diff --git a/rmgpy/data/kinetics/rules.py b/rmgpy/data/kinetics/rules.py index 5a9f05bc4a..9bfbd49e0b 100644 --- a/rmgpy/data/kinetics/rules.py +++ b/rmgpy/data/kinetics/rules.py @@ -561,12 +561,19 @@ def _get_average_kinetics(self, kinetics_list): Hence we average n, Ea, and alpha arithmetically, but we average log A (geometric average) """ + + kinetics_type = None logA = 0.0 n = 0.0 E0 = 0.0 alpha = 0.0 count = len(kinetics_list) for kinetics in kinetics_list: + if kinetics_type is None: + kinetics_type = type(kinetics) + else: + if type(kinetics) != kinetics_type: + raise KineticsError(f"Unable to average kinetics with mixed kinetics types ({kinetics_type} != {type(kinetics)})") logA += math.log10(kinetics.A.value_si) n += kinetics.n.value_si alpha += kinetics.alpha.value_si diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd index 586c4a6796..a437c82c89 100644 --- a/rmgpy/reaction.pxd +++ b/rmgpy/reaction.pxd @@ -99,6 +99,10 @@ cdef class Reaction: cpdef double get_surface_rate_coefficient(self, double T, double surface_site_density) except -2 + cpdef surface_arrhenius_to_sticking_coeff(self, double surface_site_density, Tmin=?, Tmax=?) + + cpdef sticking_coeff_to_surface_arrhenius(self, double surface_site_density, Tmin=?, Tmax=?) + cpdef fix_barrier_height(self, bint force_positive=?) cpdef reverse_arrhenius_rate(self, Arrhenius k_forward, str reverse_units, Tmin=?, Tmax=?) diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index dfb598784f..42a9a03d1a 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -729,6 +729,103 @@ def get_surface_rate_coefficient(self, T, surface_site_density): raise NotImplementedError("Can't get_surface_rate_coefficient for kinetics type {!r}".format(type(self.kinetics))) + def surface_arrhenius_to_sticking_coeff(self, surface_site_density, Tmin=None, Tmax=None): + """ + Converts `SurfaceArrhenius` kinetics to `StickingCoeff` kinetics using the provided + `surface_site_density` in SI units (mol/m^2). The reaction's kinetics type must but be + `SurfaceArrhenius`. + + Returns: + `StickingCoefficient` kinetics + """ + cython.declare(kf=StickingCoefficient) + cython.declare(sticking_coefficient=cython.double, molecular_weight_kg=cython.double) + cython.declare(Tlist=np.ndarray, klist=np.ndarray, i=cython.int, number_of_sites=cython.int) + if not isinstance(self.kinetics, SurfaceArrhenius): + raise TypeError(f'Expected a SurfaceArrhenius object but received {self.kinetics}') + # determine the adsorbate and number of surface sites + adsorbate = None + number_of_sites = 0 + for r in self.reactants: + if r.contains_surface_site(): + number_of_sites += 1 + else: + if adsorbate is None: + adsorbate = r + else: + logging.error("Error in kinetics for reaction {0!s}: " + "more than one adsorbate detected".format(self)) + raise ReactionError("More than one adsorbate detected") + + if adsorbate is None or adsorbate.contains_surface_site(): + logging.error("Problem reaction: {0!s}".format(self)) + raise ReactionError("Couldn't find the adsorbate!") + molecular_weight_kg = adsorbate.molecular_weight.value_si + + # generate temperature array to evaluate the rate coeff + 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 sticking coefficient at each temperature + klist = np.zeros_like(Tlist) + for i in range(len(Tlist)): + sticking_coefficient = self.kinetics.get_rate_coefficient(Tlist[i]) + sticking_coefficient /= math.sqrt(constants.kB * Tlist[i] / (2 * math.pi * molecular_weight_kg)) + klist[i] = sticking_coefficient + for i in range(number_of_sites): + klist *= surface_site_density + + # set the max sticking coeff to 1 + for i in np.argwhere(klist>1): + klist[i] = 1 + + # create Sticking Coeff kinetics and fit to sticking coeff array + kf = StickingCoefficient() + kf.fit_to_data(Tlist, klist, "", kf.T0.value_si) + return kf + + def sticking_coeff_to_surface_arrhenius(self, surface_site_density, Tmin=None, Tmax=None): + """ + Converts `StickingCoefficient` kinetics to `SurfaceArrhenius` kinetics using the provided + `surface_site_density` in SI units (mol/m^2). The reaction's kinetics type must but be + `StickingCoefficent`. + + Returns: + `SurfaceArrhenius` kinetics + """ + cython.declare(kf=SurfaceArrhenius) + cython.declare(Tlist=np.ndarray, klist=np.ndarray, i=cython.int) + if not isinstance(self.kinetics, StickingCoefficient): + raise TypeError(f'Expected a Sticking Coeff object but received {self.kinetics}') + + # Get the units for the reverse rate coefficient + try: + surf_reacts = [spcs for spcs in self.reactants if spcs.contains_surface_site()] + except IndexError: + raise KineticsError("Species do not have an rmgpy.molecule.Molecule." + "Cannot determine units for rate coefficient.") + n_surf = len(surf_reacts) + n_gas = len(self.reactants) - len(surf_reacts) + kunits = get_rate_coefficient_units_from_reaction_order(n_gas, n_surf) + + # generate temperature array to evaluate the rate coeff + 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 rate coeff k_f(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) + + # create Surface Arrh kinetics and fit to rate coeff array + kf = SurfaceArrhenius() + kf.fit_to_data(Tlist, klist, kunits) + return kf + def fix_diffusion_limited_a_factor(self, T): """ Decrease the pre-exponential factor (A) by the diffusion factor diff --git a/rmgpy/reactionTest.py b/rmgpy/reactionTest.py index 64e8eb002a..58ee46f618 100644 --- a/rmgpy/reactionTest.py +++ b/rmgpy/reactionTest.py @@ -137,6 +137,7 @@ class TestSurfaceReaction(unittest.TestCase): """Test surface reactions""" def setUp(self): + self.surface_site_density = 2.483e-05 # mol/m^2 for Pt111 m_h2 = Molecule().from_smiles("[H][H]") m_x = Molecule().from_adjacency_list("1 X u0 p0") m_hx = Molecule().from_smiles("[H][*]") @@ -313,6 +314,44 @@ def test_reverse_sticking_coeff_rate(self): krevrev = rxn_copy.get_rate_coefficient(T, P, surface_site_density=2.5e-5) self.assertAlmostEqual(korig / krevrev, 1.0, 0) + def test_sticking_coefficient_to_surface_arrhenius(self): + """ + Tests that the sticking_coeff_to_surface_arrhenius() method is working properly + """ + arr = self.rxn2sSC.sticking_coeff_to_surface_arrhenius(surface_site_density=self.surface_site_density) + temps = numpy.linspace(200,2000,10) + for temp in temps: + arr_k = arr.get_rate_coefficient(temp) + sc_k = self.rxn2sSC.get_rate_coefficient(temp,surface_site_density=self.surface_site_density) + self.assertAlmostEqual(arr_k,sc_k,6) + + def test_surface_arrhenius_to_sticking_coeff(self): + """ + Tests that the surface_arrhenius_to_sticking_coeff() method is working properly + """ + + sc = self.rxn2sSA.surface_arrhenius_to_sticking_coeff(surface_site_density=self.surface_site_density) + arr = self.rxn2sSA.kinetics + self.rxn2sSA.kinetics = sc + temps = numpy.linspace(200,2000,10) + for temp in temps: + arr_k = arr.get_rate_coefficient(temp) + sc_k = self.rxn2sSA.get_rate_coefficient(temp,surface_site_density=self.surface_site_density) + self.assertAlmostEqual(arr_k,sc_k,6) + self.rxn2sSA.kinetics = arr + + + # also test that the max sticking coefficient is 1 if we are converting a super fast Surface Arr rate + arr_too_fast = SurfaceArrhenius(A=(2.7e30, 'cm^3/(mol*s)'), comment="""a super fast rate""") + self.rxn2sSA.kinetics = arr_too_fast + sc = self.rxn2sSA.surface_arrhenius_to_sticking_coeff(surface_site_density=self.surface_site_density) + temps = numpy.linspace(200,2000,10) + for temp in temps: + sc_T = sc.get_sticking_coefficient(temp) + self.assertAlmostEqual(sc_T,1.0,2) + + self.rxn2sSA.kinetics = arr + class TestReaction(unittest.TestCase): """ Contains unit tests of the Reaction class.