# 溶媒和と分散補正

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

## 陰的溶媒和（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 関数の完全な詳細は
[陰的溶媒和の理論](../10-foundations/solvation-theory.md) にあります。

### 使い方

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

```python
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 安定化します（静電的溶媒和自由
エネルギー）。辞書形式は定式化を選びます:

```python
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 は不要です）:

```python
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(...)` は、エネルギー項そのものを得る幾何のみのアクセサです。
:::

## 総合例:両方の補正

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

:::{exercise}
:label: ex-solv-disp

1. 気相の反応エネルギーを計算したら水溶液の実験と食い違いました。本章のどの補正を最初に加えるのが自然で、
   どう要求しますか。
2. `mychk.dispersion(...)` は `.scf().run()` を先に要さないのに、`pcm=` は `scf(...)` の*中*に入れねばならないのは
   なぜですか。
3. 水では分散補正は ~0.0004 Ha（~0.2 kcal/mol）です。それがはるかに大きくなると予想される系を 1 つ挙げ、理由を
   述べなさい。
:::

:::{solution} ex-solv-disp
:class: dropdown

1. **陰的溶媒和** — 実験は水中、計算は真空です。`scf(...)` の呼び出しに `pcm=78.39`（または
   `pcm={"model":"iefpcm","epsilon":78.39}`）を加えます。
2. 分散（D3/D4）は原子位置と汎関数の**幾何のみ**の関数で — 電子密度に依存しない — ので波動関数が不要です。PCM の
   反応場は密度に**依存**し、SCF と*ともに*自己無撞着に解かれるので、`scf(...)` の一部でなければなりません。
3. ファンデルワールス力で結合するもの:π スタック二量体（例 ベンゼン 2 分子）、層状物質（黒鉛）、ホスト–ゲスト
   錯体、分子結晶。そこでは分散引力*こそ*が結合なので、補正は数 kcal/mol に達し、系が結合するかどうかを左右し
   えます。
:::

次は[可視化](visualization.md)で、これらの密度や軌道を確認できる図にします。