Implementing the Hartree-Fock Method in Julia

JuliaLang Advent Calendar 2022 Day 24🎄
Today I was planning to release a package I've been developing personally, but since we covered the Hartree-Fock method in a recent seminar, I decided to implement it in Julia instead.

Referencing Appendix B of "Modern Quantum Chemistry" by Szabo and Ostlund (1987), I will implement a Hartree-Fock method program (about 140 lines) in the Julia language. The original Appendix B code (about 570 lines), which I used as a reference, is distributed by the Computational Chemistry List (CCL). Since it is written in the ancient language FORTRAN, it is difficult to decipher, but I have created a program with almost the same variable names and outputs, so the correspondence should be apparent. The original code solves HeH⁺ at the HF/STO-3G level, but by changing the charges and orbital exponents, it can handle two-nucleus two-electron systems such as H₂, HeH⁺, and LiH²⁺. This time, I calculated H₂ and HeH⁺ using Julia.

Introduction

As mentioned in this article, students assigned to research labs related to computational chemistry end up studying "Modern Quantum Chemistry" by Szabo and Ostlund (1987) (commonly known as "Szabo"). I also studied Szabo thoroughly during my third year of undergraduate studies. Back then, I thought, "Hartree-Fock is outdated," but thanks to the aftereffects experience of that harsh thorough seminar, I believe I have come to be made to understand the importance of the Hartree-Fock method. The Hartree-Fock (HF) method is a type of ab initio molecular orbital method and holds a special place among the many electronic structure theories. The reasons for this include:

• The accuracy of electronic structure calculations is evaluated using the HF method as a benchmark. (The capture rate of correlation energy (the BSCE one) is set to 0% for the HF method and 100% for full CI or the exact solution.)
• It is important when studying Density Functional Theory (DFT).
• Post-HF methods (CI, MP, CC, etc.) are developed as extensions based on HF.
• Quantum Monte Carlo methods (ab initio QMC such as VMC, FN-DMC) also calculate with the HF method first and use those results for sampling.
• Even in ab initio molecular dynamics (ab initio MD), which does not use classical force fields, some form of electronic structure calculation is performed.

As long as you are in this field, you will be involved with the HF method indirectly, even if you don't use it directly. That said, general-purpose HF codes are included in packages such as:

• Gaussian
• GAMESS
• Firefly
• Molpro

Since these are available commercially or for free, I do not recommend implementing them yourself from the perspective of effort and reliability. My objectives for implementing it anew in Julia are:

1. For development of new methods
   For problems that cannot be calculated with existing packages, you may need to extend packages or develop them from scratch. Rather, I think that is what is expected of researchers in the computational field. Therefore, having a small, easily runnable program makes it easier to develop prototypes.
2. For personal learning
   For the same reasons as above, you would likely want to understand the general flow of calculations with a simple program rather than trying to extend a massive package suddenly.
3. For education
   The Appendix B code is written in the ancient language FORTRAN. Even in the 21st century, there are no students who can read or want to study FORTRAN.

These are some of the reasons. The reasons I chose Julia among the many programming languages are as follows:

• Specifically in the field of numerical calculation, the code is simpler than Python, etc.
• Since Greek characters and other symbols can be used, the correspondence between mathematical formulas and code is easy to understand.
• Validated numerics (rigorous numerics) can be used.
• Quad-precision, oct-precision, and arbitrary-precision arithmetic can be used.

While there are packages like Fermi.jl in Julia, they are not pure native implementations, so validated or arbitrary-precision arithmetic cannot be used. Since the diagonalization of the overlap matrix is ill-conditioned, a native Julia implementation of a general-purpose code would be highly promising for various applications.

Schrödinger Equation and Born-Oppenheimer Approximation

Molecules are many-body systems of electrons and nuclei. As is clear from the potential terms of the Hamiltonian, it should be noted that both electrons and nuclei are treated as point masses and point charges. A common misunderstanding is that the so-called electron cloud refers to the probability density function of the position where an electron exists, and it is not a concept regarding the size of the electron.

When solving the Schrödinger equation for a molecular system, the Born-Oppenheimer approximation (BO approximation) is often used because the masses of electrons and nuclei are significantly different. The BO approximation is sometimes explained as "treating the electrons while keeping the nuclei fixed," but this is strictly a perspective from the electrons; a different explanation is given from the perspective of the nuclei. In the first place, fixing the nuclei or treating electrons by averaging them is merely an interpretation given to the results after the approximation and does not serve as an explanation of the approximation itself. The essence of the BO approximation is to separate the motion of electrons and nuclei.

Let's look at what separation means starting from the Hamiltonian. Three Hamiltonians appear in Section 2.1 of Szabo.

$$\begin{align} \hat{H} = - \sum_{i=1}^N \frac{1}{2} \nabla_i^2 - \sum_{i=1}^N \frac{1}{2} \frac{1}{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_{A B}} \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_{A B}} + E_{\mathrm{elec}}\left(\left\{\pmb{R}_A\right\}\right) \tag{2.15} \end{align}$$

Three corresponding Schrödinger equations also appear.

$$\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}$$

The total Hamiltonian of the molecule is (2.2), and the energy we want to find is $E$ in equation (2.1). However, in standard calculation packages, equation (2.1) is not solved. Instead, equation (2.11) is solved for the motion of electrons, and equation (2.16) is solved for the motion of nuclei. Thus, in the BO approximation, the motions of electrons and nuclei are treated in separate equations. Among these, $E_\mathrm{nucl}$ in equation (2.16) corresponds to $E$ in equation (2.1), and $E \simeq E_\mathrm{nucl}$. However, since the potential part $\sum_{A=1}^M \sum_{B>A}^M \frac{Z_A Z_B}{R_{A B}} + E_{\mathrm{elec}}\left(\left\{\pmb{R}_A\right\}\right)$ of $\hat{H}_\mathrm{nucl}$ in equation (2.15) contains $E_{\mathrm{elec}}\left(\left\{\pmb{R}_A\right\}\right)$, it can be seen that to find $E_\mathrm{nucl}$ by solving equation (2.16), equation (2.10) must be solved first (multiple times for different internuclear distances). This potential term felt by the nuclei is represented by the symbol

$$\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_{A B}} + E_{\mathrm{elec}}\left(\left\{\pmb{R}_A\right\}\right) \tag{2.14} \end{align}$$

and is considered a function of the nuclear configuration $\{\pmb{R}_A\}$, called a Potential Energy Curve (PEC) or Potential Energy Surface (PES). In practice, rather than finding the entire PEC, $E_\mathrm{nucl}$ is obtained by approximating the vicinity of the minimum point as a harmonic oscillator or by adding anharmonic corrections. The laboratory I belong to develops and applies ab initio multi-component molecular orbital methods, which do not use the BO approximation in the first place, and ab initio path-integral molecular dynamics, which treats nuclear quantum and temperature effects on the PES.

Note: $\{\pmb{R}_i\}$ means the tuple $(\pmb{R}_1,\pmb{R}_2,\pmb{R}_3,\cdots)$. The issue of confusion with set notation is apparently common even outside of Szabo.

Hartree-Fock Method

You might be confused by all the symbols that have appeared, but things will become clearer if you keep in mind that $E$ and $E_\mathrm{nucl}$ do not depend on the nuclear configuration $\pmb{R}_i$, while $E_\mathrm{elec}(\{\pmb{R}_i\})$ and $E_\mathrm{tot}(\{\pmb{R}_i\})$ do depend on the nuclear configuration $\{\pmb{R}\}$. The HF method approximately calculates $E_{\mathrm{elec}}(\{\pmb{R}_A\})$ for a single given nuclear configuration $\left\{\pmb{R}_A\right\}$.

In the HF method, a single Slater determinant (3.1) is used as the trial wavefunction. By minimizing the expectation value of the energy (3.2) under the constraint (3.3), we obtain the Hartree-Fock differential equation for the spin-orbitals:

$$\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}$$

Refer to Sections 3.1 and 3.2 for the derivation. In practice, the canonical Hartree-Fock equation (3.69), which is made unique via a unitary transformation, is solved. Since these are differential equations, they are reduced to matrix equations for computation. The result is the Roothaan equations:

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

In actual implementation, this is solved by orthogonalization (3.178). Refer to Section 3.4 for the derivation and orthogonalization. Since the Fock matrix $\pmb{F}$ depends on the solution $\pmb{C}$—meaning it is self-consistent—you cannot obtain the desired solution simply by solving the eigenvalue problem once. At first glance, it might even seem impossible to solve. However, by providing appropriate initial values and solving iteratively, the resulting $\pmb{C}$ and $\pmb{\varepsilon}$ converge to constant values. This algorithm is known as the SCF procedure. This is in contrast to methods like the linear variational method, where solving the generalized eigenvalue problem $\pmb{H}\pmb{C}=E\pmb{S}\pmb{C}$ once is sufficient.

Why is Implementing the HF Method Difficult?

In the context of the HF method, the SCF procedure is often emphasized, but implementing SCF itself is easy. In Julia, specifically, it can be written in about 10 lines. Rather, the difficult parts are the handling of basis functions and the preparation of molecular integrals. Specifically:

• Writing a parser for basis set data files
• Designing appropriate data structures for basis sets
• Deriving molecular integral formulas
• Implementing molecular integrals

Most of the effort is spent on these tasks. Indeed, these parts likely take up a significant number of pages in Szabo's book. These have already been implemented in GaussianBasis.jl (a wrapper for libcint) and are used in Fermi.jl. However, as of December 2022, build errors occur, so I implemented them this time with reference to Szabo.

Implementation

I have summarized the procedure from Szabo, p. 159, Section 3.4.6, in the table below. We will proceed to implement this.

Step	Description	Equation
1	Parameter input	$$\{\pmb{R}_A\}, \{Z_A\}, N, \{\phi_\mu\}$$
2	Calculation of molecular integrals	$$S_{\mu\nu}, H_{\mu\nu}^\mathrm{core}, (\mu\nu\mid\lambda\sigma)$$
3	Calculation of $\pmb{X}$ (Equation 3.167)	$$\pmb{X} = \pmb{U} \pmb{s}^{-1/2}$$
4	Set initial values for the density matrix $\pmb{P}$	$$\boldsymbol{P}=\left(\begin{array}{cc} 0 & \cdots \\ \vdots & \ddots \end{array}\right)$$
5	Calculation of $\pmb{G}$ (Equation 3.154)	$$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	Calculation of $\pmb{F}$ (Equation 3.154)	$$\pmb{F} = \pmb{H}^\mathrm{core} + \pmb{G}$$
7	Calculation of $\pmb{F}'$ (Canonical Orthogonalization, Equation 3.177)	$$\pmb{F}' = \pmb{X}'\pmb{F}\pmb{X}$$
8	Solve Equation (3.178) to obtain $\pmb{C}'$ and $\pmb{\varepsilon}$.	$$\pmb{F}'\pmb{C}' = \pmb{C}'\pmb{\varepsilon}$$
9	Calculation of $\pmb{C}$ (Equation 3.174)	$$\pmb{C}=\pmb{X}\pmb{C}'$$
10	Calculation of $\pmb{P}$ (Equation 3.145)	$$P_{\mu\nu} = 2 \sum_a^{N/2} C_{\mu a} C_{\nu a}^*$$
11	Return to (5). Terminate if $\delta < 10^{-4}$.	$$\delta = \frac{1}{4} \sum_{\mu,v} (P_{\mu\nu}^\mathrm{new} - P_{\mu\nu}^\mathrm{old})^2$$
12	Calculation of $E_0$, etc. (Equation 3.184)	$$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

# Step (1) Parameter input
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])")
    
  # Step (2) Calculation of molecular integrals
  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))

  # One-electron integrals
  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

  # Equation (3.233), etc.
  Hcore = T + VA + VB

  # Two-electron integrals
  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

  # Step (3) Calculation of X (Equation 3.167)
  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)
    
  # Step (4) Set initial values for the density matrix.
  P = zeros(2,2)
    
  # Output matrices
  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")

  # Other initial value settings
  δ = 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")
    # Step (5) Calculation of G (Equation 3.154)
    G = [sum(λ-> sum(σ-> P[λ,σ] * (TwoE[μ,ν,σ,λ] - (1//2)*TwoE[μ,λ,σ,ν]), 1:2), 1:2) for μ=1:2, ν=1:2]
    println("G  = $G")
    # Step (6) Calculation of F (Equation 3.154)
    F = Hcore + G
    println("F  = $F")
    # Step (12) Calculation of energy, Mulliken charges, etc. (Equation 3.184)
    EN = sum(P/2 .* (Hcore + F))
    println("ELECTRONIC ENERGY = $EN")
    # Step (7) Calculation of F' (Canonical Orthogonalization, Equation 3.177)
    F° = X° * F * X
    # Step (8) Solve F'C'=C'ε (Equation 3.178)
    ε, C° = eigen(F°)
    C°[:,1] *= -1
    # Step (9) Calculation of C (Equation 3.174)
    C = X * C°
    # Step (10) Calculation of P (Equation 3.145)
    P_old = P
    P = [sum(k->2*C[μ,k]*C[ν,k], 1:1) for μ=1:2, ν=1:2]
    # Step (11) Convergence test
    δ = 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

  # Step (12) Calculation of energy, Mulliken charges, etc. (Equation 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

I was able to shorten the code from about 570 lines to about 140 lines. It's just under 70 lines if you exclude print statements and comment lines.

Calculation

Input

HF()

Output

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]

The $\mathrm{H}_2$ molecule at $R=1.4 a_0$ is calculated as follows:

Input

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

Output

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]

Results

Compare the above output with out.txt from the Computational Chemistry List (CCL). Generally, you should see that almost the same results are obtained. The discrepancies are likely due to differences in the implementation of the orthogonalization and diagonalization routines, as well as the error function erf(), etc.

HeH⁺

Author	Language	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₂

Author	Language	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$$

Points to Note

I will describe some points to note regarding this program.

• All branches based on IOP (Internal Options) were removed.
• The basis set is limited to STO-3G only. To use it with STO-nG, the upper limit of the loops needs to be changed.
• DERFOTHER(ARG) used in molecular integrals was replaced with erf(x) from SpecialFunctions.jl.
• DIAG(F,C,E) was replaced by ε, C = eigen(F) using LinearAlgebra.jl, and used for canonical orthogonalization and solving the Roothaan equations.
• Since loading LinearAlgebra.jl takes a few seconds, it would be faster to port it directly if you only run it once.
• Since it assumes a diatomic (two-nucleus) molecule, changes to the code are required for a helium atom, even if it is also a two-electron system.
• The explanation in Appendix A and the implementation in Appendix B provide formulas only for s-type orbitals.
• I used $\mu, \nu, \lambda, \sigma$ as indices for matrix elements and $i, j, k, l$ as indices for the basis functions.
• Note that $N$ is the $N$ in STO-$N$G. The number of electrons is fixed at 2.
• In Szabo's program, V1111 and TT refer to two-electron integrals, but I changed them to TwoE.
• VA and VB deviate considerably. This might be due to differences in the implementation of the error function erf(), errors in the code, or parts of the FORTRAN side being in single precision.
• Since the diagonal elements of the overlap matrix do not become exactly 1, I forced them to 1 as in Appendix B.

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

• The signs of the elements in $\pmb{X}$ or $\pmb{C}$ may differ. This is not a concern as the calculated results are the same. While the sign of an eigenvector can change depending on the algorithm, it remains a valid eigenvector. I added the following code to align them. Step (8) is similar.

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

• MULT(A,B,C,IM,M) was replaced with A*B, and MATOUT(A,IM,IN,M,N,LABEL) was replaced with something like println("A = $A").

Explanation

Step (1)

Here, I will explain only the basis functions.

STO-$n$G Basis Set

The multiple loops in molecular integrals might not make sense without understanding the concept of contracted basis sets. Please check Equation (3.283). The STO-$n$G basis set is a type of contracted basis set. Each coefficient and exponent of the STO-$n$G basis set is determined to maximize the overlap with the hydrogen 1s orbital. They are not determined to minimize the energy, so the energy can be lowered further if you optimize variationally. The reason this is not done is that STO-$n$G basis sets are intended for estimating initial structures, and the goal isn't necessarily to obtain the best variational solution. In fact, even if the energy improves (decreases), the expectation value of the nucleus-electron distance may deteriorate, which is known as the optimization dilemma. If you're interested, please read this educational paper. The STO-$n$G basis set is expressed as follows:

$$\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}$$

Specific coefficient and exponent data are distributed at sites like:

• https://www.basissetexchange.org/basis/sto-2g/format/nwchem/
• https://www.basissetexchange.org/basis/sto-3g/format/nwchem/
• https://www.basissetexchange.org/basis/sto-4g/format/nwchem/
• https://www.basissetexchange.org/basis/sto-5g/format/nwchem/
• https://www.basissetexchange.org/basis/sto-6g/format/nwchem/

(Note that this is after the $\zeta$ factor has been multiplied). In Fortran, STO-1, 2, and 3G data are stored as follows:

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/

It feels a bit strange that they are defined within the molecular integral routine, but perhaps it reflects the idea that basis functions and molecular integrals are inseparable. In Julia, this can be written as follows:

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
]

For this implementation, I designed it to pass α = [0.109818 0.405771 2.22766], c = [0.444635 0.535328 0.154329] as keyword arguments to HF(). Plotting it shows that it provides a reasonably good description, except for the cusp.

# 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-6STOnG(r,EXPON[:,1],COEF[:,1]), label="STO-1G")
plot!(0:0.01:3, r-6STOnG(r,EXPON[:,2],COEF[:,2]), label="STO-2G")
plot!(0:0.01:3, r-6STOnG(r,EXPON[:,3],COEF[:,3]), label="STO-3G")
plot!(0:0.01:3, r-6exact(r), lc="#000000", ls=:dash, label="exact")

Step (2)

The formulas for molecular integrals are given in Appendix A. Molecular integrals are only calculated at the beginning and do not change during the SCF procedure. They can be determined analytically by using Gaussian functions.

One-electron Integrals

For example, the overlap integral is

$$\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}$$

and its implementation is

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

In Julia, this can be written in just one line:

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

However, I replaced the part where the internuclear distance is input with two sets of nuclear coordinates. The design of passing only the distance to the function was rational from the viewpoint of reducing the number of internuclear distance calculations, but I chose to pass coordinates directly this time as the previous approach could hinder understanding. Furthermore, the situation would likely change assuming parallelization. The values of the functions declared as above are stored in an array as follows:

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

Since it is known that the overlap matrix is symmetric and the diagonal elements are 1, it would be sufficient to calculate only half (upper right or lower left), but I calculated everything to keep the code simple. Also, the loops were abbreviated as follows:

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

Additionally, since the diagonal elements do not become exactly 1, I added the following to obtain results close to Szabo's:

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

This doesn't mean the Julia code is wrong; it matches Szabo's results at double-precision levels. I believe this is either because the number of digits in the exponents of the basis functions is small, or because the normalization method is different.

1.0000014259978642 0.4507704116477877 # FORTRAN
1.0000014259978642 0.4507704116477877 # Julia

In fact, it can be confirmed numerically that it becomes this value.

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

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

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

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

STO-1G    	6φ|φ7
analytical	0.9999999999999997
numerical 	1.0000000000000002

STO-2G    	6φ|φ7
analytical	1.000001296011168
numerical 	1.000001296011168

STO-3G    	6φ|φ7
analytical	1.0000014259978642
numerical 	1.0000014259978638

Two-electron Integrals

As is common knowledge in our field, please verify for yourself the most important point: that the product of two-center Gaussian-type orbitals can be reduced to a one-center Gaussian-type orbital (3.211). Two-electron integrals are calculated using an eight-fold loop.

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

Since an eight-fold loop is understandably hard to read, I have abbreviated it as follows:

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

Here, $\mu, \nu, \lambda, \sigma$ are the indices for the matrix elements, and $i, j, k, l$ are the indices for the contracted basis functions. The reason for using D[μ,i] * D[ν,j] * D[σ,k] * D[λ,l] instead of c[μ,i] * c[ν,j] * c[σ,k] * c[λ,l] is to ensure normalization by considering the use of contracted basis sets and multiplying the exponents by $\zeta$ according to the atoms.

Step (3)

This step is about canonical orthogonalization. In a small basis set system like this, there is no need to perform canonical orthogonalization; you could simply solve it automatically using Cholesky decomposition by passing it as eigen(F,S). However, when trying to calculate with larger basis sets, the system often becomes ill-conditioned, and desirable eigenvalues cannot be obtained using Cholesky decomposition alone. For this reason, it is necessary to perform canonical orthogonalization to remove elements with small eigenvalues of the overlap matrix (those with high overlap or linear dependency). To do this, we need to know the eigenvalues of the overlap matrix, so we diagonalize the overlap matrix $\pmb{S}$. In other words, we solve the eigenvalue problem $\pmb{S}\pmb{U} = s\pmb{U}$. This corresponds to Equation (3.166). In Julia, after declaring using LinearAlgebra, you simply do:

s, U = eigen(S)

Since s is a vector storing the eigenvalues, we can use diagm(s) to get a matrix $\pmb{s}$ with those eigenvalues as diagonal elements. Using this, we can calculate X as:

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

Step (4)

This involves simply declaring the density matrix as a matrix with all elements set to 0.

# Step (4) Set initial values for the density matrix.
P = zeros(2,2)

Step (5)

This is the calculation of matrix $\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}$$

If we implement this directly by referencing the FORTRAN code, it looks like this:

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

However, by combining two Julia notations:

julia> [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

It can be written in a single line:

# Step (5) Equation (3.154) Calculation of G
G = [sum(λ-> sum(σ-> P[λ,σ] * (TwoE[μ,ν,σ,λ] - 1//2*TwoE[μ,λ,σ,ν]), 1:2), 1:2) for μ=1:2, ν=1:2]

While this might feel like a prime example of non-readable code, in terms of proximity to the mathematical formula, nothing beats this.

Step (6)

The calculation of the Fock matrix $\pmb{F}$ is simply adding two matrices together.

# Step (6) Equation (3.154) Calculation of F
F = Hcore + G

Steps (7-9)

What is being done in these three steps is simply solving the following generalized eigenvalue problem:

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

The code is also simple:

# Step (7) Equation (3.177) Calculation of F' (Canonical Orthogonalization)
F° = X° * F * X
# Step (8) Equation (3.178) Solve F'C'=C'ε
ε, C° = eigen(F°)
# Step (9) Equation (3.174) Calculation of C
C = X * C°

Steps (10-11)

Using the C obtained in Step (9), calculate the density matrix P. Before that, save the density matrix as P_old to use for convergence testing.

# Step (10) Equation (3.145) Calculation of 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]
# Step (11) Convergence test
δ = sqrt(sum((P-P_old).^2) / 4)

Since these calculations are wrapped in a while loop:

while δ>1e-4
# ~ Repeat here ~
end

it will terminate once δ becomes $10^{-4}$ or less.

Step (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}$$

To observe the convergence of the energy, the energy is calculated within the loop, but summarizing Step (12) looks like this:

# Step (12) Equation (3.184) Calculation of energy, Mulliken charges, etc.
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)")

For those who want to learn more

Hartree-Fock Limit

The difference between the Hartree-Fock limit energy $E_\mathrm{HF}$ and the exact solution energy $E_\mathrm{exact}$ is called the correlation energy, denoted as $E_\mathrm{corr}$.

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

It is known that the Hartree-Fock method will not converge to the exact solution no matter how much effort is applied. Methods that can incorporate this $E_\mathrm{corr}$ include post-HF methods (CI, MP, CC) and ab initio Monte Carlo methods (VMC, FN-DMC).

Plotting Potential Energy Surfaces

Since the Hartree-Fock method calculates $E_{\mathrm{elec}}(\{\pmb{R}_A\})$ for a single given nuclear configuration $\left\{\pmb{R}_A\right\}$, to obtain the potential energy curve $E_{\mathrm{tot}}(\{\pmb{R}_A\)$, calculations must be performed for varying nuclear configurations.

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 does not converge to the dissociation threshold of $\mathrm{H}_2 \rightarrow 2\mathrm{H}$ but instead tends toward a higher value. Thus, it is known that RHF cannot describe dissociation reactions where the dissociation products are open-shell systems. In other words, similar issues can occur with $\mathrm{Li}_2$. For those interested, try verifying this using packages like Gaussian16 or GAMESS. These issues can be avoided with the CI method, but they cannot be avoided with restricted calculations in the MP method or Density Functional Theory (DFT) (such as RMP2 or RB3LYP). Let's take this opportunity to experiment with various calculation methods.

References

1. A. Szabo, N. S. Ostlund, translated by K. Ohno, T. Sakai, and Y. Mochizuki, "Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory" (University of Tokyo Press, 1987)
2. Computational Chemistry List (CCL)
3. J. F. Pérez-Torres, J. Chem. Educ., 96, 4, 704–707 (2019) https://doi.org/10.1021/acs.jchemed.8b00959
   On STO-$N$G basis functions and the optimization dilemma
4. Attila Szabo, Neil S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (McGraw-Hill, 1989)
5. Jensen, F. (1999). The Journal of chemical physics, 110(14), 6601-6605. https://doi.org/10.1063/1.478567
   Convergence behavior to the HF limit
6. Sundholm, D. (1988). Chemical physics letters, 149(3), 251-256. https://doi.org/10.1016/0009-2614(88)85022-X
   Upper limit of $E_\mathrm{HF}$
7. Y. Kurokawa, H. Nakashima, and H. Nakatsuji, Phys. Rev. A 72, 062502. https://doi.org/10.1103/PhysRevA.72.062502
   Upper limit of $E_\mathrm{Exact}$
8. IUPAC GreenBook 3rd edition (2007) https://doi.org/10.1039/9781847557889
   1 $E_h = 27.21138$ eV $= 627.5095$ kcal/mol