Problem with Singlet_Carbene_Intra_Disproportionation

[JLans]

Bug Description

I am finding that RMG is generating highly negative activation energies using training reactions for the reaction family 'Singlet_Carbene_Intra_Disproportionation'.

How To Reproduce

#installing relevant modules
import rmgpy
from rmgpy.molecule import Molecule
from rmgpy import settings
from rmgpy.data.rmg import RMGDatabase
from rmgpy.rmg.model import CoreEdgeReactionModel

#defining the database
database = RMGDatabase()
database.load(
    path=settings['database.directory'],
    thermo_libraries = ['primaryThermoLibrary'],
    reaction_libraries=[],
    kinetics_depositories = ['training'],
    kinetics_families=[
        'Singlet_Carbene_Intra_Disproportionation',
    ],
    depository=False,
    solvation=False,
    surface=True,
)

#Defining molecules
adjlist="""
1  C u0 p0 c0 {2,S} {7,S} {13,S} {14,S}
2  C u0 p0 c0 {1,S} {8,S} {15,S} {16,S}
3  C u0 p0 c0 {4,S} {6,S} {11,S} {12,S}
4  C u0 p0 c0 {3,S} {10,S} {17,S} {18,S}
5  C u0 p0 c0 {10,S} {19,S} {20,S} {21,S}
6  C u0 p0 c0 {3,S} {7,D} {22,S}
7  C u0 p0 c0 {1,S} {6,D} {23,S}
8  C u0 p0 c0 {2,S} {9,D} {24,S}
9  C u0 p0 c0 {8,D} {25,S} {26,S}
10 C u0 p1 c0 {4,S} {5,S}
11 H u0 p0 c0 {3,S}
12 H u0 p0 c0 {3,S}
13 H u0 p0 c0 {1,S}
14 H u0 p0 c0 {1,S}
15 H u0 p0 c0 {2,S}
16 H u0 p0 c0 {2,S}
17 H u0 p0 c0 {4,S}
18 H u0 p0 c0 {4,S}
19 H u0 p0 c0 {5,S}
20 H u0 p0 c0 {5,S}
21 H u0 p0 c0 {5,S}
22 H u0 p0 c0 {6,S}
23 H u0 p0 c0 {7,S}
24 H u0 p0 c0 {8,S}
25 H u0 p0 c0 {9,S}
26 H u0 p0 c0 {9,S}
"""
mol1 = Molecule().from_adjacency_list(adjlist) #using the smiles for molecule 1 did not produce the same adjacency list
mol2 = Molecule().from_smiles('C=CCCC=CCCC=C')

#Get the reaction family
family = database.kinetics.families['Singlet_Carbene_Intra_Disproportionation']
#generate reactions
#reaction_list = database.kinetics.generate_reactions_from_families([mol1], [mol2], only_families=['Singlet_Carbene_Intra_Disproportionation'])
reaction_list = family.generate_reactions([mol1],[mol2]) #if family is known
reaction = reaction_list[0]

#Get the template for the reaction
template = family.retrieve_template(reaction.template)
#Add rules from training reaction
family.add_rules_from_training(thermo_database=database.thermo)
#expand rules
family.fill_rules_by_averaging_up(verbose=True)
print(family.get_kinetics_for_template(template))

#Another method
model = CoreEdgeReactionModel()
reaction_RMG = model.make_new_reaction(reaction, perform_cut=False)
print(reaction_RMG[0].kinetics)

#OUTPUT
(ArrheniusEP(A=(8.17959e+16,'s^-1'), n=-1.28029, alpha=0, E0=(-20671.3,'J/mol'), Tmin=(303.03,'K'), Tmax=(2000,'K'), comment="""Average of [From training reaction 2 used for CsJ2-C;CsJ2(CsC);CH=C]
Estimated using template [CsJ2-C;CsJ2(CsC);CH] for rate rule [CsJ2-C;CsJ2(CsC);CH3]
Euclidian distance = 1.0
family: Singlet_Carbene_Intra_Disproportionation"""), None)
Arrhenius(A=(8.17959e+16,'s^-1'), n=-1.28029, Ea=(-20.6713,'kJ/mol'), T0=(1,'K'), Tmin=(303.03,'K'), Tmax=(2000,'K'), comment="""Estimated using template [CsJ2-C;CsJ2(CsC);CH] for rate rule [CsJ2-C;CsJ2(CsC);CH3]
Euclidian distance = 1.0
family: Singlet_Carbene_Intra_Disproportionation
Ea raised from -20.7 to -20.7 kJ/mol.""")

Expected Behavior

I expected the activation energies to not be such large negative numbers.

Installation Information

Describe your installation method and system information.

• installed from source without anaconda
  • RMG-Py: 3.2.0-189-ged9f29d06
  • RMG-database: 3.2.0-5-gcfb4910cf

Additional Context

The RMG kinetics database search on the https://rmg.mit.edu/database/kinetics/search/ produces similarly large negative activation energies.
(direct link)

[JacksonBurns]

@mjohnson541 have you seen this before?

[rwest]

This is the rate rule being used in that case.
[CsJ2-C;CsJ2(CsC);CH=C]
https://rmg.mit.edu/database/kinetics/families/Singlet_Carbene_Intra_Disproportionation/rules/3/
The rule claims to have been derived from 2003_Miller_Propargyl_Recomb_High_P reaction X <=> IX
https://rmg.mit.edu/database/kinetics/libraries/2003_Miller_Propargyl_Recomb_High_P/12/
(Miller, J. A.; Klippenstein, S. J., The Recombination of Propargyl Radicals and Other Reactions on a C6H6 Potential. J. Phys. Chem. A 2003, 107, 7783-7799.)
Putting aside for a moment whether there's any reason to expect the C6 ring to represent the long chain, it's hard to see why that rate led to the given rule.
The training reaction is at https://github.com/ReactionMechanismGenerator/RMG-database/blob/64dd3b60ae6f9d5c71ee779c1c4a564099500bf0/input/kinetics/families/Singlet_Carbene_Intra_Disproportionation/training/reactions.py#L40
The rule is being derived at load time, rather than being stored in the database.

I notice the website is using RMG-database from 13 Jun 2024, which is before that family was re-trained using the new trees and blowers masel notebook ReactionMechanismGenerator/RMG-database#494 merged August 2024.
Does the problem persist with the new autogen tree?

[mjohnson541]

It seems in the database on the website the training reactions for that family all have negative barriers, but the ones in the main database are all positive?

[rwest]

It seems in the database on the website the training reactions for that family all have negative barriers, but the ones in the main database are all positive?

I'm not seeing the training reactions as having negative barriers in the website
https://rmg.mit.edu/database/kinetics/families/Singlet_Carbene_Intra_Disproportionation/training/

The training reaction is

entry(
    index = 2,
    label = "C6H6-5 <=> C6H6-6",
    degeneracy = 1.0,
    kinetics = Arrhenius(A=(1.865e+11, 's^-1'), n=0.577, Ea=(29.169, 'kcal/mol'), T0=(1, 'K')),
    rank = 5,
    shortDesc = u"""Training reaction from kinetics library: 2003_Miller_Propargyl_Recomb_High_P""",
    longDesc = 
u"""
Taken from entry: X <=> IX
""",
)

which makes it a bit of a mystery why the rule that is apparently derived from it
https://rmg.mit.edu/database/kinetics/families/Singlet_Carbene_Intra_Disproportionation/rules/3/
has that value. Is it from trying to fit Evans Polanyi?

I'm not sure how far someone should look into this.
This instance of the problem has probably gone away because the tree was re-trained.

But does it indicate there's a bug in the old style tree methods, that someone should chase down and fix?