Adding capability to have reverse surfaceArrhenius training reactions

rmgpy/data/kinetics/family.py

@@ -1232,12 +1232,25 @@ def add_rules_from_training(self, thermo_database=None, train_indices=None):
             item.template = self.get_reaction_template_labels(item)
             new_degeneracy = self.calculate_degeneracy(item)
 
+            if isinstance(entry.data, SurfaceArrhenius):
+                data = SurfaceArrheniusBEP(
+                    #  analogous to Arrhenius.to_arrhenius_ep
+                    A=deepcopy(data.A),
+                    n=deepcopy(data.n),
+                    alpha=0,
+                    E0=deepcopy(data.Ea),
+                    Tmin=deepcopy(data.Tmin),
+                    Tmax=deepcopy(data.Tmax)
+                )
+            else:
+                data = data.to_arrhenius_ep()
+
             new_entry = Entry(
                 index=index,
                 label=';'.join([g.label for g in template]),
                 item=Reaction(reactants=[g.item for g in template],
                               products=[]),
-                data=data.to_arrhenius_ep(),
+                data=data,
                 rank=entry.rank,
                 reference=entry.reference,
                 short_desc="Rate rule generated from training reaction {0}. ".format(entry.index) + entry.short_desc,
@@ -1946,7 +1959,9 @@ def _generate_reactions(self, reactants, products=None, forward=True, prod_reson
 
             # ToDo: try to remove this hard-coding of reaction family name..
             if 'adsorption' in self.label.lower() and forward:
-                if molecules_a[0].contains_surface_site() and molecules_b[0].contains_surface_site():
+                if 'Surface_Dual_Adsorption_vdW' is self.label and forward:
+                    pass
+                elif molecules_a[0].contains_surface_site() and molecules_b[0].contains_surface_site():
                     # Can't adsorb something that's already adsorbed.
                     # Both reactants either contain or are a surface site.
                     return []

rmgpy/reaction.pxd

@@ -31,6 +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
 
 cimport numpy as np
 
@@ -47,6 +48,7 @@ cdef class Reaction:
     cdef public TransitionState transition_state
     cdef public KineticsModel kinetics
     cdef public Arrhenius network_kinetics
+    cdef public SurfaceArrhenius
     cdef public bint duplicate
     cdef public float _degeneracy
     cdef public list pairs
@@ -101,6 +103,8 @@ cdef class Reaction:
 
     cpdef reverse_arrhenius_rate(self, Arrhenius k_forward, str reverse_units, Tmin=?, Tmax=?)
 
+    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 np.ndarray calculate_tst_rate_coefficients(self, np.ndarray Tlist)

rmgpy/reaction.py

@@ -53,8 +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, SurfaceArrhenius, SurfaceArrheniusBEP, StickingCoefficientBEP
+    StickingCoefficient, 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.diffusionLimited import diffusion_limiter
 from rmgpy.molecule.element import Element, element_list
 from rmgpy.molecule.molecule import Molecule, Atom
@@ -787,6 +788,29 @@ def reverse_arrhenius_rate(self, k_forward, reverse_units, Tmin=None, Tmax=None)
         kr.fit_to_data(Tlist, klist, reverse_units, kf.T0.value_si)
         return kr
 
+    def reverse_surface_arrhenius_rate(self, k_forward, reverse_units, Tmin=None, Tmax=None):
+        """
+        Reverses the given k_forward, which must be a SurfaceArrhenius 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=SurfaceArrhenius, kr=SurfaceArrhenius)
+        cython.declare(Tlist=np.ndarray, klist=np.ndarray, i=cython.int)
+        kf = k_forward
+        if not isinstance(kf, SurfaceArrhenius): # Only reverse SurfaceArrhenius rates
+            raise TypeError(f'Expected a SurfaceArrhenius 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] = kf.get_rate_coefficient(Tlist[i]) / self.get_equilibrium_constant(Tlist[i])
+        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):
         """
         Generate and return a rate coefficient model for the reverse reaction. 
@@ -800,6 +824,7 @@ def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, T
         supported_types = (
             KineticsData.__name__,
             Arrhenius.__name__,
+            SurfaceArrhenius.__name__,
             MultiArrhenius.__name__,
             PDepArrhenius.__name__,
             MultiPDepArrhenius.__name__,
@@ -825,7 +850,10 @@ def generate_reverse_rate_coefficient(self, network_kinetics=False, Tmin=None, T
             return kr
 
         elif isinstance(kf, Arrhenius):
-            return self.reverse_arrhenius_rate(kf, kunits, Tmin, Tmax)
+            if isinstance(kf, SurfaceArrhenius):
+                return self.reverse_surface_arrhenius_rate(kf, kunits, Tmin, Tmax)
+            else:
+                return self.reverse_arrhenius_rate(kf, kunits, Tmin, Tmax)
 
         elif network_kinetics and self.network_kinetics is not None:
             kf = self.network_kinetics

rmgpy/reactionTest.py

@@ -696,6 +696,35 @@ def test_generate_reverse_rate_coefficient_arrhenius(self):
             krevrev = reverse_reverse_kinetics.get_rate_coefficient(T, P)
             self.assertAlmostEqual(korig / krevrev, 1.0, 0)
 
+    def test_reverse_surface_arrhenius_rate(self):
+        """
+        Test the Reaction.reverse_surface_arrhenius_rate() method works for SurfaceArrhenius format.
+        """
+        original_kinetics = SurfaceArrhenius(
+            A=(1.195e12, 'm^2/(mol*s)'),
+            n=0.0,
+            Ea=(14.989, 'kcal/mol'),
+            T0=(1, 'K'),
+            Tmin=(300, 'K'),
+            Tmax=(2000, 'K'),
+        )
+        self.reaction2.kinetics = original_kinetics
+
+        reverse_kinetics = self.reaction2.generate_reverse_rate_coefficient()
+
+        self.reaction2.kinetics = reverse_kinetics
+        # reverse reactants, products to ensure Keq is correctly computed
+        self.reaction2.reactants, self.reaction2.products = self.reaction2.products, self.reaction2.reactants
+        reverse_reverse_kinetics = self.reaction2.generate_reverse_rate_coefficient()
+
+        # 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 = 1e5
+        for T in Tlist:
+            korig = original_kinetics.get_rate_coefficient(T, P)
+            krevrev = reverse_reverse_kinetics.get_rate_coefficient(T, P)
+            self.assertAlmostEqual(korig / krevrev, 1.0, 0)
+
     @work_in_progress
     def test_generate_reverse_rate_coefficient_arrhenius_ep(self):
         """

testing/databaseTest.py

@@ -552,6 +552,7 @@ def kinetics_check_adjlists_nonidentical(self, database):
                     species_dict[product.label] = product
 
         tst = []
+        boo = False
         # Go through all species to make sure they are nonidentical
         species_list = list(species_dict.values())
         labeled_atoms = [species.molecule[0].get_all_labeled_atoms() for species in species_list]
@@ -566,16 +567,15 @@ def kinetics_check_adjlists_nonidentical(self, database):
                 except KeyError:
                     # atom labels did not match, therefore not a match
                     continue
-                tst.append((species_list[i].molecule[0].is_isomorphic(species_list[j].molecule[0], initial_map),
-                            "Species {0} and species {1} in {2} database were found to be identical.".format(
-                                species_list[i].label, species_list[j].label, database.label)))
-
-        boo = False
-        for i in range(len(tst)):
-            if tst[i][0]:
-                logging.error(tst[i][1])
-                boo = True
-
+                m1 = species_list[i].molecule[0]
+                m2 = species_list[j].molecule[0]
+                if not m1.is_mapping_valid(m2, initial_map, equivalent=True):
+                    # the mapping is invalid so they're not isomorphic
+                    continue
+                if m1.is_isomorphic(m2, initial_map):
+                    logging.error("Species {0} and species {1} in {2} database were found to be identical.".format(
+                                species_list[i].label, species_list[j].label, database.label))
+                    boo = True
         if boo:
             raise ValueError("Error occured in databaseTest. Please check log warnings for all error messages.")
 
@@ -609,7 +609,8 @@ def kinetics_check_rate_units_are_correct(self, database, tag='library'):
                     a_factor = k.A
                     expected = copy(dimensionalities[molecularity])
                     # for each surface reactant but one, switch from (m3/mol) to (m2/mol)
-                    expected[pq.m] -= (surface_reactants - 1)
+                    for _ in range(surface_reactants - 1):
+                        expected[pq.m] -= 1
                     if pq.Quantity(1.0, a_factor.units).simplified.dimensionality != expected:
                         boo = True
                         logging.error('Reaction {0} from {1} {2}, has invalid units {3}'.format(