ArrheniusBM fitting procedure

[davidfarinajr]

Motivation or Problem

• ArrBM does not properly handle case when E0 < 0
• ArrBM does not fit to single reaction where dHrxn > 4 * w0/10 or dHrxn < -4 * w0/10
  • This is particularly a problem for exothermic reactions with large barriers (Ea from fitted rule is 0, but Ea in training reaction is >>0)

Description of Changes

• modified get_activation_energy to handle E0 < 0
  • returns DHrxn for endothermic and E0 for exothermic reactions
• modified arrbm fitting to rxns procedure
  • use get_activation_energy method when getting training rates and fitting params
  • get dHrxn only at 298K since we only use enthalpy at 298K to convert ArrBM to Arr
  • use same procedure regardless of number of training rxns
  • use avg of params (with BEP to guess E0) as initial guesses

Testing

I used this to regenerate rules for a bunch of families. It doesn't change too much, except for rules trained on single rxn (if rxn is very exothermic with high Ea)

Reviewer Tips

• Test generating BM rules with tree gen notebook
• discussion needed: not all changes in this PR are necessary. Important parts are revised get_activation_energy method and fitting to single training reactions.

[davidfarinajr]

I modified the ArrBM fit_to_data test (4efedec) to test the ArrBM fit for this reaction:

entry(
    index = 54,
    label = "CF2_r1 + C2H6 <=> C3H6F2",
    degeneracy = 1.0,
    kinetics = Arrhenius(A=(0.222791,'cm^3/(mol*s)'), n=3.59921, Ea=(320.496,'kJ/mol'), T0=(1,'K'), Tmin=(298,'K'), Tmax=(2500,'K'), comment="""Fitted to 50 data points; dA = *|/ 1.26384, dn = +|- 0.0307635, dEa = +|- 0.167414 kJ/mol"""),
    rank = 7,
    shortDesc = """M062X-D3/jun-cc-pVTZ RRHO""",
    longDesc = 
"""
Calculated with Gaussian 16 using M062X with D3 dispersion and jun-cc-pVTZ basis set
barrier = 332.008 kJ/mol

C    -0.695038    -0.080264    -0.119021
F    1.66398    -0.643084    -0.51829
F    1.066759    -0.486589    1.54559
H    0.945151    2.01382    -0.733224
H    -0.752053    1.947025    -0.291092
H    0.450849    2.144078    0.971957
C    1.305838    0.231499    0.436189
C    0.261066    1.650533    0.025447
H    -0.964565    0.081354    -1.15545
H    -1.466846    0.212255    0.582107
H    -0.37114    -1.100917    0.051497
""",
)

dHrxn = -276 kJ/mol
w0 = 519 kJ
w0/10 * 4 = 207.6 kJ/mol 
abs(dHrxn) > 4 * w0 / 10.0:
     E0 = w0 / 10.0

On the main branch, the test fails with this error

======================================================================
FAIL: test_fit_to_data (__main__.TestArrheniusBM)
Test the ArrheniusBM.fit_to_data() method.
----------------------------------------------------------------------
Traceback (most recent call last):
  File "arrheniusTest.py", line 567, in test_fit_to_data
    self.assertAlmostEqual(k, arrhenius.get_rate_coefficient(T), delta=1e-6 * k)
AssertionError: 2.893471124690076e-54 != 183.475142650761 within 2.893471124690076e-60 delta (183.475142650761 difference)

----------------------------------------------------------------------

The test passes on PR branch

[rwest]

I've just rebased this onto the official main.

I'm not sure why we'd ever have a negative E0, based on Blowers and Masel's plots? Doesn't seem to make sense. But I perhaps it results from fitting it to certain training reactions with submerged barriers?

[mjohnson541]

I have some similar fixes on the electrochem branch. After discussion we decided to simply average reactions instead of using the BM fit when E0<0.

[mjohnson541]

We decided that the BM fit shouldn't be expected to be valid for submerged barrier reactions.

[rwest]

Is it these ones, Matt? #2277

[codecov]

Codecov Report

Merging #2264 (e80a8ff) into main (e6a0040) will not change coverage.
The diff coverage is n/a.

@@           Coverage Diff           @@
##             main    #2264   +/-   ##
=======================================
  Coverage   48.14%   48.14%           
=======================================
  Files         110      110           
  Lines       30629    30629           
  Branches     7989     7989           
=======================================
  Hits        14747    14747           
  Misses      14353    14353           
  Partials     1529     1529           

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

[rwest]

I think this may address #1748 but this PR may be superseded by #2277 (and is likely incompatible with it).

[rwest]

I think #2277 (which supersedes this) is now part of #2316. Hopefully that is merged soon.

[JacksonBurns]

Closing in favor of #2316 for the time being, but if that doesn't pan out we can come back to this one.