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:46:46 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()