分散補正の理論:DFT-D3/D4

分散補正の理論:DFT-D3/D4#

Hartree-Fock もほとんどの密度汎関数も、長距離 van der Waals(London 分散)引力を捉えません — HF はそもそも相関を持たず、半局所/ハイブリッド DFT 汎関数は局所密度とその勾配から構築されるため、 離れたフラグメント上の電子間の瞬間多重極相関についての情報を一切持ちません。Grimme の DFT-D3とDFT-D4分散補正は、これを幾何依存(D3)あるいは幾何と電荷の両方に依存する(D4) 経験的エネルギー項として追加します。原子核位置(と、D4 では自己無撞着に決定される原子電荷)だけ から完全に計算され — 追加の積分もなく、SCF ループへの参加もまったくありません。本章は両モデルを 完全に導出します、crates/qc-disp/src/energy.rs、ncoord.rs、c6.rs、eeq.rs、d4.rs に直接基づき — Grimme グループの s-dftd3/dftd4/multicharge/mctc-lib の忠実な、独立に 実装された Rust 移植です。

なぜ幾何だけで相関効果を近似できるのか#

十分離れた 2 つのフラグメント間の London 分散は、物理的には相関ゆらぎ効果です — 一方のフラグメント 上の瞬間双極子が他方に相関した瞬間双極子を誘起し、それらの引力は最低次で \(1/R^6\) で減衰します (より高次の多重極が \(1/R^8\)、\(1/R^{10}\) … を寄与)。D3 と D4 の両方が利用する鍵となる経験的事実は、 この引力の主要係数 \(C_6\) が実際には元素ごとの固定定数ではないということです — それは原子の局所 的な化学環境に強く依存し(sp³ 炭素と芳香族炭素は測定可能に異なる有効分極率を持ち、したがって異なる \(C_6\) を持ちます)、その環境は純粋に幾何学的な何か — 原子がいくつの隣接原子を持ちどれだけ近いか — によってよく捉えられます。これが、分散係数を推定するために実際の電子構造計算を走らせる必要なしに 配位数を補間変数として使うことの、根拠のすべてです。

D3 配位数#

原子 \(A\) について、D3 配位数はなめらかで飽和する関数を使って実効的な隣接原子を数えます、硬い カットオフではなく:

\[ \text{CN}_A = \sum_{B\ne A}\frac{1}{1+\exp\bigl[-k_{\text{cn}}(r_{AB}^{\text{cov}}/r_{AB}-1)\bigr]}, \qquad r_{AB}^{\text{cov}} = r_A^{\text{cov}}+r_B^{\text{cov}},\qquad k_{\text{cn}}=16, \]

\(r_A^{\text{cov}}\) は表になっている D3 共有結合半径です。この計数関数は、共有結合半径の和の内側に 十分ある隣接原子については \(\approx1\)(真に結合した隣接原子)、それよりずっと外側にある隣接原子 については \(\approx0\)(無視できるほど遠い)、共有結合半径の閾値ちょうどでは正確に \(\tfrac12\) です — なめらかなシグモイドであり硬いステップではありません、これが \(\text{CN}_A\)(したがってそこから 構築される補間された \(C_6\))を幾何の微分可能な関数にするものです、解析的分散勾配に必要です。

CN 補間された \(C_6\):表になった参照値へのガウス類似度#

D3 は実行時に物理的な分極率式から \(C_6\) を計算するのではなく — 各元素の小さな参照分子について一度 オフラインで計算された、参照配位数/\(C_6\) 対の表の間を補間します。表になった参照配位数 \(\{\text{cn}^{\text{ref}}_a\}\) を持つ元素 \(Z\) について、実際の配位数 \(\text{CN}\) での参照 \(a\) への 重みは、すべての参照にわたって正規化されたガウス類似度カーネルです:

\[ w^a(\text{CN}) = \frac{\exp\bigl[-w_f(\text{CN}-\text{cn}^{\text{ref}}_a)^2\bigr]} {\sum_k\exp\bigl[-w_f(\text{CN}-\text{cn}^{\text{ref}}_k)^2\bigr]}, \]

なので、表になった配位数が原子の実際の配位数に近い参照はほぼすべての重みを得て、配位数空間で 遠く離れた参照はほとんど何も寄与しません — まさに 1 変数のなめらかな最近傍補間です。(すべての ガウス関数がゼロにアンダーフローする場合 — 極端な、表の範囲外の配位数 — 実装は 0/0 の不定形の 重みを生む代わりに、最高の表になった配位数を持つどの参照にも全重みを与えるフォールバックをします。) 原子対 \(C_6\) 係数は、両原子の参照表にわたる重み付き二重和です:

\[ C_6^{AB}(\text{CN}_A,\text{CN}_B) = \sum_{a,b}w^a(\text{CN}_A)\,w^b(\text{CN}_B)\,C_{6,ab}^{\text{ref}}(Z_A,Z_B). \]

Becke-Johnson 減衰 2 体エネルギー#

生の \(1/R^6\) と \(1/R^8\) の項は近距離で発散します、そこでは分散の描像自体が破綻し(フラグメントは もはや十分離れていない)、SCF 自身の交換反発が既にその物理を扱っているはずです。Becke-Johnson (BJ)有理減衰— qc-rs の既定 — は、汎関数ごとにわずか 2 つのフィットパラメータ(\(a_1,a_2\); \(s_6=1\) がほぼ常に、\(s_8\) は汎関数固有)を持つ汎関数非依存の数学形式を使って両方の項をなめらかに 減衰させます:

\[ E_{\text{disp}}^{(2)} = -\sum_{A>B}C_6^{AB}\left[\frac{s_6}{R_{AB}^6+R_0^6} + \frac{s_8\,\overline{R_{AB}^2}}{R_{AB}^8+R_0^8}\right], \qquad R_0 = a_1\sqrt{3\,\langle R^4\rangle_A\langle R^4\rangle_B/\langle R^2\rangle_A\langle R^2\rangle_B} + a_2, \]

\(\overline{R_{AB}^2}=3\langle R^4\rangle_A\langle R^4\rangle_B/(\langle R^2\rangle_A\langle R^2\rangle_B)\) は、\(R_0\) も構築する同じ表になった元素ごとの比(有効半径二乗の尺度 \(\langle R^4\rangle/\langle R^2\rangle\))であり、したがって減衰半径と \(C_8/C_6\) 比の両方が同じ静的な、 元素対固有の量から来ます — \(R^{-8}\) 項の相対的重みを正しくするために \(a_1,a_2,s_8\) を超える追加の フィットパラメータは不要です。有理減衰の鍵となる性質は、古い「ゼロ減衰」方式とは異なり、 \(R_{AB}\to0\) につれて \(E_{\text{disp}}^{(2)}\) が有限の非ゼロ定数に近づく(ゼロに近づくのではない) ことです — これは、非常に近距離では背後の物理的引力が実際には消えるわけではなく単に全エネルギーの 他の(反発性の)項に支配されるだけなので、分散エネルギーが人為的に消えるのではなくなめらかに飽和 すべきことを反映しています。

3 体 Axilrod-Teller-Muto 項#

2 体分散が常に全体像というわけではありません — 3 つの原子が同時に近接している場合、真の 3 体誘起 双極子相関効果(Axilrod-Teller-Muto、ATM)が、典型的には小さいものの無視できない追加の補正を 寄与します:

\[ E_{\text{disp}}^{(3)} = \sum_{A>B>C}f_{\text{damp}}(\bar R)\;C_9^{ABC}\,\frac{3\cos\theta_A\cos\theta_B \cos\theta_C+1}{(R_{AB}R_{BC}R_{CA})^3}, \qquad C_9^{ABC} = -s_9\sqrt{|C_6^{AB}C_6^{AC}C_6^{BC}|}, \]

\(\theta_A,\theta_B,\theta_C\) は 3 つの原子が形成する三角形の内角です(コンパクトな三角形配置を ほぼ直線的な配置から区別する純粋に幾何学的な角度因子)、そして零型減衰関数 \(f_{\text{damp}}(\bar R)=1/(1+6(R_0/\bar R)^{(\alpha+2)/3})\) は 2 体項が使うのと同じ対ごとの van der Waals 半径 \(R_0\) から構築され、\(\bar R=\sqrt[3]{R_{AB}^2R_{BC}^2R_{CA}^2}\)、\(s_9=1\) が標準的な 3 体スケーリングです(\(s_9=0\) とすると ATM を完全に無効化します)。\(C_9\) 係数の積の幾何平均形式 — \(-s_9\sqrt{|C_6^{AB}C_6^{AC}C_6^{BC}|}\)、別途表になった 3 体参照ではない — それ自体が、まったく 新しい 3 体参照データの集合を必要とするのではなく、既に計算済みの対ごとの \(C_6\) 値を再利用する 近似です。

検証済み例 — 水の二量体、PBE/STO-3G、D3(BJ) ありとなし:

import qc
water_dimer = """
O  -1.551007  -0.114520   0.000000
H  -1.934259   0.762503   0.000000
H  -0.599677   0.040712   0.000000
O   1.350625   0.111469   0.000000
H   1.680398  -0.373741  -0.758561
H   1.680398  -0.373741   0.758561
"""
m0 = qc.chk.new(atom=water_dimer, ao="sto-3g", unit="angstrom").scf(ref="r", xc="pbe").run()
m0.scf.energy   # -150.46818414059214

m1 = qc.chk.new(atom=water_dimer, ao="sto-3g", unit="angstrom").scf(ref="r", xc="pbe", dispersion="d3bj").run()
m1.scf.energy, m1.scf.dispersion   # (-150.46956323055147, -0.0013790899594410507)

分散補正は、この水素結合二量体について全エネルギーを \(\sim1.4\ \text{m}E_h\) 下げます(安定化させ ます)— この粗い基底では絶対値としては小さいですが、PBE 単独では完全には捉えない分子間引力に 対して定性的に正しい符号と桁数です。

D4:固定配位数から自己無撞着原子電荷へ#

D3 の配位数は化学環境の純粋に幾何学的な代理であり、構成上、幾何学的には同一だが分子全体の電荷 状態や電気陰性度バランスのために化学的には異なる 2 つの原子 — 真の電子的効果 — を区別できません。 D4 は D3 を一般化し、まず電気陰性度均衡化(EEQ)によって原子部分電荷を解き、次にそれら の電荷を(配位数に加えて)使って \(C_6\) を補間します。真の(単純化されているとはいえ)電子構造 ステップを、それ以外は高速で幾何駆動の補正である枠組みに折り込みます。

EEQ 電荷モデル自体は、分子ごとの小さな線形系を解きます — SCF ではなく、基底関数積分もなく、表に なった原子パラメータと分子幾何だけです — 固定された全電荷 \(q_{\text{tot}}\)(通常は分子電荷)の 制約下で電気陰性度均衡化エネルギーを最小化する原子電荷 \(q_A\) の集合を求めます。線形系は、自然な \(n\times n\) の電気陰性度方程式に、電荷保存制約を課す 1 つの追加行/列を加えます:

\[\begin{split} \begin{pmatrix}A & \mathbf 1\\ \mathbf 1^{\mathsf T} & 0\end{pmatrix} \begin{pmatrix}\mathbf q\\ \lambda\end{pmatrix} = \begin{pmatrix}\mathbf x\\ q_{\text{tot}}\end{pmatrix}, \qquad A_{ii} = \eta_i + \sqrt{2/\pi}/\rho_i, \qquad A_{ij} = \frac{\operatorname{erf}\bigl(r_{ij}\sqrt{\gamma_{ij}} \bigr)}{r_{ij}}, \qquad \gamma_{ij} = \frac{1}{\rho_i^2+\rho_j^2}, \end{split}\]

\(x_i=-\chi_i+k_{\text{cn}\chi,i}\sqrt{\text{CN}_i}\)(元素ごとの電気陰性度 \(\chi_i\)、異なる配位数と 結合する — D4 モデル自身の誤差関数 CN、上記の D3 のシグモイド CN とは正確な関数形式が異なります が、背後の物理的発想は同じです)、\(\eta_i\) は化学硬度パラメータ、\(\rho_i\) は各原子の電荷密度がどれ だけ「広がっている」と見なされるかを設定する電荷密度半径です(\(\operatorname{erf}(\cdot)/r\) の 非対角形はまさに、2 つのそのような広がった Gauss 型電荷分布の間の Coulomb 相互作用であり、裸の 点電荷 \(1/r\) ではありません — これが原子間距離が短くても方程式を良条件に保つものです)。Lagrange 乗数 \(\lambda\) は \(\sum_iq_i=q_{\text{tot}}\) を厳密に強制します。解かれると、D4 はこれらの電荷を — 自身の配位数と並んで — 参照 \(C_6\) 表への2 番目の補間軸として使い、各原子の実効分散係数に、 局所幾何環境と電子的電荷状態の両方への感度を与えます。

検証済み例 — D3(BJ) の代わりに D4 を使った同じ水二量体:

m2 = qc.chk.new(atom=water_dimer, ao="sto-3g", unit="angstrom").scf(ref="r", xc="pbe", dispersion="d4").run()
m2.scf.energy, m2.scf.dispersion   # (-150.46918114898963, -0.0009970083974916537)

ここでの D4 の補正(\(\sim-1.0\ \text{m}E_h\))は D3(BJ) のもの(\(\sim-1.4\ \text{m}E_h\))より絶対値 が小さいです — 2 つのモデルは異なる(ただし重なりのある)参照データセットに対して独立に較正され ており、厳密に一致することは期待されず、それぞれが同じ定性的物理を捉えることだけが期待されます: 分散は、まったく分散補正のない同じ幾何と比べて水素結合二量体を安定化させます。

分散が SCF に対してどこに位置するか#

D3 も D4 も、SCF 自己無撞着ループの完全に外側で評価されます — それらは原子核位置(D3)、あるいは 原子核位置と一度限りの EEQ 解(D4)にのみ依存し、決して収束済み電子密度自体には依存せず、Fock 行列 への項を一切寄与しません。これは、真に密度依存であり毎 SCF サイクルに参加しなければならず(線形 応答の章 に従えば)応答カーネルにも参加しなければならない PCM 溶媒和項 の構造的な対極です。分散のエネルギーは単に、(SCF 収束の後に、あるいはそれとは独立に)一度だけ加え られます — これが、分散項の解析的勾配がそれ自身の軌道応答を一切必要 としない理由でもあります:それは核座標だけの関数であり、直接微分され、Hellmann-Feynman の機微も 結合摂動解法もまったく関与しません。

練習 15

  1. D3 の配位数は、硬い「結合/非結合」カットオフではなくなめらかなシグモイド計数関数を使います。 硬いカットオフが結果の分散エネルギーを幾何の関数として不連続にする理由、そしてそれが特に幾何 最適化と解析的勾配にとって深刻な問題になる理由を一文で説明しなさい。

  2. D4 は EEQ 電荷を、配位数の代わりにではなく配位数に加えて、\(C_6\) の 2 番目の補間軸として 使います。局所幾何環境は同一だが分子的な文脈が異なる 2 つの原子について、D3 の配位数のみの 補間ではそれらを区別できないが D4 の追加された電荷軸ならできる、という具体的なシナリオを 1 つ 挙げなさい。

  3. 分散の解析的勾配が(先に導出した)解析的 Hessian や PCM 反応場 応答とは異なり CPHF/結合摂動解法を一切必要としないのはなぜですか。

本章は qc-rs の Part II「Hamiltonian 修正」の 3 つの章 — ECP、 PCM 溶媒和、そして分散 — を完結させます。それぞれ独立に切り替え可能で、 どの SCF 参照や交換相関汎関数とも組み合わさります。次章はエネルギー論から解析へと向かいます: 電荷/結合次数理論 が、収束済みの密度と MO から化学的意味をどう 読み出すかを導出します。