molecule branch: list index out of range in Molecule.isLinear

[rwest]

Running the methylformate example on @jwallen's molecule branch (jwallen/RMG-Py@d177753) it crashes pretty quickly with this:

Exploring isomer CH2CH2OH(44) in pressure-dependent network #665
Generating thermodynamics for new species...
Warning: Average RMS error in heat capacity fit to CO[CH]O[O](73) = 0.178511*R
Traceback (most recent call last):
  File "rmg.py", line 134, in <module>
    cProfile.runctx(command, global_vars, local_vars, stats_file)
  File "/usr/local/Cellar/python/2.7.3/Frameworks/Python.framework/Versions/2.7/lib/python2.7/cProfile.py", line 49, in runctx
    prof = prof.runctx(statement, globals, locals)
  File "/usr/local/Cellar/python/2.7.3/Frameworks/Python.framework/Versions/2.7/lib/python2.7/cProfile.py", line 140, in runctx
    exec cmd in globals, locals
  File "<string>", line 1, in <module>
  File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/main.py", line 315, in execute
    self.initialize(args)
  File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/main.py", line 303, in initialize
    self.reactionModel.enlarge([spec for spec in self.initialSpecies if spec.reactive])
  File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/model.py", line 580, in enlarge
    spec.generateThermoData(database)
  File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/model.py", line 87, in generateThermoData
    linear = self.molecule[0].isLinear()
  File "molecule.py", line 1015, in rmgpy.molecule.molecule.Molecule.isLinear (build/pyrex/rmgpy/molecule/molecule.c:16892)
  File "molecule.py", line 1048, in rmgpy.molecule.molecule.Molecule.isLinear (build/pyrex/rmgpy/molecule/molecule.c:16722)
IndexError: list index out of range

or, running after a make decython:

Generating thermodynamics for new species...
Warning: Average RMS error in heat capacity fit to CO[CH]O[O](73) = 0.178511*R
Traceback (most recent call last):
  File "rmg.py", line 142, in <module>
    rmg.execute(args)
  File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/main.py", line 315, in execute
    self.initialize(args)
  File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/main.py", line 303, in initialize
    self.reactionModel.enlarge([spec for spec in self.initialSpecies if spec.reactive])
  File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/model.py", line 580, in enlarge
    spec.generateThermoData(database)
  File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/model.py", line 87, in generateThermoData
    linear = self.molecule[0].isLinear()
  File "/Users/rwest/Code/rmgpy/RMG-Py/rmgpy/molecule/molecule.py", line 1048, in isLinear
    if bonds[0].isSingle() and bonds[1].isTriple():
IndexError: list index out of range

[rwest]

I ran it in ipython with --pdb and had a look around:

ipdb> bonds
[]
ipdb> atom
<Atom 'C'>
ipdb> self
Molecule(SMILES="[C].C[O].[O].[O].[H]")
ipdb> u
> /Users/rwest/Code/rmgpy/RMG-Py/rmgpy/rmg/model.py(87)generateThermoData()
     86             else:
---> 87                 linear = self.molecule[0].isLinear()
     88                 nRotors = self.molecule[0].countInternalRotors()

ipdb> self.label
'C1(OC)OO1'

that means that it was meant to be C1(OC)OO1 but by the time it got to line 1048 of isLinear(), it had no bonds! (apart from one C-O bond)

But when I play around with that molecule in isolation, both isCyclic() and isLinear() work fine.

ipdb> aaa=Molecule(SMILES=('C1(OC)OO1'))
ipdb> aaa
Molecule(SMILES="C1(OC)OO1")
ipdb> aaa.isCyclic()
True
ipdb> aaa
Molecule(SMILES="C1(OC)OO1")
ipdb> aaa.isLinear()
False
ipdb> aaa
Molecule(SMILES="C1(OC)OO1")

I have put some extra lines in that print the SMILES string at various points. The molecule is having its bonds destroyed during the call to database.thermo.getThermoData(self) on line 76 of rmg/model.py.

The first handful of (noncyclic) species get through fine, then somewhere between calling generateThermoData (line 68) and isLinear (line 87) the SMILES changes to one with missing bonds:

Starting now...
Attempting COC[O] = Molecule(SMILES="COC[O]")
Testing linearity of Molecule(SMILES="COC[O]")
Attempting CO[CH]O = Molecule(SMILES="CO[CH]O")
Testing linearity of Molecule(SMILES="CO[CH]O")
Attempting COC([O])[O] = Molecule(SMILES="COC([O])[O]")
Testing linearity of Molecule(SMILES="COC([O])[O]")
Attempting CO[CH]O[O] = Molecule(SMILES="CO[CH]O[O]")
Testing linearity of Molecule(SMILES="CO[CH]O[O]")
Warning: Average RMS error in heat capacity fit to CO[CH]O[O](73) = 0.178511*R
Attempting C1(OC)OO1 = Molecule(SMILES="C1(OC)OO1")
Testing linearity of Molecule(SMILES="[C].C[O].[O].[O].[H]")

[rwest]

The trouble is occurring in estimateThermoViaGroupAdditivity (data/thermo.py line 624)

ipdb> mmm=Molecule(SMILES="C1(OC)OO1")
ipdb> database.thermo.estimateThermoViaGroupAdditivity(mmm)
ThermoData(Tdata=([300,400,500,600,800,1000,1500],"K"), Cpdata=([91.7551,109.453,123.302,135.227,156.398,171.628,183.385],"K"), H298=(-273341,"K"), S298=(184.514,"K"), comment="""group(Cs-OsOsOsH) + gauche(Cs(RRRR)) + other(R) + group(Cs-OsHHH) + gauche(Cs(RRRR)) + other(R) + group(Os-CsCs) + gauche(Os(CsCs)) + other(R) + group(Os-OsCs) + gauche(Os(CsR)) + other(R) + group(Os-OsCs) + gauche(Os(CsR)) + other(R) + ring(Ring)""")
ipdb> mmm
Molecule(SMILES="[C].C[O].[O].[O].[H]")

Recall, it only breaks for cyclic things, I expect the trouble is around line 738 of data/thermo.py

[jwallen]

Yes, that looks like the trouble point. This should be closer (not yet tested):

                # Make a temporary structure containing only the atoms in the ring
                ringStructure = Molecule()
                newAtoms = [atom.copy() for atom in ring]
                for atom in newAtoms: ringStructure.addAtom(atom)
                for atom1 in ring:
                    for atom2 in ring:
                        if molecule.hasBond(atom1, atom2):
                            ringStructure.addBond(newAtoms[molecule.atoms.index(atom1)], newAtoms[molecule.atoms.index(atom2)], molecule.getBond(atom1, atom2).copy())

Basically the idea is that you now need to use a copy of the atoms and bonds; since they store the connectivity of the graph, they can only belong to one graph at a time.

[rwest]

No dice! The index of an atom in the original molecule may be way higher than the number of atoms in the ring, and the way addBond works has changed (you now just pass it a bond, which has already had its vertices set). I think this is what's needed (see rwest/RMG-Py@acccd69):

                # Make a temporary structure containing only the atoms in the ring
                ringStructure = Molecule()
                newAtoms = dict()
                for atom in ring:
                    newAtoms[atom] = atom.copy()
                    ringStructure.addAtom(newAtoms[atom]) # (addAtom deletes the atom's bonds)
                for atom1 in ring:
                    for atom2 in ring:
                        if molecule.hasBond(atom1, atom2):
                            ringStructure.addBond(Bond(newAtoms[atom1], newAtoms[atom2], atom1.bonds[atom2].order ))

Some functioning unit tests in the data module would help! :)