import glob
import py3Dmol
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import Draw
from rdkit import rdBase
from rdkit.Chem import rdMolAlign
print(rdBase.rdkitVersion)
import os,time
print( time.asctime())
2016.09.4 Tue Mar 28 08:30:39 2017
First let's reproduce the visualization from this notebook
m = Chem.AddHs(Chem.MolFromSmiles('COc1ccc2[C@H](O)[C@@H](COc2c1)N3CCC(O)(CC3)c4ccc(F)cc4'))
print(m)
AllChem.EmbedMultipleConfs(m,useExpTorsionAnglePrefs=True,useBasicKnowledge=True)
mb = Chem.MolToMolBlock(m)
core = m.GetSubstructMatch(Chem.MolFromSmiles('C1C(O)c2ccccc2OC1'))
print(core)
AllChem.AlignMolConformers(m,atomIds=core)
p = py3Dmol.view(width=400,height=400)
for conf in m.GetConformers():
mb = Chem.MolToMolBlock(m,confId=conf.GetId())
p.addModel(mb,'sdf')
p.setStyle({'stick':{'radius':0.1}})
p.setBackgroundColor('0xeeeeee')
p.zoomTo()
p.show()
<rdkit.Chem.rdchem.Mol object at 0x7fb5b009c440> (8, 6, 7, 5, 4, 3, 2, 12, 11, 10, 9)
Now: examples which do not work:
# create a list of all structures to be aligned
allmol = []
# conformers:
suppl = Chem.SDMolSupplier('/home/gosia/work/work_on_gitlab/icho/calcs/m1/rdkit/result_smiles_new.sdf')
# crystal:
reference = Chem.SDMolSupplier('/home/gosia/work/work_on_gitlab/icho/calcs/m1/balloon/m1_crystal.sdf')[0]
# add reference structure to the list (crystal geom):
allmol.append(reference)
for mol in suppl:
allmol.append(mol)
# find "core":
m = Chem.AddHs(Chem.MolFromSmiles('O=C1NCCNC(=O)c2nc(C(=O)NCCNC(=O)c3nc1ccc3)ccc2'))
core = m.GetSubstructMatch(Chem.MolFromSmiles('n1ccccc1'))
# align:
for mol in allmol:
AllChem.AlignMolConformers(mol,atomIds=core)
# view:
p = py3Dmol.view(width=400,height=400)
p.removeAllModels()
for mol in allmol:
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:
# 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(reference)))
Heavy Atom RMS: 3.604550104229233e-07 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
the reference structure is not aligned with the rest, yet RMSD between the first mol and reference is very small (btw i wonder why not 0), as if earlier "Align" did not work.
# find "core":
m = Chem.AddHs(Chem.MolFromSmiles('O=C1NCCNC(=O)c2nc(C(=O)NCCNC(=O)c3nc1ccc3)ccc2'))
core = m.GetSubstructMatch(Chem.MolFromSmiles('n1ccccc1'))
# create conformers, optimize their geometries and create a list of all structures to be aligned:
allmol = []
allmol.append(reference)
confs = AllChem.EmbedMultipleConfs(m,10,randomSeed=1)
# open a file to store conformers - this is because i don't know how else i can getconformers as "Mol" and not as "Conf"
w = Chem.SDWriter("test.sdf")
for conf in confs:
_ = AllChem.UFFOptimizeMolecule(m,confId=conf)
energy=AllChem.UFFGetMoleculeForceField(m,confId=conf).CalcEnergy()
m.SetProp('ENERGY', '{0:.2f}'.format(energy))
w.write(m, confId = conf)
w.close()
# load the file
suppl = Chem.SDMolSupplier('test.sdf')
for s in suppl:
allmol.append(s)
print(allmol)
[<rdkit.Chem.rdchem.Mol object at 0x7fb5a0cf5990>, <rdkit.Chem.rdchem.Mol object at 0x7fb5a0d07cb0>, <rdkit.Chem.rdchem.Mol object at 0x7fb5a0d07c10>, <rdkit.Chem.rdchem.Mol object at 0x7fb5a0c8ef30>, <rdkit.Chem.rdchem.Mol object at 0x7fb5a0c8e850>, <rdkit.Chem.rdchem.Mol object at 0x7fb5a0c9c1c0>, <rdkit.Chem.rdchem.Mol object at 0x7fb5a0c9c210>, <rdkit.Chem.rdchem.Mol object at 0x7fb5a0c9c940>, <rdkit.Chem.rdchem.Mol object at 0x7fb5a0c9c760>, <rdkit.Chem.rdchem.Mol object at 0x7fb5a0c9c4e0>, <rdkit.Chem.rdchem.Mol object at 0x7fb5a0c9c440>]
# view:
p = py3Dmol.view(width=400,height=400)
p.removeAllModels()
for mol in allmol:
mb = Chem.MolToMolBlock(mol)
p.addModel(mb,'sdf')
p.setStyle({'stick':{'radius':'0.15'}})
p.setBackgroundColor('0xeeeeee')
p.zoomTo()
p.show()
again the reference structure is not aligned