# ポスト Hartree-Fock 相関理論

Hartree-Fock は各電子に他の電子の*平均*場を与えます — その系統的誤差である**相関エネルギー**は
$E_{\text{corr}}=E_{\text{exact}}-E_{\text{HF}}$ と定義され、典型的には全エネルギーの$\sim1\%$ですが
$\gg1\ \text{kcal/mol}$ であり、化学的に決定的です。本章は Møller-Plesset 2 次摂動論(MP2)— qc-rs が
実際に実装している相関法 — をその**RI/密度フィッティング**形式で導出します。これが qc-rs が出荷する
*唯一*の形式だからです(`AGENTS.md` によれば「これらの相関法は RI-C のみ」)。`.design/25qc.lct.mp2.md`
と `crates/qc-ints-engine/src/trafo/mp2.rs` に基づきます。

## Møller-Plesset 摂動論:設定

MP2 は Hartree-Fock ハミルトニアンをゼロ次問題として、真の電子ハミルトニアンと Fock 演算子の和との差を
摂動として扱います:

$$
\hat H = \hat H_0 + \hat V, \qquad \hat H_0 = \sum_i\hat f(i), \qquad
\hat V = \hat H_{\text{elec}} - \hat H_0,
$$

$\hat f$ は Fock 演算子です(その固有関数がまさに正準 HF 軌道である — これが、任意の参照ではなく HF が
自然なゼロ次である理由です)。ゼロ次エネルギーは占有軌道エネルギーの和で、1 次補正はちょうど HF
エネルギーを回復します — $E^{(0)}+E^{(1)}=E_{\text{HF}}$ — なので**MP2 は相関を回復する最初の次数**、
すなわち 2 次補正 $E^{(2)}$ です。標準的な Rayleigh-Schrödinger 摂動論は

$$
E^{(2)} = \sum_{n\ne0}\frac{|\langle\Psi_0|\hat V|\Psi_n\rangle|^2}{E_0^{(0)}-E_n^{(0)}}.
$$

$\hat V$ が 2 電子演算子なので、$\Psi_0$ と非ゼロの行列要素を持つのは*2 電子励起*行列式
$\Psi_{ij}^{ab}$(占有軌道 $i,j$ を仮想 $a,b$ で置き換える)だけです — 単一励起は Brillouin の定理
(HF 軌道は既に単一励起-HF の結合をゼロにしている)で消え、3 重励起以上は $\hat V$ が高々 2 つの
スピン軌道が異なる行列式しか結合しないため消えます。2 電子行列要素とエネルギー分母(各仮想/占有ペアが
エネルギーギャップに $\varepsilon_a-\varepsilon_i$ を寄与する)を丁寧に展開すると、**閉殻スピン適応
MP2 エネルギー**が得られます:

$$
E_{\text{MP2}} = \sum_{i\le j}\sum_{ab}\frac{(ia|jb)\bigl[2(ia|jb)-(ib|ja)\bigr]}
{\varepsilon_i+\varepsilon_j-\varepsilon_a-\varepsilon_b},
$$

$(ia|jb)$ は化学者の記法での、占有($i,j$)・仮想($a,b$)正準 MO にわたる 2 電子積分です。これはまさに
qc-rs の作業方程式(`rimp2_energy` のドキュメントコメント)であり、$j\le i$ の占有ペアにわたって和を
取り、オフダイアゴナルペアには暗黙の因子 2 が付きます(被加数は $i\leftrightarrow j$ 対称)。

## 逆スピン/同スピン分解

分子の 2 つの項を分けると、qc-rs の API と出力全体で使われる自然な物理的分解が得られます:

$$
E_{\text{OS}} = \sum_{ijab}\frac{(ia|jb)^2}{\Delta_{ijab}}, \qquad
E_{\text{SS}} = \sum_{ijab}\frac{(ia|jb)\bigl[(ia|jb)-(ib|ja)\bigr]}{\Delta_{ijab}}, \qquad
E_{\text{MP2}} = E_{\text{OS}} + E_{\text{SS}},
$$

$\Delta_{ijab}=\varepsilon_i+\varepsilon_j-\varepsilon_a-\varepsilon_b$ が分母です。$E_{\text{OS}}$ は
Coulomb 孔(逆スピン電子対が互いを避ける)の寄与、$E_{\text{SS}}$ は Fermi 孔(同スピン、HF に組み込まれた
Pauli 排他律により既に部分的に分離されている)の寄与です — 経験的に、正準 MP2 は真の相関エネルギーに
対して $E_{\text{SS}}$ を系統的に**過大評価**し、$E_{\text{OS}}$ をそれより軽く過小評価します。これが、
qc-rs が素の MP2 と並んで実装する 2 つのスピンスケール変種の動機です:

$$
E_{\text{SCS-MP2}} = c_{\text{os}}E_{\text{OS}} + c_{\text{ss}}E_{\text{SS}}, \qquad
(c_{\text{os}},c_{\text{ss}}) = \Bigl(1.2,\ \tfrac13\Bigr)\ \text{(Grimme, 2003)},
$$

$$
E_{\text{SOS-MP2}} = c_{\text{os}}'E_{\text{OS}}, \qquad c_{\text{os}}' = 1.3\ \text{(Head-Gordon, 2004;
qc-rs の既定)},
$$

SOS-MP2 は同スピン項を単に軽く重み付けするのではなく完全に捨てます — これが、次節で示すように
$O(N^5)$ ではなく $O(N^4)$ の評価を可能にするものです。

## なぜ RI のみか:密度フィッティングされた 3 添字因子化

直接 4 添字 $(ia|jb)$ の評価は、AO→MO 積分変換だけで $O(N^5)$ かかります(4 回連続する
$O(N^4)$ スケーリングの半変換)、そして 4 中心 ERI テンソルをカルテットごとに材料化するかその場で
構築する必要があります。qc-rs は代わりに、あらゆる 2 電子積分を [RI/密度フィッティング機構](density-fitting-ri.md)
の白色化された 3 添字 MO 因子で因子化します:

$$
(ia|jb) = \sum_P B^P_{ia}\,B^P_{jb}, \qquad B^P_{pq} = \sum_Q\bigl(V^{-1/2}\bigr)_{PQ}(Q|pq),
$$

こうして各占有-仮想対積分は(はるかに小さい)補助添字 $P$ 上の内積 1 回になり、MP2 エネルギー全体は
生の 4 添字変換ではなく `nfit`$\times$`nvir` GEMM の列になります。これが `AGENTS.md` が明確に
qc-rs の相関法は**RI-C のみ**だと述べる理由です — 4 中心 MP2 変換経路は存在せず計画もされていません;
相関フィッティング用補助基底(`ric=` 引数、SCF 自身の `rijk=` 補助基底とは別物 — 両者は異なる精度/
目的のトレードオフに応えており、通常異なる公開基底セットです。例えば `cc-pvdz-jkfit` ではなく
`cc-pvdz-ri/mp2fit`)を与える必要があります。

:::{prf:algorithm} 正準 RI-MP2 ペアループ
:label: alg-rimp2-pairloop-ja

**入力:** ([RI 半変換](density-fitting-ri.md) で一度だけ構築される)白色化占有×仮想 RI 因子 $B^P_{ia}$、
正準軌道エネルギー $\varepsilon_i,\varepsilon_a$。
**出力:** $E_{\text{OS}}, E_{\text{SS}}$(したがって $E_{\text{MP2}}$、$E_{\text{SCS-MP2}}$、
$E_{\text{SOS-MP2}}$)。

1. $j\le i$ を満たす各占有ペア $(i,j)$ について(重みは $j\ne i$ なら 2、さもなくば 1):
   1. $I_{ab} = (ia|jb) = \sum_P B^P_{ia}B^P_{jb}$ を作る — 1 回の GEMM
      $(n_{\text{vir}}\times n_{\text{fit}})\cdot(n_{\text{fit}}\times n_{\text{vir}})$。
   2. $(ib|ja)$ は単に $I_{ba}$(同じ行列の添字を入れ替えたもの)— 2 回目の GEMM は不要。
   3. $E_{\text{OS}}\mathrel{+}= \text{重み}\sum_{ab}I_{ab}^2/\Delta_{ijab}$ および
      $E_{\text{SS}}\mathrel{+}= \text{重み}\sum_{ab}I_{ab}(I_{ab}-I_{ba})/\Delta_{ijab}$ を累積する。
2. $(E_{\text{OS}}, E_{\text{SS}})$ を返す。
:::

占有ペアループは驚くほど並列です — 各ペアは共有の不変な $B$ 因子と軌道エネルギーだけを読み、自分の
スレッドローカルな $(E_{\text{OS}},E_{\text{SS}})$ 累積器だけに書くので、qc-rs はこれを
`qc_hpc::thread_count()` 上に直接スレッド化し、ワーカーごとに 1 つのスクラッチバッファをループの
**外側で**確保します(コードベースの他所と同じ、ホットループ非確保の規律)。

## Laplace 変換:$O(N^4)$ での SOS-MP2

SOS-MP2 が同スピン項を完全に捨てるという決定は、漸近的により安価なアルゴリズムを可能にします。分母の
逆数を**数値 Laplace 変換**として書きます:

$$
\frac{1}{\Delta_{ijab}} \approx \sum_\alpha w_\alpha\, e^{-\Delta_{ijab}\,t_\alpha}, \qquad
\Delta_{ijab} = (\varepsilon_a-\varepsilon_i)+(\varepsilon_b-\varepsilon_j) > 0,
$$

軌道エネルギーギャップの物理的範囲にわたって $1/\Delta$ を相対誤差 $<10^{-9}$ で再現するよう校正された
少数の求積点/重み $(t_\alpha,w_\alpha)$ を使います(qc-rs はまず$\sim8$〜$12$点の minimax 最適化
フィットを試み、その最適化フィットが許容誤差に届かなければ$\sim40$点の tanh-sinh/二重指数グリッドに
フォールバックします — いずれにせよ最終エネルギーはマイクロハートリー精度です)。指数関数は 4 つの
軌道添字にわたって**分離**します:$e^{-\Delta_{ijab}t}=
e^{-(\varepsilon_a-\varepsilon_i)t/2}e^{-(\varepsilon_a-\varepsilon_i)t/2}e^{-(\varepsilon_b-\varepsilon_j)t/2}
e^{-(\varepsilon_b-\varepsilon_j)t/2}$、これにより $(i,j,a,b)$ にわたる和が**因子化**します:Laplace
スケール因子 $\tilde B^P_{ia}(t)=B^P_{ia}\,e^{-(\varepsilon_a-\varepsilon_i)t/2}$ を定義すると、

$$
E_{\text{OS}} = -\sum_\alpha w_\alpha \bigl\|G(t_\alpha)\bigr\|_F^2, \qquad
G^{PQ}(t) = \sum_{ia}\tilde B^P_{ia}(t)\,\tilde B^Q_{ia}(t),
$$

$G(t)$ は結合された占有-仮想添字にわたる単一の GEMM 縮約で構築される $n_{\text{fit}}\times
n_{\text{fit}}$ 行列であり、占有*ペア* $(i,j)$ にわたる明示的なループはまったくありません。$G(t)$ の
構築コストは求積点あたり $O(n_{\text{fit}}^2\,n_{\text{occ}}\,n_{\text{vir}})$ であり、ペアループの
$O(n_{\text{occ}}^2\,n_{\text{vir}}^2\,n_{\text{fit}})$ とは対照的です — $O(N^5)$ スケーリングの占有
ペア和を、それぞれ安価な少数($\sim40$)の求積点にわたる $O(N^4)$ スケーリングの和に交換します。これは
まさに `sos_mp2_laplace` の実装であり、直接ペアループ `rimp2_energy(...).os` と求積精度まで一致する
ことが検証されています — 2 つのアルゴリズムは*同じ*和を計算しており、単に再編成されているだけです。
分散実行では、求積点ごとに小さな $n_{\text{fit}}\times n_{\text{fit}}$ 行列 $G(t)$ だけが `Allreduce`
を必要とします — $n_{\text{occ}}\cdot n_{\text{vir}}$ にスケールするどんなテンソルよりもはるかに
安価な通信です。

## 開殻参照:UHF と ROHF の曖昧さ

**UHF** に微妙な点はありません:各スピンチャネルが自分自身の正準軌道と軌道エネルギーを持ち、同スピン
和は $\alpha\alpha$ と $\beta\beta$ の寄与に分かれ、逆スピン和は純粋な $\alpha\beta$ になります
(交換項がスピンチャネルをまたぐことはない)— `riump2_energy` は閉殻ペアループを独立な $\alpha$/
$\beta$ 軌道空間へ直接一般化したものです。

**ROHF は MP2 に対して真に曖昧**であり、これは実装の細部として片付けるより理解する価値があります:
ROHF は両スピンで共有される 1 つの**有効 Roothaan Fock 演算子**を対角化するので、その共有演算子の
固有値は per-spin の MP2 分母として意味を持ち*ません* — きちんと定義された per-spin Fock 行列
$F^\sigma$ は存在します(それが ROHF エネルギーを構築した)が、それは共有 ROHF MO 基底では対角では
ありません。2 つの弁護可能な選択肢があり、qc-rs は両方を実装します(`lct.rohf_mp2` IOP キーで選択):

| 変種 | IOP 値 | 内容 |
|---|---|---|
| **pyscf**(qc-rs の既定) | `"pyscf"` | 生の ROHF 軌道を保持;共有 ROHF MO 基底での*対角*
  $\operatorname{diag}(C^{\mathsf T}F^\sigma C)$ を per-spin 軌道エネルギーとする。PySCF の
  `mp.MP2(ROHF)`(内部では ROHF 変換軌道上の `UMP2`)と一致。 |
| **qcrs** | `"qcrs"` | **半正準化**:共有 ROHF 軌道の占有ブロック*内*と仮想ブロック*内*で別々に
  $F^\sigma$ を対角化する(占有と仮想は決して混ざらないので ROHF 密度は不変);得られた固有値が
  per-spin 軌道エネルギーになり、軌道自体は各ブロック内で回転する。 |

どちらも正当です — 両者は対角化されないオフダイアゴナル $F^\sigma_{ij}$/$F^\sigma_{ab}$ ブロックの
扱い方だけが異なります(`pyscf` では捨てる、`qcrs` では回転に折り込む)— 結果として得られる相関
エネルギーは近いものの同一ではなく、下の検証例が直接それを示します。凍結コア切り捨て
(`lct.frozen_core`、相関和から除外する最低占有軌道数を表す整数 — 既定 `0`、すなわち全電子)は、
両変種と UHF/RHF に同一に適用され、$B^P_{ia}$ を構築する前に占有ブロックの最低 `n_frozen` 行を
単に落とすだけです。

検証済み例 — 水/cc-pVDZ(相関フィッティング用補助基底として `cc-pvdz-ri/mp2fit`):

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

m = qc.chk.new(atom=water, ao="cc-pvdz", ric="cc-pvdz-ri/mp2fit", unit="angstrom").scf(ref="r").lct(method="mp2").run()
m.scf.energy, m.lct.energy, m.lct.e_corr, m.lct.e_os, m.lct.e_ss
# (-76.02679364497408, -76.23074681436064, -0.20395316938655841, -0.15237768872141289, -0.05157548066514553)

m_scs = qc.chk.new(atom=water, ao="cc-pvdz", ric="cc-pvdz-ri/mp2fit", unit="angstrom").scf(ref="r").lct(method="scs-mp2").run()
m_scs.lct.e_corr   # -0.2000450533540773 == 1.2*(-0.15237768872141289) + (1/3)*(-0.05157548066514553)

m_sos = qc.chk.new(atom=water, ao="cc-pvdz", ric="cc-pvdz-ri/mp2fit", unit="angstrom").scf(ref="r").lct(method="sos-mp2").run()
m_sos.lct.e_corr   # -0.19809099564379767 == 1.3*(-0.15237768872141289)

# 凍結コア:O 1s コア軌道を凍結
m_fc = qc.chk.new(atom=water, ao="cc-pvdz", ric="cc-pvdz-ri/mp2fit", unit="angstrom",
                   iop={"lct.frozen_core": 1}).scf(ref="r").lct(method="mp2").run()
m_fc.lct.e_corr   # -0.20161465527331127 -- 絶対値がより小さい:凍結されたコア軌道は
                  # もう仮想へ相関できない

# ROHF-MP2:CH3 ラジカル、両変種
ch3 = "C 0 0 0; H 0 1.079 0; H 0.934 -0.540 0; H -0.934 -0.540 0"
for variant in ("pyscf", "qcrs"):
    r = qc.chk.new(atom=ch3, ao="cc-pvdz", ric="cc-pvdz-ri/mp2fit", unit="angstrom", spin=2,
                    iop={"lct.rohf_mp2": variant}).scf(ref="ro").lct(method="mp2").run()
    print(variant, r.scf.energy, r.lct.e_corr)
# pyscf -39.55963721124465 -0.13088782816695832
# qcrs  -39.55963721124465 -0.13083446557503792
```

`pyscf`/`qcrs` の ROHF-MP2 相関エネルギーは、ここでは約 $5\times10^{-5}\ E_h$ で一致します — 同じ
曖昧さに対する 2 つの弁護可能な解決として予想される通り近いものの同一ではなく、まさにこれが、この
選択が黙って選ばれる内部詳細ではなく明示的で文書化された IOP キーである理由です。

## 相関がどこで止まっているか

`lct(method=...)` は、SCF 参照と CASSCF 参照のあらゆる相関法をまとめる qc-rs の傘です。今日、実際に
エネルギーを生成する実装は RI-MP2 ファミリー(`mp2`、`scs-mp2`、`sos-mp2`)だけです;`cc2`(SCF 参照の
coupled-cluster 様の方法)と CASSCF 参照の `caspt2`/`nevpt2` は依然**モック**のままです — 同じ保留
ステップの構文を受け付け、決定論的なプレースホルダーエネルギーを生成しますが、まだ物理的に意味のある
ものではありません。[マルチリファレンス展望の章](multireference-outlook.md) が、真の CASPT2/NEVPT2
実装に何が必要か、そしてなぜ CASSCF 自身のモック状態がそれを阻んでいるかを扱います。

:::{exercise}
:label: ex-post-hf-theory-ja

1. SOS-MP2 は、SCS-MP2 のように単に重み付けし直すのではなく、同スピン項を完全に捨てます。Laplace
   変換による導出を使って、$E_{\text{SS}}$ を(他ならぬ)捨てることが指数因子化のトリックと
   $O(N^5)\to O(N^4)$ の高速化を可能にする理由を 1〜2 文で説明しなさい。
2. 検証済み凍結コア例は、酸素の 1s 軌道を凍結すると $|E_{\text{corr}}|$ が絶対値で*減少*することを
   示しています($-0.20395$ 対 $-0.20161\ E_h$)。占有軌道を相関和から除去することが、なぜ(各項の
   符号を固定すれば)$|E_{\text{corr}}|$ を減少させることしかできず、決して増加させないのか説明
   しなさい。
3. `pyscf` と `qcrs` の ROHF-MP2 変種は*SCF*エネルギーでは一致しますが、MP2 相関エネルギーではわずかに
   異なります。なぜこの曖昧さは相関ステップだけに影響し、それに先立つ平均場ステップには影響しないの
   ですか。
:::

:::{solution} ex-post-hf-theory-ja
:class: dropdown

1. 同スピン項 $E_{\text{SS}}$ は反対称化された分子を持つ*交換的*な構造を持ちます、
   $(ia|jb)-(ib|ja)$ — 第 2 項は、Laplace 置換の下で独立な軌道ごとの指数因子の積に分離しない形で、
   $a$ と $b$ を 2 つの占有添字 $i,j$ を*またいで*結合します。対照的に逆スピン項の分子
   $(ia|jb)^2$ は既に同一構造の 2 つの積分の単純な積なので、指数分母は $i,a$ と $j,b$ の部分に
   きれいに因子化し、それらが $G(t)=\sum_{ia}\tilde B_{ia}\tilde B_{ia}$ 縮約へ再結合します。
   $E_{\text{SS}}$ を捨てることが、まさに和から因子化不可能な項を取り除くことです。
2. MP2 二重和の各項 $(ia|jb)[2(ia|jb)-(ib|ja)]/\Delta_{ijab}$ は、固定された占有-仮想の 4 つ組に
   対して独立に寄与し、(MP2 エネルギー式が負定値である物理的に通常の場合)各占有添字 $i$ は仮想への
   励起チャネルを*追加*するだけです。占有軌道を凍結すると、その軌道を $i$ または $j$ として含む和の
   すべての項が除去されます — 全体から既に非負の大きさの寄与を差し引くことしかできず、新しい寄与を
   加えることはできないので、軌道が活性相関空間から除去されると $|E_{\text{corr}}|$ は縮小する
   (あるいは、凍結された軌道がたまたまちょうどゼロを寄与していた場合は同じままである)ことしか
   できません。
3. ROHF の自己無撞着場は、**共有**の Roothaan 有効 Fock 演算子とその変分最小化によって完全に決まり
   ます — 明確に定義された ROHF エネルギーと ROHF 軌道/占有の集合は 1 つしかないので、両変種は
   文字通り*同じ*収束済み平均場状態から出発します(上の同一な `scf.energy` で確認済み)。曖昧さは
   ROHF が*何であるか*にあるのではなく、その定義演算子がスピン共有であるような状態に意味のある
   **per-spin の 1 粒子エネルギー**をどう割り当てるかにあります — これは、スピン分解された*ポスト*
   HF 摂動補正を構築しようとして初めて生じる問い、まさに MP2 ステップであり SCF ステップではあり
   ません。
:::

本章が基盤とする RI 因子 $B^P_{ia}$ は [密度フィッティング/RI の章](density-fitting-ri.md) で完全に
導出されます;[マルチリファレンス展望](multireference-outlook.md) は本章の「モック」な方法が終わる
ところから引き継ぎます。