cfgoldsmith / RMG-Py

Python version of the amazing Reaction Mechanism Generator (RMG).
http://reactionmechanismgenerator.github.io/RMG-Py/
Other
1 stars 4 forks source link

Unable to calculate degeneracy for bidentate species #7

Closed cfgoldsmith closed 7 years ago

cfgoldsmith commented 8 years ago

If Surface_Adsorption_Bidentate is turned on, RMG-Cat crashes. It looks like it wants to desorb one of the two sigma bonds, but then it crashes. For example:

Error: Reactant: Molecule(SMILES="O=[C]O[Ni]") Error: Reactant: Molecule(SMILES="[Ni]") Error: Product: Molecule(SMILES="O=C([Ni])O[Ni]") rmgpy.data.kinetics.common.KineticsError: Unable to calculate degeneracy for reaction <Molecule "O=[C]O[Ni]"> + <Molecule "[Ni]"> <=> <Molecule "O=C([Ni])O[Ni]"> in reaction family Surface_Adsorption_Single. Expected 1 reaction but generated 0

In this case, the reaction in question is actually in the reverse direction (I think): bidentate (or di-sigma) CO2 tries to break the X-C bond to make X-O-C.=O

cfgoldsmith commented 8 years ago

I thought I could be clever and create a forbidden structure in Surface_Adsorption_Single:

forbidden(
    label = "disigma1",
    group = 
"""
1 *1 R u1 {2,[S,D]}
2    R u0 {1,[S,D]} {3,[S,D]}
3    X u0 {2,[S,D]}
""",
    shortDesc = u"""""",
    longDesc = 
u"""
The adsorbing atom should not be adjacent to an atom that is already adsorbed.
e.g. this is not allowed:

H2C-O   <-->  H2C*-O
  | |              |
  X X           X  X
""",
)

This seems to have worked, but now I get a different error: Error: Reactant: Molecule(SMILES="O=CO[Ni]") Error: Reactant: Molecule(SMILES="[Ni]") Error: Reactant: Molecule(SMILES="[Ni]") Error: Product: Molecule(SMILES="[Ni]") Error: Product: Molecule(SMILES="O=C([Ni])O[Ni]") Traceback (most recent call last): File "/Users/claudefranklingoldsmith/code/RMG-Py/rmgpy/scoop_framework/util.py", line 108, in call return self.myfn(_args, *_kwargs) File "/Users/claudefranklingoldsmith/code/RMG-Py/rmgpy/rmg/react.py", line 104, in reactMolecules rxns = family.generateReactions(molecules) File "/Users/claudefranklingoldsmith/code/RMG-Py/rmgpy/data/kinetics/family.py", line 1351, in generateReactions reactionList.extend(self.generateReactions(reactants, forward=False)) File "/Users/claudefranklingoldsmith/code/RMG-Py/rmgpy/data/kinetics/family.py", line 1725, in generateReactions reaction.degeneracy = self.calculateDegeneracy(reaction) File "/Users/claudefranklingoldsmith/code/RMG-Py/rmgpy/data/kinetics/family.py", line 1368, in calculateDegeneracy 'but generated {2}').format(reaction, self.label, len(reactions))) KineticsError: Unable to calculate degeneracy for reaction <Molecule "O=CO[Ni]"> + <Molecule "[Ni]"> + <Molecule "[Ni]"> <=> <Molecule "[Ni]"> + <Molecule "O=C([Ni])O[Ni]"> in reaction family Surface_Adsorption_Dissociative. Expected 1 reaction but generated 0

rwest commented 8 years ago

Your analysis of the first bug seems right and your first idea seems like a good one. There may be a more universal or robust way, but that should work.

The secord bug seems to be in Surface_Adsorption_Dissociative, again it looks like it makes it in the reverse direction then can't find it in the forward. The not finding it in the forward I understand (and it's correct) but why does it make the reverse? Hmm

cfgoldsmith commented 8 years ago

I'm a little confused by the SMILES in this case. I had thought that:

SMILES="O=[C]O[Ni] is equivalent to:


O=C.
  |
  O
  |
  X

And SMILES="O=C([Ni])O[Ni] is


O=C-O
  | |
  X X

So if that is the case, then is SMILES="O=CO[Ni] equal to


O=CH
  |
  O
  |
  X

? Either way, I'm having wee trouble deciphering the second reaction.

rwest commented 8 years ago

That is indeed what each of those SMILES should mean, which would imply that the reaction vanishes/creates a hydrogen atom out of nowhere. Which doesn't make sense. Probably it made the SMILES wrong. My guess, given the Surface_Adsorption_Dissociative family, is it meant this:

O=CH
  |
  O                    O=C-O      H
  |                      | |      |
  X  + X  + X     <=>    X X  +   X

Which if you forget that one of the X is a surface site would look like this:

  OX                     OX
  |                      |
O=C-H                  O=C     H
                         |     |
        + X + X  <=>     X  +  X

which is a perfectly reasonable dissociative adsorption!

Problem once again is that in the forward direction, it doesn't let the adsorbate react because it's already adsorbed (all well and good) but in the backwards direction it doesn't realize this will be a problem.

In both cases, perhaps a hard-coded "don't let bidentates adsorbates desorb according to single or dissociative desorption reactions" rule would fix both problems. (We already have "don't let adsorbates adsorb again" hard-coded)

rwest commented 8 years ago

...or (in the meantime) you could try a forbidden structure like you did for Surface_Adsorption_Single

cfgoldsmith commented 8 years ago

OK, adding a pair of forbidden structures seems to have fixed the problem.


forbidden(
    label = "disigma1",
    group = 
"""
1 *1 R u0 {2,[S,D]}
2    R u0 {1,[S,D]} {3,[S,D]}
3    X u0 {2,[S,D]}
""",
    shortDesc = u"""""",
    longDesc = 
u"""
The adsorbing atom should not be adjacent to a molecule that is already adsorbed.
e.g. this is not allowed:

H  O=C-O   <-->  O=CH-O
|    | |              |
X    X X         X X  X
""",
)

forbidden(
    label = "disigma2",
    group = 
"""
1 *2 R u0 {2,[S,D]}
2    R u0 {1,[S,D]} {3,[S,D]}
3    X u0 {2,[S,D]}
""",
    shortDesc = u"""""",
    longDesc = 
u"""
The adsorbing atom should not be adjacent to a molecule that is already adsorbed.
e.g. this is not allowed:

H  O=C-O   <-->  O=CH-O
|    | |              |
X    X X         X X  X
""",
)

Moving forward, I'm not sure what is better (hard coding vs forbidden). Hard coding seems more straight forward and cleaner. What is your experience?

cfgoldsmith commented 8 years ago

OK, I think hard coding is the way to go. Otherwise, we have to add a lot of forbidden structures, because as we increase the chain length on the surface, we have to worry about next-nearest neighbors with surface sites, etc. I think a more robust "Don't Adsorb Anything That Already Is Adsorbed!" is the way to go.

rwest commented 8 years ago

We have already forbidden

So we need to forbid:

but not forbid Surface_Adsorption_Bidentate in the reverse direction if the adsorbate is bidentate. (forward: AB+X+X <=> AXBX)

Right?

rwest commented 8 years ago

I've decided it may be simpler to just say: "Families named 'adsorption' in reverse must make a gas phase species".

Implemented what I think is a fix, but can't make it hit my breakpoints, so perhaps my input file and/or database are not like yours. Let me know if you reproduce the problems after pulling my code.

cfgoldsmith commented 8 years ago

Excellent. I will (try to) pull the patch in the morning.