In [1]:
from rmgpy.molecule.molecule import Molecule
from rmgpy.molecule.resonance import *
from IPython.display import display
import time
In [ ]:
struct1 = Molecule().fromAdjacencyList("""
multiplicity 2
1 C u0 p0 c0 {2,S} {3,D} {7,S}
2 C u0 p0 c0 {1,S} {4,S} {8,D}
3 C u0 p0 c0 {1,D} {5,S} {11,S}
4 C u0 p0 c0 {2,S} {9,D} {10,S}
5 C u0 p0 c0 {3,S} {6,D} {15,S}
6 C u0 p0 c0 {5,D} {12,S} {16,S}
7 C u0 p0 c0 {1,S} {12,D} {18,S}
8 C u0 p0 c0 {2,D} {13,S} {19,S}
9 C u0 p0 c0 {4,D} {14,S} {22,S}
10 C u0 p0 c0 {4,S} {11,D} {23,S}
11 C u0 p0 c0 {3,S} {10,D} {24,S}
12 C u0 p0 c0 {6,S} {7,D} {17,S}
13 C u0 p0 c0 {8,S} {14,D} {20,S}
14 C u0 p0 c0 {9,S} {13,D} {21,S}
15 C u1 p0 c0 {5,S} {25,S} {26,S}
16 H u0 p0 c0 {6,S}
17 H u0 p0 c0 {12,S}
18 H u0 p0 c0 {7,S}
19 H u0 p0 c0 {8,S}
20 H u0 p0 c0 {13,S}
21 H u0 p0 c0 {14,S}
22 H u0 p0 c0 {9,S}
23 H u0 p0 c0 {10,S}
24 H u0 p0 c0 {11,S}
25 H u0 p0 c0 {15,S}
26 H u0 p0 c0 {15,S}
""")
# Kekulized form, radical on ring
struct2 = Molecule().fromAdjacencyList("""
multiplicity 2
1 C u0 p0 c0 {2,S} {3,S} {7,D}
2 C u0 p0 c0 {1,S} {4,S} {8,D}
3 C u0 p0 c0 {1,S} {5,S} {11,D}
4 C u0 p0 c0 {2,S} {9,S} {10,D}
5 C u0 p0 c0 {3,S} {6,S} {15,D}
6 C u0 p0 c0 {5,S} {12,D} {16,S}
7 C u0 p0 c0 {1,D} {12,S} {17,S}
8 C u0 p0 c0 {2,D} {13,S} {18,S}
9 C u0 p0 c0 {4,S} {14,D} {19,S}
10 C u0 p0 c0 {4,D} {11,S} {20,S}
11 C u0 p0 c0 {3,D} {10,S} {21,S}
12 C u0 p0 c0 {6,D} {7,S} {22,S}
13 C u1 p0 c0 {8,S} {14,S} {23,S}
14 C u0 p0 c0 {9,D} {13,S} {24,S}
15 C u0 p0 c0 {5,D} {25,S} {26,S}
16 H u0 p0 c0 {6,S}
17 H u0 p0 c0 {7,S}
18 H u0 p0 c0 {8,S}
19 H u0 p0 c0 {9,S}
20 H u0 p0 c0 {10,S}
21 H u0 p0 c0 {11,S}
22 H u0 p0 c0 {12,S}
23 H u0 p0 c0 {13,S}
24 H u0 p0 c0 {14,S}
25 H u0 p0 c0 {15,S}
26 H u0 p0 c0 {15,S}
""")
# Aromatic form
struct3 = Molecule().fromAdjacencyList("""
multiplicity 2
1 C u0 p0 c0 {2,B} {3,B} {7,B}
2 C u0 p0 c0 {1,B} {4,B} {8,B}
3 C u0 p0 c0 {1,B} {5,B} {11,B}
4 C u0 p0 c0 {2,B} {9,B} {10,B}
5 C u0 p0 c0 {3,B} {6,B} {15,S}
6 C u0 p0 c0 {5,B} {12,B} {16,S}
7 C u0 p0 c0 {1,B} {12,B} {18,S}
8 C u0 p0 c0 {2,B} {13,B} {19,S}
9 C u0 p0 c0 {4,B} {14,B} {22,S}
10 C u0 p0 c0 {4,B} {11,B} {23,S}
11 C u0 p0 c0 {3,B} {10,B} {24,S}
12 C u0 p0 c0 {6,B} {7,B} {17,S}
13 C u0 p0 c0 {8,B} {14,B} {20,S}
14 C u0 p0 c0 {9,B} {13,B} {21,S}
15 C u1 p0 c0 {5,S} {25,S} {26,S}
16 H u0 p0 c0 {6,S}
17 H u0 p0 c0 {12,S}
18 H u0 p0 c0 {7,S}
19 H u0 p0 c0 {8,S}
20 H u0 p0 c0 {13,S}
21 H u0 p0 c0 {14,S}
22 H u0 p0 c0 {9,S}
23 H u0 p0 c0 {10,S}
24 H u0 p0 c0 {11,S}
25 H u0 p0 c0 {15,S}
26 H u0 p0 c0 {15,S}
""")
t0 = time.time()
out1 = generateAromaticResonanceStructures(struct1)
t1 = time.time()
print t1 - t0
t0 = time.time()
out2 = generateAromaticResonanceStructures(struct2)
t1 = time.time()
print t1 - t0
t0 = time.time()
out3 = generateAromaticResonanceStructures(struct3)
t1 = time.time()
print t1 - t0
In [ ]:
for o in out1:
display(o)
In [ ]:
for o in out2:
display(o)
In [ ]:
for o in out3:
display(o)
In [3]:
mol = Molecule(SMILES='c12ccccc1c(C=[CH])ccc2')
display(mol)
print '===================='
out = generateResonanceStructures(mol)
for o in out:
display(o)
In [ ]:
print out[1].toAdjacencyList()
In [ ]:
for mol in out:
print mol.toAdjacencyList()
print '\n'
In [2]:
mol = Molecule(SMILES="C1=CC2=CC=C3C=CC4=C5C6=C(C2=C35)C1=CC=C6C=C4")
display(mol)
print '===================='
out = generateClarStructures(mol)
print len(out)
for o in out:
display(o)
In [7]:
mol.getAromaticRings()[0]
Out[7]:
In [ ]:
print out[1].toAdjacencyList()
In [8]:
mol = Molecule(SMILES="C1=CC2=CC=CC3CC=CC(=C1)C=32")
display(mol)
print '===================='
out = generateClarStructures(mol)
for o in out:
display(o)
In [ ]:
toRDKitMol(mol, removeHs=False, returnMapping=True)
In [14]:
mol = Molecule(SMILES="c1c2cccc1C(=C)C=[C]2")
display(mol)
print '===================='
out = generateResonanceStructures(mol)
for o in out:
display(o)
In [13]:
[atom.props['inRing'] for atom in out[3].atoms]
Out[13]:
In [ ]:
def getAllSimpleCyclesOfSize(self, size):
"""
Return a list of all non-duplicate monocyclic rings with length 'size'.
Naive approach by eliminating polycyclic rings that are returned by
``getAllCyclicsOfSize``.
"""
cycleList = self.getAllCyclesOfSize(size)
i = 0
#import pdb; pdb.set_trace()
while i < len(cycleList):
for vertex in cycleList[i]:
internalConnectivity = sum([1 if vertex2 in cycleList[i] else 0 for vertex2 in vertex.edges.iterkeys()])
if internalConnectivity > 2:
del cycleList[i]
break
else:
i += 1
return cycleList
In [ ]:
from unittest import *
import sys
import rmgpy.molecule.resonanceTest
tests = TestLoader().loadTestsFromModule(rmgpy.molecule.resonanceTest)
TextTestRunner(verbosity=2, stream=sys.stdout).run(tests)
In [ ]:
from rdkit import Chem
def toRDKitMol(mol, removeHs=True, returnMapping=False, sanitize=True):
"""
Convert a molecular structure to a RDKit rdmol object. Uses
`RDKit <http://rdkit.org/>`_ to perform the conversion.
Perceives aromaticity and, unless removeHs==False, removes Hydrogen atoms.
If returnMapping==True then it also returns a dictionary mapping the
atoms to RDKit's atom indices.
"""
# Sort the atoms before converting to ensure output is consistent
# between different runs
mol.sortAtoms()
atoms = mol.vertices
rdAtomIndices = {} # dictionary of RDKit atom indices
rdkitmol = Chem.rdchem.EditableMol(Chem.rdchem.Mol())
for index, atom in enumerate(mol.vertices):
rdAtom = Chem.rdchem.Atom(atom.element.symbol)
rdAtom.SetNumRadicalElectrons(atom.radicalElectrons)
if atom.element.symbol == 'C' and atom.lonePairs == 1 and mol.multiplicity == 1: rdAtom.SetNumRadicalElectrons(2)
rdkitmol.AddAtom(rdAtom)
if removeHs and atom.symbol == 'H':
pass
else:
rdAtomIndices[atom] = index
rdBonds = Chem.rdchem.BondType
orders = {1: rdBonds.SINGLE, 2: rdBonds.DOUBLE, 3: rdBonds.TRIPLE, 1.5: rdBonds.AROMATIC}
# Add the bonds
for atom1 in mol.vertices:
for atom2, bond in atom1.edges.iteritems():
index1 = atoms.index(atom1)
index2 = atoms.index(atom2)
if index1 < index2:
order = orders[bond.order]
rdkitmol.AddBond(index1, index2, order)
# Make editable mol into a mol and rectify the molecule
rdkitmol = rdkitmol.GetMol()
if sanitize:
Chem.SanitizeMol(rdkitmol)
if removeHs:
rdkitmol = Chem.RemoveHs(rdkitmol, sanitize=sanitize)
if returnMapping:
return rdkitmol, rdAtomIndices
return rdkitmol
In [ ]: