Adds resonance for radicals adjacent to sulfur lone pairs

Author: rgillis8

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)

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)

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)
 

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.