From be945f2a9fdae582173b57a2bb228121465c43f7 Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Wed, 11 Mar 2026 08:21:03 -0400 Subject: [PATCH 1/5] Delete redundant surfaceArrhenius initialization Surface arrhenius inherits from regular arrhenius, so it's actually handled further up in the code (line 331) and this should be deleted to prevent confusion --- rmgpy/reaction.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index a8ae13f3d3a..88c6c1361a6 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -392,13 +392,6 @@ def to_cantera(self, species_list=None, use_chemkin_identifier=False): rate=ct.LindemannRate() ) - elif isinstance(self.kinetics, SurfaceArrhenius): - ct_reaction = ct.Reaction( - reactants=ct_reactants, - products=ct_products, - rate=ct.InterfaceArrheniusRate() - ) - elif isinstance(self.kinetics, StickingCoefficient): ct_reaction = ct.Reaction( reactants=ct_reactants, From 7f176924d7d2dc5881844027a51e943b5cc7fe3d Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Wed, 11 Mar 2026 09:59:46 -0400 Subject: [PATCH 2/5] Add SurfaceArrheniusBM class Add SurfaceArrheniusBM class to correspond to Cantera's InterfaceBlowersMaselRate. This also implements to_cantera for the gas-phase ArrheniusBM --- rmgpy/kinetics/surface.pxd | 7 +- rmgpy/kinetics/surface.pyx | 166 +++++++++++++++++++++++++++++++++++++ rmgpy/reaction.py | 17 +++- 3 files changed, 188 insertions(+), 2 deletions(-) diff --git a/rmgpy/kinetics/surface.pxd b/rmgpy/kinetics/surface.pxd index 0a834ea5a93..b1b4376d9fa 100644 --- a/rmgpy/kinetics/surface.pxd +++ b/rmgpy/kinetics/surface.pxd @@ -28,7 +28,7 @@ cimport numpy as np from rmgpy.kinetics.model cimport KineticsModel -from rmgpy.kinetics.arrhenius cimport Arrhenius, ArrheniusEP +from rmgpy.kinetics.arrhenius cimport Arrhenius, ArrheniusEP, ArrheniusBM from rmgpy.quantity cimport ScalarQuantity, ArrayQuantity ################################################################################ @@ -80,6 +80,11 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): cdef public dict _coverage_dependence pass +################################################################################ +cdef class SurfaceArrheniusBM(ArrheniusBM): + cdef public dict _coverage_dependence + pass + ################################################################################ cdef class SurfaceChargeTransfer(KineticsModel): diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx index 7787a82e479..abfad22db63 100644 --- a/rmgpy/kinetics/surface.pyx +++ b/rmgpy/kinetics/surface.pyx @@ -785,6 +785,172 @@ cdef class SurfaceArrheniusBEP(ArrheniusEP): ################################################################################ +cdef class SurfaceArrheniusBM(ArrheniusBM): + """ + A kinetics model based on the modified Arrhenius equation, using the + Blowers-Masel equation to determine the activation energy. + Based on Blowers and Masel's 2000 paper Engineering Approximations for Activation + Energies in Hydrogen Transfer Reactions. + + It is very similar to the gas-phase :class:`ArrheniusBM`. + The only differences being the A factor has different units, + (and the catalysis community prefers to call it BEP rather than EP!) + and has a coverage_dependence parameter for coverage dependence + + The attributes are: + + ======================= ============================================================= + Attribute Description + ======================= ============================================================= + `A` The preexponential factor + `n` The temperature exponent + `w0` The average of the bond dissociation energies of the bond formed and the bond broken + `E0` The activation energy for a thermoneutral reaction + `Tmin` The minimum temperature at which the model is valid, or zero if unknown or undefined + `Tmax` The maximum temperature at which the model is valid, or zero if unknown or undefined + `Pmin` The minimum pressure at which the model is valid, or zero if unknown or undefined + `Pmax` The maximum pressure at which the model is valid, or zero if unknown or undefined + `coverage_dependence` A dictionary of coverage dependent parameters to a certain surface species with: + `a`, the coefficient for exponential dependence on the coverage, + `m`, the power-law exponent of coverage dependence, and + `E`, the activation energy dependence on coverage. + `uncertainty` Uncertainty information + `comment` Information about the model (e.g. its source) + ======================= ============================================================= + """ + + def __init__(self, A=None, n=0.0, w0=(0.0, 'J/mol'), E0=None, Tmin=None, Tmax=None, Pmin=None, Pmax=None, + coverage_dependence=None, uncertainty=None, comment=''): + KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, uncertainty=uncertainty, + comment=comment) + self.A = A + self.n = n + self.w0 = w0 + self.E0 = E0 + self.coverage_dependence = coverage_dependence + + property A: + """The preexponential factor, which has different usints from ArrheniusBM class.""" + def __get__(self): + return self._A + def __set__(self, value): + self._A = quantity.SurfaceRateCoefficient(value) + + property coverage_dependence: + """The coverage dependence parameters.""" + def __get__(self): + return self._coverage_dependence + def __set__(self, value): + self._coverage_dependence = {} + if value: + for species, parameters in value.items(): + processed_parameters = {'E': quantity.Energy(parameters['E']), + 'm': quantity.Dimensionless(parameters['m']), + 'a': quantity.Dimensionless(parameters['a'])} + self._coverage_dependence[species] = processed_parameters + + def __repr__(self): + """ + Return a string representation that can be used to reconstruct the + SurfaceArrheniusBM object. + """ + string = 'SurfaceArrheniusBM(A={0!r}, n={1!r}, w={2!r}, E0={3!r}'.format(self.A, self.n, self.w0, + self.E0) + if self.Tmin is not None: string += ', Tmin={0!r}'.format(self.Tmin) + if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) + if self.coverage_dependence: + string += ", coverage_dependence={" + for species, parameters in self.coverage_dependence.items(): + string += f"{species.to_chemkin()!r}: {{'a':{repr(parameters['a'])}, 'm':{repr(parameters['m'])}, 'E':{repr(parameters['E'])}}}," + string += "}" + if self.Pmin is not None: string += ', Pmin={0!r}'.format(self.Pmin) + if self.Pmax is not None: string += ', Pmax={0!r}'.format(self.Pmax) + if self.uncertainty is not None: string += ', uncertainty={0!r}'.format(self.uncertainty) + if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) + string += ')' + return string + + def __reduce__(self): + """ + A helper function used when pickling an SurfaceArrheniusBM object. + """ + return (SurfaceArrheniusBM, (self.A, self.n, self.w0, self.E0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, + self.uncertainty, self.coverage_dependence, self.comment)) + + cpdef SurfaceArrhenius to_arrhenius(self, double dHrxn): + """ + Return an :class:`SurfaceArrhenius` instance of the kinetics model using the + given enthalpy of reaction `dHrxn` to determine the activation energy. + + Note that despite its name it does not return a :class:`Arrhenius` object + (although :class:`SurfaceArrhenius` is a subclass of :class:`Arrhenius` + so in a way, it does). + """ + return SurfaceArrhenius( + A=self.A, + n=self.n, + Ea=(self.get_activation_energy(dHrxn) * 0.001, "kJ/mol"), + T0=(1, "K"), + Tmin=self.Tmin, + Tmax=self.Tmax, + uncertainty=self.uncertainty, + coverage_dependence=self.coverage_dependence, + comment=self.comment, + ) + + def to_cantera_kinetics(self): + """ + Converts the RMG SurfaceArrheniusBM object to a cantera InterfaceBlowersMaselRate + InterfaceBlowersMaselRate(A, b, E0 w0) where A is in units like m^2/kmol/s (depending on dimensionality) + b is dimensionless, E0 is in J/kmol, and w0 is in J/kmol + """ + import cantera as ct + + rate_units_conversion = {'1/s': 1, + 's^-1': 1, + 'm^2/(mol*s)': 1000, + 'm^4/(mol^2*s)': 1000000, + 'cm^2/(mol*s)': 1000, + 'cm^4/(mol^2*s)': 1000000, + 'm^2/(molecule*s)': 1000, + 'm^4/(molecule^2*s)': 1000000, + 'cm^2/(molecule*s)': 1000, + 'cm^4/(molecule^2*s)': 1000000, + 'cm^5/(mol^2*s)': 1000000, + 'm^5/(mol^2*s)': 1000000, + } + + if self._T0.value_si != 1: + A = self._A.value_si / (self._T0.value_si) ** self._n.value_si + else: + A = self._A.value_si + + try: + A *= rate_units_conversion[self._A.units] # convert from /mol to /kmol + except KeyError: + raise ValueError('Arrhenius A-factor units {0} not found among accepted units for converting to ' + 'Cantera Arrhenius object.'.format(self._A.units)) + + b = self._n.value_si + E = self._E0.value_si * 1000 # convert from J/mol to J/kmol + w = self._w0.value_si * 1000 # convert from J/mol to J/kmol + return ct.InterfaceBlowersMaselRate(A, b, E, w) + + def set_cantera_kinetics(self, ct_reaction, species_list): + """ + Takes in a cantera Reaction object and sets its + rate to a cantera InterfaceBlowersMaselRate object. + """ + import cantera as ct + if not isinstance(ct_reaction.rate, ct.InterfaceBlowersMaselRate): + raise TypeError("ct_reaction.rate must be an InterfaceBlowersMaselRate") + + # Set the rate parameter to a cantera InterfaceBlowersMaselRate object + ct_reaction.rate = self.to_cantera_kinetics() + + +################################################################################ + cdef class SurfaceChargeTransfer(KineticsModel): """ diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index 88c6c1361a6..a2b03ae67fd 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -55,7 +55,8 @@ PDepArrhenius, MultiArrhenius, MultiPDepArrhenius, get_rate_coefficient_units_from_reaction_order, \ SurfaceArrheniusBEP, StickingCoefficientBEP, ArrheniusChargeTransfer, ArrheniusChargeTransferBM, Marcus from rmgpy.kinetics.arrhenius import Arrhenius # Separate because we cimport from rmgpy.kinetics.arrhenius -from rmgpy.kinetics.surface import SurfaceArrhenius, StickingCoefficient, SurfaceChargeTransfer, SurfaceChargeTransferBEP # Separate because we cimport from rmgpy.kinetics.surface +from rmgpy.kinetics.surface import SurfaceArrhenius, StickingCoefficient, SurfaceChargeTransfer, SurfaceChargeTransferBEP, \ + SurfaceArrheniusBM # 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 @@ -399,6 +400,20 @@ def to_cantera(self, species_list=None, use_chemkin_identifier=False): rate=ct.StickingArrheniusRate() ) + elif isinstance(self.kinetics, ArrheniusBM): + if isinstance(self.kinetics, SurfaceArrheniusBM): + ct_reaction = ct.Reaction( + reactants=ct_reactants, + products=ct_products, + rate=ct.InterfaceBlowersMaselRate() + ) + else: + ct_reaction = ct.Reaction( + reactants=ct_reactants, + products=ct_products, + rate=ct.BlowersMaselRate() + ) + else: raise NotImplementedError(f"Unable to set cantera kinetics for {self.kinetics}") From ca1b288a319d8c51c90f9550c071396318da54a0 Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Wed, 11 Mar 2026 10:49:02 -0400 Subject: [PATCH 3/5] add basic tests for SurfaceArrheniusBM --- rmgpy/kinetics/surface.pyx | 17 +- test/rmgpy/kinetics/kineticsSurfaceTest.py | 246 ++++++++++++++++++++- 2 files changed, 250 insertions(+), 13 deletions(-) diff --git a/rmgpy/kinetics/surface.pyx b/rmgpy/kinetics/surface.pyx index abfad22db63..3d37f7da30e 100644 --- a/rmgpy/kinetics/surface.pyx +++ b/rmgpy/kinetics/surface.pyx @@ -814,15 +814,13 @@ cdef class SurfaceArrheniusBM(ArrheniusBM): `a`, the coefficient for exponential dependence on the coverage, `m`, the power-law exponent of coverage dependence, and `E`, the activation energy dependence on coverage. - `uncertainty` Uncertainty information `comment` Information about the model (e.g. its source) ======================= ============================================================= """ def __init__(self, A=None, n=0.0, w0=(0.0, 'J/mol'), E0=None, Tmin=None, Tmax=None, Pmin=None, Pmax=None, - coverage_dependence=None, uncertainty=None, comment=''): - KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, uncertainty=uncertainty, - comment=comment) + coverage_dependence=None, comment=''): + KineticsModel.__init__(self, Tmin=Tmin, Tmax=Tmax, Pmin=Pmin, Pmax=Pmax, comment=comment) self.A = A self.n = n self.w0 = w0 @@ -854,7 +852,7 @@ cdef class SurfaceArrheniusBM(ArrheniusBM): Return a string representation that can be used to reconstruct the SurfaceArrheniusBM object. """ - string = 'SurfaceArrheniusBM(A={0!r}, n={1!r}, w={2!r}, E0={3!r}'.format(self.A, self.n, self.w0, + string = 'SurfaceArrheniusBM(A={0!r}, n={1!r}, w0={2!r}, E0={3!r}'.format(self.A, self.n, self.w0, self.E0) if self.Tmin is not None: string += ', Tmin={0!r}'.format(self.Tmin) if self.Tmax is not None: string += ', Tmax={0!r}'.format(self.Tmax) @@ -865,7 +863,6 @@ cdef class SurfaceArrheniusBM(ArrheniusBM): string += "}" if self.Pmin is not None: string += ', Pmin={0!r}'.format(self.Pmin) if self.Pmax is not None: string += ', Pmax={0!r}'.format(self.Pmax) - if self.uncertainty is not None: string += ', uncertainty={0!r}'.format(self.uncertainty) if self.comment != '': string += ', comment="""{0}"""'.format(self.comment) string += ')' return string @@ -875,7 +872,7 @@ cdef class SurfaceArrheniusBM(ArrheniusBM): A helper function used when pickling an SurfaceArrheniusBM object. """ return (SurfaceArrheniusBM, (self.A, self.n, self.w0, self.E0, self.Tmin, self.Tmax, self.Pmin, self.Pmax, - self.uncertainty, self.coverage_dependence, self.comment)) + self.coverage_dependence, self.comment)) cpdef SurfaceArrhenius to_arrhenius(self, double dHrxn): """ @@ -893,7 +890,6 @@ cdef class SurfaceArrheniusBM(ArrheniusBM): T0=(1, "K"), Tmin=self.Tmin, Tmax=self.Tmax, - uncertainty=self.uncertainty, coverage_dependence=self.coverage_dependence, comment=self.comment, ) @@ -920,10 +916,7 @@ cdef class SurfaceArrheniusBM(ArrheniusBM): 'm^5/(mol^2*s)': 1000000, } - if self._T0.value_si != 1: - A = self._A.value_si / (self._T0.value_si) ** self._n.value_si - else: - A = self._A.value_si + A = self._A.value_si try: A *= rate_units_conversion[self._A.units] # convert from /mol to /kmol diff --git a/test/rmgpy/kinetics/kineticsSurfaceTest.py b/test/rmgpy/kinetics/kineticsSurfaceTest.py index 1f67453082e..75daf8c95f5 100644 --- a/test/rmgpy/kinetics/kineticsSurfaceTest.py +++ b/test/rmgpy/kinetics/kineticsSurfaceTest.py @@ -34,7 +34,7 @@ import numpy as np -from rmgpy.kinetics.surface import StickingCoefficient, SurfaceArrhenius, SurfaceChargeTransfer +from rmgpy.kinetics.surface import StickingCoefficient, SurfaceArrhenius, SurfaceChargeTransfer, SurfaceArrheniusBM from rmgpy.species import Species from rmgpy.molecule import Molecule import rmgpy.quantity as quantity @@ -770,3 +770,247 @@ def test_change_v0(self): self.surfchargerxn_reduction.change_v0(V) assert self.surfchargerxn_reduction.V0.value_si == V_i + delta assert round(abs(self.surfchargerxn_reduction.Ea.value_si- (Ea_i - (alpha *electrons * constants.F * delta))), 6) == 0 + + + +class TestSurfaceArrheniusBM: + """ + Contains unit tests of the :class:`SurfaceArrheniusBM` class. + """ + + def setup_class(self): + self.A = 1.44e18 + self.n = -0.087 + self.Ea = 63.4 + self.E0 = 63.4 + self.w0 = 1e3 + self.T0 = 1.0 + self.Tmin = 300.0 + self.Tmax = 3000.0 + s = Species().from_adjacency_list("1 X u0 p0 c0") + s.label = "X" + self.coverage_dependence = { + s: { + "a": quantity.Dimensionless(0.0), + "m": quantity.Dimensionless(-1.0), + "E": quantity.Energy(0.0, "J/mol"), + } + } + self.comment = "CH3x + Hx <=> CH4 + x + x" + self.surfarr = SurfaceArrhenius( + A=(self.A, "m^2/(mol*s)"), + n=self.n, + Ea=(self.Ea, "kJ/mol"), + T0=(self.T0, "K"), + Tmin=(self.Tmin, "K"), + Tmax=(self.Tmax, "K"), + comment=self.comment, + coverage_dependence=self.coverage_dependence, + ) + + self.surfarrBM = SurfaceArrheniusBM( + A=(self.A, "m^2/(mol*s)"), + n=self.n, + E0=(self.E0, "kJ/mol"), + w0=(self.w0, "kJ/mol"), + Tmin=(self.Tmin, "K"), + Tmax=(self.Tmax, "K"), + comment=self.comment, + coverage_dependence=self.coverage_dependence, + ) + + + def test_A(self): + """ + Test that the SurfaceArrheniusBM A property was properly set. + """ + assert abs(self.surfarrBM.A.value_si - self.A) < 1e0 + + def test_n(self): + """ + Test that the SurfaceArrheniusBM n property was properly set. + """ + assert round(abs(self.surfarrBM.n.value_si - self.n), 6) == 0 + + def test_E0(self): + """ + Test that the SurfaceArrheniusBM E0 property was properly set. + """ + assert round(abs(self.surfarrBM.E0.value_si * 0.001 - self.E0), 6) == 0 + + def test_w0(self): + """ + Test that the SurfaceArrheniusBM w0 property was properly set. + """ + assert round(abs(self.surfarrBM.w0.value_si * 0.001 - self.w0), 6) == 0 + + def test_Tmin(self): + """ + Test that the SurfaceArrheniusBM Tmin property was properly set. + """ + assert round(abs(self.surfarrBM.Tmin.value_si - self.Tmin), 6) == 0 + + def test_Tmax(self): + """ + Test that the SurfaceArrheniusBM Tmax property was properly set. + """ + assert round(abs(self.surfarrBM.Tmax.value_si - self.Tmax), 6) == 0 + + def test_comment(self): + """ + Test that the SurfaceArrheniusBM comment property was properly set. + """ + assert self.surfarrBM.comment == self.comment + + def test_coverage_dependence(self): + """ + Test that the coverage dependent parameters was properly set. + """ + for key in self.surfarrBM.coverage_dependence.keys(): + match = False + for key2 in self.coverage_dependence.keys(): + if key.is_identical(key2): + match = True + assert match + for species, parameters in self.surfarrBM.coverage_dependence.items(): + match = False + for species2 in self.coverage_dependence.keys(): + if species.is_identical(species2): + match = True + for key, value in parameters.items(): + assert value.value_si == self.coverage_dependence[species2][key].value_si + assert match + + def test_is_temperature_valid(self): + """ + Test the SurfaceArrheniusBM.is_temperature_valid() method. + """ + T_data = np.array([200, 400, 600, 800, 1000, 1200, 1400, 1600, 1800, 4000]) + valid_data = np.array([False, True, True, True, True, True, True, True, True, False], bool) + for T, valid in zip(T_data, valid_data): + valid0 = self.surfarr.is_temperature_valid(T) + assert valid0 == valid + + def test_pickle(self): + """ + Test that a SurfaceArrheniusBM object can be pickled and unpickled with no loss + of information. + """ + import pickle + + surfarrBM = pickle.loads(pickle.dumps(self.surfarrBM, -1)) + assert abs(self.surfarrBM.A.value - surfarrBM.A.value) < 1e0 + assert self.surfarrBM.A.units == surfarrBM.A.units + assert round(abs(self.surfarrBM.n.value - surfarrBM.n.value), 4) == 0 + assert round(abs(self.surfarrBM.E0.value - surfarrBM.E0.value), 4) == 0 + assert self.surfarrBM.E0.units == surfarrBM.E0.units + assert round(abs(self.surfarrBM.w0.value - surfarrBM.w0.value), 4) == 0 + assert self.surfarrBM.w0.units == surfarrBM.w0.units + assert round(abs(self.surfarrBM.Tmin.value - surfarrBM.Tmin.value), 4) == 0 + assert self.surfarrBM.Tmin.units == surfarrBM.Tmin.units + assert round(abs(self.surfarrBM.Tmax.value - surfarrBM.Tmax.value), 4) == 0 + assert self.surfarrBM.Tmax.units == surfarrBM.Tmax.units + assert self.surfarrBM.comment == surfarrBM.comment + for key in self.surfarrBM.coverage_dependence.keys(): + match = False + for key2 in surfarrBM.coverage_dependence.keys(): + if key.is_identical(key2): + match = True + assert match + for species, parameters in self.surfarrBM.coverage_dependence.items(): + match = False + for species2 in surfarrBM.coverage_dependence.keys(): + if species.is_identical(species2): + match = True + for key, value in parameters.items(): + assert value.value_si == surfarrBM.coverage_dependence[species2][key].value_si + assert match + assert dir(self.surfarrBM) == dir(surfarrBM) + + def test_repr(self): + """ + Test that an SurfaceArrheniusBM object can be reconstructed from its repr() + output with no loss of information. + """ + namespace = {} + exec("surfarrBM = {0!r}".format(self.surfarrBM), globals(), namespace) + assert "surfarrBM" in namespace + surfarrBM = namespace["surfarrBM"] + assert abs(self.surfarrBM.A.value - surfarrBM.A.value) < 1e0 + assert self.surfarrBM.A.units == surfarrBM.A.units + assert round(abs(self.surfarrBM.n.value - surfarrBM.n.value), 4) == 0 + assert round(abs(self.surfarrBM.E0.value - surfarrBM.E0.value), 4) == 0 + assert self.surfarrBM.E0.units == surfarrBM.E0.units + assert round(abs(self.surfarrBM.w0.value - surfarrBM.w0.value), 4) == 0 + assert self.surfarrBM.w0.units == surfarrBM.w0.units + assert round(abs(self.surfarrBM.Tmin.value - surfarrBM.Tmin.value), 4) == 0 + assert self.surfarrBM.Tmin.units == surfarrBM.Tmin.units + assert round(abs(self.surfarrBM.Tmax.value - surfarrBM.Tmax.value), 4) == 0 + assert self.surfarrBM.Tmax.units == surfarrBM.Tmax.units + assert self.surfarrBM.comment == surfarrBM.comment + assert [m.label for m in self.surfarrBM.coverage_dependence.keys()] == list(surfarrBM.coverage_dependence.keys()) + for species, parameters in self.surfarrBM.coverage_dependence.items(): + for key, value in parameters.items(): + assert value.value_si == surfarrBM.coverage_dependence[species.label][key].value_si + assert dir(self.surfarrBM) == dir(surfarrBM) + + def test_copy(self): + """ + Test that an SurfaceArrheniusBM object can be copied with deepcopy + with no loss of information. + """ + import copy + + surfarrBM = copy.deepcopy(self.surfarrBM) + assert abs(self.surfarrBM.A.value - surfarrBM.A.value) < 1e0 + assert self.surfarrBM.A.units == surfarrBM.A.units + assert round(abs(self.surfarrBM.n.value - surfarrBM.n.value), 4) == 0 + assert round(abs(self.surfarrBM.E0.value - surfarrBM.E0.value), 4) == 0 + assert self.surfarrBM.E0.units == surfarrBM.E0.units + assert round(abs(self.surfarrBM.w0.value - surfarrBM.w0.value), 4) == 0 + assert self.surfarrBM.w0.units == surfarrBM.w0.units + assert round(abs(self.surfarrBM.Tmin.value - surfarrBM.Tmin.value), 4) == 0 + assert self.surfarrBM.Tmin.units == surfarrBM.Tmin.units + assert round(abs(self.surfarrBM.Tmax.value - surfarrBM.Tmax.value), 4) == 0 + assert self.surfarrBM.Tmax.units == surfarrBM.Tmax.units + assert self.surfarrBM.comment == surfarrBM.comment + for key in self.surfarrBM.coverage_dependence.keys(): + match = False + for key2 in surfarrBM.coverage_dependence.keys(): + if key.is_identical(key2): + match = True + assert match + for species, parameters in self.surfarrBM.coverage_dependence.items(): + match = False + for species2 in surfarrBM.coverage_dependence.keys(): + if species.is_identical(species2): + match = True + for key, value in parameters.items(): + assert value.value_si == surfarrBM.coverage_dependence[species2][key].value_si + assert match + assert dir(self.surfarrBM) == dir(surfarrBM) + + def test_is_identical_to(self): + """ + Test that the SurfaceArrheniusBM.is_identical_to method works on itself + """ + assert self.surfarrBM.is_identical_to(self.surfarrBM) + + # def test_to_arrhenius(self): + # """ + # Test that the SurfaceArrheniusBM.to_arrhenius method works + # """ + + # surface_charge_transfer = self.surfarr.to_surface_charge_transfer(2,-2) + # assert isinstance(surface_charge_transfer, SurfaceChargeTransfer) + # surface_charge_transfer0 = SurfaceChargeTransfer( + # A = self.surfarr.A, + # n = self.surfarr.n, + # Ea = self.surfarr.Ea, + # T0 = self.surfarr.T0, + # Tmin = self.surfarr.Tmin, + # Tmax = self.surfarr.Tmax, + # electrons = -2, + # V0 = (2,'V') + # ) + # assert surface_charge_transfer.is_identical_to(surface_charge_transfer0) \ No newline at end of file From 6f81673d40dc777ffdf21d02645d94d4a53d517f Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Wed, 11 Mar 2026 14:31:39 -0400 Subject: [PATCH 4/5] update ArrheniusBM to handle surface units --- rmgpy/kinetics/arrhenius.pyx | 12 +++++++----- rmgpy/quantity.py | 21 +++++++++++++++++++++ 2 files changed, 28 insertions(+), 5 deletions(-) diff --git a/rmgpy/kinetics/arrhenius.pyx b/rmgpy/kinetics/arrhenius.pyx index ea9e462d4f8..146bc26fe4a 100644 --- a/rmgpy/kinetics/arrhenius.pyx +++ b/rmgpy/kinetics/arrhenius.pyx @@ -688,12 +688,14 @@ cdef class ArrheniusBM(KineticsModel): self.comment = 'Fitted to {0} reactions at temperatures: {1}'.format(len(rxns), Ts) # fill in parameters - A_units = ['', 's^-1', 'm^3/(mol*s)', 'm^6/(mol^2*s)'] - order = len(rxns[0].reactants) - if order != 1 and rxn.is_surface_reaction(): - raise NotImplementedError("Units not implemented for surface reactions.") - self.A = (A, A_units[order]) + if rxn.is_surface_reaction(): + A_units = quantity.SURFACERATECOEFFICIENT_SI_UNITS[rxn.kinetics.A.units] + else: + A_units = ['', 's^-1', 'm^3/(mol*s)', 'm^6/(mol^2*s)'] + order = len(rxns[0].reactants) + A_units = A_units[order] + self.A = (A, A_units) self.n = n self.w0 = (w0, 'J/mol') self.E0 = (E0, 'J/mol') diff --git a/rmgpy/quantity.py b/rmgpy/quantity.py index 5d73546d455..208147ad688 100644 --- a/rmgpy/quantity.py +++ b/rmgpy/quantity.py @@ -729,6 +729,7 @@ def __call__(self, *args, **kwargs): return quantity + Acceleration = UnitType('m/s^2') Area = UnitType('m^2') @@ -885,3 +886,23 @@ def SurfaceRateCoefficient(*args, **kwargs): # Return the Quantity or ArrayQuantity object object return quantity + +SURFACERATECOEFFICIENT_SI_UNITS = { + 's^-1': 's^-1', + 'm^3/(mol*s)': 'm^3/(mol*s)', + 'cm^3/(mol*s)': 'm^3/(mol*s)', + 'm^3/(molecule*s)': 'm^3/(mol*s)', + 'cm^3/(molecule*s)': 'm^3/(mol*s)', + 'm^2/(mol*s)': 'm^2/(mol*s)', + 'cm^2/(mol*s)': 'm^2/(mol*s)', + 'm^2/(molecule*s)': 'm^2/(mol*s)', + 'cm^2/(molecule*s)': 'm^2/(mol*s)', + 'm^5/(mol^2*s)': 'm^5/(mol^2*s)', + 'cm^5/(mol^2*s)': 'm^5/(mol^2*s)', + 'm^5/(molecule^2*s)': 'm^5/(mol^2*s)', + 'cm^5/(molecule^2*s)': 'm^5/(mol^2*s)', + 'm^4/(mol^2*s)': 'm^4/(mol^2*s)', + 'cm^4/(mol^2*s)': 'm^4/(mol^2*s)', + 'm^4/(molecule^2*s)': 'm^4/(mol^2*s)', + 'cm^4/(molecule^2*s)': 'm^4/(mol^2*s)', +} From 6dd861fb48798f0be2d05f6fbd37c5a7a5afabbd Mon Sep 17 00:00:00 2001 From: Sevy Harris Date: Wed, 11 Mar 2026 14:49:41 -0400 Subject: [PATCH 5/5] add fitting test for SurfaceArrheniusBM --- test/rmgpy/kinetics/kineticsSurfaceTest.py | 92 +++++++++++++++------- 1 file changed, 65 insertions(+), 27 deletions(-) diff --git a/test/rmgpy/kinetics/kineticsSurfaceTest.py b/test/rmgpy/kinetics/kineticsSurfaceTest.py index 75daf8c95f5..66ce14e35ec 100644 --- a/test/rmgpy/kinetics/kineticsSurfaceTest.py +++ b/test/rmgpy/kinetics/kineticsSurfaceTest.py @@ -31,14 +31,17 @@ This script contains unit tests of the :mod:`rmgpy.kinetics.surface` module. """ - +import copy +import pickle import numpy as np from rmgpy.kinetics.surface import StickingCoefficient, SurfaceArrhenius, SurfaceChargeTransfer, SurfaceArrheniusBM from rmgpy.species import Species +from rmgpy.reaction import Reaction from rmgpy.molecule import Molecule import rmgpy.quantity as quantity import rmgpy.constants as constants +import rmgpy.thermo class TestStickingCoefficient: @@ -150,7 +153,6 @@ def test_pickle(self): Test that an StickingCoefficient object can be pickled and unpickled with no loss of information. """ - import pickle stick = pickle.loads(pickle.dumps(self.stick, -1)) assert abs(self.stick.A.value - stick.A.value) < 1e0 @@ -213,7 +215,6 @@ def test_copy(self): Test that an StickingCoefficient object can be copied with deepcopy with no loss of information. """ - import copy stick = copy.deepcopy(self.stick) assert abs(self.stick.A.value - stick.A.value) < 1e0 @@ -360,7 +361,6 @@ def test_pickle(self): Test that an SurfaceArrhenius object can be pickled and unpickled with no loss of information. """ - import pickle surfarr = pickle.loads(pickle.dumps(self.surfarr, -1)) assert abs(self.surfarr.A.value - surfarr.A.value) < 1e0 @@ -423,7 +423,6 @@ def test_copy(self): Test that an SurfaceArrhenius object can be copied with deepcopy with no loss of information. """ - import copy surfarr = copy.deepcopy(self.surfarr) assert abs(self.surfarr.A.value - surfarr.A.value) < 1e0 @@ -604,7 +603,6 @@ def test_pickle(self): Test that an SurfaceChargeTransfer object can be pickled and unpickled with no loss of information. """ - import pickle surfchargerxn_reduction = pickle.loads(pickle.dumps(self.surfchargerxn_reduction, -1)) assert abs(self.surfchargerxn_reduction.A.value-surfchargerxn_reduction.A.value) < 1e0 assert self.surfchargerxn_reduction.A.units == surfchargerxn_reduction.A.units @@ -651,7 +649,6 @@ def test_copy(self): Test that an SurfaceChargeTransfer object can be copied with deepcopy with no loss of information. """ - import copy surfchargerxn_reduction = copy.deepcopy(self.surfchargerxn_reduction) assert abs(self.surfchargerxn_reduction.A.value-surfchargerxn_reduction.A.value) < 1e0 assert self.surfchargerxn_reduction.A.units == surfchargerxn_reduction.A.units @@ -819,6 +816,48 @@ def setup_class(self): coverage_dependence=self.coverage_dependence, ) + self.reaction1 = Reaction() + X = Species().from_adjacency_list('1 X u0 p0 c0\n') + OX = Species(smiles='O=*') + O2 = Species(smiles='[O][O]') + X.thermo = rmgpy.thermo.NASA( + polynomials = [ + rmgpy.thermo.NASAPolynomial(coeffs=[ + 0.000000000E+00, 0.000000000E+00, 0.000000000E+00, 0.000000000E+00, + 0.000000000E+00, 0.000000000E+00, 0.000000000E+00], Tmin=(298.0,'K'), Tmax=(1000.0, 'K')), + rmgpy.thermo.NASAPolynomial(coeffs=[ + 0.000000000E+00, 0.000000000E+00, 0.000000000E+00, 0.000000000E+00, + 0.000000000E+00, 0.000000000E+00, 0.000000000E+00], Tmin=(1000.0,'K'), Tmax=(3000.0, 'K')), + ], + Tmin = (298.0, 'K'), + Tmax = (3000.0, 'K'), + ) + OX.thermo = rmgpy.thermo.NASA( + polynomials=[ + rmgpy.thermo.NASAPolynomial(coeffs=[-2.94475701E-01, 1.44162624E-02, -2.61322704E-05, 2.19005957E-08, -6.98019420E-12, + -1.64619234E+04, -1.99445623E-01], Tmin=(298.0, 'K'), Tmax=(1000.0, 'K')), + rmgpy.thermo.NASAPolynomial(coeffs=[2.90244691E+00, -3.38584457E-04, 6.43372619E-07, -3.66326660E-10, 6.90093884E-14, + -1.70497471E+04, -1.52559728E+01], Tmin=(1000.0, 'K'), Tmax=(2000.0, 'K')), + ], + Tmin=(298.0, 'K'), + Tmax=(2000.0, 'K'), + ) + O2.thermo = rmgpy.thermo.NASA( + polynomials=[rmgpy.thermo.NASAPolynomial(coeffs=[3.53732,-0.00121571,5.31618e-06,-4.89443e-09,1.45845e-12,-1038.59,4.68368], Tmin=(100,'K'), Tmax=(1074.56,'K')), + rmgpy.thermo.NASAPolynomial(coeffs=[3.15382,0.00167804,-7.69971e-07,1.51275e-10,-1.08782e-14,-1040.82,6.16754], Tmin=(1074.56,'K'), Tmax=(5000,'K'))], + Tmin=(100,'K'), + Tmax=(5000,'K') + ) + self.reaction1.reactants = [O2, X, X] + self.reaction1.products = [OX, OX] + self.reaction1.kinetics = rmgpy.kinetics.surface.SurfaceArrhenius( + A=(1.89E21, 'cm^5/(mol^2*s)'), + n = -0.5, + Ea=(0.0, 'J/mol'), + Tmin = (200, 'K'), + Tmax = (3000, 'K'), + ) + def test_A(self): """ @@ -896,7 +935,6 @@ def test_pickle(self): Test that a SurfaceArrheniusBM object can be pickled and unpickled with no loss of information. """ - import pickle surfarrBM = pickle.loads(pickle.dumps(self.surfarrBM, -1)) assert abs(self.surfarrBM.A.value - surfarrBM.A.value) < 1e0 @@ -959,7 +997,6 @@ def test_copy(self): Test that an SurfaceArrheniusBM object can be copied with deepcopy with no loss of information. """ - import copy surfarrBM = copy.deepcopy(self.surfarrBM) assert abs(self.surfarrBM.A.value - surfarrBM.A.value) < 1e0 @@ -996,21 +1033,22 @@ def test_is_identical_to(self): """ assert self.surfarrBM.is_identical_to(self.surfarrBM) - # def test_to_arrhenius(self): - # """ - # Test that the SurfaceArrheniusBM.to_arrhenius method works - # """ - - # surface_charge_transfer = self.surfarr.to_surface_charge_transfer(2,-2) - # assert isinstance(surface_charge_transfer, SurfaceChargeTransfer) - # surface_charge_transfer0 = SurfaceChargeTransfer( - # A = self.surfarr.A, - # n = self.surfarr.n, - # Ea = self.surfarr.Ea, - # T0 = self.surfarr.T0, - # Tmin = self.surfarr.Tmin, - # Tmax = self.surfarr.Tmax, - # electrons = -2, - # V0 = (2,'V') - # ) - # assert surface_charge_transfer.is_identical_to(surface_charge_transfer0) \ No newline at end of file + + def test_to_and_from_arrhenius(self): + """Test going from SurfaceArrhenius to SurfaceArrheniusBM and back again""" + + reaction1_BM = copy.deepcopy(self.reaction1) + reaction1_BM.kinetics = rmgpy.kinetics.surface.SurfaceArrheniusBM().fit_to_reactions([self.reaction1], w0=1e6) + return_arrhenius = reaction1_BM.kinetics.to_arrhenius(dHrxn=reaction1_BM.get_enthalpy_of_reaction(298)) + assert return_arrhenius.is_similar_to(self.reaction1.kinetics) + + # not sure if .is_similar_to is wrong... because as the enthalpy of reaction changes + # assert self.reaction1.kinetics.is_similar_to(reaction1_BM.kinetics) + Ts = [500, 1000, 1500, 2000] + for T in Ts: + k1 = self.reaction1.kinetics.get_rate_coefficient(T) + dHrxn = reaction1_BM.get_enthalpy_of_reaction(T) + k2 = reaction1_BM.kinetics.get_rate_coefficient(T, dHrxn=dHrxn) + assert np.abs(np.log10(k1) - np.log10(k2)) < 0.1 + +