From 8c3bfba0ffd6e7b325dd9a0a36cc4632ecb57071 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Thu, 28 Feb 2013 14:47:06 -0500 Subject: [PATCH 01/64] Added atom types for nitrogen. --- rmgpy/molecule/atomtype.py | 183 +++++++++++++++++++++++-------------- 1 file changed, 116 insertions(+), 67 deletions(-) diff --git a/rmgpy/molecule/atomtype.py b/rmgpy/molecule/atomtype.py index 5d19e9c7b9..054ecfe26a 100644 --- a/rmgpy/molecule/atomtype.py +++ b/rmgpy/molecule/atomtype.py @@ -150,78 +150,106 @@ def isSpecificCaseOf(self, other): atomTypes = {} atomTypes['R'] = AtomType(label='R', generic=[], specific=[ 'R!H', - 'C','Cs','Cd','Cdd','Ct','CO','Cb','Cbf', 'CS', 'H', + 'C','Cs','Cd','Cdd','Ct','CO','Cb','Cbf', 'CS', + 'N','N3s','N3d','N3t','N3b','N4s','N4d','N4dd','N4t','N4b', 'O','Os','Od','Oa', 'Si','Sis','Sid','Sidd','Sit','SiO','Sib','Sibf', 'S','Ss','Sd','Sa'] ) atomTypes['R!H'] = AtomType(label='R!H', generic=['R'], specific=[ 'C','Cs','Cd','Cdd','Ct','CO','Cb','Cbf', 'CS', + 'N','N3s','N3d','N3t','N3b','N4s','N4d','N4dd','N4t','N4b', 'O','Os','Od','Oa', 'Si','Sis','Sid','Sidd','Sit','SiO','Sib','Sibf', 'S','Ss','Sd','Sa'] ) -atomTypes['C'] = AtomType('C', generic=['R','R!H'], specific=['Cs','Cd','Cdd','Ct','CO','Cb','Cbf','CS']) -atomTypes['Cs'] = AtomType('Cs', generic=['R','R!H', 'C'], specific=[]) -atomTypes['Cd'] = AtomType('Cd', generic=['R','R!H', 'C'], specific=[]) -atomTypes['Cdd'] = AtomType('Cdd', generic=['R','R!H', 'C'], specific=[]) -atomTypes['Ct'] = AtomType('Ct', generic=['R','R!H', 'C'], specific=[]) -atomTypes['CO'] = AtomType('CO', generic=['R','R!H', 'C'], specific=[]) -atomTypes['Cb'] = AtomType('Cb', generic=['R','R!H', 'C'], specific=[]) -atomTypes['Cbf'] = AtomType('Cbf', generic=['R','R!H', 'C'], specific=[]) -atomTypes['CS'] = AtomType('CS', generic=['R','R!H', 'C'], specific=[]) -atomTypes['H'] = AtomType('H', generic=['R'], specific=[]) -atomTypes['O'] = AtomType('O', generic=['R','R!H'], specific=['Os','Od','Oa']) -atomTypes['Os'] = AtomType('Os', generic=['R','R!H','O'], specific=[]) -atomTypes['Od'] = AtomType('Od', generic=['R','R!H','O'], specific=[]) -atomTypes['Oa'] = AtomType('Oa', generic=['R','R!H','O'], specific=[]) -atomTypes['Si'] = AtomType('Si', generic=['R','R!H'], specific=['Sis','Sid','Sidd','Sit','SiO','Sib','Sibf']) -atomTypes['Sis'] = AtomType('Sis', generic=['R','R!H','Si'], specific=[]) -atomTypes['Sid'] = AtomType('Sid', generic=['R','R!H','Si'], specific=[]) +atomTypes['H' ] = AtomType('H', generic=['R'], specific=[]) + +atomTypes['C' ] = AtomType('C', generic=['R','R!H'], specific=['Cs','Cd','Cdd','Ct','CO','Cb','Cbf','CS']) +atomTypes['Cs' ] = AtomType('Cs', generic=['R','R!H','C'], specific=[]) +atomTypes['Cd' ] = AtomType('Cd', generic=['R','R!H','C'], specific=[]) +atomTypes['Cdd' ] = AtomType('Cdd', generic=['R','R!H','C'], specific=[]) +atomTypes['Ct' ] = AtomType('Ct', generic=['R','R!H','C'], specific=[]) +atomTypes['CO' ] = AtomType('CO', generic=['R','R!H','C'], specific=[]) +atomTypes['Cb' ] = AtomType('Cb', generic=['R','R!H','C'], specific=[]) +atomTypes['Cbf' ] = AtomType('Cbf', generic=['R','R!H','C'], specific=[]) +atomTypes['CS' ] = AtomType('CS', generic=['R','R!H','C'], specific=[]) + +atomTypes['N' ] = AtomType('N', generic=['R','R!H'], specific=['N3s','N3d','N3t','N3b','N4s','N4d','N4dd','N4t','N4b']) +atomTypes['N3s' ] = AtomType('N3s', generic=['R','R!H','N'], specific=[]) +atomTypes['N3d' ] = AtomType('N3d', generic=['R','R!H','N'], specific=[]) +atomTypes['N3t' ] = AtomType('N3t', generic=['R','R!H','N'], specific=[]) +atomTypes['N3b' ] = AtomType('N3b', generic=['R','R!H','N'], specific=[]) +atomTypes['N4s' ] = AtomType('N4s', generic=['R','R!H','N'], specific=[]) +atomTypes['N4d' ] = AtomType('N4d', generic=['R','R!H','N'], specific=[]) +atomTypes['N4dd'] = AtomType('N4dd', generic=['R','R!H','N'], specific=[]) +atomTypes['N4t' ] = AtomType('N4t', generic=['R','R!H','N'], specific=[]) +atomTypes['N4b' ] = AtomType('N4b', generic=['R','R!H','N'], specific=[]) + +atomTypes['O' ] = AtomType('O', generic=['R','R!H'], specific=['Os','Od','Oa']) +atomTypes['Os' ] = AtomType('Os', generic=['R','R!H','O'], specific=[]) +atomTypes['Od' ] = AtomType('Od', generic=['R','R!H','O'], specific=[]) +atomTypes['Oa' ] = AtomType('Oa', generic=['R','R!H','O'], specific=[]) + +atomTypes['Si' ] = AtomType('Si', generic=['R','R!H'], specific=['Sis','Sid','Sidd','Sit','SiO','Sib','Sibf']) +atomTypes['Sis' ] = AtomType('Sis', generic=['R','R!H','Si'], specific=[]) +atomTypes['Sid' ] = AtomType('Sid', generic=['R','R!H','Si'], specific=[]) atomTypes['Sidd'] = AtomType('Sidd', generic=['R','R!H','Si'], specific=[]) -atomTypes['Sit'] = AtomType('Sit', generic=['R','R!H','Si'], specific=[]) -atomTypes['SiO'] = AtomType('SiO', generic=['R','R!H','Si'], specific=[]) -atomTypes['Sib'] = AtomType('Sib', generic=['R','R!H','Si'], specific=[]) +atomTypes['Sit' ] = AtomType('Sit', generic=['R','R!H','Si'], specific=[]) +atomTypes['SiO' ] = AtomType('SiO', generic=['R','R!H','Si'], specific=[]) +atomTypes['Sib' ] = AtomType('Sib', generic=['R','R!H','Si'], specific=[]) atomTypes['Sibf'] = AtomType('Sibf', generic=['R','R!H','Si'], specific=[]) -atomTypes['S'] = AtomType('S', generic=['R','R!H'], specific=['Ss','Sd','Sa']) -atomTypes['Ss'] = AtomType('Ss', generic=['R','R!H','S'], specific=[]) -atomTypes['Sd'] = AtomType('Sd', generic=['R','R!H','S'], specific=[]) -atomTypes['Sa'] = AtomType('Sa', generic=['R','R!H','S'], specific=[]) - -atomTypes['R' ].setActions(incrementBond=['R'], decrementBond=['R'], formBond=['R'], breakBond=['R'], incrementRadical=['R'], decrementRadical=['R']) -atomTypes['R!H' ].setActions(incrementBond=['R!H'], decrementBond=['R!H'], formBond=['R!H'], breakBond=['R!H'], incrementRadical=['R!H'], decrementRadical=['R!H']) - -atomTypes['C' ].setActions(incrementBond=['C'], decrementBond=['C'], formBond=['C'], breakBond=['C'], incrementRadical=['C'], decrementRadical=['C']) -atomTypes['Cs' ].setActions(incrementBond=['Cd','CO','Cs'], decrementBond=[], formBond=['Cs'], breakBond=['Cs'], incrementRadical=['Cs'], decrementRadical=['Cs']) -atomTypes['Cd' ].setActions(incrementBond=['Cdd','Ct'], decrementBond=['Cs'], formBond=['Cd'], breakBond=['Cd'], incrementRadical=['Cd'], decrementRadical=['Cd']) -atomTypes['Cdd' ].setActions(incrementBond=[], decrementBond=['Cd','CO','CS'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) -atomTypes['Ct' ].setActions(incrementBond=[], decrementBond=['Cd'], formBond=['Ct'], breakBond=['Ct'], incrementRadical=['Ct'], decrementRadical=['Ct']) -atomTypes['CO' ].setActions(incrementBond=['Cdd'], decrementBond=['Cs'], formBond=['CO'], breakBond=['CO'], incrementRadical=['CO'], decrementRadical=['CO']) -atomTypes['CS' ].setActions(incrementBond=['Cdd'], decrementBond=['Cs'], formBond=['CS'], breakBond=['CS'], incrementRadical=['CS'], decrementRadical=['CS']) -atomTypes['Cb' ].setActions(incrementBond=[], decrementBond=[], formBond=['Cb'], breakBond=['Cb'], incrementRadical=['Cb'], decrementRadical=['Cb']) -atomTypes['Cbf' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) - -atomTypes['H' ].setActions(incrementBond=[], decrementBond=[], formBond=['H'], breakBond=['H'], incrementRadical=['H'], decrementRadical=['H']) - -atomTypes['O' ].setActions(incrementBond=['O'], decrementBond=['O'], formBond=['O'], breakBond=['O'], incrementRadical=['O'], decrementRadical=['O']) -atomTypes['Os' ].setActions(incrementBond=['Od'], decrementBond=[], formBond=['Os'], breakBond=['Os'], incrementRadical=['Os'], decrementRadical=['Os']) -atomTypes['Od' ].setActions(incrementBond=[], decrementBond=['Os'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) -atomTypes['Oa' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) - -atomTypes['Si' ].setActions(incrementBond=['Si'], decrementBond=['Si'], formBond=['Si'], breakBond=['Si'], incrementRadical=['Si'], decrementRadical=['Si']) -atomTypes['Sis' ].setActions(incrementBond=['Sid','SiO'], decrementBond=[], formBond=['Sis'], breakBond=['Sis'], incrementRadical=['Sis'], decrementRadical=['Sis']) -atomTypes['Sid' ].setActions(incrementBond=['Sidd','Sit'], decrementBond=['Sis'], formBond=['Sid'], breakBond=['Sid'], incrementRadical=['Sid'], decrementRadical=['Sid']) -atomTypes['Sidd' ].setActions(incrementBond=[], decrementBond=['Sid','SiO'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) -atomTypes['Sit' ].setActions(incrementBond=[], decrementBond=['Sid'], formBond=['Sit'], breakBond=['Sit'], incrementRadical=['Sit'], decrementRadical=['Sit']) -atomTypes['SiO' ].setActions(incrementBond=['Sidd'], decrementBond=['Sis'], formBond=['SiO'], breakBond=['SiO'], incrementRadical=['SiO'], decrementRadical=['SiO']) -atomTypes['Sib' ].setActions(incrementBond=[], decrementBond=[], formBond=['Sib'], breakBond=['Sib'], incrementRadical=['Sib'], decrementRadical=['Sib']) -atomTypes['Sibf' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) - -atomTypes['S' ].setActions(incrementBond=['S'], decrementBond=['S'], formBond=['S'], breakBond=['S'], incrementRadical=['S'], decrementRadical=['S']) -atomTypes['Ss' ].setActions(incrementBond=['Sd'], decrementBond=[], formBond=['Ss'], breakBond=['Ss'], incrementRadical=['Ss'], decrementRadical=['Ss']) -atomTypes['Sd' ].setActions(incrementBond=[], decrementBond=['Ss'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) -atomTypes['Sa' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) + +atomTypes['S' ] = AtomType('S', generic=['R','R!H'], specific=['Ss','Sd','Sa']) +atomTypes['Ss' ] = AtomType('Ss', generic=['R','R!H','S'], specific=[]) +atomTypes['Sd' ] = AtomType('Sd', generic=['R','R!H','S'], specific=[]) +atomTypes['Sa' ] = AtomType('Sa', generic=['R','R!H','S'], specific=[]) + +atomTypes['R' ].setActions(incrementBond=['R'], decrementBond=['R'], formBond=['R'], breakBond=['R'], incrementRadical=['R'], decrementRadical=['R']) +atomTypes['R!H' ].setActions(incrementBond=['R!H'], decrementBond=['R!H'], formBond=['R!H'], breakBond=['R!H'], incrementRadical=['R!H'], decrementRadical=['R!H']) + +atomTypes['H' ].setActions(incrementBond=[], decrementBond=[], formBond=['H'], breakBond=['H'], incrementRadical=['H'], decrementRadical=['H']) + +atomTypes['C' ].setActions(incrementBond=['C'], decrementBond=['C'], formBond=['C'], breakBond=['C'], incrementRadical=['C'], decrementRadical=['C']) +atomTypes['Cs' ].setActions(incrementBond=['Cd','CO','Cs'], decrementBond=[], formBond=['Cs'], breakBond=['Cs'], incrementRadical=['Cs'], decrementRadical=['Cs']) +atomTypes['Cd' ].setActions(incrementBond=['Cdd','Ct'], decrementBond=['Cs'], formBond=['Cd'], breakBond=['Cd'], incrementRadical=['Cd'], decrementRadical=['Cd']) +atomTypes['Cdd' ].setActions(incrementBond=[], decrementBond=['Cd','CO','CS'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) +atomTypes['Ct' ].setActions(incrementBond=[], decrementBond=['Cd'], formBond=['Ct'], breakBond=['Ct'], incrementRadical=['Ct'], decrementRadical=['Ct']) +atomTypes['CO' ].setActions(incrementBond=['Cdd'], decrementBond=['Cs'], formBond=['CO'], breakBond=['CO'], incrementRadical=['CO'], decrementRadical=['CO']) +atomTypes['CS' ].setActions(incrementBond=['Cdd'], decrementBond=['Cs'], formBond=['CS'], breakBond=['CS'], incrementRadical=['CS'], decrementRadical=['CS']) +atomTypes['Cb' ].setActions(incrementBond=[], decrementBond=[], formBond=['Cb'], breakBond=['Cb'], incrementRadical=['Cb'], decrementRadical=['Cb']) +atomTypes['Cbf' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) + +atomTypes['N' ].setActions(incrementBond=['N'], decrementBond=['N'], formBond=['N'], breakBond=['N'], incrementRadical=['N'], decrementRadical=['N']) +atomTypes['N3s' ].setActions(incrementBond=['N3d','N3s'], decrementBond=[], formBond=['N3s','N4s'], breakBond=['N3s'], incrementRadical=['N3s'], decrementRadical=['N3s']) +atomTypes['N3d' ].setActions(incrementBond=['N3t'], decrementBond=['N3s'], formBond=['N3d','N4d'], breakBond=['N3d'], incrementRadical=['N3d'], decrementRadical=['N3d']) +atomTypes['N3t' ].setActions(incrementBond=[], decrementBond=['N3d'], formBond=['N4t'], breakBond=[], incrementRadical=['N3t'], decrementRadical=['N3t']) +atomTypes['N3b' ].setActions(incrementBond=[], decrementBond=[], formBond=['N3b'], breakBond=['N3b'], incrementRadical=['N3b'], decrementRadical=['N3b']) +atomTypes['N4s' ].setActions(incrementBond=['N4d'], decrementBond=[], formBond=[], breakBond=['N3s'], incrementRadical=['N4s'], decrementRadical=['N4s']) +atomTypes['N4d' ].setActions(incrementBond=['N4dd','N4t'], decrementBond=['N4s'], formBond=[], breakBond=['N3d'], incrementRadical=['N4d'], decrementRadical=['N4d']) +atomTypes['N4dd'].setActions(incrementBond=[], decrementBond=['N3d'], formBond=[], breakBond=[], incrementRadical=['N4dd'], decrementRadical=['N4dd']) +atomTypes['N4t' ].setActions(incrementBond=[], decrementBond=['N3d','N3t'], formBond=[], breakBond=['N3d','N3t'], incrementRadical=['N4t'], decrementRadical=['N4t']) +atomTypes['N4b' ].setActions(incrementBond=[], decrementBond=[], formBond=['N4b'], breakBond=['N4b'], incrementRadical=['N4b'], decrementRadical=['N4b']) + +atomTypes['O' ].setActions(incrementBond=['O'], decrementBond=['O'], formBond=['O'], breakBond=['O'], incrementRadical=['O'], decrementRadical=['O']) +atomTypes['Os' ].setActions(incrementBond=['Od'], decrementBond=[], formBond=['Os'], breakBond=['Os'], incrementRadical=['Os'], decrementRadical=['Os']) +atomTypes['Od' ].setActions(incrementBond=[], decrementBond=['Os'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) +atomTypes['Oa' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) + +atomTypes['Si' ].setActions(incrementBond=['Si'], decrementBond=['Si'], formBond=['Si'], breakBond=['Si'], incrementRadical=['Si'], decrementRadical=['Si']) +atomTypes['Sis' ].setActions(incrementBond=['Sid','SiO'], decrementBond=[], formBond=['Sis'], breakBond=['Sis'], incrementRadical=['Sis'], decrementRadical=['Sis']) +atomTypes['Sid' ].setActions(incrementBond=['Sidd','Sit'], decrementBond=['Sis'], formBond=['Sid'], breakBond=['Sid'], incrementRadical=['Sid'], decrementRadical=['Sid']) +atomTypes['Sidd'].setActions(incrementBond=[], decrementBond=['Sid','SiO'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) +atomTypes['Sit' ].setActions(incrementBond=[], decrementBond=['Sid'], formBond=['Sit'], breakBond=['Sit'], incrementRadical=['Sit'], decrementRadical=['Sit']) +atomTypes['SiO' ].setActions(incrementBond=['Sidd'], decrementBond=['Sis'], formBond=['SiO'], breakBond=['SiO'], incrementRadical=['SiO'], decrementRadical=['SiO']) +atomTypes['Sib' ].setActions(incrementBond=[], decrementBond=[], formBond=['Sib'], breakBond=['Sib'], incrementRadical=['Sib'], decrementRadical=['Sib']) +atomTypes['Sibf'].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) + +atomTypes['S' ].setActions(incrementBond=['S'], decrementBond=['S'], formBond=['S'], breakBond=['S'], incrementRadical=['S'], decrementRadical=['S']) +atomTypes['Ss' ].setActions(incrementBond=['Sd'], decrementBond=[], formBond=['Ss'], breakBond=['Ss'], incrementRadical=['Ss'], decrementRadical=['Ss']) +atomTypes['Sd' ].setActions(incrementBond=[], decrementBond=['Ss'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) +atomTypes['Sa' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) for atomType in atomTypes.values(): for items in [atomType.generic, atomType.specific, @@ -235,21 +263,23 @@ def getAtomType(atom, bonds): Determine the appropriate atom type for an :class:`Atom` object `atom` with local bond structure `bonds`, a ``dict`` containing atom-bond pairs. """ - + cython.declare(atomType=str) - cython.declare(double=cython.double, double0=cython.double, triple=cython.double, benzene=cython.double) + cython.declare(single=cython.double, double=cython.double, double0=cython.double, triple=cython.double, benzene=cython.double) atomType = '' # Count numbers of each higher-order bond type - double = 0; doubleO = 0; triple = 0; benzene = 0 + single = 0; double = 0; doubleO = 0; triple = 0; benzene = 0 for atom2, bond12 in bonds.iteritems(): - if bond12.isDouble(): - if atom2.isOxygen(): doubleO +=1 + if bond12.isSingle(): + single += 1 + elif bond12.isDouble(): + if atom2.isOxygen(): doubleO += 1 else: double += 1 elif bond12.isTriple(): triple += 1 elif bond12.isBenzene(): benzene += 1 - + # Use element and counts to determine proper atom type if atom.symbol == 'C': if double == 0 and doubleO == 0 and triple == 0 and benzene == 0: atomType = 'Cs' @@ -261,6 +291,25 @@ def getAtomType(atom, bonds): elif double == 0 and doubleO == 0 and triple == 0 and benzene == 3: atomType = 'Cbf' elif atom.symbol == 'H': atomType = 'H' + elif atom.symbol == 'N': + if double == 0 and doubleO == 0 and triple == 0 and benzene == 0 and single == 0: atomType = 'N3s' + elif double == 0 and doubleO == 0 and triple == 0 and benzene == 0 and single == 1: atomType = 'N3s' + elif double == 0 and doubleO == 0 and triple == 0 and benzene == 0 and single == 2: atomType = 'N3s' + elif double == 0 and doubleO == 0 and triple == 0 and benzene == 0 and single == 3: atomType = 'N3s' + elif double == 1 and doubleO == 0 and triple == 0 and benzene == 0 and single == 0: atomType = 'N3d' + elif double == 0 and doubleO == 1 and triple == 0 and benzene == 0 and single == 0: atomType = 'N3d' + elif double == 1 and doubleO == 0 and triple == 0 and benzene == 0 and single == 1: atomType = 'N3d' + elif double == 0 and doubleO == 1 and triple == 0 and benzene == 0 and single == 1: atomType = 'N3d' + elif double == 0 and doubleO == 0 and triple == 1 and benzene == 0 and single == 0: atomType = 'N3t' + elif double == 0 and doubleO == 0 and triple == 0 and benzene == 2 and single == 0: atomType = 'N3b' + elif double == 0 and doubleO == 0 and triple == 0 and benzene == 0 and single == 4: atomType = 'N4s' + elif double == 1 and doubleO == 0 and triple == 0 and benzene == 0 and single == 2: atomType = 'N4d' + elif double == 0 and doubleO == 1 and triple == 0 and benzene == 0 and single == 2: atomType = 'N4d' + elif double == 2 and doubleO == 0 and triple == 0 and benzene == 0 and single == 0: atomType = 'N4dd' + elif double == 0 and doubleO == 2 and triple == 0 and benzene == 0 and single == 0: atomType = 'N4dd' + elif double == 1 and doubleO == 1 and triple == 0 and benzene == 0 and single == 0: atomType = 'N4dd' + elif double == 0 and doubleO == 0 and triple == 1 and benzene == 0 and single == 1: atomType = 'N4t' + elif double == 0 and doubleO == 0 and triple == 0 and benzene == 2 and single == 1: atomType = 'N4b' elif atom.symbol == 'O': if double + doubleO == 0 and triple == 0 and benzene == 0: atomType = 'Os' elif double + doubleO == 1 and triple == 0 and benzene == 0: atomType = 'Od' @@ -277,11 +326,11 @@ def getAtomType(atom, bonds): if double + doubleO == 0 and triple == 0 and benzene == 0: atomType = 'Ss' elif double + doubleO == 1 and triple == 0 and benzene == 0: atomType = 'Sd' elif len(bonds) == 0: atomType = 'Sa' - elif atom.symbol == 'N' or atom.symbol == 'Ar' or atom.symbol == 'He' or atom.symbol == 'Ne': + elif atom.symbol == 'Ar' or atom.symbol == 'He' or atom.symbol == 'Ne': return None # Raise exception if we could not identify the proper atom type if atomType == '': raise AtomTypeError('Unable to determine atom type for atom %s.' % atom) - + return atomTypes[atomType] From 414fc012b59f3b1343d570d91a665ca8194d8524 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Thu, 7 Mar 2013 18:37:39 -0500 Subject: [PATCH 02/64] Adding attributes to the class Atom describing the number of lone electron pairs along with function to manipulate and request their values. --- rmgpy/molecule/molecule.pxd | 9 +++++++ rmgpy/molecule/molecule.py | 54 +++++++++++++++++++++++++++++++++---- 2 files changed, 58 insertions(+), 5 deletions(-) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index 2e4d3dd050..caf965e58c 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -40,6 +40,7 @@ cdef class Atom(Vertex): cdef public short charge cdef public str label cdef public AtomType atomType + cdef public short lonePairs cpdef bint equivalent(self, Vertex other) except -2 @@ -59,6 +60,14 @@ cdef class Atom(Vertex): cpdef decrementRadical(self) + cpdef setLonePairs(self, int lonePairs) + + cpdef incrementLonePairs(self) + + cpdef decrementLonePairs(self) + + cpdef updateCharge(self) + ################################################################################ cpdef object SMILEwriter diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 7d8413a1c5..e07458accb 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -60,6 +60,7 @@ class Atom(Vertex): `spinMultiplicity` ``short`` The spin multiplicity of the atom `charge` ``short`` The formal charge of the atom `label` ``str`` A string label that can be used to tag individual atoms + `lonePairs` ``short`` The number of lone electron pairs =================== =================== ==================================== Additionally, the ``mass``, ``number``, and ``symbol`` attributes of the @@ -67,7 +68,7 @@ class Atom(Vertex): e.g. ``atom.symbol`` instead of ``atom.element.symbol``. """ - def __init__(self, element=None, radicalElectrons=0, spinMultiplicity=1, charge=0, label=''): + def __init__(self, element=None, radicalElectrons=0, spinMultiplicity=1, charge=0, label='', lonePairs=0): Vertex.__init__(self) if isinstance(element, str): self.element = elements.__dict__[element] @@ -78,6 +79,7 @@ def __init__(self, element=None, radicalElectrons=0, spinMultiplicity=1, charge= self.charge = charge self.label = label self.atomType = None + self.lonePairs = lonePairs def __str__(self): """ @@ -187,10 +189,11 @@ def isSpecificCaseOf(self, other): if self.radicalElectrons == radical and self.spinMultiplicity == spin: break else: return False - for charge in atom.charge: - if self.charge == charge: break - else: - return False +# until we have charges and lone pairs in the group values we neglect them here +# for charge in atom.charge: +# if self.charge == charge: break +# else: +# return False return True def copy(self): @@ -209,6 +212,7 @@ def copy(self): a.charge = self.charge a.label = self.label a.atomType = self.atomType + a.lonePairs = self.lonePairs return a def isHydrogen(self): @@ -273,6 +277,42 @@ def decrementRadical(self): else: self.spinMultiplicity = self.radicalElectrons + 1 + def setLonePairs(self,lonePairs): + """ + Set the number of lone electron pairs. + """ + # Set the number of electron pairs + self.lonePairs = lonePairs + if self.lonePairs < 0: + raise ActionError('Unable to update Atom due to setLonePairs : Invalid lone electron pairs set "{0}".'.format(self.setLonePairs)) + + def incrementLonePairs(self): + """ + Update the lone electron pairs pattern as a result of applying a GAIN_PAIR action. + """ + # Set the new lone electron pairs count + self.lonePairs += 1 + if self.lonePairs <= 0: + raise ActionError('Unable to update Atom due to GAIN_PAIR action: Invalid lone electron pairs set "{0}".'.format(self.lonePairs)) + + def decrementLonePairs(self): + """ + Update the lone electron pairs pattern as a result of applying a LOSE_PAIR action. + """ + # Set the new lone electron pairs count + self.lonePairs -= 1 + if self.lonePairs < 0: + raise ActionError('Unable to update Atom due to LOSE_PAIR action: Invalid lone electron pairs set "{0}".'.format(self.lonePairs)) + + def updateCharge(self): + valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'He': 0, 'Ne': 0, 'Ar': 0} + orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5} + valence = valences[self.symbol] + order = 0 + for atom2, bond in self.bonds.items(): + order += orders[bond.order] + self.charge = 8 - valence - order - self.radicalElectrons - 2*self.lonePairs + def applyAction(self, action): """ Update the atom pattern as a result of applying `action`, a tuple @@ -290,6 +330,10 @@ def applyAction(self, action): for i in range(action[2]): self.incrementRadical() elif action[0].upper() == 'LOSE_RADICAL': for i in range(abs(action[2])): self.decrementRadical() + elif action[0].upper() == 'GAIN_PAIR': + for i in range(action[2]): self.incrementLonePairs() + elif action[0].upper() == 'LOSE_PAIR': + for i in range(abs(action[2])): self.decrementLonePairs() else: raise ActionError('Unable to update Atom: Invalid action {0}".'.format(action)) From 5918d4b7d8cfd791b11a514051c968a8815bb38a Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Thu, 7 Mar 2013 18:41:57 -0500 Subject: [PATCH 03/64] Removed a code section responsible for saturating the input adjacency lists with hydrogen atoms from the function fromAdjacencyList requiring that the input adjacency list have include already all hydrogen atoms. Added a new code section calculating and storing the number of lone electron pairs and electric charges for each atom which is possible only if the adjacency list contains all hydrogen atoms. --- rmgpy/molecule/adjlist.py | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index 3817423810..bc7343ae55 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -194,11 +194,10 @@ def fromAdjacencyList(adjlist, group=False): atom1.edges[atom2] = bond atom2.edges[atom1] = bond - # Add explicit hydrogen atoms to complete structure if desired + # Calculate the number of lone pair electrons requiring molecule with all hydrogen atoms present if not group: valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'He': 0, 'Ne': 0, 'Ar': 0} orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5} - newAtoms = [] for atom in atoms: try: valence = valences[atom.symbol] @@ -208,15 +207,11 @@ def fromAdjacencyList(adjlist, group=False): order = 0 for atom2, bond in atom.bonds.items(): order += orders[bond.order] - count = valence - radical - int(order) - for i in range(count): - a = Atom('H', 0, 1, 0, '') - b = Bond(atom, a, 'S') - newAtoms.append(a) - atom.bonds[a] = b - a.bonds[atom] = b - atoms.extend(newAtoms) - + lonePairs = 4 - order - radical + charge = 8 - valence - order - radical - 2*lonePairs + atom.setLonePairs(lonePairs) + atom.updateCharge() + except InvalidAdjacencyListError: print adjlist raise From 7de94672730ead6c5a7741db3e9fca866dc85d2c Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Sun, 10 Mar 2013 00:25:45 -0500 Subject: [PATCH 04/64] Added function searching for resonance isomers by radical-lone electron pair resonance and extend the electronic charge calculation for the special case of hydrogen. --- rmgpy/molecule/adjlist.py | 35 +++++++++++------- rmgpy/molecule/molecule.py | 74 +++++++++++++++++++++++++++++++++++++- 2 files changed, 96 insertions(+), 13 deletions(-) diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index bc7343ae55..87ae07a678 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -199,19 +199,30 @@ def fromAdjacencyList(adjlist, group=False): valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'He': 0, 'Ne': 0, 'Ar': 0} orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5} for atom in atoms: - try: + if not atom.isHydrogen(): + try: + valence = valences[atom.symbol] + except KeyError: + raise InvalidAdjacencyListError('Cannot add hydrogens to adjacency list: Unknown valence for atom "{0}".'.format(atom.symbol)) + radical = atom.radicalElectrons + order = 0 + for atom2, bond in atom.bonds.items(): + order += orders[bond.order] + lonePairs = 4 - order - radical + charge = 8 - valence - order - radical - 2*lonePairs + atom.setLonePairs(lonePairs) + atom.updateCharge() + else: valence = valences[atom.symbol] - except KeyError: - raise InvalidAdjacencyListError('Cannot add hydrogens to adjacency list: Unknown valence for atom "{0}".'.format(atom.symbol)) - radical = atom.radicalElectrons - order = 0 - for atom2, bond in atom.bonds.items(): - order += orders[bond.order] - lonePairs = 4 - order - radical - charge = 8 - valence - order - radical - 2*lonePairs - atom.setLonePairs(lonePairs) - atom.updateCharge() - + radical = atom.radicalElectrons + order = 0 + for atom2, bond in atom.bonds.items(): + order += orders[bond.order] + lonePairs = 1 - order - radical + charge = 2 - valence - order - radical - 2*lonePairs + atom.setLonePairs(lonePairs) + atom.updateCharge() + except InvalidAdjacencyListError: print adjlist raise diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index e07458accb..0a1d98a18e 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -311,7 +311,10 @@ def updateCharge(self): order = 0 for atom2, bond in self.bonds.items(): order += orders[bond.order] - self.charge = 8 - valence - order - self.radicalElectrons - 2*self.lonePairs + if self.isHydrogen(): + self.charge = 2 - valence - order - self.radicalElectrons - 2*self.lonePairs + else: + self.charge = 8 - valence - order - self.radicalElectrons - 2*self.lonePairs def applyAction(self, action): """ @@ -1328,6 +1331,7 @@ def generateResonanceIsomers(self): while index < len(isomers): isomer = isomers[index] newIsomers = isomer.getAdjacentResonanceIsomers() + newIsomers += isomer.getLonePairResonanceIsomers() for newIsomer in newIsomers: newIsomer.updateAtomTypes() # Append to isomer list if unique @@ -1382,6 +1386,52 @@ def getAdjacentResonanceIsomers(self): isomers.append(isomer) return isomers + + def getLonePairResonanceIsomers(self): + """ + Generate all of the resonance isomers formed by one lone electron pair radical shift. + """ + cython.declare(isomers=list, paths=list, index=cython.int, isomer=Molecule) + cython.declare(atom=Atom, atom1=Atom, atom2=Atom) + cython.declare(v1=Vertex, v2=Vertex) + + isomers = [] + + # Radicals + if self.isRadical(): + # Iterate over radicals in structure + for atom in self.vertices: + paths = self.findAllDelocalizationPathsLonePairs(atom) + for atom1, atom2 in paths: + # Adjust to (potentially) new resonance isomer + atom1.decrementRadical() + atom1.incrementLonePairs() + atom1.updateCharge() + atom2.incrementRadical() + atom2.decrementLonePairs() + atom2.updateCharge() + # Make a copy of isomer + isomer = self.copy(deep=True) + # Also copy the connectivity values, since they are the same + # for all resonance forms + for index in range(len(self.vertices)): + v1 = self.vertices[index] + v2 = isomer.vertices[index] + v2.connectivity1 = v1.connectivity1 + v2.connectivity2 = v1.connectivity2 + v2.connectivity3 = v1.connectivity3 + v2.sortingLabel = v1.sortingLabel + # Restore current isomer + atom1.incrementRadical() + atom1.decrementLonePairs() + atom1.updateCharge() + atom2.decrementRadical() + atom2.incrementLonePairs() + atom2.updateCharge() + # Append to isomer list if unique + isomers.append(isomer) + + return isomers def findAllDelocalizationPaths(self, atom1): """ @@ -1405,6 +1455,28 @@ def findAllDelocalizationPaths(self, atom1): if atom1 is not atom3 and (bond23.isDouble() or bond23.isTriple()): paths.append([atom1, atom2, atom3, bond12, bond23]) return paths + + def findAllDelocalizationPathsLonePairs(self, atom1): + """ + Find all the delocalization paths lone electron pairs next to the radical center indicated + by `atom1`. Used to generate resonance isomers. + """ + cython.declare(paths=list) + cython.declare(atom2=Atom, bond12=Bond) + + # No paths if atom1 is not a radical + if atom1.radicalElectrons <= 0: + return [] + + # Find all delocalization paths + paths = [] + for atom2, bond12 in atom1.edges.items(): + # Only single bonds are considered + if bond12.isSingle(): + # Neighboring atom must posses an lone electron pair to loose it + if atom1 is not atom2 and (atom2.lonePairs >= 1): + paths.append([atom1, atom2]) + return paths def getURL(self): """ From 77997a7aa2317445334f3e19d545fb46efb8b3ce Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Sun, 10 Mar 2013 22:31:11 -0400 Subject: [PATCH 05/64] Extended the symmetry number calculation to account for nitrogen atoms. --- rmgpy/molecule/symmetry.py | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/rmgpy/molecule/symmetry.py b/rmgpy/molecule/symmetry.py index 74b31024e1..4e4fd85426 100644 --- a/rmgpy/molecule/symmetry.py +++ b/rmgpy/molecule/symmetry.py @@ -84,6 +84,12 @@ def calculateAtomSymmetryNumber(molecule, atom): elif count == [2, 2]: symmetryNumber *= 2 elif count == [2, 1, 1]: symmetryNumber *= 1 elif count == [1, 1, 1, 1]: symmetryNumber *= 1 + elif single == 3: + # Three single bonds + if count == [3]: symmetryNumber *= 3 + elif count == [2, 1]: symmetryNumber *= 1 + elif count == [1, 1, 1]: symmetryNumber *= 1 + elif single == 2: # Two single bonds if count == [2]: symmetryNumber *= 2 @@ -99,8 +105,13 @@ def calculateAtomSymmetryNumber(molecule, atom): elif atom.radicalElectrons == 2: if single == 2: # Two single bonds - if count == [2]: symmetryNumber *= 2 - + if count == [2]: + symmetryNumber *= 2 + + for groupN in groups: + if groupN.toSMILES() == "[N+](=O)[O-]": + symmetryNumber *= 2 + return symmetryNumber ################################################################################ @@ -276,7 +287,6 @@ def calculateAxisSymmetryNumber(molecule): symmetry_broken=False end_fragments_to_remove = [] for fragment in end_fragments: # a fragment is one end of the axis - # remove the atom that was at the end of the axis and split what's left into groups terminalAtom = None for atom in terminalAtoms: @@ -298,8 +308,11 @@ def calculateAxisSymmetryNumber(molecule): end_fragments_to_remove.append(fragment) continue # next end fragment elif len(groups)==1 and terminalAtom.radicalElectrons == 0: - end_fragments_to_remove.append(fragment) - continue # next end fragment + if terminalAtom.atomType.label == 'N3d': + symmetry_broken = True + else: + end_fragments_to_remove.append(fragment) + continue # next end fragment elif len(groups)==1 and terminalAtom.radicalElectrons != 0: symmetry_broken = True elif len(groups)==2: From 6f799ee9fa2916fb63e5e64fb5ce24ec5d93439b Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 11 Mar 2013 01:54:33 -0400 Subject: [PATCH 06/64] Change removeH from True to False to enforce explicit hydrogens in databases. --- rmgpy/chemkin.py | 2 +- rmgpy/data/base.py | 7 ++++--- rmgpy/data/kinetics/common.py | 8 ++++---- rmgpy/data/kinetics/library.py | 2 +- rmgpy/data/statmech.py | 2 +- rmgpy/data/thermo.py | 3 +-- rmgpy/reaction.py | 4 ++-- rmgpy/species.py | 2 +- 8 files changed, 15 insertions(+), 15 deletions(-) diff --git a/rmgpy/chemkin.py b/rmgpy/chemkin.py index dc0a82683d..d6dc8eb8da 100644 --- a/rmgpy/chemkin.py +++ b/rmgpy/chemkin.py @@ -1286,7 +1286,7 @@ def saveSpeciesDictionary(path, species): """ with open(path, 'w') as f: for spec in species: - f.write(spec.molecule[0].toAdjacencyList(label=getSpeciesIdentifier(spec), removeH=True)) + f.write(spec.molecule[0].toAdjacencyList(label=getSpeciesIdentifier(spec), removeH=False)) f.write('\n') def saveChemkinFile(path, species, reactions, verbose = True): diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index c500db1045..9c8741048f 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -819,7 +819,8 @@ def matchNodeToStructure(self, node, structure, atoms): # Make sure labels actually point to atoms. if center is None or atom is None: return False - if isinstance(center, list): center = center[0] + if isinstance(center, list): + center = center[0] # Semantic check #1: atoms with same label are equivalent elif not atom.isSpecificCaseOf(center): return False @@ -872,7 +873,7 @@ def descendTree(self, structure, atoms, root=None): return None elif not self.matchNodeToStructure(root, structure, atoms): return None - + next = [] for child in root.children: if self.matchNodeToStructure(child, structure, atoms): @@ -1151,7 +1152,7 @@ def saveEntry(self, f, entry, name='entry'): if isinstance(entry.item, Molecule): f.write(' molecule = \n') f.write('"""\n') - f.write(entry.item.toAdjacencyList(removeH=True)) + f.write(entry.item.toAdjacencyList(removeH=False)) f.write('""",\n') elif isinstance(entry.item, Group): f.write(' group = \n') diff --git a/rmgpy/data/kinetics/common.py b/rmgpy/data/kinetics/common.py index 3f62678b86..94f1dc9a34 100644 --- a/rmgpy/data/kinetics/common.py +++ b/rmgpy/data/kinetics/common.py @@ -128,12 +128,12 @@ def sortEfficiencies(efficiencies0): if isinstance(reactant, Molecule): f.write(' reactant{0:d} = \n'.format(i+1)) f.write('"""\n') - f.write(reactant.toAdjacencyList(removeH=True)) + f.write(reactant.toAdjacencyList(removeH=False)) f.write('""",\n') elif isinstance(reactant, Species): f.write(' reactant{0:d} = \n'.format(i+1)) f.write('"""\n') - f.write(reactant.molecule[0].toAdjacencyList(label=reactant.label, removeH=True)) + f.write(reactant.molecule[0].toAdjacencyList(label=reactant.label, removeH=False)) f.write('""",\n') elif isinstance(reactant, Group): f.write(' group{0:d} = \n'.format(i+1)) @@ -146,12 +146,12 @@ def sortEfficiencies(efficiencies0): if isinstance(product, Molecule): f.write(' product{0:d} = \n'.format(i+1)) f.write('"""\n') - f.write(product.toAdjacencyList(removeH=True)) + f.write(product.toAdjacencyList(removeH=False)) f.write('""",\n') elif isinstance(reactant, Species): f.write(' product{0:d} = \n'.format(i+1)) f.write('"""\n') - f.write(product.molecule[0].toAdjacencyList(label=product.label, removeH=True)) + f.write(product.molecule[0].toAdjacencyList(label=product.label, removeH=False)) f.write('""",\n') if not isinstance(entry.item.reactants[0], Group) and not isinstance(entry.item.reactants[0], LogicNode): f.write(' degeneracy = {0:d},\n'.format(entry.item.degeneracy)) diff --git a/rmgpy/data/kinetics/library.py b/rmgpy/data/kinetics/library.py index cf9b880f73..7bd109b074 100644 --- a/rmgpy/data/kinetics/library.py +++ b/rmgpy/data/kinetics/library.py @@ -428,7 +428,7 @@ def writeArrhenius(f, arrhenius): speciesList.sort(key=lambda x: x.label) f = open(os.path.join(path, 'species.txt'), 'w') for species in speciesList: - f.write(species.molecule[0].toAdjacencyList(label=species.label, removeH=True) + "\n") + f.write(species.molecule[0].toAdjacencyList(label=species.label, removeH=False) + "\n") f.close() # Save the high-pressure limit reactions diff --git a/rmgpy/data/statmech.py b/rmgpy/data/statmech.py index d7d2a064bd..e7256101c1 100644 --- a/rmgpy/data/statmech.py +++ b/rmgpy/data/statmech.py @@ -52,7 +52,7 @@ def saveEntry(f, entry): if isinstance(entry.item, Molecule): f.write(' molecule = \n') f.write('"""\n') - f.write(entry.item.toAdjacencyList(removeH=True)) + f.write(entry.item.toAdjacencyList(removeH=False)) f.write('""",\n') elif isinstance(entry.item, Group): f.write(' group = \n') diff --git a/rmgpy/data/thermo.py b/rmgpy/data/thermo.py index f1365423dd..bd92b92d4f 100644 --- a/rmgpy/data/thermo.py +++ b/rmgpy/data/thermo.py @@ -60,7 +60,7 @@ def saveEntry(f, entry): if isinstance(entry.item, Molecule): f.write(' molecule = \n') f.write('"""\n') - f.write(entry.item.toAdjacencyList(removeH=True)) + f.write(entry.item.toAdjacencyList(removeH=False)) f.write('""",\n') elif isinstance(entry.item, Group): f.write(' group = \n') @@ -881,7 +881,6 @@ def __addGroupThermoData(self, thermoData, database, molecule, atom): in the structure `structure`, and add it to the existing thermo data `thermoData`. """ - node0 = database.descendTree(molecule, atom, None) if node0 is None: raise KeyError('Node not found in database.') diff --git a/rmgpy/reaction.py b/rmgpy/reaction.py index b65782aabe..387425fb13 100644 --- a/rmgpy/reaction.py +++ b/rmgpy/reaction.py @@ -149,10 +149,10 @@ def getURL(self): url = "http://rmg.mit.edu/database/kinetics/reaction/" for i,species in enumerate(self.reactants): - adjlist = species.molecule[0].toAdjacencyList(removeH=True) + adjlist = species.molecule[0].toAdjacencyList(removeH=False) url += "reactant{0}={1}__".format(i+1, re.sub('\s+', '%20', adjlist.replace('\n', ';'))) for i,species in enumerate(self.products): - adjlist = species.molecule[0].toAdjacencyList(removeH=True) + adjlist = species.molecule[0].toAdjacencyList(removeH=False) url += "product{0}={1}__".format(i+1, re.sub('\s+', '%20', adjlist.replace('\n', ';'))) return url.strip('_') diff --git a/rmgpy/species.py b/rmgpy/species.py index fa3b6d499b..5b25272d04 100644 --- a/rmgpy/species.py +++ b/rmgpy/species.py @@ -199,7 +199,7 @@ def toAdjacencyList(self): """ Return a string containing each of the molecules' adjacency lists. """ - output = '\n\n'.join([m.toAdjacencyList(label=self.label, removeH=True) for m in self.molecule]) + output = '\n\n'.join([m.toAdjacencyList(label=self.label, removeH=False) for m in self.molecule]) return output def hasStatMech(self): From bd662eca28d0fe37b6e6fc25b191e3eabe86b405 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Thu, 25 Apr 2013 17:21:25 -0400 Subject: [PATCH 07/64] I have added the atomic mass of nitrogen to the function loadGeometry --- rmgpy/cantherm/gaussian.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/rmgpy/cantherm/gaussian.py b/rmgpy/cantherm/gaussian.py index 54015f033a..f589a9ec9b 100644 --- a/rmgpy/cantherm/gaussian.py +++ b/rmgpy/cantherm/gaussian.py @@ -145,6 +145,8 @@ def loadGeometry(self): mass[i] = 1.00782503207 elif number[i] == 6: mass[i] = 12.0 + elif number[i] == 7: + mass[i] = 14.0030740048 elif number[i] == 8: mass[i] = 15.99491461956 else: From 1ef8c6e668cc58f5e48b049fc9cdfa819408e79a Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 6 May 2013 11:17:15 -0400 Subject: [PATCH 08/64] Added the function saveChemkinFileEdge. This function writes for debugging the molecules in the edge into a ChemkinFile similar to the function saveChemkinFile for the core molecules but with additional information. --- rmgpy/rmg/model.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index fc185f9ee1..abe6091f9f 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -1488,3 +1488,16 @@ def saveChemkinFile(self, path, verbose_path, dictionaryPath=None): saveChemkinFile(verbose_path, speciesList, rxnList, verbose = True) if dictionaryPath: saveSpeciesDictionary(dictionaryPath, speciesList) + + def saveChemkinFileEdge(self, path, verbose_path, dictionaryPath=None): + """ + Save a Chemkin file for the current model edge as well as any desired output + species and reactions to `path`. + """ + from rmgpy.chemkin import saveChemkinFile, saveSpeciesDictionaryEdge + speciesList = self.edge.species + self.outputSpeciesList + rxnList = self.edge.reactions + self.outputReactionList + saveChemkinFile(path, speciesList, rxnList, verbose = False) + saveChemkinFile(verbose_path, speciesList, rxnList, verbose = True) + if dictionaryPath: + saveSpeciesDictionaryEdge(dictionaryPath, speciesList) From 379c50c628775201b4f4567a432ad26150d18ae7 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 6 May 2013 13:01:53 -0400 Subject: [PATCH 09/64] Extension for nitrogen molecule.py and molecule.pxd - adding function isNitrogen() to confirm that an atom is a nitrogen atom - updated and renamed the functions identifying lone pair - radical electron shift resonance isomers - added functions identifying resonance isomers with nitrogen atoms having two double bonds (N4dd) and nitrogen atoms having one triple and one single bond (N4ts) - removed isRadical() in generateResonanceIsomers() because it is again used inside getAdjacentResonanceIsomers() and getLonePairRadicalIsomers() and is not required for getN4dd_N4tsResonanceIsomers() --- rmgpy/molecule/molecule.pxd | 8 ++ rmgpy/molecule/molecule.py | 181 +++++++++++++++++++++++++++++++----- 2 files changed, 164 insertions(+), 25 deletions(-) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index caf965e58c..9f48285f01 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -195,7 +195,15 @@ cdef class Molecule(Graph): cpdef list generateResonanceIsomers(self) cpdef list getAdjacentResonanceIsomers(self) + + cpdef list getLonePairRadicalResonanceIsomers(self) + + cpdef list getN4dd_N4tsResonanceIsomers(self) cpdef findAllDelocalizationPaths(self, Atom atom1) + + cpdef findAllDelocalizationPathsLonePairRadical(self, Atom atom1) + + cpdef findAllDelocalizationPathsN4dd_N4ts(self, Atom atom1) cpdef int calculateSymmetryNumber(self) except -1 diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 0a1d98a18e..88b9985895 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -235,6 +235,13 @@ def isCarbon(self): not. """ return self.element.number == 6 + + def isNitrogen(self): + """ + Return ``True`` if the atom represents a nitrogen atom or ``False`` if + not. + """ + return self.element.number == 7 def isOxygen(self): """ @@ -1326,22 +1333,23 @@ def generateResonanceIsomers(self): isomers = [self] # Iterate over resonance isomers - if self.isRadical(): - index = 0 - while index < len(isomers): - isomer = isomers[index] - newIsomers = isomer.getAdjacentResonanceIsomers() - newIsomers += isomer.getLonePairResonanceIsomers() - for newIsomer in newIsomers: - newIsomer.updateAtomTypes() - # Append to isomer list if unique - for isom in isomers: - if isom.isIsomorphic(newIsomer): - break - else: - isomers.append(newIsomer) - # Move to next resonance isomer - index += 1 + index = 0 + while index < len(isomers): + isomer = isomers[index] + newIsomers = isomer.getAdjacentResonanceIsomers() + newIsomers += isomer.getLonePairRadicalResonanceIsomers() + newIsomers += isomer.getN4dd_N4tsResonanceIsomers() + for newIsomer in newIsomers: + newIsomer.updateAtomTypes() + # Append to isomer list if unique + for isom in isomers: + if isom.isIsomorphic(newIsomer): + break + else: + isomers.append(newIsomer) + + # Move to next resonance isomer + index += 1 return isomers @@ -1387,9 +1395,9 @@ def getAdjacentResonanceIsomers(self): return isomers - def getLonePairResonanceIsomers(self): + def getLonePairRadicalResonanceIsomers(self): """ - Generate all of the resonance isomers formed by one lone electron pair radical shift. + Generate all of the resonance isomers formed by lone electron pair - radical shifts. """ cython.declare(isomers=list, paths=list, index=cython.int, isomer=Molecule) cython.declare(atom=Atom, atom1=Atom, atom2=Atom) @@ -1401,7 +1409,7 @@ def getLonePairResonanceIsomers(self): if self.isRadical(): # Iterate over radicals in structure for atom in self.vertices: - paths = self.findAllDelocalizationPathsLonePairs(atom) + paths = self.findAllDelocalizationPathsLonePairRadical(atom) for atom1, atom2 in paths: # Adjust to (potentially) new resonance isomer atom1.decrementRadical() @@ -1432,6 +1440,87 @@ def getLonePairResonanceIsomers(self): isomers.append(isomer) return isomers + + def getN4dd_N4tsResonanceIsomers(self): + """ + Generate all of the resonance isomers formed by shifts between N4dd and N4ts. + """ + cython.declare(isomers=list, paths=list, index=cython.int, isomer=Molecule) + cython.declare(atom=Atom, atom1=Atom, atom2=Atom, atom3=Atom) + cython.declare(bond12=Bond, bond13=Bond) + cython.declare(v1=Vertex, v2=Vertex) + + isomers = [] + + # Iterate over nitrogen atoms in structure + for atom in self.vertices: + paths = self.findAllDelocalizationPathsN4dd_N4ts(atom) + for atom1, atom2, atom3, bond12, bond13, direction in paths: + # from N4dd to N4ts + if direction == 1: + # Adjust to (potentially) new resonance isomer + bond12.decrementOrder() + bond13.incrementOrder() + atom2.incrementLonePairs() + atom3.decrementLonePairs() + atom1.updateCharge() + atom2.updateCharge() + atom3.updateCharge() + # Make a copy of isomer + isomer = self.copy(deep=True) + # Also copy the connectivity values, since they are the same + # for all resonance forms + for index in range(len(self.vertices)): + v1 = self.vertices[index] + v2 = isomer.vertices[index] + v2.connectivity1 = v1.connectivity1 + v2.connectivity2 = v1.connectivity2 + v2.connectivity3 = v1.connectivity3 + v2.sortingLabel = v1.sortingLabel + # Restore current isomer + bond12.incrementOrder() + bond13.decrementOrder() + atom2.decrementLonePairs() + atom3.incrementLonePairs() + atom1.updateCharge() + atom2.updateCharge() + atom3.updateCharge() + # Append to isomer list if unique + isomers.append(isomer) + + # from N4ts to N4dd + if direction == 2: + # Adjust to (potentially) new resonance isomer + bond12.decrementOrder() + bond13.incrementOrder() + atom2.incrementLonePairs() + atom3.decrementLonePairs() + atom1.updateCharge() + atom2.updateCharge() + atom3.updateCharge() + # Make a copy of isomer + isomer = self.copy(deep=True) + # Also copy the connectivity values, since they are the same + # for all resonance forms + for index in range(len(self.vertices)): + v1 = self.vertices[index] + v2 = isomer.vertices[index] + v2.connectivity1 = v1.connectivity1 + v2.connectivity2 = v1.connectivity2 + v2.connectivity3 = v1.connectivity3 + v2.sortingLabel = v1.sortingLabel + # Restore current isomer + bond12.incrementOrder() + bond13.decrementOrder() + atom2.decrementLonePairs() + atom3.incrementLonePairs() + atom1.updateCharge() + atom2.updateCharge() + atom3.updateCharge() + # Append to isomer list if unique + isomers.append(isomer) + + return isomers def findAllDelocalizationPaths(self, atom1): """ @@ -1456,9 +1545,9 @@ def findAllDelocalizationPaths(self, atom1): paths.append([atom1, atom2, atom3, bond12, bond23]) return paths - def findAllDelocalizationPathsLonePairs(self, atom1): + def findAllDelocalizationPathsLonePairRadical(self, atom1): """ - Find all the delocalization paths lone electron pairs next to the radical center indicated + Find all the delocalization paths of lone electron pairs next to the radical center indicated by `atom1`. Used to generate resonance isomers. """ cython.declare(paths=list) @@ -1467,15 +1556,57 @@ def findAllDelocalizationPathsLonePairs(self, atom1): # No paths if atom1 is not a radical if atom1.radicalElectrons <= 0: return [] - + + # In a first step we only consider nitrogen and oxygen atoms as possible radical centers + if not ((atom1.lonePairs == 0 and atom1.isNitrogen()) or(atom1.lonePairs == 2 and atom1.isOxygen())): + return [] + # Find all delocalization paths paths = [] for atom2, bond12 in atom1.edges.items(): # Only single bonds are considered if bond12.isSingle(): - # Neighboring atom must posses an lone electron pair to loose it - if atom1 is not atom2 and (atom2.lonePairs >= 1): - paths.append([atom1, atom2]) + # Neighboring atom must posses a lone electron pair to loose it + if ((atom2.lonePairs == 1 and atom2.isNitrogen()) or (atom2.lonePairs == 3 and atom2.isOxygen())) and (atom2.radicalElectrons == 0): + paths.append([atom1, atom2]) + + return paths + + def findAllDelocalizationPathsN4dd_N4ts(self, atom1): + """ + Find all the resonance structures of nitrogen atoms with two double bonds (N4dd) + and nitrogen atoms with one triple and one single bond (N4ts) + """ + cython.declare(paths=list) + cython.declare(atom2=Atom, bond12=Bond) + + # No paths if atom1 is not nitrogen + if not (atom1.isNitrogen()): + return [] + + # Find all delocalization paths + paths = [] + index_atom_2 = 0 + index_atom_3 = 0 + + for atom2, bond12 in atom1.edges.items(): + index_atom_2 = index_atom_2 + 1 + # Only double bonds are considered + if bond12.isDouble(): + for atom3, bond13 in atom1.edges.items(): + index_atom_3 = index_atom_3 + 1 + # Only double bonds are considered, at the moment we only consider non-radical nitrogen and oxygen atoms + if (bond13.isDouble() and atom3.radicalElectrons == 0 and not atom3.isOxygen() and not atom3.isCarbon() and (index_atom_2 != index_atom_3)): + paths.append([atom1, atom2, atom3, bond12, bond13, 1]) + + for atom2, bond12 in atom1.edges.items(): + # Only triple bonds are considered + if bond12.isTriple(): + for atom3, bond13 in atom1.edges.items(): + # Only single bonds are considered, at the moment we only consider negatively charged nitrogen and oxygen + if (bond13.isSingle() and ((atom3.isNitrogen() and atom3.lonePairs >= 2) or (atom3.isOxygen() and atom3.lonePairs >= 3))): + paths.append([atom1, atom2, atom3, bond12, bond13, 2]) + return paths def getURL(self): From 35b3cbea50d9febe016749576986236e321b8807 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 6 May 2013 13:09:25 -0400 Subject: [PATCH 10/64] Added functions to increment and decrement the number of lone electron pairs - also added the atom type N2d --- rmgpy/molecule/atomtype.pxd | 2 + rmgpy/molecule/atomtype.py | 111 ++++++++++++++++++++---------------- 2 files changed, 64 insertions(+), 49 deletions(-) diff --git a/rmgpy/molecule/atomtype.pxd b/rmgpy/molecule/atomtype.pxd index fe531cd90e..be8ca29cb6 100644 --- a/rmgpy/molecule/atomtype.pxd +++ b/rmgpy/molecule/atomtype.pxd @@ -36,6 +36,8 @@ cdef class AtomType: cdef public list breakBond cdef public list incrementRadical cdef public list decrementRadical + cdef public list incrementLonePair + cdef public list decrementLonePair cpdef bint isSpecificCaseOf(self, AtomType other) diff --git a/rmgpy/molecule/atomtype.py b/rmgpy/molecule/atomtype.py index 054ecfe26a..13c9d398eb 100644 --- a/rmgpy/molecule/atomtype.py +++ b/rmgpy/molecule/atomtype.py @@ -73,6 +73,8 @@ class AtomType: `breakBond` ``list`` The atom type(s) that result when an existing single bond to this atom type is broken `incrementRadical` ``list`` The atom type(s) that result when the number of radical electrons is incremented `decrementRadical` ``list`` The atom type(s) that result when the number of radical electrons is decremented + `incrementLonePair ``list`` The atom type(s) that result when the number of lone electron pairs is incremented + `decrementLonePair` ``list`` The atom type(s) that result when the number of lone electron pairs is decremented =================== =================== ==================================== """ @@ -87,6 +89,8 @@ def __init__(self, label='', generic=None, specific=None): self.breakBond = [] self.incrementRadical = [] self.decrementRadical = [] + self.incrementLonePair = [] + self.decrementLonePair = [] def __repr__(self): return '' % self.label @@ -105,6 +109,8 @@ def __reduce__(self): 'breakBond': self.breakBond, 'incrementRadical': self.incrementRadical, 'decrementRadical': self.decrementRadical, + 'incrementLonePair': self.incrementLonePair, + 'decrementLonePair': self.decrementLonePair, } return (AtomType, (), d) @@ -121,14 +127,18 @@ def __setstate__(self, d): self.breakBond = d['breakBond'] self.incrementRadical = d['incrementRadical'] self.decrementRadical = d['decrementRadical'] + self.incrementLonePair = d['incrementLonePair'] + self.decrementLonePair = d['decrementLonePair'] - def setActions(self, incrementBond, decrementBond, formBond, breakBond, incrementRadical, decrementRadical): + def setActions(self, incrementBond, decrementBond, formBond, breakBond, incrementRadical, decrementRadical, incrementLonePair, decrementLonePair): self.incrementBond = incrementBond self.decrementBond = decrementBond self.formBond = formBond self.breakBond = breakBond self.incrementRadical = incrementRadical self.decrementRadical = decrementRadical + self.incrementLonePair = incrementLonePair + self.decrementLonePair = decrementLonePair def equivalent(self, other): """ @@ -152,14 +162,14 @@ def isSpecificCaseOf(self, other): 'R!H', 'H', 'C','Cs','Cd','Cdd','Ct','CO','Cb','Cbf', 'CS', - 'N','N3s','N3d','N3t','N3b','N4s','N4d','N4dd','N4t','N4b', + 'N','N2d','N3s','N3d','N3t','N3b','N4s','N4d','N4dd','N4t','N4b', 'O','Os','Od','Oa', 'Si','Sis','Sid','Sidd','Sit','SiO','Sib','Sibf', 'S','Ss','Sd','Sa'] ) atomTypes['R!H'] = AtomType(label='R!H', generic=['R'], specific=[ 'C','Cs','Cd','Cdd','Ct','CO','Cb','Cbf', 'CS', - 'N','N3s','N3d','N3t','N3b','N4s','N4d','N4dd','N4t','N4b', + 'N','N2d','N3s','N3d','N3t','N3b','N4s','N4d','N4dd','N4t','N4b', 'O','Os','Od','Oa', 'Si','Sis','Sid','Sidd','Sit','SiO','Sib','Sibf', 'S','Ss','Sd','Sa'] @@ -176,7 +186,8 @@ def isSpecificCaseOf(self, other): atomTypes['Cbf' ] = AtomType('Cbf', generic=['R','R!H','C'], specific=[]) atomTypes['CS' ] = AtomType('CS', generic=['R','R!H','C'], specific=[]) -atomTypes['N' ] = AtomType('N', generic=['R','R!H'], specific=['N3s','N3d','N3t','N3b','N4s','N4d','N4dd','N4t','N4b']) +atomTypes['N' ] = AtomType('N', generic=['R','R!H'], specific=['N2d','N3s','N3d','N3t','N3b','N4s','N4d','N4dd','N4t','N4b']) +atomTypes['N2d' ] = AtomType('N2d', generic=['R','R!H','N'], specific=[]) atomTypes['N3s' ] = AtomType('N3s', generic=['R','R!H','N'], specific=[]) atomTypes['N3d' ] = AtomType('N3d', generic=['R','R!H','N'], specific=[]) atomTypes['N3t' ] = AtomType('N3t', generic=['R','R!H','N'], specific=[]) @@ -206,55 +217,56 @@ def isSpecificCaseOf(self, other): atomTypes['Sd' ] = AtomType('Sd', generic=['R','R!H','S'], specific=[]) atomTypes['Sa' ] = AtomType('Sa', generic=['R','R!H','S'], specific=[]) -atomTypes['R' ].setActions(incrementBond=['R'], decrementBond=['R'], formBond=['R'], breakBond=['R'], incrementRadical=['R'], decrementRadical=['R']) -atomTypes['R!H' ].setActions(incrementBond=['R!H'], decrementBond=['R!H'], formBond=['R!H'], breakBond=['R!H'], incrementRadical=['R!H'], decrementRadical=['R!H']) - -atomTypes['H' ].setActions(incrementBond=[], decrementBond=[], formBond=['H'], breakBond=['H'], incrementRadical=['H'], decrementRadical=['H']) - -atomTypes['C' ].setActions(incrementBond=['C'], decrementBond=['C'], formBond=['C'], breakBond=['C'], incrementRadical=['C'], decrementRadical=['C']) -atomTypes['Cs' ].setActions(incrementBond=['Cd','CO','Cs'], decrementBond=[], formBond=['Cs'], breakBond=['Cs'], incrementRadical=['Cs'], decrementRadical=['Cs']) -atomTypes['Cd' ].setActions(incrementBond=['Cdd','Ct'], decrementBond=['Cs'], formBond=['Cd'], breakBond=['Cd'], incrementRadical=['Cd'], decrementRadical=['Cd']) -atomTypes['Cdd' ].setActions(incrementBond=[], decrementBond=['Cd','CO','CS'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) -atomTypes['Ct' ].setActions(incrementBond=[], decrementBond=['Cd'], formBond=['Ct'], breakBond=['Ct'], incrementRadical=['Ct'], decrementRadical=['Ct']) -atomTypes['CO' ].setActions(incrementBond=['Cdd'], decrementBond=['Cs'], formBond=['CO'], breakBond=['CO'], incrementRadical=['CO'], decrementRadical=['CO']) -atomTypes['CS' ].setActions(incrementBond=['Cdd'], decrementBond=['Cs'], formBond=['CS'], breakBond=['CS'], incrementRadical=['CS'], decrementRadical=['CS']) -atomTypes['Cb' ].setActions(incrementBond=[], decrementBond=[], formBond=['Cb'], breakBond=['Cb'], incrementRadical=['Cb'], decrementRadical=['Cb']) -atomTypes['Cbf' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) - -atomTypes['N' ].setActions(incrementBond=['N'], decrementBond=['N'], formBond=['N'], breakBond=['N'], incrementRadical=['N'], decrementRadical=['N']) -atomTypes['N3s' ].setActions(incrementBond=['N3d','N3s'], decrementBond=[], formBond=['N3s','N4s'], breakBond=['N3s'], incrementRadical=['N3s'], decrementRadical=['N3s']) -atomTypes['N3d' ].setActions(incrementBond=['N3t'], decrementBond=['N3s'], formBond=['N3d','N4d'], breakBond=['N3d'], incrementRadical=['N3d'], decrementRadical=['N3d']) -atomTypes['N3t' ].setActions(incrementBond=[], decrementBond=['N3d'], formBond=['N4t'], breakBond=[], incrementRadical=['N3t'], decrementRadical=['N3t']) -atomTypes['N3b' ].setActions(incrementBond=[], decrementBond=[], formBond=['N3b'], breakBond=['N3b'], incrementRadical=['N3b'], decrementRadical=['N3b']) -atomTypes['N4s' ].setActions(incrementBond=['N4d'], decrementBond=[], formBond=[], breakBond=['N3s'], incrementRadical=['N4s'], decrementRadical=['N4s']) -atomTypes['N4d' ].setActions(incrementBond=['N4dd','N4t'], decrementBond=['N4s'], formBond=[], breakBond=['N3d'], incrementRadical=['N4d'], decrementRadical=['N4d']) -atomTypes['N4dd'].setActions(incrementBond=[], decrementBond=['N3d'], formBond=[], breakBond=[], incrementRadical=['N4dd'], decrementRadical=['N4dd']) -atomTypes['N4t' ].setActions(incrementBond=[], decrementBond=['N3d','N3t'], formBond=[], breakBond=['N3d','N3t'], incrementRadical=['N4t'], decrementRadical=['N4t']) -atomTypes['N4b' ].setActions(incrementBond=[], decrementBond=[], formBond=['N4b'], breakBond=['N4b'], incrementRadical=['N4b'], decrementRadical=['N4b']) - -atomTypes['O' ].setActions(incrementBond=['O'], decrementBond=['O'], formBond=['O'], breakBond=['O'], incrementRadical=['O'], decrementRadical=['O']) -atomTypes['Os' ].setActions(incrementBond=['Od'], decrementBond=[], formBond=['Os'], breakBond=['Os'], incrementRadical=['Os'], decrementRadical=['Os']) -atomTypes['Od' ].setActions(incrementBond=[], decrementBond=['Os'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) -atomTypes['Oa' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) - -atomTypes['Si' ].setActions(incrementBond=['Si'], decrementBond=['Si'], formBond=['Si'], breakBond=['Si'], incrementRadical=['Si'], decrementRadical=['Si']) -atomTypes['Sis' ].setActions(incrementBond=['Sid','SiO'], decrementBond=[], formBond=['Sis'], breakBond=['Sis'], incrementRadical=['Sis'], decrementRadical=['Sis']) -atomTypes['Sid' ].setActions(incrementBond=['Sidd','Sit'], decrementBond=['Sis'], formBond=['Sid'], breakBond=['Sid'], incrementRadical=['Sid'], decrementRadical=['Sid']) -atomTypes['Sidd'].setActions(incrementBond=[], decrementBond=['Sid','SiO'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) -atomTypes['Sit' ].setActions(incrementBond=[], decrementBond=['Sid'], formBond=['Sit'], breakBond=['Sit'], incrementRadical=['Sit'], decrementRadical=['Sit']) -atomTypes['SiO' ].setActions(incrementBond=['Sidd'], decrementBond=['Sis'], formBond=['SiO'], breakBond=['SiO'], incrementRadical=['SiO'], decrementRadical=['SiO']) -atomTypes['Sib' ].setActions(incrementBond=[], decrementBond=[], formBond=['Sib'], breakBond=['Sib'], incrementRadical=['Sib'], decrementRadical=['Sib']) -atomTypes['Sibf'].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) - -atomTypes['S' ].setActions(incrementBond=['S'], decrementBond=['S'], formBond=['S'], breakBond=['S'], incrementRadical=['S'], decrementRadical=['S']) -atomTypes['Ss' ].setActions(incrementBond=['Sd'], decrementBond=[], formBond=['Ss'], breakBond=['Ss'], incrementRadical=['Ss'], decrementRadical=['Ss']) -atomTypes['Sd' ].setActions(incrementBond=[], decrementBond=['Ss'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) -atomTypes['Sa' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[]) +atomTypes['R' ].setActions(incrementBond=['R'], decrementBond=['R'], formBond=['R'], breakBond=['R'], incrementRadical=['R'], decrementRadical=['R'], incrementLonePair=['R'], decrementLonePair=['R']) +atomTypes['R!H' ].setActions(incrementBond=['R!H'], decrementBond=['R!H'], formBond=['R!H'], breakBond=['R!H'], incrementRadical=['R!H'], decrementRadical=['R!H'], incrementLonePair=['R!H'], decrementLonePair=['R!H']) + +atomTypes['H' ].setActions(incrementBond=[], decrementBond=[], formBond=['H'], breakBond=['H'], incrementRadical=['H'], decrementRadical=['H'], incrementLonePair=[], decrementLonePair=[]) + +atomTypes['C' ].setActions(incrementBond=['C'], decrementBond=['C'], formBond=['C'], breakBond=['C'], incrementRadical=['C'], decrementRadical=['C'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Cs' ].setActions(incrementBond=['Cd','CO','Cs'], decrementBond=[], formBond=['Cs'], breakBond=['Cs'], incrementRadical=['Cs'], decrementRadical=['Cs'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Cd' ].setActions(incrementBond=['Cdd','Ct'], decrementBond=['Cs'], formBond=['Cd'], breakBond=['Cd'], incrementRadical=['Cd'], decrementRadical=['Cd'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Cdd' ].setActions(incrementBond=[], decrementBond=['Cd','CO','CS'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Ct' ].setActions(incrementBond=[], decrementBond=['Cd'], formBond=['Ct'], breakBond=['Ct'], incrementRadical=['Ct'], decrementRadical=['Ct'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['CO' ].setActions(incrementBond=['Cdd'], decrementBond=['Cs'], formBond=['CO'], breakBond=['CO'], incrementRadical=['CO'], decrementRadical=['CO'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['CS' ].setActions(incrementBond=['Cdd'], decrementBond=['Cs'], formBond=['CS'], breakBond=['CS'], incrementRadical=['CS'], decrementRadical=['CS'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Cb' ].setActions(incrementBond=[], decrementBond=[], formBond=['Cb'], breakBond=['Cb'], incrementRadical=['Cb'], decrementRadical=['Cb'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Cbf' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) + +atomTypes['N' ].setActions(incrementBond=['N'], decrementBond=['N'], formBond=['N'], breakBond=['N'], incrementRadical=['N'], decrementRadical=['N'], incrementLonePair=['N'], decrementLonePair=['N']) +atomTypes['N2d' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) +atomTypes['N3s' ].setActions(incrementBond=['N3d','N3s'], decrementBond=[], formBond=['N3s','N4s'], breakBond=['N3s'], incrementRadical=['N3s'], decrementRadical=['N3s'], incrementLonePair=['N3s'], decrementLonePair=['N3s']) +atomTypes['N3d' ].setActions(incrementBond=['N3t'], decrementBond=['N3s'], formBond=['N3d','N4d'], breakBond=['N3d'], incrementRadical=['N3d'], decrementRadical=['N3d'], incrementLonePair=['N3d'], decrementLonePair=['N3d']) +atomTypes['N3t' ].setActions(incrementBond=[], decrementBond=['N3d'], formBond=['N4t'], breakBond=[], incrementRadical=['N3t'], decrementRadical=['N3t'], incrementLonePair=['N3t'], decrementLonePair=['N3t']) +atomTypes['N3b' ].setActions(incrementBond=[], decrementBond=[], formBond=['N3b'], breakBond=['N3b'], incrementRadical=['N3b'], decrementRadical=['N3b'], incrementLonePair=['N3b'], decrementLonePair=['N3b']) +atomTypes['N4s' ].setActions(incrementBond=['N4d'], decrementBond=[], formBond=[], breakBond=['N3s'], incrementRadical=['N4s'], decrementRadical=['N4s'], incrementLonePair=['N4s'], decrementLonePair=['N4s']) +atomTypes['N4d' ].setActions(incrementBond=['N4dd','N4t'], decrementBond=['N4s'], formBond=[], breakBond=['N3d'], incrementRadical=['N4d'], decrementRadical=['N4d'], incrementLonePair=['N4d'], decrementLonePair=['N4d']) +atomTypes['N4dd'].setActions(incrementBond=[], decrementBond=['N3d'], formBond=[], breakBond=[], incrementRadical=['N4dd'], decrementRadical=['N4dd'], incrementLonePair=['N4d'], decrementLonePair=['N4d']) +atomTypes['N4t' ].setActions(incrementBond=[], decrementBond=['N3d','N3t'], formBond=[], breakBond=['N3d','N3t'], incrementRadical=['N4t'], decrementRadical=['N4t'], incrementLonePair=['N4t'], decrementLonePair=['N4t']) +atomTypes['N4b' ].setActions(incrementBond=[], decrementBond=[], formBond=['N4b'], breakBond=['N4b'], incrementRadical=['N4b'], decrementRadical=['N4b'], incrementLonePair=['N4b'], decrementLonePair=['N4b']) + +atomTypes['O' ].setActions(incrementBond=['O'], decrementBond=['O'], formBond=['O'], breakBond=['O'], incrementRadical=['O'], decrementRadical=['O'], incrementLonePair=['Os'], decrementLonePair=['Os']) +atomTypes['Os' ].setActions(incrementBond=['Od'], decrementBond=[], formBond=['Os'], breakBond=['Os'], incrementRadical=['Os'], decrementRadical=['Os'], incrementLonePair=['Os'], decrementLonePair=['Os']) +atomTypes['Od' ].setActions(incrementBond=[], decrementBond=['Os'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Oa' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) + +atomTypes['Si' ].setActions(incrementBond=['Si'], decrementBond=['Si'], formBond=['Si'], breakBond=['Si'], incrementRadical=['Si'], decrementRadical=['Si'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Sis' ].setActions(incrementBond=['Sid','SiO'], decrementBond=[], formBond=['Sis'], breakBond=['Sis'], incrementRadical=['Sis'], decrementRadical=['Sis'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Sid' ].setActions(incrementBond=['Sidd','Sit'], decrementBond=['Sis'], formBond=['Sid'], breakBond=['Sid'], incrementRadical=['Sid'], decrementRadical=['Sid'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Sidd'].setActions(incrementBond=[], decrementBond=['Sid','SiO'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Sit' ].setActions(incrementBond=[], decrementBond=['Sid'], formBond=['Sit'], breakBond=['Sit'], incrementRadical=['Sit'], decrementRadical=['Sit'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['SiO' ].setActions(incrementBond=['Sidd'], decrementBond=['Sis'], formBond=['SiO'], breakBond=['SiO'], incrementRadical=['SiO'], decrementRadical=['SiO'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Sib' ].setActions(incrementBond=[], decrementBond=[], formBond=['Sib'], breakBond=['Sib'], incrementRadical=['Sib'], decrementRadical=['Sib'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Sibf'].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) + +atomTypes['S' ].setActions(incrementBond=['S'], decrementBond=['S'], formBond=['S'], breakBond=['S'], incrementRadical=['S'], decrementRadical=['S'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Ss' ].setActions(incrementBond=['Sd'], decrementBond=[], formBond=['Ss'], breakBond=['Ss'], incrementRadical=['Ss'], decrementRadical=['Ss'], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Sd' ].setActions(incrementBond=[], decrementBond=['Ss'], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) +atomTypes['Sa' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) for atomType in atomTypes.values(): for items in [atomType.generic, atomType.specific, atomType.incrementBond, atomType.decrementBond, atomType.formBond, - atomType.breakBond, atomType.incrementRadical, atomType.decrementRadical]: + atomType.breakBond, atomType.incrementRadical, atomType.decrementRadical, atomType.incrementLonePair, atomType.decrementLonePair]: for index in range(len(items)): items[index] = atomTypes[items[index]] @@ -293,6 +305,7 @@ def getAtomType(atom, bonds): atomType = 'H' elif atom.symbol == 'N': if double == 0 and doubleO == 0 and triple == 0 and benzene == 0 and single == 0: atomType = 'N3s' + elif double == 1 and doubleO == 0 and triple == 0 and benzene == 0 and single == 0: atomType = 'N2d' elif double == 0 and doubleO == 0 and triple == 0 and benzene == 0 and single == 1: atomType = 'N3s' elif double == 0 and doubleO == 0 and triple == 0 and benzene == 0 and single == 2: atomType = 'N3s' elif double == 0 and doubleO == 0 and triple == 0 and benzene == 0 and single == 3: atomType = 'N3s' From 14a8a881358222fcfb51a06a2038c53afd1fff77 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 6 May 2013 13:21:26 -0400 Subject: [PATCH 11/64] Extending function isMoleculeForbidden() Because of the many new molecular structures possible by accounting for lone electron pairs, such as ions and new forms of radicals, we have to forbid, at least at the moment, certain structures until we have more thermodynamic and kinetics information for them in our libraries. - Until we have more thermodynamic data of charged molecules we will forbid them - We also forbid positively charged oxygen atoms - We forbid at least at the moment oxygen atoms with bond order larger than two and radical oxygen with a bond order of two --- rmgpy/data/base.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index 9c8741048f..b291b4784a 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -1097,6 +1097,27 @@ def isMoleculeForbidden(self, molecule): initialMap[moleculeLabeledAtoms[label]] = entryLabeledAtoms[label] if molecule.isMappingValid(entry.item, initialMap) and molecule.isSubgraphIsomorphic(entry.item, initialMap): return True + + # Until we have more thermodynamic data of charged molecules we will forbid them + # We also forbid positively charged oxygen atoms + molecule_charge = 0 + for atom in molecule.atoms: + molecule_charge += atom.charge + if atom.isOxygen() and atom.charge > 0: + return True + if molecule_charge != 0: + return True + + # We forbid at least at the moment oxygen atoms with bond order larger than two and radical oxygen with a bond order of two + orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5} + for atom2 in molecule.atoms: + order = 0 + #for bond in atom2.bonds: + for atom3, bond in atom2.bonds.items(): + order += orders[bond.order] + if (atom2.isOxygen() and order>2) or (atom2.isOxygen() and order==2 and atom2.radicalElectrons>0): + return True + return False def loadOld(self, path): From 15618c5aff1b8745eaa6ac926903b47210e8c73d Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 6 May 2013 13:40:15 -0400 Subject: [PATCH 12/64] Introduction of a new reaction family involving lone electron pairs forming single bonds Example reaction: Acatat ion (CH3CO2-) reacting with water (H2O) In both of the two contributing resonance structures of the Acetat ion one of the two oxygen atoms is connected to the carbon atom by a single bond and bears three lone electron pairs while the other oxygen atom is connected by a double bond and bears two lone electron pairs. The oxygen atom with three lone electron pairs can abstract a proton (H+) from a water molecule where one of the three lone electron pairs provides both electrons for the new single bond formed to the new hydrogen atom, resulting in an acetic acid (CH3COOH) and a hydroxide (OH-) molecule. --- rmgpy/data/kinetics/common.py | 1 + 1 file changed, 1 insertion(+) diff --git a/rmgpy/data/kinetics/common.py b/rmgpy/data/kinetics/common.py index 94f1dc9a34..898d06df35 100644 --- a/rmgpy/data/kinetics/common.py +++ b/rmgpy/data/kinetics/common.py @@ -60,6 +60,7 @@ 'SubstitutionS', 'R_Addition_CSm', '1,3_Insertion_RSR', + 'lone_electron_pair_bond', ] # The names of all of the RMG reaction families that are unimolecular From f35e603d03dec65e7f73156820894db1dc455a03 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 6 May 2013 13:46:37 -0400 Subject: [PATCH 13/64] Simple extension to draw.py to draw the lone electron pairs as angular lines around atom symbols - Needs further improvement to account for the direction of bonds to adjacent atoms and the special drawing of small molecules --- rmgpy/molecule/draw.py | 42 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 41 insertions(+), 1 deletion(-) diff --git a/rmgpy/molecule/draw.py b/rmgpy/molecule/draw.py index 08885784b9..e1be33be46 100644 --- a/rmgpy/molecule/draw.py +++ b/rmgpy/molecule/draw.py @@ -1084,7 +1084,7 @@ def __renderAtom(self, symbol, atom, x0, y0, cr, heavyFirst=True): boundingRect = [x1, y1, x2, y2] # Set color for text - if heavyAtom == 'C': cr.set_source_rgba(0.0, 0.0, 0.0, 1.0) + if heavyAtom == 'C': cr.set_source_rgba(0.0, 0.0, 0.0, 1.0) elif heavyAtom == 'N': cr.set_source_rgba(0.0, 0.0, 1.0, 1.0) elif heavyAtom == 'O': cr.set_source_rgba(1.0, 0.0, 0.0, 1.0) elif heavyAtom == 'F': cr.set_source_rgba(0.5, 0.75, 1.0, 1.0) @@ -1253,6 +1253,27 @@ def __renderAtom(self, symbol, atom, x0, y0, cr, heavyFirst=True): cr.move_to(xi, yi - extents[1]) cr.set_source_rgba(0.0, 0.0, 0.0, 1.0) cr.show_text(text) + + # Draw lone electron pairs + for i in range (atom.lonePairs): + cr.new_sub_path() + if i == 0: + x1lp = x-2 + y1lp = y-8 + x2lp = x+2 + y2lp = y-12 + elif i == 1: + x1lp = x+12 + y1lp = y-8 + x2lp = x+8 + y2lp = y-12 + elif i == 2: + x1lp = x-2 + y1lp = y-1 + x2lp = x+2 + y2lp = y+3 + self.__drawLine(cr, x1lp, y1lp, x2lp, y2lp) + elif orientation[0] == 'l' or orientation[0] == 'r': # Draw charges first text = '' @@ -1272,6 +1293,25 @@ def __renderAtom(self, symbol, atom, x0, y0, cr, heavyFirst=True): cr.arc(xi + width/2, yi + 3 * i + 1, 1, 0, 2 * math.pi) cr.set_source_rgba(0.0, 0.0, 0.0, 1.0) cr.fill() + # Draw lone electron pairs + for i in range (atom.lonePairs): + cr.new_sub_path() + if i == 0: + x1lp = x-2 + y1lp = y-8 + x2lp = x+2 + y2lp = y-12 + elif i == 1: + x1lp = x+12 + y1lp = y-8 + x2lp = x+8 + y2lp = y-12 + elif i == 2: + x1lp = x-2 + y1lp = y-1 + x2lp = x+2 + y2lp = y+3 + self.__drawLine(cr, x1lp, y1lp, x2lp, y2lp) # Update bounding rect to ensure atoms are included if boundingRect[0] < self.left: From 87b1243b0fee175708fb8743efddaef6e48f1285 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 6 May 2013 14:08:05 -0400 Subject: [PATCH 14/64] Adding GAIN_PAIR and LOSE_PAIR actions to ReactionRecipe These action will be first used in the lone_electron_pair_bond reaction family - GAIN_PAIR increases the number of lone electron pairs on a specific atom - LOSE_PAIR decreases the number of lone electron pairs on a specific atom --- rmgpy/data/kinetics/family.py | 25 ++++++++++++++++++++++++- 1 file changed, 24 insertions(+), 1 deletion(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index a647eca8c0..da1aeb56bc 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -150,6 +150,8 @@ class ReactionRecipe: BREAK_BOND `center1`, `order`, `center2` break the bond between `center1` and `center2`, which should be of type `order` GAIN_RADICAL `center`, `radical` increase the number of free electrons on `center` by `radical` LOSE_RADICAL `center`, `radical` decrease the number of free electrons on `center` by `radical` + GAIN_PAIR `center`, `pair` increase the number of lone electron pairs on `center` by `pair` + LOSE_PAIR `center`, `pair` decrease the number of lone electron pairs on `center` by `pair` ============= ============================= ================================ The actions are stored as a list in the `actions` attribute. Each action is @@ -186,6 +188,10 @@ def getReverse(self): other.addAction(['GAIN_RADICAL', action[1], action[2]]) elif action[0] == 'GAIN_RADICAL': other.addAction(['LOSE_RADICAL', action[1], action[2]]) + elif action[0] == 'LOSE_PAIR': + other.addAction(['GAIN_PAIR', action[1], action[2]]) + elif action[0] == 'GAIN_PAIR': + other.addAction(['LOSE_PAIR', action[1], action[2]]) return other def __apply(self, struct, doForward, unique): @@ -256,6 +262,23 @@ def __apply(self, struct, doForward, unique): atom.applyAction(['GAIN_RADICAL', label, 1]) elif (action[0] == 'LOSE_RADICAL' and doForward) or (action[0] == 'GAIN_RADICAL' and not doForward): atom.applyAction(['LOSE_RADICAL', label, 1]) + + elif action[0] in ['LOSE_PAIR', 'GAIN_PAIR']: + + label, change = action[1:] + change = int(change) + + # Find associated atom + atom = struct.getLabeledAtom(label) + if atom is None: + raise InvalidActionError('Unable to find atom with label "{0}" while applying reaction recipe.'.format(label)) + + # Apply the action + for i in range(change): + if (action[0] == 'GAIN_PAIR' and doForward) or (action[0] == 'LOSE_PAIR' and not doForward): + atom.applyAction(['GAIN_PAIR', label, 1]) + elif (action[0] == 'LOSE_PAIR' and doForward) or (action[0] == 'GAIN_PAIR' and not doForward): + atom.applyAction(['LOSE_PAIR', label, 1]) else: raise InvalidActionError('Unknown action "' + action[0] + '" encountered.') @@ -579,7 +602,7 @@ def loadRecipe(self, actions): self.forwardRecipe = ReactionRecipe() for action in actions: action[0] = action[0].upper() - assert action[0] in ['CHANGE_BOND','FORM_BOND','BREAK_BOND','GAIN_RADICAL','LOSE_RADICAL'] + assert action[0] in ['CHANGE_BOND','FORM_BOND','BREAK_BOND','GAIN_RADICAL','LOSE_RADICAL','GAIN_PAIR','LOSE_PAIR'] self.forwardRecipe.addAction(action) def loadForbidden(self, label, group, shortDesc='', longDesc='', history=None): From a362192de951deb0e9eb29c0723d3d6177495af5 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 6 May 2013 14:16:29 -0400 Subject: [PATCH 15/64] Adding lone electron pairs to group.py - adding attributes for lone electron pairs - adding functions __gainPair() and __losePair() and handling of the new reaction recipe actions GAIN_PAIR and LOSE_PAIR --- rmgpy/molecule/group.pxd | 5 +++++ rmgpy/molecule/group.py | 38 ++++++++++++++++++++++++++++++++++++-- 2 files changed, 41 insertions(+), 2 deletions(-) diff --git a/rmgpy/molecule/group.pxd b/rmgpy/molecule/group.pxd index 6d5c482bff..f451a77456 100644 --- a/rmgpy/molecule/group.pxd +++ b/rmgpy/molecule/group.pxd @@ -36,6 +36,7 @@ cdef class GroupAtom(Vertex): cdef public list spinMultiplicity cdef public list charge cdef public str label + cdef public list lonePairs cpdef Vertex copy(self) @@ -48,6 +49,10 @@ cdef class GroupAtom(Vertex): cpdef __gainRadical(self, short radical) cpdef __loseRadical(self, short radical) + + cpdef __gainPair(self, short radical) + + cpdef __losePair(self, short radical) cpdef applyAction(self, list action) diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py index 52e9d1dcd6..c5401d02b1 100644 --- a/rmgpy/molecule/group.py +++ b/rmgpy/molecule/group.py @@ -64,6 +64,7 @@ class GroupAtom(Vertex): `spinMultiplicity` ``list`` The allowed spin multiplicities (as short integers) `charge` ``list`` The allowed formal charges (as short integers) `label` ``str`` A string label that can be used to tag individual atoms + `lonePairs` ``list`` The number of lone electron pairs =================== =================== ==================================== Each list represents a logical OR construct, i.e. an atom will match the @@ -73,7 +74,7 @@ class GroupAtom(Vertex): order to match. """ - def __init__(self, atomType=None, radicalElectrons=None, spinMultiplicity=None, charge=None, label=''): + def __init__(self, atomType=None, radicalElectrons=None, spinMultiplicity=None, charge=None, label='', lonePairs=0): Vertex.__init__(self) self.atomType = atomType or [] for index in range(len(self.atomType)): @@ -83,6 +84,7 @@ def __init__(self, atomType=None, radicalElectrons=None, spinMultiplicity=None, self.spinMultiplicity = spinMultiplicity or [] self.charge = charge or [] self.label = label + self.lonePairs = lonePairs or [] def __reduce__(self): """ @@ -98,7 +100,7 @@ def __reduce__(self): atomType = self.atomType if atomType is not None: atomType = [a.label for a in atomType] - return (GroupAtom, (atomType, self.radicalElectrons, self.spinMultiplicity, self.charge, self.label), d) + return (GroupAtom, (atomType, self.radicalElectrons, self.spinMultiplicity, self.charge, self.label, self.lonePairs), d) def __setstate__(self, d): """ @@ -219,6 +221,34 @@ def __loseRadical(self, radical): # Set the new radical electron counts and spin multiplicities self.radicalElectrons = radicalElectrons self.spinMultiplicity = spinMultiplicity + + def __gainPair(self, pair): + """ + Update the atom group as a result of applying a GAIN_PAIR action, + where `pair` specifies the number of lone electron pairs to add. + """ + lonePairs = [] + if any([len(atomType.incrementLonePair) == 0 for atomType in self.atomType]): + raise ActionError('Unable to update GroupAtom due to GAIN_PAIR action: Unknown atom type produced from set "{0}".'.format(self.atomType)) + for lonePairs in zip(self.lonePairs): + lonePairs.append(lonePairs + pair) + # Set the new lone electron pair count + self.lonePairs = lonePairs + + def __losePair(self, pair): + """ + Update the atom group as a result of applying a LOSE_PAIR action, + where `pair` specifies the number of lone electron pairs to remove. + """ + lonePairs = [] + if any([len(atomType.decrementLonePair) == 0 for atomType in self.atomType]): + raise ActionError('Unable to update GroupAtom due to LOSE_PAIR action: Unknown atom type produced from set "{0}".'.format(self.atomType)) + for lonePairs in zip(self.lonePairs): + if lonePairs - pair < 0: + raise ActionError('Unable to update GroupAtom due to LOSE_PAIR action: Invalid lone electron pairs set "{0}".'.format(self.lonePairs)) + lonePairs.append(lonePairs - pair) + # Set the new lone electron pair count + self.lonePairs = lonePairs def applyAction(self, action): """ @@ -237,6 +267,10 @@ def applyAction(self, action): self.__gainRadical(action[2]) elif action[0].upper() == 'LOSE_RADICAL': self.__loseRadical(action[2]) + elif action[0].upper() == 'GAIN_PAIR': + self.__gainPair(action[2]) + elif action[0].upper() == 'LOSE_PAIR': + self.__losePair(action[2]) else: raise ActionError('Unable to update GroupAtom: Invalid action {0}".'.format(action)) From 355bbc1e3f60c3328fa5639c6c9e4020e104a135 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Tue, 7 May 2013 23:07:45 -0400 Subject: [PATCH 16/64] Improving performance of forbidding molecular ions - temporarily, until we have more thermo and kinetics data on ions --- rmgpy/data/base.py | 25 ++++++------------------- 1 file changed, 6 insertions(+), 19 deletions(-) diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index b291b4784a..c85ed5a997 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -1098,25 +1098,12 @@ def isMoleculeForbidden(self, molecule): if molecule.isMappingValid(entry.item, initialMap) and molecule.isSubgraphIsomorphic(entry.item, initialMap): return True - # Until we have more thermodynamic data of charged molecules we will forbid them - # We also forbid positively charged oxygen atoms - molecule_charge = 0 - for atom in molecule.atoms: - molecule_charge += atom.charge - if atom.isOxygen() and atom.charge > 0: - return True - if molecule_charge != 0: - return True - - # We forbid at least at the moment oxygen atoms with bond order larger than two and radical oxygen with a bond order of two - orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5} - for atom2 in molecule.atoms: - order = 0 - #for bond in atom2.bonds: - for atom3, bond in atom2.bonds.items(): - order += orders[bond.order] - if (atom2.isOxygen() and order>2) or (atom2.isOxygen() and order==2 and atom2.radicalElectrons>0): - return True + # Until we have more thermodynamic data of molecular ions we will forbid them + molecule_charge = 0 + for atom in molecule.atoms: + molecule_charge += atom.charge + if molecule_charge != 0: + return True return False From 3e3f1d9f31f24c107220b5d0505388c29a524e02 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Tue, 7 May 2013 23:09:44 -0400 Subject: [PATCH 17/64] Improving performance in atom symmetry number calculation - treating special case of nitrogen in NO2 groups --- rmgpy/molecule/symmetry.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/rmgpy/molecule/symmetry.py b/rmgpy/molecule/symmetry.py index 4e4fd85426..c9dd35e0f2 100644 --- a/rmgpy/molecule/symmetry.py +++ b/rmgpy/molecule/symmetry.py @@ -107,11 +107,12 @@ def calculateAtomSymmetryNumber(molecule, atom): # Two single bonds if count == [2]: symmetryNumber *= 2 - - for groupN in groups: - if groupN.toSMILES() == "[N+](=O)[O-]": - symmetryNumber *= 2 - + + if atom.isNitrogen(): + for groupN in groups: + if groupN.toSMILES() == "[N+](=O)[O-]": + symmetryNumber *= 2 + return symmetryNumber ################################################################################ From 8249d4cb23b1f9bae66a6ae6fce30e8e207551fe Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 15 May 2013 17:53:27 -0400 Subject: [PATCH 18/64] Introducing optional parameter in adjacency lists for lone electron pairs The following adjacency list example describes NO2 with the radical electron on the nitrogen atom with the following parameters where numberLoneElectronPairs is new and optional: atomIndex atomType numberRadicalElectrons numberLoneElectronPairs bonds{} 1 N 1 0 {2,D} {3,S} 2 O 0 2 {1,D} 3 O 0 3 {1,S} --- rmgpy/molecule/adjlist.py | 30 +++++++++++++++++++++++++++--- rmgpy/molecule/group.py | 2 +- 2 files changed, 28 insertions(+), 4 deletions(-) diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index b86d4d5531..eed438a445 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -142,12 +142,33 @@ def fromAdjacencyList(adjlist, group=False): elif e == '4': radicalElectrons.append(4); spinMultiplicity.append(5) index += 1 - + + # Next number defines the number of lone electron pairs (if provided) + lonePairElectrons = -1 + if not group and len(data) > index: + lpState = data[index] + if lpState[0] != '{': + if lpState == '0': + lonePairElectrons = 0 + if lpState == '1': + lonePairElectrons = 1 + if lpState == '2': + lonePairElectrons = 2 + if lpState == '3': + lonePairElectrons = 3 + if lpState == '4': + lonePairElectrons = 4 + index += 1 + else: + lonePairElectrons = -1 + else: + lonePairElectrons = -1 + # Create a new atom based on the above information if group: atom = GroupAtom(atomType, radicalElectrons, spinMultiplicity, [0 for e in radicalElectrons], label) else: - atom = Atom(atomType[0], radicalElectrons[0], spinMultiplicity[0], 0, label) + atom = Atom(atomType[0], radicalElectrons[0], spinMultiplicity[0], 0, label, lonePairElectrons) # Add the atom to the list atoms.append(atom) @@ -204,7 +225,7 @@ def fromAdjacencyList(adjlist, group=False): atom2.edges[atom1] = bond # Calculate the number of lone pair electrons requiring molecule with all hydrogen atoms present - if not group: + if not group and lonePairElectrons == -1: valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'He': 0, 'Ne': 0, 'Ar': 0} orders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5} for atom in atoms: @@ -231,6 +252,9 @@ def fromAdjacencyList(adjlist, group=False): charge = 2 - valence - order - radical - 2*lonePairs atom.setLonePairs(lonePairs) atom.updateCharge() + elif not group: + for atom in atoms: + atom.updateCharge() except InvalidAdjacencyListError: print adjlist diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py index c5401d02b1..c0aa33007c 100644 --- a/rmgpy/molecule/group.py +++ b/rmgpy/molecule/group.py @@ -74,7 +74,7 @@ class GroupAtom(Vertex): order to match. """ - def __init__(self, atomType=None, radicalElectrons=None, spinMultiplicity=None, charge=None, label='', lonePairs=0): + def __init__(self, atomType=None, radicalElectrons=None, spinMultiplicity=None, charge=None, label='', lonePairs=None): Vertex.__init__(self) self.atomType = atomType or [] for index in range(len(self.atomType)): From c019e5dcf421fdf19edc57ea8bd781f3492491e3 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 2 Aug 2013 09:46:29 -0400 Subject: [PATCH 19/64] Added functions to write the edge of the model as CHEMKIN files. --- rmgpy/chemkin.py | 18 ++++++++++++++++++ rmgpy/rmg/main.py | 16 ++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/rmgpy/chemkin.py b/rmgpy/chemkin.py index ca0ad3e6c2..abab3866b9 100644 --- a/rmgpy/chemkin.py +++ b/rmgpy/chemkin.py @@ -1350,6 +1350,24 @@ def saveSpeciesDictionary(path, species): for spec in species: f.write(spec.molecule[0].toAdjacencyList(label=getSpeciesIdentifier(spec), removeH=False)) f.write('\n') + +def saveSpeciesDictionaryEdge(path, species): + """ + Save the given list of `species` as adjacency lists in a text file `path` + on disk. + """ + with open(path, 'w') as f: + for spec in species: + f.write(spec.molecule[0].toAdjacencyList(label=getSpeciesIdentifier(spec), removeH=False)) + f.write('HF298: '+str(spec.getEnthalpy(298)/4.184/1000.0)) + f.write('\n') + f.write('S298: '+str(spec.getEntropy(298)/4.184)) + f.write('\n') + f.write('Cp300: '+str(spec.getHeatCapacity(300)/4.184)) + f.write('\n') + f.write('Cp1500: '+str(spec.getHeatCapacity(1500)/4.184)) + f.write('\n') + f.write('\n') def saveTransportFile(path, species): """ diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 547c0cda2c..e38e43a563 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -480,6 +480,8 @@ def saveEverything(self): self.saveOutputHTML() # Save a Chemkin file containing the current model core self.saveChemkinFile() + # Save a Chemkin file containing the current model edge + self.saveChemkinFileEdge() # Save the restart file if desired if self.saveRestartPeriod or self.done: self.saveRestartFile( os.path.join(self.outputDirectory,'restart.pkl'), @@ -634,6 +636,20 @@ def saveChemkinFile(self): os.unlink(latest_chemkin_path) shutil.copy2(this_chemkin_path,latest_chemkin_path) + def saveChemkinFileEdge(self): + """ + Save the current reaction edge to a Chemkin file. + """ + logging.info('Saving current edge to Chemkin file...') + this_chemkin_path = os.path.join(self.outputDirectory, 'chemkin', 'chem_edge%04i.inp' % len(self.reactionModel.edge.species)) + latest_chemkin_path = os.path.join(self.outputDirectory, 'chemkin','chem_edge.inp') + latest_chemkin_verbose_path = os.path.join(self.outputDirectory, 'chemkin', 'chem_edge_annotated.inp') + latest_dictionary_path = os.path.join(self.outputDirectory, 'chemkin','species_edge_dictionary.txt') + self.reactionModel.saveChemkinFileEdge(this_chemkin_path, latest_chemkin_verbose_path, latest_dictionary_path) + if os.path.exists(latest_chemkin_path): + os.unlink(latest_chemkin_path) + shutil.copy2(this_chemkin_path,latest_chemkin_path) + def saveRestartFile(self, path, reactionModel, delay=0): """ Save a restart file to `path` on disk containing the contents of the From ea24a4f06fd07d62c14b8ee226e566a374c81e60 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 5 Aug 2013 09:22:05 -0400 Subject: [PATCH 20/64] Modified function toAdjacencyList to print the number of lone pair electrons --- rmgpy/molecule/adjlist.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index eed438a445..b3cc06bfdd 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -314,6 +314,7 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False): atomTypes = {} atomElectronStates = {} + atomLonePairs = {} if group: for atom in atomNumbers: # Atom type(s) @@ -331,14 +332,17 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False): # Atom type atomTypes[atom] = '{0}'.format(atom.element.symbol) # Electron state(s) - atomElectronStates[atom] = '{0}'.format(getElectronState(atom.radicalElectrons, atom.spinMultiplicity)) + atomElectronStates[atom] = '{0}'.format(getElectronState(atom.radicalElectrons, atom.spinMultiplicity)) + # Lone Pair(s) + atomLonePairs[atom] = atom.lonePairs # Determine field widths atomNumberWidth = max([len(s) for s in atomNumbers.values()]) + 1 atomLabelWidth = max([len(s) for s in atomLabels.values()]) if atomLabelWidth > 0: atomLabelWidth += 1 atomTypeWidth = max([len(s) for s in atomTypes.values()]) + 1 - atomElectronStateWidth = max([len(s) for s in atomElectronStates.values()]) + atomElectronStateWidth = max([len(s) for s in atomElectronStates.values()]) + 1 + atomLonePairWidth = 1 # Assemble the adjacency list for atom in atoms: @@ -352,6 +356,8 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False): adjlist += '{0:<{1:d}}'.format(atomTypes[atom], atomTypeWidth) # Electron state(s) adjlist += '{0:<{1:d}}'.format(atomElectronStates[atom], atomElectronStateWidth) + # Lone Pair(s) + adjlist += '{0:<{1:d}}'.format(atomLonePairs[atom], atomLonePairWidth) # Bonds list atoms2 = atom.bonds.keys() From 156da3493ceb16792a9322961fc170fc7fa3a3b6 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 5 Aug 2013 09:51:55 -0400 Subject: [PATCH 21/64] A new reaction family describing the conversion of nitro (RNO2) to nitrite (RONO) compounds and vice versa. --- rmgpy/data/kinetics/common.py | 1 + 1 file changed, 1 insertion(+) diff --git a/rmgpy/data/kinetics/common.py b/rmgpy/data/kinetics/common.py index 6152c16be4..d9744367f4 100644 --- a/rmgpy/data/kinetics/common.py +++ b/rmgpy/data/kinetics/common.py @@ -82,6 +82,7 @@ 'intra_substitutionCS_isomerization', 'intra_substitutionS_cyclization', 'intra_substitutionS_isomerization', + 'intra_NO2_ONO_conversion', ] ################################################################################ From 79b36833c296651c431955010368d6773d5ec5c2 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Sun, 11 Aug 2013 13:24:34 -0400 Subject: [PATCH 22/64] Update to prevent problems with nitrates and intra_NO2_ONO_conversion. --- rmgpy/molecule/molecule.pxd | 4 ++++ rmgpy/molecule/molecule.py | 23 +++++++++++++++++++---- 2 files changed, 23 insertions(+), 4 deletions(-) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index b2666f507c..100609ba87 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -99,6 +99,8 @@ cdef class Molecule(Graph): cdef public bint implicitHydrogens cdef public int symmetryNumber + cdef public object rdMol + cdef public int rdMolConfId cdef str _fingerprint cpdef str getFingerprint(self) @@ -123,6 +125,8 @@ cdef class Molecule(Graph): cpdef str getFormula(self) + cpdef short getRadicalCount(self) + cpdef double getMolecularWeight(self) cpdef int getNumAtoms(self, str element=?) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 376edd5f5a..b5ab9542e1 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -36,6 +36,7 @@ """ import cython +import logging import os import re import element as elements @@ -60,6 +61,7 @@ class Atom(Vertex): `spinMultiplicity` ``short`` The spin multiplicity of the atom `charge` ``short`` The formal charge of the atom `label` ``str`` A string label that can be used to tag individual atoms + `coords` `lonePairs` ``short`` The number of lone electron pairs =================== =================== ==================================== @@ -294,6 +296,7 @@ def setLonePairs(self,lonePairs): self.lonePairs = lonePairs if self.lonePairs < 0: raise ActionError('Unable to update Atom due to setLonePairs : Invalid lone electron pairs set "{0}".'.format(self.setLonePairs)) + self.updateCharge() def incrementLonePairs(self): """ @@ -303,6 +306,7 @@ def incrementLonePairs(self): self.lonePairs += 1 if self.lonePairs <= 0: raise ActionError('Unable to update Atom due to GAIN_PAIR action: Invalid lone electron pairs set "{0}".'.format(self.lonePairs)) + self.updateCharge() def decrementLonePairs(self): """ @@ -312,6 +316,7 @@ def decrementLonePairs(self): self.lonePairs -= 1 if self.lonePairs < 0: raise ActionError('Unable to update Atom due to LOSE_PAIR action: Invalid lone electron pairs set "{0}".'.format(self.lonePairs)) + self.updateCharge() def updateCharge(self): valences = {'H': 1, 'C': 4, 'O': 2, 'N': 3, 'S': 2, 'Si': 4, 'He': 0, 'Ne': 0, 'Ar': 0} @@ -670,12 +675,22 @@ def getMolecularWeight(self): """ Return the molecular weight of the molecule in kg/mol. """ + cython.declare(atom=Atom, mass=cython.double) mass = 0 - H = elements.getElement('H') for atom in self.vertices: mass += atom.element.mass return mass + def getRadicalCount(self): + """ + Return the number of unpaired electrons. + """ + cython.declare(atom=Atom, radicals=cython.short) + radicals = 0 + for atom in self.vertices: + radicals += atom.radicalElectrons + return radicals + def getNumAtoms(self, element = None): """ Return the number of atoms in molecule. If element is given, ie. "H" or "C", @@ -742,7 +757,7 @@ def deleteHydrogens(self): connectivity values. If there's nothing but hydrogens, it does nothing. It destroys information; be careful with it. """ - cython.declare(atom=Atom, neighbor=Atom, hydrogens=list) + cython.declare(atom=Atom, hydrogens=list) # Check that the structure contains at least one heavy atom for atom in self.vertices: if not atom.isHydrogen(): @@ -753,7 +768,6 @@ def deleteHydrogens(self): hydrogens = [] for atom in self.vertices: if atom.isHydrogen() and atom.label == '': - neighbor = atom.edges.keys()[0] hydrogens.append(atom) # Remove the hydrogen atoms from the structure for atom in hydrogens: @@ -1338,6 +1352,7 @@ def generateResonanceIsomers(self): index = 0 while index < len(isomers): isomer = isomers[index] + newIsomers = isomer.getAdjacentResonanceIsomers() newIsomers += isomer.getLonePairRadicalResonanceIsomers() newIsomers += isomer.getN4dd_N4tsResonanceIsomers() @@ -1621,4 +1636,4 @@ def getURL(self): adjlist = self.toAdjacencyList(removeH=True) url += "{0}".format(re.sub('\s+', '%20', adjlist.replace('\n', ';'))) return url.strip('_') - + From 240b98fcd221de2b5fc5bc1633adb78d7689e7c1 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 16 Oct 2013 21:59:49 -0400 Subject: [PATCH 23/64] Wrapping RDKit for non-standard valency Because RDKit cannot handle non-standard valency like tetra-valent nitrogen RMG is falling back to previous drawing functions and replaces RDKit's toSMILE with getFormula. Temporarily fixes #156 --- rmgpy/molecule/draw.py | 563 +++++++++++++++++++++++++++++++++++-- rmgpy/molecule/molecule.py | 5 + 2 files changed, 543 insertions(+), 25 deletions(-) diff --git a/rmgpy/molecule/draw.py b/rmgpy/molecule/draw.py index 29ba822928..eb961c2812 100644 --- a/rmgpy/molecule/draw.py +++ b/rmgpy/molecule/draw.py @@ -279,36 +279,549 @@ def __generateCoordinates(self): """ atoms = self.molecule.atoms Natoms = len(atoms) + flag_nitrogen = 0 + + for atom in self.molecule.atoms: + if atom.isNitrogen(): + flag_nitrogen = 1 # Initialize array of coordinates self.coordinates = coordinates = numpy.zeros((Natoms, 2)) - # Use rdkit 2D coordinate generation: - - # Generate the RDkit molecule from the RDkit molecule, use geometry - # in order to match the atoms in the rdmol with the atoms in the - # RMG molecule (which is required to extract coordinates). - self.geometry = Geometry(None, None, self.molecule, None) - - rdmol, rdAtomIdx = self.geometry.rd_build() - AllChem.Compute2DCoords(rdmol) - - # Extract the coordinates from each atom. - for atom in atoms: - index = rdAtomIdx[atom] - point = rdmol.GetConformer(0).GetAtomPosition(index) - coordinates[index,:]= [point.x*0.6, point.y*0.6] - - # RDKit generates some molecules more vertically than horizontally, - # Especially linear ones. This will reflect any molecule taller than - # it is wide across the line y=x - ranges = numpy.ptp(coordinates, axis = 0) - if ranges[1] > ranges[0]: - temp = numpy.copy(coordinates) - coordinates[:,0] = temp[:,1] - coordinates[:,1] = temp[:,0] + if flag_nitrogen == 1: + # If there are only one or two atoms to draw, then determining the + # coordinates is trivial + if Natoms == 1: + self.coordinates[0,:] = [0.0, 0.0] + return self.coordinates + elif Natoms == 2: + self.coordinates[0,:] = [-0.5, 0.0] + self.coordinates[1,:] = [0.5, 0.0] + return self.coordinates + + if len(self.cycles) > 0: + # Cyclic molecule + backbone = self.__findCyclicBackbone() + self.__generateRingSystemCoordinates(backbone) + # Flatten backbone so that it contains a list of the atoms in the + # backbone, rather than a list of the cycles in the backbone + backbone = list(set([atom for cycle in backbone for atom in cycle])) + else: + # Straight chain molecule + backbone = self.__findStraightChainBackbone() + self.__generateStraightChainCoordinates(backbone) + + # If backbone is linear, then rotate so that the bond is parallel to the + # horizontal axis + vector0 = coordinates[atoms.index(backbone[1]),:] - coordinates[atoms.index(backbone[0]),:] + for i in range(2, len(backbone)): + vector = coordinates[atoms.index(backbone[i]),:] - coordinates[atoms.index(backbone[i-1]),:] + if numpy.linalg.norm(vector - vector0) > 1e-4: + break + else: + angle = math.atan2(vector0[0], vector0[1]) - math.pi / 2 + rot = numpy.array([[math.cos(angle), math.sin(angle)], [-math.sin(angle), math.cos(angle)]], numpy.float64) + coordinates = numpy.dot(coordinates, rot) + + # Center backbone at origin + xmin = numpy.min(coordinates[:,0]) + xmax = numpy.max(coordinates[:,0]) + ymin = numpy.min(coordinates[:,1]) + ymax = numpy.max(coordinates[:,1]) + xmid = 0.5 * (xmax + xmin) + ymid = 0.5 * (ymax + ymin) + for atom in backbone: + index = atoms.index(atom) + coordinates[index,0] -= xmid + coordinates[index,1] -= ymid + + # We now proceed by calculating the coordinates of the functional groups + # attached to the backbone + # Each functional group is independent, although they may contain further + # branching and cycles + # In general substituents should try to grow away from the origin to + # minimize likelihood of overlap + self.__generateNeighborCoordinates(backbone) + + return coordinates + + else: + + # Use rdkit 2D coordinate generation: + + # Generate the RDkit molecule from the RDkit molecule, use geometry + # in order to match the atoms in the rdmol with the atoms in the + # RMG molecule (which is required to extract coordinates). + self.geometry = Geometry(None, None, self.molecule, None) + + rdmol, rdAtomIdx = self.geometry.rd_build() + AllChem.Compute2DCoords(rdmol) + + # Extract the coordinates from each atom. + for atom in atoms: + index = rdAtomIdx[atom] + point = rdmol.GetConformer(0).GetAtomPosition(index) + coordinates[index,:]= [point.x*0.6, point.y*0.6] + + # RDKit generates some molecules more vertically than horizontally, + # Especially linear ones. This will reflect any molecule taller than + # it is wide across the line y=x + ranges = numpy.ptp(coordinates, axis = 0) + if ranges[1] > ranges[0]: + temp = numpy.copy(coordinates) + coordinates[:,0] = temp[:,1] + coordinates[:,1] = temp[:,0] + + return coordinates + + def __findCyclicBackbone(self): + """ + Return a set of atoms to use as the "backbone" of the molecule. For + cyclics this is simply the largest ring system. + """ + count = [len(set([atom for ring in ringSystem for atom in ring])) for ringSystem in self.ringSystems] + index = 0 + for i in range(1, len(self.ringSystems)): + if count[i] > count[index]: + index = i + return self.ringSystems[index] + + def __findStraightChainBackbone(self): + """ + Return a set of atoms to use as the "backbone" of the molecule. For + non-cyclics this is the largest straight chain between atoms. If carbon + atoms are present, then we define the backbone only in terms of them. + """ + # Find the terminal atoms - those that only have one explicit bond + terminalAtoms = [atom for atom in self.molecule.atoms if len(atom.bonds) == 1] + assert len(terminalAtoms) >= 2 + + # Starting from each terminal atom, find the longest straight path to + # another terminal + # The longest found is the backbone + backbone = [] + paths = [] + for atom in terminalAtoms: + paths.extend(self.__findStraightChainPaths([atom])) + + # Remove any paths that don't end in a terminal atom + # (I don't think this should remove any!) + paths = [path for path in paths if path[-1] in terminalAtoms] + + # Remove all paths shorter than the maximum + length = max([len(path) for path in paths]) + paths = [path for path in paths if len(path) == length] + + # Prefer the paths with the most carbon atoms + carbons = [sum([1 for atom in path if atom.isCarbon()]) for path in paths] + maxCarbons = max(carbons) + paths = [path for path, carbon in zip(paths, carbons) if carbon == maxCarbons] + + # At this point we could choose any remaining path, so simply choose the first + backbone = paths[0] + + assert len(backbone) > 1 + assert backbone[0] in terminalAtoms + assert backbone[-1] in terminalAtoms + + return backbone + + def __findStraightChainPaths(self, atoms0): + """ + Finds the paths containing the list of atoms `atoms0` in the + current molecule. The atoms are assumed to already be in a path, with + ``atoms0[0]`` being a terminal atom. + """ + atom1 = atoms0[-1] + paths = [] + for atom2 in atom1.bonds: + if atom2 not in atoms0: + atoms = atoms0[:] + atoms.append(atom2) + if not self.molecule.isAtomInCycle(atom2): + paths.extend(self.__findStraightChainPaths(atoms)) + if len(paths) == 0: + paths.append(atoms0[:]) + return paths + + def __generateRingSystemCoordinates(self, atoms): + """ + For a ring system composed of the given cycles of `atoms`, update the + coordinates of each atom in the system. + """ + coordinates = self.coordinates + atoms = atoms[:] + processed = [] + + # Lay out largest cycle in ring system first + cycle = atoms[0] + for cycle0 in atoms[1:]: + if len(cycle0) > len(cycle): + cycle = cycle0 + angle = - 2 * math.pi / len(cycle) + radius = 1.0 / (2 * math.sin(math.pi / len(cycle))) + for i, atom in enumerate(cycle): + index = self.molecule.atoms.index(atom) + coordinates[index,:] = [math.cos(math.pi / 2 + i * angle), math.sin(math.pi / 2 + i * angle)] + coordinates[index,:] *= radius + atoms.remove(cycle) + processed.append(cycle) + + # If there are other cycles, then try to lay them out as well + while len(atoms) > 0: + + # Find the largest cycle that shares one or two atoms with a ring that's + # already been processed + cycle = None + for cycle0 in atoms: + for cycle1 in processed: + count = sum([1 for atom in cycle0 if atom in cycle1]) + if (count == 1 or count == 2): + if cycle is None or len(cycle0) > len(cycle): cycle = cycle0 + cycle0 = cycle1 + atoms.remove(cycle) + + # Shuffle atoms in cycle such that the common atoms come first + # Also find the average center of the processed cycles that touch the + # current cycles + found = False + commonAtoms = [] + count = 0 + center0 = numpy.zeros(2, numpy.float64) + for cycle1 in processed: + found = False + for atom in cycle1: + if atom in cycle and atom not in commonAtoms: + commonAtoms.append(atom) + found = True + if found: + center1 = numpy.zeros(2, numpy.float64) + for atom in cycle1: + center1 += coordinates[cycle1.index(atom),:] + center1 /= len(cycle1) + center0 += center1 + count += 1 + center0 /= count + + if len(commonAtoms) > 1: + index0 = cycle.index(commonAtoms[0]) + index1 = cycle.index(commonAtoms[1]) + if (index0 == 0 and index1 == len(cycle) - 1) or (index1 == 0 and index0 == len(cycle) - 1): + cycle = cycle[-1:] + cycle[0:-1] + if cycle.index(commonAtoms[1]) < cycle.index(commonAtoms[0]): + cycle.reverse() + index = cycle.index(commonAtoms[0]) + cycle = cycle[index:] + cycle[0:index] + + # Determine center of cycle based on already-assigned positions of + # common atoms (which won't be changed) + if len(commonAtoms) == 1 or len(commonAtoms) == 2: + # Center of new cycle is reflection of center of adjacent cycle + # across common atom or bond + center = numpy.zeros(2, numpy.float64) + for atom in commonAtoms: + center += coordinates[self.molecule.atoms.index(atom),:] + center /= len(commonAtoms) + vector = center - center0 + center += vector + radius = 1.0 / (2 * math.sin(math.pi / len(cycle))) + + else: + # Use any three points to determine the point equidistant from these + # three; this is the center + index0 = self.molecule.atoms.index(commonAtoms[0]) + index1 = self.molecule.atoms.index(commonAtoms[1]) + index2 = self.molecule.atoms.index(commonAtoms[2]) + A = numpy.zeros((2,2), numpy.float64) + b = numpy.zeros((2), numpy.float64) + A[0,:] = 2 * (coordinates[index1,:] - coordinates[index0,:]) + A[1,:] = 2 * (coordinates[index2,:] - coordinates[index0,:]) + b[0] = coordinates[index1,0]**2 + coordinates[index1,1]**2 - coordinates[index0,0]**2 - coordinates[index0,1]**2 + b[1] = coordinates[index2,0]**2 + coordinates[index2,1]**2 - coordinates[index0,0]**2 - coordinates[index0,1]**2 + center = numpy.linalg.solve(A, b) + radius = numpy.linalg.norm(center - coordinates[index0,:]) + + startAngle = 0.0; endAngle = 0.0 + if len(commonAtoms) == 1: + # We will use the full 360 degrees to place the other atoms in the cycle + startAngle = math.atan2(-vector[1], vector[0]) + endAngle = startAngle + 2 * math.pi + elif len(commonAtoms) >= 2: + # Divide other atoms in cycle equally among unused angle + vector = coordinates[cycle.index(commonAtoms[-1]),:] - center + startAngle = math.atan2(vector[1], vector[0]) + vector = coordinates[cycle.index(commonAtoms[0]),:] - center + endAngle = math.atan2(vector[1], vector[0]) + + # Place remaining atoms in cycle + if endAngle < startAngle: + endAngle += 2 * math.pi + dAngle = (endAngle - startAngle) / (len(cycle) - len(commonAtoms) + 1) + else: + endAngle -= 2 * math.pi + dAngle = (endAngle - startAngle) / (len(cycle) - len(commonAtoms) + 1) + + count = 1 + for i in range(len(commonAtoms), len(cycle)): + angle = startAngle + count * dAngle + index = self.molecule.atoms.index(cycle[i]) + # Check that we aren't reassigning any atom positions + # This version assumes that no atoms belong at the origin, which is + # usually fine because the first ring is centered at the origin + if numpy.linalg.norm(coordinates[index,:]) < 1e-4: + vector = numpy.array([math.cos(angle), math.sin(angle)], numpy.float64) + coordinates[index,:] = center + radius * vector + count += 1 + + # We're done assigning coordinates for this cycle, so mark it as processed + processed.append(cycle) + + def __generateStraightChainCoordinates(self, atoms): + """ + Update the coordinates for the linear straight chain of `atoms` in + the current molecule. + """ + coordinates = self.coordinates + + # First atom goes at origin + index0 = self.molecule.atoms.index(atoms[0]) + coordinates[index0,:] = [0.0, 0.0] + + # Second atom goes on x-axis (for now; this could be improved!) + index1 = self.molecule.atoms.index(atoms[1]) + vector = numpy.array([1.0, 0.0], numpy.float64) + if atoms[0].bonds[atoms[1]].isTriple(): + rotatePositive = False + else: + rotatePositive = True + rot = numpy.array([[math.cos(-math.pi / 6), math.sin(-math.pi / 6)], [-math.sin(-math.pi / 6), math.cos(-math.pi / 6)]], numpy.float64) + vector = numpy.array([1.0, 0.0], numpy.float64) + vector = numpy.dot(rot, vector) + coordinates[index1,:] = coordinates[index0,:] + vector + + # Other atoms + for i in range(2, len(atoms)): + atom0 = atoms[i-2] + atom1 = atoms[i-1] + atom2 = atoms[i] + index1 = self.molecule.atoms.index(atom1) + index2 = self.molecule.atoms.index(atom2) + bond0 = atom0.bonds[atom1] + bond = atom1.bonds[atom2] + # Angle of next bond depends on the number of bonds to the start atom + numBonds = len(atom1.bonds) + if numBonds == 2: + if (bond0.isTriple() or bond.isTriple()) or (bond0.isDouble() and bond.isDouble()): + # Rotate by 0 degrees towards horizontal axis (to get angle of 180) + angle = 0.0 + else: + # Rotate by 60 degrees towards horizontal axis (to get angle of 120) + angle = math.pi / 3 + elif numBonds == 3: + # Rotate by 60 degrees towards horizontal axis (to get angle of 120) + angle = math.pi / 3 + elif numBonds == 4: + # Rotate by 0 degrees towards horizontal axis (to get angle of 90) + angle = 0.0 + elif numBonds == 5: + # Rotate by 36 degrees towards horizontal axis (to get angle of 144) + angle = math.pi / 5 + elif numBonds == 6: + # Rotate by 0 degrees towards horizontal axis (to get angle of 180) + angle = 0.0 + # Determine coordinates for atom + if angle != 0: + if not rotatePositive: angle = -angle + rot = numpy.array([[math.cos(angle), math.sin(angle)], [-math.sin(angle), math.cos(angle)]], numpy.float64) + vector = numpy.dot(rot, vector) + rotatePositive = not rotatePositive + coordinates[index2,:] = coordinates[index1,:] + vector + + def __generateNeighborCoordinates(self, backbone): + """ + Recursively update the coordinates for the atoms immediately adjacent + to the atoms in the molecular `backbone`. + """ + atoms = self.molecule.atoms + coordinates = self.coordinates + + for i in range(len(backbone)): + atom0 = backbone[i] + index0 = atoms.index(atom0) + + # Determine bond angles of all previously-determined bond locations for + # this atom + bondAngles = [] + for atom1 in atom0.bonds: + index1 = atoms.index(atom1) + if atom1 in backbone: + vector = coordinates[index1,:] - coordinates[index0,:] + angle = math.atan2(vector[1], vector[0]) + bondAngles.append(angle) + bondAngles.sort() + + bestAngle = 2 * math.pi / len(atom0.bonds) + regular = True + for angle1, angle2 in zip(bondAngles[0:-1], bondAngles[1:]): + if all([abs(angle2 - angle1 - (i+1) * bestAngle) > 1e-4 for i in range(len(atom0.bonds))]): + regular = False + + if regular: + # All the bonds around each atom are equally spaced + # We just need to fill in the missing bond locations + + # Determine rotation angle and matrix + rot = numpy.array([[math.cos(bestAngle), -math.sin(bestAngle)], [math.sin(bestAngle), math.cos(bestAngle)]], numpy.float64) + # Determine the vector of any currently-existing bond from this atom + vector = None + for atom1 in atom0.bonds: + index1 = atoms.index(atom1) + if atom1 in backbone or numpy.linalg.norm(coordinates[index1,:]) > 1e-4: + vector = coordinates[index1,:] - coordinates[index0,:] + + # Iterate through each neighboring atom to this backbone atom + # If the neighbor is not in the backbone and does not yet have + # coordinates, then we need to determine coordinates for it + for atom1 in atom0.bonds: + if atom1 not in backbone and numpy.linalg.norm(coordinates[atoms.index(atom1),:]) < 1e-4: + occupied = True; count = 0 + # Rotate vector until we find an unoccupied location + while occupied and count < len(atom0.bonds): + count += 1; occupied = False + vector = numpy.dot(rot, vector) + for atom2 in atom0.bonds: + index2 = atoms.index(atom2) + if numpy.linalg.norm(coordinates[index2,:] - coordinates[index0,:] - vector) < 1e-4: + occupied = True + coordinates[atoms.index(atom1),:] = coordinates[index0,:] + vector + self.__generateFunctionalGroupCoordinates(atom0, atom1) + + else: + + # The bonds are not evenly spaced (e.g. due to a ring) + # We place all of the remaining bonds evenly over the reflex angle + startAngle = max(bondAngles) + endAngle = min(bondAngles) + if 0.0 < endAngle - startAngle < math.pi: endAngle += 2 * math.pi + elif 0.0 > endAngle - startAngle > -math.pi: startAngle -= 2 * math.pi + dAngle = (endAngle - startAngle) / (len(atom0.bonds) - len(bondAngles) + 1) + + index = 1 + for atom1 in atom0.bonds: + if atom1 not in backbone and numpy.linalg.norm(coordinates[atoms.index(atom1),:]) < 1e-4: + angle = startAngle + index * dAngle + index += 1 + vector = numpy.array([math.cos(angle), math.sin(angle)], numpy.float64) + vector /= numpy.linalg.norm(vector) + coordinates[atoms.index(atom1),:] = coordinates[index0,:] + vector + self.__generateFunctionalGroupCoordinates(atom0, atom1) + + def __generateFunctionalGroupCoordinates(self, atom0, atom1): + """ + For the functional group starting with the bond from `atom0` to `atom1`, + generate the coordinates of the rest of the functional group. `atom0` is + treated as if a terminal atom. `atom0` and `atom1` must already have their + coordinates determined. `atoms` is a list of the atoms to be drawn, `bonds` + is a dictionary of the bonds to draw, and `coordinates` is an array of the + coordinates for each atom to be drawn. This function is designed to be + recursive. + """ + + atoms = self.molecule.atoms + coordinates = self.coordinates + + index0 = atoms.index(atom0) + index1 = atoms.index(atom1) + + # Determine the vector of any currently-existing bond from this atom + # (We use the bond to the previous atom here) + vector = coordinates[index0,:] - coordinates[index1,:] + bondAngle = math.atan2(vector[1], vector[0]) + + # Check to see if atom1 is in any cycles in the molecule + ringSystem = None + for ringSys in self.ringSystems: + if any([atom1 in ring for ring in ringSys]): + ringSystem = ringSys + + if ringSystem is not None: + # atom1 is part of a ring system, so we need to process the entire + # ring system at once + + # Generate coordinates for all atoms in the ring system + self.__generateRingSystemCoordinates(ringSystem) + + cycleAtoms = list(set([atom for ring in ringSystem for atom in ring])) + + coordinates_cycle = numpy.zeros_like(self.coordinates) + for atom in cycleAtoms: + coordinates_cycle[atoms.index(atom),:] = coordinates[atoms.index(atom),:] + + # Rotate the ring system coordinates so that the line connecting atom1 + # and the center of mass of the ring is parallel to that between + # atom0 and atom1 + center = numpy.zeros(2, numpy.float64) + for atom in cycleAtoms: + center += coordinates_cycle[atoms.index(atom),:] + center /= len(cycleAtoms) + vector0 = center - coordinates_cycle[atoms.index(atom1),:] + angle = math.atan2(vector[1] - vector0[1], vector[0] - vector0[0]) + rot = numpy.array([[math.cos(angle), -math.sin(angle)], [math.sin(angle), math.cos(angle)]], numpy.float64) + coordinates_cycle = numpy.dot(coordinates_cycle, rot) - return coordinates + # Translate the ring system coordinates to the position of atom1 + coordinates_cycle += coordinates[atoms.index(atom1),:] - coordinates_cycle[atoms.index(atom1),:] + for atom in cycleAtoms: + coordinates[atoms.index(atom),:] = coordinates_cycle[atoms.index(atom),:] + + # Generate coordinates for remaining neighbors of ring system, + # continuing to recurse as needed + self.__generateNeighborCoordinates(cycleAtoms) + + else: + # atom1 is not in any rings, so we can continue as normal + + # Determine rotation angle and matrix + numBonds = len(atom1.bonds) + angle = 0.0 + if numBonds == 2: + bond0, bond = atom1.bonds.values() + if (bond0.isTriple() or bond.isTriple()) or (bond0.isDouble() and bond.isDouble()): + angle = math.pi + else: + angle = 2 * math.pi / 3 + # Make sure we're rotating such that we move away from the origin, + # to discourage overlap of functional groups + rot1 = numpy.array([[math.cos(angle), -math.sin(angle)], [math.sin(angle), math.cos(angle)]], numpy.float64) + rot2 = numpy.array([[math.cos(angle), math.sin(angle)], [-math.sin(angle), math.cos(angle)]], numpy.float64) + vector1 = coordinates[index1,:] + numpy.dot(rot1, vector) + vector2 = coordinates[index1,:] + numpy.dot(rot2, vector) + if bondAngle < -0.5 * math.pi or bondAngle > 0.5 * math.pi: + angle = abs(angle) + else: + angle = -abs(angle) + else: + angle = 2 * math.pi / numBonds + rot = numpy.array([[math.cos(angle), -math.sin(angle)], [math.sin(angle), math.cos(angle)]], numpy.float64) + + # Iterate through each neighboring atom to this backbone atom + # If the neighbor is not in the backbone, then we need to determine + # coordinates for it + for atom, bond in atom1.bonds.iteritems(): + if atom is not atom0: + occupied = True; count = 0 + # Rotate vector until we find an unoccupied location + while occupied and count < len(atom1.bonds): + count += 1; occupied = False + vector = numpy.dot(rot, vector) + for atom2 in atom1.bonds: + index2 = atoms.index(atom2) + if numpy.linalg.norm(coordinates[index2,:] - coordinates[index1,:] - vector) < 1e-4: + occupied = True + coordinates[atoms.index(atom),:] = coordinates[index1,:] + vector + + # Recursively continue with functional group + self.__generateFunctionalGroupCoordinates(atom1, atom) def __generateAtomLabels(self): """ diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 49c0e6fc8a..183fd6b8f6 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1157,6 +1157,11 @@ def toSMILES(self): `RDKit `_ to perform the conversion. Perceives aromaticity and removes Hydrogen atoms. """ + + for atom in self.vertices: + if atom.isNitrogen(): + return self.getFormula() + rdkitmol = self.toRDKitMol() return Chem.MolToSmiles(rdkitmol) From 9d790015c71be97e2937d4f3079e7f4a3f20233b Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Thu, 24 Oct 2013 16:58:35 -0400 Subject: [PATCH 24/64] Additional Atom Energies and Bond Additivity Corrections --- rmgpy/cantherm/statmech.py | 138 ++++++++++++++++++++----------------- 1 file changed, 73 insertions(+), 65 deletions(-) diff --git a/rmgpy/cantherm/statmech.py b/rmgpy/cantherm/statmech.py index 4acd322f36..8d57110956 100644 --- a/rmgpy/cantherm/statmech.py +++ b/rmgpy/cantherm/statmech.py @@ -42,7 +42,6 @@ import rmgpy.constants as constants from rmgpy.cantherm.output import prettify from rmgpy.cantherm.gaussian import GaussianLog -from rmgpy.cantherm.molepro import MoleProLog from rmgpy.species import TransitionState from rmgpy.statmech import * @@ -196,7 +195,6 @@ def load(self): 'HinderedRotor': hinderedRotor, # File formats 'GaussianLog': GaussianLog, - 'MoleProLog': MoleProLog, 'ScanLog': ScanLog, } @@ -249,10 +247,7 @@ def load(self): except KeyError: raise InputError('Model chemistry {0!r} not found in from dictionary of energy values in species file {1!r}.'.format(self.modelChemistry, path)) if isinstance(energy, GaussianLog): - energyLog = energy; E0 = 'Gaussian' - energyLog.path = os.path.join(directory, energyLog.path) - if isinstance(energy, MoleProLog): - energyLog = energy; E0 = 'MolePro' + energyLog = energy; E0 = None energyLog.path = os.path.join(directory, energyLog.path) elif isinstance(energy, float): energyLog = None; E0 = energy @@ -292,10 +287,8 @@ def load(self): logging.debug(' Reading energy...') # The E0 that is read from the log file is without the ZPE and corresponds to E_elec - if E0 is 'Gaussian': + if E0 is None: E0 = energyLog.loadEnergy(self.frequencyScaleFactor) - elif E0 is 'MolePro': - E0 = energyLog.loadCCSDEnergy() else: E0 = E0 * constants.E_h * constants.Na # Hartree/particle to J/mol E0 = applyEnergyCorrections(E0, self.modelChemistry, atoms, bonds if self.applyBondEnergyCorrections else {}) @@ -397,7 +390,7 @@ def save(self, outputFile): logging.info('Saving statistical mechanics parameters for {0}...'.format(self.species.label)) f = open(outputFile, 'a') - numbers = {1: 'H', 6: 'C', 7: 'N', 8: 'O', 14: 'Si', 15: 'P', 16: 'S', 17: 'Cl'} + numbers = {1: 'H', 6: 'C', 7: 'N', 8: 'O', 14: 'Si', 15: 'P', 16: 'S'} conformer = self.species.conformer @@ -480,106 +473,114 @@ def applyEnergyCorrections(E0, modelChemistry, atoms, bonds): # Spin orbit correction (SOC) in Hartrees # Values taken from note 22 of http://jcp.aip.org/resource/1/jcpsa6/v109/i24/p10570_s1 and converted to hartrees # Values in millihartree are also available (with fewer significant figures) from http://jcp.aip.org/resource/1/jcpsa6/v106/i3/p1063_s1 - SOC = {'H':0.0, 'N':0.0, 'O': -0.000355, 'C': -0.000135, 'S': -0.000893, 'P': 0.0, 'Cl': -0.001338} + SOC = {'H':0.0, 'N':0.0, 'O': -0.000355, 'C': -0.000135, 'S': -0.000893, 'P': 0.0} # Step 1: Reference all energies to a model chemistry-independent basis # by subtracting out that model chemistry's atomic energies # Note: If your model chemistry does not include spin orbit coupling, you should add the corrections to the energies here if modelChemistry == 'CBS-QB3': - # 0K Energy - atomEnergies = {'H':-0.499818 , 'N':-54.520543, 'O':-74.987624, 'C':-37.785385, 'P':-340.817186, 'S': -397.657360, 'Cl': -459.683605} + atomEnergies = {'H':-0.499818 , 'N':-54.520543, 'O':-74.987624, 'C':-37.785385, 'P':-340.817186, 'S': -397.657360} elif modelChemistry == 'G3': atomEnergies = {'H':-0.5010030, 'N':-54.564343, 'O':-75.030991, 'C':-37.827717, 'P':-341.116432, 'S': -397.961110} elif modelChemistry == 'Klip_1': - atomEnergies = {'H':-0.50003976 + SOC['H'], 'N':-54.53383153 + SOC['N'], 'O':-75.00935474 + SOC['O'], 'C':-37.79266591 + SOC['C']} + atomEnergies = {'H':-0.50003976, 'N':-54.53383153, 'O':-75.00935474, 'C':-37.79266591} elif modelChemistry == 'Klip_2': #Klip QCI(tz,qz) - atomEnergies = {'H':-0.50003976 + SOC['H'], 'N':-54.53169400 + SOC['N'], 'O':-75.00714902 + SOC['O'], 'C':-37.79060419 + SOC['C']} + atomEnergies = {'H':-0.50003976, 'N':-54.53169400, 'O':-75.00714902, 'C':-37.79060419} elif modelChemistry == 'Klip_3': #Klip QCI(dz,tz) - atomEnergies = {'H':-0.50005578 + SOC['H'], 'N':-54.53128140 + SOC['N'], 'O':-75.00356581 + SOC['O'], 'C':-37.79025175 + SOC['C']} + atomEnergies = {'H':-0.50005578, 'N':-54.53128140, 'O':-75.00356581, 'C':-37.79025175} elif modelChemistry == 'Klip_2_cc': #Klip CCSD(T)(tz,qz) - atomEnergies = {'H':-0.50003976 + SOC['H'], 'O':-75.00681155 + SOC['O'], 'C':-37.79029443 + SOC['C']} + atomEnergies = {'H':-0.50003976, 'O':-75.00681155, 'C':-37.79029443} elif modelChemistry == 'CCSD(T)-F12/cc-pVDZ-F12_H-TZ': - atomEnergies = {'H':-0.499946213243 + SOC['H'], 'N':-54.526406291655 + SOC['N'], 'O':-74.995458316117 + SOC['O'], 'C':-37.788203485235 + SOC['C']} + atomEnergies = {'H':-0.499946213243, 'N':-54.526406291655, 'O':-74.995458316117, 'C':-37.788203485235} elif modelChemistry == 'CCSD(T)-F12/cc-pVDZ-F12_H-QZ': - atomEnergies = {'H':-0.499994558325 + SOC['H'], 'N':-54.526406291655 + SOC['N'], 'O':-74.995458316117 + SOC['O'], 'C':-37.788203485235 + SOC['C']} + atomEnergies = {'H':-0.499994558325, 'N':-54.526406291655, 'O':-74.995458316117, 'C':-37.788203485235} elif modelChemistry == 'CCSD(T)-F12/cc-pVDZ-F12': - atomEnergies = {'H':-0.499811124128 + SOC['H'], 'N':-54.526406291655 + SOC['N'], 'O':-74.995458316117 + SOC['O'], 'C':-37.788203485235 + SOC['C']} +# atomEnergies = {'H':-0.499811124128, 'N':-54.526406291655, 'O':-74.995458316117, 'C':-37.788203485235} + atomEnergies = {'H':-0.499811124128, 'N':-54.526406291655, 'O':-74.995458316117, 'C':-37.788203485235} elif modelChemistry == 'CCSD(T)-F12/cc-pVTZ-F12': - atomEnergies = {'H':-0.499946213243 + SOC['H'], 'N':-54.53000909621 + SOC['N'], 'O':-75.004127673424 + SOC['O'], 'C':-37.789862146471 + SOC['C']} + atomEnergies = {'H':-0.499946213243, 'N':-54.53000909621, 'O':-75.004127673424, 'C':-37.789862146471} elif modelChemistry == 'CCSD(T)-F12/cc-pVQZ-F12': - atomEnergies = {'H':-0.499994558325 + SOC['H'], 'N':-54.530515226371 + SOC['N'], 'O':-75.005600062003 + SOC['O'], 'C':-37.789961656228 + SOC['C']} + atomEnergies = {'H':-0.499994558325, 'N':-54.530515226371, 'O':-75.005600062003, 'C':-37.789961656228} elif modelChemistry == 'CCSD(T)-F12/cc-pCVDZ-F12': - atomEnergies = {'H':-0.499811124128 + SOC['H'], 'N':-54.582137180344 + SOC['N'], 'O':-75.053045547421 + SOC['O'], 'C':-37.840869118707 + SOC['C']} + atomEnergies = {'H':-0.499811124128, 'N':-54.582137180344, 'O':-75.053045547421, 'C':-37.840869118707} elif modelChemistry == 'CCSD(T)-F12/cc-pCVTZ-F12': - atomEnergies = {'H':-0.499946213243 + SOC['H'], 'N':-54.588545831900 + SOC['N'], 'O':-75.065995072347 + SOC['O'], 'C':-37.844662139972 + SOC['C']} + atomEnergies = {'H':-0.499946213243, 'N':-54.588545831900, 'O':-75.065995072347, 'C':-37.844662139972} elif modelChemistry == 'CCSD(T)-F12/cc-pCVQZ-F12': - atomEnergies = {'H':-0.499994558325 + SOC['H'], 'N':-54.589137594139 + SOC['N'], 'O':-75.067412234737 + SOC['O'], 'C':-37.844893820561 + SOC['C']} + atomEnergies = {'H':-0.499994558325, 'N':-54.589137594139, 'O':-75.067412234737, 'C':-37.844893820561} elif modelChemistry == 'CCSD(T)-F12/aug-cc-pVDZ': - atomEnergies = {'H':-0.499459066131 + SOC['H'], 'N':-54.524279516472 + SOC['N'], 'O':-74.992097308083 + SOC['O'], 'C':-37.786694171716 + SOC['C']} + atomEnergies = {'H':-0.499459066131, 'N':-54.524279516472, 'O':-74.992097308083, 'C':-37.786694171716} elif modelChemistry == 'CCSD(T)-F12/aug-cc-pVTZ': - atomEnergies = {'H':-0.499844820798 + SOC['H'], 'N':-54.527419359906 + SOC['N'], 'O':-75.000001429806 + SOC['O'], 'C':-37.788504810868 + SOC['C']} + atomEnergies = {'H':-0.499844820798, 'N':-54.527419359906, 'O':-75.000001429806, 'C':-37.788504810868} elif modelChemistry == 'CCSD(T)-F12/aug-cc-pVQZ': - atomEnergies = {'H':-0.499949526073 + SOC['H'], 'N':-54.529569719016 + SOC['N'], 'O':-75.004026586610 + SOC['O'], 'C':-37.789387892348 + SOC['C']} + atomEnergies = {'H':-0.499949526073, 'N':-54.529569719016, 'O':-75.004026586610, 'C':-37.789387892348} elif modelChemistry == 'B-CCSD(T)-F12/cc-pVDZ-F12': - atomEnergies = {'H':-0.499811124128 + SOC['H'], 'N':-54.523269942190 + SOC['N'], 'O':-74.990725918500 + SOC['O'], 'C':-37.785409916465 + SOC['C']} + atomEnergies = {'H':-0.499811124128, 'N':-54.523269942190, 'O':-74.990725918500, 'C':-37.785409916465} elif modelChemistry == 'B-CCSD(T)-F12/cc-pVTZ-F12': - atomEnergies = {'H':-0.499946213243 + SOC['H'], 'N':-54.528135889213 + SOC['N'], 'O':-75.001094055506 + SOC['O'], 'C':-37.788233578503 + SOC['C']} + atomEnergies = {'H':-0.499946213243, 'N':-54.528135889213, 'O':-75.001094055506, 'C':-37.788233578503} elif modelChemistry == 'B-CCSD(T)-F12/cc-pVQZ-F12': - atomEnergies = {'H':-0.499994558325 + SOC['H'], 'N':-54.529425753163 + SOC['N'], 'O':-75.003820485005 + SOC['O'], 'C':-37.789006506290 + SOC['C']} + atomEnergies = {'H':-0.499994558325, 'N':-54.529425753163, 'O':-75.003820485005, 'C':-37.789006506290} elif modelChemistry == 'B-CCSD(T)-F12/cc-pCVDZ-F12': - atomEnergies = {'H':-0.499811124128 + SOC['H'], 'N':-54.578602780288 + SOC['N'], 'O':-75.048064317367 + SOC['O'], 'C':-37.837592033417 + SOC['C']} + atomEnergies = {'H':-0.499811124128, 'N':-54.578602780288, 'O':-75.048064317367, 'C':-37.837592033417} elif modelChemistry == 'B-CCSD(T)-F12/cc-pCVTZ-F12': - atomEnergies = {'H':-0.499946213243 + SOC['H'], 'N':-54.586402551258 + SOC['N'], 'O':-75.062767632757 + SOC['O'], 'C':-37.842729156944 + SOC['C']} + atomEnergies = {'H':-0.499946213243, 'N':-54.586402551258, 'O':-75.062767632757, 'C':-37.842729156944} elif modelChemistry == 'B-CCSD(T)-F12/cc-pCVQZ-F12': - atomEnergies = {'H':-0.49999456 + SOC['H'], 'N':-54.587781507581 + SOC['N'], 'O':-75.065397706471 + SOC['O'], 'C':-37.843634971592 + SOC['C']} + atomEnergies = {'H':-0.49999456, 'N':-54.587781507581, 'O':-75.065397706471, 'C':-37.843634971592} elif modelChemistry == 'B-CCSD(T)-F12/aug-cc-pVDZ': - atomEnergies = {'H':-0.499459066131 + SOC['H'], 'N':-54.520475581942 + SOC['N'], 'O':-74.986992215049 + SOC['O'], 'C':-37.783294495799 + SOC['C']} + atomEnergies = {'H':-0.499459066131, 'N':-54.520475581942, 'O':-74.986992215049, 'C':-37.783294495799} elif modelChemistry == 'B-CCSD(T)-F12/aug-cc-pVTZ': - atomEnergies = {'H':-0.499844820798 + SOC['H'], 'N':-54.524927371700 + SOC['N'], 'O':-74.996328829705 + SOC['O'], 'C':-37.786320700792 + SOC['C']} + atomEnergies = {'H':-0.499844820798, 'N':-54.524927371700, 'O':-74.996328829705, 'C':-37.786320700792} elif modelChemistry == 'B-CCSD(T)-F12/aug-cc-pVQZ': - atomEnergies = {'H':-0.499949526073 + SOC['H'], 'N':-54.528189769291 + SOC['N'], 'O':-75.001879610563 + SOC['O'], 'C':-37.788165047059 + SOC['C']} - + atomEnergies = {'H':-0.499949526073, 'N':-54.528189769291, 'O':-75.001879610563, 'C':-37.788165047059} + elif modelChemistry == 'DFT_G03_b3lyp': + atomEnergies = {'H':-0.502256981529, 'N':-54.6007233648, 'O':-75.0898777574, 'C':-37.8572666349} elif modelChemistry == 'DFT_ks_b3lyp': - atomEnergies = {'H':-0.49785866 + SOC['H'], 'N':-54.45608798 + SOC['N'], 'O':-74.93566254 + SOC['O'], 'C':-37.76119132 + SOC['C']} + atomEnergies = {'H':-0.49785866, 'N':-54.45608798, 'O':-74.93566254, 'C':-37.76119132} elif modelChemistry == 'DFT_uks_b3lyp': - atomEnergies = {'H':-0.49785866 + SOC['H'], 'N':-54.45729113 + SOC['N'], 'O':-74.93566254 + SOC['O'], 'C':-37.76119132 + SOC['C']} + atomEnergies = {'H':-0.49785866, 'N':-54.45729113, 'O':-74.93566254, 'C':-37.76119132} elif modelChemistry == 'MP2_rmp2_pVDZ': - atomEnergies = {'H':-0.49927840 + SOC['H'], 'N':-54.46141996 + SOC['N'], 'O':-74.89408254 + SOC['O'], 'C':-37.73792713 + SOC['C']} + atomEnergies = {'H':-0.49927840, 'N':-54.46141996, 'O':-74.89408254, 'C':-37.73792713} elif modelChemistry == 'MP2_rmp2_pVTZ': - atomEnergies = {'H':-0.49980981 + SOC['H'], 'N':-54.49615972 + SOC['N'], 'O':-74.95506980 + SOC['O'], 'C':-37.75833104 + SOC['C']} + atomEnergies = {'H':-0.49980981, 'N':-54.49615972, 'O':-74.95506980, 'C':-37.75833104} elif modelChemistry == 'MP2_rmp2_pVQZ': - atomEnergies = {'H':-0.49994557 + SOC['H'], 'N':-54.50715868 + SOC['N'], 'O':-74.97515364 + SOC['O'], 'C':-37.76533215 + SOC['C']} + atomEnergies = {'H':-0.49994557, 'N':-54.50715868, 'O':-74.97515364, 'C':-37.76533215} + + elif modelChemistry == 'CCSD-F12/cc-pVDZ-F12': + atomEnergies = {'H':-0.499811124128, 'N':-54.524325513811, 'O':-74.992326577897, 'C':-37.786213495943} + + elif modelChemistry == 'CCSD(T)-F12/cc-pVDZ-F12_noscale': + atomEnergies = {'H':-0.499811124128, 'N':-54.526026290887, 'O':-74.994751897699, 'C':-37.787881871511} + + elif modelChemistry == 'G03_PBEPBE_6-311++g_d_p': + atomEnergies = {'H':-0.499812273282, 'N':-54.5289567564, 'O':-75.0033596764, 'C':-37.7937388736} + + elif modelChemistry == 'FCI/cc-pVDZ': +# atomEnergies = {'C':-37.760717371923} + atomEnergies = {'C':-37.789527} + elif modelChemistry == 'FCI/cc-pVTZ': + atomEnergies = {'C':-37.781266669684} + elif modelChemistry == 'FCI/cc-pVQZ': + atomEnergies = {'C':-37.787052110598} - elif modelChemistry == 'CCSD_DZ': - atomEnergies = {'H':-0.499811124 + SOC['H'], 'N':-54.52640629 + SOC['N'], 'O':-74.99545832 + SOC['O'], 'C':-37.78820349 + SOC['C']} - elif modelChemistry == 'CCSD_TZ': - atomEnergies = {'H':-0.499946213 + SOC['H'], 'N':-54.5300091 + SOC['N'], 'O':-75.00412767 + SOC['O'], 'C':-37.78986215 + SOC['C']} - elif modelChemistry == 'CCSD_QZ': - atomEnergies = {'H':-0.499994558 + SOC['H'], 'N':-54.53051523 + SOC['N'], 'O':-75.00560006 + SOC['O'], 'C':-37.78996166 + SOC['C']} - elif modelChemistry == 'CCSD_core_DZ': - atomEnergies = {'H':-0.499811124 + SOC['H'], 'N':-54.58213718 + SOC['N'], 'O':-75.05304555 + SOC['O'], 'C':-37.84086912 + SOC['C']} - - else: logging.warning('Unknown model chemistry "{0}"; not applying energy corrections.'.format(modelChemistry)) return E0 for symbol, count in atoms.items(): - if symbol in atomEnergies: E0 -= count * atomEnergies[symbol] * constants.E_h * constants.Na + if symbol in atomEnergies: E0 -= count * atomEnergies[symbol] * 4.35974394e-18 * constants.Na else: logging.warning('Ignored unknown atom type "{0}".'.format(symbol)) @@ -588,10 +589,10 @@ def applyEnergyCorrections(E0, modelChemistry, atoms, bonds): # See Gaussian thermo whitepaper at http://www.gaussian.com/g_whitepap/thermo.htm) # Note: these values are relatively old and some improvement may be possible by using newer values, particularly for carbon # However, care should be taken to ensure that they are compatible with the BAC values (if BACs are used) - atomHf = {'H': 51.63 , 'N': 112.53 ,'O': 58.99 ,'C': 169.98, 'S': 65.66, 'Cl': 28.59 } + atomHf = {'H': 51.63 , 'N': 112.53 ,'O': 58.99 ,'C': 169.98, 'S': 65.55 } # Thermal contribution to enthalpy Hss(298 K) - Hss(0 K) reported by Gaussian thermo whitepaper # This will be subtracted from the corresponding value in atomHf to produce an enthalpy used in calculating the enthalpy of formation at 298 K - atomThermal = {'H': 1.01 , 'N': 1.04, 'O': 1.04 ,'C': 0.25, 'S': 1.05, 'Cl': 1.10 } + atomThermal = {'H': 1.01 , 'N': 1.04, 'O': 1.04 ,'C': 0.25, 'S': 1.05 } # Total energy correction used to reach gas-phase reference state # Note: Spin orbit coupling no longer included in these energies, since some model chemistries include it automatically atomEnergies = {} @@ -600,19 +601,26 @@ def applyEnergyCorrections(E0, modelChemistry, atoms, bonds): for symbol, count in atoms.items(): if symbol in atomEnergies: E0 += count * atomEnergies[symbol] * 4184. - # Step 3: Bond additivity corrections + # Step 3: Bond energy corrections if modelChemistry == 'CCSD(T)-F12/cc-pVDZ-F12': - bondEnergies = { 'C-H': -0.56, 'C-C': -0.53, 'C=C': -1.90, 'C#C': -0.64, - 'O-H': -0.34, 'C-O': -0.30, 'C=O': -0.92, 'O-O': 0.03, 'N-C': -0.49, - 'N=C': -1.50, 'N#C': -3.54, 'N-O': 0.60, 'N_O': -0.17, 'N=O': -0.72, - 'N-H': -0.75, 'N-N': -1.45, 'N=N': -1.98, 'N#N': -2.05,} + bondEnergies = { 'C-H': -0.46, 'C-C': -0.68, 'C=C': -1.90, 'C#C': -3.13, + 'O-H': -0.51, 'C-O': -0.23, 'C=O': -0.69, 'O-O': -0.02, 'N-C': -0.67, + 'N=C': -1.46, 'N#C': -2.79, 'N-O': 0.74, 'N_O': -0.23, 'N=O': -0.51, + 'N-H': -0.69, 'N-N': -0.47, 'N=N': -1.54, 'N#N': -2.05,} + elif modelChemistry == 'CCSD(T)-F12/cc-pVTZ-F12': + bondEnergies = { 'C-H': -0.09, 'C-C': -0.27, 'C=C': -1.03, 'C#C': -1.79, + 'O-H': -0.06, 'C-O': -0.14, 'C=O': -0.19, 'O-O': 0.16, 'N-C': -0.18, + 'N=C': -0.41, 'N#C': -1.41, 'N-O': 0.87, 'N_O': -0.09, 'N=O': -0.23, + 'N-H': -0.01, 'N-N': -0.21, 'N=N': -0.44, 'N#N': -0.76,} + elif modelChemistry == 'CCSD(T)-F12/cc-pVQZ-F12': + bondEnergies = { 'C-H': -0.08, 'C-C': -0.26, 'C=C': -1.01, 'C#C': -1.66, + 'O-H': 0.07, 'C-O': 0.25, 'C=O': -0.03, 'O-O': 0.26, 'N-C': -0.20, + 'N=C': -0.30, 'N#C': -1.33, 'N-O': 1.01, 'N_O': -0.03, 'N=O': -0.26, + 'N-H': 0.06, 'N-N': -0.23, 'N=N': -0.37, 'N#N': -0.64,} else: - # BAC corrections from Table IX in http://jcp.aip.org/resource/1/jcpsa6/v109/i24/p10570_s1 for CBS-Q method - # H-Cl correction from CBS-QB3 enthalpy difference with Gurvich 1989, HF298=-92.31 kJ bondEnergies = { 'C-H': -0.11, 'C-C': -0.3, 'C=C': -0.08, 'C#C': -0.64, 'O-H': 0.02, 'C-O': 0.33, 'C=O': 0.55, 'N#N': -2.0, 'O=O': -0.2, - 'H-H': 1.1, 'C#N': -0.89, 'C-S': 0.43, 'S=O': -0.78, 'C-Cl': 1.29, - 'N-H': -0.42, 'C-N': -0.13, 'S-H': 0.00, 'H-Cl': 1.16 } + 'H-H': 1.1, 'C#N': -0.89, 'C-S': 0.43, 'S=O': -0.78 } for symbol, count in bonds.items(): if symbol in bondEnergies: E0 += count * bondEnergies[symbol] * 4184. From 176c21606b17850afe288936daaed6fdd8ec562b Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Thu, 24 Oct 2013 19:27:05 -0400 Subject: [PATCH 25/64] Increase atomElectronStateWidth by +1 --- rmgpy/molecule/adjlist.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index 793322db21..fe555abbe0 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -341,7 +341,7 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False): atomLabelWidth = max([len(s) for s in atomLabels.values()]) if atomLabelWidth > 0: atomLabelWidth += 1 atomTypeWidth = max([len(s) for s in atomTypes.values()]) + 1 - atomElectronStateWidth = max([len(s) for s in atomElectronStates.values()]) + atomElectronStateWidth = max([len(s) for s in atomElectronStates.values()]) + 1 atomLonePairWidth = 1 # Assemble the adjacency list From 906d5046e7c8610626e67cad7a40071c7c021e20 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Sat, 26 Oct 2013 13:19:02 -0400 Subject: [PATCH 26/64] Use OpenBabel for InChKey generation --- rmgpy/molecule/molecule.py | 53 +++++++++++++++++++++++++++++++++++++- 1 file changed, 52 insertions(+), 1 deletion(-) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 183fd6b8f6..a6a6012f91 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -40,7 +40,7 @@ import os import re import element as elements - +import openbabel from rdkit import Chem from .graph import Vertex, Edge, Graph from .group import GroupAtom, GroupBond, Group, ActionError @@ -1116,12 +1116,31 @@ def toAugmentedInChI(self): def toInChIKey(self): """ + Convert a molecular structure to an InChI Key string. Uses + `OpenBabel `_ to perform the conversion. + + or + Convert a molecular structure to an InChI Key string. Uses `RDKit `_ to perform the conversion. Removes check-sum dash (-) and character so that only the 14 + 9 characters remain. """ + + import openbabel + + + for atom in self.vertices: + if atom.isNitrogen(): + # This version does not write a warning to stderr if stereochemistry is undefined + obmol = self.toOBMol() + obConversion = openbabel.OBConversion() + obConversion.SetOutFormat('inchi') + obConversion.SetOptions('w', openbabel.OBConversion.OUTOPTIONS) + obConversion.SetOptions('K', openbabel.OBConversion.OUTOPTIONS) + return obConversion.WriteString(obmol).strip()[:-2] + if not Chem.inchi.INCHI_AVAILABLE: return "RDKitInstalledWithoutInChI" inchi = self.toInChI() @@ -1165,6 +1184,38 @@ def toSMILES(self): rdkitmol = self.toRDKitMol() return Chem.MolToSmiles(rdkitmol) + + def toOBMol(self): + """ + Convert a molecular structure to an OpenBabel OBMol object. Uses + `OpenBabel `_ to perform the conversion. + """ + cython.declare(atom=Atom, atom1=Atom, bonds=dict, atom2=Atom, bond=Bond) + cython.declare(index1=cython.int, index2=cython.int, order=cython.int) + + # Sort the atoms before converting to ensure output is consistent + # between different runs + self.sortAtoms() + + atoms = self.vertices + + obmol = openbabel.OBMol() + for atom in atoms: + a = obmol.NewAtom() + a.SetAtomicNum(atom.number) + a.SetFormalCharge(atom.charge) + orders = {'S': 1, 'D': 2, 'T': 3, 'B': 5} + for atom1 in self.vertices: + for atom2, bond in atom1.edges.iteritems(): + index1 = atoms.index(atom1) + index2 = atoms.index(atom2) + if index1 < index2: + order = orders[bond.order] + obmol.AddBond(index1+1, index2+1, order) + + obmol.AssignSpinMultiplicity(True) + + return obmol def toRDKitMol(self, removeHs=True, returnMapping=False): """ From 880460c715d7dd3e2156fa9eeeea1c3d00ef5141 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Sat, 26 Oct 2013 13:20:46 -0400 Subject: [PATCH 27/64] Do not remove hydrogen from thermoHBIcheck.txt --- rmgpy/rmg/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/rmg/model.py b/rmgpy/rmg/model.py index b8ea2b53a5..395aad0671 100644 --- a/rmgpy/rmg/model.py +++ b/rmgpy/rmg/model.py @@ -138,7 +138,7 @@ def generateThermoData(self, database, thermoClass=NASA, quantumMechanics=None): with open('thermoHBIcheck.txt','a') as f: f.write('// {0!r}\n'.format(thermo0).replace('),','),\n// ')) f.write('{0}\n'.format(molecule.toSMILES())) - f.write('{0}\n\n'.format(molecule.toAdjacencyList(removeH=True))) + f.write('{0}\n\n'.format(molecule.toAdjacencyList(removeH=False))) else: # Not too many radicals: do a direct calculation. thermo0 = quantumMechanics.getThermoData(molecule) # returns None if it fails From b67b2c232f6939f0f64e41ebc3911c30469c1836 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Tue, 29 Oct 2013 11:01:33 -0400 Subject: [PATCH 28/64] Exclude temporarily tetra-valent nitrogen from QM because RDKit does not know it --- rmgpy/qm/mopac.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/rmgpy/qm/mopac.py b/rmgpy/qm/mopac.py index 04ded33b91..d5d38a2b2c 100644 --- a/rmgpy/qm/mopac.py +++ b/rmgpy/qm/mopac.py @@ -197,6 +197,10 @@ def generateQMData(self): """ Calculate the QM data and return a QMData object, or None if it fails. """ + for atom in self.molecule.vertices: + if atom.atomType.label == 'N4s' or atom.atomType.label == 'N4d' or atom.atomType.label =='N4dd' or atom.atomType.label == 'N4t' or atom.atomType.label == 'N4b': + return None + if self.verifyOutputFile(): logging.info("Found a successful output file already; using that.") source = "QM MOPAC result file found from previous run." From 2678191334a516ab36397dca537e5f38856182a4 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 4 Nov 2013 10:40:29 -0500 Subject: [PATCH 29/64] Change N4 to N5 --- rmgpy/molecule/atomtype.py | 48 +++++++++++++++++++------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/rmgpy/molecule/atomtype.py b/rmgpy/molecule/atomtype.py index 781acc2785..db7f006105 100644 --- a/rmgpy/molecule/atomtype.py +++ b/rmgpy/molecule/atomtype.py @@ -162,14 +162,14 @@ def isSpecificCaseOf(self, other): 'R!H', 'H', 'C','Cs','Cd','Cdd','Ct','CO','Cb','Cbf', 'CS', - 'N','N2d','N3s','N3d','N3t','N3b','N4s','N4d','N4dd','N4t','N4b', + 'N','N2d','N3s','N3d','N3t','N3b','N5s','N5d','N5dd','N5t','N5b', 'O','Os','Od','Oa', 'Si','Sis','Sid','Sidd','Sit','SiO','Sib','Sibf', 'S','Ss','Sd','Sa'] ) atomTypes['R!H'] = AtomType(label='R!H', generic=['R'], specific=[ 'C','Cs','Cd','Cdd','Ct','CO','Cb','Cbf', 'CS', - 'N','N2d','N3s','N3d','N3t','N3b','N4s','N4d','N4dd','N4t','N4b', + 'N','N2d','N3s','N3d','N3t','N3b','N5s','N5d','N5dd','N5t','N5b', 'O','Os','Od','Oa', 'Si','Sis','Sid','Sidd','Sit','SiO','Sib','Sibf', 'S','Ss','Sd','Sa'] @@ -186,17 +186,17 @@ def isSpecificCaseOf(self, other): atomTypes['Cbf' ] = AtomType('Cbf', generic=['R','R!H','C'], specific=[]) atomTypes['CS' ] = AtomType('CS', generic=['R','R!H','C'], specific=[]) -atomTypes['N' ] = AtomType('N', generic=['R','R!H'], specific=['N2d','N3s','N3d','N3t','N3b','N4s','N4d','N4dd','N4t','N4b']) +atomTypes['N' ] = AtomType('N', generic=['R','R!H'], specific=['N2d','N3s','N3d','N3t','N3b','N5s','N5d','N5dd','N5t','N5b']) atomTypes['N2d' ] = AtomType('N2d', generic=['R','R!H','N'], specific=[]) atomTypes['N3s' ] = AtomType('N3s', generic=['R','R!H','N'], specific=[]) atomTypes['N3d' ] = AtomType('N3d', generic=['R','R!H','N'], specific=[]) atomTypes['N3t' ] = AtomType('N3t', generic=['R','R!H','N'], specific=[]) atomTypes['N3b' ] = AtomType('N3b', generic=['R','R!H','N'], specific=[]) -atomTypes['N4s' ] = AtomType('N4s', generic=['R','R!H','N'], specific=[]) -atomTypes['N4d' ] = AtomType('N4d', generic=['R','R!H','N'], specific=[]) -atomTypes['N4dd'] = AtomType('N4dd', generic=['R','R!H','N'], specific=[]) -atomTypes['N4t' ] = AtomType('N4t', generic=['R','R!H','N'], specific=[]) -atomTypes['N4b' ] = AtomType('N4b', generic=['R','R!H','N'], specific=[]) +atomTypes['N5s' ] = AtomType('N5s', generic=['R','R!H','N'], specific=[]) +atomTypes['N5d' ] = AtomType('N5d', generic=['R','R!H','N'], specific=[]) +atomTypes['N5dd'] = AtomType('N5dd', generic=['R','R!H','N'], specific=[]) +atomTypes['N5t' ] = AtomType('N5t', generic=['R','R!H','N'], specific=[]) +atomTypes['N5b' ] = AtomType('N5b', generic=['R','R!H','N'], specific=[]) atomTypes['O' ] = AtomType('O', generic=['R','R!H'], specific=['Os','Od','Oa']) atomTypes['Os' ] = AtomType('Os', generic=['R','R!H','O'], specific=[]) @@ -234,15 +234,15 @@ def isSpecificCaseOf(self, other): atomTypes['N' ].setActions(incrementBond=['N'], decrementBond=['N'], formBond=['N'], breakBond=['N'], incrementRadical=['N'], decrementRadical=['N'], incrementLonePair=['N'], decrementLonePair=['N']) atomTypes['N2d' ].setActions(incrementBond=[], decrementBond=[], formBond=[], breakBond=[], incrementRadical=[], decrementRadical=[], incrementLonePair=[], decrementLonePair=[]) -atomTypes['N3s' ].setActions(incrementBond=['N3d','N3s'], decrementBond=[], formBond=['N3s','N4s'], breakBond=['N3s'], incrementRadical=['N3s'], decrementRadical=['N3s'], incrementLonePair=['N3s'], decrementLonePair=['N3s']) -atomTypes['N3d' ].setActions(incrementBond=['N3t'], decrementBond=['N3s'], formBond=['N3d','N4d'], breakBond=['N3d'], incrementRadical=['N3d'], decrementRadical=['N3d'], incrementLonePair=['N3d'], decrementLonePair=['N3d']) -atomTypes['N3t' ].setActions(incrementBond=[], decrementBond=['N3d'], formBond=['N4t'], breakBond=[], incrementRadical=['N3t'], decrementRadical=['N3t'], incrementLonePair=['N3t'], decrementLonePair=['N3t']) +atomTypes['N3s' ].setActions(incrementBond=['N3d','N3s'], decrementBond=[], formBond=['N3s','N5s'], breakBond=['N3s'], incrementRadical=['N3s'], decrementRadical=['N3s'], incrementLonePair=['N3s'], decrementLonePair=['N3s']) +atomTypes['N3d' ].setActions(incrementBond=['N3t'], decrementBond=['N3s'], formBond=['N3d','N5d'], breakBond=['N3d'], incrementRadical=['N3d'], decrementRadical=['N3d'], incrementLonePair=['N3d'], decrementLonePair=['N3d']) +atomTypes['N3t' ].setActions(incrementBond=[], decrementBond=['N3d'], formBond=['N5t'], breakBond=[], incrementRadical=['N3t'], decrementRadical=['N3t'], incrementLonePair=['N3t'], decrementLonePair=['N3t']) atomTypes['N3b' ].setActions(incrementBond=[], decrementBond=[], formBond=['N3b'], breakBond=['N3b'], incrementRadical=['N3b'], decrementRadical=['N3b'], incrementLonePair=['N3b'], decrementLonePair=['N3b']) -atomTypes['N4s' ].setActions(incrementBond=['N4d'], decrementBond=[], formBond=[], breakBond=['N3s'], incrementRadical=['N4s'], decrementRadical=['N4s'], incrementLonePair=['N4s'], decrementLonePair=['N4s']) -atomTypes['N4d' ].setActions(incrementBond=['N4dd','N4t'], decrementBond=['N4s'], formBond=[], breakBond=['N3d'], incrementRadical=['N4d'], decrementRadical=['N4d'], incrementLonePair=['N4d'], decrementLonePair=['N4d']) -atomTypes['N4dd'].setActions(incrementBond=[], decrementBond=['N3d'], formBond=[], breakBond=[], incrementRadical=['N4dd'], decrementRadical=['N4dd'], incrementLonePair=['N4d'], decrementLonePair=['N4d']) -atomTypes['N4t' ].setActions(incrementBond=[], decrementBond=['N3d','N3t'], formBond=[], breakBond=['N3d','N3t'], incrementRadical=['N4t'], decrementRadical=['N4t'], incrementLonePair=['N4t'], decrementLonePair=['N4t']) -atomTypes['N4b' ].setActions(incrementBond=[], decrementBond=[], formBond=['N4b'], breakBond=['N4b'], incrementRadical=['N4b'], decrementRadical=['N4b'], incrementLonePair=['N4b'], decrementLonePair=['N4b']) +atomTypes['N5s' ].setActions(incrementBond=['N5d'], decrementBond=[], formBond=[], breakBond=['N3s'], incrementRadical=['N5s'], decrementRadical=['N5s'], incrementLonePair=['N5s'], decrementLonePair=['N5s']) +atomTypes['N5d' ].setActions(incrementBond=['N5dd','N5t'], decrementBond=['N5s'], formBond=[], breakBond=['N3d'], incrementRadical=['N5d'], decrementRadical=['N5d'], incrementLonePair=['N5d'], decrementLonePair=['N5d']) +atomTypes['N5dd'].setActions(incrementBond=[], decrementBond=['N3d'], formBond=[], breakBond=[], incrementRadical=['N5dd'], decrementRadical=['N5dd'], incrementLonePair=['N5d'], decrementLonePair=['N5d']) +atomTypes['N5t' ].setActions(incrementBond=[], decrementBond=['N3d','N3t'], formBond=[], breakBond=['N3d','N3t'], incrementRadical=['N5t'], decrementRadical=['N5t'], incrementLonePair=['N5t'], decrementLonePair=['N5t']) +atomTypes['N5b' ].setActions(incrementBond=[], decrementBond=[], formBond=['N5b'], breakBond=['N5b'], incrementRadical=['N5b'], decrementRadical=['N5b'], incrementLonePair=['N5b'], decrementLonePair=['N5b']) atomTypes['O' ].setActions(incrementBond=['O'], decrementBond=['O'], formBond=['O'], breakBond=['O'], incrementRadical=['O'], decrementRadical=['O'], incrementLonePair=['Os'], decrementLonePair=['Os']) atomTypes['Os' ].setActions(incrementBond=['Od'], decrementBond=[], formBond=['Os'], breakBond=['Os'], incrementRadical=['Os'], decrementRadical=['Os'], incrementLonePair=['Os'], decrementLonePair=['Os']) @@ -315,14 +315,14 @@ def getAtomType(atom, bonds): elif double == 0 and doubleO == 1 and triple == 0 and benzene == 0 and single == 1: atomType = 'N3d' elif double == 0 and doubleO == 0 and triple == 1 and benzene == 0 and single == 0: atomType = 'N3t' elif double == 0 and doubleO == 0 and triple == 0 and benzene == 2 and single == 0: atomType = 'N3b' - elif double == 0 and doubleO == 0 and triple == 0 and benzene == 0 and single == 4: atomType = 'N4s' - elif double == 1 and doubleO == 0 and triple == 0 and benzene == 0 and single == 2: atomType = 'N4d' - elif double == 0 and doubleO == 1 and triple == 0 and benzene == 0 and single == 2: atomType = 'N4d' - elif double == 2 and doubleO == 0 and triple == 0 and benzene == 0 and single == 0: atomType = 'N4dd' - elif double == 0 and doubleO == 2 and triple == 0 and benzene == 0 and single == 0: atomType = 'N4dd' - elif double == 1 and doubleO == 1 and triple == 0 and benzene == 0 and single == 0: atomType = 'N4dd' - elif double == 0 and doubleO == 0 and triple == 1 and benzene == 0 and single == 1: atomType = 'N4t' - elif double == 0 and doubleO == 0 and triple == 0 and benzene == 2 and single == 1: atomType = 'N4b' + elif double == 0 and doubleO == 0 and triple == 0 and benzene == 0 and single == 4: atomType = 'N5s' + elif double == 1 and doubleO == 0 and triple == 0 and benzene == 0 and single == 2: atomType = 'N5d' + elif double == 0 and doubleO == 1 and triple == 0 and benzene == 0 and single == 2: atomType = 'N5d' + elif double == 2 and doubleO == 0 and triple == 0 and benzene == 0 and single == 0: atomType = 'N5dd' + elif double == 0 and doubleO == 2 and triple == 0 and benzene == 0 and single == 0: atomType = 'N5dd' + elif double == 1 and doubleO == 1 and triple == 0 and benzene == 0 and single == 0: atomType = 'N5dd' + elif double == 0 and doubleO == 0 and triple == 1 and benzene == 0 and single == 1: atomType = 'N5t' + elif double == 0 and doubleO == 0 and triple == 0 and benzene == 2 and single == 1: atomType = 'N5b' elif atom.symbol == 'O': if double + doubleO == 0 and triple == 0 and benzene == 0: atomType = 'Os' elif double + doubleO == 1 and triple == 0 and benzene == 0: atomType = 'Od' From da085ec4555c2a8fe81f7d3ca62bd1a6bc961b56 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 4 Nov 2013 10:41:33 -0500 Subject: [PATCH 30/64] Change N4 to N5 --- rmgpy/molecule/molecule.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index a6a6012f91..be36a9f284 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1413,7 +1413,7 @@ def generateResonanceIsomers(self): newIsomers = isomer.getAdjacentResonanceIsomers() newIsomers += isomer.getLonePairRadicalResonanceIsomers() - newIsomers += isomer.getN4dd_N4tsResonanceIsomers() + newIsomers += isomer.getN5dd_N5tsResonanceIsomers() for newIsomer in newIsomers: newIsomer.updateAtomTypes() # Append to isomer list if unique @@ -1516,9 +1516,9 @@ def getLonePairRadicalResonanceIsomers(self): return isomers - def getN4dd_N4tsResonanceIsomers(self): + def getN5dd_N5tsResonanceIsomers(self): """ - Generate all of the resonance isomers formed by shifts between N4dd and N4ts. + Generate all of the resonance isomers formed by shifts between N5dd and N5ts. """ cython.declare(isomers=list, paths=list, index=cython.int, isomer=Molecule) cython.declare(atom=Atom, atom1=Atom, atom2=Atom, atom3=Atom) @@ -1529,9 +1529,9 @@ def getN4dd_N4tsResonanceIsomers(self): # Iterate over nitrogen atoms in structure for atom in self.vertices: - paths = self.findAllDelocalizationPathsN4dd_N4ts(atom) + paths = self.findAllDelocalizationPathsN5dd_N5ts(atom) for atom1, atom2, atom3, bond12, bond13, direction in paths: - # from N4dd to N4ts + # from N5dd to N5ts if direction == 1: # Adjust to (potentially) new resonance isomer bond12.decrementOrder() @@ -1563,7 +1563,7 @@ def getN4dd_N4tsResonanceIsomers(self): # Append to isomer list if unique isomers.append(isomer) - # from N4ts to N4dd + # from N5ts to N5dd if direction == 2: # Adjust to (potentially) new resonance isomer bond12.decrementOrder() @@ -1647,10 +1647,10 @@ def findAllDelocalizationPathsLonePairRadical(self, atom1): return paths - def findAllDelocalizationPathsN4dd_N4ts(self, atom1): + def findAllDelocalizationPathsN5dd_N5ts(self, atom1): """ - Find all the resonance structures of nitrogen atoms with two double bonds (N4dd) - and nitrogen atoms with one triple and one single bond (N4ts) + Find all the resonance structures of nitrogen atoms with two double bonds (N5dd) + and nitrogen atoms with one triple and one single bond (N5ts) """ cython.declare(paths=list) cython.declare(atom2=Atom, bond12=Bond) From 08ce20a584382f09a0e4b953fe6219034061d5c5 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 4 Nov 2013 10:43:18 -0500 Subject: [PATCH 31/64] Change N4 to N5 --- rmgpy/qm/mopac.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/qm/mopac.py b/rmgpy/qm/mopac.py index d5d38a2b2c..176720ec24 100644 --- a/rmgpy/qm/mopac.py +++ b/rmgpy/qm/mopac.py @@ -198,7 +198,7 @@ def generateQMData(self): Calculate the QM data and return a QMData object, or None if it fails. """ for atom in self.molecule.vertices: - if atom.atomType.label == 'N4s' or atom.atomType.label == 'N4d' or atom.atomType.label =='N4dd' or atom.atomType.label == 'N4t' or atom.atomType.label == 'N4b': + if atom.atomType.label == 'N5s' or atom.atomType.label == 'N5d' or atom.atomType.label =='N5dd' or atom.atomType.label == 'N5t' or atom.atomType.label == 'N5b': return None if self.verifyOutputFile(): From 6927c8dbf090ad705720980fd0cadb3dd0fd49ae Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 4 Nov 2013 12:36:41 -0500 Subject: [PATCH 32/64] Change N4 to N5 --- rmgpy/molecule/molecule.pxd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index a31ca6cccd..9b1d188e83 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -199,12 +199,12 @@ cdef class Molecule(Graph): cpdef list getLonePairRadicalResonanceIsomers(self) - cpdef list getN4dd_N4tsResonanceIsomers(self) + cpdef list getN5dd_N5tsResonanceIsomers(self) cpdef findAllDelocalizationPaths(self, Atom atom1) cpdef findAllDelocalizationPathsLonePairRadical(self, Atom atom1) - cpdef findAllDelocalizationPathsN4dd_N4ts(self, Atom atom1) + cpdef findAllDelocalizationPathsN5dd_N5ts(self, Atom atom1) cpdef int calculateSymmetryNumber(self) except -1 From cb2970a5cb910455c1c36576f8e1ec553d773b90 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 4 Nov 2013 20:51:57 -0500 Subject: [PATCH 33/64] Do not remove explicit hydrogen atoms --- rmgpy/data/base.py | 4 ++-- rmgpy/molecule/molecule.py | 2 +- rmgpy/qm/molecule.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/rmgpy/data/base.py b/rmgpy/data/base.py index ef6adedad2..860a5c6254 100644 --- a/rmgpy/data/base.py +++ b/rmgpy/data/base.py @@ -644,7 +644,7 @@ def getLogicNodeComponents(entry_or_item): for entry in entries: f.write(entry.label + '\n') if isinstance(entry.item, Molecule): - f.write(entry.item.toAdjacencyList(removeH=True) + '\n') + f.write(entry.item.toAdjacencyList(removeH=False) + '\n') elif isinstance(entry.item, Group): f.write(entry.item.toAdjacencyList().replace('{2S,2T}','2') + '\n') elif isinstance(entry.item, LogicOr): @@ -663,7 +663,7 @@ def comment(s): for entry in entriesNotInTree: f.write(comment(entry.label + '\n')) if isinstance(entry.item, Molecule): - f.write(comment(entry.item.toAdjacencyList(removeH=True) + '\n')) + f.write(comment(entry.item.toAdjacencyList(removeH=False) + '\n')) elif isinstance(entry.item, Group): f.write(comment(entry.item.toAdjacencyList().replace('{2S,2T}','2') + '\n')) elif isinstance(entry.item, LogicOr): diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index be36a9f284..25ecfc08db 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1691,7 +1691,7 @@ def getURL(self): # eg. http://dev.rmg.mit.edu/database/kinetics/reaction/reactant1=1%20C%200%20%7B2,S%7D;2%20O%200%20%7B1,S%7D;__reactant2=1%20C%202T;__product1=1%20C%201;__product2=1%20C%200%20%7B2,S%7D;2%20O%201%20%7B1,S%7D; url = "http://rmg.mit.edu/database/molecule/" - adjlist = self.toAdjacencyList(removeH=True) + adjlist = self.toAdjacencyList(removeH=False) url += "{0}".format(re.sub('\s+', '%20', adjlist.replace('\n', ';'))) return url.strip('_') diff --git a/rmgpy/qm/molecule.py b/rmgpy/qm/molecule.py index 97e8f98c2f..9d45bd7597 100644 --- a/rmgpy/qm/molecule.py +++ b/rmgpy/qm/molecule.py @@ -250,7 +250,7 @@ def saveThermoData(self): resultFile.write("thermoData = {0!r}\n".format(self.thermo)) resultFile.write("pointGroup = {0!r}\n".format(self.pointGroup)) resultFile.write("qmData = {0!r}\n".format(self.qmData)) - resultFile.write('adjacencyList = """\n{0!s}"""\n'.format(self.molecule.toAdjacencyList(removeH=True))) + resultFile.write('adjacencyList = """\n{0!s}"""\n'.format(self.molecule.toAdjacencyList(removeH=False))) def loadThermoData(self): """ From 5ce1b76cb40a7b314d14918d085c07e73d20a45e Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 11 Nov 2013 09:33:10 -0500 Subject: [PATCH 34/64] Exclude radicals with one heavy atom from maxRadicals limitation --- rmgpy/data/kinetics/family.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 7a2aee015f..1eb0c764fa 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1144,7 +1144,7 @@ def __generateProductStructures(self, reactantStructures, maps, forward, **optio raise ForbiddenStructureException() if len(struct.atoms) - H > maxHeavyAtoms: raise ForbiddenStructureException() - if struct.getNumberOfRadicalElectrons() > maxRadicals: + if (struct.getNumberOfRadicalElectrons() > maxRadicals) and (len(struct.atoms) - H < 1): raise ForbiddenStructureException() # Check that product structures are allowed in this family From d9441079213c682be6d9e1c1f79dd609e8b65b77 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 11 Nov 2013 10:47:36 -0500 Subject: [PATCH 35/64] Add functions to determine if biradical is singlet and is triplet --- rmgpy/molecule/molecule.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 25ecfc08db..7bbdcf8260 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1694,4 +1694,28 @@ def getURL(self): adjlist = self.toAdjacencyList(removeH=False) url += "{0}".format(re.sub('\s+', '%20', adjlist.replace('\n', ';'))) return url.strip('_') + + def isBiradicalSinglet(self): + """ + Return ``True`` if the molecule is a 1-centered biradical in singlet state, + or ``False`` otherwise. + """ + cython.declare(atom=Atom) + for atom in self.vertices: + if atom.radicalElectrons == 2: + if atom.spinMultiplicity == 1: + return True + return False + + def isBiradicalTriplet(self): + """ + Return ``True`` if the molecule is a 1-centered biradical in triplet state, + or ``False`` otherwise. + """ + cython.declare(atom=Atom) + for atom in self.vertices: + if atom.radicalElectrons == 2: + if atom.spinMultiplicity == 3: + return True + return False From c3639e097ca1eda8069f89f9aa58eef133d74df4 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 11 Nov 2013 13:27:20 -0500 Subject: [PATCH 36/64] Addding function changing triplet radical to singlet radical --- rmgpy/molecule/molecule.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 7bbdcf8260..d8b56738f0 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1718,4 +1718,15 @@ def isBiradicalTriplet(self): if atom.spinMultiplicity == 3: return True return False + + def changeTripletSinglet(self): + """ + If the molecule is a 1-centered biradical in triplet state, + change it to singlet state. + """ + cython.declare(atom=Atom) + for atom in self.vertices: + if atom.radicalElectrons == 2: + if atom.spinMultiplicity == 3: + atom.spinMultiplicity = 1 From f0d2f235458678e317a3c69b97a02792eb5fa214 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 11 Nov 2013 13:29:40 -0500 Subject: [PATCH 37/64] Function checking reverse reactions of CH2(S) for correct multiplicity --- rmgpy/data/kinetics/family.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 1eb0c764fa..00a23f08f9 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1380,9 +1380,18 @@ def __generateReactions(self, reactants, products=None, forward=True, **options) for reaction in rxnList0: products0 = reaction.products if forward else reaction.reactants + + # If the forward reaction contains CH(S) as reactant then correct product of reverse reaction from CH2(T) to CH2(S) + for product in products: + for molecule in product: + if molecule.isBiradicalSinglet() and molecule.getFormula() == 'CH2': + for product0 in products0: + if product0.isBiradicalTriplet() and product0.getFormula() == 'CH2': + product0.changeTripletSinglet() # Skip reactions that don't match the given products match = False + if len(products) == len(products0) == 1: for product in products[0]: if products0[0].isIsomorphic(product): From ab3642a4dcce5d9163c0c1a1104c2412a6d7d274 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 13 Nov 2013 17:37:40 -0500 Subject: [PATCH 38/64] Correct BAC for CCSD(T)-F12/cc-pVTZ-F12 'C-O' --- rmgpy/cantherm/statmech.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/cantherm/statmech.py b/rmgpy/cantherm/statmech.py index 8d57110956..da44c32f46 100644 --- a/rmgpy/cantherm/statmech.py +++ b/rmgpy/cantherm/statmech.py @@ -609,7 +609,7 @@ def applyEnergyCorrections(E0, modelChemistry, atoms, bonds): 'N-H': -0.69, 'N-N': -0.47, 'N=N': -1.54, 'N#N': -2.05,} elif modelChemistry == 'CCSD(T)-F12/cc-pVTZ-F12': bondEnergies = { 'C-H': -0.09, 'C-C': -0.27, 'C=C': -1.03, 'C#C': -1.79, - 'O-H': -0.06, 'C-O': -0.14, 'C=O': -0.19, 'O-O': 0.16, 'N-C': -0.18, + 'O-H': -0.06, 'C-O': 0.14, 'C=O': -0.19, 'O-O': 0.16, 'N-C': -0.18, 'N=C': -0.41, 'N#C': -1.41, 'N-O': 0.87, 'N_O': -0.09, 'N=O': -0.23, 'N-H': -0.01, 'N-N': -0.21, 'N=N': -0.44, 'N#N': -0.76,} elif modelChemistry == 'CCSD(T)-F12/cc-pVQZ-F12': From 611f32a6c66a4158db5f6f9bf3a21d50586e04b7 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Thu, 14 Nov 2013 13:34:36 -0500 Subject: [PATCH 39/64] toInChI: Fall back to openbabel of RDKit fails for nitrogen --- rmgpy/molecule/molecule.py | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index d8b56738f0..26cd76e1cc 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1094,11 +1094,26 @@ def toInChI(self): Convert a molecular structure to an InChI string. Uses `RDKit `_ to perform the conversion. Perceives aromaticity. + + or + + Convert a molecular structure to an InChI string. Uses + `OpenBabel `_ to perform the conversion. """ - if not Chem.inchi.INCHI_AVAILABLE: - return "RDKitInstalledWithoutInChI" - rdkitmol = self.toRDKitMol() - return Chem.inchi.MolToInchi(rdkitmol) + try: + if not Chem.inchi.INCHI_AVAILABLE: + return "RDKitInstalledWithoutInChI" + rdkitmol = self.toRDKitMol() + return Chem.inchi.MolToInchi(rdkitmol) + except: + pass + + # This version does not write a warning to stderr if stereochemistry is undefined + obmol = self.toOBMol() + obConversion = openbabel.OBConversion() + obConversion.SetOutFormat('inchi') + obConversion.SetOptions('w', openbabel.OBConversion.OUTOPTIONS) + return obConversion.WriteString(obmol).strip() def toAugmentedInChI(self): """ From 5d7d9be98bab14bcea9acec82bb734edcc8f54ba Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Thu, 14 Nov 2013 17:39:35 -0500 Subject: [PATCH 40/64] Make OpenBabel requirement optional RMG-Py runs without OpenBabel relying on RDKit. Tetravalent nitrogen requires OpenBabel is required as fall-back from RdKit in the drawing functions and functions like toInChI or toSMiles. --- rmgpy/molecule/molecule.py | 38 +++++++++++++++++++++----------------- 1 file changed, 21 insertions(+), 17 deletions(-) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 26cd76e1cc..b4bba197af 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -40,7 +40,10 @@ import os import re import element as elements -import openbabel +try: + import openbabel +except: + pass from rdkit import Chem from .graph import Vertex, Edge, Graph from .group import GroupAtom, GroupBond, Group, ActionError @@ -1142,24 +1145,25 @@ def toInChIKey(self): Removes check-sum dash (-) and character so that only the 14 + 9 characters remain. """ - + try: + if not Chem.inchi.INCHI_AVAILABLE: + return "RDKitInstalledWithoutInChI" + inchi = self.toInChI() + return Chem.inchi.InchiToInchiKey(inchi)[:-2] + except: + pass + import openbabel - - for atom in self.vertices: - if atom.isNitrogen(): - # This version does not write a warning to stderr if stereochemistry is undefined - obmol = self.toOBMol() - obConversion = openbabel.OBConversion() - obConversion.SetOutFormat('inchi') - obConversion.SetOptions('w', openbabel.OBConversion.OUTOPTIONS) - obConversion.SetOptions('K', openbabel.OBConversion.OUTOPTIONS) - return obConversion.WriteString(obmol).strip()[:-2] - - if not Chem.inchi.INCHI_AVAILABLE: - return "RDKitInstalledWithoutInChI" - inchi = self.toInChI() - return Chem.inchi.InchiToInchiKey(inchi)[:-2] +# for atom in self.vertices: + # if atom.isNitrogen(): + # This version does not write a warning to stderr if stereochemistry is undefined + obmol = self.toOBMol() + obConversion = openbabel.OBConversion() + obConversion.SetOutFormat('inchi') + obConversion.SetOptions('w', openbabel.OBConversion.OUTOPTIONS) + obConversion.SetOptions('K', openbabel.OBConversion.OUTOPTIONS) + return obConversion.WriteString(obmol).strip()[:-2] def toAugmentedInChIKey(self): """ From abe18ebb806e72705e363a8fbaf38c2c5b437bdd Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Tue, 19 Nov 2013 12:19:21 -0500 Subject: [PATCH 41/64] replace getformula() with SMILEwriter using openbabel for nitrogen --- rmgpy/molecule/molecule.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index b4bba197af..501a2564b8 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -525,6 +525,12 @@ def applyAction(self, action): raise ActionError('Unable to update GroupBond: Invalid action {0}.'.format(action)) ################################################################################# +try: + SMILEwriter = openbabel.OBConversion() + SMILEwriter.SetOutFormat('smi') + SMILEwriter.SetOptions("i",SMILEwriter.OUTOPTIONS) # turn off isomer and stereochemistry information (the @ signs!) +except: + pass class Molecule(Graph): """ @@ -1191,6 +1197,11 @@ def toSMARTS(self): def toSMILES(self): """ + Convert a molecular structure to an SMILES string. Uses + `OpenBabel `_ to perform the conversion. + + or + Convert a molecular structure to a canonical SMILES string. Uses `RDKit `_ to perform the conversion. Perceives aromaticity and removes Hydrogen atoms. @@ -1198,7 +1209,12 @@ def toSMILES(self): for atom in self.vertices: if atom.isNitrogen(): - return self.getFormula() + mol = self.toOBMol() + if self.getFormula() == 'H2': + return '[H][H]' + elif self.getFormula() == 'H': + return '[H]' + return SMILEwriter.WriteString(mol).strip() rdkitmol = self.toRDKitMol() From 118befdad4c1728226c6f314b735d646c06eed71 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 20 Nov 2013 15:40:37 -0500 Subject: [PATCH 42/64] Extend fingerprint generation to account for nitrogen atoms --- rmgpy/molecule/group.pxd | 1 + rmgpy/molecule/group.py | 37 +++++++++++++++++++++---------------- rmgpy/molecule/molecule.py | 7 +++++-- 3 files changed, 27 insertions(+), 18 deletions(-) diff --git a/rmgpy/molecule/group.pxd b/rmgpy/molecule/group.pxd index b059c6a627..9939b85e07 100644 --- a/rmgpy/molecule/group.pxd +++ b/rmgpy/molecule/group.pxd @@ -83,6 +83,7 @@ cdef class Group(Graph): # These read-only attribues act as a "fingerprint" for accelerating # subgraph isomorphism checks cdef public short carbonCount + cdef public short nitrogenCount cdef public short oxygenCount cdef public short sulfurCount cdef public short radicalCount diff --git a/rmgpy/molecule/group.py b/rmgpy/molecule/group.py index 4431318134..b6f83e501c 100644 --- a/rmgpy/molecule/group.py +++ b/rmgpy/molecule/group.py @@ -673,28 +673,33 @@ def updateFingerprint(self): isomorphism checks. """ cython.declare(atom=GroupAtom, atomType=AtomType) - cython.declare(carbon=AtomType, oxygen=AtomType, sulfur=AtomType) - cython.declare(isCarbon=cython.bint, isOxygen=cython.bint, isSulfur=cython.bint, radical=cython.int) + cython.declare(carbon=AtomType, nitrogen=AtomType, oxygen=AtomType, sulfur=AtomType) + cython.declare(isCarbon=cython.bint, isNitrogen=cython.bint, isOxygen=cython.bint, isSulfur=cython.bint, radical=cython.int) - carbon = atomTypes['C'] - oxygen = atomTypes['O'] - sulfur = atomTypes['S'] + carbon = atomTypes['C'] + nitrogen = atomTypes['N'] + oxygen = atomTypes['O'] + sulfur = atomTypes['S'] - self.carbonCount = 0 - self.oxygenCount = 0 - self.sulfurCount = 0 - self.radicalCount = 0 + self.carbonCount = 0 + self.nitrogenCount = 0 + self.oxygenCount = 0 + self.sulfurCount = 0 + self.radicalCount = 0 for atom in self.vertices: if len(atom.atomType) == 1: - atomType = atom.atomType[0] - isCarbon = atomType.equivalent(carbon) - isOxygen = atomType.equivalent(oxygen) - isSulfur = atomType.equivalent(sulfur) - if isCarbon and not isOxygen and not isSulfur: + atomType = atom.atomType[0] + isCarbon = atomType.equivalent(carbon) + isNitrogen = atomType.equivalent(nitrogen) + isOxygen = atomType.equivalent(oxygen) + isSulfur = atomType.equivalent(sulfur) + if isCarbon and not isNitrogen and not isOxygen and not isSulfur: self.carbonCount += 1 - elif isOxygen and not isCarbon and not isSulfur: + elif isNitrogen and not isCarbon and not isOxygen and not isSulfur: + self.nitrogenCount += 1 + elif isOxygen and not isCarbon and not isNitrogen and not isSulfur: self.oxygenCount += 1 - elif isSulfur and not isCarbon and not isOxygen: + elif isSulfur and not isCarbon and not isNitrogen and not isOxygen: self.sulfurCount += 1 if len(atom.radicalElectrons) == 1: radical = atom.radicalElectrons[0] diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 501a2564b8..e5983d2956 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -892,7 +892,7 @@ def isSubgraphIsomorphic(self, other, initialMap=None): be a :class:`Group` object, or a :class:`TypeError` is raised. """ cython.declare(group=Group, atom=Atom) - cython.declare(carbonCount=cython.short, oxygenCount=cython.short, sulfurCount=cython.short, radicalCount=cython.short) + cython.declare(carbonCount=cython.short, nitrogenCount=cython.short, oxygenCount=cython.short, sulfurCount=cython.short, radicalCount=cython.short) # It only makes sense to compare a Molecule to a Group for subgraph # isomorphism, so raise an exception if this is not what was requested @@ -901,10 +901,12 @@ def isSubgraphIsomorphic(self, other, initialMap=None): group = other # Count the number of carbons, oxygens, and radicals in the molecule - carbonCount = 0; oxygenCount = 0; sulfurCount = 0; radicalCount = 0 + carbonCount = 0; nitrogenCount = 0; oxygenCount = 0; sulfurCount = 0; radicalCount = 0 for atom in self.vertices: if atom.element.symbol == 'C': carbonCount += 1 + elif atom.element.symbol == 'N': + nitrogenCount += 1 elif atom.element.symbol == 'O': oxygenCount += 1 elif atom.element.symbol == 'S': @@ -915,6 +917,7 @@ def isSubgraphIsomorphic(self, other, initialMap=None): # needing to perform the full isomorphism check if (radicalCount < group.radicalCount or carbonCount < group.carbonCount or + nitrogenCount < group.nitrogenCount or oxygenCount < group.oxygenCount or sulfurCount < group.sulfurCount): return False From f71ccc9f53e2c96d0eaf7c8473860c250ab49e64 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 20 Nov 2013 17:56:02 -0500 Subject: [PATCH 43/64] Include NH(S) into quick fix for reverse reaction of singlet biradicals for CH2(S) see bbuesser/RMG-Py@f0d2f23 --- rmgpy/data/kinetics/family.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 5d97891288..19e9970ccd 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1384,9 +1384,9 @@ def __generateReactions(self, reactants, products=None, forward=True, **options) # If the forward reaction contains CH(S) as reactant then correct product of reverse reaction from CH2(T) to CH2(S) for product in products: for molecule in product: - if molecule.isBiradicalSinglet() and molecule.getFormula() == 'CH2': + if molecule.isBiradicalSinglet() and (molecule.getFormula() == 'CH2' or molecule.getFormula() == 'NH'): for product0 in products0: - if product0.isBiradicalTriplet() and product0.getFormula() == 'CH2': + if product0.isBiradicalTriplet() and (product0.getFormula() == 'CH2' or product0.getFormula() == 'NH'): product0.changeTripletSinglet() # Skip reactions that don't match the given products From b4d8b861dd36f279816f80d5b7d2e12248ec7356 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 20 Nov 2013 21:07:12 -0500 Subject: [PATCH 44/64] fixing failing testToAdjacencyList() The number of lone electron pairs is not defined for groups, therefore should not be printed for groups --- rmgpy/molecule/adjlist.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index fe555abbe0..18b4b46adc 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -341,8 +341,8 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False): atomLabelWidth = max([len(s) for s in atomLabels.values()]) if atomLabelWidth > 0: atomLabelWidth += 1 atomTypeWidth = max([len(s) for s in atomTypes.values()]) + 1 - atomElectronStateWidth = max([len(s) for s in atomElectronStates.values()]) + 1 - atomLonePairWidth = 1 + atomElectronStateWidth = max([len(s) for s in atomElectronStates.values()]) + atomLonePairWidth = 2 # Assemble the adjacency list for atom in atoms: @@ -356,8 +356,9 @@ def toAdjacencyList(atoms, label=None, group=False, removeH=False): adjlist += '{0:<{1:d}}'.format(atomTypes[atom], atomTypeWidth) # Electron state(s) adjlist += '{0:<{1:d}}'.format(atomElectronStates[atom], atomElectronStateWidth) - # Lone Pair(s) - adjlist += '{0:<{1:d}}'.format(atomLonePairs[atom], atomLonePairWidth) + if group == False: + # Lone Pair(s) + adjlist += '{0:>{1:d}}'.format(atomLonePairs[atom], atomLonePairWidth) # Bonds list atoms2 = atom.bonds.keys() From f313e8a6112795f94bf87aec83a13af2fe5fb598 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 20 Nov 2013 21:25:22 -0500 Subject: [PATCH 45/64] fixing the failing test testAdjacencyList() Add explicit hydrogen. --- rmgpy/molecule/moleculeTest.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/rmgpy/molecule/moleculeTest.py b/rmgpy/molecule/moleculeTest.py index 0acd8f1e6e..cef813d5df 100644 --- a/rmgpy/molecule/moleculeTest.py +++ b/rmgpy/molecule/moleculeTest.py @@ -686,12 +686,21 @@ def testAdjacencyList(self): Check the adjacency list read/write functions for a full molecule. """ molecule1 = Molecule().fromAdjacencyList(""" - 1 C 0 {2,D} - 2 C 0 {1,D} {3,S} - 3 C 0 {2,S} {4,D} - 4 C 0 {3,D} {5,S} - 5 C 1 {4,S} {6,S} - 6 C 0 {5,S} + 1 C 0 {2,D} {7,S} {8,S} + 2 C 0 {1,D} {3,S} {9,S} + 3 C 0 {2,S} {4,D} {10,S} + 4 C 0 {3,D} {5,S} {11,S} + 5 C 1 {4,S} {6,S} {12,S} + 6 C 0 {5,S} {13,S} {14,S} {15,S} + 7 H 0 {1,S} + 8 H 0 {1,S} + 9 H 0 {2,S} + 10 H 0 {3,S} + 11 H 0 {4,S} + 12 H 0 {5,S} + 13 H 0 {6,S} + 14 H 0 {6,S} + 15 H 0 {6,S} """) molecule2 = Molecule().fromSMILES('C=CC=C[CH]C') self.assertTrue(molecule1.isIsomorphic(molecule2)) From 97192bc12611d8875e41b2d2db2616accf9ce771 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 20 Nov 2013 21:27:38 -0500 Subject: [PATCH 46/64] add explicit hydrogen to setup test molecule to fix failing test --- rmgpy/molecule/moleculeTest.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/rmgpy/molecule/moleculeTest.py b/rmgpy/molecule/moleculeTest.py index cef813d5df..178430a4f0 100644 --- a/rmgpy/molecule/moleculeTest.py +++ b/rmgpy/molecule/moleculeTest.py @@ -471,7 +471,10 @@ def setUp(self): self.adjlist = """ 1 *2 C 1 {2,D} {3,S} 2 *1 O 0 {1,D} -3 C 0 {1,S} +3 C 0 {1,S} {4,S} {5,S} {6,S} +4 H 0 {3,S} +5 H 0 {3,S} +6 H 0 {3,S} """ self.molecule = Molecule().fromAdjacencyList(self.adjlist) From d9425b035aa7cea5800c05baf5a0dc565ba3cded Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 20 Nov 2013 21:40:53 -0500 Subject: [PATCH 47/64] fix testToAdjacencyList set removeH = False add number of lone electron pairs to setup molecule --- rmgpy/molecule/moleculeTest.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/rmgpy/molecule/moleculeTest.py b/rmgpy/molecule/moleculeTest.py index 178430a4f0..94d9287ec2 100644 --- a/rmgpy/molecule/moleculeTest.py +++ b/rmgpy/molecule/moleculeTest.py @@ -469,12 +469,12 @@ class TestMolecule(unittest.TestCase): def setUp(self): self.adjlist = """ -1 *2 C 1 {2,D} {3,S} -2 *1 O 0 {1,D} -3 C 0 {1,S} {4,S} {5,S} {6,S} -4 H 0 {3,S} -5 H 0 {3,S} -6 H 0 {3,S} +1 *2 C 1 0 {2,D} {3,S} +2 *1 O 0 2 {1,D} +3 C 0 0 {1,S} {4,S} {5,S} {6,S} +4 H 0 0 {3,S} +5 H 0 0 {3,S} +6 H 0 0 {3,S} """ self.molecule = Molecule().fromAdjacencyList(self.adjlist) @@ -579,7 +579,7 @@ def testToAdjacencyList(self): """ Test the Molecule.toAdjacencyList() method. """ - adjlist = self.molecule.toAdjacencyList(removeH=True) + adjlist = self.molecule.toAdjacencyList(removeH=False) self.assertEqual(adjlist.strip(), self.adjlist.strip()) def testIsomorphism(self): From 6bdd7068616d88149842fd99d27df78a74c47e1a Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 20 Nov 2013 21:46:08 -0500 Subject: [PATCH 48/64] fix testToAdjacencyList() of speciesTest.py set removeH = False --- rmgpy/speciesTest.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/rmgpy/speciesTest.py b/rmgpy/speciesTest.py index 67915c1935..1a2b8dbd22 100644 --- a/rmgpy/speciesTest.py +++ b/rmgpy/speciesTest.py @@ -102,8 +102,9 @@ def testToAdjacencyList(self): """ Test that toAdjacencyList() works as expected. """ + print self.species.toAdjacencyList() string = self.species.toAdjacencyList() - self.assertTrue(string.startswith(self.species.molecule[0].toAdjacencyList(label=self.species.label,removeH=True))) + self.assertTrue(string.startswith(self.species.molecule[0].toAdjacencyList(label=self.species.label,removeH=False))) ################################################################################ From 0714b4f46a1a4b2c659d196631f2b5fb81955681 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 20 Nov 2013 22:14:11 -0500 Subject: [PATCH 49/64] fix testSubgraphIsomorphismManyLabels() of moleculeTest.py --- rmgpy/molecule/moleculeTest.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/rmgpy/molecule/moleculeTest.py b/rmgpy/molecule/moleculeTest.py index 94d9287ec2..3116d1e945 100644 --- a/rmgpy/molecule/moleculeTest.py +++ b/rmgpy/molecule/moleculeTest.py @@ -657,16 +657,22 @@ def testSubgraphIsomorphismAgain(self): def testSubgraphIsomorphismManyLabels(self): molecule = Molecule() # specific case (species) molecule.fromAdjacencyList(""" -1 *1 C 1 {2,S} {3,S} -2 C 0 {1,S} {3,S} -3 C 0 {1,S} {2,S} +1 *1 C 1 {2,S} {3,S} {4,S} +2 C 0 {1,S} {3,S} {5,S} {6,S} +3 C 0 {1,S} {2,S} {7,S} {8,S} +4 H 0 {1,S} +5 H 0 {2,S} +6 H 0 {2,S} +7 H 0 {3,S} +8 H 0 {3,S} """) + print molecule.toAdjacencyList() group = Group() # general case (functional group) group.fromAdjacencyList(""" -1 *1 C 1 {2,S}, {3,S} -2 R 0 {1,S} -3 R 0 {1,S} +1 *1 C 1 {2,S}, {3,S} +2 R!H 0 {1,S} +3 R!H 0 {1,S} """) labeled1 = molecule.getLabeledAtoms() @@ -677,7 +683,7 @@ def testSubgraphIsomorphismManyLabels(self): self.assertTrue(molecule.isSubgraphIsomorphic(group, initialMap)) mapping = molecule.findSubgraphIsomorphisms(group, initialMap) - self.assertEqual(len(mapping), 1) + self.assertEqual(len(mapping), 2) for map in mapping: self.assertTrue(len(map) == min(len(molecule.atoms), len(group.atoms))) for key, value in map.iteritems(): From dcdce81ab3e27dfe3e03b3dcdab6a7468fd99e85 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 20 Nov 2013 22:33:30 -0500 Subject: [PATCH 50/64] fix testSetActions() in atomtypeTest.py --- rmgpy/molecule/atomtypeTest.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/rmgpy/molecule/atomtypeTest.py b/rmgpy/molecule/atomtypeTest.py index 6b1dc3ded6..80e8b87eef 100644 --- a/rmgpy/molecule/atomtypeTest.py +++ b/rmgpy/molecule/atomtypeTest.py @@ -88,7 +88,9 @@ def testSetActions(self): self.atomType.formBond, self.atomType.breakBond, self.atomType.incrementRadical, - self.atomType.decrementRadical) + self.atomType.decrementRadical, + self.atomType.incrementLonePair, + self.atomType.decrementLonePair) self.assertEqual(self.atomType.incrementBond, other.incrementBond) self.assertEqual(self.atomType.decrementBond, other.decrementBond) self.assertEqual(self.atomType.formBond, other.formBond) From 2528de2388cd99a95e1e7eb154b06cfe2f6bbed5 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Thu, 21 Nov 2013 01:45:57 -0500 Subject: [PATCH 51/64] A preliminary solution for the reverse reaction of CH2(T) + CH Introducing 1-centered biradicals to the H-abstraction family lets RMG find the special case of CH2(T) + CH = CH + CH2(T). However for the reverse reaction it recognizes that nothing happens and removes the reverse reaction from the reaction list which leads an error for empty reaction lists. --- rmgpy/data/kinetics/family.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 19e9970ccd..58a63f0d56 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1179,9 +1179,9 @@ def __createReaction(self, reactants, products, isForward): if reactants[0].isIsomorphic(products[0]): return None elif len(reactants) == len(products) == 2: - if reactants[0].isIsomorphic(products[0]) and reactants[1].isIsomorphic(products[1]): + if reactants[0].isIsomorphic(products[0]) and reactants[1].isIsomorphic(products[1]) and reactants[0].getFormula() != ('CH2' or 'CH') and reactants[1].getFormula() != ('CH2' or 'CH'): return None - elif reactants[0].isIsomorphic(products[1]) and reactants[1].isIsomorphic(products[0]): + elif reactants[0].isIsomorphic(products[1]) and reactants[1].isIsomorphic(products[0]) and reactants[0].getFormula() != ('CH2' or 'CH') and reactants[1].getFormula() != ('CH2' or 'CH'): return None # Create and return template reaction object From 65084cacf3eba4aa10dafa59eb0e6ef35b669179 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Thu, 21 Nov 2013 15:45:56 -0500 Subject: [PATCH 52/64] Example input for RMG-Py including nitrogen and lone electron pairs - uses the newly added libraries --- examples/rmg/ch3no2/input.py | 86 ++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) create mode 100644 examples/rmg/ch3no2/input.py diff --git a/examples/rmg/ch3no2/input.py b/examples/rmg/ch3no2/input.py new file mode 100644 index 0000000000..6490f928b2 --- /dev/null +++ b/examples/rmg/ch3no2/input.py @@ -0,0 +1,86 @@ +# Data sources +database( + thermoLibraries = ['KlippensteinH2O2', 'primaryThermoLibrary','DFT_QCI_thermo','CH','CHN','CHO','CHON','CN','NISTThermoLibrary','thermo_DFT_CCSDTF12_BAC','GRI-Mech3.0-N'], + reactionLibraries = [('Nitrogen_Dean_and_Bozelli',False)], + seedMechanisms = ['ERC-FoundationFuelv0.9'], + kineticsDepositories = ['training'], + kineticsFamilies = ['!Intra_Disproportionation', '!Substitution_O'], + kineticsEstimator = 'rate rules', +) + +# Constraints on generated species +generatedSpeciesConstraints( + #maximumCarbonAtoms = 7, + #maximumHydrogenAtoms = 8, + #maximumOxygenAtoms = 5, + maximumNitrogenAtoms = 2, + #maximumSiliconAtoms = 0, + #maximumSulfurAtoms = 0, + #maximumHeavyAtoms = 3, + maximumRadicalElectrons = 2, +) + +# List of species +species( + label='CH3NO2', + reactive=True, + structure=adjacencyList( + """ + 1 C 0 0 {2,S} {3,S} {4,S} {5,S} + 2 H 0 0 {1,S} + 3 H 0 0 {1,S} + 4 H 0 0 {1,S} + 5 N 0 0 {1,S} {6,D} {7,S} + 6 O 0 2 {5,D} + 7 O 0 3 {5,S} + """), +) + +species( + label='O2', + reactive=True, + structure=adjacencyList( + """ + 1 O 1 2 {2,S} + 2 O 1 2 {1,S} + """), +) + +species( + label='N2', + reactive=True, + structure=adjacencyList( + """ + 1 N 1 1 {2,T} + 2 N 1 1 {1,T} + """), +) + +# Reaction systems + +simpleReactor( + temperature=(1500,'K'), + pressure=(10.0,'bar'), + initialMoleFractions={ + "CH3NO2": 0.1, + "O2": 0.21, + "N2": 0.69, + }, + terminationConversion={ + 'CH3NO2': 0.1, + }, +) + +model( + toleranceKeepInEdge=1e-5, + toleranceMoveToCore=0.1, + toleranceInterruptSimulation=0.1, + maximumEdgeSpecies=10000 +) + +options( + units='si', + saveRestartPeriod=None, + drawMolecules=False, + generatePlots=False, +) From a575e4fe1d3c9b981cc28a95078dd64181af0269 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Thu, 21 Nov 2013 16:34:01 -0500 Subject: [PATCH 53/64] Add additional electronic states S - singlet D - doublet T - triplet Q - quartet V - quintet --- rmgpy/molecule/adjlist.py | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/rmgpy/molecule/adjlist.py b/rmgpy/molecule/adjlist.py index 18b4b46adc..b8513aae63 100644 --- a/rmgpy/molecule/adjlist.py +++ b/rmgpy/molecule/adjlist.py @@ -139,8 +139,18 @@ def fromAdjacencyList(adjlist, group=False): radicalElectrons.append(2); spinMultiplicity.append(3) elif e == '3': radicalElectrons.append(3); spinMultiplicity.append(4) + elif e == '3D': + radicalElectrons.append(3); spinMultiplicity.append(2) + elif e == '3Q': + radicalElectrons.append(3); spinMultiplicity.append(4) elif e == '4': radicalElectrons.append(4); spinMultiplicity.append(5) + elif e == '4S': + radicalElectrons.append(4); spinMultiplicity.append(1) + elif e == '4T': + radicalElectrons.append(4); spinMultiplicity.append(3) + elif e == '4V': + radicalElectrons.append(4); spinMultiplicity.append(5) index += 1 # Next number defines the number of lone electron pairs (if provided) @@ -279,10 +289,16 @@ def getElectronState(radicalElectrons, spinMultiplicity): electronState = '2S' elif radicalElectrons == 2 and spinMultiplicity == 3: electronState = '2T' - elif radicalElectrons == 3: - electronState = '3' - elif radicalElectrons == 4: - electronState = '4' + elif radicalElectrons == 3 and spinMultiplicity == 2: + electronState = '3D' + elif radicalElectrons == 3 and spinMultiplicity == 4: + electronState = '3Q' + elif radicalElectrons == 4 and spinMultiplicity == 1: + electronState = '4S' + elif radicalElectrons == 4 and spinMultiplicity == 3: + electronState = '4T' + elif radicalElectrons == 4 and spinMultiplicity == 5: + electronState = '4V' else: raise ValueError('Unable to determine electron state for {0:d} radical electrons with spin multiplicity of {1:d}.'.format(radicalElectrons, spinMultiplicity)) return electronState From 5ea2e403299281d4c6b22cb10597d0a482df849a Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Fri, 22 Nov 2013 12:55:00 -0500 Subject: [PATCH 54/64] Add function to set spin multiplicity of an atom --- rmgpy/molecule/molecule.pxd | 2 ++ rmgpy/molecule/molecule.py | 11 +++++++++++ 2 files changed, 13 insertions(+) diff --git a/rmgpy/molecule/molecule.pxd b/rmgpy/molecule/molecule.pxd index 9b1d188e83..743419fc04 100644 --- a/rmgpy/molecule/molecule.pxd +++ b/rmgpy/molecule/molecule.pxd @@ -69,6 +69,8 @@ cdef class Atom(Vertex): cpdef updateCharge(self) + cpdef setSpinMultiplicity(self, int spinMultiplicity) + ################################################################################ cpdef object SMILEwriter diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index e5983d2956..bc62ed285e 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -359,6 +359,17 @@ def applyAction(self, action): for i in range(abs(action[2])): self.decrementLonePairs() else: raise ActionError('Unable to update Atom: Invalid action {0}".'.format(action)) + + def setSpinMultiplicity(self,spinMultiplicity): + """ + Set the spin multiplicity. + """ + # Set the spin multiplicity + self.spinMultiplicity = spinMultiplicity + if self.spinMultiplicity < 0: + raise ActionError('Unable to update Atom due to spin multiplicity : Invalid spin multiplicity set "{0}".'.format(self.spinMultiplicity)) + self.updateCharge() + ################################################################################ From 59e0a7bebecddf0bf328e16722953ee7ead0e8d5 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Sun, 24 Nov 2013 10:33:22 -0500 Subject: [PATCH 55/64] Finding all possible electronic states in product structures MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The current assumptions for newly found products is that they attain the highest possible multiplicity as this in most cases the most stable species and resulting often in the fastest kinetics. However if we start with a singlet biradical reactant in the forward direction we find with this assumption only the triplet biradical as product of the reverse reaction. I became aware of this because I have added the nitrogen atom as a triradical to the three important reaction families H_abstraction, Disproportionation and R_Addition_MultipleBond. I think making it very general is the simplest solution. For that I have added the electronic states of quartet (Q) and quintet (V). I didn’t want to continue hardcode specific exceptions. Therefore I have modified the function __generateProductStructures() to find all possible electronic states for the newly found product structures and add all possible combinations of these products to a product structures list. This list can also contain spin-forbidden reactions! But instead of hardcoding the removal of them I think we should add their often very slow kinetics to the library. I have extended __generateReactions() to make a list of reactions using the new list of product structures. That way it will find the reverse reaction leading to the correct electronic state of the reactants of the forward reaction. --- rmgpy/data/kinetics/family.py | 96 +++++++++++++++++++++++------------ 1 file changed, 63 insertions(+), 33 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 58a63f0d56..045a497a16 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1081,11 +1081,10 @@ def __generateProductStructures(self, reactantStructures, maps, forward, **optio reactants are stored in the reaction family template. The `maps` parameter is a list of mappings of the top-level tree node of each *template* reactant to the corresponding *structure*. This function - returns the product structures. + returns a list of the product structures. """ - - if not forward: template = self.reverseTemplate - else: template = self.forwardTemplate + + productStructuresList = [] # Clear any previous atom labeling from all reactant structures for struct in reactantStructures: struct.clearLabeledAtoms() @@ -1098,7 +1097,8 @@ def __generateProductStructures(self, reactantStructures, maps, forward, **optio # Check that reactant structures are allowed in this family # If not, then stop for struct in reactantStructures: - if self.isMoleculeForbidden(struct): raise ForbiddenStructureException() + if self.isMoleculeForbidden(struct) and forward: + raise ForbiddenStructureException() # Generate the product structures by applying the forward reaction recipe try: @@ -1146,14 +1146,49 @@ def __generateProductStructures(self, reactantStructures, maps, forward, **optio raise ForbiddenStructureException() if (struct.getNumberOfRadicalElectrons() > maxRadicals) and (len(struct.atoms) - H < 1): raise ForbiddenStructureException() - + + # Generate other possible electronic states + electronicStrucutresList1 = [] + electronicStrucutresList2 = [] + struct01 = productStructures[0] + atoms01 = struct01.getLabeledAtoms() + for label01 in atoms01: + spin1 = struct01.getLabeledAtom(label01).spinMultiplicity + while spin1 > 0.1: + struct1 = deepcopy(productStructures[0]) + struct1.getLabeledAtom(label01).setSpinMultiplicity(spin1) + struct1.updateAtomTypes() + electronicStrucutresList1.append(struct1) + spin1 = spin1 - 2 + + if len(productStructures) == 2: + struct02 = productStructures[1] + atoms02 = struct02.getLabeledAtoms() + for label02 in atoms02: + spin2 = struct02.getLabeledAtom(label02).spinMultiplicity + while spin2 > 0.1: + struct2 = deepcopy(productStructures[1]) + struct2.getLabeledAtom(label02).setSpinMultiplicity(spin2) + struct2.updateAtomTypes() + electronicStrucutresList2.append(struct2) + spin2 = spin2 - 2 + # Check that product structures are allowed in this family # If not, then stop + forbiddenStructure = False for struct in productStructures: - struct.updateAtomTypes() - if self.isMoleculeForbidden(struct): raise ForbiddenStructureException() - - return productStructures + if self.isMoleculeForbidden(struct): + forbiddenStructure = True + + if not forbiddenStructure: + productStructuresList.append(productStructures) + + for struct1 in electronicStrucutresList1: + for struct2 in electronicStrucutresList2: + if not self.isMoleculeForbidden(struct1) and not self.isMoleculeForbidden(struct2): + productStructuresList.append([struct1,struct2]) + + return productStructuresList def isMoleculeForbidden(self, molecule): """ @@ -1179,9 +1214,9 @@ def __createReaction(self, reactants, products, isForward): if reactants[0].isIsomorphic(products[0]): return None elif len(reactants) == len(products) == 2: - if reactants[0].isIsomorphic(products[0]) and reactants[1].isIsomorphic(products[1]) and reactants[0].getFormula() != ('CH2' or 'CH') and reactants[1].getFormula() != ('CH2' or 'CH'): + if reactants[0].isIsomorphic(products[0]) and reactants[1].isIsomorphic(products[1]): return None - elif reactants[0].isIsomorphic(products[1]) and reactants[1].isIsomorphic(products[0]) and reactants[0].getFormula() != ('CH2' or 'CH') and reactants[1].getFormula() != ('CH2' or 'CH'): + elif reactants[0].isIsomorphic(products[1]) and reactants[1].isIsomorphic(products[0]): return None # Create and return template reaction object @@ -1313,13 +1348,14 @@ def __generateReactions(self, reactants, products=None, forward=True, **options) for map in mappings: reactantStructures = [molecule] try: - productStructures = self.__generateProductStructures(reactantStructures, [map], forward, **options) + productStructuresList = self.__generateProductStructures(reactantStructures, [map], forward, **options) except ForbiddenStructureException: pass else: - if productStructures is not None: - rxn = self.__createReaction(reactantStructures, productStructures, forward) - if rxn: rxnList.append(rxn) + if productStructuresList is not None: + for productStructures in productStructuresList: + rxn = self.__createReaction(reactantStructures, productStructures, forward) + if rxn: rxnList.append(rxn) # Bimolecular reactants: A + B --> products elif len(reactants) == 2 and len(template.reactants) == 2: @@ -1340,13 +1376,14 @@ def __generateReactions(self, reactants, products=None, forward=True, **options) for mapB in mappingsB: reactantStructures = [moleculeA, moleculeB] try: - productStructures = self.__generateProductStructures(reactantStructures, [mapA, mapB], forward, **options) + productStructuresList = self.__generateProductStructures(reactantStructures, [mapA, mapB], forward, **options) except ForbiddenStructureException: pass else: - if productStructures is not None: - rxn = self.__createReaction(reactantStructures, productStructures, forward) - if rxn: rxnList.append(rxn) + if productStructuresList is not None: + for productStructures in productStructuresList: + rxn = self.__createReaction(reactantStructures, productStructures, forward) + if rxn: rxnList.append(rxn) # Only check for swapped reactants if they are different if reactants[0] is not reactants[1]: @@ -1360,13 +1397,14 @@ def __generateReactions(self, reactants, products=None, forward=True, **options) for mapB in mappingsB: reactantStructures = [moleculeA, moleculeB] try: - productStructures = self.__generateProductStructures(reactantStructures, [mapA, mapB], forward, **options) + productStructuresList = self.__generateProductStructures(reactantStructures, [mapA, mapB], forward, **options) except ForbiddenStructureException: pass else: - if productStructures is not None: - rxn = self.__createReaction(reactantStructures, productStructures, forward) - if rxn: rxnList.append(rxn) + if productStructuresList is not None: + for productStructures in productStructuresList: + rxn = self.__createReaction(reactantStructures, productStructures, forward) + if rxn: rxnList.append(rxn) # If products is given, remove reactions from the reaction list that # don't generate the given products @@ -1380,14 +1418,6 @@ def __generateReactions(self, reactants, products=None, forward=True, **options) for reaction in rxnList0: products0 = reaction.products if forward else reaction.reactants - - # If the forward reaction contains CH(S) as reactant then correct product of reverse reaction from CH2(T) to CH2(S) - for product in products: - for molecule in product: - if molecule.isBiradicalSinglet() and (molecule.getFormula() == 'CH2' or molecule.getFormula() == 'NH'): - for product0 in products0: - if product0.isBiradicalTriplet() and (product0.getFormula() == 'CH2' or product0.getFormula() == 'NH'): - product0.changeTripletSinglet() # Skip reactions that don't match the given products match = False @@ -1408,7 +1438,7 @@ def __generateReactions(self, reactants, products=None, forward=True, **options) break if match: - rxnList.append(reaction) + rxnList.append(reaction) # The reaction list may contain duplicates of the same reaction # These duplicates should be combined (by increasing the degeneracy of From 9ea10088f2c7c3db691cba03a7ff881aba82f0b7 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Tue, 26 Nov 2013 00:59:11 -0500 Subject: [PATCH 56/64] Improve extension to product structure generation - improvements to obtain correct reaction degeneracy - replace while loop with if conditions --- rmgpy/data/kinetics/family.py | 193 ++++++++++++++++++++++++++++------ 1 file changed, 160 insertions(+), 33 deletions(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index 045a497a16..f312e6279b 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1097,7 +1097,7 @@ def __generateProductStructures(self, reactantStructures, maps, forward, **optio # Check that reactant structures are allowed in this family # If not, then stop for struct in reactantStructures: - if self.isMoleculeForbidden(struct) and forward: + if self.isMoleculeForbidden(struct): raise ForbiddenStructureException() # Generate the product structures by applying the forward reaction recipe @@ -1150,43 +1150,170 @@ def __generateProductStructures(self, reactantStructures, maps, forward, **optio # Generate other possible electronic states electronicStrucutresList1 = [] electronicStrucutresList2 = [] - struct01 = productStructures[0] - atoms01 = struct01.getLabeledAtoms() - for label01 in atoms01: - spin1 = struct01.getLabeledAtom(label01).spinMultiplicity - while spin1 > 0.1: - struct1 = deepcopy(productStructures[0]) - struct1.getLabeledAtom(label01).setSpinMultiplicity(spin1) - struct1.updateAtomTypes() - electronicStrucutresList1.append(struct1) - spin1 = spin1 - 2 + struct1 = productStructures[0] + struct1a = struct1.copy(True) + struct1a.updateAtomTypes() + electronicStrucutresList1.append(struct1a) + atoms1 = struct1.getRadicalAtoms() + + for atom1 in atoms1: + + radical1 = atom1.radicalElectrons + spin1 = atom1.spinMultiplicity + + if radical1 > 1 and radical1 < 4: + + if radical1 == 2 and spin1 == 3: + atom1.setSpinMultiplicity(1) + struct1a = struct1.copy(True) + struct1a.updateAtomTypes() + elif radical1 == 2 and spin1 == 1: + atom1.setSpinMultiplicity(3) + struct1a = struct1.copy(True) + struct1a.updateAtomTypes() + elif radical1 == 3 and spin1 == 4: + atom1.setSpinMultiplicity(2) + struct1a = struct1.copy(True) + struct1a.updateAtomTypes() + elif radical1 == 3 and spin1 == 2: + atom1.setSpinMultiplicity(4) + struct1a = struct1.copy(True) + struct1a.updateAtomTypes() + + for electronicStrucutres in electronicStrucutresList1: + if electronicStrucutres.isIsomorphic(struct1a): + break + else: + electronicStrucutresList1.append(struct1a) + + elif radical1 == 4: + + if spin1 == 5: + atom1.setSpinMultiplicity(3) + struct1a = struct1.copy(True) + struct1a.updateAtomTypes() + + atom1.setSpinMultiplicity(1) + struct1b = struct1.copy(True) + struct1b.updateAtomTypes() + elif spin1 == 3: + atom1.setSpinMultiplicity(5) + struct1a = struct1.copy(True) + struct1a.updateAtomTypes() + + atom1.setSpinMultiplicity(1) + struct1b = struct1.copy(True) + struct1b.updateAtomTypes() + elif spin1 == 1: + atom1.setSpinMultiplicity(5) + struct1a = struct1.copy(True) + struct1a.updateAtomTypes() + + atom1.setSpinMultiplicity(3) + struct1b = struct1.copy(True) + struct1b.updateAtomTypes() + + for electronicStrucutres in electronicStrucutresList1: + if electronicStrucutres.isIsomorphic(struct1a): + break + else: + electronicStrucutresList1.append(struct1a) + + for electronicStrucutres in electronicStrucutresList1: + if electronicStrucutres.isIsomorphic(struct1b): + break + else: + electronicStrucutresList1.append(struct1b) + if len(productStructures) == 2: - struct02 = productStructures[1] - atoms02 = struct02.getLabeledAtoms() - for label02 in atoms02: - spin2 = struct02.getLabeledAtom(label02).spinMultiplicity - while spin2 > 0.1: - struct2 = deepcopy(productStructures[1]) - struct2.getLabeledAtom(label02).setSpinMultiplicity(spin2) - struct2.updateAtomTypes() - electronicStrucutresList2.append(struct2) - spin2 = spin2 - 2 - # Check that product structures are allowed in this family - # If not, then stop - forbiddenStructure = False - for struct in productStructures: - if self.isMoleculeForbidden(struct): - forbiddenStructure = True + struct2 = productStructures[1] + struct2a = struct2.copy(True) + struct2a.updateAtomTypes() + electronicStrucutresList2.append(struct2a) + atoms2 = struct2.getRadicalAtoms() - if not forbiddenStructure: - productStructuresList.append(productStructures) + for atom2 in atoms2: + + radical2 = atom2.radicalElectrons + spin2 = atom2.spinMultiplicity + + if radical2 > 1 and radical2 < 4: + + if radical2 == 2 and spin2 == 3: + atom2.setSpinMultiplicity(1) + struct2a = struct2.copy(True) + struct2a.updateAtomTypes() + elif radical2 == 2 and spin2 == 1: + atom2.setSpinMultiplicity(3) + struct2a = struct2.copy(True) + struct2a.updateAtomTypes() + elif radical2 == 3 and spin2 == 4: + atom2.setSpinMultiplicity(2) + struct2a = struct2.copy(True) + struct2a.updateAtomTypes() + elif radical2 == 3 and spin2 == 2: + atom2.setSpinMultiplicity(4) + struct2a = struct2.copy(True) + struct2a.updateAtomTypes() + + for electronicStrucutres in electronicStrucutresList2: + if electronicStrucutres.isIsomorphic(struct2a): + break + else: + electronicStrucutresList2.append(struct2a) + + elif radical2 == 4: + + if spin2 == 5: + atom2.setSpinMultiplicity(3) + struct2a = struct2.copy(True) + struct2a.updateAtomTypes() + + atom2.setSpinMultiplicity(1) + struct2b = struct2.copy(True) + struct2b.updateAtomTypes() + elif spin2 == 3: + atom2.setSpinMultiplicity(5) + struct2a = struct2.copy(True) + struct2a.updateAtomTypes() + + atom2.setSpinMultiplicity(1) + struct2b = struct2.copy(True) + struct2b.updateAtomTypes() + elif spin2 == 1: + atom2.setSpinMultiplicity(5) + struct2a = struct2.copy(True) + struct2a.updateAtomTypes() + + atom2.setSpinMultiplicity(3) + struct2b = struct2.copy(True) + struct2b.updateAtomTypes() + + for electronicStrucutres in electronicStrucutresList2: + if electronicStrucutres.isIsomorphic(struct2a): + break + else: + electronicStrucutresList2.append(struct2a) + + for electronicStrucutres in electronicStrucutresList2: + if electronicStrucutres.isIsomorphic(struct2b): + break + else: + electronicStrucutresList2.append(struct2b) + + if len(productStructures) == 2: + + for structa in electronicStrucutresList1: + for structb in electronicStrucutresList2: + if not (self.isMoleculeForbidden(structa) or self.isMoleculeForbidden(structb)): + productStructuresList.append([structa,structb]) + elif len(productStructures) == 1: - for struct1 in electronicStrucutresList1: - for struct2 in electronicStrucutresList2: - if not self.isMoleculeForbidden(struct1) and not self.isMoleculeForbidden(struct2): - productStructuresList.append([struct1,struct2]) + for structa in electronicStrucutresList1: + if not (self.isMoleculeForbidden(structa)): + productStructuresList.append([structa]) return productStructuresList From b800aa9cdb6401e377d3cc3b068eacff839913c6 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Tue, 26 Nov 2013 01:01:48 -0500 Subject: [PATCH 57/64] Add function getRadicalAtoms to class Molecule This function returns a list of atoms of the molecules containing at least one unpaired electron. --- rmgpy/molecule/molecule.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index bc62ed285e..1085116e75 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1778,4 +1778,14 @@ def changeTripletSinglet(self): if atom.radicalElectrons == 2: if atom.spinMultiplicity == 3: atom.spinMultiplicity = 1 + + def getRadicalAtoms(self): + """ + Return the atoms in the molecule that have unpaired electrons. + """ + radicalAtomsList = [] + for atom in self.vertices: + if atom.radicalElectrons > 0: + radicalAtomsList.append(atom) + return radicalAtomsList From 76d9884ed1c126b6bb9bd1db1c38c74623c85215 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 27 Nov 2013 21:46:32 -0500 Subject: [PATCH 58/64] Correct the assertion of maximum number of radicals --- rmgpy/data/kinetics/family.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/data/kinetics/family.py b/rmgpy/data/kinetics/family.py index f312e6279b..3a5633d9d0 100644 --- a/rmgpy/data/kinetics/family.py +++ b/rmgpy/data/kinetics/family.py @@ -1144,7 +1144,7 @@ def __generateProductStructures(self, reactantStructures, maps, forward, **optio raise ForbiddenStructureException() if len(struct.atoms) - H > maxHeavyAtoms: raise ForbiddenStructureException() - if (struct.getNumberOfRadicalElectrons() > maxRadicals) and (len(struct.atoms) - H < 1): + if (struct.getNumberOfRadicalElectrons() > maxRadicals) and (len(struct.atoms) - H > 1): raise ForbiddenStructureException() # Generate other possible electronic states From 90b6b4eaa677d2f922c90115a94223b04746327e Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Wed, 27 Nov 2013 23:58:59 -0500 Subject: [PATCH 59/64] Improving stability of resonance structure generation - include the requirement for an existing lone electron pair for N5dd-N5ts resonance --- rmgpy/molecule/molecule.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 1085116e75..2f76db280b 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1720,7 +1720,7 @@ def findAllDelocalizationPathsN5dd_N5ts(self, atom1): for atom3, bond13 in atom1.edges.items(): index_atom_3 = index_atom_3 + 1 # Only double bonds are considered, at the moment we only consider non-radical nitrogen and oxygen atoms - if (bond13.isDouble() and atom3.radicalElectrons == 0 and not atom3.isOxygen() and not atom3.isCarbon() and (index_atom_2 != index_atom_3)): + if (bond13.isDouble() and atom3.radicalElectrons == 0 and atom3.lonePairs > 0 and not atom3.isOxygen() and not atom3.isCarbon() and (index_atom_2 != index_atom_3)): paths.append([atom1, atom2, atom3, bond12, bond13, 1]) for atom2, bond12 in atom1.edges.items(): From aa61f2700e001f7a73fd86bd29770426e5b813ad Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Thu, 28 Nov 2013 20:50:33 -0500 Subject: [PATCH 60/64] limit allyl resonance structure generation to monoradicals --- rmgpy/molecule/molecule.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/molecule/molecule.py b/rmgpy/molecule/molecule.py index 2f76db280b..d541bf067f 100644 --- a/rmgpy/molecule/molecule.py +++ b/rmgpy/molecule/molecule.py @@ -1662,7 +1662,7 @@ def findAllDelocalizationPaths(self, atom1): paths = [] for atom2, bond12 in atom1.edges.items(): # Vinyl bond must be capable of gaining an order - if bond12.isSingle() or bond12.isDouble(): + if (bond12.isSingle() or bond12.isDouble()) and atom1.radicalElectrons == 1: for atom3, bond23 in atom2.edges.items(): # Allyl bond must be capable of losing an order without breaking if atom1 is not atom3 and (bond23.isDouble() or bond23.isTriple()): From d71a1a068ecc558d4bfa42a66d1bb7ffbf4e4480 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Mon, 2 Dec 2013 10:34:52 -0500 Subject: [PATCH 61/64] Give the edge CHEMKIN files the same numbering as for the core ones - this prevents two edge files having the same number of species from overwriting each other --- rmgpy/rmg/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 5276e74406..432b4cd3f1 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -696,7 +696,7 @@ def saveChemkinFileEdge(self): Save the current reaction edge to a Chemkin file. """ logging.info('Saving current edge to Chemkin file...') - this_chemkin_path = os.path.join(self.outputDirectory, 'chemkin', 'chem_edge%04i.inp' % len(self.reactionModel.edge.species)) + this_chemkin_path = os.path.join(self.outputDirectory, 'chemkin', 'chem_edge%04i.inp' % len(self.reactionModel.core.species)) latest_chemkin_path = os.path.join(self.outputDirectory, 'chemkin','chem_edge.inp') latest_chemkin_verbose_path = os.path.join(self.outputDirectory, 'chemkin', 'chem_edge_annotated.inp') latest_dictionary_path = os.path.join(self.outputDirectory, 'chemkin','species_edge_dictionary.txt') From 05c51b835132b729f991ca8d88f9f74782dc7b6f Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Sat, 7 Dec 2013 11:11:59 -0500 Subject: [PATCH 62/64] Update the reaction families used by the examples - not using Substitution_O for now --- examples/rmg/1,3-hexadiene/input.py | 2 +- examples/rmg/TEOS/input.py | 2 +- examples/rmg/c3h4/input.py | 6 +++--- examples/rmg/diesel/input.py | 2 +- examples/rmg/e85/input.py | 2 +- examples/rmg/liquid_phase/input.py | 2 +- examples/rmg/methylformate/input.py | 2 +- examples/rmg/minimal/input.py | 2 +- examples/rmg/minimal_sensitivity/input.py | 2 +- 9 files changed, 11 insertions(+), 11 deletions(-) diff --git a/examples/rmg/1,3-hexadiene/input.py b/examples/rmg/1,3-hexadiene/input.py index ff5955127f..b6dbf0ca14 100644 --- a/examples/rmg/1,3-hexadiene/input.py +++ b/examples/rmg/1,3-hexadiene/input.py @@ -4,7 +4,7 @@ reactionLibraries = [], seedMechanisms = [], kineticsDepositories = ['training'], # 'all', 'default'==['training'], [], - kineticsFamilies = ['!Intra_Disproportionation'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], kineticsEstimator = 'rate rules', ) diff --git a/examples/rmg/TEOS/input.py b/examples/rmg/TEOS/input.py index 0677bf7685..44fecc5b47 100644 --- a/examples/rmg/TEOS/input.py +++ b/examples/rmg/TEOS/input.py @@ -4,7 +4,7 @@ reactionLibraries = [], seedMechanisms = [], kineticsDepositories = ['training'], - kineticsFamilies = ['!Intra_Disproportionation'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], kineticsEstimator = 'rate rules', ) diff --git a/examples/rmg/c3h4/input.py b/examples/rmg/c3h4/input.py index 1af1e146fd..ffc50b4d3b 100644 --- a/examples/rmg/c3h4/input.py +++ b/examples/rmg/c3h4/input.py @@ -1,10 +1,10 @@ # Data sources database( - thermoLibraries = ['primaryThermoLibrary', 'GRI-Mech3.0'], + thermoLibraries = ['primaryThermoLibrary', 'GRI-Mech3.0-N'], reactionLibraries = [], - seedMechanisms = [], + seedMechanisms = ['GRI-Mech3.0-N'], kineticsDepositories = ['training'], # 'all', 'default'==['training'], [], - kineticsFamilies = ['!Intra_Disproportionation'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], kineticsEstimator = 'rate rules', ) diff --git a/examples/rmg/diesel/input.py b/examples/rmg/diesel/input.py index 43389bc19f..1b3a010db4 100644 --- a/examples/rmg/diesel/input.py +++ b/examples/rmg/diesel/input.py @@ -4,7 +4,7 @@ reactionLibraries = [], seedMechanisms = [], kineticsDepositories = ['training'], # 'all', 'default'==['training'], [], - kineticsFamilies = ['!Intra_Disproportionation'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], kineticsEstimator = 'rate rules', ) diff --git a/examples/rmg/e85/input.py b/examples/rmg/e85/input.py index 95371e2951..6ffb04c275 100644 --- a/examples/rmg/e85/input.py +++ b/examples/rmg/e85/input.py @@ -4,7 +4,7 @@ reactionLibraries = [], seedMechanisms = ['GRI-Mech3.0'], kineticsDepositories = ['training'], - kineticsFamilies = ['!Intra_Disproportionation'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], kineticsEstimator = 'rate rules', ) diff --git a/examples/rmg/liquid_phase/input.py b/examples/rmg/liquid_phase/input.py index 236771e8a4..6e9f22839b 100644 --- a/examples/rmg/liquid_phase/input.py +++ b/examples/rmg/liquid_phase/input.py @@ -4,7 +4,7 @@ reactionLibraries = [], seedMechanisms = [], kineticsDepositories = ['training'], - kineticsFamilies = ['!Intra_Disproportionation'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], kineticsEstimator = 'rate rules', ) diff --git a/examples/rmg/methylformate/input.py b/examples/rmg/methylformate/input.py index 9ffd0bdb38..ccd3989236 100644 --- a/examples/rmg/methylformate/input.py +++ b/examples/rmg/methylformate/input.py @@ -4,7 +4,7 @@ reactionLibraries = [('Methylformate',False),('Glarborg/highP',False)], seedMechanisms = ['Glarborg/C2'], kineticsDepositories = ['training'], - kineticsFamilies = ['!Intra_Disproportionation'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], kineticsEstimator = 'rate rules', ) diff --git a/examples/rmg/minimal/input.py b/examples/rmg/minimal/input.py index 10839f4cba..5e5065847b 100644 --- a/examples/rmg/minimal/input.py +++ b/examples/rmg/minimal/input.py @@ -4,7 +4,7 @@ reactionLibraries = [], seedMechanisms = [], kineticsDepositories = ['training'], - kineticsFamilies = ['!Intra_Disproportionation'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], kineticsEstimator = 'rate rules', ) diff --git a/examples/rmg/minimal_sensitivity/input.py b/examples/rmg/minimal_sensitivity/input.py index bd62ddf77a..c1c57f0936 100644 --- a/examples/rmg/minimal_sensitivity/input.py +++ b/examples/rmg/minimal_sensitivity/input.py @@ -4,7 +4,7 @@ reactionLibraries = [], seedMechanisms = [], kineticsDepositories = ['training'], - kineticsFamilies = ['!Intra_Disproportionation'], + kineticsFamilies = ['!Intra_Disproportionation','!Substitution_O'], kineticsEstimator = 'rate rules', ) From 13297e344348e19d299bb1962959c9bd4c0cd434 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Sat, 7 Dec 2013 11:25:08 -0500 Subject: [PATCH 63/64] Add missing function call to calculate transport properties for seed --- rmgpy/rmg/main.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/rmgpy/rmg/main.py b/rmgpy/rmg/main.py index 00e6b15ef6..1d2d1ff6ce 100644 --- a/rmgpy/rmg/main.py +++ b/rmgpy/rmg/main.py @@ -335,6 +335,9 @@ def initialize(self, args): # DON'T generate any more reactions for the seed species at this time for seedMechanism in self.seedMechanisms: self.reactionModel.addSeedMechanismToCore(seedMechanism, react=False) + + for spec in self.reactionModel.core.species: + spec.generateTransportData(self.database) # Reaction libraries: add species and reactions from reaction library to the edge so # that RMG can find them if their rates are large enough From 32cb9301a357890b07153f37e835392bf6be79b0 Mon Sep 17 00:00:00 2001 From: Beat Buesser Date: Sat, 7 Dec 2013 11:28:58 -0500 Subject: [PATCH 64/64] Prevent RMG from crashing because of missing transport groups - add function getTransportPropertiesViaLennardJonesParameters(self,species) to transport.py - use this function as a fall back only if the library and the transport groups do not have proper information --- rmgpy/data/transport.py | 59 ++++++++++++++++++++++++++++++++++------- 1 file changed, 49 insertions(+), 10 deletions(-) diff --git a/rmgpy/data/transport.py b/rmgpy/data/transport.py index 46baa3b746..56dda0f604 100644 --- a/rmgpy/data/transport.py +++ b/rmgpy/data/transport.py @@ -274,10 +274,12 @@ def getTransportProperties(self, species): transport[0].comment = label break else: - #Transport not found in any loaded libraries, so estimate - transport = self.getTransportPropertiesViaGroupEstimates(species) - #data, library, entry = transport - #return data + try: + #Transport not found in any loaded libraries, so estimate + transport = self.getTransportPropertiesViaGroupEstimates(species) + except: + transport = self.getTransportPropertiesViaLennardJonesParameters(species) + return transport def getAllTransportProperties(self, species): @@ -368,7 +370,6 @@ def estimateCriticalPropertiesViaGroupAdditivity(self, molecule): # For transport estimation we need the atoms to already be sorted because we # iterate over them; if the order changes during the iteration then we # will probably not visit the right atoms, and so will get the transport wrong - molecule.sortVertices() if sum([atom.radicalElectrons for atom in molecule.atoms]) > 0: # radical species @@ -433,12 +434,9 @@ def estimateCriticalPropertiesViaGroupAdditivity(self, molecule): if molecule.isVertexInCycle(atom): self.__addCriticalPointContribution(groupData, self.groups['ring'], molecule, {'*':atom}) else: - self.__addCriticalPointContribution(groupData, self.groups['nonring'], molecule, {'*':atom}) + self.__addCriticalPointContribution(groupData, self.groups['nonring'], molecule, {'*':atom}) except KeyError: - logging.error("Couldn't find in any transport database:") - logging.error(molecule) - logging.error(molecule.toAdjacencyList()) - raise + raise Tb = 198.18 + groupData.Tb Vc = 17.5 + groupData.Vc @@ -493,6 +491,47 @@ def __addCriticalPointContribution(self, groupData, database, molecule, atom): groupData.structureIndex += data.structureIndex return groupData + + def getTransportPropertiesViaLennardJonesParameters(self,species): + """ + Serves as last resort if every other method to estimate Transport Properties fails. + + Generate the Lennard-Jones parameters for the species. + """ + count = sum([1 for atom in species.molecule[0].vertices if atom.isNonHydrogen()]) + + if count == 1: + sigma = (3.758e-10,"m") + epsilon = (148.6,"K") + elif count == 2: + sigma = (4.443e-10,"m") + epsilon = (110.7,"K") + elif count == 3: + sigma = (5.118e-10,"m") + epsilon = (237.1,"K") + elif count == 4: + sigma = (4.687e-10,"m") + epsilon = (531.4,"K") + elif count == 5: + sigma = (5.784e-10,"m") + epsilon = (341.1,"K") + else: + sigma = (5.949e-10,"m") + epsilon = (399.3,"K") + + shapeIndex = 1 if species.molecule[0].isLinear() else 2 + + transport = TransportData( + shapeIndex = shapeIndex, # 1 if linear, else 2 + epsilon = epsilon, + sigma = sigma, + dipoleMoment = (0, 'C*m'), + polarizability = (0, 'angstroms^3'), + rotrelaxcollnum = 0, # rotational relaxation collision number at 298 K + comment = 'Epsilon & sigma estimated with fixed Lennard Jones Parameters. This is the fallback method! Try improving transport databases!' + ) + + return (transport, None, None) class CriticalPoint: """