From 30bd07a2265c354c6e7a78262f58e5e25da01fa2 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 18 Feb 2021 19:05:30 -0500 Subject: [PATCH 1/3] added StickingCoeff to `get_rate_coefficient` reaction method So that we can call this method regardless of what type the kinetics is --- rmgpy/reaction.pxd | 2 +- rmgpy/reaction.py | 13 +++++++++++-- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd index 8a1cb6c593..8f866751df 100644 --- a/rmgpy/reaction.pxd +++ b/rmgpy/reaction.pxd @@ -95,7 +95,7 @@ cdef class Reaction: cpdef int get_stoichiometric_coefficient(self, Species spec) - cpdef double get_rate_coefficient(self, double T, double P=?) + cpdef double get_rate_coefficient(self, double T, double P=?, double surface_site_density=?) cpdef double get_surface_rate_coefficient(self, double T, double surface_site_density) except -2 diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index 8f6766ff1a..7171d1a1a8 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -657,7 +657,7 @@ def get_stoichiometric_coefficient(self, spec): if product is spec: stoich += 1 return stoich - def get_rate_coefficient(self, T, P=0): + def get_rate_coefficient(self, T, P=0, surface_site_density=0): """ Return the overall rate coefficient for the forward reaction at temperature `T` in K and pressure `P` in Pa, including any reaction @@ -666,8 +666,17 @@ def get_rate_coefficient(self, T, P=0): If diffusion_limiter is enabled, the reaction is in the liquid phase and we use a diffusion limitation to correct the rate. If not, then use the intrinsic rate coefficient. + + If the reaction has sticking coefficient kinetics, a nonzero surface site density + in `mol/m^2` must be provided """ - if diffusion_limiter.enabled: + if isinstance(self.kinetics, StickingCoefficient): + if surface_site_density <= 0: + raise ValueError("Please provide a postive surface site density in mol/m^2 " + f"for calculating the rate coefficient of {StickingCoefficient.__name__} kinetics") + else: + return self.get_surface_rate_coefficient(T, surface_site_density) + elif diffusion_limiter.enabled: try: k = self.k_effective_cache[T] except KeyError: From bcaca2f6a35db2eb7fcea6388f6d2fe11ce1112f Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 18 Feb 2021 19:11:31 -0500 Subject: [PATCH 2/3] added `reverse_sticking_coeff_rate` reaction method This method fits a SurfaceArrhenius model to the reverse rate by evaluating the forward stickingcoeff rate for a surface site density (default 2.5e-05 mol/m^2) over a temperature range --- rmgpy/reaction.pxd | 6 ++++-- rmgpy/reaction.py | 43 ++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 44 insertions(+), 5 deletions(-) diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd index 8f866751df..f3b583f4b0 100644 --- a/rmgpy/reaction.pxd +++ b/rmgpy/reaction.pxd @@ -31,7 +31,7 @@ from rmgpy.molecule.molecule cimport Atom, Molecule from rmgpy.molecule.element cimport Element from rmgpy.kinetics.model cimport KineticsModel from rmgpy.kinetics.arrhenius cimport Arrhenius -from rmgpy.kinetics.surface cimport SurfaceArrhenius +from rmgpy.kinetics.surface cimport SurfaceArrhenius, StickingCoefficient cimport numpy as np @@ -105,7 +105,9 @@ cdef class Reaction: cpdef reverse_surface_arrhenius_rate(self, SurfaceArrhenius k_forward, str reverse_units, Tmin=?, Tmax=?) - cpdef generate_reverse_rate_coefficient(self, bint network_kinetics=?, Tmin=?, Tmax=?) + cpdef reverse_sticking_coeff_rate(self, StickingCoefficient k_forward, str reverse_units, double surface_site_density, Tmin=?, Tmax=?) + + cpdef generate_reverse_rate_coefficient(self, bint network_kinetics=?, Tmin=?, Tmax=?, double surface_site_density=?) cpdef np.ndarray calculate_tst_rate_coefficients(self, np.ndarray Tlist) diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index 7171d1a1a8..aa1be07831 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -53,9 +53,9 @@ from rmgpy.exceptions import ReactionError, KineticsError from rmgpy.kinetics import KineticsData, ArrheniusBM, ArrheniusEP, ThirdBody, Lindemann, Troe, Chebyshev, \ PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, get_rate_coefficient_units_from_reaction_order, \ - StickingCoefficient, SurfaceArrheniusBEP, StickingCoefficientBEP + SurfaceArrheniusBEP, StickingCoefficientBEP from rmgpy.kinetics.arrhenius import Arrhenius # Separate because we cimport from rmgpy.kinetics.arrhenius -from rmgpy.kinetics.surface import SurfaceArrhenius # Separate because we cimport from rmgpy.kinetics.surface +from rmgpy.kinetics.surface import SurfaceArrhenius, StickingCoefficient # Separate because we cimport from rmgpy.kinetics.surface from rmgpy.kinetics.diffusionLimited import diffusion_limiter from rmgpy.molecule.element import Element, element_list from rmgpy.molecule.molecule import Molecule, Atom @@ -854,11 +854,40 @@ def reverse_surface_arrhenius_rate(self, k_forward, reverse_units, Tmin=None, Tm kr.fit_to_data(Tlist, klist, reverse_units, kf.T0.value_si) return kr - def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, Tmax=None): + def reverse_sticking_coeff_rate(self, k_forward, reverse_units, surface_site_density, Tmin=None, Tmax=None): + """ + Reverses the given k_forward, which must be a StickingCoefficient type. + You must supply the correct units for the reverse rate. + The equilibrium constant is evaluated from the current reaction instance (self). + The surface_site_density in `mol/m^2` is used to evalaute the forward rate constant. + """ + cython.declare(kf=StickingCoefficient, kr=SurfaceArrhenius) + cython.declare(Tlist=np.ndarray, klist=np.ndarray, i=cython.int) + if not isinstance(k_forward, StickingCoefficient): # Only reverse StickingCoefficient rates + raise TypeError(f'Expected a StickingCoefficient object for k_forward but received {k_forward}') + kf = k_forward + if Tmin is not None and Tmax is not None: + Tlist = 1.0 / np.linspace(1.0 / Tmax.value, 1.0 / Tmin.value, 50) + else: + Tlist = 1.0 / np.arange(0.0005, 0.0034, 0.0001) + # Determine the values of the reverse rate coefficient k_r(T) at each temperature + klist = np.zeros_like(Tlist) + for i in range(len(Tlist)): + klist[i] = \ + self.get_surface_rate_coefficient(Tlist[i], surface_site_density=surface_site_density) / \ + self.get_equilibrium_constant(Tlist[i], surface_site_density=surface_site_density) + kr = SurfaceArrhenius() + kr.fit_to_data(Tlist, klist, reverse_units, kf.T0.value_si) + return kr + + def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, Tmax=None, surface_site_density=0): """ Generate and return a rate coefficient model for the reverse reaction. Currently this only works if the `kinetics` attribute is one of several (but not necessarily all) kinetics types. + + If the reaction kinetics model is Sticking Coefficient, please provide a nonzero + surface site density in `mol/m^2` which is required to evaluate the rate coefficient. """ cython.declare(Tlist=np.ndarray, Plist=np.ndarray, K=np.ndarray, rxn=Reaction, klist=np.ndarray, i=cython.size_t, @@ -875,6 +904,7 @@ def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, T ThirdBody.__name__, Lindemann.__name__, Troe.__name__, + StickingCoefficient.__name__, ) # Get the units for the reverse rate coefficient @@ -907,6 +937,13 @@ def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, T else: return self.reverse_arrhenius_rate(kf, kunits, Tmin, Tmax) + elif isinstance(kf, StickingCoefficient): + if surface_site_density <= 0: + raise ValueError("Please provide a postive surface site density in mol/m^2 " + f"for calculating the rate coefficient of {StickingCoefficient.__name__} kinetics") + else: + return self.reverse_sticking_coeff_rate(kf, kunits, surface_site_density, Tmin, Tmax) + elif network_kinetics and self.network_kinetics is not None: kf = self.network_kinetics return self.reverse_arrhenius_rate(kf, kunits) From 7588b2a8e5fb9b92c1fe5a75897ae00e46441fe0 Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 18 Feb 2021 19:11:57 -0500 Subject: [PATCH 3/3] added test for `reverse_sticking_coeff_rate` method --- rmgpy/reactionTest.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/rmgpy/reactionTest.py b/rmgpy/reactionTest.py index 3d8b229a88..2d6f2dd742 100644 --- a/rmgpy/reactionTest.py +++ b/rmgpy/reactionTest.py @@ -36,6 +36,7 @@ import cantera as ct import numpy from external.wip import work_in_progress +from copy import deepcopy import rmgpy.constants as constants from rmgpy.kinetics import Arrhenius, ArrheniusEP, MultiArrhenius, PDepArrhenius, MultiPDepArrhenius, \ @@ -269,6 +270,29 @@ def test_equilibrium_constant_surface_kc(self): for i in range(len(Tlist)): self.assertAlmostEqual(Kclist[i] / Kclist0[i], 1.0, 4) + def test_reverse_sticking_coeff_rate(self): + """ + Test the Reaction.reverse_sticking_coeff_rate() method works for StickingCoefficient format. + """ + + original_kinetics = self.rxn2sSC.kinetics + rxn_copy = deepcopy(self.rxn2sSC) + reverse_kinetics = self.rxn2sSC.generate_reverse_rate_coefficient(surface_site_density=2.5e-5) + + rxn_copy.kinetics = reverse_kinetics + # reverse reactants, products to ensure Keq is correctly computed + rxn_copy.reactants, rxn_copy.products = rxn_copy.products, rxn_copy.reactants + reverse_reverse_kinetics = rxn_copy.generate_reverse_rate_coefficient(surface_site_density=2.5e-5) + rxn_copy.kinetics = reverse_reverse_kinetics + + # check that reverting the reverse yields the original + Tlist = numpy.arange(original_kinetics.Tmin.value_si, original_kinetics.Tmax.value_si, 200.0, numpy.float64) + P = 0 + for T in Tlist: + korig = self.rxn2sSC.get_rate_coefficient(T, P, surface_site_density=2.5e-5) + krevrev = rxn_copy.get_rate_coefficient(T, P, surface_site_density=2.5e-5) + self.assertAlmostEqual(korig / krevrev, 1.0, 0) + class TestReaction(unittest.TestCase): """ Contains unit tests of the Reaction class.