Skip to content

Commit

Permalink
Merge pull request #1600 from ReactionMechanismGenerator/resbug
Browse files Browse the repository at this point in the history 
Fix resonance structure generation and filtration bugs
  • Loading branch information
@alongd
alongd committed May 15, 2019
2 parents baed747 + a97e70b commit fd2f90d
Show file tree
Hide file tree
Showing 3 changed files with 119 additions and 103 deletions.
10 changes: 9 additions & 1 deletion rmgpy/molecule/molecule.py
View file
Expand Up @@ -1085,9 +1085,17 @@ def copy(self, deep=False):
If `deep` is ``False`` or not specified, a shallow copy is made: the
original vertices and edges are used in the new graph.
"""
other = cython.declare(Molecule)
cython.declare(g=Graph, other=Molecule, i=int, v1=Vertex, v2=Vertex)
g = Graph.copy(self, deep)
other = Molecule(g.vertices)
# Copy connectivity values and sorting labels
for i in xrange(len(self.vertices)):
v1 = self.vertices[i]
v2 = other.vertices[i]
v2.connectivity1 = v1.connectivity1
v2.connectivity2 = v1.connectivity2
v2.connectivity3 = v1.connectivity3
v2.sortingLabel = v1.sortingLabel
other.multiplicity = self.multiplicity
other.reactive = self.reactive
return other
Expand Down
125 changes: 28 additions & 97 deletions rmgpy/molecule/resonance.py
View file
Expand Up @@ -331,15 +331,6 @@ def generate_allyl_delocalization_resonance_structures(mol):
bond23.decrementOrder()
# Make a copy of structure
structure = mol.copy(deep=True)
# Also copy the connectivity values, since they are the same
# for all resonance structures
for index in xrange(len(mol.vertices)):
v1 = mol.vertices[index]
v2 = structure.vertices[index]
v2.connectivity1 = v1.connectivity1
v2.connectivity2 = v1.connectivity2
v2.connectivity3 = v1.connectivity3
v2.sortingLabel = v1.sortingLabel
# Restore current structure
atom1.incrementRadical()
atom3.decrementRadical()
Expand Down Expand Up @@ -378,15 +369,6 @@ def generate_lone_pair_multiple_bond_resonance_structures(mol):
atom3.updateCharge()
# Make a copy of structure
structure = mol.copy(deep=True)
# Also copy the connectivity values, since they are the same
# for all resonance structures
for index in xrange(len(mol.vertices)):
v1 = mol.vertices[index]
v2 = structure.vertices[index]
v2.connectivity1 = v1.connectivity1
v2.connectivity2 = v1.connectivity2
v2.connectivity3 = v1.connectivity3
v2.sortingLabel = v1.sortingLabel
# Restore current structure
atom1.incrementLonePairs()
atom3.decrementLonePairs()
Expand Down Expand Up @@ -428,15 +410,6 @@ def generate_adj_lone_pair_radical_resonance_structures(mol):
atom2.updateCharge()
# Make a copy of structure
structure = mol.copy(deep=True)
# Also copy the connectivity values, since they are the same
# for all resonance structures
for index in xrange(len(mol.vertices)):
v1 = mol.vertices[index]
v2 = structure.vertices[index]
v2.connectivity1 = v1.connectivity1
v2.connectivity2 = v1.connectivity2
v2.connectivity3 = v1.connectivity3
v2.sortingLabel = v1.sortingLabel
# Restore current structure
atom1.incrementRadical()
atom1.decrementLonePairs()
Expand Down Expand Up @@ -479,15 +452,6 @@ def generate_adj_lone_pair_multiple_bond_resonance_structures(mol):
atom2.updateCharge()
# Make a copy of structure
structure = mol.copy(deep=True)
# Also copy the connectivity values, since they are the same
# for all resonance structures
for index in xrange(len(mol.vertices)):
v1 = mol.vertices[index]
v2 = structure.vertices[index]
v2.connectivity1 = v1.connectivity1
v2.connectivity2 = v1.connectivity2
v2.connectivity3 = v1.connectivity3
v2.sortingLabel = v1.sortingLabel
# Restore current structure
if direction == 1: # The direction <increasing> the bond order
atom1.incrementLonePairs()
Expand Down Expand Up @@ -541,15 +505,6 @@ def generate_adj_lone_pair_radical_multiple_bond_resonance_structures(mol):
atom2.updateCharge()
# Make a copy of structure
structure = mol.copy(deep=True)
# Also copy the connectivity values, since they are the same
# for all resonance structures
for index in xrange(len(mol.vertices)):
v1 = mol.vertices[index]
v2 = structure.vertices[index]
v2.connectivity1 = v1.connectivity1
v2.connectivity2 = v1.connectivity2
v2.connectivity3 = v1.connectivity3
v2.sortingLabel = v1.sortingLabel
# Restore current structure
if direction == 1: # The direction <increasing> the bond order
atom1.incrementLonePairs()
Expand Down Expand Up @@ -593,15 +548,6 @@ def generate_N5dc_radical_resonance_structures(mol):
atom3.updateCharge()
# Make a copy of structure
structure = mol.copy(deep=True)
# Also copy the connectivity values, since they are the same
# for all resonance structures
for index in xrange(len(mol.vertices)):
v1 = mol.vertices[index]
v2 = structure.vertices[index]
v2.connectivity1 = v1.connectivity1
v2.connectivity2 = v1.connectivity2
v2.connectivity3 = v1.connectivity3
v2.sortingLabel = v1.sortingLabel
# Restore current structure
atom2.incrementRadical()
atom2.decrementLonePairs()
Expand Down Expand Up @@ -637,14 +583,9 @@ def generate_optimal_aromatic_resonance_structures(mol, features=None):
if not features['isCyclic']:
return []

# Copy the molecule so we don't affect the original
molecule = mol.copy(deep=True)

# First get all rings in the molecule
rings = molecule.getAllSimpleCyclesOfSize(6)

# Then determine which ones are aromatic
aromatic_bonds = molecule.getAromaticRings(rings)[1]

# Attempt to rearrange electrons to obtain a structure with the most aromatic rings
# Possible rearrangements include aryne resonance and allyl resonance
res_list = [generate_aryne_resonance_structures]
Expand All @@ -658,38 +599,37 @@ def generate_optimal_aromatic_resonance_structures(mol, features=None):

_generate_resonance_structures(kekule_list, res_list)

if len(kekule_list) > 1:
# We found additional structures, so we need to evaluate all of them
max_num = 0
mol_list = []

# Iterate through the adjacent resonance structures and keep the structures with the most aromatic rings
for mol0 in kekule_list:
aromatic_bonds = mol0.getAromaticRings()[1]
if len(aromatic_bonds) > max_num:
max_num = len(aromatic_bonds)
mol_list = [(mol0, aromatic_bonds)]
elif len(aromatic_bonds) == max_num:
mol_list.append((mol0, aromatic_bonds))
else:
# Otherwise, it is not possible to increase the number of aromatic rings by moving electrons,
# so go ahead with the inputted form of the molecule
mol_list = [(molecule, aromatic_bonds)]
# Sort all of the generated structures by number of perceived aromatic rings
mol_dict = {}
for mol0 in kekule_list:
aromatic_bonds = mol0.getAromaticRings()[1]
num_aromatic = len(aromatic_bonds)
mol_dict.setdefault(num_aromatic, []).append((mol0, aromatic_bonds))

# List of potential number of aromatic rings, sorted from largest to smallest
arom_options = sorted(mol_dict.keys(), reverse=True)

new_mol_list = []
for num in arom_options:
mol_list = mol_dict[num]
# Generate the aromatic resonance structure(s)
for mol0, aromatic_bonds in mol_list:
# Aromatize the molecule in place
result = generate_aromatic_resonance_structure(mol0, aromatic_bonds, copy=False)
if not result:
# We failed to aromatize this molecule
# This could be due to incorrect aromaticity perception by RDKit
continue

# Generate the aromatic resonance structure(s)
for mol0, aromatic_bonds in mol_list:
# Aromatize the molecule in place
success = generate_aromatic_resonance_structure(mol0, aromatic_bonds, copy=False)
if not success:
continue
for mol1 in new_mol_list:
if mol1.isIsomorphic(mol0):
break
else:
new_mol_list.append(mol0)

for mol1 in new_mol_list:
if mol1.isIsomorphic(mol0):
break
else:
new_mol_list.append(mol0)
if new_mol_list:
# We found the most aromatic resonance structures so there's no need to try smaller numbers
break

return new_mol_list

Expand Down Expand Up @@ -824,15 +764,6 @@ def generate_aryne_resonance_structures(mol):
bond.setOrderStr(new_orders[i])
# Make a copy of the molecule
new_mol = mol.copy(deep=True)
# Also copy the connectivity values, since they are the same
# for all resonance structures
for i in xrange(len(mol.vertices)):
v1 = mol.vertices[i]
v2 = new_mol.vertices[i]
v2.connectivity1 = v1.connectivity1
v2.connectivity2 = v1.connectivity2
v2.connectivity3 = v1.connectivity3
v2.sortingLabel = v1.sortingLabel
# Undo the changes to the current molecule
for i, bond in enumerate(bond_list):
bond.setOrderStr(bond_orders[i])
Expand Down
87 changes: 82 additions & 5 deletions rmgpy/molecule/resonanceTest.py
View file
Expand Up @@ -258,13 +258,37 @@ def testAromaticWithLonePairResonance(self):
molList = generate_resonance_structures(Molecule(SMILES="c1ccccc1CC=N[O]"))
self.assertEqual(len(molList), 4)

@work_in_progress
def testAromaticWithNResonance(self):
"""Test resonance structure generation for aromatic species with N5ddc <=> N5tc resonance"""
# WIP: currently generate_N5dc_resonance_structures does not apply for aromatic structures
"""Test resonance structure generation for aromatic species with lone pair resonance"""
molList = generate_resonance_structures(Molecule(SMILES="c1ccccc1CCN=[N+]=[N-]"))
self.assertEqual(len(molList), 4)
# TODO: this test cannot be run because RDKit (which checks for aromaticity) cannot process hyper-valence N
self.assertEqual(len(molList), 2)

def testAromaticWithONResonance(self):
"""Test resonance structure generation for aromatic species with heteroatoms
This test was specifically designed to recreate RMG-Py issue #1598.
Key conditions: having heteroatoms, starting with aromatic structure, and keep_isomorphic=True
"""
molList = generate_resonance_structures(Molecule().fromAdjacencyList("""
multiplicity 2
1 O u0 p2 c0 {4,S} {16,S}
2 O u1 p2 c0 {4,S}
3 O u0 p3 c-1 {4,S}
4 N u0 p0 c+1 {1,S} {2,S} {3,S} {5,S}
5 C u0 p0 c0 {4,S} {6,B} {7,B}
6 C u0 p0 c0 {5,B} {8,B} {11,S}
7 C u0 p0 c0 {5,B} {10,B} {15,S}
8 C u0 p0 c0 {6,B} {9,B} {12,S}
9 C u0 p0 c0 {8,B} {10,B} {13,S}
10 C u0 p0 c0 {7,B} {9,B} {14,S}
11 H u0 p0 c0 {6,S}
12 H u0 p0 c0 {8,S}
13 H u0 p0 c0 {9,S}
14 H u0 p0 c0 {10,S}
15 H u0 p0 c0 {7,S}
16 H u0 p0 c0 {1,S}
"""), keep_isomorphic=True)
self.assertEqual(len(molList), 1)

def testNoClarStructures(self):
"""Test that we can turn off Clar structure generation."""
Expand Down Expand Up @@ -1141,6 +1165,59 @@ def testFalseNegativePolycylicAromaticityPerception2(self):
self.assertEqual(len(out), 7)
self.assertTrue(any([m.isIsomorphic(aromatic) for m in out]))

@work_in_progress
def testInconsistentAromaticStructureGeneration(self):
"""Test an unusual case of inconsistent aromaticity perception."""
mol1 = Molecule().fromAdjacencyList("""
multiplicity 2
1 C u0 p0 c0 {2,S} {6,S} {11,S} {12,S}
2 C u0 p0 c0 {1,S} {3,D} {4,S}
3 C u0 p0 c0 {2,D} {5,S} {7,S}
4 C u0 p0 c0 {2,S} {7,D} {10,S}
5 C u1 p0 c0 {3,S} {8,S} {9,S}
6 C u0 p0 c0 {1,S} {8,D} {13,S}
7 C u0 p0 c0 {3,S} {4,D} {17,S}
8 C u0 p0 c0 {5,S} {6,D} {14,S}
9 C u0 p0 c0 {5,S} {10,D} {15,S}
10 C u0 p0 c0 {4,S} {9,D} {16,S}
11 H u0 p0 c0 {1,S}
12 H u0 p0 c0 {1,S}
13 H u0 p0 c0 {6,S}
14 H u0 p0 c0 {8,S}
15 H u0 p0 c0 {9,S}
16 H u0 p0 c0 {10,S}
17 H u0 p0 c0 {7,S}
""")

mol2 = Molecule().fromAdjacencyList("""
multiplicity 2
1 C u0 p0 c0 {2,S} {6,S} {11,S} {12,S}
2 C u0 p0 c0 {1,S} {3,D} {4,S}
3 C u0 p0 c0 {2,D} {5,S} {7,S}
4 C u0 p0 c0 {2,S} {7,D} {9,S}
5 C u1 p0 c0 {3,S} {8,S} {10,S}
6 C u0 p0 c0 {1,S} {8,D} {13,S}
7 C u0 p0 c0 {3,S} {4,D} {16,S}
8 C u0 p0 c0 {5,S} {6,D} {15,S}
9 C u0 p0 c0 {4,S} {10,D} {17,S}
10 C u0 p0 c0 {5,S} {9,D} {14,S}
11 H u0 p0 c0 {1,S}
12 H u0 p0 c0 {1,S}
13 H u0 p0 c0 {6,S}
14 H u0 p0 c0 {10,S}
15 H u0 p0 c0 {8,S}
16 H u0 p0 c0 {7,S}
17 H u0 p0 c0 {9,S}
""")

# These two slightly different adjlists should be the same structure
self.assertTrue(mol1.isIsomorphic(mol2))

# However, they give different resonance structures
res1 = generate_resonance_structures(mol1)
res2 = generate_resonance_structures(mol2)
self.assertEqual(res1, res2)


class ClarTest(unittest.TestCase):
"""
Expand Down

0 comments on commit fd2f90d

Please sign in to comment.