From a969ae4935529fe58254111f089092b2496f1662 Mon Sep 17 00:00:00 2001 From: rgillis8 Date: Thu, 15 Jun 2017 16:36:09 -0400 Subject: [PATCH] Adds resonance for radicals adjacent to sulfur lone pairs --- rmgpy/molecule/pathfinder.pxd | 4 +- rmgpy/molecule/pathfinder.py | 48 ++++++++++++++++-- rmgpy/molecule/resonance.pxd | 4 +- rmgpy/molecule/resonance.py | 93 +++++++++++++++++++++++++++++++++-- 4 files changed, 138 insertions(+), 11 deletions(-) diff --git a/rmgpy/molecule/pathfinder.pxd b/rmgpy/molecule/pathfinder.pxd index 9a816884b7a..2d1de34aea0 100644 --- a/rmgpy/molecule/pathfinder.pxd +++ b/rmgpy/molecule/pathfinder.pxd @@ -19,6 +19,8 @@ cpdef dict compute_atom_distance(list atom_indices, Molecule mol) cpdef list findAllDelocalizationPaths(Atom atom1) -cpdef list findAllDelocalizationPathsLonePairRadical(Atom atom1) +cpdef list findAllDelocalizationPathsLonePairRadicalCharge(Atom atom1) + +cpdef list findAllDelocalizationPathsLonePairRadicalDoubleBond(Atom atom1) cpdef list findAllDelocalizationPathsN5dd_N5ts(Atom atom1) diff --git a/rmgpy/molecule/pathfinder.py b/rmgpy/molecule/pathfinder.py index 87002ccbcf5..40e5a51167d 100644 --- a/rmgpy/molecule/pathfinder.py +++ b/rmgpy/molecule/pathfinder.py @@ -260,7 +260,7 @@ def findAllDelocalizationPaths(atom1): paths.append([atom1, atom2, atom3, bond12, bond23]) return paths -def findAllDelocalizationPathsLonePairRadical(atom1): +def findAllDelocalizationPathsLonePairRadicalCharge(atom1): """ Find all the delocalization paths of lone electron pairs next to the radical center indicated by `atom1`. Used to generate resonance isomers. @@ -271,9 +271,8 @@ def findAllDelocalizationPathsLonePairRadical(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())): + # In a first step we only consider nitrogen, sulfur and oxygen atoms as possible radical centers + if not ((atom1.lonePairs == 0 and atom1.isNitrogen()) or(atom1.lonePairs == 2 and atom1.isOxygen()) or (atom1.lonePairs <= 1 and atom1.isSulfur())): return [] # Find all delocalization paths @@ -282,11 +281,50 @@ def findAllDelocalizationPathsLonePairRadical(atom1): # Only single bonds are considered if bond12.isSingle(): # 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): + if ((atom2.lonePairs == 1 and atom2.isNitrogen()) or (atom2.lonePairs == 3 and atom2.isOxygen()) or (atom2.lonePairs == 2 and atom2.isSulfur())) and (atom2.radicalElectrons == 0): paths.append([atom1, atom2]) return paths +def findAllDelocalizationPathsLonePairRadicalDoubleBond(atom1): + """ + 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) + cython.declare(atom2=Atom, bond12=Bond) + + # No paths if atom1 is not a radical + if atom1.radicalElectrons <= 0: + return [] + if atom1.charge != 0: + return [] + + # Find all delocalization paths + #radical splits a lone pair on a sulfur atom + paths = [] + for atom2, bond12 in atom1.edges.items(): + if atom2.charge != 0: + continue + # Only single bonds are considered + if bond12.isSingle(): + if (atom2.lonePairs >= 1 and atom2.isSulfur()) and (atom2.radicalElectrons == 0): + paths.append([atom1, atom2, bond12, 1]) + + #sulfur radical resonates to form a lone pair and radical on the adjacent atom + #only consider sulfur radicals that have less than 2 lone pairs (i.e. can form another lone pair) + if not (atom1.lonePairs <= 1 and atom1.isSulfur()): + return paths + + #cycle through adjacent atoms until you find one that is not a radical itself and connected to atom1 by a double or triple bond + for atom2, bond12 in atom1.edges.items(): + if atom2.charge != 0: + continue + if bond12.isDouble() or bond12.isTriple(): + if atom2.radicalElectrons == 0: + paths.append([atom1, atom2, bond12, 2]) + return paths + def findAllDelocalizationPathsN5dd_N5ts(atom1): """ Find all the resonance structures of nitrogen atoms with two double bonds (N5dd) diff --git a/rmgpy/molecule/resonance.pxd b/rmgpy/molecule/resonance.pxd index 3bb51ba9a7e..a503d395043 100644 --- a/rmgpy/molecule/resonance.pxd +++ b/rmgpy/molecule/resonance.pxd @@ -11,7 +11,9 @@ cpdef list _generateResonanceStructures(list molList, list methodList, bint keep cpdef list generateAdjacentResonanceStructures(Molecule mol) -cpdef list generateLonePairRadicalResonanceStructures(Molecule mol) +cpdef list generateLonePairRadicalResonanceChargeStructures(Molecule mol) + +cpdef list generateLonePairRadicalResonanceDoubleBondStructures(Molecule mol) cpdef list generateN5dd_N5tsResonanceStructures(Molecule mol) diff --git a/rmgpy/molecule/resonance.py b/rmgpy/molecule/resonance.py index 1bc5881cabd..c5032b5875f 100644 --- a/rmgpy/molecule/resonance.py +++ b/rmgpy/molecule/resonance.py @@ -71,7 +71,8 @@ def populateResonanceAlgorithms(features=None): if features is None: methodList = [ generateAdjacentResonanceStructures, - generateLonePairRadicalResonanceStructures, + generateLonePairRadicalResonanceChargeStructures, + generateLonePairRadicalResonanceDoubleBondStructures, generateN5dd_N5tsResonanceStructures, generateAromaticResonanceStructures, generateKekuleStructure, @@ -87,7 +88,9 @@ def populateResonanceAlgorithms(features=None): if features['hasNitrogen']: methodList.append(generateN5dd_N5tsResonanceStructures) if features['hasLonePairs']: - methodList.append(generateLonePairRadicalResonanceStructures) + methodList.append(generateLonePairRadicalResonanceChargeStructures) + if features['hasLonePairs']: + methodList.append(generateLonePairRadicalResonanceDoubleBondStructures) return methodList @@ -106,6 +109,7 @@ def analyzeMolecule(mol): 'isArylRadical': False, 'hasNitrogen': False, 'hasOxygen': False, + 'hasSulfur': False, 'hasLonePairs': False, } @@ -122,6 +126,8 @@ def analyzeMolecule(mol): features['hasNitrogen'] = True if atom.isOxygen(): features['hasOxygen'] = True + if atom.isSulfur(): + features['hasSulfur'] = True if atom.lonePairs > 0: features['hasLonePairs'] = True @@ -301,7 +307,7 @@ def generateAdjacentResonanceStructures(mol): return isomers -def generateLonePairRadicalResonanceStructures(mol): +def generateLonePairRadicalResonanceChargeStructures(mol): """ Generate all of the resonance structures formed by lone electron pair - radical shifts. """ @@ -315,7 +321,7 @@ def generateLonePairRadicalResonanceStructures(mol): if mol.isRadical(): # Iterate over radicals in structure for atom in mol.vertices: - paths = pathfinder.findAllDelocalizationPathsLonePairRadical(atom) + paths = pathfinder.findAllDelocalizationPathsLonePairRadicalCharge(atom) for atom1, atom2 in paths: # Adjust to (potentially) new resonance isomer atom1.decrementRadical() @@ -348,6 +354,85 @@ def generateLonePairRadicalResonanceStructures(mol): return isomers +def generateLonePairRadicalResonanceDoubleBondStructures(mol): + """ + Generate all of the resonance structures 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) + cython.declare(v1=Vertex, v2=Vertex) + + isomers = [] + + # Radicals + if mol.isRadical(): + # Iterate over radicals in structure + for atom in mol.vertices: + paths = pathfinder.findAllDelocalizationPathsLonePairRadicalDoubleBond(atom) + for atom1, atom2, bond12, direction in paths: + #a radical adjacent to a sulfur atom forms a higher order bond and moves the radical electron to the sulfur + if direction == 1: + # Adjust to (potentially) new resonance isomer + atom1.decrementRadical() + atom2.incrementRadical() + atom2.decrementLonePairs() + bond12.incrementOrder() + atom1.updateCharge() + atom2.updateCharge() + # Make a copy of isomer + isomer = mol.copy(deep=True) + # Also copy the connectivity values, since they are the same + # for all resonance forms + for index in range(len(mol.vertices)): + v1 = mol.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() + atom2.decrementRadical() + atom2.incrementLonePairs() + bond12.decrementOrder() + atom1.updateCharge() + atom2.updateCharge() + # Append to isomer list if unique + isomer.updateAtomTypes(logSpecies=False) + isomers.append(isomer) + #now adsorb an electron from an adjacent double or triple bond to a sulfur into a lone pair + elif direction == 2: + # Adjust to (potentially) new resonance isomer + atom1.decrementRadical() + atom2.incrementRadical() + atom1.incrementLonePairs() + bond12.decrementOrder() + atom1.updateCharge() + atom2.updateCharge() + # Make a copy of isomer + isomer = mol.copy(deep=True) + # Also copy the connectivity values, since they are the same + # for all resonance forms + for index in range(len(mol.vertices)): + v1 = mol.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() + atom2.decrementRadical() + atom1.decrementLonePairs() + bond12.incrementOrder() + atom1.updateCharge() + atom2.updateCharge() + # Append to isomer list if unique + isomer.updateAtomTypes(logSpecies=False) + isomers.append(isomer) + return isomers + + def generateN5dd_N5tsResonanceStructures(mol): """ Generate all of the resonance structures formed by shifts between N5dd and N5ts.