diff --git a/arkane/pdep.py b/arkane/pdep.py index 7f0e6b22d3..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 @@ -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: diff --git a/rmgpy/pdep/configuration.pyx b/rmgpy/pdep/configuration.pyx index e44626a923..e09b6d598b 100644 --- a/rmgpy/pdep/configuration.pyx +++ b/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 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: 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 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) diff --git a/test/arkane/commonTest.py b/test/arkane/commonTest.py index 9afec8f591..336a8c3343 100644 --- a/test/arkane/commonTest.py +++ b/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): """ 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): """