from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
import glob
import py3Dmol
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import Draw
from rdkit.Chem import rdMolAlign
from rdkit import rdBase
print(rdBase.rdkitVersion)
import os,time
print( time.asctime())
2016.09.4 Tue Mar 28 15:42:19 2017
# Functions used in this notebook:
def grep_energies_from_sdf_outputs(files):
energies = []
for inp in files:
with open(inp,'r') as f:
lines = f.readlines()
for i, line in enumerate(lines):
if "M END" in line:
energies.append(float(lines[i+1]))
return energies
The crystallographic structure of "M1" (shown below) was used as a starting geometry for generating various conformers.
cm1 = open('/home/gosia/work/work_on_gitlab/icho/calcs/m1/m1_crystal.xyz','r').read()
vcm1 = py3Dmol.view(width=400,height=400)
vcm1.removeAllModels()
vcm1.addModel(cm1,'xyz')
vcm1.setStyle({'stick':{'radius':0.15,'color':'spectrum'}})
vcm1.setBackgroundColor('0xeeeeee')
vcm1.zoomTo()
vcm1.show()
We can also read it in in smiles and sdf formats (both will be used below):
# decide what is the "core" - a part of molecule, which we wish to be most aligned (rmsd-wise) among all the structures
# here: pyridine ring
m1 = Chem.AddHs(Chem.MolFromSmiles('O=C1NCCNC(=O)c2nc(C(=O)NCCNC(=O)c3nc1ccc3)ccc2'))
core_m1 = m1.GetSubstructMatch(Chem.MolFromSmiles('n1ccccc1'))
templ_m1 = Chem.SDMolSupplier('/home/gosia/work/work_on_gitlab/icho/calcs/m1/balloon/m1_crystal.sdf')
m1_crystal = templ_m1[0]
Description: notes
inps_m1_balloon = glob.glob('/home/gosia/work/work_on_gitlab/icho/calcs/m1/balloon/results_crystal/*.sdf')
e_m1_balloon = grep_energies_from_sdf_outputs(inps_m1_balloon)
Conformers:
%%html
<table>
<tr>
<td id="m1_b1" ></td>
<td id="m1_b2" ></td>
<td id="m1_b3" ></td>
</tr>
<tr>
<td> m1_b1, E = </td>
<td> m1_b2, E = </td>
<td> m1_b3, E = </td>
</tr>
<tr>
<td id="m1_b4" ></td>
<td id="m1_b5" ></td>
<td id="m1_b6" ></td>
</tr>
<tr>
<td> m1_b4, E = </td>
<td> m1_b5, E = </td>
<td> m1_b6, E = </td>
</tr>
</table>
m1_b1, E = | m1_b2, E = | m1_b3, E = |
m1_b4, E = | m1_b5, E = | m1_b6, E = |
p1_b_handles=[]
for inp in inps_m1_balloon:
f1_b = open(inp, 'r').read()
p1_b = py3Dmol.view(width=200,height=200)
p1_b.removeAllModels()
p1_b.addModel(f1_b,'sdf')
p1_b.setStyle({'stick':{}})
p1_b.setBackgroundColor('0xeeeeee')
p1_b.zoomTo()
p1_b_handles.append(p1_b)
p1_b_handles[0].insert('m1_b1')
p1_b_handles[1].insert('m1_b2')
p1_b_handles[2].insert('m1_b3')
p1_b_handles[3].insert('m1_b4')
p1_b_handles[4].insert('m1_b5')
p1_b_handles[5].insert('m1_b6')
We can now try to align all these structures (together with the crystal structure):
# write conformers as a list
allmol_m1_balloon = []
suppl_m1_balloon = Chem.SDMolSupplier('/home/gosia/work/work_on_gitlab/icho/calcs/m1/balloon/m1_crystal_out.sdf')
# add crystal structure to the list:
allmol_m1_balloon.append(m1_crystal)
for mol in suppl_m1_balloon:
allmol_m1_balloon.append(mol)
# align:
for mol in allmol_m1_balloon:
AllChem.AlignMolConformers(mol,atomIds=core_m1)
# view:
p = py3Dmol.view(width=400,height=400)
p.removeAllModels()
for mol in allmol_m1_balloon:
mb = Chem.MolToMolBlock(mol)
p.addModel(mb,'sdf')
p.setStyle({'stick':{'radius':'0.15'}})
p.setBackgroundColor('0xeeeeee')
p.zoomTo()
p.show()
# calculate RMSD /check that/:
for mol in allmol_m1_balloon:
# note that the first structure on "allmol_m1_rdkit" list is the crystal structure,
# so the RMSD value calculated for the first structure (with respect to the crystal) will be 0 or almost 0
print("Heavy Atom RMS:",AllChem.GetBestRMS(Chem.RemoveHs(mol),Chem.RemoveHs(templ_m1[0])))
Heavy Atom RMS: 3.604550104229233e-07 Heavy Atom RMS: 0.6621220451255674 Heavy Atom RMS: 0.5174150250164391 Heavy Atom RMS: 1.1744423515041509 Heavy Atom RMS: 0.527526310182849 Heavy Atom RMS: 0.8526332845365369 Heavy Atom RMS: 0.8800960083230364
# create a list of all structures to be aligned
allmol_m1_rdkit = []
suppl_m1_rdkit = Chem.SDMolSupplier('/home/gosia/work/work_on_gitlab/icho/calcs/m1/rdkit/result_smiles_new.sdf')
# add crystal structure to the list:
#allmol_m1_rdkit.append(m1_crystal)
for mol in suppl_m1_rdkit:
allmol_m1_rdkit.append(mol)
View conformers:
%%html
<table>
<tr>
<td id="m1_r1" ></td>
<td id="m1_r2" ></td>
<td id="m1_r3" ></td>
</tr>
<tr>
<td> m1_r1, E = </td>
<td> m1_r2, E = </td>
<td> m1_r3, E = </td>
</tr>
<tr>
<td id="m1_r4" ></td>
<td id="m1_r5" ></td>
<td id="m1_r6" ></td>
</tr>
<tr>
<td> m1_r4, E = </td>
<td> m1_r5, E = </td>
<td> m1_r6, E = </td>
</tr>
<tr>
<td id="m1_r7" ></td>
<td id="m1_r8" ></td>
<td id="m1_r9" ></td>
</tr>
<tr>
<td> m1_r7, E = </td>
<td> m1_r8, E = </td>
<td> m1_r9, E = </td>
</tr>
</table>
m1_r1, E = | m1_r2, E = | m1_r3, E = |
m1_r4, E = | m1_r5, E = | m1_r6, E = |
m1_r7, E = | m1_r8, E = | m1_r9, E = |
inps_m1_rdkit = glob.glob('/home/gosia/work/work_on_gitlab/icho/calcs/m1/rdkit/results_crystal_from_smiles_new/*.sdf')
p1_r_handles=[]
for inp in inps_m1_rdkit:
f1_r = open(inp, 'r').read()
p1_r = py3Dmol.view(width=200,height=200)
p1_r.removeAllModels()
p1_r.addModel(f1_r,'sdf')
p1_r.setStyle({'stick':{}})
p1_r.setBackgroundColor('0xeeeeee')
p1_r.zoomTo()
p1_r_handles.append(p1_r)
p1_r_handles[0].insert('m1_r1')
p1_r_handles[1].insert('m1_r2')
p1_r_handles[2].insert('m1_r3')
p1_r_handles[3].insert('m1_r4')
p1_r_handles[4].insert('m1_r5')
p1_r_handles[5].insert('m1_r6')
p1_r_handles[6].insert('m1_r7')
p1_r_handles[7].insert('m1_r8')
p1_r_handles[8].insert('m1_r9')
We can now align all these structures:
# align:
for mol in allmol_m1_rdkit:
AllChem.AlignMolConformers(mol,atomIds=core_m1)
# view:
p = py3Dmol.view(width=400,height=400)
p.removeAllModels()
for mol in allmol_m1_rdkit:
mb = Chem.MolToMolBlock(mol)
p.addModel(mb,'sdf')
p.setStyle({'stick':{'radius':'0.15'}})
p.setBackgroundColor('0xeeeeee')
p.zoomTo()
p.show()
# calculate RMSD /check that/:
for mol in allmol_m1_rdkit:
# note that the first structure on "allmol_m1_rdkit" list is the crystal structure,
# so the RMSD value calculated for the first structure (with respect to the crystal) will be 0 or almost 0
print("Heavy Atom RMS:",AllChem.GetBestRMS(Chem.RemoveHs(mol),Chem.RemoveHs(templ_m1[0])))
Heavy Atom RMS: 1.461877955872602 Heavy Atom RMS: 1.4214730604414139 Heavy Atom RMS: 1.8999303358326163 Heavy Atom RMS: 1.4813288602775967 Heavy Atom RMS: 1.4248638894808567 Heavy Atom RMS: 1.5474857037372614 Heavy Atom RMS: 1.4688621809782134 Heavy Atom RMS: 1.4939359811879294 Heavy Atom RMS: 1.4176054292304114
Let's align all generated conformers: //todo fixme: align all conformers!//
p = py3Dmol.view(width=400,height=400)
p.removeAllModels()
allmol = allmol_m1_balloon[1:]+allmol_m1_rdkit[1:]
m1b = Chem.MolToMolBlock(m1_crystal)
for i, mol in enumerate(allmol):
AllChem.AlignMolConformers(mol,atomIds=core_m1)
mq = Chem.MolToMolBlock(mol)
p.addModel(mq,'sdf')
p.setStyle({'stick':{'radius':'0.15'}})
p.setBackgroundColor('0xeeeeee')
p.zoomTo()
p.show()
The selection of conformers of "M1" macrocycle:
start: crystal geometry
cm7 = open('/home/gosia/work/work_on_gitlab/icho/calcs/m7/m7_crystal.xyz','r').read()
vcm7 = py3Dmol.view(width=400,height=400)
vcm7.removeAllModels()
vcm7.addModel(cm7,'xyz')
vcm7.setStyle({'stick':{}})
vcm7.setBackgroundColor('0xeeeeee')
vcm7.zoomTo()
vcm7.show()
# find "core" - a part of molecule, which we wish to be most aligned (rmsd-wise) among all the structures
m7 = Chem.AddHs(Chem.MolFromSmiles('N1C(=O)c2nc(C(=O)NCCCNC(=O)c3nc(C(=O)NCCC1)ccc3)ccc2'))
core_m7 = m7.GetSubstructMatch(Chem.MolFromSmiles('C(=O)c1nc(C=O)ccc1'))
templ_m7 = Chem.SDMolSupplier('/home/gosia/work/work_on_gitlab/icho/calcs/m7/m7_crystal.sdf')
m7_crystal = templ_m7[0]
inps_m7_balloon = glob.glob('/home/gosia/work/work_on_gitlab/icho/calcs/m7/balloon/results_crystal/*.sdf')
e_m7_balloon = grep_energies_from_sdf_outputs(inps_m7_balloon)
%%html
<table>
<tr>
<td id="m7_b1" ></td>
<td id="m7_b2" ></td>
<td id="m7_b3" ></td>
<td id="m7_b4" ></td>
</tr>
<tr>
<td> m7_b1, E = </td>
<td> m7_b2, E = </td>
<td> m7_b3, E = </td>
<td> m7_b4, E = </td>
</tr>
<tr>
<td id="m7_b5" ></td>
<td id="m7_b6" ></td>
<td id="m7_b7" ></td>
</tr>
<tr>
<td> m7_b5, E = </td>
<td> m7_b6, E = </td>
<td> m7_b7, E = </td>
</tr>
<tr>
<td id="m7_b8" ></td>
<td id="m7_b9" ></td>
<td id="m7_b10" ></td>
</tr>
<tr>
<td> m7_b8, E = </td>
<td> m7_b9, E = </td>
<td> m7_b10, E = </td>
</tr>
</table>
m7_b1, E = | m7_b2, E = | m7_b3, E = | m7_b4, E = |
m7_b5, E = | m7_b6, E = | m7_b7, E = | |
m7_b8, E = | m7_b9, E = | m7_b10, E = |
p7_b_handles=[]
for inp in inps_m7_balloon:
f7_b = open(inp, 'r').read()
p7_b = py3Dmol.view(width=200,height=200)
p7_b.removeAllModels()
p7_b.addModel(f7_b,'sdf')
p7_b.setStyle({'stick':{}})
p7_b.setBackgroundColor('0xeeeeee')
p7_b.zoomTo()
p7_b_handles.append(p7_b)
p7_b_handles[0].insert('m7_b1')
p7_b_handles[1].insert('m7_b2')
p7_b_handles[2].insert('m7_b3')
p7_b_handles[3].insert('m7_b4')
p7_b_handles[4].insert('m7_b5')
p7_b_handles[5].insert('m7_b6')
p7_b_handles[6].insert('m7_b7')
p7_b_handles[7].insert('m7_b8')
p7_b_handles[8].insert('m7_b9')
p7_b_handles[9].insert('m7_b10')
We can now align all these structures against the crystal structure:
# write conformers as a list
allmol_m7_balloon = []
suppl_m7_balloon = Chem.SDMolSupplier('/home/gosia/work/work_on_gitlab/icho/calcs/m7/balloon/m7_crystal_out.sdf')
# add crystal structure to the list:
allmol_m7_balloon.append(m7_crystal)
for mol in suppl_m7_balloon:
allmol_m7_balloon.append(mol)
# align:
for mol in allmol_m7_balloon:
AllChem.AlignMolConformers(mol,atomIds=core_m7)
p = py3Dmol.view(width=400,height=400)
p.removeAllModels()
for mol in allmol_m7_balloon:
mb = Chem.MolToMolBlock(mol)
p.addModel(mb,'sdf')
p.setStyle({'stick':{'radius':'0.15'}})
p.setBackgroundColor('0xeeeeee')
p.zoomTo()
p.show()
for mol in allmol_m7_balloon:
# note that the first structure on "allmol_m1_rdkit" list is the crystal structure,
# so the RMSD value calculated for the first structure (with respect to the crystal) will be very small
print("Heavy Atom RMS:",AllChem.GetBestRMS(Chem.RemoveHs(mol),Chem.RemoveHs(m7_crystal)))
Heavy Atom RMS: 0.0 Heavy Atom RMS: 0.5937868071467916 Heavy Atom RMS: 0.2092299897502145 Heavy Atom RMS: 2.0075607376161804 Heavy Atom RMS: 1.840950977376659 Heavy Atom RMS: 1.8774918734686699 Heavy Atom RMS: 2.0163759964267918 Heavy Atom RMS: 1.7994332303736011 Heavy Atom RMS: 1.9216906406205834 Heavy Atom RMS: 1.3635632583944308 Heavy Atom RMS: 1.380050979571361
# create a list of all structures to be aligned
inps_m7_rdkit = glob.glob('/home/gosia/work/work_on_gitlab/icho/calcs/m7/rdkit/results_crystal_from_smiles/*.sdf')
allmol_m7_rdkit = []
suppl_m7_rdkit = Chem.SDMolSupplier('/home/gosia/work/work_on_gitlab/icho/calcs/m7/rdkit/result_smiles.sdf')
# add crystal structure to the list:
#allmol_m7_rdkit.append(m7_crystal)
for mol in suppl_m7_rdkit:
allmol_m7_rdkit.append(mol)
e_m7_rdkit = grep_energies_from_sdf_outputs(inps_m7_rdkit)
View conformers:
%%html
<table>
<tr>
<td id="m7_r1" ></td>
<td id="m7_r2" ></td>
<td id="m7_r3" ></td>
</tr>
<tr>
<td> m7_r1, E = </td>
<td> m7_r2, E = </td>
<td> m7_r3, E = </td>
</tr>
<tr>
<td id="m7_r4" ></td>
<td id="m7_r5" ></td>
<td id="m7_r6" ></td>
</tr>
<tr>
<td> m7_r4, E = </td>
<td> m7_r5, E = </td>
<td> m7_r6, E = </td>
</tr>
</table>
m7_r1, E = | m7_r2, E = | m7_r3, E = |
m7_r4, E = | m7_r5, E = | m7_r6, E = |
p7_r_handles=[]
for inp in inps_m1_rdkit:
f7_r = open(inp, 'r').read()
p7_r = py3Dmol.view(width=200,height=200)
p7_r.removeAllModels()
p7_r.addModel(f7_r,'sdf')
p7_r.setStyle({'stick':{}})
p7_r.setBackgroundColor('0xeeeeee')
p7_r.zoomTo()
p7_r_handles.append(p7_r)
p7_r_handles[0].insert('m7_r1')
p7_r_handles[1].insert('m7_r2')
p7_r_handles[2].insert('m7_r3')
p7_r_handles[3].insert('m7_r4')
p7_r_handles[4].insert('m7_r5')
p7_r_handles[5].insert('m7_r6')
We can now align all these structures:
# align:
for mol in allmol_m7_rdkit:
AllChem.AlignMolConformers(mol,atomIds=core_m7)
# view:
p = py3Dmol.view(width=400,height=400)
p.removeAllModels()
for mol in allmol_m7_rdkit:
mb = Chem.MolToMolBlock(mol)
p.addModel(mb,'sdf')
p.setStyle({'stick':{'radius':'0.15'}})
p.setBackgroundColor('0xeeeeee')
p.zoomTo()
p.show()
Let's align all generated conformers: //todo fixme: align all conformers!//
p = py3Dmol.view(width=400,height=400)
p.removeAllModels()
allmol = allmol_m7_balloon[1:]+allmol_m7_rdkit[1:]
m7b = Chem.MolToMolBlock(m7_crystal)
for i, mol in enumerate(allmol):
AllChem.AlignMolConformers(mol,atomIds=core_m7)
mq = Chem.MolToMolBlock(mol)
p.addModel(mq,'sdf')
p.setStyle({'stick':{'radius':'0.15'}})
p.setBackgroundColor('0xeeeeee')
p.zoomTo()
p.show()