OpenBabelを使ってSMILES表記の分子式から3次元の分子構造を復元する方法

# はじめに

こんにちは，(株)日立製作所 研究開発グループ サービスコンピューティング研究部の露木です。

[前回の記事](https://qiita.com/MasafumiTsuyuki/private/73af68e94b14447a3d74)では，材料開発分野へのAI適用を念頭に[OpenBabel](http://openbabel.org/docs/current/Introduction/intro.html)プログラムをインストールする方法を紹介しました。今回は，そのOpenBabelを使ってSMILES表記の分子式から3次元の分子構造を復元する方法を説明します。

SMILESは情報化学 (ケモインフォマティクス，Chemoinfomatics) の分野で頻繁に用いられる分子構造の表記方法です。SMILESの特徴はたった1行の文字列で，結合を含む分子構造を表現できることです。例えば下図のような構造を持つフェノール (PhOH) 分子を，SMILESで表現すると `Oc1ccccc1` になります。しかし，SMILESには結合角や結合距離などの情報は含まれませんから，可視化やシミュレーション計算に利用するためにはカーテシアン座標やz-matrix表記の3次元分子構造を復元しなければなりません。

![out2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541219/a70684bc-0e96-7233-dbe0-f3b23dbd86ec.png)

今回説明する方法で計算した分子構造は機械学習の特徴量として用いるほか，情報化学分野のシミュレーション用の各種プログラムへの入力ファイルに利用することができます。

# 手順
## OpenBabelのインストール

最新版のOpenBabelを使いたい場合は，前回の記事で説明したようにソースコードからコンパイルしてください。だたし，今回の記事で利用する機能はある程度古い，Ubuntuのパッケージからインストールできるバージョンでも利用可能です。まだOpenBabelをインストールしていない場合は，以下のコマンドでインストールしてください。

```bash
$ sudo apt-get install openbabel
```

## OpenBabelの対応ファイル形式

もともと，OpenBabelは乱立する化学系シミュレーションプログラムの入力・出力ファイル形式を相互変換するためのツールです。OpenBabelが対応するファイル形式の一覧は `-Hall` オプションで表示できます。例えば，OpenBabel Version 2.3.2 では今回目的とするSMILESとz-matrixを含む，以下のファイル形式に対応しています。

```bash
$ obabel -Hall | grep -- --

abinit -- ABINIT Output Format [Read-only]
acesin -- ACES input format [Write-only]
acesout -- ACES output format [Read-only]
acr -- ACR format [Read-only]
adf -- ADF cartesian input format [Write-only]
adfout -- ADF output format [Read-only]
alc -- Alchemy format
arc -- Accelrys/MSI Biosym/Insight II CAR format [Read-only]
ascii -- ASCII format [Write-only]
axsf -- XCrySDen Structure Format [Read-only]
bgf -- MSI BGF format
box -- Dock 3.5 Box format
bs -- Ball and Stick format
c09out -- Crystal 09 output format [Read-only]
c3d1 -- Chem3D Cartesian 1 format
c3d2 -- Chem3D Cartesian 2 format
cac -- CAChe MolStruct format [Write-only]
caccrt -- Cacao Cartesian format
cache -- CAChe MolStruct format [Write-only]
cacint -- Cacao Internal format [Write-only]
can -- Canonical SMILES format
car -- Accelrys/MSI Biosym/Insight II CAR format [Read-only]
castep -- CASTEP format [Read-only]
ccc -- CCC format [Read-only]
cdx -- ChemDraw binary format [Read-only]
cdxml -- ChemDraw CDXML format
cht -- Chemtool format [Write-only]
cif -- Crystallographic Information File
ck -- ChemKin format
cml -- Chemical Markup Language
cmlr -- CML Reaction format
com -- Gaussian 98/03 Input [Write-only]
CONFIG -- DL-POLY CONFIG
CONTCAR -- VASP format
copy -- Copy raw text [Write-only]
crk2d -- Chemical Resource Kit diagram(2D)
crk3d -- Chemical Resource Kit 3D format
csr -- Accelrys/MSI Quanta CSR format [Write-only]
cssr -- CSD CSSR format [Write-only]
ct -- ChemDraw Connection Table format
cub -- Gaussian cube format
cube -- Gaussian cube format
dat -- Generic Output file format [Read-only]
dmol -- DMol3 coordinates format
dx -- OpenDX cube format for APBS
ent -- Protein Data Bank format
fa -- FASTA format
fasta -- FASTA format
fch -- Gaussian formatted checkpoint file format [Read-only]
fchk -- Gaussian formatted checkpoint file format [Read-only]
fck -- Gaussian formatted checkpoint file format [Read-only]
feat -- Feature format
fh -- Fenske-Hall Z-Matrix format [Write-only]
fhiaims -- FHIaims XYZ format
fix -- SMILES FIX format [Write-only]
fps -- FPS text fingerprint format (Dalke) [Write-only]
fpt -- Fingerprint format [Write-only]
fract -- Free Form Fractional format
fs -- Fastsearch format
fsa -- FASTA format
g03 -- Gaussian Output [Read-only]
g09 -- Gaussian Output [Read-only]
g92 -- Gaussian Output [Read-only]
g94 -- Gaussian Output [Read-only]
g98 -- Gaussian Output [Read-only]
gal -- Gaussian Output [Read-only]
gam -- GAMESS Output [Read-only]
gamess -- GAMESS Output [Read-only]
gamin -- GAMESS Input
gamout -- GAMESS Output [Read-only]
gau -- Gaussian 98/03 Input [Write-only]
gjc -- Gaussian 98/03 Input [Write-only]
gjf -- Gaussian 98/03 Input [Write-only]
got -- GULP format [Read-only]
gpr -- Ghemical format
gr96 -- GROMOS96 format [Write-only]
gro -- GRO format
gukin -- GAMESS-UK Input
gukout -- GAMESS-UK Output
gzmat -- Gaussian Z-Matrix Input
hin -- HyperChem HIN format
HISTORY -- DL-POLY HISTORY [Read-only]
inchi -- InChI format
inchikey -- InChIKey [Write-only]
inp -- GAMESS Input
ins -- ShelX format [Read-only]
jin -- Jaguar input format [Write-only]
jout -- Jaguar output format [Read-only]
k -- Compare molecules using InChI [Write-only]
lmpdat -- The LAMMPS data format [Write-only]
log -- Generic Output file format [Read-only]
mcdl -- MCDL format
mcif -- Macromolecular Crystallographic Info
mdl -- MDL MOL format
ml2 -- Sybyl Mol2 format
mmcif -- Macromolecular Crystallographic Info
mmd -- MacroModel format
mmod -- MacroModel format
mna -- Multilevel Neighborhoods of Atoms (MNA) [Write-only]
mol -- MDL MOL format
mol2 -- Sybyl Mol2 format
mold -- Molden format
molden -- Molden format
molf -- Molden format
molreport -- Open Babel molecule report [Write-only]
moo -- MOPAC Output format [Read-only]
mop -- MOPAC Cartesian format
mopcrt -- MOPAC Cartesian format
mopin -- MOPAC Internal
mopout -- MOPAC Output format [Read-only]
mp -- Molpro input format [Write-only]
mpc -- MOPAC Cartesian format
mpd -- MolPrint2D format [Write-only]
mpo -- Molpro output format [Read-only]
mpqc -- MPQC output format [Read-only]
mpqcin -- MPQC simplified input format [Write-only]
mrv -- Chemical Markup Language
msi -- Accelrys/MSI Cerius II MSI format [Read-only]
msms -- M.F. Sanner's MSMS input format [Write-only]
nul -- Outputs nothing [Write-only]
nw -- NWChem input format [Write-only]
nwo -- NWChem output format [Read-only]
out -- Generic Output file format [Read-only]
outmol -- DMol3 coordinates format
output -- Generic Output file format [Read-only]
pc -- PubChem format [Read-only]
pcm -- PCModel Format
pdb -- Protein Data Bank format
pdbqt -- AutoDock PDQBT format
png -- PNG 2D depiction
pos -- POS cartesian coordinates format [Read-only]
POSCAR -- VASP format
pov -- POV-Ray input format [Write-only]
  obabel -:"CC(=O)Cl acid chloride" --gen3d -O chloride.pov -xc -xf -xs -m SPF
pqr -- PQR format
pqs -- Parallel Quantum Solutions format
prep -- Amber Prep format [Read-only]
pwscf -- PWscf format [Read-only]
qcin -- Q-Chem input format [Write-only]
qcout -- Q-Chem output format [Read-only]
report -- Open Babel report format [Write-only]
res -- ShelX format [Read-only]
rsmi -- Reaction SMILES format
rxn -- MDL RXN format
sd -- MDL MOL format
sdf -- MDL MOL format
smi -- SMILES format
smiles -- SMILES format
svg -- SVG 2D depiction [Write-only]
sy2 -- Sybyl Mol2 format
t41 -- ADF TAPE41 format [Read-only]
tdd -- Thermo format
text -- Read and write raw text
therm -- Thermo format
tmol -- TurboMole Coordinate format
txt -- Title format
txyz -- Tinker XYZ format
unixyz -- UniChem XYZ format
VASP -- VASP format
vmol -- ViewMol format
xed -- XED format [Write-only]
xml -- General XML format [Read-only]
xsf -- XCrySDen Structure Format [Read-only]
xtc -- XTC format [Read-only]
xyz -- XYZ cartesian coordinates format
yob -- YASARA.org YOB format
zin -- ZINDO input format [Write-only]
```

## SMILES記法の変換

実際の変換手順はワンライナーで済む簡単なものです。例えばフェノール分子のSMILESをz-matrix形式の分子構造に変換するには，以下のコマンドを実行します。ここでは 入力形式がSMIESであることを `-i smi` オプションで指定しており，出力形式がz-matrixであることを `-o fh` オプションで指定しています。また，`--gen3D`オプションの指定により，OpenBabelは古典的な動力学法にもとづいて構造最適化を行いz-matrixの結合距離や結合角の情報を計算します。`--gen3D`オプションを指定しないと，z-matrixに必要な情報がすべて空欄となることに注意してください。本記事ではz-matrixの仕様については説明しませんが以下の出力を見れば，例えばCO結合距離は 1.364 Å と計算されたことがわかります。

```bash
$ echo 'Oc1ccccc1' | obabel -i smi -o fh --gen3D

13
O   1
C   1 1.364
C   2 1.392  1 121.397
C   3 1.396  2 119.408  1 180.3
C   4 1.396  3 120.092  2 360.0
C   5 1.395  4 120.036  3 359.9
C   2 1.391  1 117.749  3 180.1
H   1 0.971  2 108.447  3   2.2
H   3 1.085  2 121.285  1   0.4
H   4 1.087  3 120.036  2 180.0
H   5 1.087  4 120.009  3 179.9
H   6 1.087  5 120.127  4 180.0
H   7 1.086  2 120.028  1 359.8
1 molecule converted
```

また，z-matrixのみならずカーテシアン座標系へ変換することもできます。具体的には，以下のように `-o gjf` オプションを指定するとカーテシアン座標系を使ったGAUSSIANプログラムの入力ファイルが出力されます。

```bash
$ echo 'Oc1ccccc1' | obabel -i smi -o gjf --gen3D
#Put Keywords Here, check Charge and Multiplicity.



0  1
O           0.92444        -0.10001         0.10408
C           2.28120        -0.03420        -0.01566
C           2.94505         1.18747        -0.07995
C           4.33478         1.20639        -0.20936
C           5.04834         0.00790        -0.27292
C           4.37341        -1.21118        -0.20413
C           2.98497        -1.23234        -0.07408
H           0.57104         0.80441         0.09601
H           2.40334         2.12645        -0.03074
H           4.86341         2.15450        -0.26132
H           6.13034         0.02419        -0.37567
H           4.92743        -2.14489        -0.25259
H           2.45360        -2.17787        -0.02057

1 molecule converted
```

OpenBabelはあくまで簡易的な動力学法による構造最適化しかサポートしません。したがって，計算された分子構造は可視化して正しい結果を得られているか確認すべきです。OpenBabelは画像形式への変換も対応しており，簡易的な確認には便利です。画像形式で保存したいときは `-o png` とPNG形式による出力を指定し，同時に `-O out.png` と保存先のファイル名も指定します。

```bash
$ echo 'OC1ccccc1' | obabel -i smi -o png -O out.png --gen3D
```

![out.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541219/d38470c2-c6dd-dd9f-d24c-5246cd8da1c0.png)


単純に結合の有無のみを確認したい場合は， `--gen3D` オプションを削除した2次元形式での可視化が見やすく便利です。

```bash
$ echo 'OC1ccccc1' | obabel -i smi -opng -O out.png
```

![out2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541219/a70684bc-0e96-7233-dbe0-f3b23dbd86ec.png)

上記の実行例を見ていただければわかるように，OpenBabelだけでは複雑な分子構造の可視化は少々難しいものとなります。実用的な可視化を行う場合は，[avogadro](https://github.com/cryos/avogadro) のような可視化プログラムをあわせて利用するべきです。OpenBabelを使って以下のように変換したxyz形式のファイルは，avogadroへ入力することができます。

```bash
# avogadro用の入力ファイルを作成
$ echo 'Oc1ccccc1' | obabel -i smi -o xyz -O out.xyz --gen3D
```

実際に，上記コマンドで生成したファイル `out.xyz` をavogadroで開くと以下のように可視化できます。avogadroではマウス操作で視点変更できるので，OpenBabelの画像出力よりも複雑な構造もわかりやすく可視化できます。

![phoh.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/541219/a436f538-f861-fb74-6e02-19225408a293.png)

# おわりに

今回の記事でOpenBabelを使ってSMILESから簡単に分子構造を復元し，シミュレーション計算のための入力ファイルを作れるようになりました。次回以後は，作成した入力ファイルを使って第一原理計算を行い，分子の物性を求める方法を説明します。特に次回は低分子向けの第一原理計算のプログラムとして著名なGAMESSをDockerコンテナで実行可能にする方法を説明します。