Unable to calculate degeneracy for reactions involving "[CH]=C[N]" in reaction family Disproportionation

[alongd]

Apparently the Disproportionation family doesn't deal well with reactions of the type:

R.CH2R \ :CH2 + :N-CH-.CH <=> RCHR \ .CH3 + .N=C=.CH

I got errors like the following when running ethylamine pyrolysis simulations when specific reactions were added to the seed:

Unable to calculate degeneracy for reaction <Molecule "[CH2][CH]N"> + <Molecule "[CH]=C[N]"> <=> <Molecule "C[CH]N"> + <Molecule "C(=[CH])=[N]"> in reaction family Disproportionation. Expected 1 reaction but generated 2

[alongd]

A resonance structure of the C2NH species raised this degeneracy error:
C#C[N] <=>[CH]=C=[N]
So two reactions were generated

[nickvandewiele]

could you provide the exact error trace that show which methods were hit?

[mliu49]

I actually saw a similar case recently, where RMG was unable to calculate degeneracy due to resonance structures.

Traceback (most recent call last):
  File "/home/mjliu/Code/RMG-Py/rmg.py", line 165, in <module>
    rmg.execute(inputFile, output_dir, **kwargs)
  File "/home/mjliu/Code/RMG-Py/rmgpy/rmg/main.py", line 684, in execute
    forcedRecombinationProducts = self.reactionModel.enlargeEdge(unimolecularReact, bimolecularReact)
  File "/home/mjliu/Code/RMG-Py/rmgpy/rmg/model.py", line 676, in enlargeEdge
    self.processNewReactions(self.react(database, self.core.species[i]), self.core.species[i], None)
  File "/home/mjliu/Code/RMG-Py/rmgpy/rmg/model.py", line 649, in react
    reactionList.extend(database.kinetics.generateReactionsFromFamilies([moleculeA], products=None, failsSpeciesConstraints=self.failsSpeciesConstraints, only_families=only_families))
  File "/home/mjliu/Code/RMG-Py/rmgpy/data/kinetics/database.py", line 419, in generateReactionsFromFamilies
    reactionList.extend(family.generateReactions(reactants, failsSpeciesConstraints=failsSpeciesConstraints))
  File "/home/mjliu/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 1246, in generateReactions
    reactionList.extend(self.__generateReactions(reactants, forward=False, failsSpeciesConstraints=failsSpeciesConstraints))
  File "/home/mjliu/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 1478, in __generateReactions
    reaction.degeneracy = self.calculateDegeneracy(reaction)
  File "/home/mjliu/Code/RMG-Py/rmgpy/data/kinetics/family.py", line 1274, in calculateDegeneracy
    'but generated {2}').format(reaction, self.label, len(reactions)))
rmgpy.data.kinetics.common.KineticsError: Unable to calculate degeneracy for reaction <Molecule "CC12C=CC3(CCC3(C)C1)C(S)=C2S"> <=> <Molecule "C=C(C)CCC1=CC=C(C)C(S)=C1S"> in reaction family Intra_Diels_alder. Expected 1 reaction but generated 2

[alongd]

Yes, looks similar. My trace is:

Traceback (most recent call last):
  File "/home/alongd/ws/RMG-Py/rmg.py", line 165, in <module>
    rmg.execute(inputFile, output_dir, **kwargs)
  File "/home/alongd/ws/RMG-Py/rmgpy/rmg/main.py", line 566, in execute
    self.reactionModel.enlarge(objectToEnlarge)
  File "/home/alongd/ws/RMG-Py/rmgpy/rmg/model.py", line 713, in enlarge
    [corespeciesList] * familieCount
  File "/home/alongd/ws/RMG-Py/rmgpy/scoop_framework/util.py", line 153, in map_
    return map(*args, **kwargs)
  File "/home/alongd/ws/RMG-Py/rmgpy/scoop_framework/util.py", line 108, in __call__
    return self.myfn(*args, **kwargs)
  File "/home/alongd/ws/RMG-Py/rmgpy/rmg/model.py", line 1834, in reactFamily
    [familyKey] * len(combos),
  File "/home/alongd/ws/RMG-Py/rmgpy/scoop_framework/util.py", line 153, in map_
    return map(*args, **kwargs)
  File "/home/alongd/ws/RMG-Py/rmgpy/scoop_framework/util.py", line 108, in __call__
    return self.myfn(*args, **kwargs)
  File "/home/alongd/ws/RMG-Py/rmgpy/rmg/model.py", line 1855, in reactSpecies
    [familyKey] * len(combos),
  File "/home/alongd/ws/RMG-Py/rmgpy/scoop_framework/util.py", line 153, in map_
    return map(*args, **kwargs)
  File "/home/alongd/ws/RMG-Py/rmgpy/scoop_framework/util.py", line 108, in __call__
    return self.myfn(*args, **kwargs)
  File "/home/alongd/ws/RMG-Py/rmgpy/rmg/model.py", line 1872, in reactMolecules
    reactionList = family.generateReactions(list(molecules))
  File "/home/alongd/ws/RMG-Py/rmgpy/data/kinetics/family.py", line 1228, in generateReactions
    reactionList.extend(self.__generateReactions(reactants, forward=False))
  File "/home/alongd/ws/RMG-Py/rmgpy/data/kinetics/family.py", line 1458, in __generateReactions
    reaction.degeneracy = self.calculateDegeneracy(reaction)
  File "/home/alongd/ws/RMG-Py/rmgpy/data/kinetics/family.py", line 1256, in calculateDegeneracy
    'but generated {2}').format(reaction, self.label, len(reactions)))
rmgpy.data.kinetics.common.KineticsError: Unable to calculate degeneracy for reaction <Molecule "[CH2][CH]N"> + <Molecule "[CH]=C[N]"> <=> <Molecule "C[CH]N"> + <Molecule "C(=[CH])=[N]"> in reaction family Disproportionation. Expected 1 reaction but generated 2

The KineticsError rises, probably due to this resonance. The two considered reactions are:

<Molecule "N[CH][CH]"> + <Molecule "C(=[CH])[N]"> <=> <Molecule "N[CH][CH2]"> + <Molecule "C#C[N]">
<Molecule "N[CH][CH]"> + <Molecule "C(=[CH])[N]"> <=> <Molecule "N[CH][CH2]"> + <Molecule "[CH]=C=[N]">

[mliu49]

Looking at calculateDegeneracy, it's not apparent to me why the degeneracy is always expected to be 1.

[connie]

The degeneracy is not expected to be one. However the function expects only 1 reaction to come from the calculation. It expect RMG to find the multiple paths, and combine them into 1 reaction having a degeneracy = numPaths. See https://github.com/ReactionMechanismGenerator/RMG-Py/blob/master/rmgpy/data/kinetics/family.py#L1444

Of course this is a problem when you have multiple transition states, and shouldn't actually be combining the reactions into one. And also when there are resonance forms like what you guys are seeing, where RMG can't combine it into 1 reaction.

Related to:
#142
#494
#323

[connie]

Looks like the combining of reactions happens here: https://github.com/ReactionMechanismGenerator/RMG-Py/blob/master/rmgpy/data/kinetics/family.py#L1410

So it seems that we do try to combine species that are resonance isomers of each other at this point, so this should not typically cause an error. However, upon closer study I found that the reason is becausegenerateResonanceIsomers() for C#C[N] is problematic:

from rmgpy.molecule import Molecule
m = Molecule(SMILES='[CH]=C=[N]')
m.generateResonanceIsomers()
[Molecule(SMILES="[CH]=C=[N]"), Molecule(SMILES="C#C[N]"), Molecule(SMILES="[CH]C#N")]
m2 = Molecule(SMILES='C#C[N]')
m2.generateResonanceIsomers()
[Molecule(SMILES="C#C[N]")]

As you can see, the resonance isomers for [CH]=C=[N] includes C#C[N], but generating resonance isomers for C#C[N] does not produce anything. Which means our generateResonanceIsomers() function has hysteresis.

Therefore in the above case where the reaction that produces C#C[N] first, the 2nd reaction that produces [CH]=C=[N] is not seen as identical, even though had we reversed the order they would be considered identical.

[mliu49]

I think that explains my case as well. The product species in the reaction I posted has a benzene ring, so there are two possible kekule structures.

>>> m1 = Molecule(SMILES='C=C(C)CCC1=CC=C(C)C(S)=C1S')
>>> m1.generateResonanceIsomers()
[Molecule(SMILES="C=C(C)CCC1=CC=C(C)C(S)=C1S"), Molecule(SMILES="C=C(C)CCc1ccc(C)c(S)c1S"), Molecule(SMILES="C=C(C)CCC1C=CC(C)=C(S)C=1S")]
>>> m2 = Molecule(SMILES='C=C(C)CCC1C=CC(C)=C(S)C=1S')
>>> m2.generateResonanceIsomers()
[Molecule(SMILES="C=C(C)CCC1C=CC(C)=C(S)C=1S"), Molecule(SMILES="C=C(C)CCc1ccc(C)c(S)c1S")]

So starting with one of the resonance structures gives both kekule forms, while starting with the other one somehow doesn't.

[goldmanm]

overhauled in #1055