SCF 収束理論:DIIS・SOSCF・拡大 Hessian 法

SCF 収束理論:DIIS・SOSCF・拡大 Hessian 法#

SCF の章は収束制御を実務的なメニューとして示しました。本章はそれらを導出します: すべての戦略を「軌道勾配をゼロへ導く」ことの特殊例にする統一描像、CDIIS/EDIIS/ADIIS の完全な構成、安定化法（ダンピング・レベルシフト・Fermi スメアリング）、そして 1 次法が仕上げられないものを仕上げる 2 次法族（SOSCF・拡大 Hessian・TRAH）です。

統一描像:制約なし最適化としての SCF#

現在の軌道 $\mathbf C$ からユニタリ回転で到達できる新しい軌道集合をパラメータ化します、

\[ \mathbf C(\kappa) = \mathbf C\,\exp(\kappa), \qquad \kappa_{ai}=-\kappa_{ia}^{*},\quad \kappa_{ij}=\kappa_{ab}=0, \]

占有–仮想ブロック（占有を $i,j$、仮想を $a,b$）に限ります。占有どうし・仮想どうしの回転は Slater 行列式、したがってエネルギーを不変に保つので、非冗長なのは占有–仮想ブロックだけだからです。エネルギーを $\kappa$ で 2 次まで展開します:

\[ E(\kappa) = E_0 + \mathbf g^{\mathsf T}\kappa + \tfrac12\,\kappa^{\mathsf T}\mathbf H\,\kappa + \mathcal O(\kappa^3). \]

RHF では軌道勾配と軌道 Hessianは

\[ g_{ai} = 4F_{ai}, \qquad H_{ai,bj} = \delta_{ij}F_{ab} - \delta_{ab}F_{ij} + 4(ai|bj) - (ab|ij) - (aj|bi), \]

$F_{ai}$ は MO Fock 行列の占有–仮想ブロック、$(pq|rs)$ は MO 二電子積分です（UHF はスピンチャネルごとに和; ROHF は 3 ブロックの閉殻/開殻/仮想への一般化が必要で、Hessian が厳密であり続けるよう完全な交差ブロック Fock 交換子を導出します）。停留条件 $\mathbf g=0$ が SCF の収束条件そのものであり、これは AO 基底での厳密に等価な形を持ちます:交換子 $\mathbf{FDS}-\mathbf{SDF}$ は自己無撞着で消えます。この交換子はまさに、占有–仮想の勾配ブロックを MO 表現ではなく AO 表現で書いたものだからです。この 1 つの等価性が、以下のどの方法も — その簿記がどれほど異なっても — 実は「$\mathbf g\to0$ へ導くステップへの異なる近似」である理由です:

方法群	使うもの	性格
CDIIS / EDIIS / ADIIS	$\mathbf g$ のみ（サイクルをまたいで外挿）	1 次
ダンピング / レベルシフト / スメアリング	$\mathbf g$、安定化つき	1 次 + 正則化
SOSCF	$\mathbf g$ + 近似的な $\mathbf H^{-1}$（BFGS）	準 2 次
拡大 Hessian（QC-SCF）/ TRAH	$\mathbf g$ + 厳密な $\mathbf H$、Hessian ベクトル積経由	2 次

どの方法も $\mathbf H$ を陽に作りません（$\mathcal O(n_{\text{ao}}^4)$ の保存）— すべての 2 次戦略は Hessian ベクトル積 $\sigma=\mathbf H\kappa$ だけを評価し、コストは真の密度ではなく遷移密度との Fock 様縮約 1 回分です。この応答エンジン — RHF/UHF/ROHF で厳密、有効なら厳密な KS 交換相関核 $f_{xc}$ と PCM 反応場応答を含む — こそが、あらゆる 2 次法と解析的 Hessian（後の章）が築く土台の唯一の機構です。

1 次法:勾配を外挿する#

CDIIS（Pulay）#

反復部分空間の直接反転（Pulay, 1980/1982）は、過去 $n$ サイクル分の Fock 行列から、各サイクルの交換子を誤差信号として、より良い Fock 行列を外挿します。

定義 3 (DIIS 誤差ベクトル)

正規直交基底（$\mathbf X=\mathbf S^{-1/2}$）で、サイクル $i$ の誤差は $$ \mathbf e_i = \mathbf X^{\mathsf T}\big(\mathbf F_i\mathbf D_i\mathbf S - \mathbf S\mathbf D_i\mathbf F_i\big)\mathbf X, $$ 自己無撞着で厳密にゼロになります（これは AO 基底の軌道勾配そのものです）。

アルゴリズム 1 (CDIIS 外挿)

入力: 直近 $n$ 個の Fock/密度対 $\{\mathbf F_i,\mathbf D_i\}$。 出力: 次に対角化すべき外挿 $\bar{\mathbf F}$。

各履歴について誤差 $\mathbf e_i$ を作る（上の定義）。
Gram 行列 $B_{ij}=\operatorname{Tr}(\mathbf e_i^{\mathsf T}\mathbf e_j)$ を作る。
$\sum_i c_i=1$ の拘束下で $\|\sum_i c_i\mathbf e_i\|^2$ を最小化する。拘束に Lagrange 乗数 $\lambda$ を導入し Lagrangian の勾配をゼロと置くと、拡大した線形系 $$

(1)#\[\begin{pmatrix}\mathbf B & -\mathbf 1\\ -\mathbf 1^{\mathsf T} & 0\end{pmatrix}\]

(2)#\[\begin{pmatrix}\mathbf c\\ \lambda\end{pmatrix} = \begin{pmatrix}\mathbf 0\\ -1\end{pmatrix}\]

$$ が得られる。
$\bar{\mathbf F}=\sum_i c_i\mathbf F_i$ を外挿し、最新の $\mathbf F_n$ の代わりにこれを対角化する。

CDIIS は解の近くで超線形に収束します — それが全強みですが、誤差ノルムの最小化はエネルギーを下げる ことと同じではありません:履歴がまだ自己無撞着から遠いとき、外挿した $\bar{\mathbf F}$ がエネルギーを上げるステップになるのを妨げるものは構成上何もありません。実務上これは、貧弱な推定から始めると不安定・発散として現れ、下記のエネルギーベースの代替法の動機になります。さらに 2 つの実務的な安全策が重要です:悪条件の $\mathbf B$（ほぼ縮退した誤差ベクトル）は最古の履歴ベクトルの除去か対角の正則化が必要で、係数 $|c_i|$ には上限が必要です — 制約なし最小二乗解は時に、ひどく行き過ぎた外挿係数を生みえます。

EDIIS と ADIIS:制約付きエネルギー最小化#

EDIIS（Kudin, Scuseria, Cancès, 2002）は代わりにエネルギーモデルを直接最小化します。密度の凸結合 $\mathbf D(\mathbf c)=\sum_i c_i\mathbf D_i$（$c_i\ge0,\sum_i c_i=1$）に対し、HF エネルギーが密度の 2 次汎関数であることから、試行エネルギーは閉形式を持ちます:

\[ E^{\text{EDIIS}}(\mathbf c) = \sum_i c_i E_i - \tfrac12\sum_{i,j}c_ic_j\langle\mathbf F_i-\mathbf F_j,\, \mathbf D_i-\mathbf D_j\rangle, \]

$E_i$ はサイクル $i$ の全エネルギー、$\langle\cdot,\cdot\rangle$ は Frobenius 内積です。これは小さな（履歴サイズの）箱型制約付き二次計画で、アクティブセット法や射影勾配法で厳密に安価に解けます — そして拘束 $c_i\ge0$ こそが EDIIS を収束から遠くても頑健にするものです:負の係数はまさに CDIIS の行き過ぎを許すものであり、EDIIS は構造的にそれを禁じます。

ADIIS（Hu, Yang, 2010）は EDIIS の DFT 安全な改良版です:2 次エネルギー（KS-DFT では $E_{xc}$ が密度の 2 次でないため成り立たない）を仮定する代わりに、最新の反復 $n$ を参照した拡大 Roothaan–Hall モデルを組み立てます:

\[ E^{\text{ADIIS}}(\mathbf c) = 2\sum_i c_i\langle\mathbf D_i-\mathbf D_n,\,\mathbf F_n\rangle + \sum_{i,j}c_ic_j\langle\mathbf D_i-\mathbf D_n,\,\mathbf F_j-\mathbf F_n\rangle, \]

同じく $c_i\ge0,\sum c_i=1$ の下で最小化します。ADIIS と EDIIS は Hartree–Fock で厳密に一致し（2 次の仮定が厳密だから）、ADIIS は DFT でより頑健な選択です。qc-rs の algorithm="auto" は誤差ノルム $\|\mathbf e_n\|$ を監視してこれらのエネルギーベース外挿を CDIIS と混ぜます:$\|\mathbf e_n\|$ が大きい間はエネルギーベース外挿、中間域は線形補間、$\|\mathbf e_n\|$ が厳しいしきい値を下回れば純粋な CDIIS — そこで CDIIS の超線形収束は他に並ぶものがありません。

安定化法#

さらに 3 つの道具が、上記のどの戦略も置き換えずに安定させます。

ダンピングは新旧の密度を線形に混ぜます、 $$ \tilde{\mathbf D}^{(n)} = (1-\alpha)\,\mathbf D^{(n)}_{\text{new}} + \alpha\,\mathbf D^{(n-1)},\qquad 0\le\alpha<1, $$ 序盤の振動への鈍いが効果的なブレーキです（典型的な静的値は $\alpha\approx0.7$ で誤差がしきい値を下回れば解除; 動的な Hehenberger–Zerner 方式はエネルギー/勾配比から毎サイクル $\alpha$ を調整します）。

レベルシフト（Saunders–Hillier, 1973）は仮想軌道エネルギーを定数 $b$ だけ持ち上げ、小さな HOMO–LUMO ギャップ付近での偽の占有–仮想混合を抑えます。MO 基底では単に $p,q\in\text{仮想}$ について $\tilde F_{pq}=F_{pq}+b\,\delta_{pq}$;仮想空間射影子 $\mathbf Q=\mathbf S^{-1}-\mathbf D$ を使った等価な AO 基底形は、 $$ \tilde{\mathbf F} = \mathbf F + b\,\mathbf S\mathbf Q\mathbf S = \mathbf F + b\,(\mathbf S - \mathbf S\mathbf D\mathbf S), $$ これを AO Fock 行列に直接加えられます。シフトが大きすぎると収束の終盤を遅くします（安定化するはずのギャップを人為的に膨らませるため）ので、誤差が縮むにつれ解除します。

Fermi スメアリングは整数の Aufbau 占有を有限温度の Fermi–Dirac 分布で置き換えます、 $$ f_i = \Big[1+\exp\big((\varepsilon_i-\mu)/k_BT_{\text{el}}\big)\Big]^{-1}, $$ 化学ポテンシャル $\mu$ を、電子数保存 $\sum_i w_if_i=N_{\text{el}}$ から二分法/Newton 法で解きます（$f_i$ は $\mu$ について単調なので根は一意）。密度は $\mathbf D=\sum_if_i|\mathbf c_i\rangle\langle \mathbf c_i|$ となり、分数占有は真の自由エネルギー最小化なので、収束は素のエネルギーではなく電子自由エネルギー $$ A = E - T_{\text{el}}S_{\text{el}}, \qquad S_{\text{el}} = -k_B\sum_i\big[f_i\ln f_i+(1-f_i)\ln(1-f_i)\big], $$ を追うべきです。スメアリングはまさに**準縮退や金属的**な系 — 整数占有が毎サイクル軌道を出入りさせ持続的な振動を引き起こす — の修正です;よくあるアニーリングスケジュールは高温で始め、SCF が収束するにつれ $T_{\text{el}}\to0$ に冷やし、最後には通常の整数占有解を回復します。

2 次法:厳密な Hessian を使う#

1 次法が停滞する（典型的には $\|\mathbf e\|\sim10^{-3}$〜$10^{-4}$）、あるいはまったく失敗するとき、 2 次法族は $\mathbf H$ を使った明示的な Newton 風ステップを取ります。

SOSCF（Chaban–Schmidt–Gordon, 1997）は $\mathbf H$ をまったく作りません — その逆を、 (s,y) 対の履歴から BFGS 準ニュートン更新で近似し、安価な対角初期値 $H^{(0)}_{ai,ai}\approx4(\varepsilon_a-\varepsilon_i)$ から始め、$\kappa=-\mathbf H^{-1}\mathbf g$ でステップします。Hessian ベクトル積が一切不要なので最も安価な 2 次法で、解近くで超線形収束しますが遠方から始めると不安定です — これが、勾配が既に小さいときに作動する DIIS の仕上げであり、開始戦略ではない理由です。

拡大 Hessian / QC-SCF（Bacskay, 1981）は Newton 法の実際の失敗様式を直接克服します:素のステップ $\mathbf H\kappa=-\mathbf g$ は、$\mathbf H$ が負の固有値を持つ（収束から遠い場所ではよくある — まさに安定性解析が直接検査する署名です）と、そもそも降下すらしません。修正は、最低固有対について拡大固有値問題を解くことです、 $$

(3)#\[\begin{pmatrix}0 & \alpha\mathbf g^{\mathsf T}\\ \alpha\mathbf g & \mathbf H\end{pmatrix}\]

(4)#\[\begin{pmatrix}1\\ \tilde\kappa\end{pmatrix} = \mu\begin{pmatrix}1\\ \tilde\kappa\end{pmatrix}\]

;\Longrightarrow; (\mathbf H-\mu\mathbf I)\kappa=-\mathbf g, $$ 固有値 $\mu<0$ は**自動的なレベルシフト**です:$\mathbf H-\mu\mathbf I$ は構成上正定値なので、$\mathbf H$ 自体が不定でも結果のステップは降下します。最低固有対だけが必要なので、**Davidson** 反復が拡大行列の積をその場で繰り返し評価して見つけます — $\mathbf H$ を組み立てることは決してありません。

TRAH（信頼領域拡大 Hessian;Helmich-Paris, 2021）は同じ拡大問題に明示的な信頼領域を加え、 $\mathbf H$ を試行ステップ長 $\lambda\ge1$ でスケールします、 $$

(5)#\[\begin{pmatrix}0 & \mathbf g^{\mathsf T}\\ \mathbf g & \mathbf H/\lambda\end{pmatrix}\]

(6)#\[\begin{pmatrix}1\\ x\end{pmatrix} = \varepsilon\begin{pmatrix}1\\ x\end{pmatrix}\]

;\Longrightarrow; \kappa=\frac{x}{\lambda} = -(\mathbf H-\lambda\varepsilon\mathbf I)^{-1}\mathbf g, $\lambda$ を、結果の $|\kappa|$ が目標信頼半径に合うよう選び、その後、予測エネルギー低下（モデルから）と実際の低下を比べる $\rho$ テストで半径を広げるか縮めるかします。ある重要な実装上の恒等式 — 降下は $\mathbf g^{\mathsf T}\kappa=\varepsilon/\lambda\le0$ で保証され $\mathbf H$ の正定値性を必要としない — により、各試行 $\lambda$ で高価な Hessian ベクトル積を再計算せずに信頼半径を適応させられます。TRAH の適応半径が、最も難しい系（軌道的にほぼ縮退した ROHF ラジカルなど）で頑健にするもので、固定ステップは誤った電子ベイスンへ行き過ぎかねません。

収束判定#

実用の SCF コードはエネルギーだけでなく複数の基準を同時に要求します:

基準	記号	典型的な厳しいしきい値
エネルギー変化	$	\Delta E
密度 RMS 変化	$\operatorname{RMS}(\Delta\mathbf D)$	$5\times10^{-9}$
密度最大変化	$\max	\Delta\mathbf D
DIIS 誤差（交換子）	$\|\mathbf e\|$	$5\times10^{-7}$
軌道勾配	$\|\mathbf g\|$	$10^{-5}$

2 次法は自然に $\|\mathbf g\|$ だけで判定されます（それがその手法が直接ゼロへ導く量です）;1 次法は複数基準の AND から恩恵を受けます。エネルギーの平坦域だけでは、SCF が密度において実際には自己無撞着に達していないことを見逃しうるからです。

推奨される戦略の梯子#

        flowchart TD
    G["初期推定（SAD / SAP）"] --> F["遠方: ||e|| > 0.1"]
    F -->|"EDIIS（HF）/ ADIIS（DFT）<br/>必要ならダンピング/レベルシフトも"| M["中盤: 1e-3 < ||e|| <= 0.1"]
    M -->|"エネルギーモデル + CDIIS の混合"| N["近傍: ||e|| <= 1e-3"]
    N -->|"純 CDIIS（局所最速）"| T{"停滞？"}
    T -->|はい| S["SOSCF（超線形の仕上げ）"]
    T -->|いいえ、収束| D["完了"]
    N -->|"停滞 / 振動"| SM["Fermi スメアリング + T -> 0 アニール"]
    SM -->|"それでも難しい"| AH["拡大 Hessian / TRAH<br/>（大域収束保証）"]
    S --> D
    AH --> D

どの段も「$\mathbf g\to0$ への異なる近似」なので、梯子は分子がどの段を必要としても常に同じエネルギーに収束します — まさに SCF の章での実証どおりです。algorithm="xqc" はよくある場合を自動化し（有限回数の素の DIIS、未収束なら拡大 Hessian の梯子にフォールバック）、algorithm="yqc" は軌道的に縮退した開殻（ROHF/ROKS）向けに同じ考えを調整したものです。これらは近傍領域に到達するずっと前に素の 1 次 DIIS を停滞させるので、既定で素の diis ではなくこのエスカレーションを使います。

練習 4

なぜ EDIIS の係数は $c_i\ge0$ を満たさねばならず CDIIS はそうでないのですか。その拘束だけで、EDIIS が収束から遠くても頑健な理由をどう説明できますか。
拡大 Hessian の固有値方程式は、$\mathbf H$ が負の固有値を持っていても降下を保証します。この構成の中でそれを真にする唯一の量を特定し、なぜそうなのか一文で説明しなさい。
TRAH は各試行信頼半径で Hessian ベクトル積を再計算する代わりに $\mathbf H$ を $1/\lambda$ でスケールし直します。これが TRAH のマクロ反復あたりのコストにとってなぜ重要ですか。

解答練習 4

$c_i\ge0$ は $\mathbf D(\mathbf c)=\sum_ic_i\mathbf D_i$ を、それらが張る単体の内側の点である履歴密度の凸結合に制限します。CDIIS の制約なし係数は負になりえ、これは代数的には履歴密度を 超えた外挿です — まさに CDIIS を収束から遠くで不安定にする種類の行き過ぎであり、まさに $c_i\ge0$ が禁じるものです。
自動的なレベルシフト $\mu$（拡大問題の最低固有値）です。$(\mathbf H-\mu\mathbf I)\kappa=-\mathbf g$ であり $\mu$ は構成上 $\mathbf H$ の最低固有値以下なので、$\mathbf H-\mu\mathbf I$ は $\mathbf H$ 自体の符号構造によらず半正定値となり、降下方向を保証します。
Hessian ベクトル積（応答エンジン経由）はマクロ反復の高価な部分です。Krylov/Davidson 部分空間行列を $\tilde M = \tilde M_1 + \tilde M_2/\lambda$ として分離保存することで、TRAH は複数の試行信頼半径 $\lambda$ を、既に計算済みの部分空間行列上の安価な線形代数だけで再試行できます — 試行ごとに新しい $\sigma=\mathbf H\kappa$ の評価は不要です。

$\sigma=\mathbf H\kappa$ の背後の応答エンジンは本マニュアルにあと 2 回登場します:解析的 Hessian の背後の結合摂動方程式として、そして安定性解析が直接検査する演算子の固有値として。

方法群	使うもの	性格
CDIIS / EDIIS / ADIIS	\(\mathbf g\) のみ（サイクルをまたいで外挿）	1 次
ダンピング / レベルシフト / スメアリング	\(\mathbf g\)、安定化つき	1 次 + 正則化
SOSCF	\(\mathbf g\) + 近似的な \(\mathbf H^{-1}\)（BFGS）	準 2 次
拡大 Hessian（QC-SCF）/ TRAH	\(\mathbf g\) + 厳密な \(\mathbf H\)、Hessian ベクトル積経由	2 次