🎁

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

2022/12/24に公開・約38,300字

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法をベンチマークとして評価される.（相関エネルギーの取り込み率は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だとうまく計算できないはずですが, なぜか計算できてしまっています.

付録：第一原理計算とは

第一原理計算の公式な定義はIUPACのGoldbookに載っています. ab initio calculationsという項目があり, このページにはなかなか的外れな定義が書いてありましたが, 現在はちゃんとした定義に置き換えられています. 重要な部分だけ抜き出すと,

基本物理定数以外の実験値を用いない量子力学的な計算手法のこと（Methods of quantum mechanical calculations independent of any experiment other than the determination of fundamental constants.）

が第一原理計算の定義です. よく勘違いされますが, 近似を使おうとも, 基礎物理定数以外の実験値を用いなければ第一原理です. Gaussian等で実装されているHartree-Fock法では, Born-Oppenheimer近似はもちろんのこと, 平均場近似も使われていますが, 第一原理です.

第一原理計算にはIUPACによる定義がありますが, 量子化学計算という言葉の公式な定義はないと思います. ここでは, 量子化学計算の定義は分子系の電子状態を求める手法としておきましょう. 下表は全て量子化学計算の手法に分類されるものとして異論はないと思います. 第一原理のものとそうでないものに分類できます.

分類	具体例	パッケージなど
半経験的	拡張Hückel法, AM1, PM3, PM6, PM7, INDO, CNDO, MINDO, …	MOPAC, Winmostar(MOPACを内蔵)
第一原理非経験的	HF, CI, MP, CC, ECG, VMC, DMC	Gaussian, GAMESS, Firefly, Molpro, Fermi.jl, CASINOなど

よくある質問

非経験的と第一原理は違うの？
非経験的はnon-empiricalの和訳, 第一原理はab initioやfirst principleの和訳に対応します. IUPACによれば同じ意味とされます.
そもそも「経験」ってなに？
Hückel法など初期の計算手法では積分の代わりに経験的パラメータ（実験値そのものや, 実験に合うように調整されたパラメータ, 職人の勘による近似）を使っていました. これが「経験」の所以です. 一部をちゃんと計算するようになったものが半経験的手法, 全ての積分をまじめに計算するようになったのが非経験的手法・第一原理計算です.
半経験的手法って今も使われてるの？
現在も第一原理計算パッケージでは初期値推定のために使われています. 実用上も, 第一原理計算より圧倒的に速いので, 大きい分子では計算してみる価値があります.
量子化学計算は第一原理計算なの？
量子化学計算には第一原理の手法と第一原理でない手法があります. 例えばMOPAC等で実装されているAM1やPM3は半経験的分子軌道法であり, 第一原理ではありませんが, Gaussian等で実装されているHF, CI, MP, CCは第一原理計算です.
量子モンテカルロ法は第一原理計算ですか？
第一原理の量子モンテカルロ法と, そうでない量子モンテカルロ法があります. 例えば, CASINOの変分モンテカルロ法（VMC）や拡散量子モンテカルロ法（DMC）は量子化学計算かつ第一原理計算に分類されるものです. Hubbard模型などスピン系でも変分モンテカルロ法は使いますが, これは量子化学計算でも第一原理計算でもないです.
第一原理計算は変分法ですか？
そうとは限りません. HFやCIは変分法ですが, MPやCCは変分原理を破ります.
DFTは第一原理計算ですか？
第一原理計算だと思っている方が多いようですので, 第一原理計算ということにしておく方が穏便です.
量子化学計算ではない第一原理計算の具体例は？
固体物理などで, VASP等のパッケージを用いて平面波基底で結晶など周期系の計算をするものは量子化学計算の範疇ではないと思います. その中に第一原理のものがあれば, それが答えです.
分子動力学法は第一原理ですか？
いわゆる古典MDは第一原理ではありません. 力場が簡単な関数形に限定されており, その形状を決めるパラメータに実験値などが使われています.
第一原理分子動力学法は第一原理ですか？
量子分子動力学法＝第一原理分子動力学(ab initio MD)は経験的な力場を使っておらず, いわゆる古典MDと比べれば電子状態を第一原理的に扱っているという意味で, 第一原理の名を冠しているものだと思います. しかし, 核の運動を量子力学的を扱っているわけではありませんので, IUPACの定義に照らすと「量子力学的な計算ではない」という意味で, 第一原理計算には含まれないと思います.
経路積分分子動力学法は第一原理ですか？
いわゆるPIMDは核の運動は量子力学的に扱っていますが, 電子状態の第一原理計算をせずに古典力場を使っているものは第一原理ではありません. 第一原理PIMDは電子も核も第一原理的に扱っているという意味で第一原理計算に分類されてもおかしくないと思います.

参考文献

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

Discussion

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

著者: 大野周平
公開: 2022/12/24
本文更新: 2023/01/30
文章量: 約38,300字

大野周平

計算物理・計算化学の研究をしている大学院生です. 主にJulia・Fortranについての記事, たまにWeb系の記事を書いています. ご提案・ご質問等はコメント欄までお願いします. Web制作の事業は３月末に廃業予定です. 現在は案件をお受けしておりません.

バッジを贈るとは？