ポスト Hartree-Fock 相関理論

ポスト 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/密度フィッティング機構 の白色化された 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)を与える必要があります。

アルゴリズム 3 (正準 RI-MP2 ペアループ)

入力: (RI 半変換 で一度だけ構築される)白色化占有×仮想 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):

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 は依然モックのままです — 同じ保留 ステップの構文を受け付け、決定論的なプレースホルダーエネルギーを生成しますが、まだ物理的に意味のある ものではありません。マルチリファレンス展望の章 が、真の CASPT2/NEVPT2 実装に何が必要か、そしてなぜ CASSCF 自身のモック状態がそれを阻んでいるかを扱います。

練習 8

  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 相関エネルギーではわずかに 異なります。なぜこの曖昧さは相関ステップだけに影響し、それに先立つ平均場ステップには影響しないの ですか。

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