diff --git a/documentation/source/users/rmg/surfaces.rst b/documentation/source/users/rmg/surfaces.rst index 12b631b46e..1b6e755707 100644 --- a/documentation/source/users/rmg/surfaces.rst +++ b/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 ========= diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index 0f0f3e55e8..9ca55a00b9 100644 --- a/rmgpy/chemkin.pyx +++ b/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])): diff --git a/rmgpy/chemkinTest.py b/rmgpy/chemkinTest.py index 271df8d0f3..cf1c006212 100644 --- a/rmgpy/chemkinTest.py +++ b/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): diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index d80c930ecf..36d024a1bb 100644 --- a/rmgpy/data/kinetics/family.py +++ b/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() diff --git a/rmgpy/data/kinetics/library.py b/rmgpy/data/kinetics/library.py index 577db37954..220b813994 100644 --- a/rmgpy/data/kinetics/library.py +++ b/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 diff --git a/rmgpy/kinetics/surface.pxd b/rmgpy/kinetics/surface.pxd index b669794e28..8e37ac54b7 100644 --- a/rmgpy/kinetics/surface.pxd +++ b/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 diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx index f97c9137bc..2c108e9204 100644 --- a/rmgpy/kinetics/surface.pyx +++ b/rmgpy/kinetics/surface.pyx @@ -47,28 +47,34 @@ cdef class StickingCoefficient(KineticsModel): Similar to :class:`Arrhenius` but with different units for `A`. The attributes are: - =============== ============================================================= - Attribute Description - =============== ============================================================= - `A` The preexponential factor - `T0` The reference temperature - `n` The temperature exponent - `Ea` The activation energy - `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined - `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined - `Pmin` The minimum pressure at which the model is valid, or zero if unknown or undefined - `Pmax` The maximum pressure at which the model is valid, or zero if unknown or undefined - `comment` Information about the model (e.g. its source) - =============== ============================================================= + ======================= ============================================================= + Attribute Description + ======================= ============================================================= + `A` The preexponential factor + `T0` The reference temperature + `n` The temperature exponent + `Ea` The activation energy + `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined + `Pmin` The minimum pressure at which the model is valid, or zero if unknown or undefined + `Pmax` The maximum pressure at which the model is valid, or zero if unknown or undefined + `coverage_dependence` A dictionary of coverage dependent parameters to a certain surface species with: + `a`, the coefficient for exponential dependence on the coverage, + `m`, the power-law exponent of coverage dependence, and + `E`, the activation energy dependence on coverage. + `comment` Information about the model (e.g. its source) + ======================= ============================================================= """ - def __init__(self, A=None, n=0.0, Ea=None, T0=(1.0, "K"), Tmin=None, Tmax=None, Pmin=None, Pmax=None, comment=''): + def __init__(self, A=None, n=0.0, Ea=None, T0=(1.0, "K"), Tmin=None, Tmax=None, Pmin=None, Pmax=None, + coverage_dependence=None, comment=''): KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, comment=comment) self.A = A self.n = n self.Ea = Ea self.T0 = T0 + self.coverage_dependence = coverage_dependence def __repr__(self): """ @@ -78,6 +84,11 @@ cdef class StickingCoefficient(KineticsModel): string = 'StickingCoefficient(A={0!r}, n={1!r}, Ea={2!r}, T0={3!r}'.format(self.A, self.n, self.Ea, self.T0) if self.Tmin is not None: string += ', Tmin={0!r}'.format(self.Tmin) if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) + if self.coverage_dependence: + string += ", coverage_dependence={" + for species, parameters in self.coverage_dependence.items(): + string += f"{species.to_chemkin()!r}: {{'a':{parameters['a']}, 'm':{parameters['m']}, 'E':({parameters['E'].value}, '{parameters['E'].units}')}}," + string += "}" if self.Pmin is not None: string += ', Pmin={0!r}'.format(self.Pmin) if self.Pmax is not None: string += ', Pmax={0!r}'.format(self.Pmax) if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) @@ -89,7 +100,7 @@ cdef class StickingCoefficient(KineticsModel): A helper function used when pickling a StickingCoefficient object. """ return (StickingCoefficient, (self.A, self.n, self.Ea, self.T0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, - self.comment)) + self.coverage_dependence, self.comment)) property A: """The preexponential factor.""" @@ -119,6 +130,19 @@ cdef class StickingCoefficient(KineticsModel): def __set__(self, value): self._T0 = quantity.Temperature(value) + property coverage_dependence: + """The coverage dependence parameters.""" + def __get__(self): + return self._coverage_dependence + def __set__(self, value): + self._coverage_dependence = {} + if value: + for species, parameters in value.items(): + processed_parameters = {'E': quantity.Energy(parameters['E']), + 'm': quantity.Dimensionless(parameters['m']), + 'a': quantity.Dimensionless(parameters['a'])} + self._coverage_dependence[species] = processed_parameters + cpdef double get_sticking_coefficient(self, double T) except -1: """ Return the sticking coefficient (dimensionless) at temperature `T` in K. @@ -128,6 +152,7 @@ cdef class StickingCoefficient(KineticsModel): n = self._n.value_si Ea = self._Ea.value_si T0 = self._T0.value_si + stickingCoefficient = A * (T / T0) ** n * exp(-Ea / (constants.R * T)) if stickingCoefficient < 0: raise ValueError("Sticking coefficients cannot be negative, check your preexponential factor.") @@ -227,28 +252,34 @@ cdef class StickingCoefficientBEP(KineticsModel): Sticking Coefficients are between 0 and 1. The attributes are: - =============== ============================================================= - Attribute Description - =============== ============================================================= - `A` The preexponential factor - `n` The temperature exponent - `alpha` The Evans-Polanyi slope - `E0` The activation energy for a thermoneutral reaction - `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined - `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined - `Pmin` The minimum pressure at which the model is valid, or zero if unknown or undefined - `Pmax` The maximum pressure at which the model is valid, or zero if unknown or undefined - `comment` Information about the model (e.g. its source) - =============== ============================================================= + ======================= ============================================================= + Attribute Description + ======================= ============================================================= + `A` The preexponential factor + `n` The temperature exponent + `alpha` The Evans-Polanyi slope + `E0` The activation energy for a thermoneutral reaction + `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined + `Pmin` The minimum pressure at which the model is valid, or zero if unknown or undefined + `Pmax` The maximum pressure at which the model is valid, or zero if unknown or undefined + `coverage_dependence` A dictionary of coverage dependent parameters to a certain surface species with: + `a`, the coefficient for exponential dependence on the coverage, + `m`, the power-law exponent of coverage dependence, and + `E`, the activation energy dependence on coverage. + `comment` Information about the model (e.g. its source) + ======================= ============================================================= """ - def __init__(self, A=None, n=0.0, alpha=0.0, E0=None, Tmin=None, Tmax=None, Pmin=None, Pmax=None, comment=''): + def __init__(self, A=None, n=0.0, alpha=0.0, E0=None, Tmin=None, Tmax=None, Pmin=None, Pmax=None, + coverage_dependence=None, comment=''): KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, comment=comment) self.A = A self.n = n self.alpha = alpha self.E0 = E0 + self.coverage_dependence = coverage_dependence def __repr__(self): """ @@ -259,6 +290,11 @@ cdef class StickingCoefficientBEP(KineticsModel): self.E0) if self.Tmin is not None: string += ', Tmin={0!r}'.format(self.Tmin) if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) + if self.coverage_dependence: + string += ", coverage_dependence={" + for species, parameters in self.coverage_dependence.items(): + string += f"{species.to_chemkin()!r}: {{'a':{parameters['a']}, 'm':{parameters['m']}, 'E':({parameters['E'].value}, '{parameters['E'].units}')}}," + string += "}" if self.Pmin is not None: string += ', Pmin={0!r}'.format(self.Pmin) if self.Pmax is not None: string += ', Pmax={0!r}'.format(self.Pmax) if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) @@ -270,7 +306,7 @@ cdef class StickingCoefficientBEP(KineticsModel): A helper function used when pickling an StickingCoefficientBEP object. """ return (StickingCoefficientBEP, (self.A, self.n, self.alpha, self.E0, self.Tmin, self.Tmax, - self.Pmin, self.Pmax, self.comment)) + self.Pmin, self.Pmax, self.coverage_dependence, self.comment)) property A: """The preexponential factor.""" @@ -300,6 +336,19 @@ cdef class StickingCoefficientBEP(KineticsModel): def __set__(self, value): self._E0 = quantity.Energy(value) + property coverage_dependence: + """The coverage dependence parameters.""" + def __get__(self): + return self._coverage_dependence + def __set__(self, value): + self._coverage_dependence = {} + if value: + for species, parameters in value.items(): + processed_parameters = {'E': quantity.Energy(parameters['E']), + 'm': quantity.Dimensionless(parameters['m']), + 'a': quantity.Dimensionless(parameters['a'])} + self._coverage_dependence[species] = processed_parameters + cpdef double get_sticking_coefficient(self, double T, double dHrxn=0.0) except -1: """ Return the sticking coefficient (dimensionless) at @@ -309,6 +358,7 @@ cdef class StickingCoefficientBEP(KineticsModel): Ea = self.get_activation_energy(dHrxn) A = self._A.value_si n = self._n.value_si + stickingCoefficient = A * T ** n * exp(-Ea / (constants.R * T)) assert 0 <= stickingCoefficient return min(stickingCoefficient, 1.0) @@ -341,6 +391,7 @@ cdef class StickingCoefficientBEP(KineticsModel): T0=(1, "K"), Tmin=self.Tmin, Tmax=self.Tmax, + coverage_dependence=self.coverage_dependence, comment=self.comment, ) @@ -382,21 +433,35 @@ cdef class SurfaceArrhenius(Arrhenius): The attributes are: - =============== ============================================================= - Attribute Description - =============== ============================================================= - `A` The preexponential factor - `T0` The reference temperature - `n` The temperature exponent - `Ea` The activation energy - `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined - `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined - `Pmin` The minimum pressure at which the model is valid, or zero if unknown or undefined - `Pmax` The maximum pressure at which the model is valid, or zero if unknown or undefined - `uncertainty` Uncertainty information - `comment` Information about the model (e.g. its source) - =============== ============================================================= + ======================= ============================================================= + Attribute Description + ======================= ============================================================= + `A` The preexponential factor + `T0` The reference temperature + `n` The temperature exponent + `Ea` The activation energy + `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined + `Pmin` The minimum pressure at which the model is valid, or zero if unknown or undefined + `Pmax` The maximum pressure at which the model is valid, or zero if unknown or undefined + `coverage_dependence` A dictionary of coverage dependent parameters to a certain surface species with: + `a`, the coefficient for exponential dependence on the coverage, + `m`, the power-law exponent of coverage dependence, and + `E`, the activation energy dependence on coverage. + `uncertainty` Uncertainty information + `comment` Information about the model (e.g. its source) + ======================= ============================================================= """ + def __init__(self, A=None, n=0.0, Ea=None, T0=(1.0, "K"), Tmin=None, Tmax=None, Pmin=None, Pmax=None, + coverage_dependence=None, uncertainty=None, comment=''): + KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, uncertainty=uncertainty, + comment=comment) + self.A = A + self.n = n + self.Ea = Ea + self.T0 = T0 + self.coverage_dependence = coverage_dependence + property A: """The preexponential factor. @@ -406,6 +471,19 @@ cdef class SurfaceArrhenius(Arrhenius): def __set__(self, value): self._A = quantity.SurfaceRateCoefficient(value) + property coverage_dependence: + """The coverage dependence parameters.""" + def __get__(self): + return self._coverage_dependence + def __set__(self, value): + self._coverage_dependence = {} + if value: + for species, parameters in value.items(): + processed_parameters = {'E': quantity.Energy(parameters['E']), + 'm': quantity.Dimensionless(parameters['m']), + 'a': quantity.Dimensionless(parameters['a'])} + self._coverage_dependence[species] = processed_parameters + def __repr__(self): """ Return a string representation that can be used to reconstruct the @@ -414,6 +492,11 @@ cdef class SurfaceArrhenius(Arrhenius): string = 'SurfaceArrhenius(A={0!r}, n={1!r}, Ea={2!r}, T0={3!r}'.format(self.A, self.n, self.Ea, self.T0) if self.Tmin is not None: string += ', Tmin={0!r}'.format(self.Tmin) if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) + if self.coverage_dependence: + string += ", coverage_dependence={" + for species, parameters in self.coverage_dependence.items(): + string += f"{species.to_chemkin()!r}: {{'a':{parameters['a']}, 'm':{parameters['m']}, 'E':({parameters['E'].value}, '{parameters['E'].units}')}}," + string += "}" if self.Pmin is not None: string += ', Pmin={0!r}'.format(self.Pmin) if self.Pmax is not None: string += ', Pmax={0!r}'.format(self.Pmax) if self.uncertainty is not None: string += ', uncertainty={0!r}'.format(self.uncertainty) @@ -426,8 +509,7 @@ cdef class SurfaceArrhenius(Arrhenius): A helper function used when pickling a SurfaceArrhenius object. """ return (SurfaceArrhenius, (self.A, self.n, self.Ea, self.T0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, - self.uncertainty, self.comment)) - + self.coverage_dependence, self.uncertainty, self.comment)) ################################################################################ @@ -439,25 +521,40 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): It is very similar to the gas-phase :class:`ArrheniusEP`. The only differences being the A factor has different units, (and the catalysis community prefers to call it BEP rather than EP!) + and has a coverage_dependence parameter for coverage dependence The attributes are: - =============== ============================================================= - Attribute Description - =============== ============================================================= - `A` The preexponential factor - `n` The temperature exponent - `alpha` The Evans-Polanyi slope - `E0` The activation energy for a thermoneutral reaction - `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined - `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined - `Pmin` The minimum pressure at which the model is valid, or zero if unknown or undefined - `Pmax` The maximum pressure at which the model is valid, or zero if unknown or undefined - `uncertainty` Uncertainty information - `comment` Information about the model (e.g. its source) - =============== ============================================================= + ======================= ============================================================= + Attribute Description + ======================= ============================================================= + `A` The preexponential factor + `n` The temperature exponent + `alpha` The Evans-Polanyi slope + `E0` The activation energy for a thermoneutral reaction + `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined + `Pmin` The minimum pressure at which the model is valid, or zero if unknown or undefined + `Pmax` The maximum pressure at which the model is valid, or zero if unknown or undefined + `coverage_dependence` A dictionary of coverage dependent parameters to a certain surface species with: + `a`, the coefficient for exponential dependence on the coverage, + `m`, the power-law exponent of coverage dependence, and + `E`, the activation energy dependence on coverage. + `uncertainty` Uncertainty information + `comment` Information about the model (e.g. its source) + ======================= ============================================================= """ + def __init__(self, A=None, n=0.0, alpha=0.0, E0=None, Tmin=None, Tmax=None, Pmin=None, Pmax=None, uncertainty=None, + coverage_dependence=None, comment=''): + KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, uncertainty=uncertainty, + comment=comment,) + self.A = A + self.n = n + self.alpha = alpha + self.E0 = E0 + self.coverage_dependence = coverage_dependence + property A: """The preexponential factor. @@ -467,6 +564,19 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): def __set__(self, value): self._A = quantity.SurfaceRateCoefficient(value) + property coverage_dependence: + """The coverage dependence parameters.""" + def __get__(self): + return self._coverage_dependence + def __set__(self, value): + self._coverage_dependence = {} + if value: + for species, parameters in value.items(): + processed_parameters = {'E': quantity.Energy(parameters['E']), + 'm': quantity.Dimensionless(parameters['m']), + 'a': quantity.Dimensionless(parameters['a'])} + self._coverage_dependence[species] = processed_parameters + def __repr__(self): """ Return a string representation that can be used to reconstruct the @@ -476,6 +586,11 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): self.E0) if self.Tmin is not None: string += ', Tmin={0!r}'.format(self.Tmin) if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) + if self.coverage_dependence: + string += ", coverage_dependence={" + for species, parameters in self.coverage_dependence.items(): + string += f"{species.to_chemkin()!r}: {{'a':{parameters['a']}, 'm':{parameters['m']}, 'E':({parameters['E'].value}, '{parameters['E'].units}')}}," + string += "}" if self.Pmin is not None: string += ', Pmin={0!r}'.format(self.Pmin) if self.Pmax is not None: string += ', Pmax={0!r}'.format(self.Pmax) if self.uncertainty is not None: string += ', uncertainty={0!r}'.format(self.uncertainty) @@ -488,7 +603,7 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): A helper function used when pickling an SurfaceArrheniusBEP object. """ return (SurfaceArrheniusBEP, (self.A, self.n, self.alpha, self.E0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, - self.uncertainty, self.comment)) + self.uncertainty, self.coverage_dependence, self.comment)) cpdef SurfaceArrhenius to_arrhenius(self, double dHrxn): """ @@ -506,6 +621,7 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): T0=(1, "K"), Tmin=self.Tmin, Tmax=self.Tmax, - uncertainty = self.uncertainty, + uncertainty=self.uncertainty, + coverage_dependence=self.coverage_dependence, comment=self.comment, ) diff --git a/rmgpy/kinetics/surfaceTest.py b/rmgpy/kinetics/surfaceTest.py index 2cab8fec9d..f5a4a13ee0 100644 --- a/rmgpy/kinetics/surfaceTest.py +++ b/rmgpy/kinetics/surfaceTest.py @@ -36,6 +36,9 @@ import numpy as np from rmgpy.kinetics.surface import StickingCoefficient, SurfaceArrhenius +from rmgpy.species import Species +from rmgpy.molecule import Molecule +import rmgpy.quantity as quantity ################################################################################ @@ -55,6 +58,11 @@ def setUp(self): self.T0 = 1. self.Tmin = 300. self.Tmax = 3000. + s = Species().from_adjacency_list('1 X u0 p0 c0') + s.label = 'X' + self.coverage_dependence = {s: {'a': quantity.Dimensionless(0.0), + 'm': quantity.Dimensionless(-1.0), + 'E': quantity.Energy(0.0, 'J/mol')}} self.comment = 'O2 dissociative' self.stick = StickingCoefficient( A=self.A, @@ -64,6 +72,7 @@ def setUp(self): Tmin=(self.Tmin, "K"), Tmax=(self.Tmax, "K"), comment=self.comment, + coverage_dependence=self.coverage_dependence, ) def test_A(self): @@ -108,6 +117,25 @@ def test_comment(self): """ self.assertEqual(self.stick.comment, self.comment) + def test_coverage_dependence(self): + """ + Test that the coverage dependent parameters was properly set. + """ + for key in self.stick.coverage_dependence.keys(): + match = False + for key2 in self.coverage_dependence.keys(): + if key.is_identical(key2): + match = True + self.assertTrue(match) + for species, parameters in self.stick.coverage_dependence.items(): + match = False + for species2 in self.coverage_dependence.keys(): + if species.is_identical(species2): + match = True + for key, value in parameters.items(): + self.assertEqual(value.value_si, self.coverage_dependence[species2][key].value_si) + self.assertTrue(match) + def test_is_temperature_valid(self): """ Test the StickingCoefficient.is_temperature_valid() method. @@ -137,6 +165,20 @@ def test_pickle(self): self.assertAlmostEqual(self.stick.Tmax.value, stick.Tmax.value, 4) self.assertEqual(self.stick.Tmax.units, stick.Tmax.units) self.assertEqual(self.stick.comment, stick.comment) + for key in self.stick.coverage_dependence.keys(): + match = False + for key2 in stick.coverage_dependence.keys(): + if key.is_identical(key2): + match = True + self.assertTrue(match) + for species, parameters in self.stick.coverage_dependence.items(): + match = False + for species2 in stick.coverage_dependence.keys(): + if species.is_identical(species2): + match = True + for key, value in parameters.items(): + self.assertEqual(value.value_si, stick.coverage_dependence[species2][key].value_si) + self.assertTrue(match) self.assertEqual(dir(self.stick), dir(stick)) def test_repr(self): @@ -145,7 +187,7 @@ def test_repr(self): output with no loss of information. """ namespace = {} - exec('stick = {0!r}'.format(self.stick), globals(), namespace) + exec(f'stick = {self.stick!r}', globals(), namespace) self.assertIn('stick', namespace) stick = namespace['stick'] self.assertAlmostEqual(self.stick.A.value, stick.A.value, delta=1e0) @@ -160,6 +202,10 @@ def test_repr(self): self.assertAlmostEqual(self.stick.Tmax.value, stick.Tmax.value, 4) self.assertEqual(self.stick.Tmax.units, stick.Tmax.units) self.assertEqual(self.stick.comment, stick.comment) + self.assertEqual([m.label for m in self.stick.coverage_dependence.keys()], list(stick.coverage_dependence.keys())) + for species, parameters in self.stick.coverage_dependence.items(): + for key, value in parameters.items(): + self.assertEqual(value.value_si, stick.coverage_dependence[species.label][key].value_si) self.assertEqual(dir(self.stick), dir(stick)) def test_copy(self): @@ -181,6 +227,20 @@ def test_copy(self): self.assertAlmostEqual(self.stick.Tmax.value, stick.Tmax.value, 4) self.assertEqual(self.stick.Tmax.units, stick.Tmax.units) self.assertEqual(self.stick.comment, stick.comment) + for key in self.stick.coverage_dependence.keys(): + match = False + for key2 in stick.coverage_dependence.keys(): + if key.is_identical(key2): + match = True + self.assertTrue(match) + for species, parameters in self.stick.coverage_dependence.items(): + match = False + for species2 in stick.coverage_dependence.keys(): + if species.is_identical(species2): + match = True + for key, value in parameters.items(): + self.assertEqual(value.value_si, stick.coverage_dependence[species2][key].value_si) + self.assertTrue(match) self.assertEqual(dir(self.stick), dir(stick)) def test_is_identical_to(self): @@ -207,6 +267,11 @@ def setUp(self): self.T0 = 1. self.Tmin = 300. self.Tmax = 3000. + s = Species().from_adjacency_list('1 X u0 p0 c0') + s.label = 'X' + self.coverage_dependence = {s: {'a': quantity.Dimensionless(0.0), + 'm': quantity.Dimensionless(-1.0), + 'E': quantity.Energy(0.0, 'J/mol')}} self.comment = 'CH3x + Hx <=> CH4 + x + x' self.surfarr = SurfaceArrhenius( A=(self.A, 'm^2/(mol*s)'), @@ -216,6 +281,7 @@ def setUp(self): Tmin=(self.Tmin, "K"), Tmax=(self.Tmax, "K"), comment=self.comment, + coverage_dependence=self.coverage_dependence, ) def test_A(self): @@ -260,6 +326,25 @@ def test_comment(self): """ self.assertEqual(self.surfarr.comment, self.comment) + def test_coverage_dependence(self): + """ + Test that the coverage dependent parameters was properly set. + """ + for key in self.surfarr.coverage_dependence.keys(): + match = False + for key2 in self.coverage_dependence.keys(): + if key.is_identical(key2): + match = True + self.assertTrue(match) + for species, parameters in self.surfarr.coverage_dependence.items(): + match = False + for species2 in self.coverage_dependence.keys(): + if species.is_identical(species2): + match = True + for key, value in parameters.items(): + self.assertEqual(value.value_si, self.coverage_dependence[species2][key].value_si) + self.assertTrue(match) + def test_is_temperature_valid(self): """ Test the SurfaceArrhenius.is_temperature_valid() method. @@ -289,6 +374,20 @@ def test_pickle(self): self.assertAlmostEqual(self.surfarr.Tmax.value, surfarr.Tmax.value, 4) self.assertEqual(self.surfarr.Tmax.units, surfarr.Tmax.units) self.assertEqual(self.surfarr.comment, surfarr.comment) + for key in self.surfarr.coverage_dependence.keys(): + match = False + for key2 in surfarr.coverage_dependence.keys(): + if key.is_identical(key2): + match = True + self.assertTrue(match) + for species, parameters in self.surfarr.coverage_dependence.items(): + match = False + for species2 in surfarr.coverage_dependence.keys(): + if species.is_identical(species2): + match = True + for key, value in parameters.items(): + self.assertEqual(value.value_si, surfarr.coverage_dependence[species2][key].value_si) + self.assertTrue(match) self.assertEqual(dir(self.surfarr), dir(surfarr)) def test_repr(self): @@ -312,6 +411,10 @@ def test_repr(self): self.assertAlmostEqual(self.surfarr.Tmax.value, surfarr.Tmax.value, 4) self.assertEqual(self.surfarr.Tmax.units, surfarr.Tmax.units) self.assertEqual(self.surfarr.comment, surfarr.comment) + self.assertEqual([m.label for m in self.surfarr.coverage_dependence.keys()], list(surfarr.coverage_dependence.keys())) + for species, parameters in self.surfarr.coverage_dependence.items(): + for key, value in parameters.items(): + self.assertEqual(value.value_si, surfarr.coverage_dependence[species.label][key].value_si) self.assertEqual(dir(self.surfarr), dir(surfarr)) def test_copy(self): @@ -333,6 +436,20 @@ def test_copy(self): self.assertAlmostEqual(self.surfarr.Tmax.value, surfarr.Tmax.value, 4) self.assertEqual(self.surfarr.Tmax.units, surfarr.Tmax.units) self.assertEqual(self.surfarr.comment, surfarr.comment) + for key in self.surfarr.coverage_dependence.keys(): + match = False + for key2 in surfarr.coverage_dependence.keys(): + if key.is_identical(key2): + match = True + self.assertTrue(match) + for species, parameters in self.surfarr.coverage_dependence.items(): + match = False + for species2 in surfarr.coverage_dependence.keys(): + if species.is_identical(species2): + match = True + for key, value in parameters.items(): + self.assertEqual(value.value_si, surfarr.coverage_dependence[species2][key].value_si) + self.assertTrue(match) self.assertEqual(dir(self.surfarr), dir(surfarr)) def test_is_identical_to(self): diff --git a/rmgpy/reactionTest.py b/rmgpy/reactionTest.py index f0e7aaea6e..64e8eb002a 100644 --- a/rmgpy/reactionTest.py +++ b/rmgpy/reactionTest.py @@ -139,7 +139,8 @@ class TestSurfaceReaction(unittest.TestCase): def setUp(self): m_h2 = Molecule().from_smiles("[H][H]") m_x = Molecule().from_adjacency_list("1 X u0 p0") - m_hx = Molecule().from_adjacency_list("1 H u0 p0 {2,S} \n 2 X u0 p0 {1,S}") + m_hx = Molecule().from_smiles("[H][*]") + # m_hx = Molecule().from_adjacency_list("1 H u0 p0 {2,S} \n 2 X u0 p0 {1,S}") m_ch3 = Molecule().from_smiles("[CH3]") m_ch3x = Molecule().from_adjacency_list("""1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} 2 H u0 p0 c0 {1,S} @@ -148,6 +149,7 @@ def setUp(self): 5 X u0 p0 c0 {1,S}""") s_h2 = Species( + label="H2(1)", molecule=[m_h2], thermo=ThermoData(Tdata=([300, 400, 500, 600, 800, 1000, 1500, 2000], "K"), @@ -156,6 +158,7 @@ def setUp(self): H298=(0, "kcal/mol"), S298=(31.129, "cal/(mol*K)"))) s_x = Species( + label="X(2)", molecule=[m_x], thermo=ThermoData(Tdata=([300, 400, 500, 600, 800, 1000, 1500, 2000], "K"), @@ -163,6 +166,7 @@ def setUp(self): H298=(0.0, "kcal/mol"), S298=(0.0, "cal/(mol*K)"))) s_hx = Species( + label="HX(3)", molecule=[m_hx], thermo=ThermoData(Tdata=([300, 400, 500, 600, 800, 1000, 1500, 2000], "K"), @@ -171,6 +175,7 @@ def setUp(self): S298=(0.44, "cal/(mol*K)"))) s_ch3 = Species( + label="CH3(4)", molecule=[m_ch3], thermo=NASA(polynomials=[ NASAPolynomial(coeffs=[3.91547, 0.00184155, 3.48741e-06, -3.32746e-09, 8.49953e-13, 16285.6, 0.351743], @@ -184,6 +189,7 @@ def setUp(self): ) s_ch3x = Species( + label="CH3X(5)", molecule=[m_ch3x], thermo=NASA(polynomials=[ NASAPolynomial(coeffs=[-0.552219, 0.026442, -3.55617e-05, 2.60044e-08, -7.52707e-12, -4433.47, 0.692144], @@ -225,6 +231,17 @@ def setUp(self): comment="""Approximate rate""") ) + # adding coverage dependent reaction + 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={Species().from_adjacency_list('1 X u0 p0 c0'): {'E': (0.1, 'J/mol'), 'm': -1.0, 'a': 1.0}, + Species().from_adjacency_list('1 X u0 p0 c0 {2,S}\n2 H u0 p0 c0 {1,S}'): {'E': (0.2, 'J/mol'), 'm': -2.0, 'a': 2.0}})) + def test_is_surface_reaction_species(self): """Test is_surface_reaction for reaction based on Species """ self.assertTrue(self.rxn1s.is_surface_reaction()) diff --git a/rmgpy/rmg/input.py b/rmgpy/rmg/input.py index 4017004b6e..23a58c7546 100644 --- a/rmgpy/rmg/input.py +++ b/rmgpy/rmg/input.py @@ -99,12 +99,14 @@ def database( def catalyst_properties(bindingEnergies=None, surfaceSiteDensity=None, - metal=None): + metal=None, + coverageDependence=False): """ Specify the properties of the catalyst. Binding energies of C,H,O,N atoms, and the surface site density. Metal is the label of a metal in the surface metal library. Defaults to Pt(111) if not specified. + coverageDependence defaults to False if not specified. If True then coverage dependent kinetics from libraries are used. """ # Normally we wouldn't load the database until after reading the input file, # but we need to load the metal surfaces library to validate the input. @@ -140,6 +142,11 @@ def catalyst_properties(bindingEnergies=None, logging.info("Using binding energies:\n%r", rmg.binding_energies) logging.info("Using surface site density: %r", rmg.surface_site_density) + if coverageDependence: + logging.info("Coverage dependence is turned ON") + else: + logging.info("Coverage dependence is turned OFF") + rmg.coverage_dependence = coverageDependence def convert_binding_energies(binding_energies): """ @@ -509,7 +516,8 @@ def surface_reactor(temperature, termination=termination, sensitive_species=sensitive_species, sensitivity_threshold=sensitivityThreshold, - sens_conditions=sens_conditions) + sens_conditions=sens_conditions, + coverage_dependence=rmg.coverage_dependence) rmg.reaction_systems.append(system) system.log_initial_conditions(number=len(rmg.reaction_systems)) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index b986f915c1..deaf99304a 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -185,6 +185,7 @@ def clear(self): self.diffusion_limiter = None self.surface_site_density = None self.binding_energies = None + self.coverage_dependence = False self.reaction_model = None self.reaction_systems = None @@ -264,6 +265,7 @@ def load_input(self, path=None): if self.surface_site_density: self.reaction_model.surface_site_density = self.surface_site_density + self.reaction_model.coverage_dependence = self.coverage_dependence self.reaction_model.verbose_comments = self.verbose_comments self.reaction_model.save_edge_species = self.save_edge_species diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index b2b3e608c1..840d7d1550 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -649,10 +649,13 @@ def enlarge(self, new_object=None, react_edge=False, logging.info('Generating kinetics for new reactions...') for reaction in self.new_reaction_list: # If the reaction already has kinetics (e.g. from a library), - # assume the kinetics are satisfactory + # assume the kinetics are satisfactory, else generate them if reaction.kinetics is None: self.apply_kinetics_to_reaction(reaction) + if getattr(reaction.kinetics, 'coverage_dependence', None): + self.process_coverage_dependence(reaction.kinetics) + # For new reactions, convert ArrheniusEP to Arrhenius, and fix barrier heights. # self.new_reaction_list only contains *actually* new reactions, all in the forward direction. for reaction in self.new_reaction_list: @@ -865,6 +868,25 @@ def generate_thermo(self, spc, rename=False): spc.generate_energy_transfer_model() + def process_coverage_dependence(self, kinetics): + """Process the coverage dependence kinetics. + + This ensures that species that are used in coverage-dependent kinetic + expressions exist in the model. (Before this is called, they may have + only existed in a reaction libary instance). + + If self.coverage_dependence is False then + it instead removes any coverage_dependence from the kinetics. + """ + if not self.coverage_dependence: + kinetics.coverage_dependence = None + return + cov_dep = {} + for species, values in kinetics.coverage_dependence.items(): + species_in_model, is_new = self.make_new_species(species) + cov_dep[species_in_model] = values + kinetics.coverage_dependence = cov_dep + def apply_kinetics_to_reaction(self, reaction): """ retrieve the best kinetics for the reaction and apply it towards the forward @@ -1496,6 +1518,8 @@ def add_seed_mechanism_to_core(self, seed_mechanism, react=False): reversible=rxn.reversible ) r, isNew = self.make_new_reaction(rxn) # updates self.new_species_list and self.newReactionlist + if getattr(r.kinetics, 'coverage_dependence', None): + self.process_coverage_dependence(r.kinetics) if not isNew: logging.info("This library reaction was not new: {0}".format(rxn)) elif self.pressure_dependence and rxn.elementary_high_p and rxn.is_unimolecular() \ @@ -1595,6 +1619,8 @@ def add_reaction_library_to_edge(self, reaction_library): reversible=rxn.reversible ) r, isNew = self.make_new_reaction(rxn) # updates self.new_species_list and self.newReactionlist + if getattr(r.kinetics, 'coverage_dependence', None): + self.process_coverage_dependence(r.kinetics) if not isNew: logging.info("This library reaction was not new: {0}".format(rxn)) elif self.pressure_dependence and rxn.elementary_high_p and rxn.is_unimolecular() \ diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index 4fcf338c17..07821142f1 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -36,6 +36,7 @@ import logging cimport cython import numpy as np cimport numpy as np +from libc.math cimport exp import rmgpy.constants as constants cimport rmgpy.constants as constants @@ -43,7 +44,6 @@ from rmgpy.quantity import Quantity from rmgpy.quantity cimport ScalarQuantity from rmgpy.solver.base cimport ReactionSystem - cdef class SurfaceReactor(ReactionSystem): """ A reaction system consisting of a heterogeneous, isothermal, constant volume batch @@ -67,6 +67,10 @@ cdef class SurfaceReactor(ReactionSystem): cdef public np.ndarray reactions_on_surface # (catalyst surface, not core/edge surface) cdef public np.ndarray species_on_surface # (catalyst surface, not core/edge surface) + cdef public bint coverage_dependence + cdef public dict coverage_dependencies + + def __init__(self, T, P_initial, @@ -79,6 +83,7 @@ cdef class SurfaceReactor(ReactionSystem): sensitive_species=None, sensitivity_threshold=1e-3, sens_conditions=None, + coverage_dependence=False, ): ReactionSystem.__init__(self, termination, @@ -97,11 +102,14 @@ cdef class SurfaceReactor(ReactionSystem): self.initial_surface_coverages = initial_surface_coverages self.surface_volume_ratio = Quantity(surface_volume_ratio) self.surface_site_density = Quantity(surface_site_density) + self.coverage_dependence = coverage_dependence self.V = 0 # will be set from ideal gas law in initialize_model self.constant_volume = True self.sens_conditions = sens_conditions self.n_sims = n_sims + self.coverage_dependencies = {} + def convert_initial_keys_to_species_objects(self, species_dict): """ Convert the initial_gas_mole_fractions and initial_surface_coverages dictionaries @@ -166,9 +174,27 @@ cdef class SurfaceReactor(ReactionSystem): continue if spec.contains_surface_site(): species_on_surface[index] = 1 - for rxn, index in self.reaction_index.items(): + self.coverage_dependencies = {} + for rxn, rxn_index in self.reaction_index.items(): if rxn.is_surface_reaction(): - reactions_on_surface[index] = 1 + reactions_on_surface[rxn_index] = 1 + if self.coverage_dependence and rxn.kinetics.coverage_dependence: + for spec, parameters in rxn.kinetics.coverage_dependence.items(): + species_index = self.species_index[spec] + try: + list_of_coverage_deps = self.coverage_dependencies[species_index] + except KeyError: # doesn't exist yet + list_of_coverage_deps = [] + self.coverage_dependencies[species_index] = list_of_coverage_deps + list_of_coverage_deps.append((rxn_index, + parameters['a'].value_si, + parameters['m'].value_si, + parameters['E'].value_si)) + """ + self.coverage_dependencies[2] = [(3, 0.1, -1.0, 12000.0),] + means that Species with index 2 in the current simulation is used in + Reaction 3 with parameters a=0.1, m=-1, E=12 kJ/mol + """ self.species_on_surface = species_on_surface self.reactions_on_surface = reactions_on_surface @@ -206,8 +232,6 @@ cdef class SurfaceReactor(ReactionSystem): for rxn in itertools.chain(core_reactions, edge_reactions): j = self.reaction_index[rxn] - # ToDo: get_rate_coefficient should also depend on surface coverages vector - if rxn.is_surface_reaction(): """ Be careful! From here on kf and kb will now be in Volume units, @@ -344,7 +368,7 @@ cdef class SurfaceReactor(ReactionSystem): """ cdef np.ndarray[np.int_t, ndim=2] ir, ip, inet cdef np.ndarray[np.int_t, ndim=1] reactions_on_surface, species_on_surface - cdef np.ndarray[np.float64_t, ndim=1] res, kf, kr, knet, delta, equilibrium_constants + cdef np.ndarray[np.float64_t, ndim=1] res, kf, kr, knet, delta, equilibrium_constants, coverage_corrections cdef Py_ssize_t num_core_species, num_core_reactions, num_edge_species, num_edge_reactions, num_pdep_networks cdef Py_ssize_t i, j, z, first, second, third cdef double k, V, reaction_rate, surface_volume_ratio_si @@ -352,6 +376,10 @@ cdef class SurfaceReactor(ReactionSystem): cdef np.ndarray[np.float64_t, ndim=1] edge_species_rates, edge_reaction_rates, network_leak_rates cdef np.ndarray[np.float64_t, ndim=1] C cdef np.ndarray[np.float64_t, ndim=2] jacobian, dgdk + cdef list list_of_coverage_deps + cdef double surface_site_fraction, total_sites, a, m, E + + ir = self.reactant_indices ip = self.product_indices @@ -384,6 +412,8 @@ cdef class SurfaceReactor(ReactionSystem): C = np.zeros_like(self.core_species_concentrations) V = self.V # constant volume reactor + A = self.V * surface_volume_ratio_si # area + total_sites = self.surface_site_density.value_si * A # todo: double check units for j in range(num_core_species): if species_on_surface[j]: @@ -393,6 +423,26 @@ cdef class SurfaceReactor(ReactionSystem): #: surface species are in mol/m2, gas phase are in mol/m3 core_species_concentrations[j] = C[j] + # Coverage dependence + coverage_corrections = np.ones_like(kf, np.float64) + if self.coverage_dependence: + """ + self.coverage_dependencies[2] = [(3, 0.1, -1.0, 12000.0),] + means that Species with index 2 in the current simulation is used in + Reaction 3 with parameters a=0.1, m=-1, E=12 kJ/mol + """ + for i, list_of_coverage_deps in self.coverage_dependencies.items(): + # Species i, Reaction j + surface_site_fraction = N[i] / total_sites + if surface_site_fraction < 1e-15: + continue + for j, a, m, E in list_of_coverage_deps: + coverage_corrections[j] *= 10. ** (a * surface_site_fraction) *\ + surface_site_fraction ** m *\ + exp(-1 * E * surface_site_fraction / (constants.R * self.T.value_si)) + kf = kf * coverage_corrections # make a corrected copy kf, but leave the original array at self.kf unchanged + kr = kr * coverage_corrections + for j in range(ir.shape[0]): k = kf[j] if ir[j, 0] >= num_core_species or ir[j, 1] >= num_core_species or ir[j, 2] >= num_core_species: diff --git a/rmgpy/solver/surfaceTest.py b/rmgpy/solver/surfaceTest.py index b52cdb2dbc..1cc6452d21 100644 --- a/rmgpy/solver/surfaceTest.py +++ b/rmgpy/solver/surfaceTest.py @@ -28,6 +28,7 @@ ############################################################################### import unittest +import time import numpy as np @@ -83,7 +84,7 @@ def test_solve_h2(self): core_reactions = [rxn1] edge_reactions = [] - T = 1000 + T = 600 P_initial = 1.0e5 rxn_system = SurfaceReactor( T, P_initial, @@ -117,7 +118,7 @@ def test_solve_h2(self): y = np.array(y, np.float64) reaction_rates = np.array(reaction_rates, np.float64) species_rates = np.array(species_rates, np.float64) - V = constants.R * rxn_system.T.value_si * np.sum(y) / rxn_system.P_initial.value_si + total_sites = y[0,1] # Check that we're computing the species fluxes correctly for i in range(t.shape[0]): @@ -135,8 +136,8 @@ def test_solve_h2(self): # import pylab # fig = pylab.figure(figsize=(6, 6)) # pylab.subplot(2, 1, 1) - # pylab.semilogx(t, y[:, 2]) - # pylab.ylabel('Concentration (mol/m$^\\mathdefault{3 or 2}$)') + # pylab.semilogx(t, y[:, 2] / total_sites) + # pylab.ylabel('Surface coverage') # pylab.legend(['HX'], loc=4) # pylab.subplot(2, 1, 2) # pylab.semilogx(t, species_rates) @@ -147,8 +148,7 @@ def test_solve_h2(self): # pylab.tight_layout() # # pylab.show() # pylab.savefig('surfaceTestH2.pdf') - - return + # return def test_solve_ch3(self): """ @@ -280,3 +280,306 @@ def test_solve_ch3(self): # Check that we've reached equilibrium by the end self.assertAlmostEqual(reaction_rates[-1, 0], 0.0, delta=1e-2) + + def test_solve_h2_coverage_dependence(self): + """ + Test the surface batch reactor can properly apply coverage dependent parameters + with the dissociative adsorption of H2. + + Here we choose a kinetic model consisting of the dissociative adsorption reaction + H2 + 2X <=> 2 HX + We use a SurfaceArrhenius for the rate expression. + """ + h2 = Species( + molecule=[Molecule().from_smiles("[H][H]")], + thermo=ThermoData(Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"), + Cpdata=([6.955, 6.955, 6.956, 6.961, 7.003, 7.103, 7.502], "cal/(mol*K)"), + H298=(0, "kcal/mol"), + S298=(31.129, "cal/(mol*K)"))) + x = Species( + molecule=[Molecule().from_adjacency_list("1 X u0 p0")], + thermo=ThermoData(Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"), + Cpdata=([0., 0., 0., 0., 0., 0., 0.], "cal/(mol*K)"), + H298=(0.0, "kcal/mol"), + S298=(0.0, "cal/(mol*K)"))) + hx = Species( + molecule=[Molecule().from_adjacency_list("1 H u0 p0 {2,S} \n 2 X u0 p0 {1,S}")], + thermo=ThermoData(Tdata=([300, 400, 500, 600, 800, 1000, 1500], "K"), + Cpdata=([1.50, 2.58, 3.40, 4.00, 4.73, 5.13, 5.57], "cal/(mol*K)"), + H298=(-11.26, "kcal/mol"), + S298=(0.44, "cal/(mol*K)"))) + + rxn1 = Reaction(reactants=[h2, x, x], + products=[hx, 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={x: {'a': 0.0, 'm': -1.0, 'E': (0.0, 'J/mol')}})) + + rxn2 = Reaction(reactants=[h2, x, x], + products=[hx, hx], + kinetics=SurfaceArrhenius(A=(9.05e-18, 'cm^5/(mol^2*s)'), # 1e36 times slower + n=0.5, + Ea=(5.0, 'kJ/mol'), + T0=(1.0, 'K'), + coverage_dependence={x: {'a': 0.0, 'm': -1.0, 'E': (10.0, 'J/mol')}} + )) + + core_species = [h2, x, hx] + edge_species = [] + core_reactions = [rxn1] + edge_reactions = [] + + # make it slower, for benchmarking + for j in range(200): + core_species.append(hx.copy()) + for j in range(1000): + core_reactions.append(rxn2) + + T = 600 + P_initial = 1.0e5 + rxn_system = SurfaceReactor( + T, P_initial, + n_sims=1, + initial_gas_mole_fractions={h2: 1.0}, + initial_surface_coverages={x: 1.0}, + surface_volume_ratio=(1e1, 'm^-1'), + surface_site_density=(2.72e-9, 'mol/cm^2'), + coverage_dependence=True, + termination=[]) + + rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions) + + tlist = np.logspace(-13, -5, 81, dtype=np.float64) + + self.assertIsInstance(rxn1.kinetics.coverage_dependence, dict) # check to make sure coverage_dependence is still the correct type + for species, parameters in rxn1.kinetics.coverage_dependence.items(): + self.assertIsInstance(species, Species) # species should be a Species + self.assertIsInstance(parameters, dict) + self.assertIsNotNone(parameters['a']) + self.assertIsNotNone(parameters['m']) + self.assertIsNotNone(parameters['E']) + + # Integrate to get the solution at each time point + t = [] + y = [] + reaction_rates = [] + species_rates = [] + start_time = time.time() + for t1 in tlist: + rxn_system.advance(t1) + t.append(rxn_system.t) + # You must make a copy of y because it is overwritten by DASSL at + # each call to advance() + y.append(rxn_system.y.copy()) + reaction_rates.append(rxn_system.core_reaction_rates.copy()) + species_rates.append(rxn_system.core_species_rates.copy()) + run_time = time.time() - start_time + print(f"Simulation took {run_time:.3e} seconds in {self.id()}") + + # Convert the solution vectors to np arrays + t = np.array(t, np.float64) + y = np.array(y, np.float64) + reaction_rates = np.array(reaction_rates, np.float64) + species_rates = np.array(species_rates, np.float64) + total_sites = y[0,1] + + # Check that we're computing the species fluxes correctly + for i in range(t.shape[0]): + self.assertAlmostEqual(reaction_rates[i, 0], -1.0 * species_rates[i, 0], + delta=1e-6 * reaction_rates[i, 0]) + self.assertAlmostEqual(reaction_rates[i, 0], -0.5 * species_rates[i, 1], + delta=1e-6 * reaction_rates[i, 0]) + self.assertAlmostEqual(reaction_rates[i, 0], 0.5 * species_rates[i, 2], + delta=1e-6 * reaction_rates[i, 0]) + + # Check that we've reached equilibrium + self.assertAlmostEqual(reaction_rates[-1, 0], 0.0, delta=1e-2) + + def test_solve_ch3_coverage_dependence(self): + """ + Test the surface batch reactor can properly apply coverage dependent parameters + with the nondissociative adsorption of CH3 + + Here we choose a kinetic model consisting of the adsorption reaction + CH3 + X <=> CH3X + We use a sticking coefficient for the rate expression. + """ + + ch3 = Species( + molecule=[Molecule().from_smiles("[CH3]")], + thermo=NASA( + polynomials=[ + NASAPolynomial( + coeffs=[3.91547, 0.00184155, 3.48741e-06, -3.32746e-09, 8.49953e-13, 16285.6, 0.351743], + Tmin=(100, 'K'), Tmax=(1337.63, 'K')), + NASAPolynomial( + coeffs=[3.54146, 0.00476786, -1.82148e-06, 3.28876e-10, -2.22545e-14, 16224, 1.66032], + Tmin=(1337.63, 'K'), Tmax=(5000, 'K'))], + Tmin=(100, 'K'), Tmax=(5000, 'K'), E0=(135.382, 'kJ/mol'), + comment="""Thermo library: primaryThermoLibrary + radical(CH3)""" + ), + molecular_weight=(15.0345, 'amu'), + ) + + x = Species( + molecule=[Molecule().from_adjacency_list("1 X u0 p0")], + thermo=NASA(polynomials=[NASAPolynomial(coeffs=[0, 0, 0, 0, 0, 0, 0], Tmin=(298, 'K'), Tmax=(1000, 'K')), + NASAPolynomial(coeffs=[0, 0, 0, 0, 0, 0, 0], Tmin=(1000, 'K'), Tmax=(2000, 'K'))], + Tmin=(298, 'K'), Tmax=(2000, 'K'), E0=(-6.19426, 'kJ/mol'), + comment="""Thermo library: surfaceThermo""") + ) + + ch3x = Species( + molecule=[Molecule().from_adjacency_list("""1 C u0 p0 c0 {2,S} {3,S} {4,S} {5,S} +2 H u0 p0 c0 {1,S} +3 H u0 p0 c0 {1,S} +4 H u0 p0 c0 {1,S} +5 X u0 p0 c0 {1,S}""")], + thermo=NASA( + polynomials=[ + NASAPolynomial( + coeffs=[-0.552219, 0.026442, -3.55617e-05, 2.60044e-08, -7.52707e-12, -4433.47, 0.692144], + Tmin=(298, 'K'), Tmax=(1000, 'K')), + NASAPolynomial( + coeffs=[3.62557, 0.00739512, -2.43797e-06, 1.86159e-10, 3.6485e-14, -5187.22, -18.9668], + Tmin=(1000, 'K'), Tmax=(2000, 'K'))], + Tmin=(298, 'K'), Tmax=(2000, 'K'), E0=(-39.1285, 'kJ/mol'), + comment="""Thermo library: surfaceThermoNi111""") + ) + + rxn1 = Reaction(reactants=[ch3, x], + products=[ch3x], + kinetics=StickingCoefficient( + A=0.1, n=0, Ea=(0, 'kcal/mol'), T0=(1, 'K'), Tmin=(200, 'K'), Tmax=(3000, 'K'), + coverage_dependence={x: {'a': 0.0, 'm': -1.0, 'E': (0.0, 'J/mol')}}, + comment="""Exact match found for rate rule (Adsorbate;VacantSite)""" + ) + ) + core_species = [ch3, x, ch3x] + edge_species = [] + core_reactions = [rxn1] + edge_reactions = [] + + T = 800. + P_initial = 1.0e5 + rxn_system = SurfaceReactor( + T, P_initial, + n_sims=1, + initial_gas_mole_fractions={ch3: 1.0}, + initial_surface_coverages={x: 1.0}, + surface_volume_ratio=(1., 'm^-1'), + surface_site_density=(2.72e-9, 'mol/cm^2'), + coverage_dependence=True, + termination=[]) + # in chemkin, the sites are mostly occupied in about 1e-8 seconds. + + rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions) + + tlist = np.logspace(-13, -5, 81, dtype=np.float64) + + print("Surface site density:", rxn_system.surface_site_density.value_si) + + print("rxn1 rate coefficient", + rxn1.get_surface_rate_coefficient(rxn_system.T.value_si, rxn_system.surface_site_density.value_si)) + + self.assertIsInstance(rxn1.kinetics.coverage_dependence, dict) # check to make sure coverage_dependence is still the correct type + for species, parameters in rxn1.kinetics.coverage_dependence.items(): + self.assertIsInstance(species, Species) # species should be a Species + self.assertIsInstance(parameters, dict) + self.assertIsNotNone(parameters['a']) + self.assertIsNotNone(parameters['m']) + self.assertIsNotNone(parameters['E']) + + # Integrate to get the solution at each time point + t = [] + y = [] + reaction_rates = [] + species_rates = [] + t.append(rxn_system.t) + # You must make a copy of y because it is overwritten by DASSL at + # each call to advance() + y.append(rxn_system.y.copy()) + reaction_rates.append(rxn_system.core_reaction_rates.copy()) + species_rates.append(rxn_system.core_species_rates.copy()) + print("time: ", t) + print("moles:", y) + print("reaction rates:", reaction_rates) + print("species rates:", species_rates) + start_time = time.time() + for t1 in tlist: + rxn_system.advance(t1) + t.append(rxn_system.t) + # You must make a copy of y because it is overwritten by DASSL at + # each call to advance() + y.append(rxn_system.y.copy()) + reaction_rates.append(rxn_system.core_reaction_rates.copy()) + species_rates.append(rxn_system.core_species_rates.copy()) + run_time = time.time() - start_time + print(f"Simulation took {run_time:.3e} seconds in {self.id()}") + + # Convert the solution vectors to np arrays + t = np.array(t, np.float64) + y = np.array(y, np.float64) + reaction_rates = np.array(reaction_rates, np.float64) + species_rates = np.array(species_rates, np.float64) + V = constants.R * rxn_system.T.value_si * np.sum(y) / rxn_system.P_initial.value_si + + # Check that we're computing the species fluxes correctly + for i in range(t.shape[0]): + self.assertAlmostEqual(reaction_rates[i, 0], -species_rates[i, 0], + delta=1e-6 * reaction_rates[i, 0]) + self.assertAlmostEqual(reaction_rates[i, 0], -species_rates[i, 1], + delta=1e-6 * reaction_rates[i, 0]) + self.assertAlmostEqual(reaction_rates[i, 0], species_rates[i, 2], + delta=1e-6 * reaction_rates[i, 0]) + + # Check that we've reached equilibrium by the end + self.assertAlmostEqual(reaction_rates[-1, 0], 0.0, delta=1e-2) + + # Run model with Covdep off so we can test that it is actually being implemented + rxn_system = SurfaceReactor( + T, P_initial, + n_sims=1, + initial_gas_mole_fractions={ch3: 1.0}, + initial_surface_coverages={x: 1.0}, + surface_volume_ratio=(1., 'm^-1'), + surface_site_density=(2.72e-9, 'mol/cm^2'), + termination=[]) + + rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions) + + tlist = np.logspace(-13, -5, 81, dtype=np.float64) + + # Integrate to get the solution at each time point + t = [] + y_off = [] + species_rates_off = [] + t.append(rxn_system.t) + + # You must make a copy of y because it is overwritten by DASSL at + # each call to advance() + y_off.append(rxn_system.y.copy()) + species_rates_off.append(rxn_system.core_species_rates.copy()) + for t1 in tlist: + rxn_system.advance(t1) + t.append(rxn_system.t) + # You must make a copy of y because it is overwritten by DASSL at + # each call to advance() + y_off.append(rxn_system.y.copy()) + species_rates_off.append(rxn_system.core_species_rates.copy()) + run_time = time.time() - start_time + print(f"Simulation took {run_time:.3e} seconds in {self.id()}") + + # Convert the solution vectors to np arrays + t = np.array(t, np.float64) + y_off = np.array(y_off, np.float64) + species_rates_off = np.array(species_rates_off, np.float64) + + # Check that we've reached equilibrium + self.assertAlmostEqual(species_rates_off[-1, 0], 0.0, delta=1e-2) + + # Check that coverages are different + self.assertFalse(np.allclose(y,y_off)) + self.assertFalse(np.allclose(species_rates, species_rates_off)) diff --git a/rmgpy/test_data/testing_database/kinetics/libraries/surface-example/reactions.py b/rmgpy/test_data/testing_database/kinetics/libraries/surface-example/reactions.py index e423883710..f2fe64c6fe 100644 --- a/rmgpy/test_data/testing_database/kinetics/libraries/surface-example/reactions.py +++ b/rmgpy/test_data/testing_database/kinetics/libraries/surface-example/reactions.py @@ -19,6 +19,7 @@ Ea=(1.5E3, 'J/mol'), Tmin = (200, 'K'), Tmax = (3000, 'K'), + coverage_dependence = {'OX': {'E': (0.0, 'J/mol'), 'm':0.0, 'a':0.0},} ), shortDesc = u"""Default""", longDesc = u"""Made up""", diff --git a/testing/databaseTest.py b/testing/databaseTest.py index 8ea0e899f5..8be0a6e7eb 100644 --- a/testing/databaseTest.py +++ b/testing/databaseTest.py @@ -133,6 +133,12 @@ def test_kinetics(self): self.compat_func_name = test_name yield test, family_name + test = lambda x: self.kinetics_check_coverage_dependence_units_are_correct(family_name) + test_name = "Kinetics surface family {0}: check coverage dependent units are correct?".format(family_name) + test.description = test_name + self.compat_func_name = test_name + yield test, family_name + # these families have some sort of difficulty which prevents us from testing accessibility right now difficult_families = ['Diels_alder_addition', 'Intra_R_Add_Exocyclic', 'Intra_R_Add_Endocyclic', 'Retroene'] @@ -428,6 +434,32 @@ def general_check_metal_database_has_reasonable_labels(self, library): if not entry.label[0].isupper(): raise NameError('Entry {} should start with a capital letter'.format(entry.label)) + def kinetics_check_coverage_dependence_units_are_correct(self, family_name): + """Test that each surface training reaction that has coverage dependent parameters has acceptable units""" + family = self.database.kinetics.families[family_name] + training = family.get_training_depository().entries.values() + failed = False + + for entry in training: + cov_dep = entry.data.coverage_dependence + if cov_dep: + assert isinstance(cov_dep, dict) + for species, parameters in cov_dep.items(): + assert isinstance(species, str) + assert parameters['E'] + if parameters['a'].units: + "Should be dimensionless" + failed = True + logging.error(f"Entry {entry.label} has invalid units {parameters['a'].units} for a") + if parameters['m'].units: + "Should be dimensionless" + failed = True + logging.error(f"Entry {entry.label} has invalid units {parameters['m'].units} for m") + + if failed: + raise ValueError('Surface coverage dependent parameters have incorrect units.' + 'Please check log warnings for all error messages.') + def kinetics_check_training_reactions_have_surface_attributes(self, family_name): """Test that each surface training reaction has surface attributes""" family = self.database.kinetics.families[family_name]