密度フィッティングと恒等分解(RI)

密度フィッティングと恒等分解(RI)#

2 電子積分テンソル \((\mu\nu|\lambda\sigma)\) は \(O(n_{\text{ao}}^4)\) にスケールします — ほとんどの 量子化学計算における単一最大のコストです。密度フィッティング(**恒等分解(RI)**とも呼ばれる)は、 すべての 4 添字積分を、はるかに小さな 3 添字因子にわたる縮約で近似し、\(O(n_{\text{ao}}^4)\) の対象を 制御可能な小誤差で \(O(n_{\text{fit}}\cdot n_{\text{ao}}^2)\) のものに変えます。本章はこの近似自体、 qc-rs があらゆる RI バックエンドで共有する「因子 \(B\)」に変換するために使う厳密な計量の白色化、そして その単一因子が平均場 J/K とポスト HF 相関の両方にどう奉仕するかを導出します。 .design/54two-electron-integral-spine.md と .design/68ri-jk-strategies.md に基づき、AGENTS.md の統一 eri= RI 戦略記述と一致させます。

恒等分解近似#

厳密な AO 積密度 \(\rho_{\mu\nu}(\mathbf r)=\chi_\mu(\mathbf r)\chi_\nu(\mathbf r)\) は一般にどんな有限 補助基底でも厳密に展開できませんが、フィッティング関数 \(\{\varphi_P\}\) の線形結合で近似できます:

\[ \rho_{\mu\nu}(\mathbf r) \approx \tilde\rho_{\mu\nu}(\mathbf r) = \sum_P c^{\mu\nu}_P\,\varphi_P(\mathbf r). \]

Dunlap の変分 RI は、フィッティング係数 \(c^{\mu\nu}_P\) を、残差 \(\rho_{\mu\nu}-\tilde\rho_{\mu\nu}\) の 実空間ノルムではなくCoulomb 自己エネルギーを最小化するように選びます — \(J\) と \(K\) の精度を決める のは点ごとの密度誤差ではなく Coulomb エネルギーだからです。 \(\iint\frac{(\rho-\tilde\rho)(\mathbf r)(\rho-\tilde\rho)(\mathbf r')}{|\mathbf r-\mathbf r'|}\,d\mathbf r\,d\mathbf r'\) を \(c^{\mu\nu}_P\) について最小化すると、停留条件 \(\sum_Q(P|Q)c^{\mu\nu}_Q=(P|\mu\nu)\) が得られます、すなわち

\[ c^{\mu\nu}_P = \sum_Q (V^{-1})_{PQ}\,(Q|\mu\nu), \qquad V_{PQ} = (P|Q) = \iint \frac{\varphi_P(\mathbf r)\varphi_Q(\mathbf r')}{|\mathbf r-\mathbf r'|}\,d\mathbf r\,d\mathbf r', \]

\((Q|\mu\nu)\) は 1 つの補助関数と 1 つの AO 対の間の 3 中心 Coulomb 積分、\(V\) は補助基底だけにわたる 2 中心 Coulomb 計量です。このフィットをどんな 4 添字積分にも代入すると、あらゆる RI バックエンドが 実際に使う近似が得られます:

\[ (\mu\nu|\lambda\sigma) \approx \sum_{PQ}(\mu\nu|P)\,(V^{-1})_{PQ}\,(Q|\lambda\sigma). \]

これは(場当たり的な打ち切りではなく)真の変分最小化から来るので、RI 誤差はなめらかで、補助基底を 拡大することで系統的に改善可能であり、適切に対応するフィッティング基底(cc-pvdz-jkfit、 def2-svp-jkfit など、軌道基底ごとに表になっている)では経験的に微小です — これが、RI が単にメモリ 逼迫時に頼るだけの安い近似ではなく、qc-rs の推奨既定平均場経路である理由です。

白色化:計量を単一の 3 添字因子にまとめる#

\((V^{-1})_{PQ}\) を各 J/K 縮約の中で直接扱うと、Fock 構築ごとに \(n_{\text{fit}}\times n_{\text{fit}}\) 演算子を再適用することになります。qc-rs は代わりに、3 中心積分を一度だけ白色化し、計量全体を 単一の 3 添字因子に折り込みます:

\[ B^P_{\mu\nu} = \sum_Q\bigl(V^{-1/2}\bigr)_{PQ}\,(Q|\mu\nu), \qquad (\mu\nu|\lambda\sigma) \approx \sum_P B^P_{\mu\nu}\,B^P_{\lambda\sigma}, \]

これは \(V^{-1}=V^{-1/2}V^{-1/2}\) を使って \(\sum_{PQ}(\mu\nu|P)(V^{-1})_{PQ}(Q|\lambda\sigma)\) を書き 換えたものにほかなりません — 非対称な分解ではなく対称な平方根分解を選ぶことで、\(B\) が計量全体を運び、 下流のすべての縮約が \(P\) にわたる単純な内積になり、ループの中に計量が残らなくなります。qc-rs は \(V^{-1/2}\) を 2 通りの等価な方法で計算します:\(V\) の固有分解(ゼロに近い固有値を ri.metric_cutoff IOP 未満で落とす — 大きい/拡散した補助基底では計量が軽度にランク落ちしうる)、あるいは Cholesky 解 \(V=LL^{\mathsf T}\)、\(B=L^{-1}\cdot(\cdot)\)(固有分解なし、数値的により安定、その場で実行)— どちら も同じ演算子について同じ縮約済み \((\mu\nu|\lambda\sigma)\) を与えるので、どちらを使うかは物理の 選択ではなく数値安定性の選択です。

検証済み例 — 水での RI 近似誤差、そして 2 つの白色化実装の厳密な一致(ri-ram = B を保存、 ri-recomp = サイクルごとに B を再計算、どちらも同じ数式から構築):

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

m4c = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").ints(eri="4c-incore").scf(ref="r").run()
m4c.scf.energy   # -76.02679364497408 -- 厳密な 4 中心参照

mri = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").ints(eri="ri-ram").scf(ref="r").run()
mri.scf.energy   # -76.02677274798874
mri.scf.energy - m4c.scf.energy   # 2.09e-05 Ha -- RI フィッティング誤差、自動導出された jkfit で

mri2 = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").ints(eri="ri-recomp").scf(ref="r").run()
mri2.scf.energy - mri.scf.energy   # -1.4e-14 -- ri-ram と ri-recomp は浮動小数点簡約の底まで一致

厳密な 4 中心参照に対する \(2\times10^{-5}\ E_h\) のフィッティング誤差は、手動調整のない既定補助基底 としては典型的な RI-JK の値であり — 背後の軌道基底自体の精度をはるかに下回ります。明示的な rijk= 補助基底が与えられない場合、qc-rs は既定の JK フィッティング基底を未設定のままにするのではなく AO 基底から自動導出します(default_jk_aux)— 要求された厳密な方法を静かに RI へ「昇格」させる 「auto-RI」は存在しませんが、RI が明示的に要求されればいつでも妥当な補助基底の既定値があります。

白色化因子からの \(J\) と \(K\)#

\(B\) が与えられれば、Coulomb と交換の両方は密度行列との単純な縮約になり、\(B\) が常駐していれば (RAM でもディスク/RMA からストリームされていても — ストレージバックエンドはポリシーであり、この 数式を変えません)完全にローカルです:

\[ J_{\mu\nu} = \sum_P B^P_{\mu\nu}\,q_P, \qquad q_P = \sum_{\lambda\sigma}B^P_{\lambda\sigma}D_{\lambda\sigma}, \]
\[ K_{\mu\nu} = \sum_P\sum_i B^P_{\mu i}\,B^P_{\nu i}, \qquad B^P_{\mu i}=\sum_\nu B^P_{\mu\nu}C_{\nu i} \ \text{(占有 MO へ半変換済み)}. \]

\(J\) は小さな補助空間ベクトル \(q_P\) だけを必要とします(密度に対する縮約 1 回、\(B\) に対する縮約もう 1 回);\(K\) はまず \(B\) を占有 MO 添字へ半変換する必要があります(次章で 相関法が使う AO→MO 半変換とまったく類似)、その後外積の和です。どちらも完全な 4 添字テンソルには 決して触れません。これが、大規模での通常の 4 中心交換より RI が \(K\) について著しく安い理由でもあり ます — 半変換は外積の前に片方の AO 添字を \(n_{\text{occ}}\) に縮めます、\(n_{\text{ao}}^2\) を直接 \(n_{\text{ao}}^2\) と縮約するのではなく。

\(B\) を得る 2 つの方法:再計算対保存#

qc-rs の本番 RI 機構は、eri= の値で選ばれる1 つの μ 分散機構に 2 つの \(L\) ソースモード (MuRiJk/MuRiJkDistributed)です:

  • **ri-recomp**は何も保存しません:SCF サイクルごとに libcint から 3 中心 \(L=(Q|\mu\nu)\) を 再計算し、常駐する白色化 \(W=V^{-1/2}\) とその場で縮約します — 完全な \(B\) テンソルを一切材料化 しません、CPU(サイクルごとの再評価)とメモリ(小さな計量因子だけが常駐)をトレードします。 ri-recomp-disk はさらに厳しいメモリ予算向けのアウトオブコア兄弟であり、サイクルごとの中間物を チェックポイントの scratch ディレクトリにスプールします。

  • ri-ramは白色化因子 \(\tilde B=W\cdot L\) を一度だけその場で構築し、その後のすべての SCF サイクルで常駐させ続けます — サイクルごとの \(J\) と \(K\) の両方が同じ常駐 \(\tilde B\) を読み、 サイクルごとの 3 中心再計算はまったくありません。これは、完全な因子をメモリに保持するコストで サイクルあたりより高速です。

どちらも数学的に同一です — 上の検証済み例が直接示す通り浮動小数点簡約順序まで互いにビット一致します (ri-ram と ri-recomp は \(\sim10^{-14}\) だけ異なり、これは通常の浮動小数点ノイズの底であって 物理的な違いではありません)— 両者の選択は純粋なメモリ/CPU のトレードオフであり、精度のトレード オフではありません。ri-ram は大規模での incore 並みの速度の既定であり、ri-recomp は最大の基底の ためのメモリ節約の選択であり、SCF 自体が保存したか再計算したかに関わらず相関法が同じ 3 中心構成要素を 必要とするため MoTransformSource が直接再利用するものでもあります。

分散計量なしで \(B\) を分散する#

分散メモリ RI を一般に難しくする唯一の要因は、素朴な方式では計量 \(V^{-1/2}\) を、3 中心テンソルを保持 するどんな分散方式にもまたがって適用する必要がある — SCF サイクルごとの全対全通信または遠隔メモリ 行列ベクトル積 — ことです。qc-rs の設計は、\(B\) が SCF 不変(積分だけから一度構築され、すべての SCF 反復にわたって変更なく再利用される)であることを利用してこれを完全に回避し、再分散コストは構築時 に一度だけ払われます:

  1. 計量は分散配列ではなく複製配列です。 \(V\) とその因子(\(V^{-1}\)/\(V^{-1/2}\)/Cholesky 解演算子)は \(O(n_{\text{fit}}^2)\) — \(O(n_{\text{fit}}\cdot n_{\text{ao}}^2)\) の 3 中心テンソルに比べて微小 — なので、各ランクは単に自分の完全なコピーを持ちます(Workspace RMA プール の dup 配列の種類:get/put はローカルメモリコピーで MPI ゼロ)。これほど小さいものを複製 するのは無料です;これがそれ以降のすべてのステップをローカルにするものです。

  2. 白色化は一度だけ、ローカルに、1 回限りの転置の間に適用されます。 各ランクの生の 3 中心テンソル のスラブは、複製された計量因子に対するローカル GEMM で白色化されます — 反復SCF ループのどこにも 分散計量行列ベクトル積はなく、この 1 回限りの構築パス中にのみ(せいぜい)存在します。

  3. 白色化因子は補助添字 \(P\) で分散され、ノードごとに完全な AO 対を持ちます。 各ランクは \(P\) の 連続した行ブロックを所有し、自分の行については完全な \(\mu\nu\) 範囲を保持します。これがまさに \(J\)(所有行のローカル GEMM を \(D\) に対して行い、その後小さな \(O(n_{\text{ao}}^2)\) の Allreduce) と \(K\)(所有行の半変換、その後ローカル外積)の両方を各ランクで完全にローカルにするものです — SCF サイクルごとのランク間通信は Fock 行列の 1 回の小さな Allreduce だけであり、分散テンソル 縮約は決してありません。

メモリは線形に集約します:各ランクは \(B\) の自分の \(P\) 行分だけを保持するので、1 ノードの RAM に 大きすぎる \(B\) テンソルを持つ分子も、数ノードにまたがれば快適に収まります — AGENTS.md の qc-hpc/RMA 規則が記述する、主要な分散メモリターゲットである≤4ノードのメモリ集約体制です(大規模 HPC スケール アウトではありません)。マルチランク \(K\) 構築自体にはさらなる選択、ri.ram.kdist IOP キーがあります: "p-transpose"(既定 — 1 回限りの MPI_Alltoallv が \(\tilde B\) をフィット添字分散レイアウトへ転置 するので、以降のすべての SCF サイクルの \(J\)/\(K\) はランクローカルにプラス小さな \(\operatorname{Allreduce}(n_{\text{ao}}^2)\) が 1 回)または "mu-stream"(μ 分散、サイクルごとに 半変換を Gram ストリームする代わり)— サイクルごとの通信対 1 回限りの転置のトレードオフであり、 精度のトレードオフではありません。

        flowchart TD
    L["生の 3 中心 L = (Q|mu nu)<br/>(libcint, 補助殻ごと)"] --> WHITEN["白色化: B = V^-1/2 . L<br/>(複製計量、ローカル GEMM)"]
    WHITEN --> DIST["B を補助添字 P で分散<br/>(ランクごとに完全な mu-nu)"]
    DIST --> J["J: q_P = sum B^P.D (所有行)<br/>J = sum_P B^P q_P<br/>+ Allreduce(nao^2)"]
    DIST --> K["K: B を占有 MO へ半変換<br/>K = sum_P sum_i B^P_mu-i B^P_nu-i<br/>+ Allreduce(nao^2)"]
    DIST --> MO["相関: AO->MO 半変換<br/>B^P_pq (任意の MO 空間対)"]

同じ因子が相関にも奉仕する#

RI-MP2 が必要とする AO→MO 半変換は、密度行列の代わりに MO 係数と縮約する だけで、同じ白色化因子から構築されます:

\[ B^P_{pq} = \sum_{\mu\nu}C^L_{\mu p}\,B^P_{\mu\nu}\,C^R_{\nu q} \]

任意の左/右 MO 空間対について(RI-MP2 では占有×仮想ですが、同じ機構が将来のどんな相関法にも占有×占有 や仮想×仮想として奉仕します)。qc-rs の設計が意図的に利用する 1 つの機微:白色化と変換の順序は 入れ替え可能で、「先に白色化してから変換」(上記の J/K 形式)と「先に変換してから白色化」はまったく 同じ最終 \(B^P_{pq}\) を与えます。白色化は補助添字だけへの線形演算であり、軌道添字への AO→MO 縮約と 可換だからです。相関は変換してから白色化を使います:まず \(n_{\text{ao}}^2\to n_{\text{left}}\cdot n_{\text{right}}\) に縮める(占有/仮想空間が \(n_{\text{ao}}\) からかけ離れていればずっと小さな対象)、 それから計量を適用します — 完全な AO 対テンソルを先に白色化するより \(5\)〜\(20\) 倍少ない白色化 FLOP であり、対応してより小さな分散簡約です。これが、RI-MP2 の章 の \(B^P_{ia}\) が別途導出される対象ではない理由です — それは本節の \(B\) を、もう一度変換したものです。

練習 9

  1. 検証済み例は、ri-ram と ri-recomp が \(\sim10^{-14}\) で一致する一方、両者とも厳密な 4 中心 から \(\sim2\times10^{-5}\) 異なることを示しています。これらが完全に異なる種類の「誤差」— 1 つは 真の物理近似、もう 1 つは浮動小数点ノイズ — である理由と、どちらか一方だけがより大きな補助基底で 縮小する理由を一文で説明しなさい。

  2. 分子や AO 基底がどれほど大きくなっても計量 \(V\) を常に dup(完全複製、MPI ゼロ)配列として 扱えるのに対して、3 中心テンソル \(B\) にはそれができないのはなぜですか。

  3. 相関法は AO→MO 変換の後に計量を適用し(「変換してから白色化」)、平均場 J/K はその前に適用 します(「白色化してから変換」)。どちらも同じ演算の数学的に厳密な並べ替えです。相関が平均場とは 逆の順序を好む実務上の(数学上ではない)理由は何ですか。

次の 2 つの章はこの因子の上に直接築かれます:解析的微分 と 線形応答理論 はどちらも、この \(B\) と RI 整合な微分積分を必要とし、 RI-MP2 は上で示した通り、それをもう一度縮約したものにほかなりません。