# SCF: Hartree–Fock と KS-DFT
分子・基底・初期推定が揃いました。本章は、それらを軌道とエネルギーに変える計算 — **自己無撞着場（SCF）** — を
実行します。SCF は **Hartree–Fock** と **Kohn–Sham DFT** の両方を解く、以降のあらゆる手法の土台です。だから
ガイド中で最も長い章になりますが、制御は少なく、たいていの計算は 2 つ（`ref` と `xc`）だけで済みます。
## 理論の復習:SCF は何を解くのか
[Part II](../10-foundations/hartree-fock.md) より:HF も KS-DFT も、有限基底での一電子固有値問題 —
**Roothaan 方程式** — に帰着します、
$$
\mathbf{F}\,\mathbf{C} = \mathbf{S}\,\mathbf{C}\,\boldsymbol\varepsilon ,
$$
ここで Fock（または Kohn–Sham）行列 $\mathbf F$ は、それが生む軌道 $\mathbf C$ から作られる密度 $\mathbf D$ に
依存します。この循環ゆえに方程式は**反復的**に解かれます — $\mathbf F$ を作り、対角化し、密度を作り直し、変化が
止まるまで繰り返す（[初期推定の章](initial-guess.md) がそのループを始めました）。qc-rs がサイクル全体を扱います;
あなたの仕事は、**どの電子構造**を解くか（参照法、DFT なら汎関数）を指定し、時折**難しい収束をどう導くか**を
決めることです。
## 基本の実行
```python
import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
done = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run()
done.scf.energy # -76.026794 全エネルギー（hartree）
done.scf.converged # True
done.scf.ncycle # 9 使った SCF 反復数
done.scf.energy_elec # -85.220707 電子エネルギー（= 全エネルギー − 核間反発）
done.scf.energy_components # {'core': -123.1493, 'coulomb': 46.9052, 'exchange': -8.9766}
```
`.scf(...)` は保留 SCF ステップを追加し、`.run()` が実行します（関数形式 `qc.scf(mychk, ...)` と等価）。
**`energy_components`** は電子エネルギーを物理的な要素に分解します — HF では一電子の `core`（運動 + 核引力）、
クーロン反発 `coulomb`（$J$）、交換 `exchange`（$K$）;DFT ではこの最後が `xc` 項になります（後述）。
## 参照法を選ぶ:RHF, UHF, ROHF
参照法は **`ref`** とスピン多重度 **`spin` = 2S+1**（[`qc.chk.new`](molecular-input.md) で設定）で決まります。
3 つの参照法は、スピン軌道をどう制約するかで異なります（[Part II](../10-foundations/hartree-fock.md)）:
| `ref` | 閉殻（`spin=1`） | 開殻（`spin>1`） |
|---|---|---|
| `"auto"` *(既定)* | **RHF** / RKS | **UHF** / UKS |
| `"r"` | **RHF** / RKS | **ROHF** / ROKS |
| `"u"` | UHF / UKS | UHF / UKS |
| `"ro"` | **エラー** | **ROHF** / ROKS |
**`ref="auto"`（既定）** は閉殻には*制限*、開殻には*非制限*を選びます — なのでラジカルや三重項は、指定なしで
UHF になります（UHF が慣習的な開殻の既定）。`ref="ro"` は制限開殻を明示的に要求し、**閉殻ではエラー**にします
（黙って RHF を走らせない）— 誤った `charge`/`spin` に対する歯止めです。
### UHF と ROHF:トレードオフ
開殻では 2 つの参照法から選び、それらは*異なる*エネルギーを与えます:
```python
ch3 = "C 0 0 0; H 0 1.079 0; H 0.934 -0.539 0; H -0.934 -0.539 0" # メチルラジカル、二重項
u = qc.chk.new(atom=ch3, ao="cc-pvdz", unit="angstrom", spin=2).scf(ref="u").run().scf
ro = qc.chk.new(atom=ch3, ao="cc-pvdz", unit="angstrom", spin=2).scf(ref="ro").run().scf
print(u.energy, ro.energy) # -39.563802 -39.559636
```
- **UHF** は α と β 電子が*異なる*空間軌道を占めるのを許します。その自由度で**より低いエネルギー**
（−39.563802 対 −39.559636）を与えますが、波動関数はもはや純粋なスピン状態ではなく、**スピン汚染**を受けます。
これは $\langle S^2\rangle$ が厳密値（二重項なら $S(S+1)=0.75$）を超えて漂うこととして測れます:
```python
qc.prop.spin.s_squared(qc.chk.new(atom=ch3, ao="cc-pvdz", unit="angstrom", spin=2).scf(ref="u").run())
# 0.7612 （厳密な二重項 = 0.75;わずかな超過が汚染）
```
- **ROHF** は対電子に共通の空間軌道を強制するので、**スピン純粋**（$\langle S^2\rangle$ 厳密）を保ちますが、
エネルギーはわずかに高くなります。
ラジカルには **UHF/UKS** を既定に、スピン純粋な参照が重要なとき（後続の相関計算の土台にする、汚染由来の
アーティファクトを避ける等）には **ROHF/ROKS** を使います。
:::{important} 準縮退の開殻 → ROHF ではなく UHF
準縮退や非局在ホールの状態 — 例えばホールが等価なサイトに広がる対称なクラスターカチオン — は実質的に多参照で、
単一行列式の ROHF では**表現できず**、収束しません。UHF は対称性を破って（ホールを局在化して）そこへ到達でき、
正しい選択です。ROHF/ROKS は明確に定義された高スピン/二重項ラジカル向けです。
:::
## Hartree–Fock か KS-DFT か:`xc=`
*同じ* SCF 機構が両手法を走らせ、**`xc=`** が唯一のスイッチです:
- **`xc=None`**（既定）または **`xc="hf"`** → **Hartree–Fock**。
- **`xc="pbe"`**, **`"b3lyp"`**, … → その libxc 汎関数での **Kohn–Sham DFT**（LDA / GGA / meta-GGA /
ハイブリッド;理論とヤコブの梯子は [Part II](../10-foundations/dft-kohn-sham.md)）。
```python
pbe = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc="pbe").run().scf
print(pbe.energy) # -76.333409
print(dict(pbe.energy_components)) # {'core': -123.1775, 'coulomb': 46.9303, 'xc': -9.2801}
```
要素が HF の `exchange` の代わりに **`xc`** 項（交換相関）を報告している点に注目 — HF と KS 実行の唯一の構造差です。
### DFT グリッド
KS-DFT は交換相関エネルギーを分子グリッド上の**数値積分**で評価し、**`grid=`** で選びます（ORCA-DefGrid 風の
枝刈りされた段階）:
| `grid=` | 動径 × 角度 | 用途 |
|---|---|---|
| `coarse` | 50 × 194 | 下書き / 手早く |
| `medium` *(既定)* | 60 × 302 | 本番（≈ ORCA DefGrid2） |
| `fine` | 75 × 434 | 厳密 |
| `ultrafine` | 200 × 5810 | 参照 / 最高精度 |
既定の `medium` は密な `ultrafine` 参照の ≲0.4 µEh/原子 以内に較正されているので、変える必要はまれです;厳しい
エネルギー差や、meta-GGA がグリッド感受性を示すときは `fine` に上げます。（HF はグリッドなし — 解析的です。）
## 理論:最適化問題としての SCF
収束「戦略」が何を近似しているのかを理解するため、Roothaan 方程式を**軌道回転についての制約なし最適化**として
書き直します。現在の軌道 $\mathbf C$ からユニタリ回転で到達できる新しい軌道集合は
$$
\mathbf C(\kappa) = \mathbf C\,\exp(\kappa), \qquad \kappa_{ai} = -\kappa_{ia}^{*},\quad \kappa_{ii}=\kappa_{aa}=0,
$$
と書けます。ここで $\kappa$ は**占有–仮想**ブロック（占有を $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),
$$
**軌道勾配** $\mathbf g$ と**軌道 Hessian** $\mathbf H$ が現れます。RHF では、
$$
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 二電子積分です。勾配が消えること
$\mathbf g=0$ が、まさに SCF の収束条件です — そしてこれは、お馴染みの AO 交換子
$\mathbf{FDS}-\mathbf{SDF}=0$ と等価です。この交換子こそが、勾配を AO 基底で書いた占有–仮想ブロックだから
です。この 1 つの事実が、以下のすべての戦略を結びつける糸です:
:::{important} 1 つの描像、多くのアルゴリズム
**あらゆる収束戦略は「$\mathbf g\to 0$ に向かうステップを取る」ことへの、異なる近似**です。1 次法（DIIS）は
$\mathbf g$ だけを（サイクルをまたいで外挿して）使い、2 次法（SOSCF、拡大 Hessian 系）は真の Newton 様ステップ
のため $\mathbf H$ も（厳密または近似的に）使います。$\mathbf H$ を陽に作るのは $\mathcal O(n_{\text{ao}}^4)$
かかり決して行われません — 代わりにすべての 2 次法は **Hessian ベクトル積** $\sigma=\mathbf H\kappa$ だけを
評価します。これは（真の密度ではなく「遷移密度」との縮約という形で）1 回追加の Fock 構築と同じコストです。
:::
## 理論:DIIS — 反復部分空間の直接反転
**DIIS**（Pulay, 1980/1982）が既定の 1 次加速法です。現在の Fock 行列をそのまま使う代わりに、過去数サイクル分の
Fock 行列とその交換子誤差を手がかりに、より良い Fock 行列を**外挿**します。
:::{prf:algorithm} CDIIS（Pulay）ステップ
:label: alg-cdiis
**入力:** 直近 $n$ 個の Fock/密度対 $\{F_i, D_i\}$;直交化子 $X=\mathbf S^{-1/2}$。
**出力:** 対角化すべき外挿 Fock 行列 $\bar F$。
1. 各履歴 $i$ について、正規直交基底での**誤差行列**を作る:
$$
e_i = X^{\mathsf T}\big(F_i D_i S - S D_i F_i\big) X .
$$
自己無撞着で $e_i \to 0$（これは上の軌道勾配の AO 形です）。
2. $n\times n$ の Gram 行列 $B_{ij} = \operatorname{Tr}\!\big(e_i^{\mathsf T} e_j\big)$ を作る。
3. 外挿係数 $\mathbf c$（拘束 $\sum_i c_i = 1$ を乗数 $\lambda$ で課す）について、Lagrange 拡大した線形系
を解く:
$$
\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 F = \sum_i c_i F_i$ を外挿し、最新の $F_n$ ではなく*これ*を対角化する。
:::
CDIIS は誤差ノルム $\|\sum_i c_i e_i\|$ を最小化しますが、そのノルムはエネルギーではありません — 収束から
遠い場所では、外挿した $\bar F$ がエネルギーを*上げる*ステップを妨げるものは何もなく、これが CDIIS 単独では
序盤に不安定になりうる理由です。だから **`auto`** 戦略は、誤差が大きい間は**エネルギーベース**の外挿
（EDIIS/ADIIS — エネルギーそのものの 2 次モデルを、$c_i\ge0$ の制約下で最小化するので行き過ぎない）を混ぜ、
$\|e\|$ がしきい値を下回れば純粋な CDIIS に切り替えます。CDIIS の超線形局所収束が最も速いのはそこだからです。
これが「遠方はエネルギーに頑健、近傍は誤差に頑健」の橋渡しで、Part II の
[収束理論の章](../10-foundations/scf-convergence-theory.md) で詳述します。
## 理論:2 次収束 — SOSCF・拡大 Hessian・TRAH
DIIS が trailing する（典型的には $\|e\|\sim10^{-3}$〜$10^{-4}$ 付近）、あるいはまったく収束しないとき、
2 次法群は Hessian を使って $\kappa$ に明示的なステップを取ります:
- **SOSCF**（Chaban–Schmidt–Gordon）は $\mathbf H$ をまったく作らず — その*逆*を **BFGS** 準ニュートン更新
で近似します。安価な対角初期値 $H^{(0)}_{ai,ai}\approx 4(\varepsilon_a-\varepsilon_i)$ から始め、
$\kappa=-\mathbf H^{-1}\mathbf g$ でステップします。最も安価な 2 次法（Hessian ベクトル積が一切不要）で、
解に近づけば超線形収束しますが遠方では不安定になりえます — だから `soscf` は開始戦略ではなく DIIS の
*仕上げ*として最適です。
- **QC-SCF / 拡大 Hessian**（Bacskay）は Newton 法の失敗様式を直接修正します:素の Newton ステップ
$\mathbf H\kappa=-\mathbf g$ は、$\mathbf H$ が負の固有値を持つ（収束から遠い場所ではよくある）と、
そもそも降下すらしません。代わりに**拡大固有値問題**
$$
\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）は明示的な**信頼半径**を加えます:同じ拡大問題を、$\mathbf H$ を試行
ステップ長 $\lambda\ge1$ でスケールして解き、結果の $\|\kappa\|$ が目標半径に合うよう選びます — ステップが
受理されれば半径を広げ（予測と実際のエネルギー低下を比べる `ρ` テスト）、そうでなければ縮めます。これが
`trah` を最も難しい開殻系（軌道的にほぼ縮退した ROHF ラジカルなど）で頑健にするもので、固定ステップだと
誤った電子状態へ行き過ぎかねない場面です。
```{mermaid}
flowchart TD
G["初期推定密度 D₀（sad, gwh, ...）"] --> F["Fock F を構築"]
F --> E["誤差 e = FDS − SDF"]
E --> Q{"‖e‖ としきい値"}
Q -->|"大きい（遠方）"| ED["EDIIS / ADIIS
エネルギーモデルのステップ"]
Q -->|"小さい（近傍）"| CD["CDIIS
Pulay 外挿"]
Q -->|"停滞"| SO["2 次の仕上げ
SOSCF / QC-SCF / TRAH"]
ED --> F
CD --> F
SO --> F
Q -->|"収束"| D["完了:E, C, D"]
```
これはまさに `algorithm="auto"` が RHF/UHF/RKS/UKS で辿る梯子です（平衡付近ではふつう素の DIIS で十分なので、
エネルギーモデル分岐や 2 次分岐はめったに作動しません）;ROHF/ROKS が既定で `yqc` なのは、軌道的に縮退した
開殻が近傍領域に到達するずっと前に素の 1 次 DIIS を停滞させるからです。完全な導出 — Hessian ベクトル積が
なぜ $\mathbf H$ を決して材料化しないか、EDIIS/ADIIS のエネルギー汎関数、ROHF の 3 ブロック Hessian を含む —
は [SCF 収束理論](../10-foundations/scf-convergence-theory.md) にあります。
## 収束:SCF を導く
上の理論を手にした今、実務的な制御は 2 つのつまみです:**戦略**選択子 `algorithm=`（上のどの描像を使うか）と、
どれかを安定させる少数の**安定化**引数。[初期推定](initial-guess.md) と同様、**収束制御は*経路*だけを変え、
収束後の答えは変えません** — 以下のどの設定も同じエネルギーに達します。
既定の `algorithm="auto"` は、上の図の梯子を SAD 推定から辿り、平衡付近のたいていの分子を 20 サイクル未満で
収束させます。ふつうは何も変えません。難しいときは**症状で**選びます:
| 症状 | 使うもの | 方法 |
|---|---|---|
| 難しい / 振動 — 頑健な 2 次法が欲しい | QC-SCF（拡大 Hessian） | `algorithm="qc"` |
| 同上、適応的な信頼領域つき | TRAH | `algorithm="trah"` |
| 何も考えず頑健な DIIS | XQC（DIIS + 2 次のセーフティネット） | `algorithm="xqc"` |
| 終盤だけ遅い | SOSCF | `algorithm="soscf"` |
| 収束から遠くで振動 | ダンピング | `damping=0.6` |
| 準縮退 / 小ギャップ / 金属的 | Fermi スメアリング | `smearing=0.01` |
| HOMO/LUMO が毎サイクル入れ替わる | レベルシフト | `level_shift=0.3` |
理論を実地に見てみましょう:拡大 Hessian 系は、難しい開殻系で同じエネルギーにずっと少ないサイクルで到達
します。コストは、素の Fock 構築の代わりにサイクルごとの Hessian ベクトル積です:
```python
for algo in ("diis", "qc", "trah"):
c = qc.chk.new(atom=ch3, ao="cc-pvdz", unit="angstrom", spin=2).scf(ref="u", algorithm=algo).run().scf
print(algo, c.ncycle, round(c.energy, 6))
# diis 11 -39.563802
# qc 6 -39.563802
# trah 6 -39.563802
```
同じエネルギーですが、拡大 Hessian 法は 11 ではなく 6 サイクルで到達します。
### しきい値
SCF は、エネルギー変化と軌道勾配 RMS の**両方**がしきい値を下回ったら停止します:
| しきい値 | 既定 | 設定方法 |
|---|---|---|
| エネルギー `ΔE` | `1e-9` Ha | `scf(conv_tol=...)` または `conv_preset=` |
| `[F, DS]` 交換子の勾配 RMS | `1e-6` | `iop={"scf.conv_tol_grad": ...}` |
手早く締める/緩めるには **`conv_preset="tight"`** / `"loose"`、エネルギーを直接指定するには `conv_tol=` を
使います。細かなノブの背後の IOP キーは [リファレンス](../40-reference/iop.md) にあります。
## 実行結果を読む・診断する
`.run()` の後、`scf` アクセサが結果を公開します;**`run(log=...)`**（[クイックスタート](../00-intro/quickstart.md)
参照）は収束を*見る*ための実況サイクル表を流します:
```python
done = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r").run(log="stdout")
done.scf.energy, done.scf.converged, done.scf.ncycle # 見出しの数値
dict(done.scf.energy_components) # エネルギーの内訳
done.log() # 再計算なしにトランスクリプトを再表示
done.show("result") # 状態 + 結果のスナップショットを描画
```
`converged` が `False` なら、トランスクリプトを見ます:`max_cycle` でまだエネルギーが下がっているなら
「`max_cycle` を上げる、または 2 次の `algorithm` を使う」;エネルギーが*振動*しているなら「`damping`/
`level_shift` を足す、小ギャップなら `smearing`」です。
## 安定性解析
収束は SCF が*ある*停留点を見つけたことを意味するだけで、**最低**とは限りません。`stability=True` は、軌道回転
Hessian に負の固有値がないかを調べ、解が真の極小かを検査します（**internal** = *同じ*参照法のより低い解;
**external** = *より広い*参照法のより低い解、例 RHF→UHF）:
```python
h2 = qc.chk.new(atom="H 0 0 0; H 0 0 1.8", ao="cc-pvdz", unit="angstrom") # 伸ばした H2
res = h2.scf(ref="r", stability=True).run()
dict(res.scf.stability)
# {'internal_stable': True, 'internal_eigenvalue': 0.3995,
# 'external_stable': False, 'external_eigenvalue': -0.1877, 'stable': False}
```
伸ばした H₂ の RHF は **external 不安定**（外部固有値が負）です — 長い結合長では閉殻の制約が誤りで、スピンの
破れた **UHF** 解がより低くにあります。これは結合解離やジラジカルの標準的な診断で、対処は UHF で（多くは対称性を
破るため `guess(..., spin_break="mix")` とともに）再実行することです。
## 総合例と演習
```python
import qc
water = "O 0 0 0.117; H 0 0.757 -0.469; H 0 -0.757 -0.469"
# 同じ分子・基底で HF とハイブリッド汎関数
for xc in (None, "b3lyp"):
s = qc.chk.new(atom=water, ao="cc-pvdz", unit="angstrom").scf(ref="r", xc=xc).run().scf
label = "HF" if xc is None else xc
print(f"{label:6} E = {s.energy:.6f} converged={s.converged} cycles={s.ncycle}")
# HF E = -76.026794 converged=True cycles=9
# b3lyp E = -76.420349 converged=True cycles=8
```
:::{exercise}
:label: ex-scf
1. メチルラジカル CH₃ を `ref="r"`, `spin=2` で実行します。qc-rs はどのソルバを使い、結果はスピン純粋ですか。
`ref="u"` ならどうですか。
2. ある UHF 計算が二重項に対し $\langle S^2\rangle = 1.30$ を与えました。これは許容できますか。何を意味し、
何を試しますか。
3. 伸ばした結合の RHF SCF が*収束*（`converged=True`）しましたが、基底状態ではないと疑っています。これを確認
する 1 つのキーワードは何で、どんな結果が疑いを裏づけますか。
:::
:::{solution} ex-scf
:class: dropdown
1. **開殻**での `ref="r"` は **ROHF** を選び（参照表を参照）、**スピン純粋**（$\langle S^2\rangle$ 厳密）です。
`ref="u"` なら **UHF** — より低いエネルギーですが、わずかにスピン汚染します。
2. $\langle S^2\rangle = 1.30$ は厳密な二重項値 $0.75$ に対し**重い汚染**です — UHF 行列式が高スピン状態と激しく
混ざっており、エネルギーも物性も疑わしいです。**ROHF**（スピン純粋）を試すか、状態が本当に多参照か
（相関/多参照法が必要か）を確認します。
3. **`stability=True`**。負の **external** 固有値（`external_stable: False`）が、より広い参照法のより低い解の
存在を裏づけます — UHF で（`spin_break="mix"` とともに）再実行してそこへ到達します。
:::
SCF は参照波動関数とそのエネルギーを与えます。それが外す**相関**を回収するため、[次章](post-scf.md) は
ポスト SCF 法 — RI-MP2 系列 — を加えます。