SurfaceArrhenius <=> StickingCoefficient

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

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=?)

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

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.