【B3LYP】DFT初心者が知っておくべき厳選ポイントを紹介します

2024年9月20日

B3LYPはDFT法で最もよく利用される汎関数なので、「とりあえずB3LYPで」と量子化学シミュレーションを始める場合は多いかと思います。

しかし、あまりにスタンダードな汎関数であるため、初心者にとっては「DFTとB3LYPの違いが判らない」、「前例踏襲でインプットに「B3LYP/6-31G」と書いている」という方も多いのではないでしょうか？

このページでは「B3LYPは経験的なパラメータを含む」、「ハートリーフォック交換積分を混ぜている」など、ほかの人から少し突っ込まれても困らないように、B3LYPで抑えておくべき基礎的なポイントをまとめています。

【ポイント１】B3LYPは最頻出のDFTの汎関数

そもそもDFT（密度汎関数法）とは、電子密度$\rho(\boldsymbol{r})$によって電子的エネルギー$E$などを求める手法です。そのエネルギー$E$と電子密度$\rho(\boldsymbol{r})$を紐づける関係式の一つがB3LYPであり、汎関数$E_{\text{ex}}$と呼ばれるものです。

【ポイント２】B3LYPは最も古い混成汎関数

DFT法では系のエネルギーは汎関数$E_{\text{ex}}$によって、どの程度現実と近いエネルギーを与えるかがシミュレーションの精度を決めます。

これは汎関数$E_{\text{ex}}$によってHF法と現実のずれを埋め合わせるためです。

そのため、DFT法では様々な$E_{\text{ex}}$が開発・利用されてきました。

B3LYPはDFT汎関数の発展において、下のSTEP３の混成汎関数に属する汎関数であり、かつ、混成汎関数の最初期に提案されたものです。

STEP

局所密度近似（LDA）汎関数（1920年代～）

DFT法で最も歴史があり、また電子密度$\rho$を使ってエネルギーを計算するという密度汎関数法のコンセプトに最も忠実なのが局所密度近似（LDA）汎関数です。

シュレーディンガー方程式が提案された翌年の1927年にDFT法は提案され、FermiやDiracらによって間を置かずにLDA汎関数も提案されています。

STEP

一般化勾配近似（GGA）汎関数（1980年代後半～）

分子など、電子密度$\rho(\boldsymbol{r})$が急激に変化する対象の計算精度を上げるため、汎関数$E_{\text{ex}}$に$\rho(\boldsymbol{r})$だけでなく、勾配$\nabla\rho(\boldsymbol{r})$依存性も加えるのが一般化勾配近似（GGA）汎関数です。

GGA以降の汎関数は物理的な前提（負値条件など）などの条件をできるだけ多く・かつ実際の実験値と合うように作成されてきています。

STEP

混成汎関数（1990年代前半～）

GGAまでの汎関数は長距離部分で電子密度が指数的に減衰するため再現性が悪いものがあったことから、汎関数に一定の割合でハートリー・フォック交換積分を混合したのが混成汎関数です。

混成汎関数は物理的な正当性が薄いこと、パラメータの多いことなどの問題点はありますが、分子物性の再現性が高いことから、多くのDFT計算で利用されています。

【ポイント３】B3LYPは複数の汎関数を混ぜて作られている

最初期に作成された混成汎関数であるB3LYP汎関数は既存の汎関数を複数混ぜて作られたものです。

$E_\text{xc}^\text{B3LYP}=(1-a)E_{x}^\text{S}+aE_x^\text{HF}+$$b(E_x^\text{B88}-E_x^\text{S})$$+cE_c^\text{LYP}+(1-c)E_c^\text{VWN}$

混成汎関数であるB3LYPが初めて登場したことによって、小分子に対するDFT計算の平均二乗誤差が$4\text{kcal}/\text{mol}$程度以下まで抑えることができました。

有名なモデルだけどすっきりした関係式ではないんだね。

B3LYPはあくまで小分子の計算結果が合うように経験的に作られた関係式に過ぎません。

【特徴１】パラメーターが3つある

まず、B3LYPの「3」はパラメーターが３つ（$a$、$b$、$c$）あることに由来します。

$E_\text{xc}^\text{B3LYP}=(1-a)E_{x}^\text{S}+$$a$$E_x^\text{HF}+$$b$$(E_x^\text{B88}-E_x^\text{S})+$$c$$E_c^\text{LYP}+(1-c)E_c^\text{VWN}$

B3LYPではハートリー・フォック交換積分の混合割合$a=0.20$、Becke交換汎関数の混合割合$b=0.72$、LYP相関汎関数の混合割合$c=0.81$の3つのパラメータがあります。

交換汎関数と相関汎関数の意味は【特徴３】で紹介します

これらは小分子の計算誤差が小さくなるように設定されており、物理的な正当性が薄いデメリットの原因です。

【特徴２】ハートリーフォック交換積分が混ぜ合わされている

次のポイントは汎関数にハートリーフォック交換積分が混ぜ合わされている点です。

$E_\text{xc}^\text{B3LYP}=(1-a)E_{x}^\text{S}+$$aE_x^\text{HF}$$+b(E_x^\text{B88}-E_x^\text{S})$$+cE_c^\text{LYP}+(1-c)E_c^\text{VWN}$

B3LYPの登場前の汎関数は化学結合を強く見積もり過ぎるなどの問題があったことから、汎関数に一定の割合でハートリー・フォック交換積分を混合したのが混成汎関数です。

混成汎関数にはB3LYP、PBE0、M06等がありますが、ハートリー・フォック交換積分の混合割合はそれぞれ$20\%$、$25\%$、$27\%$です。

【特徴３】B3LYPは交換汎関数と相関汎関数に分けられる

B3LYPでは交換汎関数$E_x$と相関汎関数$E_c$が利用されています。

$E_\text{xc}^\text{B3LYP}=$$(1-a)E_{x}^\text{S}+$$aE_x^\text{HF}$$+b(E_x^\text{B88}-E_x^\text{S})$$+$$cE_c^\text{LYP}+(1-c)E_c^\text{VWN}$

この交換汎関数は静的電子相関の影響を含んでおり、相関汎関数は動的電子相関の影響を含んでいると解釈されています（常田貴夫著　「密度汎関数法の基礎」4章を参照）。

静的電子相関（平行スピン由来）: エネルギーが近い状態が混ざり合うことによる安定化エネルギーのことで、金属錯体や化学結合が開裂する際などでは特に重要です。
動的電子相関（反対スピン由来）: HF法が平均場近似を利用して電子間反発エネルギーを取り扱ったことにより、余分に取り込んでしまった電子相互作用のことです。

交換汎関数はBeckeの交換汎関数$E_x^\text{B88}$、相関汎関数はLee達の$E_c^\text{LYP}$などが利用されているんだね。

汎関数M06などのMinnesotaシリーズも様々な交換・相関汎関数を多数のパラメータでつなぎ合わせて作られています

【特徴４】B3LYPにはBeckeの交換汎関数が利用されている

B3LYPに利用されているBeckeの交換汎関数B88（1988年）は非常に有名です。

$E_x^{\text{B88}}=E_x^{\text{S}}-\zeta\sum_\sigma^{\alpha\beta}\int\frac{({\rho_{\sigma}})^{4/3}{(x_\sigma)}^2}{1+6\zeta{x_\sigma}\sinh^{-1}x_\sigma}$

ただし、$x^\sigma=\frac{|\nabla\rho_\sigma|}{\rho_\sigma^{4/3}}$です。

Beckeの交換汎関数B88は先ほどご紹介した一般化勾配近似（GGA）汎関数のうち、最初に提唱されたものです。$\nabla\rho$を汎関数に含めるGGAは密度でエネルギーを紐づけるというDFTのコンセプトから逸脱しますが、原子や分子におけるエネルギーの再現性がそれまでの汎関数よりも飛躍的に高められました。（平均二乗誤差が$8\text{kcal}/\text{mol}$程度以下。）

実際にどのように利用されるか

DFT計算はつまるところフォック行列$\boldsymbol{F}$を対角化することで量子化学計算を行っています。$\boldsymbol{F}=\boldsymbol{H}+\boldsymbol{J}+\boldsymbol{v}_{\text{xc}}$

このうち$\boldsymbol{v}_{\text{xc}}$行列が汎関数によって与えられるもので$E_x^\text{B88}=\int{e^\text{B88}}(\rho(\boldsymbol{r}),\nabla\rho(\boldsymbol{r}))$としたとき、${v}_{\text{xc}}(\boldsymbol{r})=\frac{\partial{e^\text{B88}}}{\partial\rho(\boldsymbol{r})}-\nabla\cdot\frac{\partial{e^\text{B88}}}{\partial\nabla\rho(\boldsymbol{r})}$と与えられます。

例えば$\frac{\partial{e^\text{B88}}}{\partial\rho(\boldsymbol{r})}$（を一部簡略化したもの）は次のように計算できます。

つまり、各グリッド点$\boldsymbol{r}$における$\rho(\boldsymbol{r})$と$\nabla\rho(\boldsymbol{r})$がわかっていれば$v(\boldsymbol{r})$もたやすく計算できるので、基底関数$\chi_\mu(\boldsymbol{r})$と$\chi_\nu(\boldsymbol{r})$による交換汎関数の表現行列の各成分は$\boldsymbol{v}^{\text{xc}}_{\mu\nu}=\int\chi_\mu(\boldsymbol{r})v_{\text{xc}}(\boldsymbol{r})\chi_\nu(\boldsymbol{r})\text{d}\boldsymbol{r}$を数値積分することで求めることができます。