Merge pull request #2090 for Coverage Dependent surface kinetics

Adding in coverage dependent effects to surface kinetics.
Goes with RMG-database Pull Request
ReactionMechanismGenerator/RMG-database#462

documentation/source/users/rmg/surfaces.rst

@@ -325,6 +325,61 @@ The following is a list of the current pre-packaged surface reaction libraries i
 
 
 
+Coverage Dependence
+===================
+
+An Arrhenius-like expression can be used to describe the effect of the coverage (:math:`\theta_{k}`) of species k on the forward rate constant for a given surface reaction. RMG has the capability to account for such coverage dependent kinetics, as specified in the following equation:
+
+.. math:: k_f = A T^b \exp \left( - \frac{E_a}{RT} \right)\prod_k 10^{a_k \theta_k} \theta_k^{m_k}\exp \left( \frac{- E_k \theta_k}{RT} \right)
+    :label: covdepeqn
+
+where the parameters :math:`a_k`, :math:`m_k`, and :math:`E_k` are the modifiers for the pre-exponential, temperature exponent, and activation energy from species k. 
+
+input file specifications
+-------------------------
+Coverage dependent parameters are applied if the "coverage_dependence" attribute is set to "True" in the catalyst properties block::
+
+    catalystProperties(
+        metal = Pt111
+        coverageDependence=True,
+    )
+
+By default, this field is false. 
+
+
+Coverage block in families and libraries
+----------------------------------------
+Coverage dependent parameters can be added to a family or library reaction by specifying a dictionary in the "coverage_dependence" field, with the names of the species serving as the keys. Nested within that dictionary is a dictionary of the coverage dependent parameters for each species. For example::
+
+    entry(
+        index = 5,
+        label = "NH_X + X <=> N_X + H_X",
+        kinetics = SurfaceArrhenius(
+            A = (7.22E20, 'cm^2/(mol*s)'), 
+            n = 0.0,
+            Ea = (5.3, 'kcal/mol'),
+            Tmin = (200, 'K'),
+            Tmax = (3000, 'K'),
+            coverage_dependence = {'N_X': {'a':0.0, 'm':0.0, 'E':(15.5, 'kcal/mol')},
+                                   'H_X': {'a':0.0, 'm':0.0, 'E':(1.0, 'kcal/mol')}},
+        ),
+        shortDesc = u"""Surface_Dissociation""",
+        longDesc = u"""
+    "The role of adsorbate–adsorbate interactions in the rate controlling step 
+    and the most abundant reaction intermediate of NH3 decomposition on Ru"
+    D.G. Vlachos et al. (2004). Catalysis Letters 96, 13–22. 
+    https://doi.org/10.1023/B:CATL.0000029523.22277.e1
+    This reaction used RMG's surface site density of Ru0001 = 2.630E-9(mol/cm^2) to calculate the A factor.
+    A = 1.9E12(1/s)/2.630E-9(mol/cm^2) = 7.22E20 cm^2/(mol*s)
+    This is R5 in Table 2 (set A)
+    """,
+        metal = "Ru",
+        facet = "0001",
+    )
+
+The species "N_X" and "H_X" are the coverage dependent species in this example, and they also happen to be species partaking in this reaction. However, any species that are listed in the species dictionary can also be used. 
+
+
 Examples
 =========
 

rmgpy/chemkin.pyx

@@ -288,6 +288,8 @@ def read_kinetics_entry(entry, species_dict, Aunits, Eunits):
             reaction.kinetics = kinetics['arrhenius high']
         elif 'sticking coefficient' in kinetics:
             reaction.kinetics = kinetics['sticking coefficient']
+        elif 'surface arrhenius' in kinetics:
+            reaction.kinetics = kinetics['surface arrhenius']
         else:
             raise ChemkinError(
                 'Unable to understand all additional information lines for reaction {0}.'.format(entry))
@@ -449,6 +451,24 @@ def _read_kinetics_line(line, reaction, species_dict, Eunits, kunits, klow_units
         # Duplicate reaction
         reaction.duplicate = True
 
+    elif 'COV' in line:
+        try:
+            k = kinetics['sticking coefficient']
+        except KeyError:
+            k = kinetics['arrhenius high']
+            k = _kinetics.SurfaceArrhenius(
+                A=(k.A.value, kunits),
+                n=k.n,
+                Ea=k.Ea,
+                T0=k.T0,
+            )
+            kinetics['surface arrhenius'] = k
+            del kinetics['arrhenius high']
+
+        tokens = case_preserved_tokens[1].split()
+        cov_dep_species = species_dict[tokens[0].strip()]
+        k.coverage_dependence[cov_dep_species] = {'a':float(tokens[1]), 'm':float(tokens[2]), 'E':(float(tokens[3]), Eunits)}
+
     elif 'LOW' in line:
         # Low-pressure-limit Arrhenius parameters
         tokens = tokens[1].split()
@@ -541,6 +561,7 @@ def _read_kinetics_line(line, reaction, species_dict, Eunits, kunits, klow_units
                                    Ea=(float(tokens[3].strip()), Eunits),
                                    T0=(1, "K"),
                                )])
+
     elif tokens[0].startswith('REV'):
         reverse_A = float(tokens[1].split()[0])
         kinetics['explicit reverse'] = line.strip()
@@ -549,6 +570,7 @@ def _read_kinetics_line(line, reaction, species_dict, Eunits, kunits, klow_units
             reaction.reversible = False
         else:
             logging.info("Ignoring explicit reverse rate for reaction {0}".format(reaction))
+
     elif line.strip() == 'STICK':
         # Convert what we thought was Arrhenius into StickingCoefficient
         k = kinetics['arrhenius high']
@@ -559,6 +581,7 @@ def _read_kinetics_line(line, reaction, species_dict, Eunits, kunits, klow_units
             T0=k.T0,
         )
         del kinetics['arrhenius high']
+
     else:
         # Assume a list of collider efficiencies
         try:
@@ -1844,6 +1867,13 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals
 
     string += '\n'
 
+    if getattr(kinetics, 'coverage_dependence', None):
+        # Write coverage dependence parameters for surface reactions
+        for species, cov_params in kinetics.coverage_dependence.items():
+            label = get_species_identifier(species)
+            string += f'    COV / {label:<41} '
+            string += f"{cov_params['a'].value:<9.3g} {cov_params['m'].value:<9.3g} {cov_params['E'].value_si/4184.:<9.3f} /\n"
+
     if isinstance(kinetics, (_kinetics.ThirdBody, _kinetics.Lindemann, _kinetics.Troe)):
         # Write collider efficiencies
         for collider, efficiency in sorted(list(kinetics.efficiencies.items()), key=lambda item: id(item[0])):

rmgpy/chemkinTest.py

@@ -44,6 +44,7 @@
 from rmgpy.species import Species
 from rmgpy.thermo import NASA, NASAPolynomial
 from rmgpy.transport import TransportData
+from rmgpy.kinetics.surface import SurfaceArrhenius, StickingCoefficient
 
 
 ###################################################
@@ -428,6 +429,81 @@ def test_mark_duplicate_reactions(self):
 
         self.assertEqual(duplicate_flags, expected_flags)
 
+    def test_read_coverage_dependence(self):
+        """Test that we can properly write coverage dependent parameters"""
+
+        folder = os.path.join(os.path.dirname(rmgpy.__file__), 'test_data/chemkin/chemkin_py')
+
+        s_x_entry = """X(1)                    X   1               G   100.000  5000.000 1554.80      1
+ 1.60299900E-01-2.52235409E-04 1.14181275E-07-1.21471653E-11 3.85790025E-16    2
+-7.08100885E+01-9.09527530E-01 7.10139498E-03-4.25619522E-05 8.98533016E-08    3
+-7.80193649E-11 2.32465471E-14-8.76101712E-01-3.11211229E-02                   4"""
+        s_hx_entry = """H*(10)                  H   1X   1          G   100.000  5000.000  952.91      1
+ 2.80339655E+00-5.41047017E-04 4.99507978E-07-7.54963647E-11 3.06772366E-15    2
+-2.34636021E+03-1.59436787E+01-3.80965452E-01 5.47228709E-03 2.60912778E-06    3
+-9.64961980E-09 4.63946753E-12-1.40561079E+03 1.01725550E+00                   4"""
+        s_h2_entry = """H2                      H   2               G   100.000  5000.000 1959.08      1
+ 2.78816619E+00 5.87640475E-04 1.59010635E-07-5.52739465E-11 4.34311304E-15    2
+-5.96144481E+02 1.12730527E-01 3.43536411E+00 2.12710383E-04-2.78625110E-07    3
+ 3.40267219E-10-7.76032129E-14-1.03135984E+03-3.90841731E+00                   4"""
+
+        s_h2 = Species().from_smiles("[H][H]")
+        s_x = Species().from_adjacency_list("1 X u0 p0")
+        s_x.label = 'X'
+        s_hx = Species().from_adjacency_list("1 H u0 p0 {2,S} \n 2 X u0 p0 {1,S}")
+        s_hx.label = 'HX'
+
+        species, thermo, formula = read_thermo_entry(s_x_entry)
+        s_x.thermo = thermo
+        species, thermo, formula = read_thermo_entry(s_hx_entry)
+        s_hx.thermo = thermo
+        species, thermo, formula = read_thermo_entry(s_h2_entry)
+        s_h2.thermo = thermo
+
+        species = [s_h2, s_x, s_hx]
+
+        self.rxn_covdep = Reaction(
+            reactants=[s_h2, s_x, s_x],
+            products=[s_hx, s_hx],
+            kinetics=SurfaceArrhenius(A=(9.05e18, 'cm^5/(mol^2*s)'),
+                                      n=0.5,
+                                      Ea=(5.0, 'kJ/mol'),
+                                      T0=(1.0, 'K'),
+                                      coverage_dependence={
+                                          s_x: {'a': 1.0, 'm': -1.0, 'E': (0.1, 'J/mol')},}))
+
+        reactions = [self.rxn_covdep]
+
+        # save_chemkin_file
+        chemkin_save_path = os.path.join(folder, 'surface', 'chem-surface_new.inp')
+        dictionary_save_path = os.path.join(folder, 'surface', 'species_dictionary_new.txt')
+        save_chemkin_file(chemkin_save_path, species, reactions, verbose=True, check_for_duplicates=True)
+        save_species_dictionary(dictionary_save_path, species, old_style=False)
+
+        # load chemkin file
+        species_load, reactions_load = load_chemkin_file(chemkin_save_path, dictionary_save_path)
+
+        # compare only labels because the objects are different
+        species_label = []
+        species_load_label = []
+        for s in range(len(species)):
+            species_label.append(species[s].label)
+            species_load_label.append(species_load[s].label)
+
+        # check the written chemkin file matches what is expected
+        self.assertTrue(all(i in species_load_label for i in species_label))
+
+        for r in range(len(reactions)):
+            for s in range(len(reactions[r].products)):
+                self.assertEqual(reactions[r].products[s].label,reactions_load[r].products[s].label)
+            for s in range(len(reactions[r].reactants)):
+                self.assertEqual(reactions[r].reactants[s].label, reactions_load[r].reactants[s].label)
+            self.assertIsNotNone(reactions[r].kinetics.coverage_dependence)
+            self.assertIsNotNone(reactions_load[r].kinetics.coverage_dependence)
+
+        # clean up
+        os.remove(chemkin_save_path)
+        os.remove(dictionary_save_path)
 
 class TestThermoReadWrite(unittest.TestCase):
 

rmgpy/data/kinetics/family.py

@@ -1211,7 +1211,8 @@ def add_rules_from_training(self, thermo_database=None, train_indices=None):
                     alpha=0,
                     E0=deepcopy(data.Ea),
                     Tmin=deepcopy(data.Tmin),
-                    Tmax=deepcopy(data.Tmax)
+                    Tmax=deepcopy(data.Tmax),
+                    coverage_dependence=deepcopy(data.coverage_dependence),
                 )
             elif isinstance(data, SurfaceArrhenius):
                 data = SurfaceArrheniusBEP(
@@ -1222,7 +1223,8 @@ def add_rules_from_training(self, thermo_database=None, train_indices=None):
                     alpha=0,
                     E0=deepcopy(data.Ea),
                     Tmin=deepcopy(data.Tmin),
-                    Tmax=deepcopy(data.Tmax)
+                    Tmax=deepcopy(data.Tmax),
+                    coverage_dependence=deepcopy(data.coverage_dependence),
                 )
             else:
                 raise NotImplementedError("Unexpected training kinetics type {} for {}".format(type(data), entry))
@@ -1303,7 +1305,8 @@ def add_rules_from_training(self, thermo_database=None, train_indices=None):
                     alpha=0,
                     E0=deepcopy(data.Ea),
                     Tmin=deepcopy(data.Tmin),
-                    Tmax=deepcopy(data.Tmax)
+                    Tmax=deepcopy(data.Tmax),
+                    coverage_dependence=deepcopy(data.coverage_dependence),
                 )
             else:
                 data = data.to_arrhenius_ep()

rmgpy/data/kinetics/library.py

@@ -434,6 +434,15 @@ def load(self, path, local_context=None, global_context=None):
             # Create a new reaction per entry
             rxn = entry.item
             rxn_string = entry.label
+
+            # Convert coverage dependence label to species label
+            if hasattr(entry.data, 'coverage_dependence'):
+                if entry.data.coverage_dependence:
+                    coverage_dependence_update = {}
+                    for key, value in entry.data.coverage_dependence.items():
+                        coverage_dependence_update[species_dict[key]] = value
+                    entry.data.coverage_dependence = coverage_dependence_update
+
             # Convert the reactants and products to Species objects using the species_dict
             reactants, products = rxn_string.split('=>')
             reversible = True

rmgpy/kinetics/surface.pxd

@@ -39,6 +39,7 @@ cdef class StickingCoefficient(KineticsModel):
     cdef public ScalarQuantity _n
     cdef public ScalarQuantity _Ea
     cdef public ScalarQuantity _T0
+    cdef public dict _coverage_dependence
 
     cpdef double get_sticking_coefficient(self, double T) except -1
 
@@ -56,7 +57,7 @@ cdef class StickingCoefficientBEP(KineticsModel):
     cdef public ScalarQuantity _n
     cdef public ScalarQuantity _alpha
     cdef public ScalarQuantity _E0
-
+    cdef public dict _coverage_dependence
     cpdef double get_sticking_coefficient(self, double T, double dHrxn=?) except -1
     cpdef double get_activation_energy(self, double dHrxn) except -1
     cpdef StickingCoefficient to_arrhenius(self, double dHrxn)
@@ -65,8 +66,10 @@ cdef class StickingCoefficientBEP(KineticsModel):
 
 ################################################################################
 cdef class SurfaceArrhenius(Arrhenius):
+    cdef public dict _coverage_dependence
     pass
 ################################################################################
 cdef class SurfaceArrheniusBEP(ArrheniusEP):
+    cdef public dict _coverage_dependence
     pass