DFT 交換相関の数値求積

DFT 交換相関の数値求積#

Kohn-Sham DFT の交換相関項 \(E_{\text{xc}}[\rho]\) は、自明な汎関数を除いてどれも閉形式を持たないため、 実用的な KS-DFT コードはすべてこれ — およびその Fock 寄与 \(V_{\text{xc}}\) — を実空間上の数値求積で 評価します:マイクロハートリー精度まで収束するのに十分細かく、しかし数百 AO の系でも手が届くほど粗い、 点と重みの分子グリッドです。本章はそのグリッド — 原子中心の動径/角度積公式、原子ごとのグリッドを 1 つの 分子グリッドに縫い合わせる Becke 空間分割、点数を精度を落とさずに抑える NWChem 式角度プルーニング、 そして AO 値・密度・XC ポテンシャル行列がブロック単位で実際にどう組み立てられるか — を導出します。 .design/27qc.dftgrid.md と qc-rs 現行の qc-grid 実装(AGENTS.md の qc-grid の項:本番の動径は Treutler-Ahlrichs M4、NWChem 5 領域プルーニング、PySCF-Bragg の Becke 分割)に基づきます。

原子中心積公式としての分子求積#

厳密な XC エネルギーは全空間にわたる 3 次元積分 \(E_{\text{xc}}=\int\rho(\mathbf r)\,\epsilon_{\text{xc}}(\mathbf r)\,d\mathbf r\) です。被積分関数はどの核でも(AO のカスプから)鋭く 尖り、それ以外ではなめらかなので、単一の大域求積ではうまく扱えません。Becke(1988)による解決策は、 各原子核に 1 つの原子中心球面グリッドを張り付け、それらをなめらかな空間分割で貼り合わせることです:

\[ \int f(\mathbf r)\,d\mathbf r \;\approx\; \sum_A\sum_i\sum_k w_{Aik}\,f(\mathbf r_{Aik}), \qquad \mathbf r_{Aik} = \mathbf R_A + r_{Ai}\,\mathbf n_k, \]

ここで \(A\) は原子、\(i\) は \(A\) を中心とする動径殻、\(k\) はその殻の単位球面上の角度方向を走ります。角度 (Lebedev)重みが \(\sum_k w_k^{\text{ang}}=1\) に規格化されていれば、単中心球面殻の重みは \(\tilde w_{Aik}=4\pi r_{Ai}^2\,w_{Ai}^{\text{rad}}\,w_k^{\text{ang}}\) であり、最終的な分子重みは 原子分割関数 \(W_A\) を折り込みます:

\[ w_{Aik} = \tilde w_{Aik}\,W_A(\mathbf r_{Aik}), \qquad \sum_A W_A(\mathbf r) = 1\ \text{(至るところ)}. \]

この単位分割条件が、収束した 1 次元動径/角度則の極限で原子グリッドの貼り合わせを厳密にするもの です — 二重計算も隙間もありません — そして qc-rs のテストスイートが数値的に検証する最初の不変量です (目標 \(10^{-12}\))。

Becke 空間分割#

Becke の分割関数は2 中心双曲座標から構成されます。点 \(\mathbf r\) と 2 つの原子核 \(A,B\) について、

\[ \mu_{AB}(\mathbf r) = \frac{|\mathbf r-\mathbf R_A| - |\mathbf r-\mathbf R_B|}{|\mathbf R_A-\mathbf R_B|} \in[-1,1], \]

これは \(A\) で \(-1\)、\(B\) で \(+1\)、垂直二等分面で \(0\) です。\(\mu_{AB}\) だけの生の切り替え関数は、両原子が どれほど異なっていても分割面をちょうど中点に置いてしまいます(水素の隣の重原子は自分の近くの空間の 分け前が不当に小さくなります)。そこで Becke はBragg-Slater 原子半径 \(R_A^{\text{BS}}\) を使って 補正します:

\[ \chi_{AB} = \frac{R_A^{\text{BS}}}{R_B^{\text{BS}}}, \qquad m_{AB} = \frac{\chi_{AB}-1}{\chi_{AB}+1}, \qquad a_{AB} = \operatorname{clip}\!\left(\frac{m_{AB}}{m_{AB}^2-1},\,-\tfrac12,\,\tfrac12\right), \]
\[ \nu_{AB} = \mu_{AB} + a_{AB}\left(1-\mu_{AB}^2\right). \]

(qc-rs はPySCF-Braggサイズ調整済み半径テーブルを使います。オリジナルの 1988 年 Becke 半径ではなく PySCF の既定と一致させています — 導出は同一で、表になっている \(R_A^{\text{BS}}\) の値だけが異なり ます。)補正済み座標 \(\nu_{AB}\) はその後、反復平滑化多項式を通されます。これは空間分割境界の導関数を 高次まで 0 にし(\(W_A\) の求積精度を損なうカスプを避ける)ためです:

\[ p_0(x)=x, \qquad p_{n+1}(x) = \tfrac32p_n(x) - \tfrac12p_n(x)^3, \qquad s_h(x) = \tfrac12\bigl(1-p_h(x)\bigr), \]

\(h\) はBecke 硬度です(既定 3 — 平滑化多項式を 3 回反復)。原子 \(A\) の非正規化セル関数は他のすべての 原子にわたる切り替え関数の積であり、全原子にわたって正規化すれば上で使った分割関数が得られます:

\[ P_A(\mathbf r) = \prod_{B\ne A} s_h\bigl(\nu_{AB}(\mathbf r)\bigr), \qquad W_A(\mathbf r) = \frac{P_A(\mathbf r)}{\sum_C P_C(\mathbf r)}. \]

対称性により、同種原子 2 つの中点では \(W_A=W_B=\tfrac12\) ちょうどになります — qc-rs のグリッド モジュールが直接検査する回帰テストの 1 つです。ゴースト原子(核を持たない基底関数)とダミー原子は、 Bragg-Slater 半径が核中心の概念であるため、既定では決して分割中心になりません;それでもグリッド上で 自分の基底関数は評価されます、自分専用の動径/角度殻を持たないだけです。

動径求積:Treutler-Ahlrichs M4 写像#

動径グリッドは \(\int_0^\infty F(r)\,dr\approx\sum_iw_i^{\text{rad}}F(r_i)\) を離散化します。qc-rs の 本番動径方式はTreutler-Ahlrichs M4 写像(Treutler & Ahlrichs, J. Chem. Phys. 102, 346 (1995))であり、PySCF の dft.radi.treutler_ahlrichs を \(\xi=1\)(元素非依存 — 原子サイズの調整は ここではなく \(R_A^{\text{BS}}\) 経由で Becke 分割の中で行われる)で忠実に移植したものです。\((-1,1)\) 上の Gauss-Chebyshev 横座標 \(x_i=\cos\theta_i\)(\(\theta_i=\tfrac{i\pi}{n+1}\))から出発し、\(L=1/\ln2\) と おくと、動径ノードと重みは

\[ r_i = -L\,(1+x_i)^{0.6}\ln\!\left(\frac{1-x_i}{2}\right), \]
\[ w_i = \frac{\pi}{n+1}\sin\theta_i \cdot L\,(1+x_i)^{0.6}\left[-\frac{0.6}{1+x_i}\ln\!\left(\frac{1-x_i}{2}\right) + \frac{1}{1-x_i}\right], \]

第 2 項は、素の \(x\) 空間 Chebyshev 則を \(r\) 空間の則に変換するのに必要な連鎖律 Jacobian \(dr/dx\) です。 この写像は、より古い Gauss-Chebyshev + Becke 有理写像方式(RadialGridKind::BeckeGaussChebyshev、 レガシーオプションとして保持)よりも DFT 精度において著しく点効率が良く — より少ない動径殻で同じ積分 電子数の許容誤差に達します — これが qc-rs の本番既定である理由です。(qc-grid クレートは両方式を 解析的 Gauss 動径モーメント \(\int_0^\infty e^{-\alpha r^2}4\pi r^2\,dr=(\pi/\alpha)^{3/2}\) に対して 検証します。)

角度求積と NWChem プルーニング#

各動径殻での単位球面上の角度積分には、生の点数ではなく代数次数(その次数までの球面調和関数に対して 厳密)で選ばれるLebedev-Laikovグリッドを使います:

Lebedev 次数

点数

用途

9

38

最内側のプルーニング領域

11

50

核近傍領域

15

86

内殻領域

29

302

medium プリセットのピーク

35

434

fine プリセットのピーク

131

5810

ultrafine(密、プルーニングなし)のピーク

すべての動径殻でピーク角度次数を使うと、密度がほぼ球対称で低次則が既に厳密に積分できる核近傍、 そして密度が指数的に小さい遠方 tail で、膨大な数の点を浪費します。qc-rs の本番プルーニングは NWChem 5 領域方式です(PySCF/gpu4pyscf 互換;以前の連続的な ORZ 式ランプは非推奨のレガシーで、 API 互換性のためだけに残されています)。Bragg 半径 \(R_A^{\text{BS}}\) の原子について、動径座標は比率 \(r/R_A^{\text{BS}}\) として表され、4 つの元素帯依存カットポイント \(\alpha_1<\alpha_2<\alpha_3<\alpha_4\) (3 つの帯:H/He、Li–Ne、それ以外 — .design/27qc.dftgrid.md の NWCHEM_ALPHAS テーブル)により 5 領域に分割され、5 つの領域は次の角度次数を取ります:

\[ \bigl[\,11,\ 15,\ \text{ladder}[j-1],\ n_{\max},\ \text{ladder}[j-1]\,\bigr], \]

\(j\) は \(n_{\max}\) が固定次数ラダー上のどこに位置するかを示すインデックスです — 物理的に重要な中間 半径の結合領域でのみ密で、核近傍と tail では軽くなります。これは、以前の非推奨ランプ \(n_{\text{ang}}(r_i)=\min\{n\in\mathcal L\mid n\ge\min(\lfloor k_sr_in_{\max}/R_A^{\text{BS}}\rfloor,\, n_{\max})\}\) とは対照的に、連続的な式ではなく離散的な元素帯/領域テーブルです — どちらも同じ発想 (必要な場所でのみ密な角度サンプリング)を表現しますが、NWChem テーブルはなめらかなランプに従う のではなく PySCF/ORCA の参照グリッドと厳密に一致するよう校正されています。

プリセットラダー:coarse → ultrafine#

qc-rs は 4 つの名前付きプリセットを公開しており、コストと精度をトレードオフする ORCA DefGrid 式の ラダーです:

プリセット

動径殻数

ピーク角度次数(点数)

プルーニング?

ultrafine に対する最悪誤差

coarse

50

23(194)

あり

\(\lesssim 8\ \mu E_h\)/原子

medium(既定)

60

29(302)

あり

\(\lesssim 0.4\ \mu E_h\)/原子

fine

75

35(434)

あり

\(\lesssim 0.2\ \mu E_h\)/原子

ultrafine

200

131(5810)

なし(密な参照)

—

プルーニングされないのは ultrafine だけです — これは他の 3 つが校正される密な参照層であり、日常の 本番作業向けの層ではありません。さらに、3d 遷移金属(Sc–Zn)はプルーニング層でより密な動径殻とより 高い角度ピークを得ます(qc-grid::tuning、元素別グリッド調整、#390)— 上記のプリセットは一般の レシピであり、調整テーブルはコンパクトな内殻密度が一般的な点数では正確に積分しづらい元素についてそれを 精密化します。

検証済み例 — 水/def2-SVP、B3LYP、ラダー全体:

import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
for g in ("coarse", "medium", "fine"):
    m = qc.chk.new(atom=water, ao="def2-svp", unit="angstrom").scf(ref="r", xc="b3lyp", grid=g).run()
    print(g, m.scf.energy, m.scf.converged)
# coarse -76.35812209974867 True
# medium -76.35812511536525 True
# fine   -76.35812496750043 True

ここでの coarse→medium のエネルギーシフトは \(\approx3\times10^{-6}\ E_h\)(3 原子分)で、上記の 校正済み \(\lesssim8\to\lesssim0.4\ \mu E_h\)/原子の許容誤差と整合しています。medium→fine は 9 桁目 を変えるだけで、化学的精度をはるかに超えています;medium が qc-rs の既定なのはまさにこの理由です。

グリッド上の AO 値と密度#

各 AO は(GGA/meta-GGA では勾配も)保持されるすべてのグリッド点で評価される必要があります。原子 \(A\) を 中心とする縮約 Cartesian Gauss 関数は

\[ \chi_{\kappa,lmn}(\mathbf r) = \sum_p c_{\kappa p}\, x_A^l y_A^m z_A^n\, e^{-\alpha_p r_A^2}, \qquad \mathbf r_A = \mathbf r - \mathbf R_A, \]

勾配は通常の積の法則で得られます。例えば \(\partial_x(x^ly^mz^ne^{-\alpha r^2}) = (lx^{l-1}y^mz^n - 2\alpha x^{l+1}y^mz^n)e^{-\alpha r^2}\); 球面 AO は、qc-rs の他所と同じ solid-harmonic 変換 \(\chi^{\text{sph}}_s=\sum_cC^{(l)}_{sc}\chi^{\text{cart}}_c\) を(値にも各勾配成分にも)適用します。 制限・非制限密度とその勾配は、AO 値と密度行列から直接得られます:

\[ \rho(\mathbf r) = \sum_{\mu\nu}D_{\mu\nu}\chi_\mu(\mathbf r)\chi_\nu(\mathbf r), \qquad \nabla\rho^\sigma(\mathbf r) = \sum_{\mu\nu}D^\sigma_{\mu\nu}\bigl[(\nabla\chi_\mu)\chi_\nu + \chi_\mu(\nabla\chi_\nu)\bigr], \]

libxc の GGA 縮約勾配不変量は \(\sigma_{\alpha\alpha}=\nabla\rho^\alpha\cdot\nabla\rho^\alpha\)、 \(\sigma_{\alpha\beta}=\nabla\rho^\alpha\cdot\nabla\rho^\beta\)、 \(\sigma_{\beta\beta}=\nabla\rho^\beta\cdot\nabla\rho^\beta\) です(制限 GGA は \(\sigma=\nabla\rho\cdot\nabla\rho\) だけが必要)。重要な回帰テストは、密度を重なり整合の電子数に対して 積分し直します:\(\sum_gw_g\rho_g\approx\operatorname{Tr}[DS]\) — ここでの不一致は、XC 汎関数とは無関係 に、グリッドか AO 評価器のバグを意味します。

AO 評価は分子全体にわたって密に行われることは決してありません:各グリッドブロックはアクティブ AO — ブロックのどこかで大きさがカットオフ \(\tau_\phi\) を超える AO、\(\max_g|\chi_\mu(\mathbf r_g)|>\tau_\phi\) — をスクリーニングし、そのアクティブな部分集合だけが評価され、密度へ縮約され、後で \(V_{\text{xc}}\) へ縮約し戻されます。空間的に広がった分子では、これによりブロックあたりの AO バッファは 系全体のサイズによらず小さく保たれます — 2 電子積分に使われる Schwarz 殻カルテットスクリーニングと 同じ局所性の論法です。

\(E_{\text{xc}}\) と \(V_{\text{xc}}\) の組み立て#

libxc の粒子あたりエネルギー密度 \(\epsilon_{\text{xc}}\) とその汎関数微分(点ごとに評価、 SCF の章 がラッパーを扱う)が与えられれば、XC エネルギーは単純な重み付き和です:

\[ E_{\text{xc}} = \int\rho(\mathbf r)\,\epsilon_{\text{xc}}(\mathbf r)\,d\mathbf r \approx \sum_g w_g\,\rho_g\,\epsilon_{\text{xc},g}. \]

Fock 寄与は、離散化されたエネルギーを密度行列で微分することで得られます。制限 LDA では、

\[ V^{\text{xc}}_{\mu\nu} = \sum_g w_g\,v_{\rho,g}\,\chi_{\mu g}\chi_{\nu g}, \]

制限 GGA は積勾配恒等式 \(\nabla(\chi_\mu\chi_\nu)=(\nabla\chi_\mu)\chi_\nu+\chi_\mu(\nabla\chi_\nu)\) を 使って勾配項を加えます:

\[ V^{\text{xc}}_{\mu\nu} = \sum_g w_g\Bigl[v_{\rho,g}\chi_{\mu g}\chi_{\nu g} + 2v_{\sigma,g}\,\nabla\rho_g\cdot\nabla(\chi_{\mu g}\chi_{\nu g})\Bigr]. \]

非制限 GGA の形は、交差勾配不変量 \(\sigma_{\alpha\beta}\) を通じて両スピンチャネルを結合します:

\[ V^{\alpha}_{\mu\nu} = \sum_g w_g\Bigl[v_{\rho^\alpha,g}\chi_{\mu g}\chi_{\nu g} + \bigl(2v_{\sigma_{\alpha\alpha},g}\nabla\rho^\alpha_g + v_{\sigma_{\alpha\beta},g}\nabla\rho^\beta_g\bigr) \cdot\nabla(\chi_{\mu g}\chi_{\nu g})\Bigr], \]

\(\beta\) は全体で \(\alpha\leftrightarrow\beta\) を交換して得ます。これは、SCF の章 が説明する通りに、Coulomb/交換項と並んで組み立て可能な Fock ビルダーにはまり込みます:純粋な RKS は \(F=H+J[D]+V_{\text{xc}}[D]\)、大域ハイブリッドは \(-a_xK[D^\sigma]\) を加え、範囲分離ハイブリッドは さらに \(K\) を短距離/長距離に分けます — 本章のグリッドと libxc の機構は、どの参照/ハイブリッドの 組み合わせであっても、そのFock 式に \(V_{\text{xc}}\) と \(E_{\text{xc}}\) をまさに供給するものです。

        flowchart TD
    SPEC["grid= プリセット<br/>(coarse/medium/fine/ultrafine)"] --> RAD["Treutler-Ahlrichs<br/>原子ごとの動径殻"]
    RAD --> PRUNE["NWChem 5 領域<br/>角度プルーニング"]
    PRUNE --> PART["Becke 空間分割<br/>W_A(r), sum_A W_A = 1"]
    PART --> BLOCK["グリッドブロック<br/>(workspace レーンサイズ)"]
    BLOCK --> AOSCREEN["AO スクリーニング<br/>+ AO/勾配評価"]
    AOSCREEN --> DENS["rho, grad-rho<br/>(+ meta-GGA の tau)"]
    DENS --> LIBXC["libxc: eps_xc, v_rho, v_sigma, ..."]
    LIBXC --> ACC["E_xc + V_xc を<br/>FockAccumulator へ加算"]

性能モデル:Workspace・スクリーニング・スレッド化#

このパイプラインのすべての量はグリッド点数にスケールし、それは原子数に線形にスケールします — 2 電子 積分テンソルとは異なり、グリッドの作業量は組み合わせ的には爆発しませんが、素朴に扱えば中規模 DFT ジョブの実時間を支配します。qc-rs の性能規律は、コードベースの他所と同じ Workspace の規則 を反映します:

  • ホットループでの確保なし。 グリッドエンジンは、デバッグ/エクスポートモード (store_points=True)以外では自分の点/重みを所有しません;本番コードは run スコープの Workspace を一度だけ埋め、MolecularGridView/GridBlock スライスを再借用するので、 ブロックを反復しても何も確保しません。

  • ブロックローカルな AO スクリーニングは、ブロックあたりのアクティブ AO バッファを、分子全体の AO 数ではなく、そのブロック局所の AO 密度に比例させます — ERI に対する Schwarz スクリーニングを 効果的にするのと同じ局所性の論法です。ブロックサイズ \(B\)、アクティブ AO 数 \(N_{\text{AO,active}}\) に対して、LDA ワーカーレーンはおおよそ \(M_{\text{lane}}^{\text{LDA}}\sim8B(N_{\text{AO,active}}+c_{\text{LDA}})\) バイトを、AO 勾配も必要な GGA レーンはおおよそ \(M_{\text{lane}}^{\text{GGA}}\sim8B(4N_{\text{AO,active}}+c_{\text{GGA}})\) バイトを必要とします;workspace プランナーは、各ワーカーレーンと簡約バッファがプロセスのメモリ 予算に収まる最大のブロックサイズ \(B\) を選びます。

  • スレッド化:グリッドブロックループは、今日 SCF サイクルの中で真にスレッド並列な唯一の部分 (LDA・GGA・meta-GGA — \(\tau\) と Laplacian を含む — の RKS・UKS とも)です — 各ワーカースレッドは 互いに素なブロックとプライベート/タイル化された部分 \(V_{\text{xc}}\) を持ち、ビット再現性のため 最後に固定のスレッド ID 順で簡約されます。ブロックごとの縮約自体は単一スレッド GEMM (qc_array::blas::gemm_st)なので、外側のスレッド並列グリッドループと内側のマルチスレッド BLAS 呼び出しとの間に調停すべき入れ子のオーバーサブスクリプションはありません — スレッド予算は 1 つで、 BlasThreadGuard は不要です(仮に各ワーカーの内側でマルチスレッド BLAS を走らせる設計だったら 必要だったでしょう)。

2 つのエネルギー式:RKS 対 UKS#

\(E_{\text{xc}}\) と \(V_{\text{xc}}^\sigma\) が分かれば、純粋な(非ハイブリッド)KS-DFT の電子エネルギーは HF と同じトレース構造に従い、厳密交換を XC 汎関数で置き換えます:

\[ \text{RKS:}\quad E_{\text{elec}} = \operatorname{Tr}[DH] + \tfrac12\operatorname{Tr}[DJ] + E_{\text{xc}}, \qquad F = H + J[D] + V_{\text{xc}}[D], \]
\[ \text{UKS:}\quad E_{\text{elec}} = \operatorname{Tr}\bigl[(D^\alpha{+}D^\beta)H\bigr] + \tfrac12\operatorname{Tr}\bigl[(D^\alpha{+}D^\beta)J\bigr] + E_{\text{xc}}, \qquad F^\sigma = H + J[D^\alpha{+}D^\beta] + V_{\text{xc}}^\sigma. \]

純粋汎関数ではどちらの式にも交換項がまったくありません — これは単に異なる数値係数ではなく Hartree-Fock との根本的な構造の違いであり、純粋な GGA/LDA 実行が 2 電子交換縮約のコストを一切 払わない理由であり、ハイブリッド(xc="b3lyp"、xc="pbe0")がこのグリッド機構と SCF の章 の厳密交換 \(K\) 構築の両方を必要とする理由です。

練習 7

  1. 水/def2-SVP/B3LYP は coarse・medium・fine でそれぞれエネルギー \(-76.358122100\)、 \(-76.358125115\)、\(-76.358124968\ E_h\) を与えます。校正済みの最悪許容誤差(\(\lesssim8\)、 \(\lesssim0.4\)、\(\lesssim0.2\ \mu E_h\)/原子、3 原子)を使って、観測された coarse\(\to\)medium シフト(\(\approx3\times10^{-6}\ E_h\))はこれらの許容誤差と整合していますか。また、はるかに小さい medium\(\to\)fine シフトは medium が精度ラダー上のどこに位置するかについて何を示しますか。

  2. Becke 分割は Bragg-Slater 半径を使って \(\mu_{AB}\) を \(\nu_{AB}\) に補正します。水素に結合した 重原子について、生の \(\mu_{AB}\)(素朴な垂直二等分境界)を使うのがなぜ悪い選択なのか一文で 説明しなさい。

  3. ゴースト原子は基底関数を持ちますが決して分割中心にはなりません。もしゴースト原子に Bragg-Slater 半径を与えて中心として扱ったら、Becke 分割の正規化 \(\sum_AW_A=1\) に何が問題を起こすでしょうか。

ここでのブロックローカルな AO/密度機構は、線形応答理論(基底状態密度 ではなく遷移密度を評価する)と、解析的微分(同じ \(s_h(\nu_{AB})\) 切り替え関数から導かれる Becke 重み微分 \(\partial W_A/\partial R_C\) が DFT 核勾配の一部になる)に そのまま再登場します。