線形応答:CPHF/CPKS

線形応答:CPHF/CPKS#

一見無関係に見える驚くほど多くの問いかけ — 「この SCF 解は本当に極小なのか、それとも鞍点なのか」 「停留点近傍でエネルギー曲線は 2 次でどう振る舞うか(Hessian)」「どの Newton ステップが SCF を最も 速く収束させるか」「この分子はどの光子エネルギーで光を吸収するか」— はすべて同じ線形代数の対象 に帰着します:試行軌道回転への電子(軌道)Hessian の作用、\(\sigma=H\kappa\) です。qc-rs はこの作用を 一度だけ、単一の注入可能な FockResponseKernel として構築し、上記のあらゆる消費者はそれに適用 される少数の汎用ソルバのどれかの薄いラッパになります。本章はこの統一対象とその 4 つの消費者を導出 します。.design/51scf-linear-response.md に基づきます。

統一対象:\(A\) と \(B\)#

SCF 収束の章 と同じ方法で試行軌道回転をパラメータ化します、 \(\mathbf C(\kappa)=\mathbf C\exp(\kappa)\)、非冗長な占有-仮想ブロックに限定して。エネルギーの \(\kappa\) での 2 次展開はそこで軌道 Hessian \(H_{ai,bj}\) を定義しました;線形応答理論はまさにその同じ 対象(とその近縁)に標準的な名前を与え、それが SCF 収束をはるかに超えるものを支配していることを 示すだけです。閉殻(RHF/RKS)singlet 回転について、実正準 MO で占有 \(i,j\)・仮想 \(a,b\) とすると:

\[ A_{ia,jb} = \delta_{ij}\delta_{ab}(\varepsilon_a-\varepsilon_i) + 2(ia|jb) - (ij|ab)\ \bigl[+\,2f^{xc}_{ia,jb}\bigr], \qquad B_{ia,jb} = 2(ia|bj) - (ib|aj)\ \bigl[+\,2f^{xc}_{ia,bj}\bigr], \]

角括弧の \(f^{xc}\) 項(交換相関核、\(E_{\text{xc}}\) の第 2 汎関数微分)は KS 参照では存在し純粋 HF では存在しません。\(A,B\) の 2 つの組み合わせがいたるところで再登場します:\(A+B\) はまさに SCF 収束の章の実軌道回転 Hessianそのものです — そこで拡大 Hessian/QC-SCF と TRAH 法を駆動するのと 同じ対象です — 一方 \(A-B\) は別の物理的問い(三重項/スピンフリップ応答、Coulomb 様の項を落とし交換 だけを残す)を支配します。ある問題がどちらの組み合わせを必要とするか、そして固有値問題か線形解法か は以下にまとめます:

消費者

必要な演算子

解く問題

内部安定性

\(A+B\)

最低固有値 \(<0\) ⇒ 不安定(実軌道回転方向)

外部/三重項安定性

\(A-B\)(または三重項カーネル)

最低固有値 \(<0\) ⇒ スピンフリップ/対称性破れ不安定性

QC-SCF/拡大 Hessian/TRAH

\(A+B\)

拡大(境界付き)行列の最低固有対 — 既に導出済み

CPHF/CPKS(核/場応答)

\(A+B\)

摂動固有の右辺 \(b\) に対する線形解法 \((A+B)U=-b\)

TDA(Tamm-Dancoff 励起状態)

\(A\) のみ

最低数固有対、\(AX=\omega X\)

フル TDDFT/RPA

\(A,B\) 結合

\((A-B)(A+B)Z=\omega^2Z\)(Hermite 形)

どの行も同じ背後の \(A\)/\(B\) 作用を必要とします — 組み合わせとソルバだけが異なります。これが構造上の 要点のすべてです:\(\sigma=(A{+}B)\kappa\)(または \(A\kappa\)、あるいは \((A{-}B)\kappa\))を一度だけ、 単一の再利用可能なプリミティブとして構築すれば、すべての消費者は薄いラッパになります。

\(A\) や \(B\) を一切組み立てずに \(\sigma\) を構築する#

\(A\) や \(B\) を明示的に \((n_{\text{occ}}n_{\text{vir}})\times(n_{\text{occ}}n_{\text{vir}})\) の密行列 として構築することは、まさに qc-hpc Workspace 規則 が材料化を禁じる \(O(N^4)\) スケーリングの対象です。代わりに \(\sigma\) 構築は、\(A+B\) が試行ベクトル \(X\) に及ぼす作用を 単一の応答 Fock 構築 — 単一添字の遷移密度に対する Fock 様の縮約 — として評価し、通常の SCF Fock 構築が既に使うのとまったく同じ \(J/K\)(そして KS では同じ XC グリッド)機構を再利用します:

アルゴリズム 4 (1 回の応答 Fock 構築による \(\sigma=(A{+}B)X\))

入力: 試行ベクトル \(X\)(形状 \(n_{\text{vir}}\times n_{\text{occ}}\))、占有/仮想 MO 係数 \(C_o,C_v\)、軌道エネルギー \(\varepsilon\)。 出力: \(\sigma=(A+B)X\)、\(X\) と同じ形状。

  1. AO 遷移密度 \(P = C_v\,X\,C_o^{\mathsf T}\) を作る(一般に対称である必要はない;\(A+B\) の組み合わせ に限っては対称化された \(P+P^{\mathsf T}\) を使う — これが \(A+B\) を 2 回ではなく単一の対称応答 Fock 構築に崩壊させるものである)。

  2. 注入された FockResponseKernel で応答 Fock \(G[P]\) を構築する:HF/ハイブリッドでは \(2J[P]-a_xK[P]\)、KS では \(+2f^{xc}[P]\) を加える(同じ DFT 求積グリッド で 評価されるが、libxc の 1 次ではなく 2 次汎関数微分を使う)。

  3. 占有-仮想ブロックへ変換し直す:\(\sigma_{\text{2e}} = C_v^{\mathsf T}G[P]\,C_o\)。

  4. 対角(軌道エネルギー差)部分を加える:\(\sigma = \sigma_{\text{2e}} + (\varepsilon_a-\varepsilon_i)X_{ia}\)。

1 回の \(\sigma\) 構築のコストはちょうど 1 回の Fock 構築です — 同じ \(J/K\) エンジン(incore/RI/ direct、SCF 自身が使ったものと同じ)、同じ Schwarz スクリーニング、KS では同じスレッド化 XC グリッド ループ。これが qc-rs の応答基盤における最も重要な設計上の決定です:\(J/K\) や XC グリッドループへの 将来のあらゆる高速化が、安定性解析・QC-SCF・CPHF・TDA/TDDFT に自動的に恩恵を与えます、別個の最適化 作業なしで。すべてが同一の Fock 構築機構を経由するからです — qc-rs は意図的に応答エンジンに独自の 積分経路を持たせません。複数の同時試行ベクトル(Davidson のブロック反復、CPHF の複数右辺 — 摂動ごとに 1 つ)は自然にバッチ化されます:複数の遷移密度を束ね、多数の別個の Fock 構築の代わりに 1 回のバッチ Fock 構築を行い、積分/グリッドアクセスをそれらすべてで償却します。

参照カバレッジ:1 つの抽象化、3 つの回転構造#

応答形式は 6 つの参照(RHF/UHF/ROHF × HF/KS)をカバーしなければなりませんが、きれいに 2 つの独立な 軸に分解します — スピン構造(非冗長な回転空間と \(\sigma\) の組み立て方を決める)とカーネル族 (\(f^{xc}\) を加えるかどうか)— なので qc-rs は 6 つの別個の経路ではなく、3 つの回転構造と 2 つの カーネル族の直積を実装します:

参照

非冗長回転空間

スピンチャネル

カーネル

RHF/RKS

占有-仮想(空間軌道)

1(singlet/triplet はカーネル選択で分離)

RestrictedSinglet/Triplet(RKS では +\(f^{xc}\))

UHF/UKS

(\(\alpha\) 占有-仮想) \(\oplus\) (\(\beta\) 占有-仮想)

2(\(\alpha\)/\(\beta\) は \(J\) と \(f^{xc}_{\sigma\sigma'}\) を通じて結合)

Unrestricted(UKS では +\(f^{xc}_{\sigma\sigma'}\))

ROHF/ROKS

(閉殻-開殻) \(\oplus\) (閉殻-仮想) \(\oplus\) (開殻-仮想)

1 セット、3 空間

制約付きスピン軌道 Hessian + ROHF 結合演算子

Kohn-Sham 化は 3 つすべてで構造的に同一です:各 HF 参照カーネルは単に \(f^{xc}\) 項を得るだけで (UKS では同スピン・交差スピンブロック)、背後の回転構造はまったく変わりません — RKS は RHF に \(f^{xc}\) を足したもの、UKS は UHF に \(f^{xc}_{\sigma\sigma'}\) を足したもの、ROKS は ROHF に \(f^{xc}\) を足したものです。ROHF/ROKS は真に最も難しいケースです(1 つや 2 つの分離した回転空間では なく 3 つの結合した回転空間)、半占有軌道が \(f^{xc}\) にどう入るかに特に注意が必要ですが、4 つ目の アーキテクチャパターンを必要としません — 同じ \(\sigma\) 構築抽象化の内側でより注意深い会計を必要 とするだけです。

外部ポテンシャル(ECP、PCM、点電荷/溶媒)は特殊ケースというより組織原理として組み込まれます: \(\sigma\) 構築は「各密度依存項の応答」の和であり、密度に依存しない外部ポテンシャル(固定 ECP や点 電荷)は \(\sigma\) にまったく何も追加で寄与しません — そのエネルギー寄与は、通常の SCF Fock 構築が 使った軌道エネルギーと収束済み密度に既に焼き込まれています;真に密度依存の環境項(どんな密度 (遷移密度を含む)にも応答する PCM 反応場のような)だけが自身の応答寄与を \(\sigma\) 構築に加える必要 があります。

4 つの汎用ソルバ、1 つの共有エンジン#

\(\sigma\) 構築が不透明なクロージャとして存在すれば、すべての消費者は少数の再利用可能なソルバ集合の どれかをそれに呼び出すことに帰着します — これらのソルバはどれも \(\sigma\) が RHF から来たのか UKS から来たのか、あるいは他の何かから来たのかを知りも気にもしません:

  • davidson_lowest — 対角前処理付きブロック Davidson 法による、大きな疎作用演算子の最低 \(k\) 固有対。安定性解析(\(k=1\)、\(A+B\) または \(A-B\) について)、TDA(\(k=\)欲しい励起状態数、\(A\) のみに ついて)、そして拡大 Hessian ソルバの内部で使われます。

  • augmented_hessian_step — 拡大 Hessian/QC-SCF 固有値問題 (境界付き行列の最低固有対、自動レベルシフトを与える)、davidson_lowest の上に構築されます。

  • linear_solve — 対称正定値系 \((A+B)X=\text{rhs}\) に対する前処理付き共役勾配法。これが CPHF/CPKS ソルバであり、摂動が核変位(解析的 Hessian に供給される) であれ、外部電場(分極率)であれ、他のどんな単一添字摂動であれ、同一に使われます。

  • rpa_solve — フル TDDFT/RPA(TDA の単純な \(A\) のみの固有値問題とは対照的に)のための結合 \((A-B)(A+B)Z=\omega^2Z\) 固有値問題の Hermite 化簡約。

内部安定性解析は davidson_lowest(apply_aplusb, diag_precond, k=1) です:\(A+B\) の最低固有値が負なら、 現在の SCF 解は軌道回転空間で真の極小ではなく鞍点であることを意味し、その固有ベクトルはより低い解へ 軌道を回転させる方向を与えます(SOSCF が使うのと同じ Cayley 回転機構)。 外部/三重項安定性は \(A-B\) または専用の三重項カーネルに対する同じ呼び出しです。

検証済み例 — 水/STO-3G 安定性解析、内部(実軌道)と外部(スピンフリップ)両チャネル:

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="sto-3g", unit="angstrom").scf(ref="r", stability=True).run()
m.scf.energy, m.scf.converged
# (-74.96294665653868, True)
m.scf.stability
# {'internal_stable': True, 'internal_eigenvalue': 0.5236790238088829,
#  'external_stable': True, 'external_eigenvalue': 0.36925930146991304, 'stable': True}

ここでは内部(\(A+B\))と外部(\(A-B\))の両方の最低固有値が正です — この制限閉殻解は、実軌道回転の 意味でもスピン対称性破れの意味でも、真の局所極小です。まさに本節の 2 つの固有値の行が答える 2 つの 問いです。

CPHF/CPKS:安定性/QC-SCF と何が違うか#

安定性解析と QC-SCF は \(A+B\) の固有値だけを必要とします — 「どの方向に、もしあれば、エネルギーが 下がるか」と問います。CPHF/CPKS は違う問いを発します:「特定の外部摂動(核変位、電場)が与えられた とき、軌道は 1 次でどう応答するか」。これは線形解法であり固有値問題ではありません、摂動への応答 \(U\) は(参照が安定である、すなわち \(A+B\) がゼロ固有値を持たないと仮定すれば)極値問題ではなく一意に 決まるからです:

\[ (A+B)_{ai,bj}\,U_{bj} = -b_{ai}, \]

\(b\) は摂動されている対象の微分積分から構築される摂動固有の右辺です。核座標摂動については、この \(b\) は 解析的勾配の章 が導入したのと同じ 1 次微分(勾配レベル)積分から 構築されます — これはまさにその章が自身の範囲外だとフラグを立てた「真に軌道応答を必要とする」場合 であり、次章、解析的 Hessian がこの右辺と結果の折り返しを完全に導出 します。外部電場については、\(b\) は代わりに AO 双極子積分から構築されます — 同じソルバ、同じ \(A+B\) 演算子、異なる摂動固有の右辺です;これが静的双極子分極率計算が最終的に帰着するものです。いずれに せよ、linear_solve は安定性解析と QC-SCF が既に使うのと同一の apply_aplusb クロージャに対して 呼ばれます — これが「CPHF は SCF 収束機構を再利用する」がゆるいアナロジーではなく文字通りのコード 共有の事実である意味です。

励起状態応答が今日どこに立っているか#

TDA とフル TDDFT/RPA は、この設計の中でアーキテクチャ的に考慮されています(それぞれ davidson_lowest(apply_a, ...) と rpa_solve)— 本章が導出する同じ \(\sigma\) 構築抽象化が、まさに 実際の実装が消費するはずのものです。しかし本稿執筆時点で qc.td(...) は qc-rs ではモックの ステップです:同じ保留ステップの構文を受け付け、プレースホルダーの converged=True 結果を記録します が、背後に実際の固有値解法はまだありません。これは ポスト HF 相関の章 の cc2/caspt2/nevpt2 に関する正直な記述を反映しています — アーキテクチャ上の継ぎ目は実在し重要な 役割を担っています(内部安定性解析は既に TDA/TDDFT が必要とするのとまったく同じ apply_a/ apply_aplusb 機構を行使しています)が、励起状態ソルバ自体はまだそこに配線されていません。

練習 11

  1. 検証済み安定性の例は、\(A+B\) について 1 つ、\(A-B\) について 1 つの2 つの固有値を報告します。 制限閉殻参照が真の極小であると確信するために、なぜ片方だけではなく両方を確認する必要があるのか 一文で説明しなさい。

  2. CPHF の右辺 \(b\) は核変位摂動と電場摂動でまったく異なりますが、両方とも同じ linear_solve(apply_aplusb, ...) 呼び出しで解かれます。左辺演算子 \(A+B\) の何が摂動非依存に しているのですか、そしてそれがなぜここでのコード再利用を可能にするのですか。

  3. 固定(密度非依存)の ECP や点電荷外部ポテンシャルは、応答カーネル \(\sigma\) を構築する際には完全に 無視できるのに、通常の SCF Fock 行列やエネルギーを構築する際にはまったく無視されないのはなぜ ですか。

次章はこの機構を直接引き継ぎます:解析的 Hessian は構造的に「\(3N_{\text{atom}}\) 回の CPHF 解と 1 回の折り返し縮約」にほかならず — 本章のすべての式がそこに、核座標摂動に特化して、 再登場します。