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

[SCF の章](../20-guide/scf.md) は収束制御を実務的なメニューとして示しました。本章はそれらを導出します:
すべての戦略を「軌道勾配をゼロへ導く」ことの特殊例にする統一描像、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 行列を外挿します。

:::{prf:definition} DIIS 誤差ベクトル
:label: def-diis-error

正規直交基底（$\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 基底の軌道勾配*そのもの*です）。
:::

:::{prf:algorithm} CDIIS 外挿
:label: alg-cdiis-full

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

1. 各履歴について誤差 $\mathbf e_i$ を作る（上の定義）。
2. Gram 行列 $B_{ij}=\operatorname{Tr}(\mathbf e_i^{\mathsf T}\mathbf e_j)$ を作る。
3. $\sum_i c_i=1$ の拘束下で $\|\sum_i c_i\mathbf e_i\|^2$ を最小化する。拘束に Lagrange 乗数 $\lambda$ を
   導入し Lagrangian の勾配をゼロと置くと、拡大した線形系
   $$
   \begin{pmatrix}\mathbf B & -\mathbf 1\\ -\mathbf 1^{\mathsf T} & 0\end{pmatrix}
   \begin{pmatrix}\mathbf c\\ \lambda\end{pmatrix} = \begin{pmatrix}\mathbf 0\\ -1\end{pmatrix}
   $$
   が得られる。
4. $\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$ が負の固有値を持つ（収束から遠い場所ではよくある — まさに
[安定性解析](../20-guide/scf.md) が直接検査する署名です）と、そもそも降下すらしません。修正は、
最低固有対について**拡大固有値問題**を解くことです、
$$
\begin{pmatrix}0 & \alpha\mathbf g^{\mathsf T}\\ \alpha\mathbf g & \mathbf H\end{pmatrix}
\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$ でスケールします、
$$
\begin{pmatrix}0 & \mathbf g^{\mathsf T}\\ \mathbf g & \mathbf H/\lambda\end{pmatrix}
\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|$ | $10^{-8}$ |
| 密度 RMS 変化 | $\operatorname{RMS}(\Delta\mathbf D)$ | $5\times10^{-9}$ |
| 密度最大変化 | $\max|\Delta\mathbf D|$ | $10^{-7}$ |
| DIIS 誤差（交換子） | $\|\mathbf e\|$ | $5\times10^{-7}$ |
| 軌道勾配 | $\|\mathbf g\|$ | $10^{-5}$ |

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

## 推奨される戦略の梯子

```{mermaid}
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 の章](../20-guide/scf.md) での実証どおりです。`algorithm="xqc"` はよくある場合を
自動化し（有限回数の素の DIIS、未収束なら拡大 Hessian の梯子にフォールバック）、`algorithm="yqc"` は軌道的に
縮退した開殻（ROHF/ROKS）向けに同じ考えを調整したものです。これらは近傍領域に到達するずっと前に素の 1 次
DIIS を停滞させるので、既定で素の `diis` ではなくこのエスカレーションを使います。

:::{exercise}
:label: ex-scf-theory

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

:::{solution} ex-scf-theory
:class: dropdown

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

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