ケムインフォマティクスに虚空投げ

Energy Minimizationのパラメータはsteps以外にもstep_sizeとconvergenceがあります。それぞれdefault値は以下のとおりです。

step_size = 2 pm
convergence = 0.01 KJ/mol/nm

MMTKでは単位の記述はUnitsモジュールを利用すると便利です。

　from MMTK import Units

具体的にEMのパラメータ設定は、以下のように記述することができます。

　minimizer(steps = 200, step_size = 0.04*Units.Ang, convergence = 0.02*Units.kJ/(Units.mol*Units.nm))



人気ブログランキング(クリックして応援してね)

MMTKを使ってEnergy Minimizationを行ってみます。以下にそのコードを示します。

EM.py


それでは、順番にみていきます。

UniverseはMMTKにおけるシミュレーションの基盤となるクラスです。topology、boundary condition, force fieldなどの設定を行うことができます。

　universe = InfiniteUniverse(Amber94ForceField())

分子(タンパク質)の読み込みは、Proteinクラスで行いました。'bala1'は、MMTKにはじめから含まれているもので、その実体は、$Python23\Lib\site-packages\MMTK\Database\PDB\bALA1.pdbです。ファイル名を指定してPDBを入力することもできます。

　universe.protein = Protein('bala1')

今回は、Energy Minimizationを2つの最適化手法で実施しています。
まずは、SteepestDescent(SD)です。StandardLogOutputでPotential Energy等の情報が50step毎に表示されます。SDのステップ数は200です。

　minimizer1 = SteepestDescentMinimizer(universe, actions = 　　[StandardLogOutput(50)])
　minimizer1(steps = 200)

次に、ConjugateGradientで最適化です。SDと使い方はほぼ同じです。

　minimizer2 = ConjugateGradientMinimizer(universe, actions = [StandardLogOutput(50)])
　minimizer2(steps = 100)

最後に、EM後の構造を表示します。拡張子(.pdb)と関連付けされたソフトウェアが起動し、構造を表示させます。

　view(universe)

プログラムを実行すると、以下のような結果が表示されると思います。

SteepestDescent
Step 0
Potential energy: -80.205828, Gradient norm: 809.115225

Step 50
Potential energy: -110.589675, Gradient norm: 62.273830

Step 100
Potential energy: -111.893701, Gradient norm: 64.541604

Step 150
Potential energy: -112.424998, Gradient norm: 67.033451

Step 200
Potential energy: -112.810568, Gradient norm: 62.190661

ConjugateGradient
Step 0
Potential energy: -112.810568, Gradient norm: 62.190661

Step 50
Potential energy: -116.103856, Gradient norm: 16.267230

Step 100
Potential energy: -116.504514, Gradient norm: 21.127653




人気ブログランキング(クリックして応援してね)

CDKのAPIドキュメントの構築方法をメモしたいと思います。といいましても、たった一行なのですが....　ただ、antは事前にインストールする必要があります。

CDKのルートディレクトに移動して、以下のコマンドを実行します。

　> ant -f javadoc.xml

実行後、doc\api以下にHTMLファイルが生成されています。index.htmlを見るといつもCDKのサイトで見ているAPIドキュメントが構築されています。



ネット上で閲覧できない古いバージョンのAPIを調べたいときにも役にたちますね。



人気ブログランキング(クリックして応援してね)

OpenBabelのWindows GUI版を使ってみました。なかなか使い勝手もよく便利だと思います。そういえば、BabelWinというGUI版のbabelもありましたね。私はコマンドライン版に慣れきってしまっているため使うことはほとんどないですが、GUI版は、バイナリで配布されているため、インストール作業も必要なく便利なのは確かだと思います。





人気ブログランキング(クリックして応援してね)

The Molecular Modelling Toolkit (MMTK)は、分子モデリングに関連する計算を行うためのPythonの拡張ライブラリです。MMTKは私にとって、Pythonを覚えるきっかになったライブラリです。ここ3-4年全く使っていなかったのですが、最近?Windows版とMacOS X版のbinary installersが公開されたようです。はじめはWindows版と言ってもCygwin上で動くのかなと思ったのですが、そうではありませんでした。MMTKは、MD,MMを手軽に体験できる便利なツールですので、興味のある方はぜひ使ってみてください。

それでは、Windows(XP)にMMTKをインストールします。

1. Pythonのインストール
　Windows版のPython 2.3.4日本語環境用インストーラを利用します。

2. pywin32インストール
　pywin32-209.win32-py2.3.exeをインストールします。

3. netCDFのインストール
　NetCDF 3.6.1のprecompileされたWindows DLLを取得し、展開後、ファイルをC:\WINDOSW\SYSTEM32以下にコピーします。

4. Numericのインストール
　Numeric-24.2.win32-py2.3.exeをインストールします。

5. ScientificPythonのインストール
　ScientificPython-2.4.6.win32-py2.3.exeをインストールします。

6. MMTKのインストール
　やっとたどり着きました。MMTK-2.4.4.win32-py2.3.exeをインストールします。

注意すべきところは、各ソフトウェアのバージョンで、組合せによってはMMTKが動作しない場合があります。今のところは、私の環境では、上記バージョンの組合せで動作しています。MMTKのサイトにExamplesがありますので、まずはこれを利用して動作確認することができます。



人気ブログランキング(クリックして応援してね)

これまでは、emacsやIDLEを利用してPythonのプログラムを書いていましたが、EclipseのプラグインであるPyDevが便利だと聞いたので試してみました。

�@インストール
　http://sourceforge.net/projects/pydev/ からorg.python.pydev.feature-1_2_2.zipを取得して、適当な場所で展開します。そして、featuresとplugins下のファイルをeclipseのfeaturesとplugins下にコピーします。

�Aインストールの確認
　eclipseを起動し、[Help]→[About Eclipse SDF]→[Plug-in Details]でpydevがインストールされているか確認できます。

�B設定
　[Window]→[Preferences]で左側のTreeから[Pydev]→[Interpreter Python]を選択します。そして、右上の[New]ボタンを押してPythonの実行ファイルの場所を設定します。あとは、通常通り[File]→[New]→[Project]でPydev Projectが作成できます。

試しにfrownsをimportして使ってみました。



おお！コード補完されました。



人気ブログランキング(クリックして応援してね)

MOBYのクライアントの1つであるRemoraを使ってみました。通常では様々なWebサービスを個々に使って解析しなければならないタスクを、下図のようなワークフローを構築することにより、一度に解析することができます。



MOBYのように複数のWebサービスを連携して利用できる仕組みがあるとRemoraのFrequently Asked Workflows (FAW)のように利用頻度の高いワークフローをテンプレートとして再利用しやすいですし、研究の効率化にも繋がると思います。

ただ、サイトを見る限りでは、MOBYを使いこなすのは、私には大変そうだなという印象を受けました。日本にもユーザーグループが設立されているようですし、使いこなせれば便利なのは間違いないと思いますので、これから学んでいこうと思ってます。



人気ブログランキング(クリックして応援してね)

SDF形式では、'M　END'と'$$$$'の間に様々な情報が付加されていることがあります。例えば以下のような情報が挙げられます。

....
M END
＞ ＜CAS_Number＞
XXXXXX

＞ ＜Catolog_Number＞
XXXXXX

＞ ＜logP＞
XXX
$$$$

CDKではこれら付加情報も分子情報と共に読み込んでいますので、次のような方法で具体的な値を得ることができます。

　　fr = new FileReader(new File(filename));
　　imr = new IteratingMDLReader(fr,
　　　　　DefaultChemObjectBuilder.getInstance());
　　while (imr.hasNext()){
　　　mol = (IMolecule)imr.next();
　　　String mol_id = (String)mol.getProperty("CAS_Number");
　　　String cat_num = (String)mol.getProperty("Catalog_Number");
　　　String logP = (String)mol.getProperty("logP");
　　　System.out.println(mol_id+" "+cat_num+" "+logP);
　　}

内部的にはHashtableに'CAS_Number'はキー、'XXXXXX'は値として格納されています。getProperties()メソッドを用いれてば、一度に全てのキーと値の情報を取得することができ便利です。

　Hashtable ht = (Hashtable)mol.getProperties();
　System.out.println(ht);



人気ブログランキング(クリックして応援してね)

CDKにおけるMDL形式(mol,sdf)ファイルの入力方法をメモしたいと思います。CDKのバージョンは、cdk-20060718を用いて行います。

まずは、１分子のみ(mol形式)の入力方法です。分子の入力確認をするために、とりあえずは、getAtomCount()メソッドで原子数を表示させています。
　
　fr = new FileReader(new File(filename));
　mr = new MDLReader(fr);
　mol = (IMolecule)mr.read(new Molecule());
　//原子数の表示
　System.out.println(mol.getAtomCount());　

次は、複数の分子が入っているファイル(sdf形式)から一分子づつ入力する方法です。

　fr = new FileReader(new File(filename));
　imr = new IteratingMDLReade(fr,
　　　　　　　　DefaultChemObjectBuilder.getInstance());
　while (imr.hasNext()){
　　mol = (IMolecule)imr.next();
　　//原子数の表示
　　System.out.println(mol.getAtomCount());
　}

一分子づつ入力しますので、大規模なsdf形式の入力に適しています。

最後は、sdf形式を一度に入力する方法です。

　fr = new FileReader(new File(filename));
　mr = new MDLReader(fr);
　chemfile = (IChemFile)mr.read(new ChemFile());　←�@
　//原子数の表示
　for(int i=0;i＜chemfile.getChemSequenceCount();i++){
　　IChemSequence cs = chemfile.getChemSequence(i);
　　　for(int j=0;j＜cs.getChemModelCount();j++){
　　　　IChemModel cm = cs.getChemModel(j);
　　　　ISetOfMolecules sms = cm.getSetOfMolecules();
　　　　　for(int k=0;k＜sms.getMoleculeCount();k++){
　　　　　　IMolecule mol = sms.getMolecule(k);
　　　　　　System.out.println(mol.getAtomCount());
　　　　　}
　　　}
　}

上記ソースの�@の行のみで複数の分子(sdf形式)を一度に入力しています。入力された各分子の原子数を知るためには、"//原子の表示" 以下の作業が必要となります。ChemSequenceクラス, ChemModelクラス, Moleculeクラスの関係が理解できると思います。



人気ブログランキング(クリックして応援してね)

先日の記事では、隣接・結合行列を取得するためにAdjacencyMatrixクラスとConnectionMatrixクラスを用いました。CDKでは、これらクラスを利用しなくても簡単に隣接・結合行列を作成できます。ソースコードを載せておきますので、メソッドの使い方等の参考になればと思います。

Matrix.java




人気ブログランキング(クリックして応援してね)