🤞
Julia言語による対称直交化と正準直交化の実装

2025/01/06に公開
A. ザボ, N.S.オストランド著, 大野公男, 阪井健男, 望月祐志訳『新しい量子化学 上　電子構造の理論入門』(東京大学出版会, 1987) の3.4.5と3.5.2で解説されている対称直交化と正準直交化をJuliaで実装します.

 一般固有値問題から標準固有値問題へ量子力学に関わる分野では, よく一般固有値問題
HC=ESC\pmb{H} \pmb{C} = E \pmb{S} \pmb{C}HC=ESCに遭遇します. これを標準固有値問題
HC=EC\pmb{H} \pmb{C} = E \pmb{C}HC=ECの形に変形する手続きとして対称直交化や正準直交化が知られています. 例えば, Hartree-Fock法ではRoothaan方程式
FC=SCε(3.159)\pmb{F} \pmb{C} = \pmb{S} \pmb{C} \pmb{\varepsilon} \tag{3.159}FC=SCε(3.159)を前述の手法によって直交化し,
F′C′=C′ε(3.178)\pmb{F}' \pmb{C}' = \pmb{C}' \pmb{\varepsilon} \tag{3.178}F′C′=C′ε(3.178)に変形してから解きます. 直交化の具体的な説明はA. ザボ, N.S.オストランド著, 大野公男, 阪井健男, 望月祐志訳『新しい量子化学 上　電子構造の理論入門』(東京大学出版会, 1987) の3.4.5や3.5.2にあります. 対称直交化と正準直交化では, 重なり行列 S\pmb{S}S を単位行列にするための行列 X\pmb{X}X が異なります.
Xsymmetric=S−1/2=Us−1/2U†(3.167)
\pmb{X}_\mathrm{symmetric}
= \pmb{S}^{-1/2}
= \pmb{U} \pmb{s}^{-1/2} \pmb{U}^{\dagger}
\tag{3.167}
Xsymmetric​=S−1/2=Us−1/2U†(3.167)Xcanonical=Us−1/2(3.169)
\pmb{X}_\mathrm{canonical}
= \pmb{U} \pmb{s}^{-1/2}
\tag{3.169}
Xcanonical​=Us−1/2(3.169)いずれの場合でも
X†SX=1(3.165)
\pmb{X}^{\dagger} \pmb{S} \pmb{X} = \pmb{1} \tag{3.165}
X†SX=1(3.165)となることを実際に計算してみましょう.

 重なり行列の対角化3.5.2節に下記の具体的な計算例があります. 重なり行列 S の固有値 s と固有ベクトルを並べた行列 U を使って X を計算し, これらの数値と一致するか確かめます.
S=(1.00.45080.45081.0)(3.251)
\pmb{S} = \left(\begin{array}{ll}
    1.0 & 0.4508 \\
    0.4508 & 1.0
\end{array}\right)
\tag{3.251}
S=(1.00.4508​0.45081.0​)(3.251)重なり行列の定義
julia> S = [1.0  0.4508; 0.4508  1.0]
2×2 Matrix{Float64}:
 1.0     0.4508
 0.4508  1.0
s, U = eigen(S) で重なり行列を対角化してsとUが得られます. 実際に教科書の値と一致していることが確かめられました.
U=([2]−1/2[2]−1/2[2]−1/2−[2]−1/2)(3.259)
\pmb{U}
= \left(\begin{array}{lr}
    {[2]^{-1/2}} & {[2]^{-1/2}} \\
    {[2]^{-1/2}} & -[2]^{-1/2}
\end{array}\right)
\tag{3.259}
U=([2]−1/2[2]−1/2​[2]−1/2−[2]−1/2​)(3.259)s−1/2=(s1−1/200s2−1/2)=(0.83020.00.01.3493)(3.260)
s^{-1 / 2}
= \left(\begin{array}{cc}
    s_1^{-1/2} & 0 \\
    0 & s_2^{-1/2}
\end{array}\right)
= \left(\begin{array}{ll}
    0.8302 & 0.0 \\
    0.0 & 1.3493
\end{array}\right)
\tag{3.260}
s−1/2=(s1−1/2​0​0s2−1/2​​)=(0.83020.0​0.01.3493​)(3.260)パッケージの読み込み, 重なり行列の対角化, 並び替え
julia> using LinearAlgebra
julia> s, U = eigen(S)
julia> s = sort(s, rev=true)
julia> U = U[:, sortperm(s)]
2×2 Matrix{Float64}:
 0.707107  -0.707107
 0.707107   0.707107

julia> diagm(s .^ (-1/2))
2×2 Matrix{Float64}:
 0.830226  0.0
 0.0       1.34938

 対称直交化Juliaでは S ^ (-1/2) が直ちに計算できてしまいます. これは U * diagm(s .^ (-1/2)) * U' と一致しており, 多少の誤差が入るものの X' * S * X が実際に単位行列になっていることも確かめられます.
Xsymmetric=S−1/2=Us−1/2U†=(1.0898−0.2596−0.25961.0898)(3.261)
\pmb{X}_\mathrm{symmetric}
= \pmb{S}^{-1/2}
= \pmb{U}
  \pmb{s}^{-1/2}
  \pmb{U}^{\dagger}
= \left(\begin{array}{rr}
    1.0898 & -0.2596 \\
    -0.2596 & 1.0898
\end{array}\right)
\tag{3.261}
Xsymmetric​=S−1/2=Us−1/2U†=(1.0898−0.2596​−0.25961.0898​)(3.261)対称直交化
julia> X = S ^ (-1/2)
2×2 Symmetric{Float64, Matrix{Float64}}:
  1.0898    -0.259578
 -0.259578   1.0898

julia> X = U * diagm(s .^ (-1/2)) * U'
2×2 Matrix{Float64}:
  1.0898    -0.259578
 -0.259578   1.0898

julia> X' * S * X
2×2 Matrix{Float64}:
  1.0          -7.50288e-17
 -1.29604e-16   1.0

 正準直交化diagm(s .^ (-1/2)) は単に s の要素の平方根の逆数 si−1/2s_i^{-1/2}si−1/2​ を並べた行列です. X' * S * X が実際に単位行列になっていることが確かめられます.
Xcanonical=Us−1/2=(0.58710.95410.5871−0.9541)(3.262)
\pmb{X}_\mathrm{canonical}
= \pmb{U}
  \pmb{s}^{-1/2}
= \left(\begin{array}{rr}
    0.5871 & 0.9541 \\
    0.5871 & -0.9541
\end{array}\right)
\tag{3.262}
Xcanonical​=Us−1/2=(0.58710.5871​0.9541−0.9541​)(3.262)正準直交化
julia> X = U * diagm(s .^ (-1/2))
2×2 Matrix{Float64}:
 0.587058  -0.954157
 0.587058   0.954157

julia> X' * S * X
2×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0

 まとめ対称直交化と正準直交化をJuliaで実装し, どちらの手法でも X†SX=1\pmb{X}^{\dagger} \pmb{S} \pmb{X} = \pmb{1}X†SX=1 となっていることを確かめました. どちらを選んでも良いように思われますが, 正準直交化は基底の線形従属性が高い場合への対策が可能であるため, 量子化学計算で実際に用いられています.

 余談直観的なイメージはA. ザボ, N.S.オストランド著, 大野公男, 阪井健男, 望月祐志訳『新しい量子化学 上　電子構造の理論入門』(東京大学出版会, 1987)の図3.7を参照してください. 中井 浩巳 『手で解く量子化学Ⅰ』 (丸善出版, 2022) の4.3にも解説があり, 演習問題も用意されています. 将来的にはTwoBody.jlのRayleigh–Ritz法への正準直交化の導入を計画しています. LAPACKにおける一般固有値問題のソルバdsygvなどは上記の手順とは異なり, dpotrfでCholesky分解を行ってdsygstで標準固有値問題に帰着してからdsyevで解いています.
https://gist.github.com/ohno/d8887dbae0fdd07eddb1b4d3e1c92e6b
Discussion

ログインするとコメントできます