From 4782a7c59399c5db81266a05b9d334edab32f1ee Mon Sep 17 00:00:00 2001 From: Richard West Date: Tue, 15 Jul 2025 14:57:11 -0400 Subject: [PATCH 01/10] Add note to docstring about adding energy_correction The property E0 when accessed from a Configuration, includes the energy_correction which was added in 169a34ec00c1e4275e38a690a49fc1d86bdd8b4c --- rmgpy/pdep/configuration.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/pdep/configuration.pyx b/rmgpy/pdep/configuration.pyx index e44626a923..f6854ca63a 100644 --- a/rmgpy/pdep/configuration.pyx +++ b/rmgpy/pdep/configuration.pyx @@ -77,7 +77,7 @@ cdef class Configuration(object): return string property E0: - """The ground-state energy of the configuration in J/mol.""" + """The ground-state energy of the configuration in J/mol. Applies the energy_correction.""" def __get__(self): return sum([float(spec.conformer.E0.value_si) for spec in self.species]) + self.energy_correction From f60ac9c5e7916fe7fe2b02b8dcf1416d5a210f60 Mon Sep 17 00:00:00 2001 From: Richard West Date: Tue, 15 Jul 2025 15:01:51 -0400 Subject: [PATCH 02/10] Report energy_correction in __repr__ of a Configuration. This could help when reporting Configuration objects, and at least makes debugging less misleading. Without this offset, all your energies are wrong. See https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2791. But this alone doesn't fix the issue since Configuration objects are not printed inte Arkane input files using __repr__. We should generally have better unit tests that __repr__ reproduces things. --- rmgpy/pdep/configuration.pyx | 1 + 1 file changed, 1 insertion(+) diff --git a/rmgpy/pdep/configuration.pyx b/rmgpy/pdep/configuration.pyx index f6854ca63a..e09b6d598b 100644 --- a/rmgpy/pdep/configuration.pyx +++ b/rmgpy/pdep/configuration.pyx @@ -73,6 +73,7 @@ cdef class Configuration(object): if self.sum_states is not None: string += 'sum_states={0}, '.format(self.sum_states) string += 'active_k_rotor={0}, '.format(self.active_k_rotor) string += 'active_j_rotor={0}, '.format(self.active_j_rotor) + if self.energy_correction != 0.0: string += 'energy_correction={0}, '.format(self.energy_correction) string += ')' return string From e089a374e2a3786298ed1630a5476387693bbfa6 Mon Sep 17 00:00:00 2001 From: Richard West Date: Wed, 16 Jul 2025 14:21:21 -0400 Subject: [PATCH 03/10] Use lowest energy Configuration for energy_correction. As noted in https://github.com/ReactionMechanismGenerator/RMG-Py/pull/2791 the energy_correction was not actually an average of the Configuration energies because some Configurations (reactants and products) are bimolecular. This corrects the calculation, but also switches to using the MINIMUM rather than the average, so that the lowest should end up at zero. --- rmgpy/rmg/pdep.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/rmgpy/rmg/pdep.py b/rmgpy/rmg/pdep.py index 54d62df106..7e5050daf0 100644 --- a/rmgpy/rmg/pdep.py +++ b/rmgpy/rmg/pdep.py @@ -788,25 +788,21 @@ def update(self, reaction_model, pdep_settings, requires_rms=False): # Log the network being updated logging.info("Updating {0!s}".format(self)) - E0 = [] # Generate states data for unimolecular isomers and reactants if necessary for isomer in self.isomers: spec = isomer.species[0] if not spec.has_statmech(): spec.generate_statmech() - E0.append(spec.conformer.E0.value_si) for reactants in self.reactants: for spec in reactants.species: if not spec.has_statmech(): spec.generate_statmech() - E0.append(spec.conformer.E0.value_si) # Also generate states data for any path reaction reactants, so we can # always apply the ILT method in the direction the kinetics are known for reaction in self.path_reactions: for spec in reaction.reactants: if not spec.has_statmech(): spec.generate_statmech() - E0.append(spec.conformer.E0.value_si) # While we don't need the frequencies for product channels, we do need # the E0, so create a conformer object with the E0 for the product # channel species if necessary @@ -814,12 +810,12 @@ def update(self, reaction_model, pdep_settings, requires_rms=False): for spec in products.species: if spec.conformer is None: spec.conformer = Conformer(E0=spec.get_thermo_data().E0) - E0.append(spec.conformer.E0.value_si) - # Use the average E0 as the reference energy (`energy_correction`) for the network + # Use the lowest E0 as the reference energy (`energy_correction`) for the network # The `energy_correction` will be added to the free energies and enthalpies for each # configuration in the network. - energy_correction = -np.array(E0).mean() + energy_correction = -min(sum(spec.conformer.E0.value_si for spec in stationary_point.species) + for stationary_point in self.reactants + self.isomers + self.products) for spec in self.reactants + self.products + self.isomers: spec.energy_correction = energy_correction self.energy_correction = energy_correction From 104442ddfd66230270e78abe3efdc6bbf0d3c024 Mon Sep 17 00:00:00 2001 From: Richard West Date: Wed, 16 Jul 2025 16:20:18 -0400 Subject: [PATCH 04/10] Remove energy_correction from TS E0 values when writing Arkane input file The E0 values for TS have been modified using the energy_correction, but the E0 values for species have not been. We can't change the species, because it needs to be applied to the Configuration (some are bimolecular, some are unimolecular), so to be consistent we must remove it from the TS. This leaves both values, and a comment. It's a bit ugly, and I think there should be a better fix later, but this I think solves an immediate bug, and now the network_X_Y.py files written during an RMG job are at least internally consistent. (They have a different definition of 0 from inside the RMG job, but they're self-consistent.) --- arkane/pdep.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/arkane/pdep.py b/arkane/pdep.py index 7f0e6b22d3..db18a1d706 100644 --- a/arkane/pdep.py +++ b/arkane/pdep.py @@ -699,7 +699,10 @@ def save_input_file(self, path): f.write(' label = {0!r},\n'.format(ts.label)) if ts.conformer is not None: if ts.conformer.E0 is not None: - f.write(' E0 = {0!r},\n'.format(ts.conformer.E0)) + if self.network.energy_correction: + f.write(f' E0 = ({ts.conformer.E0.value_si * 0.001:.3f} - {self.network.energy_correction * 0.001:.3f}, "kJ/mol"), # removing the applied energy_correction\n') + else: + f.write(' E0 = {0!r},\n'.format(ts.conformer.E0)) if len(ts.conformer.modes) > 0: f.write(' modes = [\n') for mode in ts.conformer.modes: From a25dda7e069b559ad8697e8a010938fb0a6ab268 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 17 Jul 2025 14:09:23 -0400 Subject: [PATCH 05/10] Fix diffmodels.py exception handling for surface mechanisms. When comparing mechanisms that are NOT surface models, it would do the load_chemkin_file as part of an `except` block, which then means if you get an error during `load_chemkin_file` the stack trace and error message are all confused and it looks like a KeyError: 'surface_path1' which is misleading. --- rmgpy/tools/diffmodels.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/rmgpy/tools/diffmodels.py b/rmgpy/tools/diffmodels.py index e93609a77e..6d24a2d686 100644 --- a/rmgpy/tools/diffmodels.py +++ b/rmgpy/tools/diffmodels.py @@ -361,17 +361,16 @@ def execute(chemkin1, species_dict1, thermo1, chemkin2, species_dict2, thermo2, chemkin2 = chemkin_gas2 kwargs['surface_path2'] = chemkin_surface2 - try: + if 'surface_path1' in kwargs and 'surface_path2' in kwargs: surface_path1 = kwargs['surface_path1'] surface_path2 = kwargs['surface_path2'] model1.species, model1.reactions = load_chemkin_file( chemkin1, species_dict1, thermo_path=thermo1, surface_path=surface_path1) model2.species, model2.reactions = load_chemkin_file( chemkin2, species_dict2, thermo_path=thermo2, surface_path=surface_path2) - except KeyError: - if 'surface_path1' in kwargs or 'surface_path2' in kwargs: - logging.warning('Please specify 2 surface input files if you are comparing a surface mechanism') - + elif 'surface_path1' in kwargs or 'surface_path2' in kwargs: + raise ValueError('Please specify 2 surface input files if you are comparing a surface mechanism') + else: model1.species, model1.reactions = load_chemkin_file(chemkin1, species_dict1, thermo_path=thermo1) model2.species, model2.reactions = load_chemkin_file(chemkin2, species_dict2, thermo_path=thermo2) From 717010c57914670469af25e1f4be1cc42ebf4802 Mon Sep 17 00:00:00 2001 From: Richard West Date: Thu, 17 Jul 2025 14:11:16 -0400 Subject: [PATCH 06/10] Use the Tmin,Tmax,Pmin,Pmax if specified in an Arkane job. For the Chebyshev fitting it used to ignore the specified limits and just use the lowest and highest points from the grid. Chebyshev grid points aren't on the boundaries of the valid range, so you'd end up with not the range you specified. One consequence is network.py files created during an RMG job would not reproduce the fits if you later ran them in Arkane. Without this change, an RMG job would create TCHEB/ 300.000 2500.000 / PCHEB/ 0.010 98.692 / but then when you run the network1_1.py file through Arkane you would get TCHEB/ 302.558 2335.460 / PCHEB/ 0.012 78.776 / --- arkane/pdep.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/arkane/pdep.py b/arkane/pdep.py index db18a1d706..8e61304b59 100644 --- a/arkane/pdep.py +++ b/arkane/pdep.py @@ -132,9 +132,11 @@ def __init__(self, network, if Tlist is not None: self.Tlist = Tlist - self.Tmin = (np.min(self.Tlist.value_si), "K") - self.Tmax = (np.max(self.Tlist.value_si), "K") self.Tcount = len(self.Tlist.value_si) + if self.Tmin is None: + self.Tmin = (np.min(self.Tlist.value_si), "K") + if self.Tmax is None: + self.Tmax = (np.max(self.Tlist.value_si), "K") else: self.Tlist = None @@ -143,9 +145,11 @@ def __init__(self, network, self.Pcount = Pcount if Plist is not None: self.Plist = Plist - self.Pmin = (np.min(self.Plist.value_si) * 1e-5, "bar") - self.Pmax = (np.max(self.Plist.value_si) * 1e-5, "bar") self.Pcount = len(self.Plist.value_si) + if self.Pmin is None: + self.Pmin = (np.min(self.Plist.value_si) * 1e-5, "bar") + if self.Pmax is None: + self.Pmax = (np.max(self.Plist.value_si) * 1e-5, "bar") else: self.Plist = None From e279681b925820c5100e839edfbf2841a538c230 Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 18 Jul 2025 10:18:34 -0400 Subject: [PATCH 07/10] Fix string formatting error. the {reaction:s} was causing a TypeError: unsupported format string passed to TemplateReaction.__format__ but {reaction} should do the intended conversion. This line is seldom hit. --- rmgpy/rmg/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index d677755b07..0eb946dee3 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -501,7 +501,7 @@ def make_new_reaction(self, forward, check_existing=True, generate_thermo=True, reactants = [self.make_new_species(reactant, generate_thermo=generate_thermo)[0] for reactant in forward.reactants] products = [self.make_new_species(product, generate_thermo=generate_thermo)[0] for product in forward.products] except: - logging.error(f"Error when making species in reaction {forward:s} from {forward.family:s}") + logging.error(f"Error when making species in reaction {forward} from {forward.family}") raise if forward.specific_collider is not None: From d6767248c2728c44179a5305a36407df3b978cf0 Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 18 Jul 2025 11:29:28 -0400 Subject: [PATCH 08/10] Fix unit tests for Arkane reading of P and T limits. Because the T limits were being read incorrectly, a test of the rate coefficient that was supposedly at Tmax (but was in fact at 700K) was failing once the T limit reading was fixed. Also the pressure limit reading test was using .value instead of .value_si so that broke also. --- test/arkane/commonTest.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/test/arkane/commonTest.py b/test/arkane/commonTest.py index 9afec8f591..e415ae4d9d 100644 --- a/test/arkane/commonTest.py +++ b/test/arkane/commonTest.py @@ -110,9 +110,11 @@ def setup_class(cls): arkane = Arkane() job_list = arkane.load_input_file(os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "..", "arkane", "data", "methoxy.py")) pdepjob = job_list[-1] + cls.pdepjob = pdepjob cls.kineticsjob = job_list[0] pdepjob.active_j_rotor = True network = pdepjob.network + cls.network = network cls.Nisom = len(network.isomers) cls.Nreac = len(network.reactants) cls.Nprod = len(network.products) @@ -169,11 +171,12 @@ def test_temperatures_units(self): """ assert str(self.TmaxUnits) == "K" - def test_temperatures_value(self): + def test_temperatures_limits(self): """ - Test the temperature value. + Test the temperature limits. """ - assert self.TminValue == 450.0 + assert self.pdepjob.Tmin.value_si == 450.0 + assert self.pdepjob.Tmax.value_si == 1200.0 def test_temperatures_list(self): """ @@ -181,11 +184,12 @@ def test_temperatures_list(self): """ assert np.array_equal(self.TlistValue, np.array([450, 500, 678, 700])) - def test_min_pressure_value(self): + def test_pressure_limits(self): """ - Test the minimum pressure value. + Test the pressure limits. """ - assert "%0.7f" % self.PminValue == str(0.0101325) + assert self.pdepjob.Pmin.value_si == 1013.25 # Pa + assert self.pdepjob.Pmax.value_si == 101325000.0 # Pa def test_pressure_count(self): """ @@ -234,8 +238,8 @@ def test_calculate_tst_rate_coefficient(self): """ Test the calculation of the high-pressure limit rate coef for one of the kinetics jobs at Tmin and Tmax. """ - assert "%0.7f" % self.kineticsjob.reaction.calculate_tst_rate_coefficient(self.TminValue) == str(46608.5904933) - assert "%0.5f" % self.kineticsjob.reaction.calculate_tst_rate_coefficient(self.Tmaxvalue) == str(498796.64535) + assert "%0.7f" % self.kineticsjob.reaction.calculate_tst_rate_coefficient(450.0) == str(46608.5904933) + assert "%0.5f" % self.kineticsjob.reaction.calculate_tst_rate_coefficient(700.0) == str(498796.64535) def test_tunneling(self): """ From d5c011368852439b67ffc26a67eca13f912e0bbf Mon Sep 17 00:00:00 2001 From: Richard West Date: Fri, 18 Jul 2025 11:37:21 -0400 Subject: [PATCH 09/10] Clean up remaining tests in arkane/commonTest.py This was harder to read, parse, understand, or maintain. Also added a few extra asserts (units etc.) --- test/arkane/commonTest.py | 47 +++++++++++++-------------------------- 1 file changed, 16 insertions(+), 31 deletions(-) diff --git a/test/arkane/commonTest.py b/test/arkane/commonTest.py index e415ae4d9d..336a8c3343 100644 --- a/test/arkane/commonTest.py +++ b/test/arkane/commonTest.py @@ -115,61 +115,44 @@ def setup_class(cls): pdepjob.active_j_rotor = True network = pdepjob.network cls.network = network - cls.Nisom = len(network.isomers) - cls.Nreac = len(network.reactants) - cls.Nprod = len(network.products) - cls.Npath = len(network.path_reactions) - cls.PathReaction2 = network.path_reactions[2] - cls.TminValue = pdepjob.Tmin.value - cls.Tmaxvalue = pdepjob.Tmax.value - cls.TmaxUnits = pdepjob.Tmax.units - cls.TlistValue = pdepjob.Tlist.value - cls.PminValue = pdepjob.Pmin.value - cls.Pcount = pdepjob.Pcount - cls.Tcount = pdepjob.Tcount - cls.GenTlist = pdepjob.generate_T_list() - cls.PlistValue = pdepjob.Plist.value - cls.maximum_grain_size_value = pdepjob.maximum_grain_size.value - cls.method = pdepjob.method - cls.rmgmode = pdepjob.rmgmode # test Arkane's interactions with the network module def test_num_isom(self): """ Test the number of isomers identified. """ - assert self.Nisom == 2 + assert len(self.network.isomers) == 2 def test_num_reac(self): """ Test the number of reactants identified. """ - assert self.Nreac == 1 + assert len(self.network.reactants) == 1 def test_num_prod(self): """ Test the number of products identified. """ - assert self.Nprod == 1 + assert len(self.network.products) == 1 def test_n_path_reactions(self): """ Test the whether or not RMG mode is turned on. """ - assert self.Npath == 3 + assert len(self.network.path_reactions) == 3 def test_path_reactions(self): """ Test a path reaction label """ - assert str(self.PathReaction2) == "CH2OH <=> methoxy" + assert str(self.network.path_reactions[2]) == "CH2OH <=> methoxy" # test Arkane's interactions with the pdep module def test_temperatures_units(self): """ Test the Temperature Units. """ - assert str(self.TmaxUnits) == "K" + assert str(self.pdepjob.Tmax.units) == "K" def test_temperatures_limits(self): """ @@ -182,7 +165,7 @@ def test_temperatures_list(self): """ Test the temperature list. """ - assert np.array_equal(self.TlistValue, np.array([450, 500, 678, 700])) + assert np.array_equal(self.pdepjob.Tlist.value_si, np.array([450, 500, 678, 700])) def test_pressure_limits(self): """ @@ -195,43 +178,45 @@ def test_pressure_count(self): """ Test the number pressures specified. """ - assert self.Pcount == 7 + assert self.pdepjob.Pcount == 7 def test_temperature_count(self): """ Test the number temperatures specified. """ - assert self.Tcount == 4 + assert self.pdepjob.Tcount == 4 def test_pressure_list(self): """ Test the pressure list. """ - assert np.array_equal(self.PlistValue, np.array([0.01, 0.1, 1, 3, 10, 100, 1000])) + assert np.array_equal(self.pdepjob.Plist.value, np.array([0.01, 0.1, 1, 3, 10, 100, 1000])) + assert self.pdepjob.Plist.units == 'atm' def test_generate_temperature_list(self): """ Test the generated temperature list. """ - assert list(self.GenTlist) == [450.0, 500.0, 678.0, 700.0] + assert list(self.pdepjob.generate_T_list()) == [450.0, 500.0, 678.0, 700.0] def test_maximum_grain_size_value(self): """ Test the max grain size value. """ - assert self.maximum_grain_size_value == 0.5 + assert self.pdepjob.maximum_grain_size.value == 0.5 + assert self.pdepjob.maximum_grain_size.units == 'kcal/mol' def test_method(self): """ Test the master equation solution method chosen. """ - assert self.method == "modified strong collision" + assert self.pdepjob.method == "modified strong collision" def test_rmg_mode(self): """ Test the whether or not RMG mode is turned on. """ - assert self.rmgmode == False + assert self.pdepjob.rmgmode == False # Test Arkane's interactions with the kinetics module def test_calculate_tst_rate_coefficient(self): From 116301dc3b9989e5b348cc189199a9fbd6c1c3c7 Mon Sep 17 00:00:00 2001 From: Richard West Date: Sat, 11 Oct 2025 22:59:04 -0400 Subject: [PATCH 10/10] Change rounding precision in a unit test. It was spuriously failing on the github runner for continuous integration. Unrelated to changes on this pull request. --- test/rmgpy/reactionTest.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/rmgpy/reactionTest.py b/test/rmgpy/reactionTest.py index 4e1cc0cc50..887b834a5e 100644 --- a/test/rmgpy/reactionTest.py +++ b/test/rmgpy/reactionTest.py @@ -2936,9 +2936,9 @@ def test_multi_arrhenius(self): # Check that the reaction string is the same assert repr(converted_rxn) == repr(ct_rxn) # Check that the Arrhenius rates are identical - assert round(abs(converted_rxn.rate.pre_exponential_factor - ct_rxn.rate.pre_exponential_factor), 3) == 0 - assert round(abs(converted_rxn.rate.temperature_exponent - ct_rxn.rate.temperature_exponent), 7) == 0 - assert round(abs(converted_rxn.rate.activation_energy - ct_rxn.rate.activation_energy), 7) == 0 + assert round((converted_rxn.rate.pre_exponential_factor / ct_rxn.rate.pre_exponential_factor), 7) == 1 + assert round((converted_rxn.rate.temperature_exponent / ct_rxn.rate.temperature_exponent), 7) == 1 + assert round((converted_rxn.rate.activation_energy / ct_rxn.rate.activation_energy), 7) == 1 def test_pdep_arrhenius(self): """