結晶構造解析のcifデータを元に量子化学計算のPsi4の入力ファイルを生成してみたよ

今回の内容は、次の記事の計算に用いたインプットファイルの作成法についてです。

結晶構造解析データを元に量子化学計算を行いたい

量子化学では分子のエネルギーを計算的に求めることができます。計算で求めたエネルギーを使って、分子の最安定構造や、電子軌道の形、分子の振動、分子間の相互作用といった種々の情報を得ることができます。このように量子化学計算は、化学においてとても有用で、なくてはならないものですが、あくまでも理論計算なので、量子化学計算で求めた結果がどれだけ実験的な事実と一致するかを確認する必要があります。化学という学問において、実験と計算は車の両輪のような関係で、どちらか片方ではうまく機能しません。

量子化学計算の入力ファイルとして実験から得られた分子構造を用いるということは一般によく用いられる手法です。今回は、実験的に得られた結晶構造解析データであるcifファイルから二つの分子の構造を抽出し、量子化学計算プログラムであるPsi4で分子間相互作用計算を行うdatファイルを作るところまでをまとめます。Psi4での計算の仕方や結果の見方は上で紹介した記事をご参考ください。

手順

(1)cifファイルを準備します
(2)結晶構造描画ソフトMercuryでcifファイルを開きます
(3)分子間相互作用計算を行いたい二つの分子をそれぞれ表示しxyz形式で保存
(4)Pythonで二つのxyz形式のファイルをdatファイルに変換

上で紹介した記事では、薬効物質であるアスピリンの二量体間の分子間相互作用を計算するために、Crystallography Open Database(COD)¹という結晶構造のデータベースからアスピリンの結晶構造を無料でダウンロードしています。
(2)においてMercuryというソフトを用いていますが、Cambridge Crystallographic Data Center²が配布している無料版のものでできます。
計算したい分子のxyz形式のファイルさえできて仕舞えば、Psi4の計算に用いるdatファイルはコピペして作っても難しくないですが、これをPythonプログラムで変換します。

xyz形式のファイルをdatファイルに変換するプログラム

手順の(3)で作った二つのxyz形式のファイルのパスをそれぞれ変数file1とfile2に代入しています。実際に用いる場合はここを適切なパスに変更します。

#xyz形式のファイルを指定
file1 = "xyz形式のファイルのパス1"
file2 = "xyz形式のファイルのパス2"

files = [file1, file2]
mols = []

for i in range(len(files)):
    with open(files[i], "r") as f:
        lines = f.readlines()
        
    for j in range(len(lines)):
        lines[j] = "     " + lines[j]

    mol = ''.join(lines[2:])
    
    mols.append(mol)

psi4_text = "memory 600 mb\n"\
            "\n"\
            "molecule_dimer {\n"\
            "     0 1\n"\
            + mols[0] +\
            "     --\n"\
            "     0 1\n"\
            + mols[1] +\
            "     units angstrom\n"\
            "}\n"\
            "\n"\
            "set globals {\n"\
            "    basis          3-21g\n"\
            "}\n"\
            "\n"\
            "set sapt {\n"\
            "    print          1\n"\
            "}\n"\
            "\n"\
            "energy('sapt0')\n"
        
with open("input.dat", 'wb') as a_file:
    a_file.write(psi4_text.encode('ASCII'))

Mercuryでcifファイルをxyz形式にするところが自動化されてないので、わざわざpythonプログラムを使わず、手作業で書き換えてもそんなに作業量は変わらないかもしれませんね。

1. Crystallography Open Database
   http://www.crystallography.net/cod/ ↩

2. Cambridge Crystallographic Data Center
   https://www.ccdc.cam.ac.uk/ ↩