Python で量子化学計算 - ご注文はリード化合物ですか？〜医薬化学録にわ〜

明けましておめでとうございます。
今年も、創薬化学、ケモインフォマティクスを中心にブログを更新していこうと思います。

量子化学計算といえば、Gaussian や GAMESS を利用している方も多いと思いますが、Python には量子化学計算をしてくれるライブラリがあります。
psicode.org

conda install psi4 psi4-rt python=3.7 -c psi4

で使えるようになります。

使い方ですが、基本的な使い方はこちらのブログにまとまっています（この記事もこちらを参考にして書いています）。
future-chem.com

また、psi4 を使いやすくした psikit なるものもあります（この記事もこちらを参考にして書いています）。
github.com

これだけではこの記事の新規性がないので、いくつかの分子をまとめて量子化学計算できるコードを作成したので、これを元に psi4 の基本的な使い方をまとめます。

import numpy as np
import pandas as pd
from rdkit import Chem
from rdkit.Chem.rdDistGeom import ETKDGv3, EmbedMolecule
from rdkit.Chem.rdForceFieldHelpers import MMFFHasAllMoleculeParams, MMFFOptimizeMolecule
from rdkit.Chem.Draw import IPythonConsole
import psi4

import datetime
import time

# ハードウェア側の設定
psi4.set_num_threads(nthread = 2)
psi4.set_memory("1GB")

# 入力する分子
smiles_list = ["C", "CC", "CCC", "CCCC", "CCCCC"]

for smiles in smiles_list:
    t = datetime.datetime.fromtimestamp(time.time())
    psi4.set_output_file("../result/psi4_log/{}{}{}_{}{}_{}.log".format(t.year,
                                                  t.month,
                                                  t.day,
                                                  t.hour,
                                                  t.minute, smiles))

    # SMILES から三次元構造を発生させる
    mol = Chem.MolFromSmiles(smiles)
    mol = Chem.AddHs(mol)
    params = ETKDGv3()
    params.randomSeed = 1
    EmbedMolecule(mol, params)

    # MMFF で構造最適化
    MMFFOptimizeMolecule(mol)
    conf = mol.GetConformer()

    # Psi4 に入力可能な形式に変換
    mol_input = "0 1"

    for atom in mol.GetAtoms():
        mol_input += "\n " + atom.GetSymbol() + " " + str(conf.GetAtomPosition(atom.GetIdx()).x)\
        + " " + str(conf.GetAtomPosition(atom.GetIdx()).y)\
        + " " + str(conf.GetAtomPosition(atom.GetIdx()).z)

    geometry = psi4.geometry(mol_input)

    # 汎関数、基底関数を設定
    level = "b3lyp/6-31G*"

    # 構造最適化計算(単位： au)
    energy, wave_function = psi4.optimize(level, molecule = geometry, return_wfn = True)

    LUMO_idx = wave_function.nalpha()
    HOMO_idx = LUMO_idx - 1

    # HOMO を表示（単位： au）
    homo = wave_function.epsilon_a_subset("AO", "ALL").np[HOMO_idx]
    lumo = wave_function.epsilon_a_subset("AO", "ALL").np[LUMO_idx]

    with open("../result/psi4_summary.txt", "a") as writer:
        output_string = "{}\nEnergy: {}\nHOMO: {}\nLUMO: {}\n\n".format(smiles, str(energy), str(homo), str(lumo))
        writer.write(output_string)

最初に、ハードウェアの設定をします。スレッドやメモリを大きくすれば、計算速度も上がると思いますが、手元のマシンのスペックと相談することになります。
出力するログファイルも設定できますが、当然、分子ごとに異なるので、for 文の中で名前を決めてあげます。

量子化学計算では当然三次元構造が必要なので、SMILES から立体配座を作成します。
水素原子を付加した後、ETKDG というアルゴリズムで構造を発生させます。イメージ的には、distance geometry 法に近いと思いますが、アルゴリズムを色々改良しているようです。
コンフォメーションを発生させてから、MMFF94 で構造最適化をしていますが、やった方がいいのか、厳密な検討はしていません（通常の DG 法は最適化する）。
余談ですが、MMFF94S もありますが、結晶構造の場合、こちらを使った方がいいらしいです（厳密な比較はしていない）。

立体配座を作成したら、いよいよ、psi4 にインプットするデータを作成します。データは、Gaussian や GAMESS へのインプットファイルをイメージすると分かりやすいです。最初の行に形式電荷などの情報、その下に元素記号と座標を書いていきます。String として作成したら、psi4.geometry で psi の Molecule オブジェクトに変換しましょう。今回は XYZ 形式（座標をそのまま使う）でインプットファイルを作成していますが、Z-matrix 形式にした方が計算が収束しやすい気がします（case by case ですが）。

計算の基本設定は簡単で、汎関数/基底関数を決めて、psi4.optimize に引数として入力すると計算をしてくれます。psi4.optimize はデフォルトだとエネルギーしか返してくれないので、return_wfn を True にして、色々な情報が入っている wave_function オブジェクトを取り出しましょう。
一点計算を行いたい時は、optimize を energy に変更すれば同じようにできます。

wave_function.nalpha でHOMO の軌道の番号が取れますが、1 始まりなので、wave_function.epsilon_a_subset('AO', 'ALL') で軌道のエネルギーを取り出すときは、番号が一つずれる（0 始まり）なので、気をつけてください。
epsilon_a_subset の引数を変えることで、他の情報も取り出せそうですが、psi4 のドキュメントが思ったより充実していないので、試行錯誤する必要がありそうです・・・

計算結果はログファイルにまとめられていますが、分量が多く、見にくいので、エネルギーを txt ファイルに出力しておくと、結果を確認しやすいです。

データの単位ですが、RDKit で構造を発生させた場合、座標の単位はオングストローム(10^-10 m) です。
psi4 へのインプットファイルも同様です。
今回は出力していませんが、psi4 の出力座標の単位はボーアになっているので、長さのスケールには注意してください。

また、エネルギーの単位は原子単位（ハートリー）になっています。
単位はちゃんとドキュメントに書いていて欲しい・・・

psi4 もまだまだ理解していないことも多いので、随時勉強したことを紹介していければと思います。
とりあえず原子レベルでの情報を取り出せるようになりたい。