In [1]:
import rmgpy
from rmgpy.molecule.resonance import *
from rmgpy.molecule.molecule import Molecule
from rmgpy.species import Species
from rmgpy.data.rmg import RMGDatabase
import mipshell
In [2]:
molecule = Molecule().fromSMILES('C1=CC=C2C(C=CC3C2=CC=C2C=CC=CC=32)=C1')
In [3]:
molecule
Out[3]:
In [4]:
molecule.isAromatic()
Out[4]:
In [5]:
SSSR = molecule.getAromaticSSSR()
SSSR
Out[5]:
In [8]:
atoms = set()
for ring in SSSR:
atoms.update(ring)
atoms = list(atoms)
# Get list of bonds involving the ring atoms, ignoring bonds to hydrogen
bonds = set()
for atom in atoms:
bonds.update([atom.bonds[key] for key in atom.bonds.keys() if key.isNonHydrogen()])
bonds = list(bonds)
a = len(SSSR)
b = len(bonds)
# Connectivity matrix which indicates which rings and bonds each atom is in
# Part of equality constraint Ax=b
c = []
for atom in atoms:
inRing = [1 if atom in ring else 0 for ring in SSSR]
inBond = [1 if atom in [bond.atom1, bond.atom2] else 0 for bond in bonds]
c.append(inRing + inBond)
z = [1] * len(SSSR) + [0] * len(bonds)
In [31]:
class clar(mipshell.Problem):
def model(self, a, b, c, constraints=None):
"""
a is number of rings
b is number of bonds
c is constraint matrix
"""
x = mipshell.VarVector([a + b], 'x', type=mipshell.BIN)
mipshell.maximize(mipshell.sum_(x[i] for i in range(a)))
for row in c:
mipshell.sum_(row[i]*x[i] for i in range(a + b)) == 1
if constraints is not None:
for cons in constraints:
mipshell.sum_(cons[0]*x[i] for i in range(a + b)) == cons[1]
def getSolution(self):
x = [var.val for var in self.vars]
return x
In [32]:
problem = clar('test')
problem.model(a, b, c)
problem.optimize()
In [33]:
problem.getSolution()
Out[33]:
In [41]:
lp = glpk.LPX()
lp.obj.maximize = True
lp.rows.add(len(atoms))
for r in lp.rows:
r.bounds = (1, 1)
lp.cols.add(len(SSSR) + len(bonds))
for c in lp.cols:
c.bounds = (0, 1)
c.kind = bool
lp.obj[:] = z
lp.matrix = a
lp.intopt()
print lp.status
In [29]:
lp.nbin
Out[29]:
In [33]:
lp.obj.value
Out[33]:
In [38]:
[c.value for c in lp.cols]
Out[38]:
In [40]:
lp.kind
Out[40]:
In [ ]: