## SMILES 記法
化学構造をコンピュータで表現するには、例えば ChemDraw のようなソフトウェアを用いることができます。ただ、プログラムから扱うためには、もっとシンプルなデータ構造をもつ表現のほうが適するでしょう。SMILES は化学構造を英数字のみからなる文字列リテラルの形式で表現します。
Wikipedia では以下のように説明されています。
SMILES記法SMILES記法 とは、分子の化学構造をASCII符号の英数字で文字列化した、構造の曖昧性の無い表記方法である。SMILES文字列は多くの種類の分子エディタにおいてインポート可能で、二次元の図表あるいは三次元のモデルとして表示することができる。
## 分子構造生成の実装例
以下では、ケモインフォマティクスの Python ライブラリとして著名な RDKit を用いた、SMILES 文字列をソースとする分子構造生成を実装します（コードはエッセンス部分だけ載せます）。
### SMILES の読み込み
RDKit の MolFromSmiles メソッドを使うと、SMILES を RDKit の MOL オブジェクトに変換できます。
from rdkit import Chem
smiles = "CCO" # ethanol
mol = Chem.MolFromSmiles(smiles)
### 立体構造生成
SMILES から作成した MOL オブジェクトは立体構造情報を含まないので、これを付加してやる必要があります。
RDKit では ETKDG 法 1 を使用した3次元構造生成を行うことができます。これを用いて複数の構造をサンプリングし、分子力学法 2 による構造最適化を行った際のエネルギーがもっとも安定になるものを正解構造として採用します。
from rdkit import Chem
from rdkit.Chem import AllChem
def generate3Dmol(rdmol, nconfs=100, seed=1):
# ETKDG-v3 で複数の立体配座を生成
params = AllChem.ETKDGv3()
params.randomSeed = seed
confs = AllChem.EmbedMultipleConfs(rdmol, numConfs=nconfs, params=params)
# 分子力学法 (MMFF) で最適化
for cid in confs:
AllChem.MMFFOptimizeMolecule(rdmol, confId=cid)
# 最適化を実行、エネルギーを計算
props = []
for cid in confs:
mmff = AllChem.MMFFGetMoleculeForceField(rdmol, AllChem.MMFFGetMoleculeProperties(rdmol), confId=cid)
props.append((cid, mmff.CalcEnergy()))
# 最適化後のエネルギーが低い順にソート
props.sort(key=lambda x: x[1])
# 最もエネルギーの低い構造の conformer id を取り出す
best_cid = props[0][0]
# MOL オブジェクトを生成
mol_cid = Chem.Mol(rdmol)
best_conf = rdmol.GetConformer(best_cid)
mol_cid.RemoveAllConformers()
mol_cid.AddConformer(best_conf, assignId=True)
return mol_cid
from rdkit import Chem
smiles = "CCO" # ethanol
mol = Chem.MolFromSmiles(smiles)
# 1000 構造を生成して最適な構造を取り出す
best_mol = generate3Dmol(mol, nconfs=1000)
### 3次元座標の付加
import numpy as np
def get_coords(rdmol):
coords = []
conf = rdmol.GetConformer()
for atom in rdmol.GetAtoms():
idx = atom.GetIdx()
symbol = atom.GetSymbol()
pos = conf.GetAtomPosition(idx)
coords.append([symbol, pos.x, pos.y, pos.z])
return coords
from rdkit import Chem
smiles = "CCO" # ethanol
mol = Chem.MolFromSmiles(smiles)
# 1000 構造を生成して最適な構造を取り出す
best_mol = generate3Dmol(mol, nconfs=1000)
# 座標を取得
coords = get_coords(best_mol)
### xyz ファイルとしてエクスポート
xyz は化学構造フォーマットの一つで、各原子の元素記号と3次元座標の情報のみを含むシンプルなフォーマットです。
def export_xyz(coords, mol_name):
# 各原子のシンボルと座標
rows = ["\t".join([c[0], f"{c[1]:.6f}", f"{c[2]:.6f}", f"{c[3]:.6f}"]) for c in coords]
# 原子数
natoms = f"{len(coords)}"
# ブロックを結合
content = "\n".join([natoms, mol_name] + rows)
# エクスポート
with open(f"{mol_name}.xyz", "w") as f:
f.writelines(content)
from rdkit import Chem
smiles = "CCO" # ethanol
mol = Chem.MolFromSmiles(smiles)
# 1000 構造を生成して最適な構造を取り出す
best_mol = generate3Dmol(mol, nconfs=1000)
# 座標を取得
coords = get_coords(best_mol)
# xyz ファイルをエクスポート
export_xyz(coords, "ethanol")
## CLI の実装
本記事で作成したプログラムを CLI ツールとして実装しました。PyPI にも公開してありますので、以下でインストールが可能です。
pip install smiles2mol
ソースコードは以下に配置してあります。
### 実行例
半導体材料として用いられる3,4-エチレンジオキシチオフェン (EDOT) を作成してみます。SMILES は C1COC2=CSC=C2O1 です。
smiles2mol -s "C1COC2=CSC=C2O1" -o edot.mol
以下のような構造が得られます。エチレンジオキシ基の配座がいす形で生成されていることがわかります。
EDOT
15 16 0 0 0 0 0 0 0 0999 V2000
-1.6694 0.5428 0.4018 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.5836 -0.6722 -0.5322 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.4141 -1.4708 -0.2719 O 0 0 0 0 0 0 0 0 0 0 0 0
0.6653 -0.6727 -0.0933 C 0 0 0 0 0 0 0 0 0 0 0 0
1.9520 -1.1484 -0.1064 C 0 0 0 0 0 0 0 0 0 0 0 0
3.0861 0.1228 0.1237 S 0 0 0 0 0 0 0 0 0 0 0 0
1.8400 1.2993 0.2584 C 0 0 0 0 0 0 0 0 0 0 0 0
0.6003 0.7230 0.1440 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.5499 1.4324 0.2332 O 0 0 0 0 0 0 0 0 0 0 0 0
-2.5793 1.1140 0.1920 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.6995 0.2265 1.4521 H 0 0 0 0 0 0 0 0 0 0 0 0
-2.4592 -1.3144 -0.3939 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.5548 -0.3560 -1.5826 H 0 0 0 0 0 0 0 0 0 0 0 0
2.2862 -2.1657 -0.2513 H 0 0 0 0 0 0 0 0 0 0 0 0
2.0799 2.3394 0.4263 H 0 0 0 0 0 0 0 0 0 0 0 0
1 11 1 0 0 0 0
2 12 1 0 0 0 0
2 3 1 0 0 0 0
2 1 1 0 0 0 0
3 4 1 0 0 0 0
4 8 1 0 0 0 0
5 4 2 0 0 0 0
5 6 1 0 0 0 0
6 7 1 0 0 0 0
7 15 1 0 0 0 0
8 9 1 0 0 0 0
8 7 2 0 0 0 0
9 1 1 0 0 0 0
10 1 1 0 0 0 0
13 2 1 0 0 0 0
14 5 1 0 0 0 0
M END
## おわりに
RDKit を利用した構造生成は分子力学法に基づくため、ほとんどの分子で量子化学計算を使った構造最適化のほうがよい結果を示すはずです（量子化学計算では分子力学法とは異なって、電子をあらわに考慮するからです）。
この手法は量子化学計算のための初期構造生成に用いるのがよいかと思います（初期構造作成のために、GaussView のような GUI ソフトウェアを経由する手間は省けるかもしれません）。
また、錯体や無機分子のような重原子を含む分子や、タンパク質のような巨大分子には適用しにくく、あくまで有機小分子に限ります。
#### 脚注
-
ディスタンスジオメトリ法による構造生成アルゴリズムに、多重結合などの化学的な構造制約と CCDC の経験的な二面角分布の効果を加えた方法 (S. Riniker, G. A. Landrum, J. Chem. Inf. Model. (2015).↗)。 例えば、経験的に「いす型」構造をとることが知られているシクロヘキサンなんかもうまく生成してくれます。 ↩
-
分子を振動子とみなし、原子間の相互作用や分子振動の効果を表現するポテンシャルエネルギー関数（分子力場）を計算して構造安定性を評価する手法。RDKit では Merck 社が開発した MMFF という分子力場を使用できます (T. A. Halgren, J. Comp. Chem. (1996)↗)。 ↩