k(T,P) values from partial networks don't always satisfy thermo consistence

[jwallen]

We are using the "activated species algorithm" of Matheu, thinking that many small partial networks will be faster in the long run than fewer large partial networks. One tradeoff of this is that the k(T,P) values for the forward and reverse directions come from different partial networks (and are treated as irreversible), so there is no guarantee that thermodynamic consistency is satisfied, although they should be close.

An example of a case where they are not close appears early in the 1,3-hexadiene example, and involves the reaction

C=[CH](13) + [CH]=C[CH2](50) <=> C=CC=C[CH2](5)

As expected two partial networks are created containing this path reaction: one with C=[CH](13) + [CH]=C[CH2](50) as the source and one with C=CC=C[CH2](5) as the source. The two partial networks (123 and 16, respectively) return k(1000 K, 1 bar) = 5.8e13 cm^3/mol_s and 2.0e27 cm^3/mol_s for the forward direction as given above (note that the latter was computed via the equilibrium constant). The former is close to the high-P rate of 7.23e13 cm^3/mol*s, while the latter exceeds it by 14 orders of magnitude.

[rwest]

Hmm. 10^14 times faster is indeed "not close"!

[jwallen]

The two relevant partial networks contain all of the same path reactions. The only difference is in how they treat C=[CH](13) + [CH]=C[CH2](50); in 123 they are a reactant channel, while in 16 they are a product channel. The consequence of this is that the two networks compute k(E) for this reaction in different ways; in 123 the ILT method is applied to the forward kinetics (since we know the frequencies of the reactant channel), while in 16 the kinetics are first reversed and the ILT method applied to the reverse kinetics.

To fix the thermodynamic consistence issue, we should always compute k(E) in the same way for each partial network. Ultimately I think it would be better to always apply the ILT method in the direction the kinetics are known. Often this will be the exothermic direction, which generally has weaker T dependence and is therefore safer to apply the ILT method to. (Note that this is different from RRKM theory, which must always be applied to situations with unimolecular reactants.)

[jwallen]

Incidentally, I wonder if we are getting an accurate estimate for the thermo of [CH]=C[CH2](50). RMG-Py seems to think the most stable resonance form is actually [CH]C=C, by 6 kcal/mol over [CH]=C[CH2].

[rwest]

Always computing k(E) the same way would ensure both directions gave the same k(T,P) and solve the thermo consistence issue, and sounds like a good idea, but I can't help feeling that 10^14 is such a large error it must be a bug, and by doing as you suggest we just hide the bug. As you said - the two directions ought to be close. Can 10^14 be explained some non-buggy way? Is the method really that sensitive to our [poor] estimates? (If so, how can we use it?)

[jwallen]

In this particular case, the issue seems to be one of barrier heights vs. activation energies:

• In the forward direction, the high-P kinetics are kf(T) = 7.23e13 cm^3/mol*s. Note that the activation energy is 0 kJ/mol, which makes sense for an R_Recombination reaction.
• Fitting the reverse high-P kinetics by computing kr at several temperatures and fitting an Arrhenius expression gives kr(T) = 4.9e17 * T^0.47 * exp(-493 kJ/mol / RT) s^-1.
• RMG predicts that the enthalpy of reaction for this reaction is -497 kJ/mol at 298 K, consistent with the fitted Ea for the reverse kinetics.
• However, the ground-state energy for the minima on the PES are computed at (or at least near) 0 K, for which the energy difference is only -200 kJ/mol.

To make matters more interesting, the inverse Laplace transform method assumes that the activation energy is equivalent to the barrier height. Clearly this makes a huge difference.

[jwallen]

Aha, it is in fact a bug. In particular, the ground-state energy E0 for C=[CH] was not being set because its thermo came from GRI-Mech and was already a MultiNASA. The fix is to also set E0 for this case, by converting back to Wilhoit just long enough to extrapolate to near 0 K.

The fixes I've made toward computing k(E) the same way in both partial networks should still be used, I think, but they are not the source of a thermo inconsistency of 14 orders of magnitude.

[rwest]

I'm hitting the new exception you introduced to check for inconsistencies. Seem to be off by ~2.5

CO2(8) + H(17) -> OjCHO(55)
877.012194043 310.991929235
OjCHO(55) -> CO2(8) + H(17)
154662008113.0 54843748592.3

rmgpy.rmg.pdep.PressureDependenceError: Forward and reverse PDepReactions generated from networks 256 and 819 do not satisfy thermodynamic consistency.
  File "/home/rwest/code/RMG-Py/rmgpy/rmg/model.py", line 667, in enlarge

Networks: http://rmg.mit.edu/measure/networks/2475ef5c78b9aa3909320d397b8b47cc http://rmg.mit.edu/measure/networks/96863605511aa70a4fa0af5d6309bbf2

[jwallen]

The only difference between these two networks is that, in the first, CO2(8) + H(17) is treated as a product channel, while in the second it is a reactant channel. Both partial networks estimate k(E) for the forward and reverse directions in exactly the same way; however, the first discards the association direction (because product channels are treated as infinite sinks). This causes the k(T,P) value for OjCHO(55) -> CO2(8) + H(17) from the first network to be a (slight but significant) overestimate because, in that network, CO2(8) + H(17) can never reassociate back to OjCHO(55).

The solution, therefore, is to never treat configurations consisting of all core species as product channels. I should have a commit for this pushed shortly.

[rwest]

Now off by ~70:

CO2(8) + H(17) -> OjCHO(55)
557.623081103 7.85591671697
OjCHO(55) -> CO2(8) + H(17)
98337407483.4 1385399043.79
rmgpy.rmg.pdep.PressureDependenceError: Forward and reverse PDepReactions generated from networks 256 and 819 do not satisfy thermodynamic consistency.
  File "/home/rwest/code/RMG-Py/rmgpy/rmg/model.py", line 667, in enlarge

The measure input files for networks 256 and 819 are identical, so current guess is that the results are the same but the Chebyshev fit is bad, at least at the T,P at which the thermo equivalence is tested.

[jwallen]

The issue is one of barrier height vs. activation energy. The reaction CO2(8) + H(17) -> OjCHO(55) is originally estimated to have an activation energy of 13.8 kJ/mol. However, CO2 is very stable, so the above reaction is actually endothermic by 53.3 kJ/mol at 298 K. We were already adjusting the activation energy to 53.3 kJ/mol (a 40 kJ/mol change!) because we believe that we are overpredicting the rate of k(T) otherwise.

However, the approximate E0 values for the reactants and products say that the reaction is endothermic by 59.6 kJ/mol at ~0 K, so even the adjusted barrier was still negative from the point of view of the pressure dependence algorithm. This caused the computed k(E) values in both networks 256 and 819 to be off from the macroscopic equilibrium enought to trigger the thermodynamic consistency error.

The short-term fix, which I am currently testing, is simply to adjust the barrier using the enthalpy of reaction at 0 K instead of 298 K. Not sure about the long-term fix; probably need to consult with Bill.

[jwallen]

Still having this issue. From 1,3-hexadiene example (networks 95 and 321):

C(=C[CH2])C=C(5) <=> [H](10) + C1C=CC=C1(312)
334.649706171 5.77219991874
[H](10) + C1C=CC=C1(312) -> C(=C[CH2])C=C(5)
177248.009737 3057.25936265

Network #95 (source is C(=C[CH2])C=C(5)): http://rmg.mit.edu/measure/networks/ae545df2a0c5d908cb22ba185db91170 Network #321 (source is [H](10) + C1C=CC=C1(312)): http://rmg.mit.edu/measure/networks/8b52d42aa191455f4b15830a707c6bac

These networks are different in that they have a different number of explored isomers.

[jwallen]

I'm no longer convinced that we should expect the forward and reverse directions to be "close" to thermodynamic consistence. In the above example, the offending net reaction is actually a well-skipping reaction, and in competition with many other well-skipping reactions as well as collisional stabilization. If one or more of these competing reactions is to an infinite sink, then the rate of reaction to that sink will be overestimated (perhaps greatly so if the reaction is uphill in energy), if for no other reason than the negelecting of the reverse process. This means that the well-skipping reaction we are interested in will be slower than otherwise expected (perhaps significantly so). This is what seems to be happening in the most recent system.

My current proposed solution is to forget the thermo consistence check, and instead keep the net reaction from the network that is more extensively explored. If both networks are explored similarly, then keep the net reaction that is faster -- not necessarily the one in the exothermic direction, but the one that predicts the higher k(T, P) when comparing in the same direction. Again, need to discuss this with Bill.

[rwest]

Sounds good.

[rwest]

With the latest version running the methylformate example with Reservoir State pdep it ends:

Using 336 grains from -189.25 to 1212.39 kJ/mol in steps of 4.18 kJ/mol
Calculating densities of states for Network #577...
Calculating phenomenological rate coefficients for Network #577...
Error: For path reaction CH2O(14) + OH(19) <=> [CH2]OO(85):
Error:     Expected kf(292.578 K) = 3.30873e-19
Error:       Actual kf(292.578 K) = 1.79156e-18
Error:     Expected Keq(292.578 K) = 1.26848e-33
Error:       Actual Keq(292.578 K) = 3.91513e-33

Traceback (most recent call last):
  File "/home/rwest/code/RMG-Py/rmg.py", line 134, in <module>
    cProfile.runctx(command, global_vars, local_vars, stats_file)
  File "/usr/lib/python2.6/cProfile.py", line 49, in runctx
    prof = prof.runctx(statement, globals, locals)
  File "/usr/lib/python2.6/cProfile.py", line 140, in runctx
    exec cmd in globals, locals
  File "<string>", line 1, in <module>
  File "/home/rwest/code/RMG-Py/rmgpy/rmg/main.py", line 336, in execute
    self.reactionModel.enlarge(object)
  File "/home/rwest/code/RMG-Py/rmgpy/rmg/model.py", line 628, in enlarge
    self.updateUnimolecularReactionNetworks(database)
  File "/home/rwest/code/RMG-Py/rmgpy/rmg/model.py", line 1220, in updateUnimolecularReactionNetworks
    network.update(self, database, self.pressureDependence)
  File "/home/rwest/code/RMG-Py/rmgpy/rmg/pdep.py", line 440, in update
    K, p0 = self.calculateRateCoefficients(Tlist, Plist, Elist, method)
  File "/home/rwest/code/RMG-Py/rmgpy/measure/network.py", line 642, in calculateRateCoefficients
    Kij, Gnj, Fim = self.calculateMicrocanonicalRates(Elist, densStates0, T)
  File "/home/rwest/code/RMG-Py/rmgpy/measure/network.py", line 491, in calculateMicrocanonicalRates
    raise NetworkError('Invalid k(E) values computed for path reaction "{0}" do not satisfy'.format(rxn))
rmgpy.measure.network.NetworkError: Invalid k(E) values computed for path reaction "CH2O(14) + OH(19) <=> [CH2]OO(85)" do not satisfy

[github-actions[bot]]

This issue is being automatically marked as stale because it has not received any interaction in the last 90 days. Please leave a comment if this is still a relevant issue, otherwise it will automatically be closed in 30 days.