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 を導く#

SCF ループは反復的なので、数百ではなく数サイクルで収束させる加速戦略が必要です。初期推定と同様、収束制御は経路だけを変え、収束後の答えは変えません — 以下のどの設定も同じエネルギーに達します。

既定の algorithm="auto" は SAD 推定から DIIS（direct inversion of the iterative subspace — 過去のサイクルから次の Fock 行列を外挿）を走らせ、平衡付近のたいていの分子を 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`

2 次戦略（qc/trah）は厳密な軌道 Hessian を作り Newton ステップを取ります — 1 サイクルあたりの仕事は多いが、難しい系ではサイクル数がずっと少なくなります:

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 系列 — を加えます。