Merge pull request #165 from sudarshanv01/bugfix_numbers_solver

BUGFIX: Numbers solver

Author: vossjo

catmap/mappers/min_resid_mapper.py

@@ -593,7 +593,6 @@ def minresid_iteration(isMapped):
                                                 [resid,sol_pt,cvg])
 
                                     checked_dirs[k] = 1
-
                                 if sol_cvgs:
                                     self._coverage = sol_cvgs
                                     if self.use_numbers_solver:

catmap/solvers/integrated_rate_control_solver.py

@@ -45,12 +45,12 @@ def get_integrated_DRC_rate(self,params,*args,**kwargs):
             for dname in self.descriptor_names:
                 Ef = self.species_definitions[dname]['formation_energy'][ref_idx]
                 desc_energies.append(Ef)
-            
+
             ref_dict = self.scaler.get_free_energies(desc_energies)
             self.reference_free_energies = ref_dict
             self._static_rate_controls = DRCs
             self._initialized = True
-        
+
         Gs = {}
         for key,G in zip(self.adsorbate_names + self.transition_state_names, params):
             Gs[key] = G

catmap/solvers/mean_field_solver.py

@@ -20,7 +20,7 @@ def __init__(self,reaction_model=None):
                             "stagnated or diverging (residual = ${resid})."+\
                             " Assuming Jacobian is 0.",
                             }
-        # Create a log file of the solver 
+        # Create a log file of the solver
         self.outer_solver_log = 'outer_solver.log'
 
     def get_rxn_rates(self,coverages,rate_constants):
@@ -110,7 +110,7 @@ def get_selectivity(self,rxn_parameters,weights=None):
         tofs = self.get_turnover_frequency(rxn_parameters)
         if weights is None:
             weights = [1]*len(tofs) #use for weighted selectivity (e.g. % carbon)
-            
+
         if self.products is None:
             self.products = [g for g,r in zip(self.gas_names,tofs) if r >0]
         if self.reactants is None:
@@ -241,7 +241,7 @@ def get_rxn_order(self,rxn_parameters,epsilon=1e-10):
                 dP = (new_p[i] - current_Ps[i])/current_Ps[i]
                 DRC_i.append(float(dTOF/dP))
             DRC.append(DRC_i)
-        self.gas_pressures = current_Ps 
+        self.gas_pressures = current_Ps
         self._rxn_order = DRC
         return DRC
 
@@ -299,12 +299,12 @@ def rate_equation_term(self,species_list,rate_string,d_wrt=None,specific_ads_nam
             adsorbate_names = specific_ads_names
         else:
             adsorbate_names = self.adsorbate_names
-        
+
         # Separate the species_list into gas and sites
         gas_idxs = [self.gas_names.index(gas)
                 for gas in species_list if gas in self.gas_names]
         #allows for multiple site types
-        sites = [s for s in species_list if s in self.site_names] 
+        sites = [s for s in species_list if s in self.site_names]
 
         # Create an adsorbate list by making sure that the adsorbates do not
         # include elements already present in the sites list.
@@ -314,9 +314,9 @@ def rate_equation_term(self,species_list,rate_string,d_wrt=None,specific_ads_nam
                 ads_idxs.append(adsorbate_names.index(ads))
         if len(gas_idxs+ads_idxs+sites) != len(species_list):
             raise ValueError('Undefined species in '+','.join(species_list))
-        
+
         if not d_wrt:
-            # If the derivative is not required, simply 
+            # If the derivative is not required, simply
             # return the rate string based on the pressures and coverages
             # of the relevant species
             for iden in gas_idxs:
@@ -346,17 +346,17 @@ def rate_equation_term(self,species_list,rate_string,d_wrt=None,specific_ads_nam
             # Also get the site in which the adsorbate is present
             d_site = self.species_definitions[d_wrt]['site']
 
-            if (d_wrt not in sites # Not differentiating with respect to a site 
-                and d_idx in ads_idxs # Differentiating with respect to an adsorbate 
-                and d_site not in sites # Site not present in the rate expression 
-                and self.fix_x_star == True): # Fixing x_star
+            if (d_wrt not in sites # Not differentiating with respect to a site
+                and d_idx in ads_idxs # Differentiating with respect to an adsorbate
+                and d_site not in sites # Site not present in the rate expression
+                and (self.use_numbers_solver == False or self.fix_x_star == True)): # Fixing x_star
                 # Scenario 1: Differentiating with respect to an adsorbate
-                # but there is no site term in the rate equation. 
+                # but there is no site term in the rate equation.
                 multiplier = ads_idxs.count(d_idx) #get order
                 ads_idxs.remove(d_idx) #reduce order by 1
                 if multiplier != 1:
                     rate_string = str(multiplier)+'*'+rate_string
-            
+
             elif (len(ads_idxs)>0 # Adsorbate is present in the rate expression
                   and d_site not in sites # Site not present in the rate expression
                   and d_wrt in self.site_names # Differentiating with respect to a site
@@ -367,7 +367,7 @@ def rate_equation_term(self,species_list,rate_string,d_wrt=None,specific_ads_nam
                 # of each adsorbate with respect to the site.
                 temp_rate_string = ''
                 for ie, ads_idx in enumerate(ads_idxs):
-                    # get the other adsorbates by getting 
+                    # get the other adsorbates by getting
                     # all other index of ads_idxs expect ie
                     other_ads_idxs = []
                     for je in range(len(ads_idxs)):
@@ -376,7 +376,7 @@ def rate_equation_term(self,species_list,rate_string,d_wrt=None,specific_ads_nam
                     # the rate string will just be the coverage of
                     # the other adsorbates multiplied to each other
                     # multiplied by a -1 because the derivative with
-                    # respect to the site is negative of that of the 
+                    # respect to the site is negative of that of the
                     # adsorbate.
                     if other_ads_idxs:
                         temp_rate_string += '*'.join(['theta['+str(i)+']' for i in other_ads_idxs])
@@ -387,26 +387,26 @@ def rate_equation_term(self,species_list,rate_string,d_wrt=None,specific_ads_nam
                 rate_string += '*-1*('+temp_rate_string + ')'
                 ads_idxs = [] #no more adsorbates
 
-            elif ( d_site in sites # Site present in the rate expression 
-                  and d_idx not in ads_idxs # Adsorbates not present in the rate expression 
-                  and self.fix_x_star == True): # Fixing x_star
-                # Scenario 3: There is an empty site in 
+            elif ( d_site in sites # Site present in the rate expression
+                  and d_idx not in ads_idxs # Adsorbates not present in the rate expression
+                  and (self.use_numbers_solver == False or self.fix_x_star == True)): # Fixing x_star
+                # Scenario 3: There is an empty site in
                 # the rate equation but there is no adsorbate term
                 multiplier = sites.count(d_site)
-                # Account for the fact that d/dsite does not have 
+                # Account for the fact that d/dsite does not have
                 # a -1 pre-multiplier while d/dads does
                 if d_wrt not in self.site_names:
-                    multiplier = -1 * multiplier # dsite/dtheta_* is negative 
+                    multiplier = -1 * multiplier # dsite/dtheta_* is negative
                 sites.remove(d_site) # reduce the order of site by 1
                 rate_string = str(multiplier)+'*'+rate_string
 
-            elif ( d_site in sites # Site present in the rate expression 
+            elif ( d_site in sites # Site present in the rate expression
                    and d_idx in ads_idxs # Adsorbates present in the rate expression
                    and self.fix_x_star == True): # Fixing x_star
                 # Scenario 4: Both an empty site and adsorbate in the rate equation
                 # need to use chain rule.
                 ads_mult = ads_idxs.count(d_idx)
-                # Account for the fact that d/dsite doesnt have 
+                # Account for the fact that d/dsite doesnt have
                 # a -1 pre-multiplier while d/dads does
                 if d_wrt not in self.site_names:
                     site_mult = -1 * sites.count(d_site) #negative 1 to account for d_site/d_ads
@@ -438,18 +438,18 @@ def rate_equation_term(self,species_list,rate_string,d_wrt=None,specific_ads_nam
                         [ads_str]*(ads_mult-1))+ ')'
                 rate_string += '*'+mult_rule
 
-            elif self.fix_x_star == False:
-                # Scenario 5: If x_star is not fixed, then 
-                # the partial derivatives are to be taken 
-                # independently for each site and adsorbate   
+            elif (self.use_numbers_solver == True and self.fix_x_star == False):
+                # Scenario 5: If x_star is not fixed, then
+                # the partial derivatives are to be taken
+                # independently for each site and adsorbate
                 if d_idx in ads_idxs and d_wrt not in self.site_names:
-                    # Taking a derivative with respect to an 
+                    # Taking a derivative with respect to an
                     # adsorbate coverage
                     multiplier = ads_idxs.count(d_idx) #get order
                     ads_idxs.remove(d_idx) #reduce order by 1
                     if multiplier != 1:
                         rate_string = str(multiplier)+'*'+rate_string
-                elif d_wrt in sites: 
+                elif d_wrt in sites:
                     # Taking the derivative with respect to a site
                     multiplier = sites.count(d_site)
                     sites.remove(d_site) #reduce order by 1
@@ -463,9 +463,9 @@ def rate_equation_term(self,species_list,rate_string,d_wrt=None,specific_ads_nam
             else:
                 # No dependence on either the adsorbate or the site
                 return '0'
-            
-            # Create the rate string for the derivative by taking into 
-            # account the derivative done by string manipulation in the 
+
+            # Create the rate string for the derivative by taking into
+            # account the derivative done by string manipulation in the
             # previous step
             for iden in gas_idxs:
                 rate_string += '*p['+str(iden)+']'
@@ -577,7 +577,7 @@ def rate_equations(self):
 
         if self.use_numbers_solver:
             # In case we use the numbers solver, the steady state
-            # equations needs to be n+1 dimensional 
+            # equations needs to be n+1 dimensional
             adsorbate_names = self.adsorbate_names + self.site_names
             # Remove 'g' if adsorbate_names if present
             if 'g' in adsorbate_names:
@@ -609,7 +609,7 @@ def rate_equations(self):
             surface_names = list(surface_names)
             surface_names.remove('g')
             surface_names = tuple(surface_names)
-            species_iterate = self.adsorbate_names + surface_names 
+            species_iterate = self.adsorbate_names + surface_names
         else:
             species_iterate = self.adsorbate_names
 
@@ -630,7 +630,7 @@ def rate_equations(self):
             dcdt_strings.append(dcdt_str)
 
         all_strings = rate_strings + dcdt_strings
-        
+
         return all_strings
 
     def jacobian_equations(self,adsorbate_interactions=True):
@@ -683,7 +683,7 @@ def jacobian_equations(self,adsorbate_interactions=True):
             # Standard CatMAP solver needs only an nxn Jacobian matrix
             adsorbate_names = self.adsorbate_names
 
-        if self.DEBUG: 
+        if self.DEBUG:
             with open('adsorbate_names.log','w') as f:
                 f.write(str(adsorbate_names))
 
@@ -693,10 +693,10 @@ def jacobian_equations(self,adsorbate_interactions=True):
         for i, ads_i in enumerate(adsorbate_names):
             # The second index is for the column
             for j, ads_j in enumerate(adsorbate_names):
-                # Populate the Jacobian matrix based on the 
+                # Populate the Jacobian matrix based on the
                 # required dimensions
                 J_str = 'J['+str(i)+']['+str(j)+'] = 0'
-                # Iterate over the elementary reactions such that 
+                # Iterate over the elementary reactions such that
                 # we get how much of the species are present
                 for k, rxn in enumerate(self.elementary_rxns):
                     # Get the number of ads_i in the forward and reverse
@@ -705,7 +705,7 @@ def jacobian_equations(self,adsorbate_interactions=True):
                     rxnCounts = [-1.0*rxn[0].count(ads_i),
                                   1.0*rxn[-1].count(ads_i)]
                     rxnOrder = [o for o in rxnCounts if o]
-                    # In the following, the drdx term is constructed 
+                    # In the following, the drdx term is constructed
                     # and then multiplied by the reaction order to get
                     # the df/dx term.
                     if rxnOrder:
@@ -714,7 +714,7 @@ def jacobian_equations(self,adsorbate_interactions=True):
                         rRate_string = self.rate_equation_term(rxn[-1], 'kr['+str(k)+']', ads_j, adsorbate_names)
                         if adsorbate_interactions == True:
                             if ads_j not in self.site_names:
-                                # If one of the "adsorbates" is an empty site then 
+                                # If one of the "adsorbates" is an empty site then
                                 # there should be no inclusion of adsorbate-adsorbate interaction
                                 dfRate_string = self.rate_equation_term(rxn[0],'(kfkBT['+str(k)+'])*dEf['+str(k)+']['+str(j)+']')
                                 drRate_string = self.rate_equation_term(rxn[-1],'(krkBT['+str(k)+'])*dEr['+str(k)+']['+str(j)+']')