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

@@ -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)
 

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)