計算手法とエネルギー・最適化構造の関係:コンフォメーション探索における注意点

01_計算化学

RDKitを用いたコンフォメーション探索:MMFFによる配座異性体生成とクラスタリング」という記事では,

  1. RDKitに実装されているディスタンス・ジオメトリー法による配座発生
  2. MMFFによる構造最適化と同一構造の除去
  3. 取得した多数のコンフォマーをDBSCANを用いてクラスタリング

という一連の作業を行い,コンフォマーを取得しました.

この過程では一貫してMMFFを用いてエネルギー評価を行ってきました.今回は用いる計算レベルによってエネルギーの値,ひいては最適化構造にどのような影響を与えるかを見ていきたいと思います.

本記事では量子化学計算エンジンとしてPsi4を用いています。本ブログのPsi4関係の記事は「Psi4の使い方」という記事にまとめています。参照してみてください。
【まとめ】Psi4の使い方
本ブログではPythonのオープンソースライブラリであるPsi4を用いて,量子化学計算の初心者を対象にした記事を多数用意しています。 記事の内容を一つずつ習得することで,これから計算化学を始めたい人が知識・経験ゼロの状態から徐々にステップアップしていくことで,量子化学計算の基本コ...

計算手法とポテンシャルエネルギー曲面の形

計算化学に馴染みのない人にありがちなのが「計算した構造だから正しい」と思い込んでしまうことです.

実際にはほぼ全ての分子についてシュレディンガー方程式は厳密には解けませんので,計算にあたって何らかの近似を行っています.どのような近似を行うかによって計算手法それぞれに特徴があり,ある手法では考慮できない相互作用なども存在します.

特に近年の量子化学計算における主流である密度汎関数法(DFT)では,分子軌道法とは異なり精度を上げるための理論的な道筋が明確でないこともあり,非常に多くの汎関数が開発されています.その様相は「Density Functional Theory Zoo」などと揶揄される程です.

一昔前は,「有機低分子はとりあえずB3LYPで最適化すれば大丈夫」といった風潮もありましたが,近年ではB3LYPでうまく評価できない系に関する知見も蓄積されてきたため,そのような「幻想」は消滅したように思います.一方で決定的なものは存在しないため,汎関数の選択は困難なままです.

構造最適化ができない例

まずはやや極端な例ですが,アルゴン二量体を取り上げてみます.このような希ガス二量体は分散力のみで安定化されているため,分散力を考慮できない計算手法では最適化構造を見出すことができません.

ポテンシャルエネルギー曲面の図示

アルゴン–アルゴン間の距離を0.9から5.4Åまで変化させながらエネルギー計算を行っていくことでポテンシャルエネルギー曲面を描いてみましょう.

  • HF
  • MP2
  • B3LYP
  • B3LYP-D3(BJ)
  • wB97XD

の5つについて,cc-pVTZを基底関数系として計算することにします.

import psi4
print(psi4.__version__) # 1.3.2

# CPUとRAMは各自の環境に合わせて設定してください
psi4.set_num_threads(10)
psi4.set_memory('14GB')

ar_geom = '''
0 1
Ar  0 0 0
Ar  0 0 {}
'''

def calc_ar2_energy(level):
    theory = level + '/cc-pvtz'
    calc_energy = []
    for i in np.arange(0.9, 5.5, 0.1):
        ar2 = psi4.geometry(ar_geom.format(i))
        calc_energy.append(psi4.energy(theory, molecule=ar2))
    return calc_energy


hf = calc_ar2_energy('hf')
mp2 = calc_ar2_energy('mp2')
b3lyp = calc_ar2_energy('b3lyp')
b3lyp_d3bj = calc_ar2_energy('b3lyp-d3bj')
wb97x_d = calc_ar2_energy('wb97x-d')

得られたエネルギーをグラフ化したものが次のものになります.このままではエネルギーの極小値があるのかよくわかりません.

Ar dimer energy

エネルギー極小値を求める

scipy.optimize.fminbound(f, x1, x2)

scipyには色々な最適化手法が実装されています.ここでは区間x1からx2の間の最小値を求めるoptimize.fminboundを使うことで,先ほどのエネルギー極小値を求めてみます.

下のコードではwB97X-Dの場合の極小値を求めています.

from scipy import optimize

def wb97xd_energy(r):
    ar2 = psi4.geometry(ar_geom.format(r))
    return psi4.energy('wb97x-d/cc-pvtz', molecule=ar2)

optimize.fminbound(wb97xd_energy, x1=1.0, x2=5.0)
# 4.191968507707002

先ほどのグラフを実際に拡大してみると,確かに4.2Å付近で極小値を取っていることが確認できます.

Min wb97xd

同様にしてMP2とB3LYPで極小値を求めてみると,

  • MP2: 3.93358892716643,
  • B3LYP: 4.99999521993633

となりました.B3LYPは分散力を考慮できないため,区間[1.0,5.0]でエネルギーは単調減少となります.そのため4Å付近の極小値が存在せず構造最適化ができません.

コンフォマー間のエネルギーの大小が異なる例

続いて構造最適化は可能なものの,コンフォマー間のエネルギーの大小が異なる例を見てみます.「A sobering assessment of small‐molecule force field methods for low energy conformer predictions」という論文でもとりあげられているように,こちらのケースは非常に多く存在します.

1,2-ジフルオロエタンを例に,anti型とgauche型のエネルギーをいくつかの計算レベルで比べてみることとします.この分子は実験・理論の両面からよく研究されており,gauch型配座の方が安定であることがわかっています.

Psi4を用いた制限付き構造最適化

psi4.set_options({‘optking__fixed_dihedral’: ‘x1 x2 x3 x4 angle’})

まずは「RDKitを用いて制約付きで立体構造を生成する」で行ったような二面角をスキャンしながらの構造最適化を行ってみます.psi4の最適化モジュール「optking」のオプションで「fixed_dihedral」を用います.

注意点としては,このオプションでは指定する原子番号のカウントが1から始めます.pythonのように0からでない点に注意が必要です.

下記のコードではRDKitのコンフォマーオブジェクトの座標からPsi4のMoleculeオブジェクトを作成するためにgenXyzAngという関数を定義しています.

その上で,

  1. RDKitを用いて分子を構築し,19のコンフォマーを発生
  2. 各コンフォマーに対して,0から180°までの二面角を対応させる
  3. RDKitのコンフォマーオブジェクトの座標を用いてPsi4のMoleculeオブジェクトを作成
  4. Psi4で二面角を固定しながらHF/STO-3Gで構造最適化

という処理を行っています.

from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Draw, PandasTools, Descriptors
from rdkit.Chem.Draw import IPythonConsole
print(rdBase.rdkitVersion)
# 2019.09.3

mol = Chem.MolFromSmiles('FCCF')
mol = Chem.AddHs(mol)
confids = AllChem.EmbedMultipleConfs(mol, 19, AllChem.ETKDGv2())
len(confids) # 19

def genXyzAng(m, angle, confId=-1):
    conf = m.GetConformer(confId)
    Chem.rdMolTransforms.SetDihedralDeg(conf, 0,1,2,3, angle)

    xyz = '0 1'
    for atom, (x,y,z) in zip(m.GetAtoms(), conf.GetPositions()):
        xyz += '\n'
        xyz += '{}\t{}\t{}\t{}'.format(atom.GetSymbol(), x, y, z)

    return xyz


psi4.set_output_file('FCCF_HF.txt')

HF_FCCF = []
for i, angle in zip(confids, range(0,181,10)):
    conf = mol.GetConformer(i)
    Chem.rdMolTransforms.SetDihedralDeg(conf, 0, 1, 2, 3, angle)
    pmol = psi4.geometry(genXyzAng(mol, angle, i))
    # psi4による二面角の指定:1,2,3,4であることに注意
    psi4.set_options({'optking__fixed_dihedral': '1 2 3 4 {}'.format(angle)})
    psi4.set_options({'geom_maxiter': 200})
        
    e = psi4.optimize('hf/sto-3g', molecule=pmol)
    HF_FCCF.append(e)

得られたエネルギーをプロットしたものが下のグラフです.これから,HF/STO-3Gレベルにおいては

  • anti型とgauche型2つのコンフォマーがあること
  • この計算レベルではgauche型の方が安定であること
  • gauche型は70°あたりが最適構造であること

などが見てとれます.

Scan HF STO 3G

Psi4を用いて様々な計算レベルで構造最適化

続いて二面角60°と180°を初期構造として様々な計算手法で構造最適化を行い,anti型とgauche型のエネルギー差を比べてみます.

具体的には

  • HF, MP2, B3LYP, B3LYP-D3(BJ)
  • STO-3G, 6-31G(d), 6-31+G(d,p), 6-311+G(2d,p)

の組合せで4 x 4通りの計算レベルについて調べてみます.

また

  • GetSubstructMatchを利用して二面角の指定範囲をmatchとして取得,設定
  • 最適化したPsi4のMoleculeオブジェクトの座標をRDKitのコンフォマーへと移すpsi2conf関数を定義

することで,一般性の向上と最適化後の二面角の取得を行っています.最適化後の座標を含んだMOLファイルも保存しておきます.

なおpsi4.Moleculeオブジェクトから(x,y,z)座標を取り出す際,単位がBohrになっていますのでpsi4.constants.bohr2angstromsを用いて変換しています.

from datetime import datetime
import itertools

smiles = 'FCCF'
pattern = Chem.MolFromSmarts('FCCF')
theory = ['hf', 'mp2', 'b3lyp', 'b3lyp-d3bj']
basis = ['sto-3g', '6-31G(d)', '6-31+G(d,p)', '6-311+G(2d,p)']
angles = [60, 180]


def mol2conf(sm, angle):
    # sm: smiles, angle: dihedral angle to be set in the conformer
    mol = Chem.MolFromSmiles(sm)
    mol = Chem.AddHs(mol)
    _ = AllChem.EmbedMolecule(mol)
    match = mol.GetSubstructMatch(pattern)
    conf = mol.GetConformer(-1)
    Chem.rdMolTransforms.SetDihedralDeg(conf, *match, angle)

    return mol, conf


def conf2xyz(mol, conf):
    # mol: rdkit mol, conf: conformer of mol
    xyz = '0 1'
    for atom, (x, y, z) in zip(mol.GetAtoms(), conf.GetPositions()):
        xyz += '\n'
        xyz += '{}\t{}\t{}\t{}'.format(atom.GetSymbol(), x, y, z)

    return xyz


def psi4_optimize(th, ba, xyz):
    # th: theory, ba: basis set, xyz: coordinate
    psi4.set_module_options('optking', {'geom_maxiter': 200})
    psimol = psi4.geometry(xyz)
    energy = psi4.optimize(th + '/' + ba, molecule=psimol)

    return energy, psimol


def psi2conf(psimol, conf):
    # mol: psi4 mol, conf: conformer of rdkit mol
    for i in range(psimol.natom()):
        x = psimol.x(i) * psi4.constants.bohr2angstroms
        y = psimol.y(i) * psi4.constants.bohr2angstroms
        z = psimol.z(i) * psi4.constants.bohr2angstroms
        conf.SetAtomPosition(i, (x, y, z))


results = []
for (angle, th, ba) in itertools.product(angles, theory, basis):
    print('--')
    print('{}/{} with {}deg'.format(th, ba, angle))
    mol, conf = mol2conf(smiles, angle)
    geom = conf2xyz(mol, conf)
    match = mol.GetSubstructMatch(pattern)

    now = datetime.now()
    log_file = '{}_{}_{}.txt'.format(datetime.strftime(now, '%y%m%d_%H%M'), th, ba[:5])
    psi4.set_output_file(log_file)

    energy, psimol = psi4_optimize(th, ba, geom)
    psi2conf(psimol, conf)
    final_angle = Chem.rdMolTransforms.GetDihedralDeg(conf, *match)

    mol = Chem.RemoveHs(mol)
    Chem.MolToMolFile(mol, '{}_{}_{}_{}deg.txt'.format(th, ba[:5], smiles, angle))

    results.append((smiles, angle, th, ba, energy, final_angle))

with open('./result.txt', 'w') as f:
    print(results, file=f)

計算レベル毎のgauch型とanti型のエネルギー差は以下のようになりました.なお最適化構造の二面角は,anti型では179.5から180.0°とほとんど変化がなかったのに対し,gauche型は69.0から72.4°と計算手法による違いが大きかったです.

HF/STO-3Gではgauche型の方が安定になっていますが,より大きい基底関数系を用いるとanti型の方が安定になっています.この結果は他の計算レベルでの結果を踏まえると,

  • HFを用いるとうまく評価できないエラー
  • 最小基底系では評価できないエラー

がキャンセルしあって偶然出てきた結果であると結論付けられそうです.

計算手法 エネルギー差 (kcal/mol)
HF/STO-3G -0.3290318484135456
HF/6-31G(d) 0.45126090174887473
HF/6-31+G(d,p) 0.1304722850608613
HF/6-311+G(2d,p) -0.30009113129279874
MP2/STO-3G -0.3584640475723512
MP2/6-31G(d) -0.21332636897446614
MP2/6-31+G(d,p) -0.4833678449151578
MP2/6-311+G(2d,p) -0.8516391230649175
B3LYP/STO-3G -0.3677167585435006
B3LYP/6-31G(d) -0.3728175376350604
B3LYP/6-31+G(d,p) -0.6567267035116128
B3LYP/6-311+G(2d,p) -0.8938457404183217
B3LYP-D3(BJ)/STO-3G -0.373095195655657
B3LYP-D3(BJ)/6-31G(d) -0.3839035364556454
B3LYP-D3(BJ)/6-31+G(d,p) -0.6666511601584338
B3LYP-D3(BJ)/6-311+G(2d,p) -0.9046106549686137

エネルギー計算と構造最適化を異なる計算レベルで行う

一般に構造最適化はエネルギー計算を複数回(通常は数十回以上)行うため,エネルギー計算よりもはるかに計算量が多くなります.そのため,構造最適化を高レベルで行うことが困難なことが多々あります.

そのような場合に,

  1. 構造最適化(と振動数計算)を低レベルの計算手法で行い
  2. 最適化構造に対して高レベルのエネルギー計算を行う

という2段階でエネルギーを評価することが行われます.

例えば下の囲みで示す

計算レベル1//計算レベル2

のように「//」で区切られた表記方法は,

  1. 計算レベル2で最適化・振動数計算
  2. 計算レベル2で最適化した構造に対して,計算レベル1でエネルギー計算

という2段階でエネルギーを評価する計算手法を意味しています.ある程度の大きさの有機化合物を対象とした計算化学の論文でよく見る表記方法です.

この場合,自由エネルギーなどの算出には,計算レベル2の振動数計算で得られた補正項を計算レベル1で得たエネルギーに加えることで対応しています.

もちろん計算レベル1と計算レベル2では最適化構造は異なるため,このような補正は厳密には正しくありませんが,振動数は計算手法を変えてもあまり変化しない(系統的に少しズレるだけ)ことなどを理由に広く用いられているようです.

なおIR振動数を予想するために掛ける係数のデータベースなども古くから存在します.

終わりに

今回は「計算手法とエネルギー・最適化構造の関係:コンフォメーション探索における注意点」という話題について,

  • 計算方法によっては構造最適化できない場合がある
  • 計算手法によっては異性体間のエネルギーの大小が逆転する場合がある

ということを具体例をもとに見てきました.またpythonから利用可能な量子化学計算パッケージであるpsi4の利用方法の復習も行いました.

今回見たように,計算の実行にあたっては興味のある系をきちんと理解した上で,適切な計算手法を選択することが大切になります.

もっとも量子化学計算は計算量が多いため,いたずらに高レベルの手法を用いると計算終了までに時間がかかり,スループットが悪くなってしまいます.低コストで高精度の計算手法を目指して,

などの開発が活発に行われているのも納得でしょう.

次回は,シュレディンガー方程式の解である波動関数から得られる電子密度などの各種値について見ていきたいと思います.

>>次の記事:「計算化学における電荷:Psi4を用いた電子密度解析

コメント

タイトルとURLをコピーしました