🎁

Julia言語によるHartree-Fock法の実装

2022/12/24に公開

Julia

JuliaLang Advent Calendar 2022 24日目です🎄
今日は個人で開発しているパッケージを公開する予定でしたが, 先日のゼミでHartree-Fock法を扱ったので急遽Juliaで実装することにしました.

ザボ・オストランド著『新しい量子化学上』(1987)の付録Bを参考に, Julia言語でHartree-Fock法のプログラム（約140行）を実装します. 参考にした付録Bのコード（約570行）は Computational Chemistry List (CCL) で配布されています. これは古代言語FORTRANで書かれたコードなので解読が難しいですが, ほぼ同じ変数名・出力のプログラムを作ったので, 対応関係は察せられると思います. 元々はHeH⁺をHF/STO-3Gで解くコードですが, 電荷や軌道指数を変えれば H₂, HeH⁺, LiH²⁺ などの2核2電子系が扱えます. 今回はH₂とHeH⁺をJuliaで計算しました.

はじめに

こちらの記事の通り, 計算化学に関わる研究室に配属された学生はザボ・オストランド著『新しい量子化学上』(1987)（通称「ザボ」）を勉強することになります. 筆者も学部3年生の時にみっちりザボを勉強しました. 学部3年生の当時は「Hartree-Fockなんて時代遅れだ」なんて思っていましたが, その~~過酷な~~丁寧なゼミの~~後遺症~~おかげで, Hartree-Fock法の重要性が~~解らせられ~~理解できたと思っています. Hartree-Fock法（HF法）は第一原理分子軌道法の一種あり, 数ある電子状態理論の中でも特別な存在です. その理由は

電子状態計算の精度はHF法をベンチマークとして評価される.（相関エネルギー（BSCEの方）の取り込み率はHF法を0%, full CIまたは厳密解を100%とする）
密度汎関数法（DFT）を勉強する際に重要.
post-HF法（CI, MP, CC等）はHFをベースとして拡張されたもの.
量子モンテカルロ法（VMC, FN-DMC等のab initio QMC）でも最初にHF法で計算して, その結果を利用してサンプリングする.
古典力場を用いない第一原理分子動力学法（abinitio MD）でも, 何らかの電子状態計算を行っている.

などです. この分野に居る以上, 直接的にはHF法を使わないにしても, 間接的にはHF法に関わることになります. とはいえ, HF法の汎用コードは

等のパッケージに含まれており, 有償・無償で手に入りますので, 労力や信頼性の観点からも, 自分で実装することはお勧めしません. Juliaで新たに実装した目的は

新規手法の開発のため
既存のパッケージで計算できないような問題については, パッケージを拡張したり, 1から開発することになることがあります. むしろ, それが計算系の研究者に求められていることだと思います. そこで, まずは簡単に動かせる, 小さなプログラムがあるとプロトタイプを開発しやすくなります.
自身の学習のため
上記と同じ理由から, 巨大なパッケージをいきなり拡張しようとするよりは, 簡単なプログラムで大まかな計算の流れを理解しておきたいと思うはずです.
教育のため
付録Bのコードは古代言語FORTRANで書かれています. 21世紀にもなってFORTRANを読める・勉強したい学生はいません.

などです. 数あるプログラミング言語の中からJuliaを選んだ理由は以下の通りです.

数値計算の分野に限って言えば, Python等よりもコードがシンプルになる.
ギリシャ文字等が使えるので, 数式とコードの対応がわかりやすい.
精度保証付き演算が利用できる.
四倍精度, 八倍精度, 任意精度演算が利用できる.

JuliaでもFermi.jl等のパッケージがありますが, 純粋なネイティブ実装ではないので, 精度保証付き演算や任意精度演算は使えません. 重なり行列の対角化が悪条件なので, もし汎用コードをJuliaでネイティブ実装できれば応用が期待できます.

Schrödinger方程式とBorn-Oppenheimer近似

分子は電子と原子核の多体系です. ハミルトニアンのポテンシャル項を見れば明らかですが, 電子と原子核はどちらも質点・点電荷として扱われることに注意してください. よくある勘違いですが, いわゆる電子雲とは, 電子が存在する位置の確率密度関数について言及したものであって, 電子の大きさに関する概念ではありません.

分子系のSchrödinger方程式を解くとき, 電子と原子核の質量が大きく異なることから, Born-Oppenheimer近似（BO近似）がよく用いられます. BO近似は「原子核を止めて電子を扱うこと」などと説明されることがありますが, これはあくまで電子から見た話であって, 原子核から見ればまた別の説明が与えられます. そもそも, 核を止めるとか, 電子を平均化して扱うというのは, 近似をした後の結果に与えられる解釈に過ぎず, 近似そのものの説明にはなっていません. BO近似の本質は電子と原子核の運動を分離することです.

分離とは何なのか, ハミルトニアンから見ていきましょう. ザボの2.1節にはハミルトニアンが3つ登場します.

\begin{align} \hat{H} = - \sum_{i=1}^N \frac{1}{2} \nabla_i^2 - \sum_{i=1}^N \frac{1}{2 M_A} \nabla_A^2 - \sum_{i=1}^N \sum_{A=1}^M \frac{Z_A}{r_{i A}} + \sum_{i=1}^N \sum_{j>i}^N \frac{1}{r_{i j}} + \sum_{A=1}^M \sum_{B>A}^M \frac{Z_A Z_B}{R_{i j}} \tag{2.2} \end{align}

\begin{align} \hat{H}_\mathrm{elec} = - \sum_{i=1}^N \frac{1}{2} \nabla_i^2 - \sum_{i=1}^N \sum_{A=1}^M \frac{Z_A}{r_{i A}} + \sum_{i=1}^N \sum_{j>i}^N \frac{1}{r_{i j}} \tag{2.10} \end{align}

\begin{align} \hat{H}_\mathrm{nucl} = - \sum_{i=1}^N \frac{1}{2 M_A} \nabla_A^2 + \sum_{A=1}^M \sum_{B>A}^M \frac{Z_A Z_B}{R_{i j}} + E_{\mathrm{elec}}\left(\left\{\pmb{R}_A\right\}\right) \tag{2.15} \end{align}

それぞれに対応するSchrödinger方程式も3つ登場します.

\begin{align} \hat{H} \varPsi = E \varPsi \tag{2.1} \end{align}

\begin{align} \hat{H}_\mathrm{elec} \varPsi_\mathrm{elec} = E_\mathrm{elec} \varPsi_\mathrm{elec} \tag{2.11} \end{align}

\begin{align} \hat{H}_\mathrm{nucl} \varPsi_\mathrm{nucl} = E_\mathrm{nucl} \varPsi_\mathrm{nucl} \tag{2.16} \end{align}

分子の全ハミルトニアンは(2.2)であり, 求めたいエネルギーは式(2.1)の $E$ です. しかし, 普通の計算パッケージでは式(2.1)は解きません. 代わりに, 電子の運動について式(2.11)を, 原子核の運動について式(2.16)を解きます. このように, BO近似では電子と原子核の運動を別の方程式に分けて扱います. このうち, 式(2.16)の $E_\mathrm{nucl}$ が式(2.1)の $E$ に相当するものであり, $E \simeq E_\mathrm{nucl}$ です. しかし, 式(2.15) $\hat{H}_\mathrm{nucl}$ のポテンシャル部分 $\sum_{A=1}^M \sum_{B>A}^M \frac{Z_A Z_B}{R_{i j}} + E_{\mathrm{elec}}\left(\left\{\pmb{R}_A\right\}\right)$ には $E_{\mathrm{elec}}\left(\left\{\pmb{R}_A\right\}\right)$ が入ってしまっていますので, 式(2.16)を解いて $E_\mathrm{nucl}$ を求めるには, 先に式(2.10)を（核間距離を変えて何度も）解かなければいけないことがわかります. この, 原子核が感じるポテンシャル項を

\begin{align} E_{\mathrm{tot}}\left(\left\{\pmb{R}_A\right\}\right) := \sum_{A=1}^M \sum_{B>A}^M \frac{Z_A Z_B}{R_{i j}} + E_{\mathrm{elec}}\left(\left\{\pmb{R}_A\right\}\right) \tag{2.14} \end{align}

という記号で表し, 核配置 $\{\pmb{R}_i\}$ の関数と考えてポテンシャルエネルギー曲線(Potential Energy Curve, PEC)とかポテンシャルエネルギー曲面(Potential Energy Surface, PES)と呼んだりします. 実用上はPECを全て求めるというよりも, 極小点付近を調和振動子に近似したり, 非調和補正を加えたりして $E_\mathrm{nucl}$ を求めます. 筆者の所属する研究室では, BO近似をそもそも用いない第一原理多成分分子軌道法や, PES上で核の量子効果・温度効果を扱う第一原理経路積分分子動力学法を開発・応用しています.

※ $\{\pmb{R}_i\}$ は組 $(\pmb{R}_1,\pmb{R}_2,\pmb{R}_3,\cdots)$ を意味します. 集合の記号と紛らわしい問題はザボ以外でもよくあるそうです.

Hartree-Fock法

いろいろ記号が登場して混乱したと思いますが, $E$ や $E_\mathrm{nucl}$ は核配置 $\pmb{R}_i$ に依存せず, $E_\mathrm{elec}(\{\pmb{R}_i\})$ や $E_\mathrm{tot}(\{\pmb{R}_i\})$ は核配置 $\{\pmb{R}\}$ に依存するという点を踏まえれば, 整理できると思います. HF法は与えられた核配置1点 $\left\{\pmb{R}_A\right\}$ に対して $E_{\mathrm{elec}}(\{\pmb{R}_A\})$ を近似的に求めます.

HF法では試行波動関数に1つのSlater行列式(3.1)を用い, 束縛条件(3.3)の下でエネルギーの期待値(3.2)を最小化するようなスピン軌道を定める方程式, Hartree-Fock微分方程式

\begin{align} h(1) \chi_a(1) + \sum_{b \neq a}\left[\int d x_2\left|\chi_b(2)\right|^2 r_{12}^{-1}\right] \chi_a(1) - \sum_{b \neq a}\left[\int d x_2 \chi_b^*(2) \chi_a(2) r_{12}^{-1}\right] \chi_b(1) = \varepsilon_a \chi_a(1) \tag{3.4} \end{align}

を得ます. 導出は3.1節と3.2節を読んでください. 実際には, ユニタリ変換によって一意性を与えられた正準Hartree-Fock方程式(3.69)が解かれます. これらは微分方程式ですので, コンピュータで解くために行列方程式に帰着します. その結果がRoothaan方程式

\begin{align} \pmb{F} \pmb{C} = \pmb{S} \pmb{C} \pmb{\varepsilon} \tag{3.159} \end{align}

です. 実際には, これを直交化して(3.178)を解きます. 導出・直交化については3.4節を読んでください. このFock行列 $\pmb{F}$ が解である $\pmb{C}$ に依存してしまっている＝自己無撞着（じこむどうちゃく, self-consistent）であることから, ただ固有値問題を解くだけでは所望の解が得られません. そもそも, 解きようがないようにも思われます. しかし, 適当な初期値を与えて繰り返し解いていくと, 得られる $\pmb{C}$ と $\pmb{\varepsilon}$ が一定の値に収束します. このアルゴリズムをSCFの手続きといいます. 線形変分法のように一般化固有値問題 $\pmb{H}\pmb{C}=E\pmb{S}\pmb{C}$ を一度だけ解けば終わりであることとは対照的です.

HF法の実装は何が難しいのか？

HF法というと, SCFの手続きが強調されますが, SCFの実装は簡単です. 特に, Juliaだと10行程度で書けてしまいます. どちらかといえば, 基底関数の扱いと分子積分の用意が難しいところです. 具体的には

基底関数系データファイルのパーサーを自作する
基底関数系の適切なデータ構造を設計する
分子積分の公式の導出
分子積分の実装

などに労力のほとんどが費やされます. 実際, ザボでもかなりのページ数を割いて解説されている部分だと思われます. この辺りはGaussianBasis.jl（libcintのラッパー）で既に実装されており, Fermi.jlで使われています. しかし, 2022年12月現在, ビルドでエラーが出ますので, 今回はザボを参考に実装しました.

実装

ザボ p.159 3.4.6 を下表にまとめました. これを実装していきます.

手順	説明	数式
1	パラメータ入力	$\{\pmb{R}_A\}, \{Z_A\}, N, \{\phi_\mu\}$
2	分子積分の計算	$S_{\mu\nu}, H_{\mu\nu}^\mathrm{core}, (\mu\nu\mid\lambda\sigma)$
3	式(3.167), $\pmb{X}$ の計算	$\pmb{X} = \pmb{U} \pmb{s}^{-1/2}$
4	密度行列 $\pmb{P}$ の初期値設定	$\boldsymbol{P}=\left(\begin{array}{cc} 0 & \cdots \\ \vdots & \ddots \end{array}\right)$
5	式(3.154), $\pmb{G}$ の計算	$G_{\mu\nu} = \sum_{\lambda,\sigma} P_{\lambda\sigma} \left[ (\mu\nu\mid\sigma\lambda) - \frac{1}{2}(\mu\lambda\mid\sigma\nu) \right]$
6	式(3.154), $\pmb{F}$ の計算	$\pmb{F} = \pmb{H}^\mathrm{core} + \pmb{G}$
7	式(3.177), $\pmb{F}'$ の計算（正準直交化）	$\pmb{F}' = \pmb{X}'\pmb{F}\pmb{X}$
8	式(3.178)を解いて $\pmb{C}'$ と $\pmb{\varepsilon}$ を得る.	$\pmb{F}'\pmb{C}' = \pmb{C}'\pmb{\varepsilon}$
9	式(3.174), $\pmb{C}$ の計算	$\pmb{C}=\pmb{X}\pmb{C}'$
10	式(3.145), $\pmb{P}$ の計算	$P_{\mu\nu} = 2 \sum_a^{N/2} C_{\mu a} C_{\nu a}^*$
11	(5)に戻る. $\delta<10^{-4}$ なら終了.	$\delta = \frac{1}{4} \sum_{\mu,v} (P_{\mu\nu}^\mathrm{new} - P_{\mu\nu}^\mathrm{old})^2$
12	式(3.184), $E_0$ 等の計算	$E_0 = \frac{1}{2} \sum_{\mu,v} P_{\mu\nu}\left(H_{\mu \nu}^\mathrm{core} + F_{\mu\nu} \right)$

# Shuhei Ohono https://github.com/ohno
# Licence: MIT

# using Pkg
# Pkg.add("SpecialFunctions")
using SpecialFunctions
using LinearAlgebra

# 手順(1)  パラメータ入力
function HF(; R=[0.0 0.0 0.0; 0.0 0.0 1.4632], Z=[2,1], ζ=[2.0925,1.24], α=[0.109818 0.405771 2.22766], c=[0.444635 0.535328 0.154329],  A=[α;α], D=[c;c])
  A[1,:] *= ζ[1]^2
  A[2,:] *= ζ[2]^2
  D[1,:] .*= ((2*A[1,:]/π).^(3/4))
  D[2,:] .*= ((2*A[2,:]/π).^(3/4))
  println("STO-3G FOR ATOMIC NUMBERS $(Z[1]) AND $(Z[2])")
    
  # 手順(2)  分子積分の計算
  F0(x) = x<1e-6 ? 1-x/3 : sqrt(π/x)*erf(sqrt(x))/2
  integral_S(α, β, RA, RB) = (π/(α+β))^(3//2) * exp(-α*β*sum((RA-RB).^2)/(α+β))
  integral_T(α, β, RA, RB) = α*β/(α+β) * (3-2*α*β*sum((RA-RB).^2)/(α+β)) * (π/(α+β))^(3//2) * exp(-α*β*sum((RA-RB).^2)/(α+β))
  integral_V(α, β, RA, RB, RC, ZC) = -2*π/(α+β) * F0((α+β)*sum(((α*RA+β*RB)/(α+β)-RC).^2)) * exp(-α*β*sum((RA-RB).^2)/(α+β)) * ZC
  integral_TwoE(α, β, γ, δ, RA, RB, RC, RD) = 2*(π^(5//2))/((α+β)*(γ+δ)*sqrt(α+β+γ+δ)) * exp(-α*β/(α+β)*sum((RA-RB).^2)-γ*δ/(γ+δ)*sum((RC-RD).^2)) * F0((α+β)*(γ+δ)/(α+β+γ+δ)*sum(((α*RA+β*RB)/(α+β)-(γ*RC+δ*RD)/(γ+δ)).^2))

  # 一電子積分
  S = zeros(2,2)
  T = zeros(2,2)
  VA = zeros(2,2)
  VB = zeros(2,2)
  for i=1:3, j=1:3
    for μ=1:2, ν=1:2
       S[μ,ν] += D[μ,i] * D[ν,j] * integral_S(A[μ,i], A[ν,j], R[μ,:], R[ν,:])
       T[μ,ν] += D[μ,i] * D[ν,j] * integral_T(A[μ,i], A[ν,j], R[μ,:], R[ν,:])
      VA[μ,ν] += D[μ,i] * D[ν,j] * integral_V(A[μ,i], A[ν,j], R[μ,:], R[ν,:], R[1,:], Z[1])
      VB[μ,ν] += D[μ,i] * D[ν,j] * integral_V(A[μ,i], A[ν,j], R[μ,:], R[ν,:], R[2,:], Z[2])
    end
  end
  S[1,1] = 1
  S[2,2] = 1

  # 式(3.233)など
  Hcore = T + VA + VB

  # 二電子積分
  TwoE = zeros(2,2,2,2)
  for μ=1:2, ν=1:2, σ=1:2, λ=1:2
    for i=1:3, j=1:3, k=1:3, l=1:3
      TwoE[μ,ν,σ,λ] += D[μ,i] * D[ν,j] * D[σ,k] * D[λ,l] * integral_TwoE(A[μ,i], A[ν,j], A[σ,k], A[λ,l], R[μ,:], R[ν,:], R[σ,:], R[λ,:])
    end
  end

  # 手順(3) 式(3.167), Xの計算
  s, U = eigen(S)
  s = sort(s, rev=true)
  U = U[:, sortperm(s)]
  U[:,2] .*= -1
  X = U * (diagm(s) ^ (-1//2))
  X°= transpose(X)
    
  # 手順(4)  密度行列に対する初期値を設定する.
  P = zeros(2,2)
    
  # 行列の出力
  println("\nR\tζ1\tζ2\tS12\tT11")
  println(norm(R[1,:]-R[2,:]), "\t", ζ[1], "\t", ζ[2], "\t", S[1,2], "\t", T[1,1])
  println("T12\tT22\tV11A\tV12A\tV22A")
  println(T[1,2], "\t", T[2,2], "\t", VA[1,1], "\t", VA[1,2], "\t", VA[2,2])
  println("V11B\tV12B\tV22B\tV1111\tV2111")
  println(VB[1,1], "\t", VB[1,2], "\t", VB[2,2], "\t", TwoE[1,1,1,1], "\t", TwoE[2,1,1,1])
  println("V2121\tV2211\tV2221\tV2222")
  println(TwoE[2,1,2,1], "\t", TwoE[2,2,1,1], "\t", TwoE[2,2,2,1], "\t", TwoE[2,2,2,2])  
  println("\nS = $S")
  println("X = $X")
  println("Hcore = $Hcore\n")
  for μ=1:2, ν=1:2, σ=1:2, λ=1:2
    println("($μ,$ν,$σ,$λ)  $(TwoE[μ,ν,σ,λ])")
  end
  println("\nP  = $P")

  # その他の初期値設定
  δ = Inf
  ITER = 0
  EN = 0.0
  G = zeros(2,2)
  F = zeros(2,2)
  ε = zeros(2,2)
  C = zeros(2,2)
  P = zeros(2,2)
  while δ>1e-4  # CRIT = 1e-4
    ITER = ITER+1
    println("\nSTART OF ITERATION NUMBER = $ITER")
    # 手順(5)  式(3.154)  Gの計算
    G = [sum(λ-> sum(σ-> P[λ,σ] * (TwoE[μ,ν,σ,λ] - (1//2)*TwoE[μ,λ,σ,ν]), 1:2), 1:2) for μ=1:2, ν=1:2]
    println("G  = $G")
    # 手順(6)  式(3.154)  Fの計算
    F = Hcore + G
    println("F  = $F")
    # 手順(12)  式(3.184)  エネルギー・マリケン電荷などの計算
    EN = sum(P/2 .* (Hcore + F))
    println("ELECTRONIC ENERGY = $EN")
    # 手順(7)  式(3.177)  F'の計算（正準直交化）
    F° = X° * F * X
    # 手順(8)  式(3.178)  F'C'=C'εを解く
    ε, C° = eigen(F°)
    C°[:,1] *= -1
    # 手順(9)  式(3.174)  Cの計算
    C = X * C°
    # 手順(10)  式(3.145)  Pの計算
    P_old = P
    P = [sum(k->2*C[μ,k]*C[ν,k], 1:1) for μ=1:2, ν=1:2]
    # 手順(11)  収束判定
    δ = sqrt(sum((P-P_old).^2) / 4)
    println("F' = $F°")
    println("C' = $C°")
    println("ε  = $ε")
    println("C  = $C")
    println("P  = $P")
    println("δ (CONVERGENCE OF DENSITY MATRIX) = $δ")
    if (ITER>25)  # MAXIT = 25
      println("NO CONVERGENCE = SCF")
      break
    end
  end

  # 手順(12)  式(3.184)  エネルギー・マリケン電荷などの計算
  EN = sum(P/2 .* (Hcore + F))
  ENT = EN + Z[1] * Z[2] / norm(R[1,:]-R[2,:])
  println("\nELECTRONIC ENERGY = $EN")
  println("TOTAL ENERGY = $ENT\n")
  println("G  = $G")
  println("F  = $F")
  println("ε  = $ε")
  println("C  = $C")
  println("P  = $P")
  println("\nPS = $(P * S)")
  return ENT
end

約570行から約140行程度まで短縮できました. print文やコメント行を省けば70行弱です.

計算

入力

HF()

出力

STO-3G FOR ATOMIC NUMBERS 2 AND 1

R	ζ1	ζ2	S12	T11
1.4632	2.0925	1.24	0.4507704116477877	2.164312561949012
T12	T22	V11A	V12A	V22A
0.1670128658640491	0.7600329435650853	-4.139827204696602	-1.1029125678873686	-1.2652458734233516
V11B	V12B	V22B	V1111	V2111
-0.6772300819326643	-0.4113054702897494	-1.2266154680582526	1.307151607555482	0.43727932526541685
V2121	V2211	V2221	V2222
0.17726712195066152	0.6057033663335998	0.31179457036897396	0.7746083600328786

S = [1.0 0.4507704116477877; 0.4507704116477877 1.0]
X = [0.587064281204403 0.954131072232266; 0.587064281204403 -0.954131072232266]
Hcore = [-2.6527447246802547 -1.347205172313069; -1.3472051723130691 -1.731828397916519]

(1,1,1,1)  1.307151607555482
(1,1,1,2)  0.43727932526541696
(1,1,2,1)  0.43727932526541696
(1,1,2,2)  0.6057033663335999
(1,2,1,1)  0.437279325265417
(1,2,1,2)  0.17726712195066158
(1,2,2,1)  0.17726712195066155
(1,2,2,2)  0.311794570368974
(2,1,1,1)  0.43727932526541685
(2,1,1,2)  0.17726712195066152
(2,1,2,1)  0.17726712195066152
(2,1,2,2)  0.3117945703689739
(2,2,1,1)  0.6057033663335998
(2,2,1,2)  0.311794570368974
(2,2,2,1)  0.31179457036897396
(2,2,2,2)  0.7746083600328786

P  = [0.0 0.0; 0.0 0.0]

START OF ITERATION NUMBER = 1
G  = [0.0 0.0; 0.0 0.0]
F  = [-2.6527447246802547 -1.347205172313069; -1.3472051723130691 -1.731828397916519]
ELECTRONIC ENERGY = -0.0
F' = [-2.4397325070831126 -0.5158386381846956; -0.5158386381846952 -1.5386669016152497]
C' = [0.9104452936060882 0.4136295049322813; 0.413629504932281 -0.9104452936060881]
ε  = [-2.674086050382983, -1.3043133583153792]
C  = [0.9291466749147285 -0.6258570361991943; 0.13983314881885123 1.1115112521952]
P  = [1.7266270870101923 0.259851010535784; 0.259851010535784 0.03910661901718999]
δ (CONVERGENCE OF DENSITY MATRIX) = 0.8828668530136912

START OF ITERATION NUMBER = 2
G  = [1.2623300126571688 0.37400402478104683; 0.3740040247810466 0.9889530230082917]
F  = [-1.3904147120230859 -0.9732011475320221; -0.9732011475320226 -0.7428753749082273]
ELECTRONIC ENERGY = -4.141862876133926
F' = [-1.4060434198411986 -0.362710270326269; -0.36271027032626846 -0.17013631078031075]
C' = [0.9649913879927128 0.26228156835717215; 0.26228156835717176 -0.9649913879927127]
ε  = [-1.5046269095429181, -0.07155282107859129]
C  = [0.8167629696037698 -0.7667521273196227; 0.31626098151699245 1.0747044081211563]
P  = [1.3342034970319372 0.5166205168672434; 0.5166205168672434 0.2000420168601829]
δ (CONVERGENCE OF DENSITY MATRIX) = 0.27917630406864175

START OF ITERATION NUMBER = 3
G  = [1.2013462807310709 0.3038061649826069; 0.30380616498260676 0.9284329239279128]
F  = [-1.4513984439491838 -1.043399007330462; -1.0433990073304624 -0.8033954739886062]
ELECTRONIC ENERGY = -4.226491899128991
F' = [-1.4963056517218072 -0.36296996790009084; -0.36296996790009034 -0.152937775785973]
C' = [0.9694747623841949 0.24519111953761485; 0.24519111953761452 -0.9694747623841949]
ε  = [-1.5881048570658307, -0.06113857044194933]
C  = [0.8030884703111406 -0.7810630461867005; 0.33519953873863295 1.068948942884806]
P  = [1.2899021822933756 0.5383897696292174; 0.5383897696292174 0.22471746154118458]
δ (CONVERGENCE OF DENSITY MATRIX) = 0.029661780077229388

START OF ITERATION NUMBER = 4
G  = [1.1946701848800272 0.29716257964642057; 0.29716257964642045 0.921870489419819]
F  = [-1.4580745398002275 -1.0500425926666483; -1.0500425926666486 -0.8099579084967]
ELECTRONIC ENERGY = -4.227522925343759
F' = [-1.505447587902247 -0.3630336337410123; -0.36303363374101166 -0.15289349529356425]
C' = [0.9698136576776785 0.2438472254954788; 0.2438472254954784 -0.9698136576776782]
ε  = [-1.5967277463620375, -0.06161333683377371]
C  = [0.8020051725696232 -0.7821753489063077; 0.3366807431238953 1.0684833412246901]
P  = [1.2864245936568623 0.5400393949798973; 0.5400393949798973 0.22670784558091672]
δ (CONVERGENCE OF DENSITY MATRIX) = 0.002318284869555417

START OF ITERATION NUMBER = 5
G  = [1.1941478316400083 0.29665158100169653; 0.2966515810016964 0.9213575616263323]
F  = [-1.4585968930402464 -1.0505535913113724; -1.0505535913113726 -0.8104708362901867]
ELECTRONIC ENERGY = -4.2275292681003185
F' = [-1.5061566175000514 -0.3630389132754954; -0.3630389132754951 -0.15290558836387533]
C' = [0.9698390836541676 0.24374608061842734; 0.2437460806184272 -0.9698390836541674]
ε  = [-1.5973978488695526, -0.06166435699437396]
C  = [0.8019235937822431 -0.7822589871650616; 0.33679217527649846 1.0684482223943566]
P  = [1.286162900529256 0.5401631831109376; 0.5401631831109376 0.22685793865495132]
δ (CONVERGENCE OF DENSITY MATRIX) = 0.00017439769686151862

START OF ITERATION NUMBER = 6
G  = [1.1941085339307838 0.29661318951207055; 0.2966131895120704 0.9213189761538769]
F  = [-1.458636190749471 -1.0505919828009984; -1.0505919828009986 -0.8105094217626421]
ELECTRONIC ENERGY = -4.227529304014095
F' = [-1.506209922337154 -0.3630393122251443; -0.36303931222514385 -0.15290658995087186]
C' = [0.9698409901720303 0.24373849466618105; 0.24373849466618078 -0.9698409901720301]
ε  = [-1.5974482349709667, -0.06166827731705905]
C  = [0.801917475038031 -0.7822652596745968; 0.33680053251778796 1.0684455880206865]
P  = [1.2861432735427423 0.5401724652562576; 0.5401724652562576 0.2268691974085311]
δ (CONVERGENCE OF DENSITY MATRIX) = 1.3079512369232684e-5

ELECTRONIC ENERGY = -4.227525525484362
TOTAL ENERGY = -2.86065838497042

G  = [1.1941085339307838 0.29661318951207055; 0.2966131895120704 0.9213189761538769]
F  = [-1.458636190749471 -1.0505919828009984; -1.0505919828009986 -0.8105094217626421]
ε  = [-1.5974482349709667, -0.06166827731705905]
C  = [0.801917475038031 -0.7822652596745968; 0.33680053251778796 1.0684455880206865]
P  = [1.2861432735427423 0.5401724652562576; 0.5401724652562576 0.2268691974085311]

PS = [1.5296370380671058 1.1199277981091527; 0.6424383867623044 0.4703629619328946]

$R=1.4 a_0$ での $\mathrm{H}_2$ 分子は以下のように計算します.

入力

HF(R=[0 0 0; 0 0 1.4], Z=[1,1], ζ=[1.24,1.24])

出力

STO-3G FOR ATOMIC NUMBERS 1 AND 1

R	ζ1	ζ2	S12	T11
1.4	1.24	1.24	0.659319009323656	0.7600329435650853
T12	T22	V11A	V12A	V22A
0.2364553122909982	0.7600329435650853	-1.2266154680582526	-0.5974182717025165	-0.6538282202290772
V11B	V12B	V22B	V1111	V2111
-0.6538282202290772	-0.5974182717025164	-1.2266154680582526	0.7746083600328786	0.4441090500220031
V2121	V2211	V2221	V2222
0.29702947357791903	0.5696777252566694	0.44410905002200307	0.7746083600328786

S = [1.0 0.659319009323656; 0.659319009323656 1.0]
X = [0.5489339075446148 1.2114655016152216; 0.5489339075446148 -1.2114655016152216]
Hcore = [-1.1204107447222444 -0.9583812311140347; -0.9583812311140347 -1.1204107447222444]

(1,1,1,1)  0.7746083600328786
(1,1,1,2)  0.44410905002200307
(1,1,2,1)  0.4441090500220031
(1,1,2,2)  0.5696777252566694
(1,2,1,1)  0.4441090500220032
(1,2,1,2)  0.29702947357791903
(1,2,2,1)  0.29702947357791903
(1,2,2,2)  0.4441090500220031
(2,1,1,1)  0.4441090500220031
(2,1,1,2)  0.29702947357791903
(2,1,2,1)  0.29702947357791903
(2,1,2,2)  0.4441090500220032
(2,2,1,1)  0.5696777252566694
(2,2,1,2)  0.4441090500220031
(2,2,2,1)  0.44410905002200307
(2,2,2,2)  0.7746083600328786

P  = [0.0 0.0; 0.0 0.0]

START OF ITERATION NUMBER = 1
G  = [0.0 0.0; 0.0 0.0]
F  = [-1.1204107447222444 -0.9583812311140347; -0.9583812311140347 -1.1204107447222444]
ELECTRONIC ENERGY = -0.0
F' = [-1.2527982649241152 0.0; 0.0 -0.475604797574814]
C' = [-1.0 0.0; -0.0 1.0]
ε  = [-1.2527982649241152, -0.475604797574814]
C  = [-0.5489339075446148 1.2114655016152216; -0.5489339075446148 -1.2114655016152216]
P  = [0.6026568697043994 0.6026568697043994; 0.6026568697043994 0.6026568697043994]
δ (CONVERGENCE OF DENSITY MATRIX) = 0.6026568697043994

START OF ITERATION NUMBER = 2
G  = [0.7548736629110914 0.3644955517065228; 0.3644955517065228 0.7548736629110914]
F  = [-0.3655370818111531 -0.5938856794075119; -0.5938856794075119 -0.3655370818111531]
ELECTRONIC ENERGY = -1.8310009829233074
F' = [-0.5782027179991921 0.0; 0.0 0.6702710272828108]
C' = [-1.0 0.0; -0.0 1.0]
ε  = [-0.5782027179991921, 0.6702710272828108]
C  = [-0.5489339075446148 1.2114655016152216; -0.5489339075446148 -1.2114655016152216]
P  = [0.6026568697043994 0.6026568697043994; 0.6026568697043994 0.6026568697043994]
δ (CONVERGENCE OF DENSITY MATRIX) = 0.0

ELECTRONIC ENERGY = -1.8310009829233074
TOTAL ENERGY = -1.1167152686375932

G  = [0.7548736629110914 0.3644955517065228; 0.3644955517065228 0.7548736629110914]
F  = [-0.3655370818111531 -0.5938856794075119; -0.5938856794075119 -0.3655370818111531]
ε  = [-0.5782027179991921, 0.6702710272828108]
C  = [-0.5489339075446148 1.2114655016152216; -0.5489339075446148 -1.2114655016152216]
P  = [0.6026568697043994 0.6026568697043994; 0.6026568697043994 0.6026568697043994]

PS = [0.9999999999999996 0.9999999999999996; 0.9999999999999996 0.9999999999999996]

結果

上記の出力はComputational Chemistry List (CCL)のout.txtと見比べてください. 概ね, 同じ結果が得られてると思います. ズレの結果は直交化・対角化のルーチンと, 誤差関数erf()の実装の違い等だと思われます.

HeH⁺

作者	言語	ELECTRONIC ENERGY	TOTAL ENERGY
Szabo	FORTRAN	$-4.227~529~132~03$	$-2.860~661~991~52$
Ohno	Julia	$-4.227~525~525~48$	$-2.860~658~384~97$

H₂

作者	言語	ELECTRONIC ENERGY	TOTAL ENERGY
Szabo	FORTRAN	$-1.831~000~883~01$	$-1.116~715~168~72$
Ohno	Julia	$-1.831~000~982~92$	$-1.116~715~268~64$

注意点

このプログラムの注意点について述べます.

IOPによる分岐は全て削除した.
基底関数もSTO-3Gのみに限定した. STO-nGで使うにはループの上限を変える必要がある.
分子積分で使われるDERFOTHER(ARG)はSpecialFunctions.jlのerf(x)に置き換えた.
DIAG(F,C,E)はLinearAlgebra.jlを使てε,C = eigen(F)に置き換え, 正準直交化やRoothann方程式を解く際に用いた.
LinearAlgebra.jlのロードに数秒かかるので, たった1回だけ実行するならそのまま移植した方が早い.
二原子（原子核が2つの）分子を前提としているので, 同じ2電子系でもヘリウム原子ではコードの変更が必要.
付録Aの解説, 付録Bの実装はs型軌道のみの公式を与えている.
$\mu,\nu,\lambda,\sigma$ は行列要素の添え字, $i,j,k,l$ は基底関数の添え字とした
$N$ はSTO- $N$ Gの $N$ なので注意. 電子数は2で固定である.
ザボのプログラムではV1111やTTが二電子積分だが, TwoEに変えた.
VAやVBがかなりズレる. 誤差関数erf()の実装の違いか, コードに誤りがあるか, FORTRAN側が単精度になってしまっている部分があるか.
重なり行列の対角要素が1にならないので, 付録Bと同じように強制的に1にしている.

S[1,1] = 1
S[2,2] = 1

$\pmb{X}$ や $\pmb{C}$ の要素の符号が異なることがある. 計算結果は同じなので気にしなくてよい. アルゴリズムによって固有ベクトルの符号が変わることがあるものの, 固有ベクトルであることには変わりない. ここでは揃えるために下記のコードを加えた. 手続き(8)も同様.

s = sort(s, rev=true)
U = U[:, sortperm(s)]
U[2,:] .*= -1

MULT(A,B,C,IM,M)はA*B, MATOUT(A,IM,IN,M,N,LABEL)はprintln("A = $A")のように置き換えた.

解説

手続き(1)

ここでは基底関数についてのみ解説します.

STO- $n$ G基底関数系

分子積分の多重ループは短縮基底という概念を知らないと意味が分からないと思います. 式(3.283)を確認しておいてください. STO- $n$ G基底関数系も短縮基底の一種です. STO- $n$ G基底関数系の各係数・指数は水素の1s軌道と重なりを最大化するように決定されています. エネルギーを最小化するわけではありませんので, 変分的に最適化すればもっとエネルギーは下げられます. そうしないのは, STO- $n$ G基底関数系は初期構造の推定に使うためのものであって, 変分的に良い解が得られれば良いというものではないからです. 実際, エネルギーが改善しても（下がっても）, 核-電子間距離の期待値などは悪化することがあり, 最適化のジレンマと呼ばれています. 興味がある方はこちらの教育的な論文などを読んでください. STO- $n$ G基底関数系は下記のような表式です.

\phi_\mathrm{1s}^{\mathrm{STO}-\mathrm{NG}}(\pmb{r})=\sum_{i=1}^N C_i \phi_{1 \mathrm{~s}}^{\mathrm{GF}}\left(\alpha_i, \pmb{r}\right)

\phi_\mathrm{1s}^{\mathrm{GF}}(\alpha, \pmb{r})=\left(\frac{2 \alpha}{\pi}\right)^{3 / 4} e^{-\alpha r^2}

具体的な係数・指数のデータは

等で配布されています（ζの因子が掛けられた後なので注意）. FortranではSTO-1,2,3Gのデータを下記のように格納しています.

DATA COEF,EXPON/1.0D0,2*0.0D0,
0.678914D0,0.430129D0,0.0D0,
0.444635D0,0.535328D0,0.154329D0

0.270950D0,2*0.0D0,
0.151623D0,0.851819D0,0.0D0,
0.109818D0,0.405771D0,2.22766D0/

分子積分のルーチンの中で定義されていて違和感を覚えますが, 基底関数と分子積分は切り離せないという思想の表れでしょうか. これはJuliaで書くと

COEF = [
1.0  0.678914  0.444635
0.0  0.430129  0.535328
0.0  0.0       0.154329
]

EXPON = [
0.270950  0.151623  0.109818
0.0       0.851819  0.405771
0.0       0.0       2.22766
]

といった具合です. 今回はHF()のキーワード引数としてα = [0.109818 0.405771 2.22766], c = [0.444635 0.535328 0.154329]のように渡す仕様にしました. 描写してみると, カスプ以外はそこそこ良い記述となっていることがわかります.

# using Pkg
# Pkg.add("Plots")
using Plots

COEF = [
1.0  0.678914  0.444635
0.0  0.430129  0.535328
0.0  0.0       0.154329
]

EXPON = [
0.270950  0.151623  0.109818
0.0       0.851819  0.405771
0.0       0.0       2.22766
]

exact(r) = 1/sqrt(pi)*exp(-r)
STOnG(r,α,C) = sum(C .* exp.(-α.*r^2) .* (2*α./pi).^(3/4))

plot(xlims=(0,3), xlabel="\$r\$", ylabel="\$\\psi_{1s}(r)\$", size=(500,400))
plot!(0:0.01:3, r->STOnG(r,EXPON[:,1],COEF[:,1]), label="STO-1G")
plot!(0:0.01:3, r->STOnG(r,EXPON[:,2],COEF[:,2]), label="STO-2G")
plot!(0:0.01:3, r->STOnG(r,EXPON[:,3],COEF[:,3]), label="STO-3G")
plot!(0:0.01:3, r->exact(r), lc="#000000", ls=:dash, label="exact")

手続き(2)

分子積分の公式は付録Aに載っています. 分子積分は最初に計算するだけで, SCFの手続きの中では変化しません. Gauss関数を用いることで, 解析的に求めることができます.

一電子積分

例えば重なり積分は

\begin{align} (A \mid B) = [\pi /(\alpha+\beta)]^{3 / 2} \exp \left[-\alpha \beta /(\alpha+\beta)\left|\boldsymbol{R}_A-\boldsymbol{R}_B\right|^2\right] \tag{A.9} \end{align}

であり, この実装は

C*********************************************************************
      FUNCTION S(A,B,RAB2)
C
C CALCULATES OVERLAPS FOR UN-NORMALIZED PRIMITIVES
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DATA PI/3.1415926535898D0/
      S=(PI/(A+B))**1.5D0*DEXP(-A*B*RAB2/(A+B))
      RETURN
      END

です. これは, Juliaでは

integral_S(α, β, RA, RB) = (π/(α+β))^(3//2) * exp(-α*β*sum((RA-RB).^2)/(α+β))

の一行だけで書けます. ただし, 核間距離を入力するところを, 核の座標を2つ入力するように置き換えました. 関数側には距離だけ渡すという設計は, 核間距離の計算回数を減らせるという観点では合理的でしたが, 理解の妨げになるので今回は座標を渡すだけにしました. また, 並列化を前提にすると, 事情が変わると思われます. 以上のように宣言した関数の値を下記のように配列に格納します.

S = zeros(2,2)
for μ=1:2
  forν=1:2
    for i=1:3
      for j=1:3
        S[μ,ν] += D[μ,i] * D[ν,j] * integral_S(A[μ,i], A[ν,j], R[μ,:], R[ν,:])
      end
    end
  end
end

なお, 重なり行列は対称行列かつ対角要素は1だとわかっているので, 右上か左下の半分だけ計算すれば良いのですが, コードをシンプルにするために全て計算しています. また, ループは

for μ=1:2, ν=1:2
  for i=1:3, j=1:3
    S[μ,ν] += D[μ,i] * D[ν,j] * integral_S(A[μ,i], A[ν,j], R[μ,:], R[ν,:])
  end
end

のように略記しました. なお, 対角要素がぴったり1にならないので, ザボと近い結果を得るために

S[1,1] = 1
S[2,2] = 1

を加えておきます. Juliaのコードが誤りということではなく, ザボと倍精度レベルで一致します. 基底関数の指数の桁数が少ないためか, あるいは規格化の方法が間違っているためだと考えています.

1.0000014259978642 0.4507704116477877 # FORTRAN
1.0000014259978642 0.4507704116477877 # Julia

実際, 数値的にもこの値になることが確かめられます.

STOnG(r, C, α) = sum(C .* exp.(-α.*r^2))

for n in 1:3
  println("\nSTO-$(n)G    \t<φ|φ>")
  dr = 0.0001

  C = COEF[:,n]
  α = EXPON[:,n]
  C .*= ((2*α/π).^(3//4))
  println("analytical\t", sum(i->sum(j-> C[i] * C[j] * (π/(α[i]+α[j]))^(3//2) ,1:n),1:n))
  println("numerical \t", sum(r->(STOnG(r,C,α))^2 * 4*π*r^2*dr, 0.0:dr:100.0))

#   C = COEF[:,n]
#   α = EXPON[:,n]
#   C /= sqrt(sum(i->sum(j-> C[i] * C[j] * (π/(α[i]+α[j]))^(3//2) ,1:n),1:n))
#   println("analytical\t", sum(i->sum(j-> C[i] * C[j] * (π/(α[i]+α[j]))^(3//2) ,1:n),1:n))
#   println("numerical \t", sum(r->(STOnG(r,C,α))^2 * 4*π*r^2*dr, 0.0:dr:100.0))
end

STO-1G    	<φ|φ>
analytical	0.9999999999999997
numerical 	1.0000000000000002

STO-2G    	<φ|φ>
analytical	1.000001296011168
numerical 	1.000001296011168

STO-3G    	<φ|φ>
analytical	1.0000014259978642
numerical 	1.0000014259978638

二電子積分

我々の分野では常識ですが, 2中心のガウス型軌道の積が1中心のガウス型軌道に直せる（3.211）という最重要事項は自分で確認しておいてください. 二電子積分は8重のループで計算します.

integral_TwoE(α, β, γ, δ, RA, RB, RC, RD) = 2*(π^(5//2))/((α+β)*(γ+δ)*sqrt(α+β+γ+δ)) * exp(-α*β/(α+β)*sum((RA-RB).^2)-γ*δ/(γ+δ)*sum((RC-RD).^2)) * F0((α+β)*(γ+δ)/(α+β+γ+δ)*sum(((α*RA+β*RB)/(α+β)-(γ*RC+δ*RD)/(γ+δ)).^2))

TwoE = zeros(2,2,2,2)
for μ=1:2
  for ν=1:2
    for σ=1:2
      for λ=1:2
        for i=1:3
          for j=1:3
            for k=1:3
              for l=1:3
                TwoE[μ,ν,σ,λ] += D[μ,i] * D[ν,j] * D[σ,k] * D[λ,l] * integral_TwoE(A[μ,i], A[ν,j], A[σ,k], A[λ,l], R[μ,:], R[ν,:], R[σ,:], R[λ,:])
              end
            end
          end
        end
      end
    end
  end
end

さすがに8重は読みにくいので, 下記のように略記しました.

for μ=1:2, ν=1:2, σ=1:2, λ=1:2
  for i=1:3, j=1:3, k=1:3, l=1:3
    TwoE[μ,ν,σ,λ] += D[μ,i] * D[ν,j] * D[σ,k] * D[λ,l] * integral_TwoE(A[μ,i], A[ν,j], A[σ,k], A[λ,l], R[μ,:], R[ν,:], R[σ,:], R[λ,:])
  end
end

$\mu,\nu,\lambda,\sigma$ は行列要素の添え字, $i,j,k,l$ は短縮基底の添え字です. ここでD[μ,i] * D[ν,j] * D[σ,k] * D[λ,l]かc[μ,i] * c[ν,j] * c[σ,k] * c[λ,l]ではない理由は, 指数に $\zeta$ をかけて原子に合わせて変えたり, 短縮基底を用いていること等を考慮して規格化するです.

手順(3)

正準直交化についてです. 今回のような小さい基底関数系で正準直交化を行う必要はなく, eigen(F,S)のように渡してしまえばコレスキー分解で直交化して勝手に解いてくれます. しかし, より大きな基底関数系で計算しようとすると, 悪条件となり, 単純なコレスキー分解のみでは正常な固有値が求められないことが多くなります. このため, 正準直交化を行い, 重なり行列の固有値が小さいもの（重なりが大きいもの・線形従属性が大きいもの）を取り除くという操作が必要になります. この時に重なり行列の固有値を知る必要がありますので, 重なり行列 $\pmb{S}$ を対角化します. つまり固有値問題 $\pmb{S}\pmb{U} = s\pmb{U}$ を解けばよいわけです. これは式(3.166)に相当します. Juliaではusing LiearAlgebra.jlを宣言してから

s, U = eigen(S)

とするだけです. この時のsは固有値を格納したベクトルなので, diagm(s)を用いることで, その固有値を対角要素に並べた行列 $\pmb{s}$ が利用できます. これを用いて

X = U * (diagm(s) ^ (-1//2))

とすれば, Xが得られます.

手順(4)

密度行列を全ての要素が0の行列として宣言するだけです.

# 手順(4)  密度行列に対する初期値を設定する.
P = zeros(2,2)

手順(5)

行列 $\pmb{G}$ の計算です.

G_{\mu\nu} = \sum_{\lambda,\sigma} P_{\lambda\sigma} \left[ (\mu\nu\mid\sigma\lambda) - \frac{1}{2}(\mu\lambda\mid\sigma\nu) \right] \tag{3.154}

これをFORTRANのコードを参考にそのまま実装すると

function FORMG(G,P,TT)
  for I in 1:2
    for J in 1:2
      G[I,J] = 0.0
      for K = 1:2
        for L = 1:2
          G[I,J] += P[K,L]*(TT[I,J,K,L]-0.5*TT[I,L,K,J])
        end
      end
    end
  end
end

のようになりますが, Juliaの2つの記法

julai> [10*i+j for i=1:3, j=1:4]
3×4 Matrix{Int64}:
 11  12  13  14
 21  22  23  24
 31  32  33  34

julia> sum(i->i^2, 1:5)
55

を組み合わせると, 1行で書けます.

# 手順(5)  式(3.154)  Gの計算
G = [sum(λ-> sum(σ-> P[λ,σ] * (TwoE[μ,ν,σ,λ] - 1//2*TwoE[μ,λ,σ,ν]), 1:2), 1:2) for μ=1:2, ν=1:2]

リーダブルではないコードの代表例のようになっている気もしますが, 数式に近いという観点ではこれに勝るものはありません.

手順(6)

Fock行列 $\pmb{F}$ の計算は単に2つの行列を足すだけです.

# 手順(6)  式(3.154)  Fの計算
F = Hcore + G

手順(7-9)

この3つの手順でやっていることは, 下記の一般固有値問題を解くということだけです.

\begin{align} \pmb{F} \pmb{C} = \pmb{S} \pmb{C} \pmb{\varepsilon} \tag{3.159} \end{align}

コードもシンプルです.

# 手順(7)  式(3.177)  F'の計算（正準直交化）
F° = X° * F * X
# 手順(8)  式(3.178)  F'C'=C'εを解く
ε, C° = eigen(F°)
# 手順(9)  式(3.174)  Cの計算
C = X * C°

手順(10-11)

手順(9)で得られたCを使って密度行列Pを計算します. その前に密度行列をP_oldとして保存しておき, 収束判定で使います.

# 手順(10)  式(3.145)  Pの計算
P_old = P
P = [sum(k->2*C[i,k]*C[j,k], 1:Int(2/2)) for i=1:2, j=1:2]
# 手順(11)  収束判定
δ = sqrt(sum((P-P_old).^2) / 4)

これらの計算はWhile文

while δ>1e-4
～ここを繰り返す～
end

に挟まれていますので, δが $10^{-4}$ 以下なら終了します.

手順(12)

E_0 = \frac{1}{2} \sum_{\mu,v} P_{\mu\nu}\left(H_{\mu \nu}^\mathrm{core} + F_{\mu\nu} \right) \tag{3.184}

エネルギーの収束の様子を見るためにループ内でもエネルギーを計算していますが, この手順(12)をまとめてみると下記のようになります.

# 手順(12)  式(3.184)  エネルギー・マリケン電荷などの計算
EN = sum(P/2 .* (Hcore + F))
ENT = EN + Z[1] * Z[2] / norm(R[1,:]-R[2,:])
println("\nELECTRONIC ENERGY = $EN")
println("TOTAL ENERGY = $ENT\n")
println("\nPS = $(P * S)")

もっと学びたい人へ

Hartree-Fock極限

Hartree-Fock極限のエネルギー $E_\mathrm{HF}$ と厳密解のエネルギー $E_\mathrm{exact}$ の差を相関エネルギー(correlation energy)といい, $E_\mathrm{corr}$ で表します.

E_\mathrm{corr} := E_\mathrm{exact} - E_\mathrm{HF}

Hartree-Fock法では, いくら頑張っても厳密解には収束しないことが知られています. この $E_\mathrm{corr}$ を取り込める手法がpost-HF法（CI,MP,CC）や第一原理モンテカルロ法（VMC, FN-DMC）などです.

ポテンシャルエネルギー曲面の描写

Hartree-Fock法が求めているのは与えられた核配置1点 $\left\{\pmb{R}_A\right\}$ に対する $E_{\mathrm{elec}}\left\{\pmb{R}_A\right\}$ ですので, ポテンシャルエネルギー曲線 $E_{\mathrm{tot}}\left(\left\{\pmb{R}_A\right\}\right)$ を求めるには, 核配置を変えて計算しなくてはいけません.

X = 0.2:0.1:10.0
Y = []
for x in X
    E = HF(R=[0 0 0; 0 0 x], Z=[1,1], ζ=[1.24,1.24])
    push!(Y,E)
end

plot(X, Y, label="", xlabel="\$R\$", ylabel="\$E_\\mathrm{tot}\\{R\\}\$", lw=2)

RHFは $\mathrm{H}_2 \rightarrow 2\mathrm{H}$ の解離閾値に漸近せず, より高い値に向かっていきます. このように, RHFは解離生成物が開殻系になるような解離反応を記述できないことが知られています. つまり $\mathrm{Li}_2$ でも同様の問題が起きうるということです. 余力のある人はGaussian16やGAMESS等のパッケージを用いて検証してみましょう. これら問題はCI法では回避できますが, MP法やDFT法の制限付き計算（RMP2やRB3LYPなど）では回避できません. これを機にいろいろな計算手法を試してみましょう.

参考文献

A.ザボ, N.S.オストランド著, 大野公男, 阪井健男, 望月祐志訳『新しい量子化学上　電子構造の理論入門』(東京大学出版会, 1987)
Computational Chemistry List (CCL)
J. F. Pérez-Torres, J. Chem. Educ., 96, 4, 704–707 (2019) https://doi.org/10.1021/acs.jchemed.8b00959
STO- $N$ G基底関数・最適化のジレンマについて
Attila Szabo, Neil S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (McGraw-Hill, 1989)
Jensen, F. (1999). The Journal of chemical physics, 110(14), 6601-6605. https://doi.org/10.1063/1.478567
HF極限への収束の様子
Sundholm, D. (1988). Chemical physics letters, 149(3), 251-256. https://doi.org/10.1016/0009-2614(88)85022-X
$E_\mathrm{HF}$ の上限
Y. Kurokawa, H. Nakashima, and H. Nakatsuji, Phys. Rev. A 72, 062502. https://doi.org/10.1103/PhysRevA.72.062502
$E_\mathrm{Exact}$ の上限
IUPAC GreenBook 3rd edition (2007) https://doi.org/10.1039/9781847557889
1𝐸_h=27.21138 eV=627.5095 kcal/mol