SMILESからGaussianのインプットファイルを作成する
はじめに
材料物性のデータベースはMaterials Project、Polyinfoなど最近でこそ色々公開されてきていますが、実用面ではまだまだ十分なわけではなく、自分自身でシミュレーション等を行って物性を取得する必要があるのが現状です。材料系のシミュレーション技術は対象とする時空間スケールに応じた手法が色々ありますが、ほとんどのシミュレータは化学構造や計算条件を記載したテキストファイルを元にシミュレーションを実行するというスタイルがとられます。
大量データをシミュレーションで生成しようとすると・・・
通常、インプットファイルはGUI上で原子を配置したり結合を作ったりして作成します。こういったGUIは、日本だとクロスアビリティから発売されているWinmostarがメジャーでしょうか。
さて、例えば1000件の有機物分子の物性を何らか取得したいとします。そうなるとインプットファイルの作成もGUIで原子をポチポチ触りながら1000種類の分子を作ることになってしまいます。・・・大変すぎますね。ということで、今回はRDKitを使ってGaussianのインプットファイルを作成してみます。ここで紹介する方法は、分子の座標以外の部分を書き換えればNWChemやMOPACなどのGaussian以外の量子化学計算ソルバーに対しても、そのまま利用可能です。
1. SMILESから3次元座標を生成
下記のコードでSMILESから3次元構造が埋め込まれRDKit Molオブジェクトを作成できます。大まかな流れとしては、
1. SMILESからRDKitのMolオブジェクトを作成
2. 水素を付加
3. 3次元座標の埋め込み
4. 分子力場で構造最適化
というものです。
from rdkit import Chem
from rdkit.Chem import AllChem
smiles = "N(=N/c1ccccc1)/c2ccccc2"
mol = Chem.MolFromSmiles(smiles)
mol = Chem.AddHs(mol)
AllChem.EmbedMolecule(mol)
AllChem.MMFFOptimizeMolecule(mol)
Gaussianのインプットファイルは計算条件を記載したセクションと、分子の情報が記載されたセクションに分かれています。まずは後者の分子の情報が記載されたセクションを作成してみます。このセクションは、系内の原子について
- 元素シンボル
- 座標
をそれぞれ指定する必要があるので、上記のコードから元素シンボルと座標を取り出す必要があります。各原子の元素シンボルは次のように取り出せます。
for atom in mol.GetAtoms():
print(atom.GetSymbol())
各原子の座標は次のように取り出せます。
for c in mol.GetConformers():
print(c.GetPositions())
これらを組み合わせることで次のように元素シンボルと座標をpandas dataframeで取得することができます。
import pandas as pd
atoms = pd.DataFrame(columns=['element', 'x', 'y', 'z'])
elements = []
for atom in mol.GetAtoms():
elements.append(atom.GetSymbol())
for c in mol.GetConformers():
positions = c.GetPositions()
atoms['element'] = elements
atoms['x'] = positions[:, 0]
atoms['y'] = positions[:, 1]
atoms['z'] = positions[:, 2]
atomsの中身はこんな感じです(長いので先頭十行だけ表示)。これで分子の座標のセクションは完成です！
element x y z
0 N -0.720965 -0.809274 -2.190419
1 N 0.527900 -0.968200 -2.180961
2 C 1.289714 -0.602906 -1.054229
3 C 1.856365 0.673989 -0.964581
4 C 2.647289 1.024660 0.131455
5 C 2.893721 0.091931 1.138021
6 C 2.358346 -1.192089 1.045468
7 C 1.567841 -1.539652 -0.051860
8 C -1.383721 -0.262658 -1.074494
9 C -1.902405 -1.097945 -0.078075
10 C -2.596961 -0.561362 1.007968
2. 計算条件のセクションを作る
一番シンプルなやり方は下記のようにヘッダーの部分をそのまま文字列として入力してしまう方法です。計算対象によってファイル名を変えたい場合や電荷やスピン多重度を変えたい場合は変数として代入してもいいかもしれないですね。
header = """
!%chk=temp
#p b3lyp/6-31g* scf=tight
azobenzene
0 1
"""
print(header)
3. Gaussianのインプットファイルとして書き出す
これで準備は整ったのであとは書き出すだけです。次のように出力できます
with open('azobenzene.gjf', 'w') as gin:
gin.write(header)
for v in atoms.values:
gin.write(' '.join(map(str, v)) + '\n')
生成したazobenzen.gjfの中身は次のようになっています（これも先頭のみ）。これで完成です！
!%chk=temp
#p b3lyp/6-31g* scf=tight
azobenzene
0 1
N -0.7209650995770261 -0.8092744269939769 -2.1904190657860143
N 0.5278996647584542 -0.9682002770105621 -2.180961084549009
C 1.2897141129687846 -0.6029063245123301 -1.0542287084888842
C 1.8563645290676343 0.6739893731151152 -0.96458126992599
C 2.6472891939674357 1.0246598525781698 0.13145485439940768
C 2.8937212695795465 0.09193147487459612 1.1380211103823847
C 2.3583458549200627 -1.192089269166475 1.0454680102159424
C 1.5678408543620515 -1.5396521646615697 -0.051859684093346965
C -1.3837208753186618 -0.26265811671665446 -1.074493625713662
C -1.9024045685996003 -1.0979451348037084 -0.07807508001816853
おわりに
上記では最初にsmilesを指定していましたが、予め準備しておいたSMILESのリストに対してループさせたりPubChemやChemSpiderのAPIと繋げたりすると、計算資源の有る限り自動でシミュレーションを行うことができてとても楽ちんです。
おまけ：py3Dmolでの可視化
一応構造がちゃんとできているか確認しましょう。こーいうときはpy3Dmolで可視化できて・・・
import py3Dmol
v = py3Dmol.view(width=400, height=400)
mb = Chem.MolToMolBlock(mol)
v.addModel(mb, 'sdf')
v.setBackgroundColor('0xeeeeee')
v.setStyle({'stick':{}})
v.zoomTo()
v.show()
ちゃんとアゾベンゼンになってますね。ちなみにjupyter notebook上ではグリグリ動かすこともできます。
Discussion