解析的微分:核勾配

解析的微分:核勾配#

あらゆる幾何最適化・振動解析・第一原理分子動力学ステップは、各原子核にはたらく力 \(-\partial E/\partial X_A\) を必要とします。これを有限差分で計算すると、追加で \(6N_{\text{atom}}\) 回の SCF エネルギー(各 Cartesian 座標を前後に摂動する)がかかり数値的に ノイズが多くなります;解析的に計算すると、エネルギー構築が既に使ったのと同じ微分積分クラスへの おおよそもう 1 回分のパスだけで済み、機械精度が得られます。本章は完全な解析的 RHF/UHF/ROHF 勾配を 導出します — なぜ結合摂動解法が不要なのか、ROHF の正しさに関する機微を含むあらゆる SCF 参照に対する 厳密なエネルギー重み付き密度、そして 2 電子勾配のカルテット/RI 構造 — .design/74analytic-nuclear-gradients.md に基づきます。勾配/幾何最適化の章 がコンパクトな使い方ダイジェスト を持ちます;本章は導出を完全に持ちます。

全微分と鍵となる単純化#

変分的 SCF エネルギー \(E(\{C\},\{\varepsilon\};X)\) について、原子核 Cartesian 座標 \(X_A\) に関する 全微分は、原理的には 3 種類の項を持ちます:固定 AO 基底・密度でハミルトニアンの演算子が \(X_A\) に どう依存するか、AO 基底関数自体がどう動くか(各 GTO は原子核を中心とするので)、そして収束済み MO 係数が幾何摂動にどう応答するか。最後のもの — 軌道応答 — が一般に勾配を高価にするもの(結合摂動 Hartree-Fock/Kohn-Sham、次章で完全に導出)です。しかしエネルギー勾配 に限っては、変分原理の直接的帰結として厳密にゼロになります:

\[ \frac{\partial E}{\partial X_A}\bigg|_{\text{via }C} = \sum_{\mu\nu}\frac{\partial E}{\partial D_{\mu\nu}}\cdot\underbrace{\frac{\partial D_{\mu\nu}}{\partial C}\frac{\partial C}{\partial X_A}}_{\text{軌道応答}} = 0, \]

SCF 収束で \(\partial E/\partial C=0\) だからです — これが変分エネルギー汎関数に適用された Hellmann-Feynman 定理そのものです。落とし穴は、HF/KS 軌道が単に固定演算子の固有関数であるだけでなく、それ自体 幾何(AO 重なり \(S\) を通じて)に依存する直交規格性拘束 \(C^{\mathsf T}SC=I\) を満たすことです。その 拘束を微分することが、軌道回転応答は完全に消えるにもかかわらず、幾何依存の項 — Pulay 力 — を再び持ち込むものです。生き残るすべての項を集めると(RHF 形式;UHF/ROHF/KS は同じ骨格の per-spin/ 有効一般化):

\[ \frac{\partial E}{\partial X_A} = \underbrace{\sum_{\mu\nu}D_{\mu\nu}\frac{\partial h_{\mu\nu}}{\partial X_A}}_{\text{コア Hamiltonian}} + \underbrace{\tfrac12\sum_{\mu\nu\lambda\sigma}D_{\mu\nu}D_{\lambda\sigma} \frac{\partial(\mu\nu|\lambda\sigma)}{\partial X_A}\bigg|_{JK}}_{\text{2 電子 }(J-c_xK)} - \underbrace{\sum_{\mu\nu}W_{\mu\nu}\frac{\partial S_{\mu\nu}}{\partial X_A}}_{\text{Pulay}} + \frac{\partial E_{\text{xc}}}{\partial X_A} + \frac{\partial E_{\text{ECP}}}{\partial X_A} + \frac{\partial E_{\text{PCM}}}{\partial X_A} + \frac{\partial V_{nn}}{\partial X_A}. \]

実務上の帰結:収束した HF あるいは純粋/ハイブリッド KS 勾配は、収束済み密度 \(D\) ともう 1 つの対象 — エネルギー重み付き密度 \(W\) — だけを必要とします — 反復的な結合摂動解法はまったく不要です。 (Z ベクトル/CPHF 機構が必要になるのは変分的でない量の勾配 — MP2 やその他のポスト SCF 相関 エネルギー — だけであり、これは本章が扱う素の SCF 勾配の範囲外であり、代わりに 線形応答理論 が本質的になる場面です。)

エネルギー重み付き(Lagrangian)密度 \(W\)#

\(W\) はまさに、Pulay 項が軌道直交規格性拘束の幾何依存性を吸収するために必要とする対象です — それは 重なり項において、Hamiltonian 項での \(D\) が果たす役割を果たします。

RHF:\(W = C_{\text{occ}}^{\mathsf T}\operatorname{diag}(2\varepsilon_{\text{occ}})\,C_{\text{occ}}\) — 占有 MO から \(D\) を構築するのと構造的に同一ですが、各軌道を占有数ではなく軌道エネルギーで 重み付けます。

UHF:\(W = \sum_\sigma C_{\text{occ}}^{\sigma\mathsf T}\operatorname{diag}(\varepsilon^\sigma_{\text{occ}})\, C_{\text{occ}}^\sigma\) — 同じ構成を、各スピンチャネルで独立に行います。

ROHF/ROKS — 正しさの地雷。 ROHF は両スピンで共有される 1 つの有効 Roothaan Fock 演算子を対角 化します(これはROHF-MP2 を曖昧にしたのと同じ事実です)ので、 orbital_energies として保存される固有値は、その有効演算子の固有値であって真のスピン分解された 拘束に対する真の Lagrange 乗数ではありません — この固有値で RHF/UHF の式を素朴に再利用すると、 もっともらしく見えるが間違った勾配を与えます(開殻がゼロという自明な場合にのみ正しい)。正しい 構成は代わりに、収束したスピン Fock 行列 \(F_\alpha,F_\beta\)(収束済み ROHF 密度から既に得られる) と per-spin AO 射影子から直接 \(W\) を再構築します:

\[ W = W_\alpha + W_\beta, \qquad W_\sigma = P_\sigma\,F_\sigma\,P_\sigma, \qquad P_\alpha = C_{\text{occ}}^{\alpha\mathsf T}C_{\text{occ}}^\alpha,\quad P_\beta = C_{\text{occ}}^{\beta\mathsf T}C_{\text{occ}}^\beta, \]

\(P_\alpha\) は \(\alpha\) 占有空間(単一占有かつ二重占有軌道)への射影、\(P_\beta\) は \(\beta\) 占有 空間(二重占有のみ)への射影です — これが正準的構成です(PySCF の grad/rohf.py:make_rdm1e と 一致;以前の設計草案が推測して後で修正しなければならなかった、より単純な MO 基底ブロック則では ありません)。qc-rs 自身の設計ノートによれば、これは「勾配全体の中で最も微妙な正しさのポイント」 です — ここでの誤りは、妥当に見える(桁数も定性的方向も正しい)にもかかわらず数値的に正しくない 勾配を生むので、信頼される前に有限差分比較で検証されます。

1 電子(コア Hamiltonian)と Pulay 勾配#

\(h=T+V_{ne}\) なので、\(\partial h/\partial X = \partial T/\partial X + \partial V_{ne}/\partial X\) です。 運動エネルギー項は純粋に AO 基底中心の微分です(運動エネルギー演算子自体は核座標を持たない)。核 引力項は物理的に異なる 2 つの部分に分かれます:AO 基底中心部分(bra/ket Gauss 関数がどう動く か)と、真にHellmann-Feynmanな部分 — 演算子 \(1/|\mathbf r-\mathbf R_B|\) 自体が核 \(B\) の位置に 依存するので、それを \(R_B\) について微分すると、その核の有効電荷 \(Z_{\text{eff},B}\) でスケールされた、 核そのものを中心とする項が得られます:

\[ \sum_{\mu\nu}D_{\mu\nu}\frac{\partial h_{\mu\nu}}{\partial X_A} = \underbrace{\sum_{\mu\nu\in A}D_{\mu\nu} \left(\frac{\partial T_{\mu\nu}}{\partial X_A}+\frac{\partial V_{ne,\mu\nu}^{\text{basis}}}{\partial X_A}\right)}_{\text{AO 中心、原子 }A\text{ 上の殻へ散布}} + \underbrace{Z_{\text{eff},A}\sum_{\mu\nu}D_{\mu\nu}\frac{\partial}{\partial R_A}\frac{1}{|\mathbf r-\mathbf R_A|}\bigg|_{\mu\nu}}_{\text{Hellmann-Feynman、核ごとに 1 項}}. \]

Pulay 項は \(W\) が手元にあれば構造的に単純です — 単に \(-\sum_{\mu\nu}W_{\mu\nu}\, \partial S_{\mu\nu}/\partial X_A\) であり、Hamiltonian 項が \(D\) と縮約するのとまったく同じ方法で \(W\) と縮約される AO 中心の重なり微分です。ここでのすべての AO 中心寄与は自然に殻対にわたるループと して構成され、各対の寄与を、微分される殻を所有する原子へ散布します;Hellmann-Feynman 部分は代わり に原子核にわたるループです、演算子自体が基底関数ではなく核の上にあるからです。古典的な核反発勾配 \(\partial V_{nn}/\partial X_A\)(点電荷間の初等的な Coulomb の法則の微分)を加えれば 1 電子側は 完成します。

2 電子勾配:なぜ Hellmann-Feynman 項がないのか#

化学者の記法では、\((\mu\nu|\lambda\sigma)=\iint\chi_\mu(\mathbf r_1)\chi_\nu(\mathbf r_1)\,r_{12}^{-1}\,\chi_\lambda(\mathbf r_2)\chi_\sigma(\mathbf r_2)\,d\mathbf r_1\,d\mathbf r_2\) です。 演算子 \(r_{12}^{-1}\) は核座標をまったく持ちません — \(V_{ne}\) の \(1/|\mathbf r-\mathbf R_B|\) とは 異なり、この演算子が乗る核が存在しません。したがって 2 電子勾配は、\(\mu,\nu,\lambda,\sigma\) の 4 つの AO 中心を通じてだけ幾何に依存し、Hellmann-Feynman 項をまったく持ちません — すべての寄与が 基底中心(Pulay 型)の微分です。これが、2 電子勾配が基底中心の微分積分だけを必要とし、裸の核を中心 とする積分を決して必要としない構造的理由です。

厳密な 4 中心経路では、微分は最初の bra 殻上で取られ(libcint の int2e_ip1)、GTO が \(\chi(\mathbf r-\mathbf A)\) であることから \(\partial\chi/\partial\mathbf A=-\nabla\chi\) です — 符号 反転を伴う電子座標勾配です。カルテットの完全な微分は 4 つの中心すべてを順番に微分する必要があります が、置換対称性と並進対称性が残りのほとんどを無料で回復します:カルテットのどの殻が「最初の」 (微分される)役割を演じるかを回転させることで残りの中心をカバーし、分子全体の恒等式 \(\sum_A\mathbf F_A=0\)(並進不変性 — 自由分子にはたらく正味の力はゼロになる)は最後の中心の寄与を 直接計算せずに他から推論することを可能にします。エネルギー側の Coulomb/交換パターンは変更なく勾配 縮約に引き継がれます — \(J\) は両方の添字対に対して全密度と縮約し、\(K\) は片方の添字を交差させてスピン 密度と縮約します — エネルギー側 Fock 構築と同じ \(\tfrac12\)、\(a_x\)(厳密交換分率)、範囲分離重み スカラーが組み立て時に適用されます。

RI 2 電子勾配:追加の必須部分#

RI 因子化された 2 電子エネルギーは \((\mu\nu|\lambda\sigma)_{\text{RI}}=\sum_{PQ}(\mu\nu|P)\, V^{-1}_{PQ}\,(Q|\lambda\sigma)\) です(RI の章で完全に導出)。この因子化 された形を微分するには、3 中心積分 \((\mu\nu|P)\) と2 中心計量 \(V_{PQ}=(P|Q)\) の両方を微分 する必要があります — そして \(V^{-1}\) 自体が幾何に依存するので、その微分は連鎖律を通じて考慮され なければなりません。フィッティングされた電荷 \(\gamma_P=\sum_QV^{-1}_{PQ}c_Q\)(\(c_Q=\sum_{\mu\nu} (Q|\mu\nu)D_{\mu\nu}\)、RI の章の常駐白色化 \(B\) 因子の \(q_P\) そのもの、新しい 対象は不要)を定義すると、Coulomb エネルギーは \(E_J=\tfrac12\gamma^{\mathsf T}V\gamma\)(\(c=V\gamma\) を使用)であり、\(V^{-1}\) 連鎖律恒等式 \(\partial(K^{-1}R)=K^{-1}(\partial R - \partial K\cdot q)\) (PCM の章も自身の応答微分にこの恒等式を使う)を通じてそれを微分すると、 明示的な \(\partial V^{-1}/\partial X\) を残さないきれいな 2 項の結果が得られます:

\[ \frac{\partial E_J}{\partial X} = \sum_P\gamma_P\sum_{\mu\nu}\frac{\partial(P|\mu\nu)}{\partial X}D_{\mu\nu} \;-\; \tfrac12\sum_{PQ}\gamma_P\gamma_Q\,\frac{\partial(P|Q)}{\partial X}. \]

類似の交換勾配式はもう 1 つのフィッティングされた対象、2 重に MO 変換された計量因子 \(G_{PQ}=\sum_{ij}\tilde d_{P,i,j}\tilde d_{Q,j,i}\)(意図的に \(i,j\) について対称ではない — 一方の 添字は占有重み付けを持ち、他方は持ってはならない)を必要としますが、構造は同じです: \(\partial(P|\mu\nu)/\partial X\) からの AO 中心項と \(\partial(P|Q)/\partial X\) からの計量項です。

3 つの異なる微分積分クラスとその原子散布先:

積分

微分対象

散布先

int3c2e_ip1

\((P|\mu\nu)\)、AO 中心

原子(\(\mu\))、原子(\(\nu\))

int3c2e_ip2

\((P|\mu\nu)\)、補助中心

原子(\(P\))

int2c2e_ip1

\((P|Q)\)、補助中心

原子(\(P\))、原子(\(Q\))

補助中心散布は、RI 勾配が持ち 4 中心勾配が持たない、真に新しい構造的部分です:補助フィッティング 関数自体が原子中心の Gauss 関数なので、真の物理的な力を運びます — そして決定的に、この補助基底 応答は、一部のコードがトグル可能な近似として扱うのとは異なり、qc-rs にとってオプションではありません。 qc-rs の既定 JK フィッティング基底(default_jk_aux)が原子中心である以上、補助中心項を省略すると 力の実在部分が落ち、有限差分検査に完全に失敗します;仮想的な幾何非依存の補助基底だけがそれらを省略 できるでしょう。

4 中心の場合との構造的な違いは直接指摘する価値があります:4 中心では、\(\sum_A\mathbf F_A=0\) は 内部一貫性(置換/転置の帳簿つけ)だけの検査です、あらゆる寄与が構成上既にバランスしているからです。 RI では、同じ和ゼロの恒等式が真の欠落した補助応答項の検出器です — \((P|\mu\nu)\) の補助中心微分 や計量微分を落とすと、並進不変性が微妙にではなく測定可能なほど破れます。そして RI 自体が近似で ある以上、RI 勾配は厳密な 4 中心勾配と RI フィッティング誤差(RI 対 4 中心エネルギーで見られる のと同じ \(10^{-5}\) オーダーのずれ)までしか一致すべきではありません — 浮動小数点精度までではあり ません。

検証済み例 — 水/cc-pVDZ、厳密な 4 中心勾配を RI 勾配と比較し、両方の並進不変性を確認:

import qc, numpy as np
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").scf(ref="r").run()
g4c = np.array(m4c.scf.gradient)
g4c.sum(axis=0)   # [-1.4e-29, 2.2e-15, -8.0e-15] -- 並進不変性は ~1e-15 まで

mri = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").ints(eri="ri-ram").scf(ref="r").run()
gri = np.array(mri.scf.gradient)
gri.sum(axis=0)   # [1.5e-27, 1.5e-15, -4.4e-16] -- こちらも数値ノイズまでゼロ

np.abs(g4c - gri).max()   # 2.05e-05 -- エネルギーで見られる RI フィッティング誤差スケールと一致し、
                          # 浮動小数点レベルの一致ではない

両方の勾配は数値ノイズ(\(\sim10^{-15}\))まで並進不変であり、4c/RI の差(\(\sim2\times10^{-5}\))は 浮動小数点ノイズの底ではなくまさに RI 近似誤差のスケールに位置します — 設計ノートが予測する通りの 2 つの異なる診断挙動です。

解析的勾配が必要としないもの(そして何が必要とするか)#

本章の導出全体は、エネルギーが軌道の停留汎関数であるという事実にかかっています — それが Hellmann-Feynman を通じて軌道応答項をゼロに崩壊させたものです。2 つの重要な帰結が続き、どちらも 本章の範囲の境界を示すので明示的に述べる価値があります:

  • DFT XC・ECP・PCM の各寄与はそれぞれ自身の勾配項を加えます(\(\partial E_{\text{xc}}/\partial X\) など)が、どれも軌道応答を再導入しません — これらはすべて同じ停留 \(E\) への変分的寄与のまま であり、同じ Hellmann-Feynman の崩壊がそれぞれに適用されます。XC 項はグリッド重み自体の幾何微分 (DFT グリッドの章で導出された Becke 分割重み微分 \(\partial W_A/\partial R_C\))に加えて AO-原子追従微分 \(\partial\chi/\partial X=-\nabla\chi\) を必要とします;ECP 項は 自身の Pulay 部分を持たずに \(\partial\langle\chi|U_{\text{ECP}}|\chi\rangle/\partial X\) を必要と します(ECP は真に 2 中心の演算子行列であり、各中心が直接微分されます、直交規格性拘束ではあり ません);PCM は同じテッセラポテンシャル積分微分に加えて空洞/テッセラ幾何と応答の微分を必要と します。

  • 変分的でない量の勾配はまったく別の問題です。 MP2(および他のポスト SCF 相関エネルギー)は 参照軌道に関して停留ではありません — 参照行列式は相関エネルギーではなくSCFエネルギーに 対して最適化されたものだからです — なのでその勾配は本当に軌道応答を必要とします、Z ベクトル/ 結合摂動方程式を通じて。これは本章が導出する SCF 勾配の範囲から明示的に外れており、代わりに 次章の主題です。

練習 10

  1. 検証済み例は、4c/RI 勾配の差(\(\sim2\times10^{-5}\))が、どちらの勾配自身の並進不変性からのずれ (\(\sim10^{-15}\))よりおよそ 3 桁大きいことを示しています。これら 2 つの数値がまったく異なるもの を測っている理由、そしてなぜどちらも他方を縮小させると期待すべきでないのかを説明しなさい。

  2. ROHF の orbital_energies は有効 Roothaan Fock 演算子の固有値であり、それを RHF/UHF 式の \(W\) に直接使うと、もっともらしく見えるが間違った勾配を与えます。有限差分検査なしにはこの誤りが 見逃されやすい理由(なぜ間違った勾配が依然「妥当に見える」のか)を一文で説明しなさい。

  3. 2 電子勾配には Hellmann-Feynman 項がありませんが、1 電子核引力勾配にはあります。どちらの演算子 も形式的には 2 つの AO 中心の間で作用します。これを説明する \(r_{12}^{-1}\) と \(1/|\mathbf r-\mathbf R_B|\) の間の 1 つの構造的違いは何ですか。

ここで導出したエネルギー重み付き密度 \(W\) と Pulay 力の構造は、解析的 Hessian の章 で一般化されて再登場します — 2 次核微分は軌道応答を実際に必要とします(次に導出する 結合摂動方程式)、まさに停留量の 2 次微分はそれ自体もはや停留ではないからです。