Fix E0 energy_correction in pressure dependence calculations.

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
 
@@ -699,7 +703,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:

rmgpy/pdep/configuration.pyx

@@ -73,11 +73,12 @@ 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
 
     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 
 

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:

rmgpy/rmg/pdep.py

@@ -788,38 +788,34 @@ 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
         for products in self.products:
             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

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)
 

test/arkane/commonTest.py

@@ -110,132 +110,121 @@ 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.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
+        cls.network = network
 
     # 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_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):
         """
         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_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):
         """
         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):
         """
         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):
         """