溶媒和と分散補正

溶媒和と分散補正#

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

陰的溶媒和（PCM）#

理論#

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

使い方#

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 補正は、これを安価な幾何のみの原子対和として加え戻します、

\[ E_{\text{disp}} = -\sum_{A<B} \Big( s_6\frac{C_6^{AB}}{R_{AB}^6} + s_8\frac{C_8^{AB}}{R_{AB}^8} \Big) f_{\text{damp}}(R_{AB}), \]

表化された分散係数 \(C_n^{AB}\) と、短距離（汎関数が既に相互作用を記述する領域）で補正を切るダンピング関数を伴い ます。パラメータは汎関数に依存します。

使い方#

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}")

練習 7

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

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

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

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