Merge pull request #2080 to enable calculation of reverse rate for St…

…ickingCoefficient rates.

Merge pull request #2080 from ReactionMechanismGenerator/reverse_sticking_coeff
Added StickingCoefficient kinetics handling to reaction get_rate_coefficient method so that we can use this method for all types of kinetics.
Added a method to reverse StickingCoefficient kinetics with default surface site density.

Will be useful for rendering results on the website, and for reversing training reactions when
they are published in the reverse direction.

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
 
@@ -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
 
@@ -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
@@ -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:
@@ -845,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,
@@ -866,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
@@ -898,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)

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.