From e099db7da7399e313f60e2e345351cbf916f1cd7 Mon Sep 17 00:00:00 2001 From: mazeau Date: Tue, 2 Mar 2021 15:34:01 -0500 Subject: [PATCH 01/43] Adding in coverage parameters to SurfaceArrhenius Adding coverage dependent parameters to surface arrhenius kinetics added in coverage dependence in kinetics model Adding _coverage_dependence to SurfaceArrhenius and SurfaceArrhenius Cython declarations adding variable names to coverage documentation and making it more readable Fixing the updated docstrings for coverage_dependence Adding in coverage_dependence to SurfaceArrhenius and SurfaceArrheniusBEP --- rmgpy/data/kinetics/family.py | 9 +- rmgpy/kinetics/surface.pxd | 5 +- rmgpy/kinetics/surface.pyx | 196 +++++++++++++++++++++++----------- 3 files changed, 144 insertions(+), 66 deletions(-) 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/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..965037c3a5 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: + `E`, the activation energy dependence on coverage, + `m`, the power-law exponent of coverage dependence, and + `a`, the coefficient for exponential dependence on the 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): """ @@ -89,7 +95,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 +125,13 @@ 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, dict): + self._coverage_dependence = dict + cpdef double get_sticking_coefficient(self, double T) except -1: """ Return the sticking coefficient (dimensionless) at temperature `T` in K. @@ -128,6 +141,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 +241,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: + `E`, the activation energy dependence on coverage, + `m`, the power-law exponent of coverage dependence, and + `a`, the coefficient for exponential dependence on the 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): """ @@ -270,7 +290,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 +320,13 @@ 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, dict): + self._coverage_dependence = dict + cpdef double get_sticking_coefficient(self, double T, double dHrxn=0.0) except -1: """ Return the sticking coefficient (dimensionless) at @@ -309,6 +336,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 +369,7 @@ cdef class StickingCoefficientBEP(KineticsModel): T0=(1, "K"), Tmin=self.Tmin, Tmax=self.Tmax, + coverage_dependence=self.coverage_dependence, comment=self.comment, ) @@ -382,21 +411,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: + `E`, the activation energy dependence on coverage, + `m`, the power-law exponent of coverage dependence, and + `a`, the coefficient for exponential dependence on the 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 +449,13 @@ 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, dict): + self._coverage_dependence = dict + def __repr__(self): """ Return a string representation that can be used to reconstruct the @@ -426,8 +476,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 +488,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: + `E`, the activation energy dependence on coverage, + `m`, the power-law exponent of coverage dependence, and + `a`, the coefficient for exponential dependence on the 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 +531,13 @@ 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, dict): + self._coverage_dependence = dict + def __repr__(self): """ Return a string representation that can be used to reconstruct the @@ -488,7 +559,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 +577,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, ) From 21a4da7bffad800cfd26ca2d4c9c74d5c80e45fb Mon Sep 17 00:00:00 2001 From: mazeau Date: Thu, 11 Mar 2021 19:55:00 -0500 Subject: [PATCH 02/43] Use coverage parameters to modify the reaction rates using coverage dependent parameters multiply kr by cov parameters moving cov_dep and co to the init WIP: Export cantera files with coverage dependent parameters make cov dict into cantera readable format for export WIP: Adding unittests to test cantera export Make sure units of coverage dependence are si --- rmgpy/solver/surface.pyx | 44 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 42 insertions(+), 2 deletions(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index 4fcf338c17..bd82515ccc 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -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 dict cov_dep + cdef public dict cov_dep_index_smiles + cdef public dict current_surface_coverages + def __init__(self, T, P_initial, @@ -79,6 +83,9 @@ cdef class SurfaceReactor(ReactionSystem): sensitive_species=None, sensitivity_threshold=1e-3, sens_conditions=None, + cov_dep={}, + cov_dep_index_smiles={}, + current_surface_coverages={}, ): ReactionSystem.__init__(self, termination, @@ -101,6 +108,9 @@ cdef class SurfaceReactor(ReactionSystem): self.constant_volume = True self.sens_conditions = sens_conditions self.n_sims = n_sims + self.cov_dep = cov_dep + self.cov_dep_index_smiles = cov_dep_index_smiles + self.current_surface_coverages = current_surface_coverages def convert_initial_keys_to_species_objects(self, species_dict): """ @@ -161,16 +171,20 @@ cdef class SurfaceReactor(ReactionSystem): #: 1 if it's on a surface, 0 if it's in the gas phase reactions_on_surface = np.zeros((self.num_core_reactions + self.num_edge_reactions), np.int) species_on_surface = np.zeros((self.num_core_species), np.int) + cov_dep_index_smiles = {} for spec, index in self.species_index.items(): if index >= self.num_core_species: continue if spec.contains_surface_site(): species_on_surface[index] = 1 + # map for species index to SMILES for coverage dependence + cov_dep_index_smiles[index] = spec.smiles for rxn, index in self.reaction_index.items(): if rxn.is_surface_reaction(): reactions_on_surface[index] = 1 self.species_on_surface = species_on_surface self.reactions_on_surface = reactions_on_surface + self.cov_dep_index_smiles = cov_dep_index_smiles # Set initial conditions self.set_initial_conditions() @@ -206,8 +220,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, @@ -219,6 +231,10 @@ cdef class SurfaceReactor(ReactionSystem): rxn.get_surface_rate_coefficient(self.T.value_si, self.surface_site_density.value_si )) + + if rxn.kinetics.coverage_dependence: + self.cov_dep[j] = rxn.kinetics.coverage_dependence # dict + else: if not warned and rxn.kinetics.is_pressure_dependent(): logging.warning("Pressure may be varying, but using initial pressure to evaluate k(T,P) expressions!") @@ -359,6 +375,8 @@ cdef class SurfaceReactor(ReactionSystem): kf = self.kf # are already 'per m3 of reactor' even for surface reactions kr = self.kb # are already 'per m3 of reactor' even for surface reactions + cov_dep = self.cov_dep # load in coverage dependent parameters + cov_dep_index_smiles = self.cov_dep_index_smiles # load in species index to SMILES inet = self.network_indices knet = self.network_leak_coefficients @@ -384,10 +402,16 @@ 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 + current_surface_coverages = {} for j in range(num_core_species): if species_on_surface[j]: C[j] = (N[j] / V) / surface_volume_ratio_si + surface_site_fraction = N[j] / total_sites + smiles = cov_dep_index_smiles[j] + current_surface_coverages[smiles] = surface_site_fraction else: C[j] = N[j] / V #: surface species are in mol/m2, gas phase are in mol/m3 @@ -395,6 +419,14 @@ cdef class SurfaceReactor(ReactionSystem): for j in range(ir.shape[0]): k = kf[j] + if j in cov_dep: + cov_dep_j = cov_dep[j] + for smiles, parameters in cov_dep_j.items(): + if smiles in current_surface_coverages: + theta = current_surface_coverages[smiles] + if theta >= 0.1: + k *= 10. ** (parameters['a'] * theta) * theta ** parameters['m'] * np.exp(-1 * \ + parameters['E'].value_si * theta / (constants.R * self.T.value_si)) if ir[j, 0] >= num_core_species or ir[j, 1] >= num_core_species or ir[j, 2] >= num_core_species: reaction_rate = 0.0 elif ir[j, 1] == -1: # only one reactant @@ -404,6 +436,14 @@ cdef class SurfaceReactor(ReactionSystem): else: # three reactants!! (really?) reaction_rate = k * C[ir[j, 0]] * C[ir[j, 1]] * C[ir[j, 2]] k = kr[j] + if j in cov_dep: + cov_dep_j = cov_dep[j] + for smiles, parameters in cov_dep_j.items(): + if smiles in current_surface_coverages: + theta = current_surface_coverages[smiles] + if theta >= 0.1: + k *= 10. ** (parameters['a'] * theta) * theta ** parameters['m'] * np.exp(-1 * \ + parameters['E'].value_si * theta / (constants.R * self.T.value_si)) if ip[j, 0] >= num_core_species or ip[j, 1] >= num_core_species or ip[j, 2] >= num_core_species: pass elif ip[j, 1] == -1: # only one reactant From 9cf442462aa7272f08b9d2ecf8e45fdf917b5774 Mon Sep 17 00:00:00 2001 From: Emily Mazeau Date: Fri, 12 Mar 2021 14:21:42 -0500 Subject: [PATCH 03/43] Fix the coverage_dependence setter functions. Now enforces units --- rmgpy/kinetics/surface.pyx | 40 ++++++++++++++++++++++++++++++-------- 1 file changed, 32 insertions(+), 8 deletions(-) diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx index 965037c3a5..d72f71a4c5 100644 --- a/rmgpy/kinetics/surface.pyx +++ b/rmgpy/kinetics/surface.pyx @@ -129,8 +129,14 @@ cdef class StickingCoefficient(KineticsModel): """The coverage dependence parameters.""" def __get__(self): return self._coverage_dependence - def __set__(self, dict): - self._coverage_dependence = dict + 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: """ @@ -324,8 +330,14 @@ cdef class StickingCoefficientBEP(KineticsModel): """The coverage dependence parameters.""" def __get__(self): return self._coverage_dependence - def __set__(self, dict): - self._coverage_dependence = dict + 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: """ @@ -453,8 +465,14 @@ cdef class SurfaceArrhenius(Arrhenius): """The coverage dependence parameters.""" def __get__(self): return self._coverage_dependence - def __set__(self, dict): - self._coverage_dependence = dict + 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): """ @@ -535,8 +553,14 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): """The coverage dependence parameters.""" def __get__(self): return self._coverage_dependence - def __set__(self, dict): - self._coverage_dependence = dict + 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): """ From 548ce1cc0c56d26f38132930714345e2d4f61bba Mon Sep 17 00:00:00 2001 From: mazeau Date: Fri, 12 Mar 2021 14:21:42 -0500 Subject: [PATCH 04/43] Adding in coverage dep to __repr__ in surface kinetics --- rmgpy/kinetics/surface.pyx | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx index d72f71a4c5..3e4dfcf81a 100644 --- a/rmgpy/kinetics/surface.pyx +++ b/rmgpy/kinetics/surface.pyx @@ -84,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 is not None: + string += ", coverage_dependence={" + for species, parameters in self.coverage_dependence.items(): + string += f"'{species}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}" + 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) @@ -285,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 is not None: + string += ", coverage_dependence={" + for species, parameters in self.coverage_dependence.items(): + string += f"'{species}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}" + 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) @@ -482,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 is not None: + string += ", coverage_dependence={" + for species, parameters in self.coverage_dependence.items(): + string += f"'{species}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}" + 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) @@ -571,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 is not None: + string += ", coverage_dependence={" + for species, parameters in self.coverage_dependence.items(): + string += f"'{species}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}" + 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) From 3494321d0ee78ea9768b102cafc73d3903b3f719 Mon Sep 17 00:00:00 2001 From: mazeau Date: Wed, 17 Mar 2021 13:32:31 -0400 Subject: [PATCH 05/43] Adding in new unit tests to the solver surface test test that coverage dependence is doing what it should be doing Richard: Fix some surface solver tests details. This is meant to represent CH3X, not HX. Doesn't make a difference to what we're testing, but better to be right. --- rmgpy/solver/surfaceTest.py | 234 ++++++++++++++++++++++++++++++++++++ 1 file changed, 234 insertions(+) diff --git a/rmgpy/solver/surfaceTest.py b/rmgpy/solver/surfaceTest.py index b52cdb2dbc..31549b63f1 100644 --- a/rmgpy/solver/surfaceTest.py +++ b/rmgpy/solver/surfaceTest.py @@ -280,3 +280,237 @@ 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={'*': {'E': (0.0, 'J/mol'), 'm': -1.0, 'a': 0.0}})) + + core_species = [h2, x, hx] + edge_species = [] + core_reactions = [rxn1] + edge_reactions = [] + + T = 1000 + 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'), + 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, str) # species should be a string + self.assertIsInstance(parameters, dict) + self.assertIsNotNone(parameters['E']) + self.assertIsNotNone(parameters['m']) + self.assertIsNotNone(parameters['a']) + + # Integrate to get the solution at each time point + t = [] + y = [] + reaction_rates = [] + species_rates = [] + 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()) + + # 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], -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={'*': {'E': (0.0, 'J/mol'), 'm': -1.0, 'a': 0.0}}, + 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'), + 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, str) # species should be a string + self.assertIsInstance(parameters, dict) + self.assertIsNotNone(parameters['E']) + self.assertIsNotNone(parameters['m']) + self.assertIsNotNone(parameters['a']) + + # 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) + 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()) + + # 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) From 18a1810a4cae22f555805e6a5815f60ea5014b79 Mon Sep 17 00:00:00 2001 From: mazeau Date: Fri, 19 Mar 2021 18:39:50 -0400 Subject: [PATCH 06/43] Adding in coverage dependence to kinetics unittests --- rmgpy/kinetics/surfaceTest.py | 39 +++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/rmgpy/kinetics/surfaceTest.py b/rmgpy/kinetics/surfaceTest.py index 2cab8fec9d..d4031aeb94 100644 --- a/rmgpy/kinetics/surfaceTest.py +++ b/rmgpy/kinetics/surfaceTest.py @@ -55,6 +55,7 @@ def setUp(self): self.T0 = 1. self.Tmin = 300. self.Tmax = 3000. + self.coverage_dependence = {'*': {'E': (0.0, 'J/mol'), 'm': -1.0, 'a': 0.0}} self.comment = 'O2 dissociative' self.stick = StickingCoefficient( A=self.A, @@ -64,6 +65,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 +110,12 @@ 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. + """ + self.assertEqual(self.stick.coverage_dependence, self.coverage_dependence) + def test_is_temperature_valid(self): """ Test the StickingCoefficient.is_temperature_valid() method. @@ -137,6 +145,10 @@ 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) + self.assertEqual(self.stick.coverage_dependence.keys(), 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][key].value_si) self.assertEqual(dir(self.stick), dir(stick)) def test_repr(self): @@ -160,6 +172,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(self.stick.coverage_dependence.keys(), 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][key].value_si) self.assertEqual(dir(self.stick), dir(stick)) def test_copy(self): @@ -181,6 +197,10 @@ 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) + self.assertEqual(self.stick.coverage_dependence.keys(), 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][key].value_si) self.assertEqual(dir(self.stick), dir(stick)) def test_is_identical_to(self): @@ -207,6 +227,7 @@ def setUp(self): self.T0 = 1. self.Tmin = 300. self.Tmax = 3000. + self.coverage_dependence = {'*': {'E': (0.0, 'J/mol'), 'm': -1.0, 'a': 0.0}} self.comment = 'CH3x + Hx <=> CH4 + x + x' self.surfarr = SurfaceArrhenius( A=(self.A, 'm^2/(mol*s)'), @@ -216,6 +237,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 +282,12 @@ 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. + """ + self.assertEqual(self.surfarr.coverage_dependence, self.coverage_dependence) + def test_is_temperature_valid(self): """ Test the SurfaceArrhenius.is_temperature_valid() method. @@ -289,6 +317,9 @@ 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 species, parameters in self.surfarr.coverage_dependence.items(): + for key, value in parameters.items(): + self.assertEqual(value.value_si, surfarr.coverage_dependence[species][key].value_si) self.assertEqual(dir(self.surfarr), dir(surfarr)) def test_repr(self): @@ -312,6 +343,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(self.surfarr.coverage_dependence.keys(), 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][key].value_si) self.assertEqual(dir(self.surfarr), dir(surfarr)) def test_copy(self): @@ -333,6 +368,10 @@ 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) + self.assertEqual(self.surfarr.coverage_dependence.keys(), 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][key].value_si) self.assertEqual(dir(self.surfarr), dir(surfarr)) def test_is_identical_to(self): From bb2c693cd3f77d53c2cac6d1fd680598541f8224 Mon Sep 17 00:00:00 2001 From: mazeau Date: Fri, 19 Mar 2021 19:13:48 -0400 Subject: [PATCH 07/43] Adding in coverage dependent database tests --- testing/databaseTest.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) 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] From b0be206758f45b8a85db23a8abf95e77f081abbb Mon Sep 17 00:00:00 2001 From: Chris B <56306881+ChrisBNEU@users.noreply.github.com> Date: Sun, 21 Mar 2021 20:50:51 -0400 Subject: [PATCH 08/43] Add species labels in reactionTest --- rmgpy/reactionTest.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/rmgpy/reactionTest.py b/rmgpy/reactionTest.py index f0e7aaea6e..ed7cb3b5ad 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], From a16a2be35137c4f933f94111822ef2a0ed67944e Mon Sep 17 00:00:00 2001 From: Emily Mazeau Date: Sat, 26 Jun 2021 00:46:50 -0400 Subject: [PATCH 09/43] adding a coverage-dependent reaction to reactionTest --- rmgpy/reactionTest.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/rmgpy/reactionTest.py b/rmgpy/reactionTest.py index ed7cb3b5ad..64e8eb002a 100644 --- a/rmgpy/reactionTest.py +++ b/rmgpy/reactionTest.py @@ -231,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()) From 35c7b0bcafa5a37bed39e55b96f519ed7bb7aa57 Mon Sep 17 00:00:00 2001 From: mazeau Date: Mon, 5 Apr 2021 13:38:19 -0400 Subject: [PATCH 10/43] Changing coverage dependence SMILES to species objects Look for coverage_dependence attribute --- rmgpy/data/kinetics/library.py | 9 +++++++++ rmgpy/kinetics/surface.pyx | 8 ++++---- rmgpy/kinetics/surfaceTest.py | 5 +++-- rmgpy/rmg/model.py | 6 ++++++ rmgpy/solver/surface.pyx | 32 ++++++++++++++++---------------- 5 files changed, 38 insertions(+), 22 deletions(-) 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.pyx b/rmgpy/kinetics/surface.pyx index 3e4dfcf81a..2c5156f957 100644 --- a/rmgpy/kinetics/surface.pyx +++ b/rmgpy/kinetics/surface.pyx @@ -87,7 +87,7 @@ cdef class StickingCoefficient(KineticsModel): if self.coverage_dependence is not None: string += ", coverage_dependence={" for species, parameters in self.coverage_dependence.items(): - string += f"'{species}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}" + string += f"'{species!r}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," 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) @@ -293,7 +293,7 @@ cdef class StickingCoefficientBEP(KineticsModel): if self.coverage_dependence is not None: string += ", coverage_dependence={" for species, parameters in self.coverage_dependence.items(): - string += f"'{species}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}" + string += f"'{species!r}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," 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) @@ -495,7 +495,7 @@ cdef class SurfaceArrhenius(Arrhenius): if self.coverage_dependence is not None: string += ", coverage_dependence={" for species, parameters in self.coverage_dependence.items(): - string += f"'{species}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}" + string += f"'{species!r}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," 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) @@ -589,7 +589,7 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): if self.coverage_dependence is not None: string += ", coverage_dependence={" for species, parameters in self.coverage_dependence.items(): - string += f"'{species}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}" + string += f"'{species!r}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," 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) diff --git a/rmgpy/kinetics/surfaceTest.py b/rmgpy/kinetics/surfaceTest.py index d4031aeb94..27d8232ba3 100644 --- a/rmgpy/kinetics/surfaceTest.py +++ b/rmgpy/kinetics/surfaceTest.py @@ -36,6 +36,7 @@ import numpy as np from rmgpy.kinetics.surface import StickingCoefficient, SurfaceArrhenius +from rmgpy.species import Species ################################################################################ @@ -55,7 +56,7 @@ def setUp(self): self.T0 = 1. self.Tmin = 300. self.Tmax = 3000. - self.coverage_dependence = {'*': {'E': (0.0, 'J/mol'), 'm': -1.0, 'a': 0.0}} + self.coverage_dependence = {Species().from_adjacency_list('1 X u0 p0 c0'): {'E': (0.0, 'J/mol'), 'm': -1.0, 'a': 0.0}} self.comment = 'O2 dissociative' self.stick = StickingCoefficient( A=self.A, @@ -227,7 +228,7 @@ def setUp(self): self.T0 = 1. self.Tmin = 300. self.Tmax = 3000. - self.coverage_dependence = {'*': {'E': (0.0, 'J/mol'), 'm': -1.0, 'a': 0.0}} + self.coverage_dependence = {Species().from_adjacency_list('1 X u0 p0 c0'): {'E': (0.0, 'J/mol'), 'm': -1.0, 'a': 0.0}} self.comment = 'CH3x + Hx <=> CH4 + x + x' self.surfarr = SurfaceArrhenius( A=(self.A, 'm^2/(mol*s)'), diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index b2b3e608c1..835440d4d4 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -887,6 +887,12 @@ def apply_kinetics_to_reaction(self, reaction): reaction.reverse = None reaction.kinetics = kinetics + if hasattr(kinetics, 'coverage_dependence'): + if kinetics.coverage_dependence: + for species, values in kinetics.coverage_dependence.items(): + species_in_model = self.make_new_species(species) + kinetics.coverage_dependence[species_in_model] = values + def generate_kinetics(self, reaction): """ Generate best possible kinetics for the given `reaction` using the kinetics database. diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index bd82515ccc..b1aaa04a58 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -68,7 +68,7 @@ cdef class SurfaceReactor(ReactionSystem): cdef public np.ndarray species_on_surface # (catalyst surface, not core/edge surface) cdef public dict cov_dep - cdef public dict cov_dep_index_smiles + cdef public dict cov_dep_index_species cdef public dict current_surface_coverages def __init__(self, @@ -84,7 +84,7 @@ cdef class SurfaceReactor(ReactionSystem): sensitivity_threshold=1e-3, sens_conditions=None, cov_dep={}, - cov_dep_index_smiles={}, + cov_dep_index_species={}, current_surface_coverages={}, ): ReactionSystem.__init__(self, @@ -109,7 +109,7 @@ cdef class SurfaceReactor(ReactionSystem): self.sens_conditions = sens_conditions self.n_sims = n_sims self.cov_dep = cov_dep - self.cov_dep_index_smiles = cov_dep_index_smiles + self.cov_dep_index_species = cov_dep_index_species self.current_surface_coverages = current_surface_coverages def convert_initial_keys_to_species_objects(self, species_dict): @@ -171,20 +171,20 @@ cdef class SurfaceReactor(ReactionSystem): #: 1 if it's on a surface, 0 if it's in the gas phase reactions_on_surface = np.zeros((self.num_core_reactions + self.num_edge_reactions), np.int) species_on_surface = np.zeros((self.num_core_species), np.int) - cov_dep_index_smiles = {} + cov_dep_index_species = {} for spec, index in self.species_index.items(): if index >= self.num_core_species: continue if spec.contains_surface_site(): species_on_surface[index] = 1 - # map for species index to SMILES for coverage dependence - cov_dep_index_smiles[index] = spec.smiles + # map for species index to species for coverage dependence + cov_dep_index_species[index] = spec for rxn, index in self.reaction_index.items(): if rxn.is_surface_reaction(): reactions_on_surface[index] = 1 self.species_on_surface = species_on_surface self.reactions_on_surface = reactions_on_surface - self.cov_dep_index_smiles = cov_dep_index_smiles + self.cov_dep_index_species = cov_dep_index_species # Set initial conditions self.set_initial_conditions() @@ -376,7 +376,7 @@ cdef class SurfaceReactor(ReactionSystem): kf = self.kf # are already 'per m3 of reactor' even for surface reactions kr = self.kb # are already 'per m3 of reactor' even for surface reactions cov_dep = self.cov_dep # load in coverage dependent parameters - cov_dep_index_smiles = self.cov_dep_index_smiles # load in species index to SMILES + cov_dep_index_species = self.cov_dep_index_species # load in species index to species inet = self.network_indices knet = self.network_leak_coefficients @@ -410,8 +410,8 @@ cdef class SurfaceReactor(ReactionSystem): if species_on_surface[j]: C[j] = (N[j] / V) / surface_volume_ratio_si surface_site_fraction = N[j] / total_sites - smiles = cov_dep_index_smiles[j] - current_surface_coverages[smiles] = surface_site_fraction + species = cov_dep_index_species[j] + current_surface_coverages[species] = surface_site_fraction else: C[j] = N[j] / V #: surface species are in mol/m2, gas phase are in mol/m3 @@ -421,9 +421,9 @@ cdef class SurfaceReactor(ReactionSystem): k = kf[j] if j in cov_dep: cov_dep_j = cov_dep[j] - for smiles, parameters in cov_dep_j.items(): - if smiles in current_surface_coverages: - theta = current_surface_coverages[smiles] + for species, parameters in cov_dep_j.items(): + if species in current_surface_coverages: + theta = current_surface_coverages[species] if theta >= 0.1: k *= 10. ** (parameters['a'] * theta) * theta ** parameters['m'] * np.exp(-1 * \ parameters['E'].value_si * theta / (constants.R * self.T.value_si)) @@ -438,9 +438,9 @@ cdef class SurfaceReactor(ReactionSystem): k = kr[j] if j in cov_dep: cov_dep_j = cov_dep[j] - for smiles, parameters in cov_dep_j.items(): - if smiles in current_surface_coverages: - theta = current_surface_coverages[smiles] + for species, parameters in cov_dep_j.items(): + if species in current_surface_coverages: + theta = current_surface_coverages[species] if theta >= 0.1: k *= 10. ** (parameters['a'] * theta) * theta ** parameters['m'] * np.exp(-1 * \ parameters['E'].value_si * theta / (constants.R * self.T.value_si)) From a1982430188a0f0aecfbde6596f1daa357003d82 Mon Sep 17 00:00:00 2001 From: mazeau Date: Sun, 2 May 2021 12:31:47 -0400 Subject: [PATCH 11/43] surfaceTest can check for coverage dependence paramaters --- rmgpy/kinetics/surfaceTest.py | 74 +++++++++++++++++++++++++++-------- 1 file changed, 58 insertions(+), 16 deletions(-) diff --git a/rmgpy/kinetics/surfaceTest.py b/rmgpy/kinetics/surfaceTest.py index 27d8232ba3..d402147931 100644 --- a/rmgpy/kinetics/surfaceTest.py +++ b/rmgpy/kinetics/surfaceTest.py @@ -37,6 +37,7 @@ from rmgpy.kinetics.surface import StickingCoefficient, SurfaceArrhenius from rmgpy.species import Species +from rmgpy.molecule import Molecule ################################################################################ @@ -146,10 +147,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) - self.assertEqual(self.stick.coverage_dependence.keys(), stick.coverage_dependence.keys()) + 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(): - for key, value in parameters.items(): - self.assertEqual(value.value_si, stick.coverage_dependence[species][key].value_si) + 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): @@ -158,7 +169,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) @@ -173,10 +184,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(self.stick.coverage_dependence.keys(), stick.coverage_dependence.keys()) + 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][key].value_si) + self.assertEqual(value.value_si, stick.coverage_dependence[species.label][key].value_si) self.assertEqual(dir(self.stick), dir(stick)) def test_copy(self): @@ -198,10 +209,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) - self.assertEqual(self.stick.coverage_dependence.keys(), stick.coverage_dependence.keys()) + 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(): - for key, value in parameters.items(): - self.assertEqual(value.value_si, stick.coverage_dependence[species][key].value_si) + 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): @@ -318,9 +339,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(): - for key, value in parameters.items(): - self.assertEqual(value.value_si, surfarr.coverage_dependence[species][key].value_si) + 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): @@ -344,10 +376,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(self.surfarr.coverage_dependence.keys(), surfarr.coverage_dependence.keys()) + 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][key].value_si) + self.assertEqual(value.value_si, surfarr.coverage_dependence[species.label][key].value_si) self.assertEqual(dir(self.surfarr), dir(surfarr)) def test_copy(self): @@ -369,10 +401,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) - self.assertEqual(self.surfarr.coverage_dependence.keys(), surfarr.coverage_dependence.keys()) + 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(): - for key, value in parameters.items(): - self.assertEqual(value.value_si, surfarr.coverage_dependence[species][key].value_si) + 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): From 2fcea33d4e772b41a5c4347075e10c92f13248a8 Mon Sep 17 00:00:00 2001 From: mazeau Date: Thu, 6 May 2021 18:42:06 -0400 Subject: [PATCH 12/43] adding in coverage dependence to one of the testing examples --- .../kinetics/libraries/surface-example/reactions.py | 1 + 1 file changed, 1 insertion(+) 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""", From f62cedf7a3035097acbb1acb6086e3a8e7ea681e Mon Sep 17 00:00:00 2001 From: mazeau Date: Fri, 7 May 2021 16:40:03 -0400 Subject: [PATCH 13/43] use species labels in the __repr__ --- rmgpy/kinetics/surface.pyx | 8 +++---- rmgpy/kinetics/surfaceTest.py | 41 +++++++++++++++++++++++++++++++---- 2 files changed, 41 insertions(+), 8 deletions(-) diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx index 2c5156f957..542cdc99c8 100644 --- a/rmgpy/kinetics/surface.pyx +++ b/rmgpy/kinetics/surface.pyx @@ -87,7 +87,7 @@ cdef class StickingCoefficient(KineticsModel): if self.coverage_dependence is not None: string += ", coverage_dependence={" for species, parameters in self.coverage_dependence.items(): - string += f"'{species!r}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," + string += f"{species.to_chemkin()!r}: {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," 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) @@ -293,7 +293,7 @@ cdef class StickingCoefficientBEP(KineticsModel): if self.coverage_dependence is not None: string += ", coverage_dependence={" for species, parameters in self.coverage_dependence.items(): - string += f"'{species!r}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," + string += f"{species.to_chemkin()!r}: {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," 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) @@ -495,7 +495,7 @@ cdef class SurfaceArrhenius(Arrhenius): if self.coverage_dependence is not None: string += ", coverage_dependence={" for species, parameters in self.coverage_dependence.items(): - string += f"'{species!r}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," + string += f"{species.to_chemkin()!r}: {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," 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) @@ -589,7 +589,7 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): if self.coverage_dependence is not None: string += ", coverage_dependence={" for species, parameters in self.coverage_dependence.items(): - string += f"'{species!r}': {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," + string += f"{species.to_chemkin()!r}: {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," 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) diff --git a/rmgpy/kinetics/surfaceTest.py b/rmgpy/kinetics/surfaceTest.py index d402147931..8625600d47 100644 --- a/rmgpy/kinetics/surfaceTest.py +++ b/rmgpy/kinetics/surfaceTest.py @@ -38,6 +38,7 @@ from rmgpy.kinetics.surface import StickingCoefficient, SurfaceArrhenius from rmgpy.species import Species from rmgpy.molecule import Molecule +import rmgpy.quantity as quantity ################################################################################ @@ -57,7 +58,10 @@ def setUp(self): self.T0 = 1. self.Tmin = 300. self.Tmax = 3000. - self.coverage_dependence = {Species().from_adjacency_list('1 X u0 p0 c0'): {'E': (0.0, 'J/mol'), 'm': -1.0, 'a': 0.0}} + s = Species().from_adjacency_list('1 X u0 p0 c0') + s.label = 'X' + self.coverage_dependence = {s: {'E': quantity.Energy(0.0, 'J/mol'), 'm': quantity.Dimensionless(-1.0), + 'a': quantity.Dimensionless(0.0)}} self.comment = 'O2 dissociative' self.stick = StickingCoefficient( A=self.A, @@ -116,7 +120,20 @@ def test_coverage_dependence(self): """ Test that the coverage dependent parameters was properly set. """ - self.assertEqual(self.stick.coverage_dependence, self.coverage_dependence) + 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): """ @@ -249,7 +266,10 @@ def setUp(self): self.T0 = 1. self.Tmin = 300. self.Tmax = 3000. - self.coverage_dependence = {Species().from_adjacency_list('1 X u0 p0 c0'): {'E': (0.0, 'J/mol'), 'm': -1.0, 'a': 0.0}} + s = Species().from_adjacency_list('1 X u0 p0 c0') + s.label = 'X' + self.coverage_dependence = {s: {'E': quantity.Energy(0.0, 'J/mol'), 'm': quantity.Dimensionless(-1.0), + 'a': quantity.Dimensionless(0.0)}} self.comment = 'CH3x + Hx <=> CH4 + x + x' self.surfarr = SurfaceArrhenius( A=(self.A, 'm^2/(mol*s)'), @@ -308,7 +328,20 @@ def test_coverage_dependence(self): """ Test that the coverage dependent parameters was properly set. """ - self.assertEqual(self.surfarr.coverage_dependence, self.coverage_dependence) + 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): """ From b683d1b323fb255aa8c22716cac7d80564db1251 Mon Sep 17 00:00:00 2001 From: mazeau Date: Wed, 3 Mar 2021 15:26:47 -0500 Subject: [PATCH 14/43] Writing coverage dependencies in Chemkin output --- rmgpy/chemkin.pyx | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index 0f0f3e55e8..3d03de8c6f 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -1789,6 +1789,10 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals kinetics.Ea.value_si / 4184. ) string += '\n STICK' + if kinetics.coverage_dependence is not None: + for species, cov_params in kinetics.coverage_dependence.items(): + string += f'\n COV / {species.label:<41}' + string += f"{cov_params['E'].value_si:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['a'].value:<6.3f} /" elif isinstance(kinetics, _kinetics.Arrhenius): conversion_factor = kinetics.A.get_conversion_factor_from_si_to_cm_mol_s() if not isinstance(kinetics, _kinetics.SurfaceArrhenius): @@ -1805,6 +1809,10 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals kinetics.n.value_si, kinetics.Ea.value_si / 4184. ) + if isinstance(kinetics, _kinetics.SurfaceArrhenius) and kinetics.coverage_dependence: + for species, cov_params in kinetics.coverage_dependence.items(): + string += f'\n COV / {species.label:<41}' + string += f"{cov_params['E'].value_si:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['a'].value:<6.3f} /" elif isinstance(kinetics, (_kinetics.Lindemann, _kinetics.Troe)): arrhenius = kinetics.arrheniusHigh conversion_factor = arrhenius.A.get_conversion_factor_from_si_to_cm_mol_s() From 268a3e6af54fd57b6e18601a9013b578ebfb3119 Mon Sep 17 00:00:00 2001 From: Richard West Date: Sat, 26 Jun 2021 00:36:10 -0400 Subject: [PATCH 15/43] Can read Coverage Dependent kinetics from Chemkin files For now, only reactions with COV are turned into SurfaceArrhenius objects, all others are left as Arrhenius. I'm not quite sure what's happening around lines 467-473 (Emily's changes, I'm just committing them) --- rmgpy/chemkin.pyx | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index 3d03de8c6f..defd44f4a2 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -207,6 +207,7 @@ def read_kinetics_entry(entry, species_dict, Aunits, Eunits): kinetics.update({ 'chebyshev coefficients': [], 'efficiencies': {}, + 'coverage dependence': {} }) # Note that the subsequent lines could be in any order @@ -288,6 +289,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 +452,25 @@ def _read_kinetics_line(line, reaction, species_dict, Eunits, kunits, klow_units # Duplicate reaction reaction.duplicate = True + elif 'COV' in line: + k = kinetics['arrhenius high'] + kinetics['surface arrhenius'] = _kinetics.SurfaceArrhenius( + A=(k.A.value, kunits), + n=k.n, + Ea=k.Ea, + T0=k.T0, + ) + del kinetics['arrhenius high'] + + tokens = tokens[1].split() + kinetics['coverage dependence'][species_dict[tokens[0].strip()]] = {'E': (tokens[1].strip(), Eunits), 'm': tokens[2].strip(), 'a':tokens[3].strip()} + try: + # is a sticking coefficient + kinetics['sticking coefficient'][species_dict[tokens[0].strip()]] = {'E': (tokens[1].strip(), Eunits), 'm': tokens[2].strip(), 'a':tokens[3].strip()} + except KeyError: + # is a not a sticking coefficient + kinetics['surface arrhenius'].coverage_dependence[species_dict[tokens[0].strip()]] = {'E': (tokens[1].strip(), 'J/mol'), 'm': tokens[2].strip(), 'a':tokens[3].strip()} + elif 'LOW' in line: # Low-pressure-limit Arrhenius parameters tokens = tokens[1].split() From 9cdd49b31525bd45665816eea9ed9075de963d51 Mon Sep 17 00:00:00 2001 From: mazeau Date: Tue, 1 Jun 2021 16:16:49 -0400 Subject: [PATCH 16/43] Turn coverage dependence on/off through the catalyst block in the input file --- rmgpy/rmg/input.py | 12 ++++++++++-- rmgpy/rmg/main.py | 2 ++ rmgpy/rmg/model.py | 23 ++++++++++++++++++---- rmgpy/solver/surface.pyx | 41 ++++++++++++++++++++++++---------------- 4 files changed, 56 insertions(+), 22 deletions(-) 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 835440d4d4..ebdf5ed8c3 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -652,6 +652,9 @@ def enlarge(self, new_object=None, react_edge=False, # assume the kinetics are satisfactory if reaction.kinetics is None: self.apply_kinetics_to_reaction(reaction) + elif hasattr(reaction.kinetics, 'coverage_dependence'): + if reaction.kinetics.coverage_dependence: + self.apply_coverage_dependence_to_reaction(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. @@ -865,6 +868,14 @@ def generate_thermo(self, spc, rename=False): spc.generate_energy_transfer_model() + def apply_coverage_dependence_to_reaction(self, kinetics): + """Apply the coverage dependence kinetics""" + 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 @@ -888,10 +899,8 @@ def apply_kinetics_to_reaction(self, reaction): reaction.kinetics = kinetics if hasattr(kinetics, 'coverage_dependence'): - if kinetics.coverage_dependence: - for species, values in kinetics.coverage_dependence.items(): - species_in_model = self.make_new_species(species) - kinetics.coverage_dependence[species_in_model] = values + if reaction.kinetics.coverage_dependence: + self.apply_coverage_dependence_to_reaction(kinetics) def generate_kinetics(self, reaction): """ @@ -1502,6 +1511,9 @@ 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 hasattr(r.kinetics, 'coverage_dependence'): + if r.kinetics.coverage_dependence: + self.apply_coverage_dependence_to_reaction(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() \ @@ -1601,6 +1613,9 @@ 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 hasattr(r.kinetics, 'coverage_dependence'): + if r.kinetics.coverage_dependence: + self.apply_coverage_dependence_to_reaction(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 b1aaa04a58..495d4c6d18 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -67,6 +67,7 @@ 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 cov_dep cdef public dict cov_dep_index_species cdef public dict current_surface_coverages @@ -83,6 +84,7 @@ cdef class SurfaceReactor(ReactionSystem): sensitive_species=None, sensitivity_threshold=1e-3, sens_conditions=None, + coverage_dependence=False, cov_dep={}, cov_dep_index_species={}, current_surface_coverages={}, @@ -104,6 +106,7 @@ 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 @@ -312,6 +315,7 @@ cdef class SurfaceReactor(ReactionSystem): surface_volume_ratio_si = self.surface_volume_ratio.value_si # 1/m total_surface_sites = V * surface_volume_ratio_si * self.surface_site_density.value_si # total surface sites in reactor + coverage_dependence = self.coverage_dependence for spec, coverage in self.initial_surface_coverages.items(): i = self.get_species_index(spec) @@ -368,6 +372,8 @@ 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 bint coverage_dependence + cdef dict cov_dep, cov_dep_index_species ir = self.reactant_indices ip = self.product_indices @@ -375,6 +381,7 @@ cdef class SurfaceReactor(ReactionSystem): kf = self.kf # are already 'per m3 of reactor' even for surface reactions kr = self.kb # are already 'per m3 of reactor' even for surface reactions + coverage_dependence = self.coverage_dependence # turn coverage dependence on/off cov_dep = self.cov_dep # load in coverage dependent parameters cov_dep_index_species = self.cov_dep_index_species # load in species index to species @@ -419,14 +426,15 @@ cdef class SurfaceReactor(ReactionSystem): for j in range(ir.shape[0]): k = kf[j] - if j in cov_dep: - cov_dep_j = cov_dep[j] - for species, parameters in cov_dep_j.items(): - if species in current_surface_coverages: - theta = current_surface_coverages[species] - if theta >= 0.1: - k *= 10. ** (parameters['a'] * theta) * theta ** parameters['m'] * np.exp(-1 * \ - parameters['E'].value_si * theta / (constants.R * self.T.value_si)) + if coverage_dependence: + if j in cov_dep: + cov_dep_j = cov_dep[j] + for species, parameters in cov_dep_j.items(): + if species in current_surface_coverages: + theta = current_surface_coverages[species] + if theta >= 1e-15: + k *= 10. ** (parameters['a'].value_si * theta) * theta ** parameters['m'].value_si *\ + np.exp(-1 * parameters['E'].value_si * theta / (constants.R * self.T.value_si)) if ir[j, 0] >= num_core_species or ir[j, 1] >= num_core_species or ir[j, 2] >= num_core_species: reaction_rate = 0.0 elif ir[j, 1] == -1: # only one reactant @@ -436,14 +444,15 @@ cdef class SurfaceReactor(ReactionSystem): else: # three reactants!! (really?) reaction_rate = k * C[ir[j, 0]] * C[ir[j, 1]] * C[ir[j, 2]] k = kr[j] - if j in cov_dep: - cov_dep_j = cov_dep[j] - for species, parameters in cov_dep_j.items(): - if species in current_surface_coverages: - theta = current_surface_coverages[species] - if theta >= 0.1: - k *= 10. ** (parameters['a'] * theta) * theta ** parameters['m'] * np.exp(-1 * \ - parameters['E'].value_si * theta / (constants.R * self.T.value_si)) + if coverage_dependence: + if j in cov_dep: + cov_dep_j = cov_dep[j] + for species, parameters in cov_dep_j.items(): + if species in current_surface_coverages: + theta = current_surface_coverages[species] + if theta >= 1e-15: + k *= 10. ** (parameters['a'].value_si * theta) * theta ** parameters['m'].value_si *\ + np.exp(-1 * parameters['E'].value_si * theta / (constants.R * self.T.value_si)) if ip[j, 0] >= num_core_species or ip[j, 1] >= num_core_species or ip[j, 2] >= num_core_species: pass elif ip[j, 1] == -1: # only one reactant From 127d7257c0250385433e5e38c6e548bdae531a8b Mon Sep 17 00:00:00 2001 From: mazeau Date: Fri, 4 Jun 2021 17:36:28 -0400 Subject: [PATCH 17/43] Do not write coverage dependent parameters if coverage dependence is turned off save_chemkin_file has coverage_dependence flag, which is passed on to write_kinetics_entry --- rmgpy/chemkin.pyx | 39 ++++++++++++++++++++++++--------------- 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index defd44f4a2..e86e95067b 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -1710,7 +1710,8 @@ def write_reaction_string(reaction, java_library=False): ################################################################################ -def write_kinetics_entry(reaction, species_list, verbose=True, java_library=False, commented=False): +def write_kinetics_entry(reaction, species_list, verbose=True, java_library=False, commented=False, + coverage_dependence=False): """ Return a string representation of the reaction as used in a Chemkin file. Use `verbose = True` to turn on kinetics comments. @@ -1742,7 +1743,8 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals specific_collider=reaction.specific_collider, reversible=reaction.reversible, kinetics=kinetics) - string += write_kinetics_entry(new_reaction, species_list, verbose, java_library, commented) + string += write_kinetics_entry(new_reaction, species_list, verbose, java_library, commented, + coverage_dependence=coverage_dependence) string += "DUPLICATE\n" if commented: @@ -1811,10 +1813,12 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals kinetics.Ea.value_si / 4184. ) string += '\n STICK' - if kinetics.coverage_dependence is not None: - for species, cov_params in kinetics.coverage_dependence.items(): - string += f'\n COV / {species.label:<41}' - string += f"{cov_params['E'].value_si:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['a'].value:<6.3f} /" + if coverage_dependence: + if kinetics.coverage_dependence is not None: + for species, cov_params in kinetics.coverage_dependence.items(): + label = get_species_identifier(species) + string += f'\n COV / {label:<41}' + string += f"{cov_params['E'].value_si:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['a'].value:<6.3f} /" elif isinstance(kinetics, _kinetics.Arrhenius): conversion_factor = kinetics.A.get_conversion_factor_from_si_to_cm_mol_s() if not isinstance(kinetics, _kinetics.SurfaceArrhenius): @@ -1832,9 +1836,11 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals kinetics.Ea.value_si / 4184. ) if isinstance(kinetics, _kinetics.SurfaceArrhenius) and kinetics.coverage_dependence: - for species, cov_params in kinetics.coverage_dependence.items(): - string += f'\n COV / {species.label:<41}' - string += f"{cov_params['E'].value_si:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['a'].value:<6.3f} /" + if coverage_dependence is True: + for species, cov_params in kinetics.coverage_dependence.items(): + label = get_species_identifier(species) + string += f'\n COV / {label:<41}' + string += f"{cov_params['E'].value_si:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['a'].value:<6.3f} /" elif isinstance(kinetics, (_kinetics.Lindemann, _kinetics.Troe)): arrhenius = kinetics.arrheniusHigh conversion_factor = arrhenius.A.get_conversion_factor_from_si_to_cm_mol_s() @@ -2101,7 +2107,7 @@ def save_transport_file(path, species): )) -def save_chemkin_file(path, species, reactions, verbose=True, check_for_duplicates=True): +def save_chemkin_file(path, species, reactions, verbose=True, check_for_duplicates=True, coverage_dependence=False): """ Save a Chemkin input file to `path` on disk containing the provided lists of `species` and `reactions`. @@ -2149,7 +2155,7 @@ def save_chemkin_file(path, species, reactions, verbose=True, check_for_duplicat global _chemkin_reaction_count _chemkin_reaction_count = 0 for rxn in reactions: - f.write(write_kinetics_entry(rxn, species_list=species, verbose=verbose)) + f.write(write_kinetics_entry(rxn, species_list=species, verbose=verbose, coverage_dependence=coverage_dependence)) # Don't forget to mark duplicates! f.write('\n') f.write('END\n\n') @@ -2159,7 +2165,7 @@ def save_chemkin_file(path, species, reactions, verbose=True, check_for_duplicat def save_chemkin_surface_file(path, species, reactions, verbose=True, check_for_duplicates=True, - surface_site_density=None): + surface_site_density=None, coverage_dependence=False): """ Save a Chemkin *surface* input file to `path` on disk containing the provided lists of `species` and `reactions`. @@ -2209,7 +2215,8 @@ def save_chemkin_surface_file(path, species, reactions, verbose=True, check_for_ global _chemkin_reaction_count _chemkin_reaction_count = 0 for rxn in reactions: - f.write(write_kinetics_entry(rxn, species_list=species, verbose=verbose)) + f.write(write_kinetics_entry(rxn, species_list=species, verbose=verbose, + coverage_dependence=coverage_dependence)) f.write('\n') f.write('END\n\n') f.close() @@ -2302,11 +2309,13 @@ def save_chemkin(reaction_model, path, verbose_path, dictionary_path=None, trans # We should already have marked everything as duplicates by now so use check_for_duplicates=False save_chemkin_file(gas_path, gas_species_list, gas_rxn_list, verbose=False, check_for_duplicates=False) save_chemkin_surface_file(surface_path, surface_species_list, surface_rxn_list, verbose=False, - check_for_duplicates=False, surface_site_density=reaction_model.surface_site_density) + check_for_duplicates=False, surface_site_density=reaction_model.surface_site_density, + coverage_dependence=reaction_model.coverage_dependence) logging.info('Saving annotated version of Chemkin files...') save_chemkin_file(gas_verbose_path, gas_species_list, gas_rxn_list, verbose=True, check_for_duplicates=False) save_chemkin_surface_file(surface_verbose_path, surface_species_list, surface_rxn_list, verbose=True, - check_for_duplicates=False, surface_site_density=reaction_model.surface_site_density) + check_for_duplicates=False, surface_site_density=reaction_model.surface_site_density, + coverage_dependence=reaction_model.coverage_dependence) else: # Gas phase only From 1c2d51d37995d3a4f39f4417a441b18621e98811 Mon Sep 17 00:00:00 2001 From: mazeau Date: Wed, 12 May 2021 16:20:14 -0400 Subject: [PATCH 18/43] Adding in chemkin test to check we can read and write coverage dependent parameters --- rmgpy/chemkinTest.py | 77 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 77 insertions(+) diff --git a/rmgpy/chemkinTest.py b/rmgpy/chemkinTest.py index 271df8d0f3..702dd52a08 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,82 @@ 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: {'E': (0.1, 'J/mol'), 'm': -1.0, 'a': 1.0},})) + + 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, + coverage_dependence=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): From 884c4e58bf36eb811afff846e9194a2a2ef8d2f2 Mon Sep 17 00:00:00 2001 From: Richard West Date: Sun, 27 Jun 2021 23:13:35 -0400 Subject: [PATCH 19/43] Turn on coverage dependence in the coverage-dependence solver tests. This defaults to False if not specified as true. The fact that this wasn't noticed, tells us that the test isn't really testing that we do evaluate coverage-dependence correctly. --- rmgpy/solver/surfaceTest.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/rmgpy/solver/surfaceTest.py b/rmgpy/solver/surfaceTest.py index 31549b63f1..121c0be0db 100644 --- a/rmgpy/solver/surfaceTest.py +++ b/rmgpy/solver/surfaceTest.py @@ -331,6 +331,7 @@ def test_solve_h2_coverage_dependence(self): 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) @@ -452,6 +453,7 @@ def test_solve_ch3_coverage_dependence(self): 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. From 93979ea7c9960e18c6950a07c25ea201f0164e55 Mon Sep 17 00:00:00 2001 From: Emily Mazeau Date: Tue, 29 Jun 2021 22:50:11 -0400 Subject: [PATCH 20/43] Correcting ordering of COV parameters in Chemkin files. Should be (a, m, E). --- rmgpy/chemkin.pyx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index e86e95067b..e68167247e 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -463,13 +463,13 @@ def _read_kinetics_line(line, reaction, species_dict, Eunits, kunits, klow_units del kinetics['arrhenius high'] tokens = tokens[1].split() - kinetics['coverage dependence'][species_dict[tokens[0].strip()]] = {'E': (tokens[1].strip(), Eunits), 'm': tokens[2].strip(), 'a':tokens[3].strip()} + kinetics['coverage dependence'][species_dict[tokens[0].strip()]] = {'E': (tokens[3].strip(), Eunits), 'm': tokens[2].strip(), 'a':tokens[1].strip()} try: # is a sticking coefficient - kinetics['sticking coefficient'][species_dict[tokens[0].strip()]] = {'E': (tokens[1].strip(), Eunits), 'm': tokens[2].strip(), 'a':tokens[3].strip()} + kinetics['sticking coefficient'][species_dict[tokens[0].strip()]] = {'E': (tokens[3].strip(), Eunits), 'm': tokens[2].strip(), 'a':tokens[1].strip()} except KeyError: # is a not a sticking coefficient - kinetics['surface arrhenius'].coverage_dependence[species_dict[tokens[0].strip()]] = {'E': (tokens[1].strip(), 'J/mol'), 'm': tokens[2].strip(), 'a':tokens[3].strip()} + kinetics['surface arrhenius'].coverage_dependence[species_dict[tokens[0].strip()]] = {'E': (tokens[3].strip(), 'J/mol'), 'm': tokens[2].strip(), 'a':tokens[1].strip()} elif 'LOW' in line: # Low-pressure-limit Arrhenius parameters @@ -1818,7 +1818,7 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals for species, cov_params in kinetics.coverage_dependence.items(): label = get_species_identifier(species) string += f'\n COV / {label:<41}' - string += f"{cov_params['E'].value_si:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['a'].value:<6.3f} /" + string += f"{cov_params['a'].value:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['E'].value_si:<6.3f} /" elif isinstance(kinetics, _kinetics.Arrhenius): conversion_factor = kinetics.A.get_conversion_factor_from_si_to_cm_mol_s() if not isinstance(kinetics, _kinetics.SurfaceArrhenius): @@ -1840,7 +1840,7 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals for species, cov_params in kinetics.coverage_dependence.items(): label = get_species_identifier(species) string += f'\n COV / {label:<41}' - string += f"{cov_params['E'].value_si:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['a'].value:<6.3f} /" + string += f"{cov_params['a'].value:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['E'].value_si:<6.3f} /" elif isinstance(kinetics, (_kinetics.Lindemann, _kinetics.Troe)): arrhenius = kinetics.arrheniusHigh conversion_factor = arrhenius.A.get_conversion_factor_from_si_to_cm_mol_s() From 9bb42c18cfef9bc4f0d8728af29d1c888fa84167 Mon Sep 17 00:00:00 2001 From: Richard West Date: Wed, 30 Jun 2021 23:59:18 -0400 Subject: [PATCH 21/43] Chemkin reading COV: Do parameters in order a, m, E. also, turning the forgotten 'J/mol' into Eunits ! --- rmgpy/chemkin.pyx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index e68167247e..54e979b77a 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -463,13 +463,13 @@ def _read_kinetics_line(line, reaction, species_dict, Eunits, kunits, klow_units del kinetics['arrhenius high'] tokens = tokens[1].split() - kinetics['coverage dependence'][species_dict[tokens[0].strip()]] = {'E': (tokens[3].strip(), Eunits), 'm': tokens[2].strip(), 'a':tokens[1].strip()} + kinetics['coverage dependence'][species_dict[tokens[0].strip()]] = {'a':tokens[1].strip(), 'm': tokens[2].strip(), 'E': (tokens[3].strip(), Eunits)} try: # is a sticking coefficient - kinetics['sticking coefficient'][species_dict[tokens[0].strip()]] = {'E': (tokens[3].strip(), Eunits), 'm': tokens[2].strip(), 'a':tokens[1].strip()} + kinetics['sticking coefficient'][species_dict[tokens[0].strip()]] = {'a':tokens[1].strip(), 'm': tokens[2].strip(), 'E': (tokens[3].strip(), Eunits)} except KeyError: # is a not a sticking coefficient - kinetics['surface arrhenius'].coverage_dependence[species_dict[tokens[0].strip()]] = {'E': (tokens[3].strip(), 'J/mol'), 'm': tokens[2].strip(), 'a':tokens[1].strip()} + kinetics['surface arrhenius'].coverage_dependence[species_dict[tokens[0].strip()]] = {'a':tokens[1].strip(), 'm': tokens[2].strip(), 'E': (tokens[3].strip(), Eunits)} elif 'LOW' in line: # Low-pressure-limit Arrhenius parameters From 9e17c6e1b7a5e6ffac3e278cc4220dd349a8a443 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 1 Jul 2021 00:07:38 -0400 Subject: [PATCH 22/43] Chemkin reading COV parameters: turn them into float() values --- rmgpy/chemkin.pyx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index 54e979b77a..72f61f7b91 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -463,13 +463,13 @@ def _read_kinetics_line(line, reaction, species_dict, Eunits, kunits, klow_units del kinetics['arrhenius high'] tokens = tokens[1].split() - kinetics['coverage dependence'][species_dict[tokens[0].strip()]] = {'a':tokens[1].strip(), 'm': tokens[2].strip(), 'E': (tokens[3].strip(), Eunits)} + kinetics['coverage dependence'][species_dict[tokens[0].strip()]] = {'a':float(tokens[1]), 'm':float(tokens[2]), 'E':(float(tokens[3]), Eunits)} try: # is a sticking coefficient - kinetics['sticking coefficient'][species_dict[tokens[0].strip()]] = {'a':tokens[1].strip(), 'm': tokens[2].strip(), 'E': (tokens[3].strip(), Eunits)} + kinetics['sticking coefficient'][species_dict[tokens[0].strip()]] = {'a':float(tokens[1]), 'm':float(tokens[2]), 'E':(float(tokens[3]), Eunits)} except KeyError: # is a not a sticking coefficient - kinetics['surface arrhenius'].coverage_dependence[species_dict[tokens[0].strip()]] = {'a':tokens[1].strip(), 'm': tokens[2].strip(), 'E': (tokens[3].strip(), Eunits)} + kinetics['surface arrhenius'].coverage_dependence[species_dict[tokens[0].strip()]] = {'a':float(tokens[1]), 'm':float(tokens[2]), 'E': (float(tokens[3]), Eunits)} elif 'LOW' in line: # Low-pressure-limit Arrhenius parameters From fec8618e4d36037d9c3534e9f46d2dbc4518e446 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 1 Jul 2021 00:12:02 -0400 Subject: [PATCH 23/43] Chemkin reading COV parameters: should preserve case of species names. If we're even going to bother writing this code we really ought to test it --- rmgpy/chemkin.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index 72f61f7b91..9062dbc666 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -462,7 +462,7 @@ def _read_kinetics_line(line, reaction, species_dict, Eunits, kunits, klow_units ) del kinetics['arrhenius high'] - tokens = tokens[1].split() + tokens = case_preserved_tokens[1].split() kinetics['coverage dependence'][species_dict[tokens[0].strip()]] = {'a':float(tokens[1]), 'm':float(tokens[2]), 'E':(float(tokens[3]), Eunits)} try: # is a sticking coefficient From e670ceeda2e506bfca3144da3045584a379ab58b Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 1 Jul 2021 00:27:28 -0400 Subject: [PATCH 24/43] Put coverage-dependence paramaters in order a, m, E, consistently. Dictionaries are order-preserving, and humans' brains are, so we'll have fewer bugs if we're consistent with ourselves and other software's conventions. Chemkin and Cantera use (a, m, E) ordering. Now we do, throughout. --- rmgpy/chemkinTest.py | 2 +- rmgpy/kinetics/surface.pyx | 24 ++++++++++++------------ rmgpy/kinetics/surfaceTest.py | 10 ++++++---- rmgpy/solver/surfaceTest.py | 12 ++++++------ 4 files changed, 25 insertions(+), 23 deletions(-) diff --git a/rmgpy/chemkinTest.py b/rmgpy/chemkinTest.py index 702dd52a08..6b8bf1ba19 100644 --- a/rmgpy/chemkinTest.py +++ b/rmgpy/chemkinTest.py @@ -470,7 +470,7 @@ def test_read_coverage_dependence(self): Ea=(5.0, 'kJ/mol'), T0=(1.0, 'K'), coverage_dependence={ - s_x: {'E': (0.1, 'J/mol'), 'm': -1.0, 'a': 1.0},})) + s_x: {'a': 1.0, 'm': -1.0, 'E': (0.1, 'J/mol')},})) reactions = [self.rxn_covdep] diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx index 542cdc99c8..96998e8de6 100644 --- a/rmgpy/kinetics/surface.pyx +++ b/rmgpy/kinetics/surface.pyx @@ -59,9 +59,9 @@ cdef class StickingCoefficient(KineticsModel): `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: - `E`, the activation energy dependence on coverage, + `a`, the coefficient for exponential dependence on the coverage, `m`, the power-law exponent of coverage dependence, and - `a`, the coefficient for exponential dependence on the coverage + `E`, the activation energy dependence on coverage. `comment` Information about the model (e.g. its source) ======================= ============================================================= @@ -87,7 +87,7 @@ cdef class StickingCoefficient(KineticsModel): if self.coverage_dependence is not None: string += ", coverage_dependence={" for species, parameters in self.coverage_dependence.items(): - string += f"{species.to_chemkin()!r}: {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," + 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) @@ -264,9 +264,9 @@ cdef class StickingCoefficientBEP(KineticsModel): `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: - `E`, the activation energy dependence on coverage, + `a`, the coefficient for exponential dependence on the coverage, `m`, the power-law exponent of coverage dependence, and - `a`, the coefficient for exponential dependence on the coverage + `E`, the activation energy dependence on coverage. `comment` Information about the model (e.g. its source) ======================= ============================================================= @@ -293,7 +293,7 @@ cdef class StickingCoefficientBEP(KineticsModel): if self.coverage_dependence is not None: string += ", coverage_dependence={" for species, parameters in self.coverage_dependence.items(): - string += f"{species.to_chemkin()!r}: {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," + 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) @@ -445,9 +445,9 @@ cdef class SurfaceArrhenius(Arrhenius): `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: - `E`, the activation energy dependence on coverage, + `a`, the coefficient for exponential dependence on the coverage, `m`, the power-law exponent of coverage dependence, and - `a`, the coefficient for exponential dependence on the coverage + `E`, the activation energy dependence on coverage. `uncertainty` Uncertainty information `comment` Information about the model (e.g. its source) ======================= ============================================================= @@ -495,7 +495,7 @@ cdef class SurfaceArrhenius(Arrhenius): if self.coverage_dependence is not None: string += ", coverage_dependence={" for species, parameters in self.coverage_dependence.items(): - string += f"{species.to_chemkin()!r}: {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," + 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) @@ -537,9 +537,9 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): `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: - `E`, the activation energy dependence on coverage, + `a`, the coefficient for exponential dependence on the coverage, `m`, the power-law exponent of coverage dependence, and - `a`, the coefficient for exponential dependence on the coverage + `E`, the activation energy dependence on coverage. `uncertainty` Uncertainty information `comment` Information about the model (e.g. its source) ======================= ============================================================= @@ -589,7 +589,7 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): if self.coverage_dependence is not None: string += ", coverage_dependence={" for species, parameters in self.coverage_dependence.items(): - string += f"{species.to_chemkin()!r}: {{'E':({parameters['E'].value}, '{parameters['E'].units}'), 'm':{parameters['m']}, 'a':{parameters['a']}}}," + 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) diff --git a/rmgpy/kinetics/surfaceTest.py b/rmgpy/kinetics/surfaceTest.py index 8625600d47..f5a4a13ee0 100644 --- a/rmgpy/kinetics/surfaceTest.py +++ b/rmgpy/kinetics/surfaceTest.py @@ -60,8 +60,9 @@ def setUp(self): self.Tmax = 3000. s = Species().from_adjacency_list('1 X u0 p0 c0') s.label = 'X' - self.coverage_dependence = {s: {'E': quantity.Energy(0.0, 'J/mol'), 'm': quantity.Dimensionless(-1.0), - 'a': quantity.Dimensionless(0.0)}} + 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, @@ -268,8 +269,9 @@ def setUp(self): self.Tmax = 3000. s = Species().from_adjacency_list('1 X u0 p0 c0') s.label = 'X' - self.coverage_dependence = {s: {'E': quantity.Energy(0.0, 'J/mol'), 'm': quantity.Dimensionless(-1.0), - 'a': quantity.Dimensionless(0.0)}} + 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)'), diff --git a/rmgpy/solver/surfaceTest.py b/rmgpy/solver/surfaceTest.py index 121c0be0db..0e91a1171b 100644 --- a/rmgpy/solver/surfaceTest.py +++ b/rmgpy/solver/surfaceTest.py @@ -315,7 +315,7 @@ def test_solve_h2_coverage_dependence(self): n=0.5, Ea=(5.0, 'kJ/mol'), T0=(1.0, 'K'), - coverage_dependence={'*': {'E': (0.0, 'J/mol'), 'm': -1.0, 'a': 0.0}})) + coverage_dependence={'*': {'a': 0.0, 'm': -1.0, 'E': (0.0, 'J/mol')}})) core_species = [h2, x, hx] edge_species = [] @@ -342,9 +342,9 @@ def test_solve_h2_coverage_dependence(self): for species, parameters in rxn1.kinetics.coverage_dependence.items(): self.assertIsInstance(species, str) # species should be a string self.assertIsInstance(parameters, dict) - self.assertIsNotNone(parameters['E']) - self.assertIsNotNone(parameters['m']) self.assertIsNotNone(parameters['a']) + self.assertIsNotNone(parameters['m']) + self.assertIsNotNone(parameters['E']) # Integrate to get the solution at each time point t = [] @@ -435,7 +435,7 @@ def test_solve_ch3_coverage_dependence(self): products=[ch3x], kinetics=StickingCoefficient( A=0.1, n=0, Ea=(0, 'kcal/mol'), T0=(1, 'K'), Tmin=(200, 'K'), Tmax=(3000, 'K'), - coverage_dependence={'*': {'E': (0.0, 'J/mol'), 'm': -1.0, 'a': 0.0}}, + coverage_dependence={'*': {'a': 0.0, 'm': -1.0, 'E': (0.0, 'J/mol')}}, comment="""Exact match found for rate rule (Adsorbate;VacantSite)""" ) ) @@ -470,9 +470,9 @@ def test_solve_ch3_coverage_dependence(self): for species, parameters in rxn1.kinetics.coverage_dependence.items(): self.assertIsInstance(species, str) # species should be a string self.assertIsInstance(parameters, dict) - self.assertIsNotNone(parameters['E']) - self.assertIsNotNone(parameters['m']) self.assertIsNotNone(parameters['a']) + self.assertIsNotNone(parameters['m']) + self.assertIsNotNone(parameters['E']) # Integrate to get the solution at each time point t = [] From 27978dd0512dde542b39970edd392a0db9026b00 Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 2 Jul 2021 22:31:34 -0400 Subject: [PATCH 25/43] Chemkin reading COV parameters. Fixes for reading sticking coefficients. If there's a STICK line before a COV line we can cope. --- rmgpy/chemkin.pyx | 33 +++++++++++++++++---------------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index 9062dbc666..d5f1feb82a 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -207,7 +207,6 @@ def read_kinetics_entry(entry, species_dict, Aunits, Eunits): kinetics.update({ 'chebyshev coefficients': [], 'efficiencies': {}, - 'coverage dependence': {} }) # Note that the subsequent lines could be in any order @@ -453,23 +452,22 @@ def _read_kinetics_line(line, reaction, species_dict, Eunits, kunits, klow_units reaction.duplicate = True elif 'COV' in line: - k = kinetics['arrhenius high'] - kinetics['surface arrhenius'] = _kinetics.SurfaceArrhenius( - A=(k.A.value, kunits), - n=k.n, - Ea=k.Ea, - T0=k.T0, - ) - del kinetics['arrhenius high'] - - tokens = case_preserved_tokens[1].split() - kinetics['coverage dependence'][species_dict[tokens[0].strip()]] = {'a':float(tokens[1]), 'm':float(tokens[2]), 'E':(float(tokens[3]), Eunits)} try: - # is a sticking coefficient - kinetics['sticking coefficient'][species_dict[tokens[0].strip()]] = {'a':float(tokens[1]), 'm':float(tokens[2]), 'E':(float(tokens[3]), Eunits)} + k = kinetics['sticking coefficient'] except KeyError: - # is a not a sticking coefficient - kinetics['surface arrhenius'].coverage_dependence[species_dict[tokens[0].strip()]] = {'a':float(tokens[1]), 'm':float(tokens[2]), 'E': (float(tokens[3]), Eunits)} + 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 @@ -563,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() @@ -571,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'] @@ -581,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: From ad508c90f5416735095b14d9936b1f5db59bd6d6 Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 2 Jul 2021 22:33:10 -0400 Subject: [PATCH 26/43] Chemkin COV writing: correct units for the E parameter The E should be in the same units as Ea, which is kcal/mol. That's 4184 times smaller. --- rmgpy/chemkin.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index d5f1feb82a..115fb317d6 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -1841,7 +1841,7 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals for species, cov_params in kinetics.coverage_dependence.items(): label = get_species_identifier(species) string += f'\n COV / {label:<41}' - string += f"{cov_params['a'].value:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['E'].value_si:<6.3f} /" + string += f"{cov_params['a'].value:<9.3g} {cov_params['m'].value:<6.3f} {cov_params['E'].value_si/4184.:<6.3f} /" elif isinstance(kinetics, (_kinetics.Lindemann, _kinetics.Troe)): arrhenius = kinetics.arrheniusHigh conversion_factor = arrhenius.A.get_conversion_factor_from_si_to_cm_mol_s() From b12ded625793b5b30b6b562a566f58dc75944d53 Mon Sep 17 00:00:00 2001 From: mazeau Date: Tue, 6 Jul 2021 15:48:36 -0400 Subject: [PATCH 27/43] update __repr__ string writing do not write if the coverage dependence dictionary is empty --- rmgpy/kinetics/surface.pyx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx index 96998e8de6..2c108e9204 100644 --- a/rmgpy/kinetics/surface.pyx +++ b/rmgpy/kinetics/surface.pyx @@ -84,7 +84,7 @@ 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 is not None: + 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}')}}," @@ -290,7 +290,7 @@ 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 is not None: + 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}')}}," @@ -492,7 +492,7 @@ 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 is not None: + 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}')}}," @@ -586,7 +586,7 @@ 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 is not None: + 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}')}}," From dd5f1192b5462ebe1cbf39a4d452efdee938ff9d Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 8 Jul 2021 13:21:19 -0400 Subject: [PATCH 28/43] Revert "Do not write coverage dependent parameters if coverage dependence is turned off" This reverts commit 7841f2496620c1c36556adec886e182fd457da79. Had to manually tweak a few things to resolve conflicts. The point is to remove the extra code from the chemkin-writing methods. If we turn off coverage, the parameters should not be in the model. If the parameters are in the model, they should be written by chemkin. --- rmgpy/chemkin.pyx | 41 +++++++++++++++++------------------------ rmgpy/chemkinTest.py | 3 +-- 2 files changed, 18 insertions(+), 26 deletions(-) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index 115fb317d6..a73a276ddc 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -1711,8 +1711,7 @@ def write_reaction_string(reaction, java_library=False): ################################################################################ -def write_kinetics_entry(reaction, species_list, verbose=True, java_library=False, commented=False, - coverage_dependence=False): +def write_kinetics_entry(reaction, species_list, verbose=True, java_library=False, commented=False): """ Return a string representation of the reaction as used in a Chemkin file. Use `verbose = True` to turn on kinetics comments. @@ -1744,8 +1743,7 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals specific_collider=reaction.specific_collider, reversible=reaction.reversible, kinetics=kinetics) - string += write_kinetics_entry(new_reaction, species_list, verbose, java_library, commented, - coverage_dependence=coverage_dependence) + string += write_kinetics_entry(new_reaction, species_list, verbose, java_library, commented) string += "DUPLICATE\n" if commented: @@ -1814,12 +1812,11 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals kinetics.Ea.value_si / 4184. ) string += '\n STICK' - if coverage_dependence: - if kinetics.coverage_dependence is not None: - for species, cov_params in kinetics.coverage_dependence.items(): - label = get_species_identifier(species) - string += f'\n COV / {label:<41}' - string += f"{cov_params['a'].value:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['E'].value_si:<6.3f} /" + if kinetics.coverage_dependence is not None: + for species, cov_params in kinetics.coverage_dependence.items(): + label = get_species_identifier(species) + string += f'\n COV / {label:<41}' + string += f"{cov_params['a'].value:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['E'].value_si:<6.3f} /" elif isinstance(kinetics, _kinetics.Arrhenius): conversion_factor = kinetics.A.get_conversion_factor_from_si_to_cm_mol_s() if not isinstance(kinetics, _kinetics.SurfaceArrhenius): @@ -1837,11 +1834,10 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals kinetics.Ea.value_si / 4184. ) if isinstance(kinetics, _kinetics.SurfaceArrhenius) and kinetics.coverage_dependence: - if coverage_dependence is True: - for species, cov_params in kinetics.coverage_dependence.items(): - label = get_species_identifier(species) - string += f'\n COV / {label:<41}' - string += f"{cov_params['a'].value:<9.3g} {cov_params['m'].value:<6.3f} {cov_params['E'].value_si/4184.:<6.3f} /" + for species, cov_params in kinetics.coverage_dependence.items(): + label = get_species_identifier(species) + string += f'\n COV / {label:<41}' + string += f"{cov_params['a'].value:<9.3g} {cov_params['m'].value:<6.3f} {cov_params['E'].value_si/4184.:<6.3f} /" elif isinstance(kinetics, (_kinetics.Lindemann, _kinetics.Troe)): arrhenius = kinetics.arrheniusHigh conversion_factor = arrhenius.A.get_conversion_factor_from_si_to_cm_mol_s() @@ -2108,7 +2104,7 @@ def save_transport_file(path, species): )) -def save_chemkin_file(path, species, reactions, verbose=True, check_for_duplicates=True, coverage_dependence=False): +def save_chemkin_file(path, species, reactions, verbose=True, check_for_duplicates=True): """ Save a Chemkin input file to `path` on disk containing the provided lists of `species` and `reactions`. @@ -2156,7 +2152,7 @@ def save_chemkin_file(path, species, reactions, verbose=True, check_for_duplicat global _chemkin_reaction_count _chemkin_reaction_count = 0 for rxn in reactions: - f.write(write_kinetics_entry(rxn, species_list=species, verbose=verbose, coverage_dependence=coverage_dependence)) + f.write(write_kinetics_entry(rxn, species_list=species, verbose=verbose)) # Don't forget to mark duplicates! f.write('\n') f.write('END\n\n') @@ -2166,7 +2162,7 @@ def save_chemkin_file(path, species, reactions, verbose=True, check_for_duplicat def save_chemkin_surface_file(path, species, reactions, verbose=True, check_for_duplicates=True, - surface_site_density=None, coverage_dependence=False): + surface_site_density=None): """ Save a Chemkin *surface* input file to `path` on disk containing the provided lists of `species` and `reactions`. @@ -2216,8 +2212,7 @@ def save_chemkin_surface_file(path, species, reactions, verbose=True, check_for_ global _chemkin_reaction_count _chemkin_reaction_count = 0 for rxn in reactions: - f.write(write_kinetics_entry(rxn, species_list=species, verbose=verbose, - coverage_dependence=coverage_dependence)) + f.write(write_kinetics_entry(rxn, species_list=species, verbose=verbose)) f.write('\n') f.write('END\n\n') f.close() @@ -2310,13 +2305,11 @@ def save_chemkin(reaction_model, path, verbose_path, dictionary_path=None, trans # We should already have marked everything as duplicates by now so use check_for_duplicates=False save_chemkin_file(gas_path, gas_species_list, gas_rxn_list, verbose=False, check_for_duplicates=False) save_chemkin_surface_file(surface_path, surface_species_list, surface_rxn_list, verbose=False, - check_for_duplicates=False, surface_site_density=reaction_model.surface_site_density, - coverage_dependence=reaction_model.coverage_dependence) + check_for_duplicates=False, surface_site_density=reaction_model.surface_site_density) logging.info('Saving annotated version of Chemkin files...') save_chemkin_file(gas_verbose_path, gas_species_list, gas_rxn_list, verbose=True, check_for_duplicates=False) save_chemkin_surface_file(surface_verbose_path, surface_species_list, surface_rxn_list, verbose=True, - check_for_duplicates=False, surface_site_density=reaction_model.surface_site_density, - coverage_dependence=reaction_model.coverage_dependence) + check_for_duplicates=False, surface_site_density=reaction_model.surface_site_density) else: # Gas phase only diff --git a/rmgpy/chemkinTest.py b/rmgpy/chemkinTest.py index 6b8bf1ba19..cf1c006212 100644 --- a/rmgpy/chemkinTest.py +++ b/rmgpy/chemkinTest.py @@ -477,8 +477,7 @@ def test_read_coverage_dependence(self): # 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, - coverage_dependence=True) + 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 From 61f7c88df966bbfbc84e302b842b8265aee341c5 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 8 Jul 2021 13:34:46 -0400 Subject: [PATCH 29/43] Chemkin writing coverage dependence: fix bugs - An empty dictionary shouldn't be printed (even though it's not 'None') - We need sticking coefficients cov parameters to be in the right units too! (we previously fixed from J/mol to kcal/mol in one place, but not both) --- rmgpy/chemkin.pyx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index a73a276ddc..d507484499 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -1812,11 +1812,11 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals kinetics.Ea.value_si / 4184. ) string += '\n STICK' - if kinetics.coverage_dependence is not None: + if kinetics.coverage_dependence: for species, cov_params in kinetics.coverage_dependence.items(): label = get_species_identifier(species) string += f'\n COV / {label:<41}' - string += f"{cov_params['a'].value:<9.3e} {cov_params['m'].value:<6.3f} {cov_params['E'].value_si:<6.3f} /" + string += f"{cov_params['a'].value:<9.3g} {cov_params['m'].value:<6.3g} {cov_params['E'].value_si/4184:<6.3f} /" elif isinstance(kinetics, _kinetics.Arrhenius): conversion_factor = kinetics.A.get_conversion_factor_from_si_to_cm_mol_s() if not isinstance(kinetics, _kinetics.SurfaceArrhenius): @@ -1837,7 +1837,7 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals for species, cov_params in kinetics.coverage_dependence.items(): label = get_species_identifier(species) string += f'\n COV / {label:<41}' - string += f"{cov_params['a'].value:<9.3g} {cov_params['m'].value:<6.3f} {cov_params['E'].value_si/4184.:<6.3f} /" + string += f"{cov_params['a'].value:<9.3g} {cov_params['m'].value:<6.3g} {cov_params['E'].value_si/4184.:<6.3f} /" elif isinstance(kinetics, (_kinetics.Lindemann, _kinetics.Troe)): arrhenius = kinetics.arrheniusHigh conversion_factor = arrhenius.A.get_conversion_factor_from_si_to_cm_mol_s() From 78f87d8267f97847138e6be7afb7b154fd0c5184 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 8 Jul 2021 13:44:23 -0400 Subject: [PATCH 30/43] Chemkin writing coverage dependence: reduce code duplication Protects us from bugs where we fix things in one place and not the other. --- rmgpy/chemkin.pyx | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/rmgpy/chemkin.pyx b/rmgpy/chemkin.pyx index d507484499..9ca55a00b9 100644 --- a/rmgpy/chemkin.pyx +++ b/rmgpy/chemkin.pyx @@ -1812,11 +1812,6 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals kinetics.Ea.value_si / 4184. ) string += '\n STICK' - if kinetics.coverage_dependence: - for species, cov_params in kinetics.coverage_dependence.items(): - label = get_species_identifier(species) - string += f'\n COV / {label:<41}' - string += f"{cov_params['a'].value:<9.3g} {cov_params['m'].value:<6.3g} {cov_params['E'].value_si/4184:<6.3f} /" elif isinstance(kinetics, _kinetics.Arrhenius): conversion_factor = kinetics.A.get_conversion_factor_from_si_to_cm_mol_s() if not isinstance(kinetics, _kinetics.SurfaceArrhenius): @@ -1833,11 +1828,6 @@ def write_kinetics_entry(reaction, species_list, verbose=True, java_library=Fals kinetics.n.value_si, kinetics.Ea.value_si / 4184. ) - if isinstance(kinetics, _kinetics.SurfaceArrhenius) and kinetics.coverage_dependence: - for species, cov_params in kinetics.coverage_dependence.items(): - label = get_species_identifier(species) - string += f'\n COV / {label:<41}' - string += f"{cov_params['a'].value:<9.3g} {cov_params['m'].value:<6.3g} {cov_params['E'].value_si/4184.:<6.3f} /" elif isinstance(kinetics, (_kinetics.Lindemann, _kinetics.Troe)): arrhenius = kinetics.arrheniusHigh conversion_factor = arrhenius.A.get_conversion_factor_from_si_to_cm_mol_s() @@ -1877,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])): From c1099110f15542b8c0aafb95e66ef60e64002150 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 8 Jul 2021 14:42:21 -0400 Subject: [PATCH 31/43] Renaming apply_coverage_depndence_to_reaction to process_coverage_dependence. This better describes what it does. Also expanded docstring. Also removed it from INSIDE the apply_kinetics_to_reaction method and removed if from and elif block, so that it is duplicated in less code (but still called in the same way, I think) --- rmgpy/rmg/model.py | 31 +++++++++++++++---------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index ebdf5ed8c3..108a1c5ae8 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -649,12 +649,12 @@ 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) - elif hasattr(reaction.kinetics, 'coverage_dependence'): - if reaction.kinetics.coverage_dependence: - self.apply_coverage_dependence_to_reaction(reaction.kinetics) + + 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. @@ -868,8 +868,13 @@ def generate_thermo(self, spc, rename=False): spc.generate_energy_transfer_model() - def apply_coverage_dependence_to_reaction(self, kinetics): - """Apply the coverage dependence kinetics""" + 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). + """ cov_dep = {} for species, values in kinetics.coverage_dependence.items(): species_in_model, is_new = self.make_new_species(species) @@ -898,10 +903,6 @@ def apply_kinetics_to_reaction(self, reaction): reaction.reverse = None reaction.kinetics = kinetics - if hasattr(kinetics, 'coverage_dependence'): - if reaction.kinetics.coverage_dependence: - self.apply_coverage_dependence_to_reaction(kinetics) - def generate_kinetics(self, reaction): """ Generate best possible kinetics for the given `reaction` using the kinetics database. @@ -1511,9 +1512,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 hasattr(r.kinetics, 'coverage_dependence'): - if r.kinetics.coverage_dependence: - self.apply_coverage_dependence_to_reaction(r.kinetics) + 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() \ @@ -1613,9 +1613,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 hasattr(r.kinetics, 'coverage_dependence'): - if r.kinetics.coverage_dependence: - self.apply_coverage_dependence_to_reaction(r.kinetics) + 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() \ From 6494c17b80fbc2756710e1117dd6b2eadc1c9f5d Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 8 Jul 2021 14:44:07 -0400 Subject: [PATCH 32/43] process_coverage_dependence now REMOVES parameters if turned off. I'm still not very sure about the `coverage_dependence=False` settings, and what should be turned off by it. --- rmgpy/rmg/model.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index 108a1c5ae8..840d7d1550 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -874,7 +874,13 @@ def process_coverage_dependence(self, 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) From bc99a4e2a71b8cbceb770afe32a3bec8d1f1f934 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 8 Jul 2021 20:45:03 -0400 Subject: [PATCH 33/43] SurfaceReactor: Remove unused (and mutable) default arguments. Putting dictionaries as a default argument is a famous gotcha https://docs.python-guide.org/writing/gotchas/ I wondered if it's being done on purpose for some efficiency reason, to share memory, but decided no. a risk of error, and is not needed anyway as we never initialize the object with those parameters anyway. --- rmgpy/solver/surface.pyx | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index 495d4c6d18..440ec6f44d 100644 --- a/rmgpy/solver/surface.pyx +++ b/rmgpy/solver/surface.pyx @@ -85,9 +85,6 @@ cdef class SurfaceReactor(ReactionSystem): sensitivity_threshold=1e-3, sens_conditions=None, coverage_dependence=False, - cov_dep={}, - cov_dep_index_species={}, - current_surface_coverages={}, ): ReactionSystem.__init__(self, termination, @@ -111,9 +108,10 @@ cdef class SurfaceReactor(ReactionSystem): self.constant_volume = True self.sens_conditions = sens_conditions self.n_sims = n_sims - self.cov_dep = cov_dep - self.cov_dep_index_species = cov_dep_index_species - self.current_surface_coverages = current_surface_coverages + + self.cov_dep = {} + self.cov_dep_index_species = {} + self.current_surface_coverages = {} def convert_initial_keys_to_species_objects(self, species_dict): """ From 7eda4c3a3845465304b5e50d0bc7012f9acb822b Mon Sep 17 00:00:00 2001 From: Richard West Date: Mon, 12 Jul 2021 16:10:22 -0400 Subject: [PATCH 34/43] CovDep test: report timing of the simulations. I'm hoping to speed it up later. --- rmgpy/solver/surfaceTest.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/rmgpy/solver/surfaceTest.py b/rmgpy/solver/surfaceTest.py index 0e91a1171b..16f37872c6 100644 --- a/rmgpy/solver/surfaceTest.py +++ b/rmgpy/solver/surfaceTest.py @@ -28,6 +28,7 @@ ############################################################################### import unittest +import time import numpy as np @@ -351,6 +352,7 @@ def test_solve_h2_coverage_dependence(self): y = [] reaction_rates = [] species_rates = [] + start_time = time.time() for t1 in tlist: rxn_system.advance(t1) t.append(rxn_system.t) @@ -359,6 +361,8 @@ def test_solve_h2_coverage_dependence(self): 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) @@ -489,6 +493,7 @@ def test_solve_ch3_coverage_dependence(self): 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) @@ -497,6 +502,8 @@ def test_solve_ch3_coverage_dependence(self): 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) From 6dee14584a753e600e3bbe2f41efa6be3692022d Mon Sep 17 00:00:00 2001 From: Richard West Date: Mon, 12 Jul 2021 16:12:24 -0400 Subject: [PATCH 35/43] CovDep test: plot graphs. They look the same :-/ --- rmgpy/solver/surfaceTest.py | 50 ++++++++++++++++++++++++------------- 1 file changed, 33 insertions(+), 17 deletions(-) diff --git a/rmgpy/solver/surfaceTest.py b/rmgpy/solver/surfaceTest.py index 16f37872c6..7cbc94e26f 100644 --- a/rmgpy/solver/surfaceTest.py +++ b/rmgpy/solver/surfaceTest.py @@ -132,23 +132,22 @@ def test_solve_h2(self): # Check that we've reached equilibrium self.assertAlmostEqual(reaction_rates[-1, 0], 0.0, delta=1e-2) - # # Visualize the simulation results - # 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.legend(['HX'], loc=4) - # pylab.subplot(2, 1, 2) - # pylab.semilogx(t, species_rates) - # pylab.legend(['H2', 'X', 'HX'], loc=4) - # pylab.xlabel('Time (s)') - # pylab.ylabel('Rate (mol/m$^\\mathdefault{3 or 2}$*s)') - # # fig.subplots_adjust(left=0.21, bottom=0.10, right=0.95, top=0.95, wspace=0.20, hspace=0.35) - # pylab.tight_layout() - # # pylab.show() - # pylab.savefig('surfaceTestH2.pdf') - + # Visualize the simulation results + 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.legend(['HX'], loc=4) + pylab.subplot(2, 1, 2) + pylab.semilogx(t, species_rates) + pylab.legend(['H2', 'X', 'HX'], loc=4) + pylab.xlabel('Time (s)') + pylab.ylabel('Rate (mol/m$^\\mathdefault{3 or 2}$*s)') + # fig.subplots_adjust(left=0.21, bottom=0.10, right=0.95, top=0.95, wspace=0.20, hspace=0.35) + pylab.tight_layout() + # pylab.show() + pylab.savefig('surfaceTestH2.pdf') return def test_solve_ch3(self): @@ -383,6 +382,23 @@ def test_solve_h2_coverage_dependence(self): # Check that we've reached equilibrium self.assertAlmostEqual(reaction_rates[-1, 0], 0.0, delta=1e-2) + # Visualize the simulation results + 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.legend(['HX'], loc=4) + pylab.subplot(2, 1, 2) + pylab.semilogx(t, species_rates) + pylab.legend(['H2', 'X', 'HX'], loc=4) + pylab.xlabel('Time (s)') + pylab.ylabel('Rate (mol/m$^\\mathdefault{3 or 2}$*s)') + # fig.subplots_adjust(left=0.21, bottom=0.10, right=0.95, top=0.95, wspace=0.20, hspace=0.35) + pylab.tight_layout() + # pylab.show() + pylab.savefig('surfaceTestH2covdep.pdf') + def test_solve_ch3_coverage_dependence(self): """ Test the surface batch reactor can properly apply coverage dependent parameters From 7425b4e66b5a369bf302231c8d3625b51b2b32b8 Mon Sep 17 00:00:00 2001 From: Richard West Date: Mon, 12 Jul 2021 16:41:04 -0400 Subject: [PATCH 36/43] CovDep test: species key should be a Species --- rmgpy/solver/surfaceTest.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/rmgpy/solver/surfaceTest.py b/rmgpy/solver/surfaceTest.py index 7cbc94e26f..a2494706d7 100644 --- a/rmgpy/solver/surfaceTest.py +++ b/rmgpy/solver/surfaceTest.py @@ -315,7 +315,7 @@ def test_solve_h2_coverage_dependence(self): n=0.5, Ea=(5.0, 'kJ/mol'), T0=(1.0, 'K'), - coverage_dependence={'*': {'a': 0.0, 'm': -1.0, 'E': (0.0, 'J/mol')}})) + coverage_dependence={x: {'a': 0.0, 'm': -1.0, 'E': (0.0, 'J/mol')}})) core_species = [h2, x, hx] edge_species = [] @@ -340,7 +340,7 @@ def test_solve_h2_coverage_dependence(self): 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, str) # species should be a string + self.assertIsInstance(species, Species) # species should be a Species self.assertIsInstance(parameters, dict) self.assertIsNotNone(parameters['a']) self.assertIsNotNone(parameters['m']) @@ -455,7 +455,7 @@ def test_solve_ch3_coverage_dependence(self): products=[ch3x], kinetics=StickingCoefficient( A=0.1, n=0, Ea=(0, 'kcal/mol'), T0=(1, 'K'), Tmin=(200, 'K'), Tmax=(3000, 'K'), - coverage_dependence={'*': {'a': 0.0, 'm': -1.0, 'E': (0.0, 'J/mol')}}, + coverage_dependence={x: {'a': 0.0, 'm': -1.0, 'E': (0.0, 'J/mol')}}, comment="""Exact match found for rate rule (Adsorbate;VacantSite)""" ) ) @@ -488,7 +488,7 @@ def test_solve_ch3_coverage_dependence(self): 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, str) # species should be a string + self.assertIsInstance(species, Species) # species should be a Species self.assertIsInstance(parameters, dict) self.assertIsNotNone(parameters['a']) self.assertIsNotNone(parameters['m']) From bbbba5e7e85b9a02433b565460ca0c7a4676c33c Mon Sep 17 00:00:00 2001 From: Richard West Date: Mon, 12 Jul 2021 16:45:01 -0400 Subject: [PATCH 37/43] CovDep test: change T, and plot coverage. At 600 K the eqm coverage is much higher, so you get a more noticeable coverage dependence. By plotting coverage instead of concentration it's easier to see what's happening. --- rmgpy/solver/surfaceTest.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/rmgpy/solver/surfaceTest.py b/rmgpy/solver/surfaceTest.py index a2494706d7..4c72f69be0 100644 --- a/rmgpy/solver/surfaceTest.py +++ b/rmgpy/solver/surfaceTest.py @@ -84,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, @@ -118,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]): @@ -136,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) @@ -322,7 +322,7 @@ def test_solve_h2_coverage_dependence(self): core_reactions = [rxn1] edge_reactions = [] - T = 1000 + T = 600 P_initial = 1.0e5 rxn_system = SurfaceReactor( T, P_initial, @@ -368,7 +368,7 @@ def test_solve_h2_coverage_dependence(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]): @@ -386,8 +386,8 @@ def test_solve_h2_coverage_dependence(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) From 4670244b2e9023cd2c04abd74a29043e692a2856 Mon Sep 17 00:00:00 2001 From: Richard West Date: Mon, 12 Jul 2021 23:46:19 -0400 Subject: [PATCH 38/43] CovDep Test: make the surface reactor test slower. Add a bunch of redundant species and reactions, to slow it down a bit, for benchmarking purposes. Simulation took 3.550e+00 seconds in rmgpy.solver.surfaceTest.SurfaceReactorCheck.test_solve_h2_coverage_dependence --- rmgpy/solver/surfaceTest.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/rmgpy/solver/surfaceTest.py b/rmgpy/solver/surfaceTest.py index 4c72f69be0..53dc09050c 100644 --- a/rmgpy/solver/surfaceTest.py +++ b/rmgpy/solver/surfaceTest.py @@ -317,11 +317,26 @@ def test_solve_h2_coverage_dependence(self): 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( From cf8146b9a208096d6029c2c21ae817fc1cce0080 Mon Sep 17 00:00:00 2001 From: Richard West Date: Mon, 12 Jul 2021 23:37:14 -0400 Subject: [PATCH 39/43] Refactor the coverage dependence calculations in surface solver. Main goal was to optimize everything happening inside the loops in the "residual" function (which gets called a lot). Simulation took 4.541e-01 seconds in rmgpy.solver.surfaceTest.SurfaceReactorCheck.test_solve_h2_coverage_dependence (previous run was 3.550e+00, so we're about 8x faster) --- rmgpy/solver/surface.pyx | 95 +++++++++++++++++++++------------------- 1 file changed, 49 insertions(+), 46 deletions(-) diff --git a/rmgpy/solver/surface.pyx b/rmgpy/solver/surface.pyx index 440ec6f44d..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 @@ -68,9 +68,8 @@ cdef class SurfaceReactor(ReactionSystem): cdef public np.ndarray species_on_surface # (catalyst surface, not core/edge surface) cdef public bint coverage_dependence - cdef public dict cov_dep - cdef public dict cov_dep_index_species - cdef public dict current_surface_coverages + cdef public dict coverage_dependencies + def __init__(self, T, @@ -109,9 +108,7 @@ cdef class SurfaceReactor(ReactionSystem): self.sens_conditions = sens_conditions self.n_sims = n_sims - self.cov_dep = {} - self.cov_dep_index_species = {} - self.current_surface_coverages = {} + self.coverage_dependencies = {} def convert_initial_keys_to_species_objects(self, species_dict): """ @@ -172,20 +169,34 @@ cdef class SurfaceReactor(ReactionSystem): #: 1 if it's on a surface, 0 if it's in the gas phase reactions_on_surface = np.zeros((self.num_core_reactions + self.num_edge_reactions), np.int) species_on_surface = np.zeros((self.num_core_species), np.int) - cov_dep_index_species = {} for spec, index in self.species_index.items(): if index >= self.num_core_species: continue if spec.contains_surface_site(): species_on_surface[index] = 1 - # map for species index to species for coverage dependence - cov_dep_index_species[index] = spec - 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 - self.cov_dep_index_species = cov_dep_index_species # Set initial conditions self.set_initial_conditions() @@ -232,10 +243,6 @@ cdef class SurfaceReactor(ReactionSystem): rxn.get_surface_rate_coefficient(self.T.value_si, self.surface_site_density.value_si )) - - if rxn.kinetics.coverage_dependence: - self.cov_dep[j] = rxn.kinetics.coverage_dependence # dict - else: if not warned and rxn.kinetics.is_pressure_dependent(): logging.warning("Pressure may be varying, but using initial pressure to evaluate k(T,P) expressions!") @@ -313,7 +320,6 @@ cdef class SurfaceReactor(ReactionSystem): surface_volume_ratio_si = self.surface_volume_ratio.value_si # 1/m total_surface_sites = V * surface_volume_ratio_si * self.surface_site_density.value_si # total surface sites in reactor - coverage_dependence = self.coverage_dependence for spec, coverage in self.initial_surface_coverages.items(): i = self.get_species_index(spec) @@ -362,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 @@ -370,8 +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 bint coverage_dependence - cdef dict cov_dep, cov_dep_index_species + cdef list list_of_coverage_deps + cdef double surface_site_fraction, total_sites, a, m, E + + ir = self.reactant_indices ip = self.product_indices @@ -379,9 +387,6 @@ cdef class SurfaceReactor(ReactionSystem): kf = self.kf # are already 'per m3 of reactor' even for surface reactions kr = self.kb # are already 'per m3 of reactor' even for surface reactions - coverage_dependence = self.coverage_dependence # turn coverage dependence on/off - cov_dep = self.cov_dep # load in coverage dependent parameters - cov_dep_index_species = self.cov_dep_index_species # load in species index to species inet = self.network_indices knet = self.network_leak_coefficients @@ -410,29 +415,36 @@ cdef class SurfaceReactor(ReactionSystem): A = self.V * surface_volume_ratio_si # area total_sites = self.surface_site_density.value_si * A # todo: double check units - current_surface_coverages = {} for j in range(num_core_species): if species_on_surface[j]: C[j] = (N[j] / V) / surface_volume_ratio_si - surface_site_fraction = N[j] / total_sites - species = cov_dep_index_species[j] - current_surface_coverages[species] = surface_site_fraction else: C[j] = N[j] / V #: 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 coverage_dependence: - if j in cov_dep: - cov_dep_j = cov_dep[j] - for species, parameters in cov_dep_j.items(): - if species in current_surface_coverages: - theta = current_surface_coverages[species] - if theta >= 1e-15: - k *= 10. ** (parameters['a'].value_si * theta) * theta ** parameters['m'].value_si *\ - np.exp(-1 * parameters['E'].value_si * theta / (constants.R * self.T.value_si)) if ir[j, 0] >= num_core_species or ir[j, 1] >= num_core_species or ir[j, 2] >= num_core_species: reaction_rate = 0.0 elif ir[j, 1] == -1: # only one reactant @@ -442,15 +454,6 @@ cdef class SurfaceReactor(ReactionSystem): else: # three reactants!! (really?) reaction_rate = k * C[ir[j, 0]] * C[ir[j, 1]] * C[ir[j, 2]] k = kr[j] - if coverage_dependence: - if j in cov_dep: - cov_dep_j = cov_dep[j] - for species, parameters in cov_dep_j.items(): - if species in current_surface_coverages: - theta = current_surface_coverages[species] - if theta >= 1e-15: - k *= 10. ** (parameters['a'].value_si * theta) * theta ** parameters['m'].value_si *\ - np.exp(-1 * parameters['E'].value_si * theta / (constants.R * self.T.value_si)) if ip[j, 0] >= num_core_species or ip[j, 1] >= num_core_species or ip[j, 2] >= num_core_species: pass elif ip[j, 1] == -1: # only one reactant From fa20deaadbf83402921df790d17627a159e09ab7 Mon Sep 17 00:00:00 2001 From: Richard West Date: Mon, 12 Jul 2021 23:54:58 -0400 Subject: [PATCH 40/43] CovDep test: stop plotting the graphs. These are helpful for manual inspection and helped me confirm it's working, but they're not needed for the automated test suite or continuous integration. --- rmgpy/solver/surfaceTest.py | 66 ++++++++++++++++++------------------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/rmgpy/solver/surfaceTest.py b/rmgpy/solver/surfaceTest.py index 53dc09050c..73a0d00cc4 100644 --- a/rmgpy/solver/surfaceTest.py +++ b/rmgpy/solver/surfaceTest.py @@ -132,23 +132,23 @@ def test_solve_h2(self): # Check that we've reached equilibrium self.assertAlmostEqual(reaction_rates[-1, 0], 0.0, delta=1e-2) - # Visualize the simulation results - import pylab - fig = pylab.figure(figsize=(6, 6)) - pylab.subplot(2, 1, 1) - 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) - pylab.legend(['H2', 'X', 'HX'], loc=4) - pylab.xlabel('Time (s)') - pylab.ylabel('Rate (mol/m$^\\mathdefault{3 or 2}$*s)') - # fig.subplots_adjust(left=0.21, bottom=0.10, right=0.95, top=0.95, wspace=0.20, hspace=0.35) - pylab.tight_layout() - # pylab.show() - pylab.savefig('surfaceTestH2.pdf') - return + # # Visualize the simulation results + # import pylab + # fig = pylab.figure(figsize=(6, 6)) + # pylab.subplot(2, 1, 1) + # 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) + # pylab.legend(['H2', 'X', 'HX'], loc=4) + # pylab.xlabel('Time (s)') + # pylab.ylabel('Rate (mol/m$^\\mathdefault{3 or 2}$*s)') + # # fig.subplots_adjust(left=0.21, bottom=0.10, right=0.95, top=0.95, wspace=0.20, hspace=0.35) + # pylab.tight_layout() + # # pylab.show() + # pylab.savefig('surfaceTestH2.pdf') + # return def test_solve_ch3(self): """ @@ -397,22 +397,22 @@ def test_solve_h2_coverage_dependence(self): # Check that we've reached equilibrium self.assertAlmostEqual(reaction_rates[-1, 0], 0.0, delta=1e-2) - # Visualize the simulation results - import pylab - fig = pylab.figure(figsize=(6, 6)) - pylab.subplot(2, 1, 1) - 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) - pylab.legend(['H2', 'X', 'HX'], loc=4) - pylab.xlabel('Time (s)') - pylab.ylabel('Rate (mol/m$^\\mathdefault{3 or 2}$*s)') - # fig.subplots_adjust(left=0.21, bottom=0.10, right=0.95, top=0.95, wspace=0.20, hspace=0.35) - pylab.tight_layout() - # pylab.show() - pylab.savefig('surfaceTestH2covdep.pdf') + # # Visualize the simulation results + # import pylab + # fig = pylab.figure(figsize=(6, 6)) + # pylab.subplot(2, 1, 1) + # 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) + # pylab.legend(['H2', 'X', 'HX'], loc=4) + # pylab.xlabel('Time (s)') + # pylab.ylabel('Rate (mol/m$^\\mathdefault{3 or 2}$*s)') + # # fig.subplots_adjust(left=0.21, bottom=0.10, right=0.95, top=0.95, wspace=0.20, hspace=0.35) + # pylab.tight_layout() + # # pylab.show() + # pylab.savefig('surfaceTestH2covdep.pdf') def test_solve_ch3_coverage_dependence(self): """ From a00934b7e74914bc527241e01a630b593e74df20 Mon Sep 17 00:00:00 2001 From: Chris B <56306881+ChrisBNEU@users.noreply.github.com> Date: Wed, 21 Jul 2021 16:52:47 -0400 Subject: [PATCH 41/43] added check that covdep makes a difference changed CH3 coverage unittest to run an additional simulation without coverage dependence. the solution arrays from each are compared. --- rmgpy/solver/surfaceTest.py | 64 +++++++++++++++++++++++++++---------- 1 file changed, 47 insertions(+), 17 deletions(-) diff --git a/rmgpy/solver/surfaceTest.py b/rmgpy/solver/surfaceTest.py index 73a0d00cc4..8448255c9f 100644 --- a/rmgpy/solver/surfaceTest.py +++ b/rmgpy/solver/surfaceTest.py @@ -397,23 +397,6 @@ def test_solve_h2_coverage_dependence(self): # Check that we've reached equilibrium self.assertAlmostEqual(reaction_rates[-1, 0], 0.0, delta=1e-2) - # # Visualize the simulation results - # import pylab - # fig = pylab.figure(figsize=(6, 6)) - # pylab.subplot(2, 1, 1) - # 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) - # pylab.legend(['H2', 'X', 'HX'], loc=4) - # pylab.xlabel('Time (s)') - # pylab.ylabel('Rate (mol/m$^\\mathdefault{3 or 2}$*s)') - # # fig.subplots_adjust(left=0.21, bottom=0.10, right=0.95, top=0.95, wspace=0.20, hspace=0.35) - # pylab.tight_layout() - # # pylab.show() - # pylab.savefig('surfaceTestH2covdep.pdf') - def test_solve_ch3_coverage_dependence(self): """ Test the surface batch reactor can properly apply coverage dependent parameters @@ -554,3 +537,50 @@ def test_solve_ch3_coverage_dependence(self): # 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'), + coverage_dependence=False, + 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.array_equal(y,y_off)) + self.assertFalse(np.array_equal(species_rates, species_rates_off)) From 99ab2817edde23e9e73a38a9727f67691f6661c8 Mon Sep 17 00:00:00 2001 From: Chris B <56306881+ChrisBNEU@users.noreply.github.com> Date: Wed, 21 Jul 2021 18:01:16 -0400 Subject: [PATCH 42/43] updated documentation with coverage dependence added instructions on input file construction and updating the database. --- documentation/source/users/rmg/surfaces.rst | 55 +++++++++++++++++++++ 1 file changed, 55 insertions(+) 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 ========= From 5918981394fd0829447e15f78692a3648a6f9945 Mon Sep 17 00:00:00 2001 From: Chris B <56306881+ChrisBNEU@users.noreply.github.com> Date: Thu, 22 Jul 2021 15:01:59 -0400 Subject: [PATCH 43/43] updated covdep vs non covdep check to 'all close' --- rmgpy/solver/surfaceTest.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/rmgpy/solver/surfaceTest.py b/rmgpy/solver/surfaceTest.py index 8448255c9f..1cc6452d21 100644 --- a/rmgpy/solver/surfaceTest.py +++ b/rmgpy/solver/surfaceTest.py @@ -546,7 +546,6 @@ def test_solve_ch3_coverage_dependence(self): initial_surface_coverages={x: 1.0}, surface_volume_ratio=(1., 'm^-1'), surface_site_density=(2.72e-9, 'mol/cm^2'), - coverage_dependence=False, termination=[]) rxn_system.initialize_model(core_species, core_reactions, edge_species, edge_reactions) @@ -582,5 +581,5 @@ def test_solve_ch3_coverage_dependence(self): self.assertAlmostEqual(species_rates_off[-1, 0], 0.0, delta=1e-2) # Check that coverages are different - self.assertFalse(np.array_equal(y,y_off)) - self.assertFalse(np.array_equal(species_rates, species_rates_off)) + self.assertFalse(np.allclose(y,y_off)) + self.assertFalse(np.allclose(species_rates, species_rates_off))