SCF: Hartree–Fock と KS-DFT

SCF: Hartree–Fock と KS-DFT#

分子・基底・初期推定が揃いました。本章は、それらを軌道とエネルギーに変える計算 — 自己無撞着場（SCF） — を実行します。SCF は Hartree–Fock と Kohn–Sham DFT の両方を解く、以降のあらゆる手法の土台です。だからガイド中で最も長い章になりますが、制御は少なく、たいていの計算は 2 つ（ref と xc）だけで済みます。

理論の復習:SCF は何を解くのか#

Part II より: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$ を作り、対角化し、密度を作り直し、変化が止まるまで繰り返す（初期推定の章がそのループを始めました）。qc-rs がサイクル全体を扱います; あなたの仕事は、どの電子構造を解くか（参照法、DFT なら汎関数）を指定し、時折難しい収束をどう導くかを決めることです。

基本の実行#

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 で設定）で決まります。 3 つの参照法は、スピン軌道をどう制約するかで異なります（Part II）:

`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 つの参照法から選び、それらは異なるエネルギーを与えます:

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$）を超えて漂うこととして測れます:
```
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 を使います。

重要

準縮退の開殻 → 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）。

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 つの事実が、以下のすべての戦略を結びつける糸です:

重要

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 行列を外挿します。

アルゴリズム 1 (CDIIS（Pulay）ステップ)

入力: 直近 $n$ 個の Fock/密度対 $\{F_i, D_i\}$;直交化子 $X=\mathbf S^{-1/2}$。 出力: 対角化すべき外挿 Fock 行列 $\bar F$。

各履歴 $i$ について、正規直交基底での誤差行列を作る: $$ e_i = X^{\mathsf T}\big(F_i D_i S - S D_i F_i\big) X . $$ 自己無撞着で $e_i \to 0$（これは上の軌道勾配の AO 形です）。
$n\times n$ の Gram 行列 $B_{ij} = \operatorname{Tr}\!\big(e_i^{\mathsf T} e_j\big)$ を作る。
外挿係数 $\mathbf c$（拘束 $\sum_i c_i = 1$ を乗数 $\lambda$ で課す）について、Lagrange 拡大した線形系を解く: $$

(1)#\[\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}. $$
$\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 の収束理論の章で詳述します。

理論: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{split} \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, \end{split}\]

を、その最低固有対について解きます。固有値 $\mu<0$ が自動的なレベルシフトとして働くので、 $\mathbf H-\mu\mathbf I$ は常に正定値になり、$\mathbf H$ 自体が不定でもステップは常に降下します。最低固有対だけが必要で、拡大行列の積をその場で評価する Davidson 反復で見つけます（$\mathbf H$ は決して作りません）。
TRAH（信頼領域拡大 Hessian）は明示的な信頼半径を加えます:同じ拡大問題を、$\mathbf H$ を試行ステップ長 $\lambda\ge1$ でスケールして解き、結果の $\|\kappa\|$ が目標半径に合うよう選びます — ステップが受理されれば半径を広げ（予測と実際のエネルギー低下を比べる ρ テスト）、そうでなければ縮めます。これが trah を最も難しい開殻系（軌道的にほぼ縮退した ROHF ラジカルなど）で頑健にするもので、固定ステップだと誤った電子状態へ行き過ぎかねない場面です。

        flowchart TD
    G["初期推定密度 D₀（sad, gwh, ...）"] --> F["Fock F を構築"]
    F --> E["誤差 e = FDS − SDF"]
    E --> Q{"‖e‖ としきい値"}
    Q -->|"大きい（遠方）"| ED["EDIIS / ADIIS<br/>エネルギーモデルのステップ"]
    Q -->|"小さい（近傍）"| CD["CDIIS<br/>Pulay 外挿"]
    Q -->|"停滞"| SO["2 次の仕上げ<br/>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 収束理論にあります。

収束:SCF を導く#

上の理論を手にした今、実務的な制御は 2 つのつまみです:戦略選択子 algorithm=（上のどの描像を使うか）と、どれかを安定させる少数の安定化引数。初期推定と同様、収束制御は経路だけを変え、収束後の答えは変えません — 以下のどの設定も同じエネルギーに達します。

既定の 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 ベクトル積です:

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 キーはリファレンスにあります。

実行結果を読む・診断する#

.run() の後、scf アクセサが結果を公開します;run(log=...)（クイックスタート参照）は収束を見るための実況サイクル表を流します:

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）:

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") とともに）再実行することです。

総合例と演習#

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

練習 4

メチルラジカル CH₃ を ref="r", spin=2 で実行します。qc-rs はどのソルバを使い、結果はスピン純粋ですか。 ref="u" ならどうですか。
ある UHF 計算が二重項に対し $\langle S^2\rangle = 1.30$ を与えました。これは許容できますか。何を意味し、何を試しますか。
伸ばした結合の RHF SCF が収束（converged=True）しましたが、基底状態ではないと疑っています。これを確認する 1 つのキーワードは何で、どんな結果が疑いを裏づけますか。

解答練習 4

開殻での ref="r" は ROHF を選び（参照表を参照）、スピン純粋（$\langle S^2\rangle$ 厳密）です。 ref="u" なら UHF — より低いエネルギーですが、わずかにスピン汚染します。
$\langle S^2\rangle = 1.30$ は厳密な二重項値 $0.75$ に対し重い汚染です — UHF 行列式が高スピン状態と激しく混ざっており、エネルギーも物性も疑わしいです。ROHF（スピン純粋）を試すか、状態が本当に多参照か（相関/多参照法が必要か）を確認します。
stability=True。負の external 固有値（external_stable: False）が、より広い参照法のより低い解の存在を裏づけます — UHF で（spin_break="mix" とともに）再実行してそこへ到達します。

SCF は参照波動関数とそのエネルギーを与えます。それが外す相関を回収するため、次章はポスト SCF 法 — RI-MP2 系列 — を加えます。