Skip to content

Remove incorrectly calculated and applied energy_correction from pdep.py #2791

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from

Conversation

donerancl
Copy link
Contributor

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.

Motivation or Problem

A clear and concise description of what what you're trying to fix or improve. Please reference any issues that this addresses.

Description of Changes

A clear and concise description of what what you've changed or added.

Testing

A clear and concise description of testing that you've done or plan to do.

Reviewer Tips

Suggestions for verifying that this PR works or other notes for the reviewer.

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.
Copy link

Regression Testing Results

WARNING:root:Initial mole fractions do not sum to one; normalizing.
WARNING:root:Initial mole fractions do not sum to one; normalizing.
WARNING:root:Initial mole fractions do not sum to one; normalizing.
⚠️ One or more regression tests failed.
Please download the failed results and run the tests locally or check the log to see why.

Detailed regression test results.

Regression test aromatics:

Reference: Execution time (DD:HH:MM:SS): 00:00:01:09
Current: Execution time (DD:HH:MM:SS): 00:00:01:09
Reference: Memory used: 2748.31 MB
Current: Memory used: 2731.22 MB

aromatics Passed Core Comparison ✅

Original model has 15 species.
Test model has 15 species. ✅
Original model has 11 reactions.
Test model has 11 reactions. ✅

aromatics Failed Edge Comparison ❌

Original model has 106 species.
Test model has 106 species. ✅
Original model has 358 reactions.
Test model has 358 reactions. ✅

Non-identical thermo! ❌
original: C=CC1C=CC2=CC1C=C2
tested: C=CC1C=CC2=CC1C=C2

Hf(300K) S(300K) Cp(300K) Cp(400K) Cp(500K) Cp(600K) Cp(800K) Cp(1000K) Cp(1500K)
83.22 82.78 35.48 45.14 53.78 61.40 73.58 82.20 95.08
83.22 84.16 35.48 45.14 53.78 61.40 73.58 82.20 95.08

Identical thermo comments:
thermo: Thermo group additivity estimation: group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cs-(Cds-Cds)(Cds-Cds)CsH) + group(Cds-Cds(Cds-Cds)(Cds-Cds)) + group(Cds- CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-CdsCsH) + group(Cds-Cds(Cds-Cds)H) + group(Cds-Cds(Cds-Cds)H) + group(Cds-CdsHH) + Estimated bicyclic component: polycyclic(s3_5_6_ane) - ring(Cyclohexane) - ring(Cyclopentane) + ring(1,3-Cyclohexadiene) + ring(Cyclopentadiene)

Observables Test Case: Aromatics Comparison

✅ All Observables varied by less than 0.500 on average between old model and new model in all conditions!

aromatics Passed Observable Testing ✅

Regression test liquid_oxidation:

Reference: Execution time (DD:HH:MM:SS): 00:00:02:22
Current: Execution time (DD:HH:MM:SS): 00:00:02:25
Reference: Memory used: 2876.39 MB
Current: Memory used: 2855.80 MB

liquid_oxidation Passed Core Comparison ✅

Original model has 37 species.
Test model has 37 species. ✅
Original model has 241 reactions.
Test model has 241 reactions. ✅

liquid_oxidation Failed Edge Comparison ❌

Original model has 214 species.
Test model has 214 species. ✅
Original model has 1590 reactions.
Test model has 1593 reactions. ❌
The tested model has 3 reactions that the original model does not have. ❌
rxn: CC(CC(C)OO)O[O](90) + CC(CCCOO)O[O](108) <=> oxygen(1) + CC([O])CC(C)OO(110) + CC([O])CCCOO(122) origin: Peroxyl_Disproportionation
rxn: CC(CC(C)OO)O[O](90) + CC(CCCOO)O[O](108) <=> oxygen(1) + CC(=O)CC(C)OO(100) + CC(O)CCCOO(152) origin: Peroxyl_Termination
rxn: CC(CC(C)OO)O[O](90) + CC(CCCOO)O[O](108) <=> oxygen(1) + CC(=O)CCCOO(115) + CC(O)CC(C)OO(143) origin: Peroxyl_Termination

Non-identical kinetics! ❌
original:
rxn: CCC(CC)O[O](37) + CCCCCO[O](36) <=> oxygen(1) + CCC([O])CC(67) + CCCCC[O](69) origin: Peroxyl_Disproportionation
tested:
rxn: CCC(CC)O[O](37) + CCCCCO[O](35) <=> oxygen(1) + CCC([O])CC(69) + CCCCC[O](67) origin: Peroxyl_Disproportionation

k(1bar) 300K 400K 500K 600K 800K 1000K 1500K 2000K
k(T): 3.54 4.28 4.73 5.02 5.39 5.62 5.91 6.06
k(T): 7.83 7.49 7.23 7.02 6.68 6.42 5.95 5.61

kinetics: Arrhenius(A=(3.2e+12,'cm^3/(mol*s)'), n=0, Ea=(4.064,'kcal/mol'), T0=(1,'K'), comment="""Estimated from node Root_Ext-5R-R_7R!H->C_N-7C-inRing_Ext-5R-R""")
kinetics: Arrhenius(A=(3.18266e+20,'cm^3/(mol*s)'), n=-2.694, Ea=(0,'kcal/mol'), T0=(1,'K'), comment="""Estimated from node Root_Ext-5R-R_7R!H->C_N-7C-inRing""")
kinetics: Estimated from node Root_Ext-5R-R_7R!H->C_N-7C-inRing_Ext-5R-R
kinetics: Estimated from node Root_Ext-5R-R_7R!H->C_N-7C-inRing

Observables Test Case: liquid_oxidation Comparison

✅ All Observables varied by less than 0.100 on average between old model and new model in all conditions!

liquid_oxidation Passed Observable Testing ✅

Regression test nitrogen:

Reference: Execution time (DD:HH:MM:SS): 00:00:01:28
Current: Execution time (DD:HH:MM:SS): 00:00:01:28
Reference: Memory used: 2862.44 MB
Current: Memory used: 2858.19 MB

nitrogen Passed Core Comparison ✅

Original model has 41 species.
Test model has 41 species. ✅
Original model has 360 reactions.
Test model has 360 reactions. ✅

nitrogen Failed Edge Comparison ❌

Original model has 133 species.
Test model has 133 species. ✅
Original model has 983 reactions.
Test model has 983 reactions. ✅

Non-identical thermo! ❌
original: O1[C]=N1
tested: O1[C]=N1

Hf(300K) S(300K) Cp(300K) Cp(400K) Cp(500K) Cp(600K) Cp(800K) Cp(1000K) Cp(1500K)
116.46 53.90 11.62 12.71 13.49 13.96 14.14 13.85 13.58
141.64 58.66 12.26 12.27 12.09 11.96 12.26 12.72 12.15

thermo: Thermo group additivity estimation: group(O2s-CdN3d) + group(N3d-OCd) + group(Cd-HN3dO) + ring(Cyclopropene) + radical(CdJ-NdO)
thermo: Thermo group additivity estimation: group(O2s-CdN3d) + group(N3d-OCd) + group(Cd-HN3dO) + ring(oxirene) + radical(CdJ-NdO)

Non-identical kinetics! ❌
original:
rxn: NCO(66) <=> O1[C]=N1(126) origin: Intra_R_Add_Endocyclic
tested:
rxn: NCO(66) <=> O1[C]=N1(126) origin: Intra_R_Add_Endocyclic

k(1bar) 300K 400K 500K 600K 800K 1000K 1500K 2000K
k(T): -49.54 -33.65 -24.16 -17.85 -10.01 -5.35 0.80 3.82
k(T): -66.25 -46.19 -34.19 -26.21 -16.28 -10.36 -2.54 1.31

kinetics: Arrhenius(A=(6.95187e+18,'s^-1'), n=-1.628, Ea=(88.327,'kcal/mol'), T0=(1,'K'), comment="""Estimated from node Backbone0_N-2R!H-inRing_N-1R!H-inRing_Sp-2R!H-1R!H""")
kinetics: Arrhenius(A=(6.95187e+18,'s^-1'), n=-1.628, Ea=(111.271,'kcal/mol'), T0=(1,'K'), comment="""Estimated from node Backbone0_N-2R!H-inRing_N-1R!H-inRing_Sp-2R!H-1R!H""")
Identical kinetics comments:
kinetics: Estimated from node Backbone0_N-2R!H-inRing_N-1R!H-inRing_Sp-2R!H-1R!H

Observables Test Case: NC Comparison

✅ All Observables varied by less than 0.200 on average between old model and new model in all conditions!

nitrogen Passed Observable Testing ✅

Regression test oxidation:

Reference: Execution time (DD:HH:MM:SS): 00:00:02:29
Current: Execution time (DD:HH:MM:SS): 00:00:02:28
Reference: Memory used: 2738.84 MB
Current: Memory used: 2718.27 MB

oxidation Passed Core Comparison ✅

Original model has 59 species.
Test model has 59 species. ✅
Original model has 694 reactions.
Test model has 694 reactions. ✅

oxidation Passed Edge Comparison ✅

Original model has 230 species.
Test model has 230 species. ✅
Original model has 1526 reactions.
Test model has 1526 reactions. ✅

Observables Test Case: Oxidation Comparison

✅ All Observables varied by less than 0.500 on average between old model and new model in all conditions!

oxidation Passed Observable Testing ✅

Regression test sulfur:

Reference: Execution time (DD:HH:MM:SS): 00:00:00:56
Current: Execution time (DD:HH:MM:SS): 00:00:00:56
Reference: Memory used: 2834.41 MB
Current: Memory used: 2824.60 MB

sulfur Passed Core Comparison ✅

Original model has 27 species.
Test model has 27 species. ✅
Original model has 74 reactions.
Test model has 74 reactions. ✅

sulfur Failed Edge Comparison ❌

Original model has 89 species.
Test model has 89 species. ✅
Original model has 227 reactions.
Test model has 227 reactions. ✅
The original model has 1 reactions that the tested model does not have. ❌
rxn: O(4) + SO2(15) (+N2) <=> SO3(16) (+N2) origin: primarySulfurLibrary
The tested model has 1 reactions that the original model does not have. ❌
rxn: O(4) + SO2(15) (+N2) <=> SO3(16) (+N2) origin: primarySulfurLibrary

Observables Test Case: SO2 Comparison

✅ All Observables varied by less than 0.100 on average between old model and new model in all conditions!

sulfur Passed Observable Testing ✅

Regression test superminimal:

Reference: Execution time (DD:HH:MM:SS): 00:00:00:39
Current: Execution time (DD:HH:MM:SS): 00:00:00:39
Reference: Memory used: 2923.50 MB
Current: Memory used: 2893.71 MB

superminimal Passed Core Comparison ✅

Original model has 13 species.
Test model has 13 species. ✅
Original model has 21 reactions.
Test model has 21 reactions. ✅

superminimal Passed Edge Comparison ✅

Original model has 18 species.
Test model has 18 species. ✅
Original model has 28 reactions.
Test model has 28 reactions. ✅

Regression test RMS_constantVIdealGasReactor_superminimal:

Reference: Execution time (DD:HH:MM:SS): 00:00:02:30
Current: Execution time (DD:HH:MM:SS): 00:00:02:26
Reference: Memory used: 3471.36 MB
Current: Memory used: 3448.94 MB

RMS_constantVIdealGasReactor_superminimal Passed Core Comparison ✅

Original model has 13 species.
Test model has 13 species. ✅
Original model has 19 reactions.
Test model has 19 reactions. ✅

RMS_constantVIdealGasReactor_superminimal Passed Edge Comparison ✅

Original model has 13 species.
Test model has 13 species. ✅
Original model has 19 reactions.
Test model has 19 reactions. ✅

Observables Test Case: RMS_constantVIdealGasReactor_superminimal Comparison

✅ All Observables varied by less than 0.100 on average between old model and new model in all conditions!

RMS_constantVIdealGasReactor_superminimal Passed Observable Testing ✅

Regression test RMS_CSTR_liquid_oxidation:

Reference: Execution time (DD:HH:MM:SS): 00:00:06:25
Current: Execution time (DD:HH:MM:SS): 00:00:06:15
Reference: Memory used: 3447.05 MB
Current: Memory used: 3466.62 MB

RMS_CSTR_liquid_oxidation Passed Core Comparison ✅

Original model has 37 species.
Test model has 37 species. ✅
Original model has 202 reactions.
Test model has 202 reactions. ✅

RMS_CSTR_liquid_oxidation Passed Edge Comparison ✅

Original model has 248 species.
Test model has 248 species. ✅
Original model has 2057 reactions.
Test model has 2057 reactions. ✅

Observables Test Case: RMS_CSTR_liquid_oxidation Comparison

✅ All Observables varied by less than 0.100 on average between old model and new model in all conditions!

RMS_CSTR_liquid_oxidation Passed Observable Testing ✅

Regression test fragment:

Reference: Execution time (DD:HH:MM:SS): 00:00:00:42
Current: Execution time (DD:HH:MM:SS): 00:00:00:42
Reference: Memory used: 2661.14 MB
Current: Memory used: 2665.01 MB

fragment Passed Core Comparison ✅

Original model has 10 species.
Test model has 10 species. ✅
Original model has 2 reactions.
Test model has 2 reactions. ✅

fragment Passed Edge Comparison ✅

Original model has 33 species.
Test model has 33 species. ✅
Original model has 47 reactions.
Test model has 47 reactions. ✅

Observables Test Case: fragment Comparison

✅ All Observables varied by less than 0.100 on average between old model and new model in all conditions!

fragment Passed Observable Testing ✅

Regression test RMS_constantVIdealGasReactor_fragment:

Reference: Execution time (DD:HH:MM:SS): 00:00:03:09
Current: Execution time (DD:HH:MM:SS): 00:00:03:04
Reference: Memory used: 3631.22 MB
Current: Memory used: 3623.94 MB

RMS_constantVIdealGasReactor_fragment Passed Core Comparison ✅

Original model has 10 species.
Test model has 10 species. ✅
Original model has 2 reactions.
Test model has 2 reactions. ✅

RMS_constantVIdealGasReactor_fragment Passed Edge Comparison ✅

Original model has 27 species.
Test model has 27 species. ✅
Original model has 24 reactions.
Test model has 24 reactions. ✅

Observables Test Case: RMS_constantVIdealGasReactor_fragment Comparison

✅ All Observables varied by less than 0.100 on average between old model and new model in all conditions!

RMS_constantVIdealGasReactor_fragment Passed Observable Testing ✅

Regression test minimal_surface:

Reference: Execution time (DD:HH:MM:SS): 00:00:00:45
Current: Execution time (DD:HH:MM:SS): 00:00:00:45
Reference: Memory used: 2836.48 MB
Current: Memory used: 2833.31 MB

minimal_surface Passed Core Comparison ✅

Original model has 11 species.
Test model has 11 species. ✅
Original model has 3 reactions.
Test model has 3 reactions. ✅

minimal_surface Passed Edge Comparison ✅

Original model has 38 species.
Test model has 38 species. ✅
Original model has 38 reactions.
Test model has 38 reactions. ✅

Observables Test Case: minimal_surface Comparison

✅ All Observables varied by less than 0.500 on average between old model and new model in all conditions!

minimal_surface Passed Observable Testing ✅

beep boop this comment was written by a bot 🤖

@donerancl
Copy link
Contributor Author

donerancl commented May 13, 2025

@JacksonBurns @alongd @mjohnson541 would you be willing to review this PR? It is rather straightforward if we are willing to not apply energy_correction in RMG's pdep mode. This correction is only intended to change the reference E0 and only the relative E0 is important, so I think it is not necessary.

@mjohnson541
Copy link
Contributor

My understanding is that the intent of this correction is exclusively so our energy diagrams are in the vicinity of zero. The choice of the correction is arbitrary, but intended to adjust the y-axis location in energy diagrams while keeping relative energies the same.

Generally, I hate having this done in the vicinity of legitimate physical calculations and I think if we're going to shift the diagram we should do that when we draw the diagram. However, historically other people have had much stronger opinions about our energy diagrams than I do. So I would wait for @alongd's opinion.

I don't mind getting rid of it altogether, but adding the species energy correction back in might be the easiest way to restore the desired functionality.

@rwest
Copy link
Member

rwest commented May 14, 2025

Oh, boy, now I feel bad. I came across this problem (I think) back in December, and was half way to opening an issue and/or pull request, then got distracted.

I was debugging a PDep network and as commonly done I we would find the networkX_Y.py file and make Arkane draw the PDF. Getting tired of this in in 9a6bd59 I made RMG save the PDF of the network automatically as it crashes. Strangely, the two outputs were different!
Something like these two.

Screenshot 2025-05-14 at 2 24 41 PM

I got about as far as @donerancl in realizing it's the inconsistently applied "correction".

I think I determined that during running of the RMG job, the correction was correctly applied, but it was inconsistently (and incorrectly) saved in the arkane input files. The conformers would have the correction applied, but not the TS, or something like that.

I kept the PDFs and windows and tabs open for about a month hoping to get around to finishing my investigation, but eventually had to reboot the laptop. I have one more stashed commit from February when I changed branch to do something else, where around line 26 of rmgpy/pdep/configuration.pyx I added

        if self.energy_correction != 0.0: string += 'energy_correction={0}, '.format(self.energy_correction)

I can't remember if that fixed it.

Based on when David (I think) added the correction, I seem to recall it helped numerical issues (you end up taking exponentials of this energy when evaluating a partition function?) and was not just to make the PES figures look nice.
That gives me reason to pause before approving the current PR. We should check several scenarios to see if the results change and it's still numerically stable and doesn't give overflows.
eg. the commit message ot 169a34e.
see also 7f26746

@rwest
Copy link
Member

rwest commented May 15, 2025

Would this be a good solution?

  1. do not store an energy_correction. remove the current one entirely
  2. just before evaluating equilibrium constants, partition functions, etc., whatever needs energy values and was running into overflow errors, re-normalize or shift energies at that point, but don't save the results.
  3. just before drawing PES diagrams, re-normalize or shift energies at that point, perhaps so the entrance channel, or the lowest well, is zero, but don't save the results.

I'm slightly worried that point 2 would be prone to bugs. I haven't looked.

@mjohnson541
Copy link
Contributor

I think that's definitely the most readable way to handle this, but it does require identifying all locations where this might happen and some of the algorithmic fixes may be much more than couple line fixes. I suspect we have object stored vectors that may become 0/Inf/NaN without the correction which may require more involved refactoring.

It would be good to fix that, but given the benefits/work cost perhaps its most practical now to ensure the correction is applied to all the energies correctly and add some detailed comments explaining what is going on.

@donerancl
Copy link
Contributor Author

@rwest @mjohnson541 Thanks for your discussion -- just returned from vacation. I think you are right that we need to find everywhere the energy_correction may be applied. I used this pull request branch to generate a mechanism for NH3 combustion and the RMG-Py run did fine, but I tried to run some of the network.py files in the pdep directory with Arkane and ran into many numerical errors.

Also, somewhat related, I am noticing that many of the networks from the network.py could be merged (share isomers), which makes me think RMG isn't updating its networks as intended when new species/reactions are added. I'm still investigating what's happening here but should have another PR coming in the near future.

@donerancl
Copy link
Contributor Author

notes for my reference on places where energy_correction is applied...
in rmgpy/rmg/pdep.py

    def update_configurations(self, reaction_model):
        """
        Sort the reactants and products of each of the network's path reactions
        into isomers, reactant channels, and product channels. You must pass 
        the current `reaction_model` because some decisions on sorting are made
        based on which species are in the model core. 
        """
        ...
        if self.energy_correction:
            for spec in self.reactants + self.products + self.isomers:
                spec.energy_correction = self.energy_correction

in rmgpy/pdep/configuration.pyx

cpdef class Configuration(object):
    ...
    def __init__(self, *species):
         ...
        self.energy_correction = 0.0
    ...
    property E0:
        """The ground-state energy of the configuration in J/mol."""
        def __get__(self):
            return sum([float(spec.conformer.E0.value_si) for spec in self.species]) + self.energy_correction 
    ...
    cpdef double get_enthalpy(self, double T) except 100000000:
        """
        Return the enthalpy in kJ/mol at the specified temperature `T` in K.
        """
        cdef double h = 0.0
        for spec in self.species:
            h += spec.get_enthalpy(T)
        h += self.energy_correction
        return h
    ...
    cpdef double get_free_energy(self, double T) except 100000000:
        """
        Return the Gibbs free energy in kJ/mol at the specified temperature
        `T` in K.
        """
        cdef double g = 0.0
        for spec in self.species:
            g += spec.get_free_energy(T)
        g += self.energy_correction
        return g        

in rmgpy/pdep/configuration.pxd

cdef class Configuration(object):

    cdef public list species
    cdef public np.ndarray e_list
    cdef public np.ndarray dens_states
    cdef public np.ndarray sum_states
    cdef public bint active_j_rotor
    cdef public bint active_k_rotor
    cdef public float energy_correction

in arkane/explorer.py

    def execute(self, output_file, plot, file_format='pdf', print_summary=True, species_list=None,
                thermo_library=None, kinetics_library=None):
        """Execute an ExplorerJob"""
        ...
        for network in self.networks:
            for rxn in network.path_reactions:
                if rxn.transition_state is None:
                    E0 = sum([spec.conformer.E0.value_si for spec in rxn.reactants])
                    if rxn.network_kinetics is None:
                        E0 += rxn.kinetics.Ea.value_si
                    else:
                        E0 += rxn.network_kinetics.Ea.value_si
                    if network.energy_correction:
                        E0 += network.energy_correction
                    rxn.transition_state = rmgpy.species.TransitionState(conformer=Conformer(E0=(E0 * 0.001, "kJ/mol")))
       

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants