Open Babel ライブラリ練習ノート

Open Babel のライブラリ版である libopenbabel を使う練習をしました.
これはすなわち,自作ソフトウェアに Open Babel の機能を組み込むことが可能になるということです.
このページでは,ファイル形式の変換,水素付加,および分子力学計算の機能について,C++ 言語でサンプルコードを作成し,動作確認しています.
C/C++ なら system() 関数を使って実行ファイルを呼び出すという方法もありますが,ライブラリからしか呼び出せない関数を使うかどうかが分岐点だと思います.

このページでは,下記の順序で OpenBabel オブジェクトを利用しています.

  1. OpenBabel::OBConversion:ファイル形式の変換を担うクラスです.コンストラクタ引数として入出力ファイルストリームを与えています
  2. OpenBabel::OBMol:水素付加など編集用のクラスです.OBConversion で読み取った分子構造を OBMol に渡しています
  3. OpenBabel::OBForceField:分子力学計算用のクラスです.エネルギー極小化とコンフォメーション探索をおこないます.入出力構造として OBMol を渡しています

分子モデリングソフト Builcule のドキュメントOpenBabel による分子力学計算は,このページで紹介する機能とファイル出力の場合分け(系全体か分子ごとか)を組み合わせて実装しています.


インフォメーション

目次(ページ内リンク)


練習 1:ファイル形式の変換
練習 2:水素付加
練習 3:エネルギー極小化
練習 4:コンフォメーション探索

ウェブサイト

下は参考にしているページで,開発者向けのチュートリアルページや API のページにサンプルコードが掲載されています.
そのサンプルコードを改変するところから練習を始めます.


Developer:Cpp Tutorial - Open Babel
Open Babel: Getting Started(Open Babel: API Documentation 内のページ)

Debian パッケージ

このページは,Debian GNU/Linux バージョン 11 "bullseye" で作成しました.
バージョン 12 "bookworm" では.libopenbabel-dev のバージョンは変わらなかったので,内容はほぼそのままにしてあります.
下に使用している主なパッケージを示します.

パッケージ名(バージョン)
g++ (12.2.0)GNU C++ コンパイラ
libopenbabel-dev (3.1.1)Open Babel のヘッダファイルと静的ライブラリのパッケージ(共有ライブラリのパッケージも同時にインストールされます)
libopenbabel-doc (3.1.1)Open Babel の API ドキュメント./usr/share/doc/libopenbabel-dev/html/ にインストールされます

練習 1:ファイル形式の変換

ソースコード

"Developer:Cpp Tutorial - Open Babel" の最初のコードを改変したものを下に記します.例外処理等を削除して,入出力のファイル形式を変更しています.

Open Babel といえばこれ,というコマンド,
obabel [-i<input-type>] <infilename> [-o<output-type>] -O<outfilename>
を使うプログラムです.OBConversion というクラスを使います.

  1. 入力用と出力用のファイルストリームへのポインタを引数にしてファイル形式変換用クラスオブジェクト OBConversion をインスタンス化
  2. OBConversion のメンバ関数 SetInAndOutFormats() 関数で,入力用と出力用のファイル形式を設定
  3. OBConversion のメンバ関数 Convert() を呼ぶ

#include <openbabel/obconversion.h>

int main(int argc, char **argv) {
 std::ifstream ifs(argv[1]);  //ここでは第一引数を入力ファイル名とし,その名前で入力ファイルストリームを作成する
 std::ofstream ofs(argv[2]);  //ここでは第二引数を出力ファイル名とし,その名前で出力ファイルストリームを作成する

 OpenBabel::OBConversion conv(&ifs, &ofs);  //ファイル形式変換用のクラスオブジェクト.引数は,入力および出力ファイルストリーム
 conv.SetInAndOutFormats("XYZ","PDB");  //入出力ファイル形式を設定.第一引数が入力,第二引数が出力.変数にする場合は const char *.
 conv.Convert();  //変換.ここでは破棄しているけれど,変換したファイルの数が戻り値

 return 0;
}


これと同様のコードを,当サイトで開発中の Builcule の中に記述したら,ファイルに原子の情報が出力されませんでした.
"Open Babel: Getting Started" に例示されているように,ファイル入出力には OpenBabel::OBMol オブジェクトを使います.
すなわち,
OpenBabel::OBConversion conv(&ifs, &ofs);
conv.SetInAndOutFormats("XYZ", "PDB");
OpenBabel::OBMol mol;
conv.Read(&mol);
conv.Write(&mol);

ビルドと実行

上のコードを main.cc という名前で保存しました.

コンパイル

ヘッダファイルがインストールされたディレクトリ /usr/include/openbabel3/ を -I オプションで与えています.
~$ g++ -c main.cc -I/usr/include/openbabel3/

リンク

生成する実行ファイルの名前を "babe" としています.
~$ g++ -o babe main.o -lopenbabel

実行

第一引数は変換元の(すでに存在する)ファイルの名前,第二引数は変換先の(これから作成する)ファイルの名前です.
methane.xyz は,Builcule で作成しました.
~$ ./babe methane.xyz methane.pdb
とすると,下に示す入力ファイルから出力ファイルが生成しました.成功したようです.

入力ファイル:methane.xyz

5
This file was creatwd with Builcule.
C         -0.00000        0.00000       -0.00000
H          1.07000        0.00000       -0.00000
H         -0.35667        0.00000        1.00881
H         -0.35667       -0.87365       -0.50440
H         -0.35667        0.87365       -0.50440

出力ファイル:methane.pdb

COMPND    This file was creatwd with Builcule.
AUTHOR    GENERATED BY OPEN BABEL 3.1.1
HETATM    1  C   UNL     1       0.000   0.000   0.000  1.00  0.00           C
HETATM    2  H   UNL     1       1.070   0.000   0.000  1.00  0.00           H
HETATM    3  H   UNL     1      -0.357   0.000   1.009  1.00  0.00           H
HETATM    4  H   UNL     1      -0.357  -0.874  -0.504  1.00  0.00           H
HETATM    5  H   UNL     1      -0.357   0.874  -0.504  1.00  0.00           H
CONECT    1    4    5    2    3
CONECT    2    1
CONECT    3    1
CONECT    4    1
CONECT    5    1
MASTER        0    0    0    0    0    0    0    0    5    0    5    0
END

練習 2:水素付加

このセクションでの作業は,Open Babel のコマンドで書くと,
~$ obabel -ipdb in.pdb -O inh.pdb -p 7.4
と同等です.

Protein Data Bank などのデータバンクからタンパク質などの構造をダウンロードする場合,その分子構造が X-線結晶回折法で決定されたものであれば水素原子が含まれません.
このようなときには必要に応じて,水素原子を編集する手段が必要となります.
libopenbabel の OBMol クラスの説明を見ると,水素原子を編集する関数がいくつか用意されていることが判ります.
下に列挙します.

ソースコード

練習 1 と同じく例外処理等は省いています.

ここでは上記のなかから,
bool AddHydrogens(bool polaronly = false, bool correctForPH = false, double pH = 7.4)
関数を選択して練習しています.

ファイル入出力には上と同じく OBConversion インスタンスを利用します.
練習 1 と違う点は,ファイルの入出力に OpenBabel::OBConversion::Read(OpenBabel::OBMol*) と OpenBabel::OBConversion::Write(OpenBabel::OBMol*) という関数を使っている点です.


#include <openbabel/mol.h>
#include <openbabel/obconversion.h>

int main(int argc,char **argv) {
 std::ifstream ifs(argv[1]);  //第一引数を入力ファイル名とし,入力ファイルストリームを作成する
 std::ofstream ofs(argv[2]);  //第二引数を出力ファイル名とし,出力ファイルストリームを作成する

 OpenBabel::OBConversion conv(&ifs, &ofs);
 conv.SetInAndOutFormats("pdb","pdb");  //入出力ファイル形式は PDB 形式に固定
 OpenBabel::OBMol mol;

 conv.Read(&mol);  //OBConversion でファイルを開いて,OBMol に分子構造のデータを渡している
 mol.AddHydrogens(false, true, 7.4);  //OBMol 内で水素付加.pH による影響を補正して,すべての水素を付加
 conv.Write(&mol);  //OBConversion は OBMol から水素付加されたデータを入手してファイルに出力している

 return 0;
}

ビルドと実行

libbabel で WHN に水素付加

ビルドは上と同じです.
~$ g++ -c main.cc -I/usr/include/openbabel3/
~$ g++ -o babe main.o -lopenbabel

まず Builcule で水素を有しない XYZ 形式のトリペプチド Trp-His-Asn を作成し,in.pdb として保存しました.

実行コマンドは,
~$ ./babe in.pdb inh.pdb です. in.pdb に対して,水素付加された構造が,inh.pdb として出力されます.

画像は OpenBabel で水素付加したペプチドを,Builcule で表示したものです.
この構造と中性でのペプチドの構造とでは多少差異があります.

水素付加に関しては,Open Babel をプログラムに組み込む場合には,入力構造に注意する必要がありそうです.


練習 3:エネルギー極小化

Obminimize - Open Babel にコマンドの解説があります.
このセクションでの作業は,Open Babel のコマンドで書くと,
~$ obminimize -n 10000 -cg -c 1e-6 -ff UFF inh.pdb > out.pdb
と同等です.
言葉で書くと,obminimize コマンドを使い,最大ステップ数 1000,極小値探索アルゴリズムは共役勾配法,収束条件 1e-6,力場 UFF(Universal force field)の条件で inh.pdb のエネルギーを極小化せよ」となります.
作成時点でファイル名を指定して保存するオプションを見つけられなかったので,標準出力をリダイレクトしています.

力場については,Molecular Mechanics and Force Fields — Open Babel 3.0.1 documentation に,「有機分子には Generalized Amber Force Field (GAFF) または MMFF94 Force Field (MMFF94) を,その他の種類の分子には Universal Force Field (UFF) を使用することをお勧めします.」と書いてあります.ここでは UFF を使ってみました.

Builcule で上のセクションで作成した inh.pdb の誤り 3 箇所を修正し,改めて in.pdb としました.

ソースコード

Open Babel の API ドキュメント(libopenbabel-doc パッケージ)の OBForceField Class Reference ページに,プログラムの断片があります.
それと上記プログラムを組み合わせて下のコードを作成しました,


#include <openbabel/forcefield.h>
#include <openbabel/mol.h>
#include <openbabel/obconversion.h>
#include <iostream>

int main(int argc,char **argv) {
 std::ifstream ifs(argv[1]);
 std::ofstream ofs(argv[2]);

 OpenBabel::OBConversion conv(&ifs, &ofs);
 conv.SetInAndOutFormats("pdb","pdb");
 OpenBabel::OBMol mol;

 conv.Read(&mol);

//力場を発生.ここでは UFF 力場を使用
 OpenBabel::OBForceField* pFF
     = OpenBabel::OBForceField::FindForceField("UFF");

//力場に分子を置く.パラメータが見つからないときなどには,その旨が表示される(これは,UFF ではなく MMFF94 を使ってみると確認できる)
 pFF->SetLogFile(&std::cerr);
 pFF->SetLogLevel(OBFF_LOGLVL_LOW);  //表示させるのは LOGLVL_LOW だけど
 //怖いので,今回は 1 箇所例外処理を記述
 if (!pFF->Setup(mol))
  std::cerr << "ERROR: could not setup force field." << std::endl;

//勾配降下法でエネルギー極小化.10000 ステップ経過しても収束しなければ終了(ここの例では 1600 ステップあまりで収束しました)
 pFF->ConjugateGradients(10000);

 //エネルギー極小化後の構造を OpenBabel::OBMol に渡す
 pFF->GetCoordinates(mol);
 //OpenBabel::OBMol の構造を OpenBabel::OBConversion に渡してファイル出力
 conv.Write(&mol);

 return 0;
}

ビルドと実行

libbabel で WHN をエネルギー極小化

ビルドは上と同じです.
~$ g++ -c main.cc -I/usr/include/openbabel3/
~$ g++ -o babe main.o -lopenbabel

出力ファイル名を out.pdb とします.
実行コマンドは,
~$ ./babe inh.pdb out.pdb
です.

エネルギー極小化前後の構造を Builcule で開きました.
青色がエネルギー極小化前,赤色が極小化後です.


練習 4:コンフォメーション探索

このセクションでの作業は,Open Babel のコマンドで書くと,
~$ obminimize -n 10000 -cg -c 1e-6 -ff MMFF94 dp.pdb > tmp.pdb #初期構造のエネルギー極小化
~$ obconformer 250 100 tmp.pdb > tmp2.pdb #コンフォメーション探索本体
~$ obminimize -n 100 -cg -c 1e-6 -ff MMFF94 tmp2.pdb > dp2.pdb #得られたコンフォーマーのエネルギー極小化
~$ rm tmp.pdb tmp2.pdb #中間ファイルの削除
とおそらくほぼ同等です.

練習3 のコードを少し書き換えるだけで,安定なコンフォメーションを探索するコードに変更できます.
基本的には,OpenBabel::OBForceField クラスの関数 pFF->WeightedRotorSearch() を呼び出すだけです.
変更した箇所にコメントをつけました.


#include <openbabel/forcefield.h>
#include <openbabel/mol.h>
#include <openbabel/obconversion.h>
#include <iostream>

int main(int argc,char **argv) {
 std::ifstream ifs(argv[1]);
 std::ofstream ofs(argv[2]);

 OpenBabel::OBConversion conv(&ifs, &ofs);
 conv.SetInAndOutFormats("pdb","pdb");
 OpenBabel::OBMol mol;

 conv.Read(&mol);

 //力場は, obconfomer() と同じ MMFF94 とする
 OpenBabel::OBForceField* pFF
     = OpenBabel::OBForceField::FindForceField("MMFF94");
 pFF->SetLogFile(&std::cerr);
 pFF->SetLogLevel(OBFF_LOGLVL_LOW);
 if (!pFF->Setup(mol))
  std::cerr << "ERROR: could not setup force field." << std::endl;

 //初期構造のエネルギーを極小化
 //ドキュメントに「最良の結果を得るためにはこの構造を最小にする必要があります」と書いてあるので記述
 //なお,obconformer() ではこの処理をおこなっていないので,手動では obconformer() を呼び出す前に obominimize() を実行せよ,ということらしい
 pFF->ConjugateGradients(10000);

 //コンフォメーション探索本体.引数は発生させるコンフォーマー数,エネルギー極小化のステップ数.Obconformer - Open Babel に書いてある値をそのまま使用
 pFF->WeightedRotorSearch(250, 100);

 //念のため得られたコンフォーマーのエネルギーを極小化.ステップ数は,obconformer と同じ
 pFF->ConjugateGradients(100);

 pFF->GetCoordinates(mol);

 conv.Write(&mol);

 return 0;
}

ビルドと実行

libbabel で リン酸無水物のコンフォメーション探索

ビルドは上と同じです.
~$ g++ -c main.cc -I/usr/include/openbabel3/
~$ g++ -o babe main.o -lopenbabel

入力ファイル名を dp.pdb,出力ファイル名を dp2.pdb とします.
実行コマンドは,
~$ ./babe dp.pdb dp2.pdb
です.

リン酸の無水物を試料として動作確認をおこないました.
コンフォメーション探索前後の構造を Builcule で開きました.
上が入力構造,下が得られた安定コンフォーマーです.

上では MMFF94 力場を使っていますが,力場を変更すると得られる安定コンフォーマーも異なります,


参考書の検索