diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 17edff2a89..0d4b1701ac 100644 --- a/rmgpy/data/kinetics/family.py +++ b/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 [] diff --git a/rmgpy/reaction.pxd b/rmgpy/reaction.pxd index 5c5fe76cfd..a266d0756a 100644 --- a/rmgpy/reaction.pxd +++ b/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) diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index 9b3dc0ae98..ebe9ea22a2 100644 --- a/rmgpy/reaction.py +++ b/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 diff --git a/rmgpy/reactionTest.py b/rmgpy/reactionTest.py index 56903b55e2..2ab07830c2 100644 --- a/rmgpy/reactionTest.py +++ b/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): """ diff --git a/testing/databaseTest.py b/testing/databaseTest.py index d36c2cf000..ed0a65393d 100644 --- a/testing/databaseTest.py +++ b/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(