Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix resonance structure generation and filtration bugs #1600

Merged
merged 4 commits into from May 15, 2019
Merged

Fix resonance structure generation and filtration bugs #1600

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7ce6b7d
Add failing unit test for interesting resonance bug
mliu49 May 13, 2019
fa588a2
Fix optimal arom res struct generation to try all structures
mliu49 May 13, 2019
8289d3c
Copy connectivity values and sorting labels in Molecule.copy
mliu49 May 14, 2019
a97e70b
Add unit test for heteroatom+aromatic resonance structure generation
mliu49 May 14, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
10 changes: 9 additions & 1 deletion rmgpy/molecule/molecule.py
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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