From a3baf03459118c00c22da8cb64b6eec02f9d15a3 Mon Sep 17 00:00:00 2001 From: Anna Doner <47831564+donerancl@users.noreply.github.com> Date: Tue, 13 May 2025 13:52:48 -0400 Subject: [PATCH] Remove incorrect energy_correction from pdep.py After finding many submerged barriers when running RMG in pdep mode, I discovered that the energy_correction is calculated incorrectly *and* only applied to the transition state. The code calculates the energy_correction by averaging all species E0s. However, it should average the total energy for each stationary point (some of which are bimolecular), and I believe it would be more intuitive to use the minimum energy on the surface. A simple solution would be to remove the energy_correction altogether (what I've done here for now). If it is desired, we need to apply energy_correction to all stationary points, not just transition states, so that barriers are not submerged. The energy_correction could be calculated as follows... ``` stationary_point_E0s = [sum([spec.conformer.E0.value_si for spec in stationary_point]) for stationary_point in self.reactants + self.isomers + self.products])] energy_correction = -np.array(stationary_point_E0s).min() ``` I would need to track down where to apply it to the wells though. --- rmgpy/rmg/pdep.py | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/rmgpy/rmg/pdep.py b/rmgpy/rmg/pdep.py index 0e1398e51a7..bd0a2ab09d6 100644 --- a/rmgpy/rmg/pdep.py +++ b/rmgpy/rmg/pdep.py @@ -816,14 +816,6 @@ def update(self, reaction_model, pdep_settings): 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 - # The `energy_correction` will be added to the free energies and enthalpies for each - # configuration in the network. - energy_correction = -np.array(E0).mean() - for spec in self.reactants + self.products + self.isomers: - spec.energy_correction = energy_correction - self.energy_correction = energy_correction - # Determine transition state energies on potential energy surface # In the absence of any better information, we simply set it to # be the reactant ground-state energy + the activation energy @@ -851,9 +843,9 @@ def update(self, reaction_model, pdep_settings): 'type "{2!s}".'.format(rxn, self.index, rxn.kinetics.__class__)) rxn.fix_barrier_height(force_positive=True) if rxn.network_kinetics is None: - E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.kinetics.Ea.value_si + energy_correction + E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.kinetics.Ea.value_si else: - E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.network_kinetics.Ea.value_si + energy_correction + E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants]) + rxn.network_kinetics.Ea.value_si rxn.transition_state = rmgpy.species.TransitionState(conformer=Conformer(E0=(E0 * 0.001, "kJ/mol"))) # Set collision model