C++ による RDKit 入門

RDKit は Python/C++ ベースのケムインフォマティクスソフトウェア環境です．
このページでは，C++ を使い，RDKit で分子を構築したり，構築された分子からいくつかの情報を取得するような，ごく基本的な方法を練習しています．
教材は，RDKit のドキュメント Getting Started with the RDKit in C++ およびサンプルプログラムです．

インフォメーション

ウェブサイト

RDKit: Open-Source Cheminformatics Software は Usefl Links からなる簡素なページです．そこで見つかる下記ドキュメントで勉強しました．

The RDKit Documentation：RDKit のドキュメントのトップページです
Getting Started with the RDKit in C++："The RDKit Documentation" の直下にある，RDKit を C++ で始めるためのドキュメントです
RDKit: Main Page："The RDKit Documentation" の直下にある，C++ API リファレンスのページです
rdkit/Docs/Book/C++Examples：GitHub にあるサンプルプログラムです．example1〜17 があります

目次（ページ内リンク）

留意すべき点
例１ RDKit::SmilesToMol クラスを使った分子の構築
例２ RDKit::MolFileToMol 関数を使った分子の読み込み
例3 RDKit::SDMolSupplier クラスを使った複数の分子の読み込み
例4 表記法と相互変換の例
例5 原子と結合へのアクセス

留意すべき点

上 2 個は，ビルドオプションとして与えるパスはここを設定する，というメモです．
下 4 個は，"Getting Started with the RDKit in C++" 冒頭の内容を改変して引用しています．

Debian パッケージのインストール先 (1)
librdkit-dev のヘッダファイルは，/usr/include/rdkit/ 以下にインストールされます．コンパイルの I オプションではここを指定すればよさそうです．

Debian パッケージのインストール先 (2)
librdkit のライブラリは，/usr/lib/libRDKit*** というファイル名でインストールされます．ビルドの l オプションでライブラリの名前がわからなくなったらここを探します．
librdkit1 (20253.1)なら "ls /usr/lib/libRDKit*.2025.03.1" とでもすれば，何というライブラリがインストールされているか判ります．

Python より C++ の方が API は不完全です
RDKit の高度な機能の多くは Python で開発されており，C++ へのバックポートは需要に応じておこなわれます．そのため，コンフォーマー間の RMS 差の計算（GetConformerRMSMatrix）など，C++ では利用できない機能もあります．

コンパイル
2018 年 3 月リリースの時点で，RDKit はモダン C++ を使用しています．これはリリース時点で C++17 までを意味します．Clang と gcc では，フラグ -std=c++17 を使用してください．

メモリリーク
RDKit では，多くの関数が分子へのポインタを返し，そのようなポインタを引数として受け取ります．std::shared_ptr や std::unique_ptr といったスマートポインタがメモリリークを防ぐのに役立ちます．

Python を使うか C++ を使うか
独自のプログラミングロジックを多用し，より複雑なことをおこなうのであれば，C++ の方がよりパフォーマンスが向上する可能性が高いでしょう．

例１ RDKit::SmilesToMol クラスを使った分子の構築

最初のプログラムでは，RDKit::SmilesToMol を使って分子を構築し，構築された分子から情報を取得します．
構築されるべき分子クラスは 2 個あります．RDKit::ROMol と RDKit::RWMol です。
ROMol は読み取り専用分子です（Read-Only）．分子を編集する必要がある場合は RWMol（Read-Write）を使います．
どちらも GraphMol.h で宣言されています。

ソースファイル

下のソースコードは，「インフォメーション」セクションで紹介した，GitHub にあるサンプルプログラムの example1.cpp を参考にして作成したものです．
RDKit::ROMol への共有ポインタを生成し，その原子数を端末に出力します．


#include <GraphMol/SmilesParse/SmilesParse.h>
#include <GraphMol/RWMol.h>
#include <iostream>
#include <memory>

int main() {
 std::string smiles = "Cc1ccccc1"; // トルエンの SMILES 形式

 RDKit::ROMol *mol_0 = RDKit::SmilesToMol(smiles);  //SmilesToMol は RDKit::ROMol*を返す
 std::shared_ptr<RDKit::ROMol> mol_1(mol_0);  //スマートポインタに所有権を移す

 std::shared_ptr<RDKit::ROMol> mol_2(RDKit::SmilesToMol(smiles));  //上 2 行と同じことを 1 行で

 std::cout << "原子数：" << mol_1->getNumAtoms() << std::endl;
 std::cout << "原子数：" << mol_2->getNumAtoms() << std::endl;

 return 0;
}

コード中の SMILES 形式 "Cc1ccccc1" について．
大文字の "C" は脂肪族炭素，小文字の "C" は芳香族炭素を表します，同じ数字（ここでは "1"）どうしは環を成していることを示します．
6 個の芳香族炭素が環を形成していて，それに脂肪族炭素が 1 個結合しているわけだから，メチルベンゼン == トルエンです．

スマートポインタを生成する定番の std::make_shared は，使わないようです．方法はあるのかもしれませんが，調べていません．

冒頭に示した C++ API リファレンスのページの RDKit::ROMol Class Reference を見るともう少し理解が深まります．
例えば，

RDKit::SmilesToMol()：SMILES 形式の文字列（ここでは "Cc1ccccc1"）を引数にして RWMol へのポインタ（生ポインタ）を返す
RDKit::ROMol::getNumAtoms()：ROMol クラスの公開関数．水素以外の原子数を unsigned int 型で返す

コンパイル

「留意すべき点」で記したように，

-I/usr/include/rdkit でヘッダファイルがインストールされたディレクトリを指定
C++17 を使用

すればコンパイルできます．すなわち，
~/wk$ g++ exp1.cc -c -I/usr/include/rdkit -std=c++17

リンク

l オプションで指定するライブラリ名は，「留意すべき点」で触たようにインクルードしているヘッダファイルから推定できそうです．
~/wk$ g++ exp1.o -o exp1 -lRDKitSmilesParse
とやってみたらリンクできました．ここでは exp1 という実行ファイルが生成します．

実行

実行ファイルのあるディレクトリで
~/wk$ ./exp1
とすると，下のように原子数 7 と出力されました．SMILES 形式から分子が構築できたようです．

原子数：7
原子数：7

例２ RDKit::MolFileToMol 関数を使った分子の読み込み

このセクションでは，MDL MOL 形式の分子ファイルを読み込んで，RDKit::ROMol インスタンスを構築する練習をします．
MOL 形式の分子ファイルには 1 個の分子が記述されています．

分子ファイル

読み込みに使ったファイルを下に示します．
エタンを Mol 形式で記述したファイルです．本サイトで開発している分子モデリングソフト Builcule で作成して，テキストエディタで 3 行めまでを編集しました．
3 行めまでは何か書いてあればよさそうです．ただし行数は 3 行でないと，後述のサンプルコードで RDKit::FileParseException という例外が発生しました．
ファイルは，作業ディレクトリ ~/wk に ethane.mol という名前で保存してあるとします．


エタン  //化合物名
 OpenBabel12152307473D  //作成したプログラム
 //コメント行
  8  7  0  0  0  0  0  0  0  0999 V2000
   -0.6290   -0.4450    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6290   -1.1030   -0.9310 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6290   -1.1030    0.9310 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.6290    0.4450    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.6290    1.1030   -0.9310 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.6290    1.1030    0.9310 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.5600    0.2140    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.5600   -0.2140    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  1  3  1  0  0  0  0
  1  4  1  0  0  0  0
  1  7  1  0  0  0  0
  4  5  1  0  0  0  0
  4  6  1  0  0  0  0
  4  8  1  0  0  0  0
M  END

ソースファイル

例１と同じくサンプルプログラムの example1.cpp を抜粋，改変したものです．
これを作業ディレクトリ ~/wk に exp2.cc という名前で保存するとします．


#include <iostream>
#include <GraphMol/GraphMol.h>
#include <GraphMol/FileParsers/FileParsers.h>

int main(int argc, char **argv) {
 std::string mol_file = "./ethane.mol";  //ファイルパス

 RDKit::ROMol *mol0 = RDKit::MolFileToMol(mol_file);  //RDKit::MolFileToMol 関数の戻り値は RWMol へのポインタ
 std::shared_ptr<RDKit::ROMol> mol1(mol0);  //スマートポインタに所有権を移す

 //std::shared_ptr<RDKit::RWMol> mol2(RDKit::MolFileToMol(mol_file));  //上 2 行と同じことを 1 行で

 //分子情報の出力
 std::cout << "Number of atoms(atomic number > 1) " << mol1->getNumAtoms(true) << std::endl;  //水素以外の原子数
 std::cout << "Number of atoms " << mol1->getNumAtoms(false) << std::endl;  //水素を含めた原子数
 std::cout << "Number of heavy atoms " << mol1->getNumHeavyAtoms() << std::endl; //水素以外の原子数
 std::cout << "Number of bonds(atomic number > 1) " << mol1->getNumBonds(true) << std::endl;  //水素が関与しない結合数
 std::cout << "Number of bonds " << mol1->getNumBonds(false) << std::endl;  //水素の関与も含めた総結合数

 return 0;
}

C++ API リファレンスから RDKit::MolFileToMol() の関数宣言を引用すると，
RWMol *RDKit::v1::MolFileToMol(const std::string &fName, bool sanitize = true, bool removeHs = true, bool strictParsing = true)
となっています．戻り値は RWMol * なので，読み書き可能です．生ポインタなので，コード内では共有ポインタに所有権を移しています．
引数の意味は下のとおりです．

const std::string &fName
bool ssanitize：true だと，芳香環を検出結合型を設定したり，結合次数が明示的でない場合に設定したりします
bool sremoveHs：true だと，分子が sanitize されている場合に限り，水素を除去します
bool sstrictParsing：false だと，分子の構造に関するパーサが緩くなるそうです

また，C++ API リファレンスの RDKit::ROMol Class Reference で目にした下記の関数を使ってみました．

getNumHeavyAtoms()
getNumBonds(bool onlyHeavy=1)

ビルド

~/wk$ g++ exp2.cc -c -I/usr/include/rdkit -std=c++17
~/wk$ g++ exp2.o -o exp2 -lRDKitGraphMol -lRDKitFileParsers

実行

~/wk$ ./exp2
とすると，下のように出力されました．
ファイルが正しく読み込まれたようです．

Number of atoms(atomic number > 1) 2
Number of atoms 8
Number of heavy atoms 2
Number of bonds(atomic number > 1) 1
Number of bonds 7

例3 RDKit::SDMolSupplier クラスを使った複数の分子の読み込み

このセクションでは，RDKit::SDMolSupplier クラスを使って，SDF 形式の分子ファイルを開く練習をします．
SDF 形式は，複数の MOL 形式の分子を一つのファイルにまとめたものです．
ファイルパスを引数にして RDKit::SDMolSupplier インスタンスを作成し，個々の分子を while ループで取得します．

分子ファイル

下が練習に使った SDF ファイルです．
アンモニアと水分子を別々に MOL 形式で作成し，テキストエディタで 1 ファイルにまとめたものです．複数の分子を区別するには，"$$$$" を区切り行とします．
ファイルは，作業ディレクトリ ~/wk に mols.sdf という名前で保存してあるとします．


ammonia
 OpenBabel12152314473D
 //コメント行
  4  3  0  0  0  0  0  0  0  0999 V2000
    0.0000    0.0000    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.0000   -0.6580   -0.9310 H   0  0  0  0  0  0  0  0  0  0  0  0
    0.0000   -0.6580    0.9310 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.9310    0.6580    0.0000 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  1  3  1  0  0  0  0
  1  4  1  0  0  0  0
M  END
$$$$
H2O
 OpenBabel12152314523D
 //コメント行
  3  2  0  0  0  0  0  0  0  0999 V2000
    0.3100    0.0000    0.3100 O   0  0  0  0  0  0  0  0  0  0  0  0
    0.3100   -0.6580   -0.6210 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6210    0.6580    0.3100 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  1  3  1  0  0  0  0
M  END

ソースファイル

下のコードは，サンプルプログラムの example2.cpp の一部を改変したものです．
ここでは，これを作業ディレクトリ ~/wk に exp3.cc という名前で保存します．


#include <iostream>
#include <GraphMol/GraphMol.h>
#include <GraphMol/FileParsers/MolSupplier.h>

int main(int argc, char **argv) {
 RDKit::SDMolSupplier mol_supplier("./mols.sdf");  //引数は SDF ファイルのパス

 //mol_supplier から 1 個ずつ分子を取得して原子数を表示
 std::unique_ptr<RDKit::ROMol> mol;
 while (!mol_supplier.atEnd()) {
  mol.reset(mol_supplier.next());
  std::cout << mol->getProp<std::string>("_Name") << " has "
    << mol->getNumAtoms(false) << " atoms." << std::endl;
 }

 return 0;
}

RDKit::SDMolSupplier には 3 種類のコンストラクタが記されていました．

SDMolSupplier()
SDMolSupplier(const std::string &fileName, const MolFileParserParams &params=MolFileParserParams())
SDMolSupplier(std::istream *inStream, bool takeOwnership=true, const MolFileParserParams &params=MolFileParserParams())

です．ここでは 2 番めの，const std::string & を第一引数とするコンストラクタを使っています．

また，std::unique_ptr<RWMol> を返す関数は 2 個ありました．

std::unique_ptr<RWMol> next() override
std::unique_ptr<RWMol> operator[](unsigned int idx)

です．ここでは std::unique_ptr<RWMol> next() を利用しています．

ビルド

リンクオプションに，-lRDKitRDGeneral も必要ですが，どこに書いてあったのか分からなくなりました．違う書き方があるのかもしれません．

~/wk$ g++ exp3.cc -c -I/usr/include/rdkit -std=c++17
~/wk$ g++ exp3.o -o exp3 -lRDKitGraphMol -lRDKitFileParsers -lRDKitRDGeneral

実行

~/wk$ ./exp3
とすると，下のように出力されました．
ファイルが正しく読み込まれたようです．

ammonia has 4 atoms.
H2O has 3 atoms.

例4 表記法と相互変換の例

サンプルプログラムの example3.cpp では，いくつかの表記法と相互変換が例示されています．
そのなかから抜粋して試してみます．

ソースコード

下にコードを示します．
新たに，RDKit::MolToSmiles() と RDKit::MolToMolBlock() という関数を使っています．


#include <iostream>
#include <GraphMol/FileParsers/FileParsers.h>
#include <GraphMol/SmilesParse/SmilesWrite.h>
#include <GraphMol/SmilesParse/SmilesParse.h>
#include <GraphMol/MolOps.h>

int main(int argc, char **argv) {
 //ピリジンを SMILES 形式の文字列から作成
 std::shared_ptr<RDKit::ROMol> mol1(RDKit::SmilesToMol("c1cccnc1"));  //出力するだけなので，RDKit::ROMol で充分
 std::cout << "ピリジン１：" << RDKit::MolToSmiles(*mol1) << std::endl;

 //ピリジンをケクレ化した（芳香環ではなく単結合と二重結合から成る）SMILES 形式の文字列から作成
 std::shared_ptr<RDKit::RWMol> mol2(RDKit::SmilesToMol("C1=CC=CN=C1"));  //後で編集（ケクレ化と名前の設定）するので，RDKit::RWMol 
 std::cout << "ピリジン２：" << (*mol2) << std::endl;

 //RDKit::RWMol インスタンスをケクレ化
 RDKit::MolOps::Kekulize(*mol2);
 std::cout << "ピリジン３：" << RDKit::MolToSmiles(*mol2) << std::endl;
 std::cout << std::endl;

 //ブロック化 
 mol2->setProp("_Name", "pyridine");
 std::cout << RDKit::MolToMolBlock(*mol2) << std::endl;

 return 0;
}

RDKit::MolToSmiles や RDKit::MolToMolBlock() の関数宣言は，RDKit: RDKit Namespace Reference にあります．


std::string RDKit::MolToSmiles(const ROMol &mol,  //分子
 bool doIsomericSmiles = true,  //SMILES に立体化学と同位体情報を含めるか
 bool doKekule = false,  //ケクレ化をするか
 int rootedAtAtom = -1,  //SMILES を指定された原子から始める
 bool canonical = true,  //SMILES を正規化するか
 bool allBondsExplicit = false,  //すべての結合を明示するか
 bool allHsExplicit = false,  //すべての水素を明示するか
 bool doRandom = false,  //SMILES 形式の最初の原子をランダムに選択するか
 bool ignoreAtomMapNumbers = false  //分子を正規化する際に SMILES 形式に与えられた原子マップ番号を無視するか
)


std::string RDKit::MolToMolBlock(const ROMol &mol,  //分子
 boolincludeStereo = true,  //立体化学情報を含めるか
 int confId = -1, //使用するコンフォーマーを選択．-1 は分子内で見つかる最初のコンフォーマー
 bool kekulize = true,  //ケクレ化をするか
 bool forceV3000 = false  //ForthV3000 形式で出力するか．これは高分子用の形式．原子数が 999 以上だと自動的に実行される
)

ビルド

~/wk$ g++ exp4.cc -c -I/usr/include/rdkit -std=c++17
~/wk$ g++ exp4.o -o exp4 -lRDKitSmilesParse -lRDKitFileParsers -lRDKitGraphMol -lRDKitRDGeneral

実行結果

~/wk$ ./exp4
とすると，下のように出力されました．

ピリジン１：c1ccncc1
ピリジン２：c1ccncc1
ピリジン３：C1=CC=NC=C1

pyridine
     RDKit          2D

  6  6  0  0  0  0  0  0  0  0999 V2000
    1.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7500   -1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.7500   -1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.5000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.7500    1.2990    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.7500    1.2990    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  2  0
  2  3  1  0
  3  4  2  0
  4  5  1  0
  5  6  2  0
  6  1  1  0
M  END

例5 原子と結合へのアクセス

構築した分子の各原子あるいは各結合にアクセスする練習をします．
分子，原子，結合と言ったクラスには多くのメンバ関数やメンバ変数があります．ここではそのごく一部を使います．API リファレンスへのリンクを張っておきます．

サンプルプログラムの example6.cpp からいくつかの例を抜粋し改変してコードを作成しました．

ソースコード


#include <iostream>
#include <memory>
#include <GraphMol/SmilesParse/SmilesParse.h>
//#include <GraphMol/GraphMol.h>  //example6.cpp ではインクルードしていたけど，省略してもビルドできました
//#include <GraphMol/AtomIterators.h>  //同上

int main() {
 //SMILES 形式の文字列からピリジンを作成
 std::shared_ptr<RDKit::ROMol> mol(RDKit::SmilesToMol("n1ccccc1"));
 //おまけ．ブタジエン
 //std::shared_ptr<RDKit::ROMol> mol(RDKit::SmilesToMol("C=CC=C"));

 std::cout  << "各原子の原子番号；";
 for(auto atom: mol->atoms())
  std::cout << atom->getAtomicNum() << " ";
 std::cout << std::endl;

 //インデックスでも各々の原子へのポインタを取得できる
 const RDKit::Atom *atom = mol->getAtomWithIdx(0);
 std::cout << "インデックス 0 の原子の原子番号：" << atom->getAtomicNum() << std::endl;

 std::cout  << "各結合のタイプ：";
 for(auto bond : mol->bonds())
  std::cout << bond->getBondType() << " ";
 std::cout << std::endl;

 //インデックスでも結合へのポインタを取得できる
 const RDKit::Bond *bond = mol->getBondWithIdx(0);
 std::cout << "インデックス 0 の結合は芳香族か：" << bond->getIsAromatic() << " " << std::endl;

 return 0;
}

ビルド

$ g++ exp5.cc -c -I/usr/include/rdkit -std=c++17
$ g++ exp5.o -o exp5 -lRDKitSmilesParse -lRDKitGraphMol

実行結果

~/wk$ ./exp5
とすると，下のように出力されました．芳香族の結合には，タイプ 12 が与えられています．

各原子の原子番号；7 6 6 6 6 6 
インデックス 0 の原子の原子番号：7
各結合のタイプ：12 12 12 12 12 12 
インデックス 0 の結合は芳香族か：1

ブタジエンの場合は下のように出力されました．ブタジエンの結合のタイプは 2 と 1 でした．芳香族ではないという判定です．

各原子の原子番号；6 6 6 6 
インデックス 0 の原子の原子番号：6
各結合のタイプ：2 1 2 
インデックス 0 の結合は芳香族か：0

パッケージ名(バージョン)	注
g++ (14.2.0)	GNU C++ コンパイラ
librdkit-dev (202503.1)	RDKit（Collection of cheminformatics and machine-learning software）の開発ファイル
librdkit1t64 (202503.1)	RDKit（Collection of cheminformatics and machine-learning software）の共有ライブラリ
libboost1.88-dev	Boost が提供するライブラリ
rdkit-doc(202503.1)	念の為にインストールしておきました．/usr/share/doc/rdkit/html/ 以下にインストールされます．RDKit のコピーのようです

C++ による RDKit 入門

インフォメーション

ウェブサイト

関連パッケージのインストール

目次（ページ内リンク）

留意すべき点

例１ RDKit::SmilesToMol クラスを使った分子の構築

ソースファイル

コンパイル

リンク

実行

例２ RDKit::MolFileToMol 関数を使った分子の読み込み

分子ファイル

ソースファイル

ビルド

実行

例3 RDKit::SDMolSupplier クラスを使った複数の分子の読み込み

分子ファイル

ソースファイル

ビルド

実行

例4 表記法と相互変換の例

ソースコード

ビルド

実行結果

例5 原子と結合へのアクセス

ソースコード

ビルド

実行結果