構造最適化の計算レベルと反応エネルギーの関係

2020年8月16日公開

はじめに

DFT法をはじめとした量子化学計算によって化学反応の熱力学量（エネルギー）を求めることが日常的に行われるようになってきています．これらの値を求めるためには分子構造が必要ですが，通常，構造もDFT法などで構造最適化計算を行って得ます．

構造最適化計算は基底関数を小さめにして低コストな方法で行い，高精度なエネルギーを得るために基底関数を大きくして全電子エネルギーを求め直すことがしばしば行われています．しかし，コンピュータが高速化し，さらに様々な基底関数や汎関数（あるいは電子相関理論とその近似法）が提案されているなかで，どの程度までコストを抑えた手法で構造最適化を行えばよいのか，悩みの多いところです．

本稿では，様々な計算手法で構造最適化計算を行い，その構造を用いてDLPNO-CEPA/1/(ma-)def2-SVPで反応のΔEを求めて比較することで，構造最適化の計算レベルの影響を調べました．

結果と考察

有機分子の反応

まずはスキーム1に示すD-フルクトースからヒドロキシメチルフルフラール（HMF）と水3分子が生成する反応について計算を行いました．

スキーム1．D-フルクトースからのHMF生成．

計算結果を表1にまとめます．SCS-MP2/def2-TZVPで構造最適化計算を行った結果を基準としたΔEの差分（ΔΔE）や同じくSCS-MP2/def2-TZVPの構造を基準とした構造のRMSD値も記載してあります．ちなみに構造レベルは構造最適化の計算レベルのことで，ΔE，ΔΔEの値は構造最適化に用いた手法ではなく，DLPNO-CEPA/1/def2-SVPのシングルポイント計算で求めたものです．

表1．D-フルクトースからのHMF生成のΔE，ΔΔEと各分子の構造のRMSD（SCS-MP2/def2-TZVP基準）．
構造レベル ΔE [kcal/mol] ΔΔE [kcal/mol] RMSD (D-フルクトース) [Å] RMSD (水) [Å] RMSD (HMF) [Å] 平均RMSD [Å]
SCS-MP2/def2-TZVP 39.78 - - - - -
SCS-MP2/def2-SVP 40.09 -0.32 0.043 0.010 0.014 0.022
PBE-D4/def2-SVP 39.30 0.47 0.112 0.012 0.031 0.052
TPSS-D3(BJ)/def2-SVP 39.22 0.55 0.111 0.011 0.027 0.049
B3LYP-D4/det2-SVP 40.63 -0.85 0.084 0.007 0.020 0.037
PBE0-D3(BJ)/def2-SVP 40.06 -0.28 0.089 0.008 0.024 0.040
PBE0-D4/def2-SVP 40.08 -0.30 0.089 0.008 0.025 0.040
PBE-D4/def2-TZVP 39.84 -0.06 0.057 0.005 0.035 0.032
B3LYP-D3(BJ)/det2-TZVP 40.32 -0.55 0.034 0.004 0.022 0.020
B3LYP-D4/det2-TZVP 40.32 -0.54 0.035 0.004 0.021 0.020
PBE0-D3(BJ)/def2-TZVP 40.85 -1.07 0.038 0.003 0.034 0.025
PBE0-D4/def2-TZVP 40.87 -1.10 0.038 0.003 0.033 0.024
PBEh-3c 39.98 -0.21 0.042 0.003 0.035 0.027
B97-3c 40.35 -0.58 0.035 0.004 0.075 0.038

この計算データから求めたΔEのmax-minは1.65 kcal/molで標準偏差は0.46でした．炭素の混成状態が変わったり，分子内に水素結合が多数あったりして構造最適化の計算レベルの影響がありそうな反応ですが，意外と汎関数・基底関数の依存性はないようです．強いて言えばPBE0-D3(BJ)またはPBE0-D4とdef2-TZVPを用いた場合にSCS-MP2/def2-TZVPとの差が1.1 kcal/mol程度あります．この2つの手法については，むしろdef2-SVPを使ったほうがエネルギーの差異は小さくなっています．ただ構造は平均的にはdef2-TZVPのほうが良いようで，構造とエネルギーの関係は見えてきません．

有機分子の異性化

続いてスキーム2に示すISOL24の反応#10について検討しました．計算結果は表2にまとめます．

スキーム2．ISOL24の反応#10．
表2. ISOL24の反応#10のΔE，ΔΔEと各分子の構造のRMSD（SCS-MP2/def2-TZVP基準）．
構造レベル ΔE [kcal/mol] ΔΔE [kcal/mol] RMSD (化合物3) [Å] RMSD (化合物4) [Å] 平均RMSD [Å]
SCS-MP2/def2-TZVP 5.24 - - - -
SCS-MP2/def2-SVP 5.22 0.024 0.050 0.022 0.036
PBE-D4/def2-SVP 4.92 0.32 0.068 0.094 0.081
TPSS-D3(BJ)/def2-SVP 5.11 0.13 0.049 0.093 0.071
B3LYP-D4/det2-SVP 5.22 0.017 0.047 0.082 0.065
PBE0-D3(BJ)/def2-SVP 5.09 0.15 0.054 0.082 0.068
PBE0-D4/def2-SVP 5.07 0.17 0.057 0.080 0.069
PBE-D4/def2-TZVP 5.01 0.23 0.066 0.069 0.067
B3LYP-D3(BJ)/det2-TZVP 5.24 0.00 0.050 0.070 0.060
B3LYP-D4/det2-TZVP 5.21 0.033 0.054 0.070 0.062
PBE0-D3(BJ)/def2-TZVP 5.15 0.090 0.064 0.066 0.065
PBE0-D4/def2-TZVP 5.13 0.11 0.067 0.065 0.066
PBEh-3c 5.36 -0.12 0.036 0.110 0.073
B97-3c 5.28 -0.040 0.033 0.072 0.052

この反応のΔEのmax-minは0.45 kcal/mol，標準偏差は0.12でした．この反応はDFTだとΔEに関しては定性的にも誤った結果を与える汎関数が多いのですが，構造はどの方法でも問題なく求まるようです．そもそも特殊な構造をもつ化合物ではありませんので当然かもしれませんが．この計算においては構造最適化の計算レベルは反応のエネルギーに与える影響は小さく，重要なのは全電子エネルギーの計算方法であると考えられます．

有機金属の反応

有機系ではほとんど差が出ないのでそろそろ飽きて来ましたが，遷移金属が入った反応を検討しました．スキーム3に示すPd(dmpe)錯体（Pdゼロ価錯体）へのクロロメタンの酸化的付加です．計算結果は表3にまとめます．

スキーム3．Pd(dmpe)錯体へのクロロメタンの酸化的付加．
表3．酸化的付加のΔE，ΔΔEと各分子の構造のRMSD（SCS-MP2/def2-TZVP基準）．
構造レベル ΔE [kcal/mol] ΔΔE [kcal/mol] RMSD (Pd(dmpe)) [Å] RMSD (MeCl) [Å] RMSD (PdCl(dmpe)Me) [Å] 平均RMSD [Å]
SCS-MP2/def2-TZVP -47.06 - - - - -
SCS-MP2/def2-SVP -47.31 0.25 0.061 0.008 0.032 0.033
PBE-D4/def2-SVP -47.64 0.58 0.161 0.015 0.046 0.074
TPSS-D3(BJ)/def2-SVP -47.74 0.68 0.046 0.013 0.152 0.070
B3LYP-D4/det2-SVP -47.65 0.59 0.052 0.011 0.155 0.073
PBE0-D3(BJ)/def2-SVP -47.49 0.43 0.141 0.007 0.037 0.062
PBE0-D4/def2-SVP -47.49 0.43 0.141 0.007 0.036 0.061
PBE-D4/def2-TZVP -47.63 0.57 0.100 0.007 0.031 0.046
B3LYP-D3(BJ)/det2-TZVP -47.45 0.39 0.096 0.004 0.038 0.046
B3LYP-D4/det2-TZVP -47.43 0.37 0.091 0.004 0.037 0.044
PBE0-D3(BJ)/def2-TZVP -47.46 0.40 0.079 0.004 0.030 0.038
PBE0-D4/def2-TZVP -47.49 0.43 0.077 0.004 0.030 0.037
PBEh-3c -47.33 0.27 0.107 0.007 0.038 0.051
B97-3c -47.65 0.59 0.121 0.007 0.033 0.054

スキーム3の酸化的付加のΔEのmax-minは0.68 kcal/mol，標準偏差は0.15で，再びどの汎関数，基底関数でもほぼ同じ結果になりました．構造のRMSDが0.16とか0.14程度の計算結果もありますが，エネルギー的には1 kcal/mol以下の差異です．

反応の活性化エネルギー

それでは遷移状態ならばどうでしょう．計算したのはスキーム4に示したトリエチルアミンとクロロエタンとが反応して活性錯合体が形成される際の活性化エネルギーです．なお，計算にあたっては塩素原子のみにdiffuse関数を含むma-def2-XVP（X = S or TZ）を用いました．計算結果を表4にまとめます．

スキーム4．トリエチルアミンとクロロエタンの反応．
表4．トリエチルアミンとクロロエタンの反応のΔE‡，ΔΔE‡と各分子の構造のRMSD（SCS-MP2/def2-TZVP基準）．
構造レベル ΔE‡ [kcal/mol] ΔΔE‡ [kcal/mol] RMSD (Et3N) [Å] RMSD (EtCl) [Å] RMSD (TS) [Å] 平均RMSD [Å]
SCS-MP2/(ma-)def2-TZVP 30.59 - - - - -
SCS-MP2/(ma-)def2-SVP 30.57 0.027 0.017 0.087 0.008 0.037
PBE-D4/(ma-)def2-SVP 30.00 0.60 0.062 0.093 0.020 0.058
TPSS-D3(BJ)/(ma-)def2-SVP 29.91 0.68 0.051 0.094 0.018 0.054
B3LYP-D3(BJ)/(ma-)det2-SVP 30.39 0.20 0.051 0.093 0.016 0.053
PBE0-D3(BJ)/(ma-)def2-SVP 30.91 -0.32 0.112 0.078 0.009 0.066
PBE0-D4/(ma-)def2-SVP 30.54 0.054 0.036 0.083 0.009 0.043
PBE-D4/(ma-)def2-TZVP 29.56 1.03 0.044 0.054 0.010 0.036
B3LYP-D3(BJ)/(ma-)det2-TZVP 30.39 0.20 0.040 0.050 0.007 0.032
B3LYP-D4/(ma-)det2-TZVP 30.33 0.26 0.045 0.059 0.008 0.037
PBE0-D3(BJ)/(ma-)def2-TZVP 30.46 0.14 0.030 0.032 0.005 0.022
PBE0-D4/(ma-)def2-TZVP 30.50 0.091 0.031 0.035 0.005 0.024
PBEh-3c(+ma) 30.44 0.16 0.035 0.034 0.006 0.025
B97-3c(+ma) 30.05 0.55 0.046 0.031 0.008 0.028

遷移状態だと多少分布が広がり，ΔE‡のmax-minは1.34 kcal/mol，標準偏差は0.35となりました．エネルギー的に一番差異があるのはΔΔE‡ = 1.03 kcal/molのPBE-D4/(ma-)def2-TZVPですが，構造のRMSDは特に悪いわけではありません．そして1 kcal/mol程度の差異はエネルギー計算の方法に何を選択するかで変わってくるような差異です．

構造の平均RMSDとΔΔEの分布

構造を求める手法と反応のエネルギーの間には明確な関係はなさそうですが，最後に平均RMSDとΔΔEがどのように分布しているかを確認しておきます．

図1．構造の平均RMSDとΔΔEの分布．

0.40で構造最適化の計算レベルがエネルギー計算に与える影響はやはり小さいようです．

結論

様々な計算方法で構造最適化を行い，得られた構造を用いてDLPNO-CEPA/1/(ma-)def2-SVPレベルで反応のΔEまたはΔE‡を求めて比較しました．また，構造最適化計算で得られた構造のRMSDの比較も行いました．今回試した反応においては構造のRMSDは最大で0.16程度ありましたが，この程度の構造の差では反応のエネルギーに対する影響はほとんどありませんでした．

反応の各種エネルギーを議論するにあたっては，今回試した限りでは構造最適化の計算レベルの影響は小さく，分散力補正したDFTとsplit-valence（double-zeta）基底関数で十分であると考えられます．一方で全電子エネルギーを求める方法の影響は大きいことが知られています．ORCAが使えるならばDLPNO-CCSD(T)がベストでしょう．

計算方法

全ての計算にはORCA 4.2.1を利用しました．構造最適化ではRI近似またはRIJCOSX近似を適用し，補助基底関数にはdef2/Jを用いました．DFTの積分グリッドはGrid6，COSXの積分グリッドはGridX7を使用しました．TightOptキーワードを用いて構造最適化を行いました．得られた構造を元にNumFreqで振動数計算を行って，平衡構造（もしくは遷移状態）であることを確認しました．反応のΔEの計算は各レベルで得られた構造を元に，DLPNO-CEPA/1/def2-SVPレベルで計算を行いました．SCF計算にはdef2/Jを補助基底関数を用いたRIJCOSX近似を適用し，COSXの積分グリッドはGridX8としました．アニオンを含む計算では，アニオン部分にのみdiffuse関数（minimally augmented def2 basis set）を使用しました．構造のRMSDはCalculate Root-mean-square deviation (RMSD) of Two Molecules Using Rotationを用いて求めました．

© 2020 化ケ屋