書記が物理やるだけ#417 ハートリー-フォック法の計算-2
ハートリー-フォック法のいくつかの方法を比較してみる。
問題
解答
制限ハートリー＝フォック法（RHF）は、全ての分子軌道（MO）において α スピンと β スピンの電子が同一の空間軌道を共有することを前提とする。このため、適用範囲は電子が全て対になっている閉殻系に限定される。
非制限ハートリー＝フォック法（UHF）は、α スピンと β スピンの電子がそれぞれ独立した異なる空間軌道を持つことを許容する。これにより、閉殻系に加えて、ラジカルや三重項などの開殻系にも適用可能である。
制限開殻ハートリー＝フォック法（ROHF）では、閉殻部分のMOはRHFと同様に制限するが、開殻部分のMOについては、電子が片方のスピンのみで占有する共通の空間軌道を用いる。これにより、ROHFはスピン純粋性を保ったまま開殻系を記述できる。
RHFは、計算コストが最も低く、閉殻系に対してはスピン純粋な波動関数を与えることから、安定な閉殻分子の計算において最も効率的な手法として選択される。
一方、開殻分子を扱う際には、計算が比較的容易でスピン分極効果を取り込めるUHFが広く用いられる。しかし、UHFではスピン汚染が発生する可能性があるため、そのチェックと評価が必須となる。
スピン汚染とは、UHF法によって得られた波動関数が、スピン演算子 S^2 の固有関数ではないために、本来求めたい多重度（例：二重項や三重項）以外の高次スピン状態（例：四重項、五重項など）が混入してしまう現象を指す。この汚染が発生する原因は、UHFが α 軌道と β 軌道を独立させることによって、波動関数のスピン対称性が破綻することにある。結果として、計算されたエネルギーやスピン密度分布といった物性値の信頼性が低下する可能性がある。
スピン汚染を厳密に避け、より信頼性の高い結果を得たい場合には、ROHFが選択される。
水分子カチオンの計算を以下に示す。
import numpy as np
from pyscf import gto, scf
# --- 分子の定義（水分子カチオン: H2O+） ---
# H2O+ は開殻系 (9電子系、スピン多重度 2)
# 電子数 9 (5α, 4β) -> Nα ≠ Nβ
mol = gto.Mole()
mol.atom = """
O 0.000000 0.000000 0.117300
H 0.000000 0.757200 -0.469200
H 0.000000 -0.757200 -0.469200
"""
mol.basis = '6-31g*'
mol.charge = 1 # 総電荷 +1
mol.spin = 1 # スピン多重度 2 (2S+1=2) -> S=1/2 -> mol.spin=1
mol.build()
print(f"--- 分子情報 ---")
print(f"全電子数: {mol.nelectron}")
Nalpha = (mol.nelectron + mol.spin) // 2
Nbeta = (mol.nelectron - mol.spin) // 2
print(f"alpha電子数 (Nα): {Nalpha}")
print(f"beta電子数 (Nβ): {Nbeta}")
print("----------------\n")
# 理論的なスピン角運動量二乗の期待値 (S=1/2 の場合、S^2 = S(S+1) = 0.5 * 1.5 = 0.75)
Ideal_S2 = mol.spin * (mol.spin + 2) / 4
# ----------------------------------------------------
# 1. UHF法 (Unrestricted Hartree-Fock: 非制限)
# ----------------------------------------------------
print("--- 1. UHF法（非制限）の実行 ---")
mf_uhf = scf.UHF(mol)
E_uhf = mf_uhf.kernel()
print(f'UHF 総エネルギー: {E_uhf:.10f} Hartrees')
# UHFの波動関数とエネルギーの確認
E_alpha = mf_uhf.mo_energy[0]
E_beta = mf_uhf.mo_energy[1]
# ★ スピン汚染の確認 (S^2期待値)
S2_uhf = mf_uhf.spin_square()[0]
print(f" スピン汚染 S^2 期待値: {S2_uhf:.6f} (理論値: {Ideal_S2:.2f})")
print(f" α-HOMO ({Nalpha}番目の軌道): {E_alpha[Nalpha-1]:.6f} Hartrees")
print(f" β-HOMO ({Nbeta}番目の軌道): {E_beta[Nbeta-1]:.6f} Hartrees")
# UHFのMulliken電荷解析
populations_uhf, charges_uhf = mf_uhf.mulliken_pop()
print("\n UHF Mulliken atomic charges:")
for i, symbol in enumerate(mol.elements):
print(f" 原子 {i} ({symbol}): {charges_uhf[i]:.4f}")
# ----------------------------------------------------
# 2. ROHF法 (Restricted Open-shell Hartree-Fock: 制限開殻)
# ----------------------------------------------------
print("\n--- 2. ROHF法（制限開殻）の実行 ---")
mf_rohf = scf.ROHF(mol)
E_rohf = mf_rohf.kernel()
print(f'ROHF 総エネルギー: {E_rohf:.10f} Hartrees')
# ROHFの波動関数とエネルギーの確認
E_mo = mf_rohf.mo_energy
Nocc_open = Nalpha # 開殻部分を含めた占有軌道数 (αの数)
# ★ スピン純粋性の確認 (ROHFは定義上スピン純粋: 理論値に近いことを確認)
S2_rohf = mf_rohf.spin_square()[0]
print(f" スピン汚染 S^2 期待値: {S2_rohf:.6f} (理論値: {Ideal_S2:.2f})")
# ROHFではαとβは同じ空間軌道を持つため、エネルギーはMO全体のリストから取得
print(f" ROHF HOMO-α/開殻軌道 (軌道 {Nocc_open-1}): {E_mo[Nocc_open-1]:.6f} Hartrees")
# ROHFの場合、HOMO-βは閉殻部分の最も高い軌道エネルギーを指す
print(f" ROHF 閉殻部 HOMO (軌道 {Nbeta-1}): {E_mo[Nbeta-1]:.6f} Hartrees")
# ROHFのMulliken電荷解析
populations_rohf, charges_rohf = mf_rohf.mulliken_pop()
print("\n ROHF Mulliken atomic charges:")
for i, symbol in enumerate(mol.elements):
print(f" 原子 {i} ({symbol}): {charges_rohf[i]:.4f}")
# ----------------------------------------------------
# 3. RHF法 (比較用: H2O+には不適切な手法)
# ----------------------------------------------------
print("\n--- 3. RHF法（制限：不適切）の実行 ---")
# 開殻系に対してRHFを強制すると、通常は高エネルギーまたは計算が破綻する
mf_rhf = scf.RHF(mol)
E_rhf = mf_rhf.kernel()
print(f'RHF 総エネルギー: {E_rhf:.10f} Hartrees')
print("\n--- 全ての計算が完了しました ---")結果は以下の通りである。
--- 分子情報 ---
全電子数: 9
alpha電子数 (Nα): 5
beta電子数 (Nβ): 4
----------------
--- 1. UHF法（非制限）の実行 ---
converged SCF energy = -75.610498235251 <S^2> = 0.7565988 2S+1 = 2.006588
UHF 総エネルギー: -75.6104982353 Hartrees
スピン汚染 S^2 期待値: 0.756599 (理論値: 0.75)
α-HOMO (5番目の軌道): -1.100760 Hartrees
β-HOMO (4番目の軌道): -1.047942 Hartrees
** Mulliken pop alpha | beta **
pop of 0 O 1s 0.99900 | 0.99900
pop of 0 O 2s 0.51224 | 0.47068
pop of 0 O 3s 0.46678 | 0.47874
pop of 0 O 2px 0.68312 | 0.00000
pop of 0 O 2py 0.50519 | 0.46571
pop of 0 O 2pz 0.57267 | 0.54507
pop of 0 O 3px 0.31549 | 0.00000
pop of 0 O 3py 0.26692 | 0.26501
pop of 0 O 3pz 0.30687 | 0.32014
pop of 0 O 3dxy 0.00000 | 0.00000
pop of 0 O 3dyz 0.01138 | 0.01186
pop of 0 O 3dz^2 0.00223 | 0.00181
pop of 0 O 3dxz 0.00138 | 0.00000
pop of 0 O 3dx2-y2 0.00107 | 0.00352
pop of 1 H 1s 0.19088 | 0.22134
pop of 1 H 2s -0.01305 | -0.00211
pop of 2 H 1s 0.19088 | 0.22134
pop of 2 H 2s -0.01305 | -0.00211
In total 5.00000 | 4.00000
** Mulliken atomic charges ( Nelec_alpha | Nelec_beta ) **
charge of 0O = -0.20588 ( 4.64435 3.56153 )
charge of 1H = 0.60294 ( 0.17782 0.21924 )
charge of 2H = 0.60294 ( 0.17782 0.21924 )
UHF Mulliken atomic charges:
原子 0 (O): -0.2059
原子 1 (H): 0.6029
原子 2 (H): 0.6029
--- 2. ROHF法（制限開殻）の実行 ---
converged SCF energy = -75.6060375138485
ROHF 総エネルギー: -75.6060375138 Hartrees
スピン汚染 S^2 期待値: 0.750000 (理論値: 0.75)
ROHF HOMO-α/開殻軌道 (軌道 4): -0.702718 Hartrees
ROHF 閉殻部 HOMO (軌道 3): -1.074610 Hartrees
** Mulliken pop **
pop of 0 O 1s 1.99799
pop of 0 O 2s 0.98468
pop of 0 O 3s 0.94723
pop of 0 O 2px 0.67958
pop of 0 O 2py 0.97157
pop of 0 O 2pz 1.11807
pop of 0 O 3px 0.31905
pop of 0 O 3py 0.53253
pop of 0 O 3pz 0.62671
pop of 0 O 3dxy 0.00000
pop of 0 O 3dyz 0.02325
pop of 0 O 3dz^2 0.00379
pop of 0 O 3dxz 0.00137
pop of 0 O 3dx2-y2 0.00370
pop of 1 H 1s 0.41130
pop of 1 H 2s -0.01606
pop of 2 H 1s 0.41130
pop of 2 H 2s -0.01606
** Mulliken atomic charges **
charge of 0O = -0.20952
charge of 1H = 0.60476
charge of 2H = 0.60476
ROHF Mulliken atomic charges:
原子 0 (O): -0.2095
原子 1 (H): 0.6048
原子 2 (H): 0.6048
--- 3. RHF法（制限：不適切）の実行 ---
converged SCF energy = -75.6060375138485
RHF 総エネルギー: -75.6060375138 Hartrees
--- 全ての計算が完了しました ---(1)総エネルギーとスピン分極効果
総エネルギーの比較において、UHF法が最も低い総エネルギー（−75.6105 Hartrees）を与えている。これは、UHF法が α 軌道と β 軌道を独立させることで、スピン分極効果という形で電子相関の一部を間接的に取り込み、ROHFよりも柔軟な波動関数を持つためである。これは、より低いエネルギーを与えるという変分原理の観点から理論的に期待される結果である。
(2)スピン純粋性とスピン汚染
ROHF法は、その定義通り ⟨S^2⟩=0.750000 を与え、完璧なスピン純粋性を保っている。これが開殻計算におけるROHFの最大の優位性である。
一方で、UHF法は ⟨S^2⟩=0.756599 であり、理論値 0.75 から僅かに逸脱している。これは、目的とする二重項状態以外の高次スピン状態（この場合は四重項など）が波動関数に混入していることを示しており、スピン汚染が発生している証拠である。
(3)軌道エネルギーと制限性の違い
UHF法では、α 軌道と β 軌道が独立して最適化されるため、それぞれのHOMO（最高被占有分子軌道）エネルギーが異なり、非制限的な性質を示している。
一方、ROHF法では、開殻軌道（α のみで占有）と、それより内側の閉殻軌道（α と β で占有）のエネルギーが異なる。これは、閉殻部分には制限を課しつつ、開殻部分の軌道を特別に扱うというROHFの制限開殻の定義を反映した結果である。
本記事のもくじはこちら：
