Aromaticity Detection

RMG relies solely on RDKit for aromaticity detection. After creating a molecule, we will convert it to an RDKit molecule, then query each of the bonds to determine if RDKit considers them aromatic. If a molecule has a six-membered ring (RMG constraint) and each of the six bonds are considered by RDKit to be aromatic, then RMG considers the ring to be aromatic.

Essentially, we assume that RDKit is always correct, which is not true. A notable example is when rings have double bonded substituents. Because RDKit uses an atom-centered counting method to perform a Huckel rule check, such molecules are considered aromatic. However, we don't think they should be. This notebook tests aromaticity perception by RDKit and Open Babel for such cases.

Interesting link: http://blueobelisk.shapado.com/questions/aromaticity-perception-differences

RDKit



In [ ]:

    
from rdkit import Chem
from rdkit.Chem import Draw

m = Chem.MolFromSmiles('c1ccccc1CCN=[N+]=[N-]')
Draw.MolToImage(m)



In [ ]:

    
print 'Atoms\n====='
for atom in m.GetAtoms():
    print atom.GetIdx(), atom.GetSymbol(), atom.IsInRing(), atom.GetIsAromatic()

print 'Bonds\n====='
for bond in m.GetBonds():
    print bond.GetIdx(), bond.GetBondTypeAsDouble(), bond.GetIsAromatic()

OpenBabel

http://python.zirael.org/e-openbabel8.html



In [ ]:

    
import pybel



In [ ]:

    
m2 = pybel.readstring('smi', 'C1(=C)C=CC=CC1=C')

print 'Atoms\n====='
for atom in m2.atoms:
    print atom.OBAtom.IsAromatic()

print 'Rings\n====='
for ring in m2.sssr:
    print ring.IsAromatic()



In [ ]:

    
import rmgpy.molecule.generator as generator
from rmgpy.molecule import Molecule
from IPython.display import display

mol = Molecule().fromAdjacencyList("""
1  C u0 p0 c0 {2,D} {6,S} {7,S}
2  C u0 p0 c0 {1,D} {3,S} {8,S}
3  C u0 p0 c0 {2,S} {4,D} {9,S}
4  C u0 p0 c0 {3,D} {5,S} {10,S}
5  C u0 p0 c0 {4,S} {6,D} {11,S}
6  C u0 p0 c0 {1,S} {5,D} {12,S}
7  O u1 p2 c0 {1,S}
8  H u0 p0 c0 {2,S}
9  H u0 p0 c0 {3,S}
10 H u0 p0 c0 {4,S}
11 H u0 p0 c0 {5,S}
12 H u0 p0 c0 {6,S}
""")
display(mol)
mol.getBond(mol.atoms[0], mol.atoms[6]).order = 2
mol.atoms[6].radicalElectrons = 0
print mol.toAdjacencyList()



In [ ]:

    
obmol, obAtomIndices = generator.toOBMol(mol, returnMapping=True)



In [ ]:

    
obAtomIndices



In [ ]:

    
dir(obmol)



In [ ]: