From fec029fe75b9d366560eea823d0cfac765db729b Mon Sep 17 00:00:00 2001 From: davidfarinajr Date: Thu, 18 Feb 2021 19:11:31 -0500 Subject: [PATCH] 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 | 39 ++++++++++++++++++++++++++++++++++++--- 2 files changed, 40 insertions(+), 5 deletions(-) diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd index 8f866751dfa..f3b583f4b03 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 7171d1a1a81..18645a8d929 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,7 +854,32 @@ 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). + """ + cython.declare(kf=StickingCoefficient, kr=SurfaceArrhenius) + cython.declare(Tlist=np.ndarray, klist=np.ndarray, i=cython.int) + kf = k_forward + if not isinstance(kf, StickingCoefficient): # Only reverse StickingCoefficient rates + raise TypeError(f'Expected a StickingCoefficient object for k_forward but received {kf}') + 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 @@ -875,6 +900,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 +933,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)