In [ ]:
import rmgpy
from rmgpy.data.rmg import RMGDatabase
from rmgpy.reaction import Reaction
from rmgpy.molecule import Molecule, Bond, PeriodicSystem
from rmgpy.molecule.resonance import *
import numpy as np
import time
In [ ]:
databasePath = rmgpy.settings['database.directory']
database = RMGDatabase()
database.load(
path = databasePath,
thermoLibraries = [],
reactionLibraries = [],
seedMechanisms = [],
kineticsFamilies = ['R_Addition_MultipleBond', 'Intra_R_Add_Exocyclic', 'Diels_alder_addition'],
)
In [ ]:
reactants = [
Molecule().fromAdjacencyList("""1 C u0 p0 c0 {2,B} {6,B} {7,S}
2 C u0 p0 c0 {1,B} {3,B} {8,S}
3 C u0 p0 c0 {2,B} {4,B} {9,S}
4 C u0 p0 c0 {3,B} {5,B} {10,S}
5 C u0 p0 c0 {4,B} {6,B} {11,S}
6 C u0 p0 c0 {1,B} {5,B} {12,S}
7 H u0 p0 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}
"""),
Molecule(SMILES='[H]')
]
rxns = database.kinetics.families['R_Addition_MultipleBond'].generateReactions(reactants)
rxns[0]
In [ ]:
reactants = [
Molecule().fromAdjacencyList("""1 C u0 p0 c0 {2,B} {3,B} {6,S}
2 C u0 p0 c0 {1,B} {4,B} {5,S}
3 C u0 p0 c0 {1,B} {7,B} {11,S}
4 C u0 p0 c0 {2,B} {8,B} {14,S}
5 C u0 p0 c0 {2,S} {9,D} {15,S}
6 C u0 p0 c0 {1,S} {10,D} {18,S}
7 C u0 p0 c0 {3,B} {8,B} {12,S}
8 C u0 p0 c0 {4,B} {7,B} {13,S}
9 C u0 p0 c0 {5,D} {10,S} {16,S}
10 C u0 p0 c0 {6,D} {9,S} {17,S}
11 H u0 p0 c0 {3,S}
12 H u0 p0 c0 {7,S}
13 H u0 p0 c0 {8,S}
14 H u0 p0 c0 {4,S}
15 H u0 p0 c0 {5,S}
16 H u0 p0 c0 {9,S}
17 H u0 p0 c0 {10,S}
18 H u0 p0 c0 {6,S}
"""),
Molecule(SMILES='[H]')
]
rxns = database.kinetics.families['R_Addition_MultipleBond'].generateReactions(reactants)
rxns[0]
In [ ]:
rxns
In [ ]:
reactants = [
Molecule().fromAdjacencyList("""multiplicity 2
1 C u0 p0 c0 {2,S} {3,S} {12,S} {13,S}
2 C u0 p0 c0 {1,S} {9,S} {10,S} {11,S}
3 C u0 p0 c0 {1,S} {4,B} {5,B}
4 C u0 p0 c0 {3,B} {6,B} {14,S}
5 C u0 p0 c0 {3,B} {8,B} {18,S}
6 C u0 p0 c0 {4,B} {7,B} {15,S}
7 C u0 p0 c0 {6,B} {8,B} {16,S}
8 C u0 p0 c0 {5,B} {7,B} {17,S}
9 C u1 p0 c0 {2,S} {19,S} {20,S}
10 H u0 p0 c0 {2,S}
11 H u0 p0 c0 {2,S}
12 H u0 p0 c0 {1,S}
13 H u0 p0 c0 {1,S}
14 H u0 p0 c0 {4,S}
15 H u0 p0 c0 {6,S}
16 H u0 p0 c0 {7,S}
17 H u0 p0 c0 {8,S}
18 H u0 p0 c0 {5,S}
19 H u0 p0 c0 {9,S}
20 H u0 p0 c0 {9,S}
""")
]
rxns = database.kinetics.families['Intra_R_Add_Exocyclic'].generateReactions(reactants)
rxns[0]
In [ ]:
rxns
In [ ]:
reactants = [
Molecule().fromAdjacencyList("""1 C u0 p0 c0 {4,B} {5,B} {7,S}
2 C u0 p0 c0 {3,B} {5,B} {8,S}
3 C u0 p0 c0 {2,B} {6,B} {9,S}
4 C u0 p0 c0 {1,B} {6,B} {10,S}
5 C u0 p0 c0 {1,B} {2,B} {18,S}
6 C u0 p0 c0 {3,B} {4,B} {23,S}
7 C u0 p0 c0 {1,S} {12,D} {17,S}
8 C u0 p0 c0 {2,S} {13,D} {19,S}
9 C u0 p0 c0 {3,S} {14,D} {22,S}
10 C u0 p0 c0 {4,S} {11,D} {24,S}
11 C u0 p0 c0 {10,D} {12,S} {15,S}
12 C u0 p0 c0 {7,D} {11,S} {16,S}
13 C u0 p0 c0 {8,D} {14,S} {20,S}
14 C u0 p0 c0 {9,D} {13,S} {21,S}
15 H u0 p0 c0 {11,S}
16 H u0 p0 c0 {12,S}
17 H u0 p0 c0 {7,S}
18 H u0 p0 c0 {5,S}
19 H u0 p0 c0 {8,S}
20 H u0 p0 c0 {13,S}
21 H u0 p0 c0 {14,S}
22 H u0 p0 c0 {9,S}
23 H u0 p0 c0 {6,S}
24 H u0 p0 c0 {10,S}
"""),
Molecule(SMILES='C=C')
]
rxns = database.kinetics.families['Diels_alder_addition'].generateReactions(reactants)
rxns
In [ ]:
rxns[1]
In [ ]:
mol = Molecule(SMILES='c12ccccc1cccc2')
out = generateResonanceStructures(mol)
#print out[0].toAdjacencyList()
out[1]
In [ ]:
out = Molecule(SMILES='C1=CC=C2C=C3C=CC=CC3=CC2=C1').generateResonanceIsomers()
print out[1].toAdjacencyList()
In [ ]:
out[1]
In [ ]:
mol = Molecule().fromAdjacencyList("""1 C u0 p0 c0 {2,B} {6,B} {7,S}
2 C u0 p0 c0 {1,B} {3,B} {8,S}
3 C u0 p0 c0 {2,B} {4,B} {9,S}
4 C u0 p0 c0 {3,B} {5,B} {10,S}
5 C u0 p0 c0 {4,B} {6,B} {11,S}
6 C u0 p0 c0 {1,B} {5,B} {12,S}
7 H u0 p0 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}
13 C u0 p0 c0 {14,D} {15,S} {16,S}
14 C u0 p0 c0 {13,D} {17,S} {18,S}
15 H u0 p0 c0 {13,S}
16 H u0 p0 c0 {13,S}
17 H u0 p0 c0 {14,S}
18 H u0 p0 c0 {14,S}
""")
mol.addBond(Bond(mol.atoms[0], mol.atoms[12], order='S'))
mol.addBond(Bond(mol.atoms[3], mol.atoms[13], order='S'))
bond = mol.getBond(mol.atoms[12], mol.atoms[13])
bond.order = 'S'
bond = mol.getBond(mol.atoms[0], mol.atoms[1])
bond.order = 'B-'
bond = mol.getBond(mol.atoms[1], mol.atoms[2])
bond.order = 'B+'
bond = mol.getBond(mol.atoms[2], mol.atoms[3])
bond.order = 'B-'
print mol.toAdjacencyList()
t0 = time.time()
mol.redistributeAromaticElectrons()
t1 = time.time()
print t1 - t0
mol.update()
mol
In [ ]:
mol = Molecule().fromAdjacencyList("""1 C u0 p0 c0 {2,S} {3,B} {4,B}
2 C u0 p0 c0 {1,S} {5,B} {6,B}
3 C u0 p0 c0 {1,B} {8,B} {15,S}
4 C u0 p0 c0 {1,B} {9,B} {16,S}
5 C u0 p0 c0 {2,B} {10,B} {18,S}
6 C u0 p0 c0 {2,B} {12,B} {22,S}
7 C u0 p0 c0 {8,B} {9,B} {13,S}
8 C u0 p0 c0 {3,B} {7,B} {14,S}
9 C u0 p0 c0 {4,B} {7,B} {17,S}
10 C u0 p0 c0 {5,B} {11,B} {19,S}
11 C u0 p0 c0 {10,B} {12,B} {20,S}
12 C u0 p0 c0 {6,B} {11,B} {21,S}
13 H u0 p0 c0 {7,S}
14 H u0 p0 c0 {8,S}
15 H u0 p0 c0 {3,S}
16 H u0 p0 c0 {4,S}
17 H u0 p0 c0 {9,S}
18 H u0 p0 c0 {5,S}
19 H u0 p0 c0 {10,S}
20 H u0 p0 c0 {11,S}
21 H u0 p0 c0 {12,S}
22 H u0 p0 c0 {6,S}
23 C u0 p0 c0 {24,D} {25,S} {26,S}
24 C u0 p0 c0 {23,D} {27,S} {28,S}
25 H u0 p0 c0 {23,S}
26 H u0 p0 c0 {23,S}
27 H u0 p0 c0 {24,S}
28 H u0 p0 c0 {24,S}
""")
mol.addBond(Bond(mol.atoms[2], mol.atoms[22], order='S'))
mol.addBond(Bond(mol.atoms[4], mol.atoms[23], order='S'))
bond = mol.getBond(mol.atoms[22], mol.atoms[23])
bond.decrementOrder()
bond = mol.getBond(mol.atoms[0], mol.atoms[2])
bond.decrementOrder()
bond = mol.getBond(mol.atoms[0], mol.atoms[1])
bond.incrementOrder()
bond = mol.getBond(mol.atoms[1], mol.atoms[4])
bond.decrementOrder()
print mol.toAdjacencyList()
t0 = time.time()
mol.redistributeAromaticElectrons()
t1 = time.time()
print t1 - t0
mol.update()
mol
In [ ]:
mol = Molecule().fromAdjacencyList("""1 C u0 p0 c0 {2,B} {6,B} {7,S}
2 C u0 p0 c0 {1,B} {3,B} {8,S}
3 C u0 p0 c0 {2,B} {4,B} {9,S}
4 C u0 p0 c0 {3,B} {5,B} {10,S}
5 C u0 p0 c0 {4,B} {6,B} {11,S}
6 C u0 p0 c0 {1,B} {5,B} {12,S}
7 H u0 p0 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}
13 H u1 p0 c0
""")
mol.atoms[0].applyAction(['CHANGE_BOND', '*1', '-1', '*2'])
mol.atoms[1].applyAction(['CHANGE_BOND', '*1', '-1', '*2'])
bond = mol.getBond(mol.atoms[0], mol.atoms[1])
bond.decrementOrder()
mol.addBond(Bond(mol.atoms[0], mol.atoms[12], order='S'))
mol.atoms[1].radicalElectrons += 1
mol.atoms[12].radicalElectrons -= 1
print mol.toAdjacencyList()
t0 = time.time()
mol.redistributeAromaticElectrons()
t1 = time.time()
print t1 - t0
mol.update()
mol
In [ ]:
def redistributeAromaticElectrons(mol):
"""
Redistributes pi-electrons around aromatic rings after incrementing or decrementing a benzene bond.
"""
from rmgpy.molecule import PeriodicSystem
import numpy as np
rings = mol.getAllCyclesOfSize(6)
invalidRings = []
# We want to identify previously aromatic rings that have been modified
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 and not valid:
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 invalidRings:
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:
v -= bondOrders[bond.order]
v -= atom.radicalElectrons
b.append(v)
for index, bond in enumerate(endo):
if bond.order == 'B-':
# Must become single bond
new_a = [0] * len(endo)
new_a[index] = 1
a.append(new_a)
b.append(1)
elif bond.order == 'B+':
# Must become double bond
new_a = [0] * len(endo)
new_a[index] = 1
a.append(new_a)
b.append(2)
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):
bond.order = bondOrdersInv[x[index]]
In [ ]:
def redistributeAromaticElectrons(mol):
sssr = mol.getAromaticRings()
if sssr == []:
return
rings = [ring for ring in sssr if len(ring) == 6]
invalidRings = []
for ring in rings:
endoBonds = set()
exoBonds = set()
valid = True
for atom1 in ring:
for atom2, bond in atom1.bonds.iteritems():
if atom2 in ring:
endoBonds.add(bond)
if bond.order == 'B+' or bond.order == 'B-':
valid = False
else:
exoBonds.add(bond)
if not valid:
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 invalidRings:
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:
v -= bondOrders[bond.order]
v -= atom.radicalElectrons
b.append(v)
for index, bond in enumerate(endo):
if bond.order == 'B-':
# Must become single bond
new_a = [0] * len(endo)
new_a[index] = 1
a.append(new_a)
b.append(1)
elif bond.order == 'B+':
# Must become double bond
new_a = [0] * len(endo)
new_a[index] = 1
a.append(new_a)
b.append(2)
a = np.array(a)
b = np.array(b)
x = np.rint(np.dot(np.linalg.pinv(a), b)).tolist()
import pdb; pdb.set_trace()
for index, bond in enumerate(endo):
bond.order = bondOrdersInv[x[index]]
In [ ]: