ちょっと改変。

SDFには最小のエネルギーを持つ配座だけを
テキストには全部のログを残すようにしてみた。

#! usr/env/python
import sys
from rdkit import Chem
from rdkit.Chem import AllChem
sdf = Chem.SDMolSupplier(sys.argv[1])
name = "conf_"+str(sys.argv[1]).rsplit(".sdf")[0]+".txt"
f = open( name , "w")
f.write("molid\tconfId\tmin_energy\n")

mols = [ Chem.AddHs(m) for m in sdf ] 

out = Chem.SDWriter("genconf.sdf")

for molid, mol in enumerate(mols):
    ids = AllChem.EmbedMultipleConfs(mol, numConfs = 10, clearConfs = True,
            pruneRmsThresh = 0.5)
    e_list = []
    for cid in ids:
        AllChem.UFFOptimizeMolecule(mol, maxIters = 1000  ,confId = cid)
        ff = AllChem.UFFGetMoleculeForceField(mol, confId = cid)
        ff.Minimize()
        energy=ff.CalcEnergy()
        e_list.append(energy)
        f.write("%s\t%s\t%s\n"%(molid, cid, energy))
    min_cid = e_list.index(min(e_list))
    #print min(e_list)
    mol.SetProp("molid","%s"%molid)
    mol.SetProp("Min_energy","%s" % min(e_list))
    mol.SetProp("confId","%s"%min_cid)
    out.write(mol, confId = min_cid)
out.close()
f.close()


lioni$ time python conformation_v2.py mol100.sdf 
real	1m36.280s
user	1m36.107s

sys	0m0.121s

lion$ wc conf_mol100.txt 
 1001  3003 18811 conf_mol100.txt

lioni$ head -21 conf_mol100.txt 

molid	confId	min_energy
0	0	64.3825031596
0	1	66.8213267669
0	2	62.302429729
0	3	66.0064341133
0	4	64.3867981839
0	5	67.427604423
0	6	67.2825997948
0	7	67.2140755789
0	8	67.9395352633
0	9	60.378773236
1	0	72.1261194316
1	1	67.5461349343
1	2	73.0026727272
1	3	69.8917224767
1	4	65.6810320174
1	5	70.1950915266
1	6	70.2749205787
1	7	69.8981659667
1	8	72.9899634079
1	9	69.898166115

ちゃんとテキストの方には100分子で10配座のデータがきました。
時間がちょっとかかるのでもっと荒くていいならiter =200くらいでいいかなぁ。

広告

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中