From aa882682c34461b45ce8dbe276d0aa4160935d7a Mon Sep 17 00:00:00 2001 From: Emily Mazeau Date: Wed, 5 Feb 2020 18:34:22 -0500 Subject: [PATCH 1/6] adding in a way to get reverse rates for surface arrhenius reactions --- rmgpy/reaction.pxd | 4 ++++ rmgpy/reaction.py | 32 ++++++++++++++++++++++++++++++-- 2 files changed, 34 insertions(+), 2 deletions(-) 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 From 6854b13cca3ad298ebbab01d7b5a03e2bb127839 Mon Sep 17 00:00:00 2001 From: Emily Mazeau Date: Wed, 5 Feb 2020 18:48:53 -0500 Subject: [PATCH 2/6] don't overwrite surface arrhenius data with arrheniusEP --- rmgpy/data/kinetics/family.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 17edff2a89..66db046ffa 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, From e8e292875bb412c8690b0a0584d2dc31c700d5b9 Mon Sep 17 00:00:00 2001 From: Emily Mazeau Date: Thu, 6 Feb 2020 22:41:03 -0500 Subject: [PATCH 3/6] Adding a hard-coded exception for Surface Dual Adsorption vdW --- rmgpy/data/kinetics/family.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 66db046ffa..0d4b1701ac 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1959,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 [] From eb65ec0ba9a39e480828038126e683d2dbfc9e63 Mon Sep 17 00:00:00 2001 From: Emily Mazeau Date: Tue, 18 Feb 2020 17:57:34 -0500 Subject: [PATCH 4/6] Fixing a kinetic units test in databaseTest Surface families that have only one surface species and nothing else, start with units of inverse seconds, and so you can't decrease the power of metres by zero, because there is no power of meters to decrease. Instead of subtracting zero from a dictionary value which didn't exist, we now subtract one from it zero times, thus avoiding the key error! --- testing/databaseTest.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/testing/databaseTest.py b/testing/databaseTest.py index d36c2cf000..714a69039c 100644 --- a/testing/databaseTest.py +++ b/testing/databaseTest.py @@ -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( From 8c78e91e558a03078ef6a071851f223cb2e36a2b Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 20 Feb 2020 14:14:42 -0500 Subject: [PATCH 5/6] Fix isomorphism checking in databaseTest. Previously, if there was an non-isomorphism present in the labeled atoms, then the mapping would be invalid, but this is not (currently) checked by the is_isomorphic method, so two non-isomorphic molecules would be found isomorphic, and would fail the test. Now we ensure that the mapping is valid, before passing it to is_isomorphic. This should work whatever the outcome of the discussion around PR #1894 --- testing/databaseTest.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/testing/databaseTest.py b/testing/databaseTest.py index 714a69039c..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.") From d87147c7040eaed70f3a15f5b03bd598ddda709e Mon Sep 17 00:00:00 2001 From: Emily Mazeau Date: Thu, 20 Feb 2020 15:23:16 -0500 Subject: [PATCH 6/6] adding a test to check that we can generate reverse rate coefficients for surface arrhenius --- rmgpy/reactionTest.py | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) 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): """