RDKitで配座を発生させてエネルギーを計算してみる

RDKit2012_09はしばらくβ版で確認していましたが正式版がリリースされました。

今回はMCSを実装した点が今までと変わった点のように思います。

残念ながらまだwin64bit用のバイナリは無いようです。取りあえず自分のMBAではビルドもできて動作も確認できました。

MCS関連はまた今度試すとして今回は、SDFから配座を発生させて、UFFによる最適化後のエネルギーを計算させるスクリプトを書きました。


#! usr/bin/python 
import sys 
from rdkit import Chem 
from rdkit.Chem import AllChem 
sdf = Chem.SDMolSupplier( sys.argv[1] )
mols = [ Chem.AddHs(m) for m in sdf ] 
out = Chem.SDWriter("genconf.sdf")
 
for mol in mols:
    ids = AllChem.EmbedMultipleConfs(mol, numConfs = 10, clearConfs = True,
          pruneRmsThresh = 0.5)
     for cid in ids:
         AllChem.UFFOptimizeMolecule(mol, maxIters = 1000 ,confId = cid)
         ff = AllChem.UFFGetMoleculeForceField(mol, confId = cid)
         ff.Minimize()
         energy=ff.CalcEnergy()
         mol.SetProp("UFF_Energy", "%s"%energy)
         mol.SetProp("confId", "%s" % cid)
         out.write(mol,confId = cid)
 
out.close()

AllChem.EmbedMultipleConfsで分子に対して10個の配座を発生させます。

AllChem.UFFOptimizeMoleculeにて構造を最適化します(1000回に設定)0なら最適化が終わっていて

1、がかえってくる場合は1000以上の数値を設定した方がいいと思われますが今回次のステップにはここのクライテリアを抜いてしまいました。

ff = AllChem.UFFGetMoleculeForceField(mol, confId = cid)
ff.Minimize()
energy=ff.CalcEnergy()
にて最適化後のエネルギーを受け取ります。

今回は出力を全部にしたのですが、最低のものだけを取り出したいのです、、、

リストのminでは数値が出てくるけどインデックスがとれないのでどうしたらいいんだろう。

list.index(min_energy)で行けそうだから後で追加しよう。

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s