Aromatic Oxygen (and Nitrogen, for that matter)

rmgpy/molecule/atomtype.py

@@ -170,7 +170,7 @@ def isSpecificCaseOf(self, other):
     'H','He',
     'C','Cs','Cd','Cdd','Ct','CO','Cb','Cbf','CS',
     'N','N1d','N3s','N3d','N3t','N3b','N5s','N5d','N5dd','N5t','N5b',
-    'O','Os','Od','Oa','Ot',
+    'O','Os','Od','Oa','Ot','Ob',
     'Ne',
     'Si','Sis','Sid','Sidd','Sit','SiO','Sib','Sibf',
     'S','Ss','Sd','Sa',
@@ -181,7 +181,7 @@ def isSpecificCaseOf(self, other):
     'Val4','Val5','Val6','Val7',
     'C','Cs','Cd','Cdd','Ct','CO','Cb','Cbf','CS',
     'N','N1d','N3s','N3d','N3t','N3b','N5s','N5d','N5dd','N5t','N5b',
-    'O','Os','Od','Oa','Ot',
+    'O','Os','Od','Oa','Ot','Ob',
     'Ne',
     'Si','Sis','Sid','Sidd','Sit','SiO','Sib','Sibf',
     'S','Ss','Sd','Sa',
@@ -195,7 +195,7 @@ def isSpecificCaseOf(self, other):
     'N','N1d','N3s','N3d','N3t','N3b','N5s','N5d','N5dd','N5t','N5b']
 )
 atomTypes['Val6'] = AtomType(label='Val6', generic=['R','R!H'], specific=[
-    'O','Os','Od','Oa','Ot',
+    'O','Os','Od','Oa','Ot','Ob',
     'S','Ss','Sd','Sa']
 )
 atomTypes['Val7'] = AtomType(label='Val7', generic=['R','R!H'], specific=[
@@ -228,11 +228,12 @@ def isSpecificCaseOf(self, other):
 atomTypes['N5t' ] = AtomType('N5t',  generic=['R','R!H','N','Val5'],  specific=[])
 atomTypes['N5b' ] = AtomType('N5b',  generic=['R','R!H','N','Val5'],  specific=[])
 
-atomTypes['O'   ] = AtomType('O',    generic=['R','R!H','Val6'],      specific=['Os','Od','Oa','Ot'])
+atomTypes['O'   ] = AtomType('O',    generic=['R','R!H','Val6'],      specific=['Os','Od','Oa','Ot','Ob'])
 atomTypes['Os'  ] = AtomType('Os',   generic=['R','R!H','O','Val6'],  specific=[])
 atomTypes['Od'  ] = AtomType('Od',   generic=['R','R!H','O','Val6'],  specific=[])
 atomTypes['Oa'  ] = AtomType('Oa',   generic=['R','R!H','O','Val6'],  specific=[])
 atomTypes['Ot'  ] = AtomType('Ot',   generic=['R','R!H','O','Val6'],  specific=[])
+atomTypes['Ob'  ] = AtomType('Ob',   generic=['R','R!H','O','Val6'],  specific=[])
 
 atomTypes['Ne'  ] = AtomType('Ne',   generic=['R','R!H'],      specific=[])
 
@@ -293,6 +294,8 @@ def isSpecificCaseOf(self, other):
 atomTypes['Od'  ].setActions(incrementBond=[],               decrementBond=['Os'],           formBond=[],            breakBond=[],            incrementRadical=[],       decrementRadical=[],       incrementLonePair=['Od'],  decrementLonePair=['Od'])
 atomTypes['Oa'  ].setActions(incrementBond=[],               decrementBond=[],               formBond=[],            breakBond=[],            incrementRadical=[],       decrementRadical=[],       incrementLonePair=[],      decrementLonePair=[])
 atomTypes['Ot'  ].setActions(incrementBond=[],               decrementBond=['Od'],           formBond=[],            breakBond=[],            incrementRadical=[],       decrementRadical=[],       incrementLonePair=['Ot'],  decrementLonePair=['Ot'])
+atomTypes['Ob'  ].setActions(incrementBond=[],               decrementBond=[],               formBond=['Ob'],        breakBond=['Ob'],        incrementRadical=['Ob'],   decrementRadical=['Ob'],   incrementLonePair=['Ob'],  decrementLonePair=['Ob']) # not sure if we should have incrementLonePair and dercementLonePair options here?
+
 
 atomTypes['Ne'  ].setActions(incrementBond=[],               decrementBond=[],               formBond=[],            breakBond=[],            incrementRadical=['Ne'],   decrementRadical=['Ne'],   incrementLonePair=[],    decrementLonePair=[])
 
@@ -379,6 +382,7 @@ def getAtomType(atom, bonds):
     elif atom.symbol == 'O':
         if   double + doubleO == 0 and triple == 0 and benzene == 0: atomType = 'Os'
         elif double + doubleO == 1 and triple == 0 and benzene == 0: atomType = 'Od'
+        elif double + doubleO == 0 and triple == 0 and benzene == 2: atomType = 'Ob'
         elif len(bonds) == 0:                                        atomType = 'Oa'
         elif double + doubleO == 0 and triple == 1 and benzene == 0: atomType = 'Ot'
     elif atom.symbol == 'Ne':

rmgpy/molecule/atomtypeTest.py

@@ -257,6 +257,21 @@ def testOtherTypes(self):
         self.assertEqual(self.atomType(self.mol6, 0), 'Ar')
         self.assertEqual(self.atomType(self.mol7, 0), 'He')
         self.assertEqual(self.atomType(self.mol8, 0), 'Ne')
+
+    def testFuranOxygen(self):
+        """
+        Test the O in aromatic resonance form of furan is Ob
+        """
+        furan = Molecule().fromSMILES('C1=COC=C1')
+        aromatics = furan.getAromaticResonanceIsomers()
+        self.assertEqual(len(aromatics), 1)
+        aromatic = aromatics[0]
+        aromatic.updateAtomTypes()
+        #self.assertTrue(aromatic.isAromatic())
+        for atom in aromatic.atoms:
+            if atom.isOxygen():
+                type = getAtomType(atom, aromatic.getBonds(atom))
+                self.assertEqual(type.label, 'Ob')
 
 ################################################################################
 

rmgpy/molecule/molecule.py

@@ -1577,26 +1577,26 @@ def isLinear(self):
     def isAromatic(self):
         """ 
         Returns ``True`` if the molecule is aromatic, or ``False`` if not.  
-        Iterates over the SSSR's and searches for rings that consist solely of Cb 
-        atoms.  Assumes that aromatic rings always consist of 6 atoms. 
-        In cases of naphthalene, where a 6 + 4 aromatic system exists,
-        there will be at least one 6 membered aromatic ring so this algorithm
-        will not fail for fused aromatic rings.
+        Iterates over the SSSR's and searches for rings that consist solely of Xb 
+        atoms, where X could be anything (i.e. Cb, Ob, N3b, N5b).
+        If at least one ring of 'b' atoms is found, then it's Aromatic.
+        Be sure to call updateAtomTypes() before using this.
         """
-        cython.declare(SSSR=list, vertices=list, polycyclicVertices=list)
+        cython.declare(SSSR=list, vertices=list, polycyclicVertices=list, label=str)
         SSSR = self.getSmallestSetOfSmallestRings()
         if SSSR:
             for cycle in SSSR:
-                if len(cycle) == 6:
                     for atom in cycle:
-                        #print atom.atomType.label
-                        if atom.atomType.label == 'Cb' or atom.atomType.label == 'Cbf':
-                            continue                        
-                        # Go onto next cycle if a non Cb atomtype was discovered in this cycle
-                        break 
+                        label = atom.atomType.label
+                        #print 'in isAromatic cycle: '+label
+                        if label[-1] == 'b' or label[-2:] == 'bf':
+                            continue
+                        # Go on to next cycle if a non-b atomtype was discovered in this cycle
+                        break
                     else:
-                        # Molecule is aromatic when all 6 atoms are type 'Cb'
-                        return True    
+                        # All n atoms in this ring are some type of 'b'
+                        return True 
+            # exhausted all rings without finding an aromatic one
         return False
 
     def countInternalRotors(self):
@@ -1896,24 +1896,18 @@ def getAromaticResonanceIsomers(self):
             aromatic = False
             rings = molecule.getSmallestSetOfSmallestRings()            
             for ring0 in rings:
-                # In RMG, only 6-member rings can be considered aromatic, so ignore all other rings                
                 aromaticBonds = []
-                if len(ring0) == 6:
-                    # Figure out which atoms and bonds are aromatic and reassign appropriately:
-                    for i, atom1 in enumerate(ring0):
-                        if not atom1.isCarbon():
-                            # all atoms in the ring must be carbon in RMG for our definition of aromatic
-                            break
-                        for atom2 in ring0[i+1:]:
-                            if molecule.hasBond(atom1, atom2):
-                                if str(rdkitmol.GetBondBetweenAtoms(rdAtomIndices[atom1],rdAtomIndices[atom2]).GetBondType()) == 'AROMATIC':
-                                    aromaticBonds.append(molecule.getBond(atom1, atom2))
-                if len(aromaticBonds) == 6:
+                # Figure out which atoms and bonds are aromatic and reassign appropriately:
+                for i, atom1 in enumerate(ring0):
+                    for atom2 in ring0[i+1:]:
+                        if molecule.hasBond(atom1, atom2):
+                            if str(rdkitmol.GetBondBetweenAtoms(rdAtomIndices[atom1],rdAtomIndices[atom2]).GetBondType()) == 'AROMATIC':
+                                aromaticBonds.append(molecule.getBond(atom1, atom2))
+                if len(aromaticBonds) == len(ring0):
                     aromatic = True
-                    # Only change bonds if there are all 6 are aromatic.  Otherwise don't do anything
+                    # Only change bonds if there are all aromatic.  Otherwise don't do anything
                     for bond in aromaticBonds:
                         bond.order = 'B'
-
             if aromatic:              
                 isomers.append(molecule)
 

rmgpy/molecule/moleculeTest.py

@@ -986,6 +986,14 @@ def testAromaticBenzene(self):
         m = Molecule().fromSMILES('C1=CC=CC=C1')
         isomers = m.generateResonanceIsomers()
         self.assertTrue(any(isomer.isAromatic() for isomer in isomers))
+
+    def testAromaticFuran(self):
+        """
+        Test the Molecule.isAromatic() method for Furan.
+        """
+        m = Molecule().fromSMILES('C1=COC=C1')
+        isomers = m.generateResonanceIsomers()
+        self.assertTrue(any(isomer.isAromatic() for isomer in isomers))
 
     def testAromaticNaphthalene(self):
         """