In [ ]:

    

from rmgpy.molecule import Molecule
from rmgpy.molecule.resonance import *
from rmgpy.molecule import PeriodicSystem

import numpy as np
import itertools
from IPython.display import display
import pdb


    

In [ ]:

    

mol = Molecule(SMILES='[C]1=CC2=CC=C3C=CC4=CC=C5C=CC6=CC=C1C1=C2C3=C4C5=C61')
out = generateAromaticResonanceStructures(mol)
display(out[0])
out[0].kekulize()
display(out[0])


    

In [ ]:

    

def kekulize(self):

    rings = self.getAllCyclesOfSize(6)
    
    # Identify aromatic rings
    aromaticRings = []
    for ring in rings:
        endoBonds = set()
        exoBonds = set()
        aromatic = True
        for atom1 in ring:
            for atom2, bond in atom1.bonds.iteritems():
                if atom2 in ring:
                    if bond.order[0] != 'B':
                        aromatic = False
                        break
                    endoBonds.add(bond)
                else:
                    exoBonds.add(bond)
            if not aromatic:
                break
        if aromatic:
            aromaticRings.append([ring, list(endoBonds), list(exoBonds), -1, -1])

    updateKekulizationDOF(aromaticRings)
    resolvedRings = []
    itercount = 0
    maxiter = 2 * len(aromaticRings)
    while aromaticRings and itercount < maxiter:
        ringData = aromaticRings.pop()
        successful = kekulizeRing(*ringData)
        if successful:
            resolvedRings.append(ringData)
            updateKekulizationDOF(aromaticRings)
        else:
            print 'Unsuccessful kekulization!'
            aromaticRings.append(ringData)
            updateKekulizationDOF(aromaticRings)
        itercount += 1
    
    try:
        self.updateAtomTypes()
    except:
        pdb.set_trace()

def updateKekulizationDOF(aromaticRings):
    for ringData in aromaticRings:
        endoDOF = 6
        for bond in ringData[1]:
            if bond.order != 'B':
                # Remove one dof for each bond that is not aromatic
                endoDOF -= 1
        exoDOF = 6
        for bond in ringData[2]:
            if bond.order != 'B':
                # Remove one dof for each bond that is not aromatic
                exoDOF -= 1
        for atom in ringData[0]:
            exoDOF -= atom.radicalElectrons
            exoDOF -= atom.lonePairs
        ringData[3] = endoDOF
        ringData[4] = exoDOF
    
    aromaticRings.sort(key=lambda x: x[3] + x[4], reverse=True)
    
    
def kekulizeRing(ring, endoBonds, exoBonds, endoDOF, exoDOF):
    print 'Solving ring'
    
    bondOrders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5}
    bondOrdersInv = {1: 'S', 2: 'D', 3: 'T'}
    valences = PeriodicSystem.valences

    
    # Make a copy of the list
    unresolved = list(endoBonds)
    resolved = []
    # Check if any bonds are predetermined
    for bond in endoBonds:
        if bond.order[0] != 'B':
            # Bond has already been assigned, so mark as resolved
            resolved.append(bond)
            unresolved.remove(bond)
        elif bond.order == 'B+':
            # Bond was incremented, so it must be a double bond
            bond.order = 'D'
            resolved.append(bond)
            unresolved.remove(bond)
        elif bond.order == 'B-':
            # Bond was decremented, so it must be a single bond
            bond.order = 'S'
            resolved.append(bond)
            unresolved.remove(bond)
    
    # Check status
    if len(resolved) == len(endoBonds):
        return True
    
    itercount = 0
    maxiter = 2 * len(endoBonds)
    while unresolved and itercount < maxiter:
        bond0 = unresolved.pop()
        # Check available valence for both atoms
        doublePossible = True
        doubleRequired = False
        localEndoDOF = 2
        localExoDOF = 2
        for atom in [bond0.atom1, bond0.atom2]:
            occupied = 0
            uncertain = 0
            for bond in atom.bonds.itervalues():
                if bond.order[0] != 'B':
                    occupied += bondOrders[bond.order]
                    if bond is not bond0:
                        if bond in endoBonds:
                            localEndoDOF -= 1
                        else:
                            localExoDOF -= 1
                elif bond.order == 'B':
                    # The atom has a benzene bond, so at least one electron is occupied, but there is a second uncertain electron
                    occupied += 1
                    uncertain += 1
                else:
                    ValueError('Unexpected bond order {0}.'.format(bond.order))
            occupied += atom.radicalElectrons
            localExoDOF -= atom.radicalElectrons
            occupied += 2 * atom.lonePairs
            localExoDOF -= atom.lonePairs
            available = valences[atom.element.symbol] - occupied
            if available < 0:
                ValueError('Atom {0} cannot have negative available valence.'.format(atom))
            elif available == 0:
                doublePossible = False
                break
            elif available == 1 and uncertain == 1:
                doubleRequired = True
                break
        print doublePossible, doubleRequired
        if doublePossible and doubleRequired:
            bond0.order = 'D'
            resolved.append(bond0)
            endoDOF -= 1
        elif doublePossible and not doubleRequired:
            if (endoDOF == 6 and exoDOF == 0) or (endoDOF == 6 and localExoDOF == 0) or (localEndoDOF == 1 and localExoDOF == 2) or endoDOF == 1:
                # Go ahead an assume this bond is double
                bond0.order = 'D'
                resolved.append(bond0)
                endoDOF -= 1
            else:
                # Come back to this bond later
                print 'Postponing resolution of this bond with endoDOF={0} and exoDOF={1}...'.format(endoDOF, exoDOF)
                unresolved.insert(0, bond0)
        else:
            bond0.order = 'S'
            resolved.append(bond0)
            endoDOF -= 1
            
        itercount += 1
        
    if unresolved:
        print 'Unable to solve ring...'
        return False
        
    print '-'.join([bond.order for bond in endoBonds])
            
    return True


    

In [ ]: