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

OpenBabelでエネルギー計算ができるようになって、さらに便利になりました。そこで、超シンプルな配座探索プログラム(obconfopt)を作ってみました。

基本的には、ランダムに配座を発生させ、最安定な配座を探すというランダム探索法を採用しています。

　Step 1. 初期配座の入力
　↓
　Step 2. torsion角をランダムに設定し、新しい配座の生成。
　↓
　Step 3. Energy Minimization
　↓
　Step 4. 前に発生させた配座と異なっていれば保持する。
　↓
　Step 5. 終了判定

終了するまでは、Step2からStep5を繰り返します。
ここで、Step2はobrotateを、Step3はobenergyとobminimizeを、そしてStep4はobfitを参考にすることができます。簡単ですので、興味のある方は、作ってみてください。

例としてペンタンを使ってテストしました。


実用レベルで使うには、環構造を考慮したり、ヒューリスティックなアルゴリズムを採用したりとまだまだ改良する部分は多いですが、全くのゼロからobconfoptに相当するプログラムを書くのは、相当大変ですので、OpenBablの便利さを再認識しました。

Structure-CDKを使って分子を描画してみます。Structure-CDKではSwingPainterというとても便利なクラスがあります。

　JFrame frame = new JFrame();
　frame.setTitle("TestStructCDK");
　frame.setSize(400,400);
　Container con = frame.getContentPane();

　SwingPainter painter = new SwingPainter();
　painter.setAtomContainer(mol);
　con.add(painter);
　frame.setVisible(true);

こんな感じで表示されます。

CDKのRenderer2DModelクラスを使って、分子の描画を行ってみます。まずは、JPanelクラスを継承したTestRender2DPanelを作ります。

TestRender2DPanel.java


描画用のPanelができれば後は簡単です。ここでは、JFrameにgetContentPane()を使ってTestRender2DPanelを追加しています。

　JFrame frame = new JFrame();
　frame.setTitle("Render2DPanel");
　frame.setSize(400,400);
　Container con = frame.getContentPane();
　TestRender2DPanel r2dp = new TestRender2DPanel(mol);
　con.add(r2dp);
　frame.setVisible(true);

こんな感じで表示されます。

OpenBabelのOBAtomIeratorを使って、分子の中心座標を計算してみます。

vector3 calcMass(OBMol& mol){

　vector3 center(0,0,0);
　OBAtomIterator i;
　OBAtom *atom;

　for(atom = mol.BeginAtom(i);atom;atom = mol.NextAtom(i)){
　　center += atom->GetVector();
　}

　center /= (float) mol.NumAtoms();

　return center;
}

マクロ定義がされていますので、for文を次のように書くこともできます。

　FOR_ATOMS_OF_MOL(atom,mol)
　{
　　center += atom->GetVector();
　}

ライブラリによって、原子が配列の0から格納されていたり、1から格納されていたりするので、Iteratorを使うことによって、コーディングミスを減らすことができると思います。

OpenBabel2.1.0になって、エネルギー計算を筆頭に有益な機能がいくつか追加されているようです。

今回は、logPとTPSAを計算してみます。

ob_exp1.cpp


私の環境では、コンパイルを次のようにして行いました。
　
　$ g++ -I/usr/local/include/openbabel-2.0 -L/usr/local/lib ob_exp1.cpp -lopenbabel

実行

　$ ./a.out
LogP 1.3922
TPSA 20.23

あと、PDBをfasta形式で出力できるようになってました。

　$ babel -ipdb input.pdb -ofsa output.fsa

単純な変換ですが、便利だと思います。

SD Fileでは、[Data_Header]、[Data]の組み合わせで、任意のデータが付加されていたり、また、付加したりすることができます。
CDKでは、これらデータを全て取得したい場合は、
　
　Hashtable ht = mol.getProperties();

特定の[Data_Header]のみ取得したい場合は(ここでは、"descriptors_mw")、

　String data = (String)mol.getProperty("descriptors_mw");

さらに、Tile行を取得したい場合は、

　String data = (String)mol.getProperty("Title");

のようにして処理できます。

また、データを付加したい場合は、

　ht.put("TEST DATA","9999.9"); 　　　　
　mw.setSdFields(ht);

となります。


TestIteratingMDLReader.java


　

QuteMolですが、Real Time Ambient Occlusionを使って分子を描画することが特徴の１つのようです[IEEE Trans Vis Comput Graph. 2006 Sep-Oct;12(5):1237-44]。




確かに、マウスでグルグル回しても描画レベルを保っています。すごい。

CDK-Tavernaは、CDKをTavernaから利用するためのpluginです。通常はスクリプトで書くような処理を、ワークフローを使うことにより簡便に処理することができます。プログラム全体を視覚的に理解できるところが気に入っています。

CDK-Tavernaを使ってSD FileをSmilesに変換してみます。ワークフローは下のようになります。

出力はこんな感じです。


Tavernaの使い方は、とても分かりやすいので、複雑なワークフローの作成も楽だと思います。
それとTavernaからBiomobyやSoaplabなども使えるみたいです。

HOSE(Hierarchical Organization of Spherical Environments)コード[1]は、部分構造を表現するコードの１つです。ある原子を中心に球形の層を第1層、第2層、第3層...と考慮することによりコード生成を行います。

CDKでは、HOSECodeGeneratorクラスを用いることにより、HOSEコードを生成させることができます。この例では、Caffeic acidを用い、9番目の原子を中心に第3層までコード化しています。

　HOSECodeGenerator hcg = new HOSECodeGenerator();
　try {
　　String HOSE_CODE = hcg.getHOSECode(mol, mol.getAtomAt(9), 3);
　　System.out.println(HOSE_CODE);
　} catch (CDKException e) {
　　e.printStackTrace();
　}

出力(HOSEコード):
　C-2;=CC(C,=OO/=CC,,/)



HOSEコードは、C13 NMRスペクトルの予測に用いられていることで有名ですね。

[1] Anal. Chim. Acta. (1978) 103:355-365



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

VPythonを利用すれば、簡単にMD等のシミュレーションにおける計算過程をリアルタイムに可視化できるので、学習や教育目的には適していると思います。

この↓サイトでは、箱の中のボールの動きをシミュレーションしています。

　http://wwwmaths.anu.edu.au/comptlsci/Tutorial-Gas/tute-gas.html

単原子分子のMDなんかVPythonで可視化するにはお手軽でいいかも。

こちら↓のサイトも参考になりました。

　http://lurbano-5.memphis.edu/GeoMod/index.php/Main_Page


とりあえずは、分子描画用のクラスを作ってみようと思ってます。



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