Blowers Masel Arrhenius trained and evaluated with T-varying ∆Hr(T) ?

[rwest]

Bug Description

Here

RMG-Py/rmgpy/data/kinetics/family.py

Line 3502 in 3e37656

kinetics = kinetics.to_arrhenius(rxn.get_enthalpy_of_reaction(T))

an ArrheniusBM is converted to an Arrhenius using the ∆Hrxn evaluated at the temperature of interest, something like 1000K, rather than at 298K.

Looking at the code used to fit the BM expressions here

RMG-Py/rmgpy/kinetics/arrhenius.pyx

Line 616 in 91cdd1d

xdata.append([T, rxn.get_enthalpy_of_reaction(T)])

perhaps this is correct, or at least consistent (we shouldn't change one without changing both.)

Skimming the Blowers and Masel paper https://doi.org/10.1002/aic.690461015 I can't see it explicitly say what temperature it's evaluated at (it just says "∆Hr is the heat of reaction", and later "[we have] looked up the heats of reaction [for the 151 reactions in Table 1]") but given that it doesn't talk about barrier heights changing at different temperatures, just barrier heights, it doesn't make sense to me for ∆Hr to be a function of temperature. If Ea is f(∆Hr) and Ea independent of T then ∆Hr is independent of T. I would guess they looked up the 298K value, but as they don't actually give any of the values in the paper I'm not sure how to check. If we're re-training our own parameters then probably its up to us anyway - as long as we're consistent.

To my knowledge (and ability to search) all other calls of to_arrhenius in RMG, whether called on instances of ArrheniusEP, SurfaceArrheniusBEP, StickingCoefficientBEP, or ArrheniusBM, all use ∆Hr298.

My gut says ∆Hr shouldn't be a function of Temperature, and so the BM fitting and evaluating segments linked above are wrong. But I've run out of time to dig deeper.

[mjohnson541]

My understanding is that ∆Hr is a relatively weak function of temperature. When fitting we evaluate ∆Hr(T) to be thorough in fitting, but we don't lose much accuracy when we convert to Arrhenius using ∆Hr(T=1000K)(T=1000K because that's what the trees are currently fitted at).

[rwest]

As we set about implementing Blowers Masel rate expressions in Cantera (Cantera/cantera#987 ) we have come across this issue again. Should the Blowers Masel expression be evaluated using ∆Hr298 or ∆Hr(T) ?

Yes, Matt's right ∆Hrxn is a pretty weak function of T, but when @12Chao converts all H abstractions (except those with negative barriers) into Blowers Masel expressions, fitted at 298K but evaluated at T, the ignition delay of n-heptane changes about 3%. So not very big, but not zero.

I guess we have a few questions. 1) what did Blowers and Masel intend? 2) what actually works best? 3) what will users expect? (we have chance to influence this, as for now we are the only users).

(Blowers' PhD thesis is at http://hdl.handle.net/2142/82466 but restricted access)

[mjohnson541]

For (1) the paper (https://aiche.onlinelibrary.wiley.com/doi/abs/10.1002/aic.690461015) almost comically doesn't even have the word temperature in the whole document. In general it talks about dHrxn like it isn't temperature dependent. As far as what they're actually using they say: "In order to do so we have taken all of the hydrogen transfer reactions listed in the NIST tables Westley 1980., and looked up the heats of reaction." Which makes me think they're using dHrxn298 or dHrxn0.

[mjohnson541]

One thing I particularly don't like about evaluating the BM using dHrxn(T) is that it causes programmatic branching in the rate evaluation which makes it more expensive to differentiate in the reverse-mode and much harder to generate symbolic representations for. Of course NASA polynomials have the same issue, but in an ideal world I'd like to move away from them as well.

[rwest]

Another problem with making Ea a function of ∆Hrxn(T) is that you cannot recover an Arrhenius expression. Ea is not just a function of the reaction. I think the ∆H should be a function of other state variables: concentration (in the case of non-ideal solvents), surface coverage (for co-adsorbates in heterogeneous catalysis), etc. so that these effects that we have learned to predict for thermochemistry are also altering kinetics. Which might also make automatic differentiation tricky. But I think I'm veering towards saying it shouldn't be a function of T. (hmm. but what about temperature-dependent solvation? or crossing into a supercritical regime? there are all sorts of counterexamples and problems)

But in the context of Cantera, for example, where you have a manager keeping track of species enthalpies, that is already at the state Temperature and it's in easier to query the current ∆Hrxn(T) than a reference ∆Hrxn.

[mjohnson541]

So I think we might want to think about how we format things this way:

RMG est format => RMG output => Simulator format

We use RMG to generate the rate/thermochemistry estimates in some programmatic format, we write the RMG mechanism to a file in some format and when the simulator reads the file it compiles it into a third format.

In the first case we're generating the estimates so the format is whatever our estimators spit out and we really only care that the format allows the estimators to perform well. In the simulator format we care that we've represented the estimation format properly, but we also care a lot about our ability to evaluate it optimally. The RMG output can be anything as long as it has enough information to generate simulator format.

Following the solvation example. We could have RMG output that is solvent specific. We could also output something much more similar to the RMG est format and let the simulator interpret that to whatever format it thinks is most convenient. These put different burdens on different software, but I think neither of these prevents the simulator from compiling its own kinetics format at the end.

[rwest]

@sevyharris and I just spent a while getting very confused as to why our Blowers-Masel fit was not fitting the training reaction when evaluated with H(T) but only with H(1000K). Turns out it's because we had a single training reaction, and this line

RMG-Py/rmgpy/kinetics/arrhenius.pyx

Line 601 in 336273d

T = 1000.0

added in #1978
shows that a single training reaction is fitted using H(1000K) not H(T).

Having been thinking about this quite a bit (in the context of adding BM to Cantera) I'm leaning towards the conclusion that Blowers and Masel intended that it should be H(298K) not H(T). In any case, I think we should be consistent!

[github-actions]

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.