

In [ ]:

    
def kekulize(self):
    from rmgpy.molecule import PeriodicSystem
    import numpy as np
    import itertools

    rings = self.getAllCyclesOfSize(6)

    # We want to identify all aromatic rings, especially ones that have been modified
    aromaticRings = []
    invalidRings = []
    for ring in rings:
        endoBonds = set()
        exoBonds = set()
        aromatic = True
        valid = True
        for atom1 in ring:
            for atom2, bond in atom1.bonds.iteritems():
                if atom2 in ring:
                    if bond.order[0] != 'B':
                        aromatic = False
                        break
                    elif bond.order == 'B+' or bond.order == 'B-':
                        valid = False
                    endoBonds.add(bond)
                else:
                    exoBonds.add(bond)
            if not aromatic:
                break
        if aromatic:
            if valid:
                aromaticRings.append((ring, list(endoBonds), list(exoBonds)))
            else:
                invalidRings.append((ring, list(endoBonds), list(exoBonds)))

    bondOrders = {'S': 1, 'D': 2, 'T': 3, 'B': 1.5}
    bondOrdersInv = {1: 'S', 2: 'D', 3: 'T'}
    valences = PeriodicSystem.valences

    for ring, endo, exo in itertools.chain(invalidRings, aromaticRings):
        print 'Solving ring...'
        exoAromatic = sum([1 for bond in exo if bond.order == 'B'])
        print 'There are {} exocyclic aromatic bonds.'.format(exoAromatic)
        endoNonArom = sum([bondOrders[bond.order] for bond in endo if bond.order[0] != 'B'])
        print 'Endocyclic Bond:' + '-'.join([bond.order for bond in endo])

        a = []
        b = []
        for atom in ring:
            a.append([1 if atom in [bond.atom1, bond.atom2] else 0 for bond in endo])
            v = valences[atom.element.symbol]
            for bond in atom.bonds.itervalues():
                if bond in exo:
                    if exoAromatic == 1:
                        v -= np.ceil(bondOrders[bond.order])
                    else:
                        v -= np.floor(bondOrders[bond.order])
            v -= atom.radicalElectrons
            b.append(v)
        constraints = []
        for index, bond in enumerate(endo):
            if bond.order == 'B-':
                # Must become single bond
                new_a = [0] * len(endo)
                new_a[index] = 1
                constraints.append((new_a, 1))
            elif bond.order == 'B+':
                # Must become double bond
                new_a = [0] * len(endo)
                new_a[index] = 1
                constraints.append((new_a, 2))
            elif bond.order != 'B':
                # Bond has already been altered, so keep the new order
                new_a = [0] * len(endo)
                new_a[index] = 1
                constraints.append((new_a, bondOrders[bond.order]))
        if len(constraints) == 0:
            # Move on to next ring and come back later
            print 'No constraints, picking arbitrary constraint.'
            # Pick an arbitrary bond to constrain
            index = 0
            new_a = [0] * len(endo)
            new_a[index] = 1
            constraints.append((new_a, 2))
        for new_a, new_b in constraints:
            a.append(new_a)
            b.append(new_b)

        a = np.array(a)
        b = np.array(b)

        x = np.rint(np.dot(np.linalg.pinv(a), b)).tolist()

        for index, bond in enumerate(endo):
            try:
                bond.order = bondOrdersInv[x[index]]
            except KeyError:
                pdb.set_trace()