溶媒和と分散補正

溶媒和と分散補正#

気相の HF や DFT 計算には 2 つの物理効果が欠けており、qc-rs はどちらにも安価で標準的な補正を持ちます。 溶媒和は分子を真空ではなく溶媒中に置きます。分散は多くの汎関数が外す長距離ファンデルワールス引力を 加えます。両者は独立で — どちらか、両方、どちらもなし、が選べ — どちらも自動的に勾配へ 流れ込みます。

陰的溶媒和（PCM）#

理論#

化学の多くは真空ではなく溶液中で起こります。溶媒分子を 1 つずつ明示的にモデル化するのは高価です;陰的（連続体） モデルは、溶媒を 1 つの数 — 誘電率 \(\varepsilon\)（水で ≈ 78.4、真空で ≈ 1）— で特徴づけられる分極可能な 誘電連続体で置き換えます。溶質はその連続体に彫られた分子形の空洞に座り、その電荷分布が誘電体を分極させ、 それが今度は溶質に働き返す反応場を作ります。**分極連続体モデル（PCM）**はその相互分極を自己無撞着に解きます。 反応場が溶質密度に依存し逆もまた然りなので、PCM は追加の Fock 項として SCF ループに直接組み込まれます。

具体的には、空洞表面は（GePol により）小さな面積要素 — テッセラ — に離散化され、反応場は みかけ表面電荷（ASC） \(\mathbf q\)（テッセラごとに 1 値）として捉えられます。これは各テッセラでの溶質の 静電ポテンシャル \(\mathbf v\) に対する境界積分方程式を解いて求めます。qc-rs が実装する 2 つの定式化は、 誘電境界条件をどう課すかで異なります:

\[ \textbf{IEF-PCM（既定）:}\quad \mathbf T(\varepsilon)\,\mathbf q = -\mathbf R_\infty\,\mathbf v, \qquad \mathbf T(\varepsilon) = \big[2\pi f(\varepsilon)\,\mathbf I - \mathbf D\mathbf A\big]\mathbf S, \quad f(\varepsilon)=\frac{\varepsilon+1}{\varepsilon-1}, \quad \mathbf R_\infty = 2\pi\mathbf I - \mathbf D\mathbf A, \]
\[ \textbf{C-PCM（COSMO スケーリング）:}\quad \mathbf S_{\text{scaled}}\,\mathbf q = -\mathbf v, \qquad \mathbf S_{\text{scaled}} = \mathbf S\Big/\frac{\varepsilon-1}{\varepsilon+k}, \]

\(\mathbf S\) と \(\mathbf D\) は誘電体の Green 関数から作られる単層・二重層境界積分演算子（既定の外側真空 では \(1/r\)）、\(\mathbf A\) はテッセラ面積の対角行列です。\(\mathbf q\) が分かれば分極エネルギーは単に \(E_{\text{pol}}=\tfrac12\,\mathbf v\cdot\mathbf q\)。\(\mathbf v\) は現在の SCF 密度から作られるので、この 線形系は毎 SCF サイクル解かれます — これこそ PCM が収束後に適用する一発補正ではなく追加の Fock 項として 実装されている理由です。空洞構築と Green 関数の完全な詳細は 陰的溶媒和の理論 にあります。

使い方#

scf(...) に pcm= を渡します:素の数は誘電率、辞書はモデルを選びます:

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

gas = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run()
sol = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", pcm=78.39).run()

print(f"gas   E = {gas.scf.energy:.6f}")     # -76.026794
print(f"water E = {sol.scf.energy:.6f}")     # -76.036749
print(f"ΔG_solv = {(sol.scf.energy - gas.scf.energy)*627.509:.3f} kcal/mol")   # -6.247

溶媒和したエネルギーはより低い — ここでは誘電体が分子を −6.25 kcal/mol 安定化します（静電的溶媒和自由 エネルギー）。辞書形式は定式化を選びます:

scf(..., pcm={"model": "iefpcm", "epsilon": 78.39})   # IEF-PCM（既定）
scf(..., pcm={"model": "cpcm",   "epsilon": 78.39})   # C-PCM

  • IEF-PCM（積分方程式形式、既定）と C-PCM（導体様）が 2 つの標準定式化で、ごくわずかに異なるエネルギーを 与えます（ここでは −76.036749 対 −76.036812）。

  • cuda ビルドでは "device": "cuda"（または "auto"）を加えると AO↔空洞結合を GPU で走らせ、エネルギーは 同一です。

PCM は HF と KS-DFT に合成でき、解析的勾配も含まれるので、溶媒中で構造最適化できます。

分散補正（DFT-D3 / D4）#

理論#

London 分散 — 瞬間的に誘起された双極子間の弱い引力的ファンデルワールス力 — は純粋な長距離相関効果です。 Hartree–Fock はこれを完全に外し、多くの半局所 DFT 汎関数もほとんど捉えないので、π スタック・層状物質・ホスト– ゲスト錯体・分子結晶を過小結合します。DFT-D 補正は、これを安価な幾何のみの原子対和として加え戻します。

既定の Becke–Johnson（有理）ダンピングつき D3 は、正確には

\[ E_{\text{disp}}^{\text{D3(BJ)}} = -\sum_{A<B} C_6^{AB} \left( \frac{s_6}{R_{AB}^6 + R_0^6} + \frac{s_8\,\big(C_8^{AB}/C_6^{AB}\big)}{R_{AB}^8 + R_0^8} \right), \qquad R_0 = a_1\sqrt{C_8^{AB}/C_6^{AB}} + a_2 . \]

物理的に異なる 3 つの要素がこれに入ります:

  • \(C_6^{AB}\) は固定の原子定数ではありません — 各原子の配位数（滑らかな距離ベースの切替関数による 連続的な近傍数）から補間されます。だから \(\text{sp}^3\) 環境の炭素と芳香族炭素は異なる \(C_6\) を得ます。 この局所環境への感受性が D3 の「3」の由来です。

  • \(C_8^{AB}/C_6^{AB} = 3\sqrt{\langle r^4\rangle_A/\langle r^2\rangle_A}\cdot\sqrt{\langle r^4\rangle_B/\langle r^2\rangle_B}\) は表化された元素ごとの半径から。

  • \(s_6,s_8,a_1,a_2\) は 4 つの Becke–Johnson パラメータで、DFT 汎関数ごとに再フィットされています （dispersion(...) に渡す汎関数名で引かれます）。

有理（BJ）ダンピングの美点は、分母の \(+R_0^6\) そのものにあります:指数型の切替関数と違い、 \(R_{AB}\to0\) でも決して発散しないので、短く強く結合した接触でも補正は有限（かつ小さい、\(C_6/R_0^6\) は 有界なので）に保たれ、別途の短距離カットオフ処理は不要です。d3zero は代わりに古い指数型の切替形式を使います （だから下の表でコンパクトな水分子への寄与がほぼゼロなのです）;d4 は配位数のみの \(C_6\) を、幾何依存の EEQ 部分電荷で置き換え、局所環境への感受性をさらに精密化します。

使い方#

mychk.dispersion(functional, method=...) は SCF エネルギーに加えるべき分散エネルギーを返します （幾何のみなので SCF は不要です）:

mol = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom")

mol.dispersion("pbe")                      # -0.00035947   D3(BJ)、既定
mol.dispersion("pbe", method="d3zero")     # -0.00000363   D3 ゼロダンピング
mol.dispersion("pbe", method="d3bjatm")    # -0.00035947   D3(BJ) + ATM 3 体項
mol.dispersion("pbe", method="d4")         # -0.00019523   D4（EEQ 電荷 + ATM）

total = mol.scf(ref="r", xc="pbe").run().scf.energy + mol.dispersion("pbe")   # 分散補正した全エネルギー

method=

補正

d3bj (既定)

Becke–Johnson ダンピングの D3 — 標準的な選択

d3zero

ゼロダンピングの D3（古いダンピング形式）

d3bjatm

D3(BJ) + Axilrod–Teller–Muto 3 体項

d4

D4 — 幾何依存電荷（EEQ）+ ATM

寄与するのは実原子だけ — ゴーストとダミー原子は除外されます。水の分散エネルギーはごく小さい（分散結合の 接触がない）です;補正が効くのは、より大きな弱結合の集合体で、数 kcal/mol になりえます。ここで d3bjatm が 表示桁で d3bj に等しいのは、水の 3 体（ATM）項が無視できるためで — 密で分極しやすい系では大きくなります。

Tip

SCF ワークフロー内の分散 scf(...) に dispersion=... を渡すと（例 scf(ref="r", xc="pbe", dispersion="d3bj")）、補正が報告される scf.energy と勾配に自動的に折り込まれます — 分散補正した構造最適化にはこれが欲しい形です。上の単体 mychk.dispersion(...) は、エネルギー項そのものを得る幾何のみのアクセサです。

総合例:両方の補正#

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"

gas   = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc="pbe").run().scf.energy
solv  = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc="pbe", pcm=78.39).run().scf.energy
disp  = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").dispersion("pbe")

print(f"gas-phase PBE      : {gas:.6f}")
print(f"+ solvation (water): {solv:.6f}   (ΔG_solv = {(solv-gas)*627.509:.2f} kcal/mol)")
print(f"+ dispersion (D3BJ): {gas + disp:.6f}")

練習 8

  1. 気相の反応エネルギーを計算したら水溶液の実験と食い違いました。本章のどの補正を最初に加えるのが自然で、 どう要求しますか。

  2. mychk.dispersion(...) は .scf().run() を先に要さないのに、pcm= は scf(...) の中に入れねばならないのは なぜですか。

  3. 水では分散補正は ~0.0004 Ha（~0.2 kcal/mol）です。それがはるかに大きくなると予想される系を 1 つ挙げ、理由を 述べなさい。

次は可視化で、これらの密度や軌道を確認できる図にします。