Hartree-Fock Method

Introduction

Understanding the electronic structure of atoms and molecules is fundamental to our grasp of solid-state physics and condensed matter. The Hartree-Fock (HF) method, a cornerstone of quantum chemistry, provides an approximate solution to the Schrödinger equation for many-electron systems.

The Schrödinger Equation for Many-Electron Systems

The time-independent Schrödinger equation for a system with \(N\) electrons and \(M\) nuclei is given by:

where \(\hat{H}\) is the Hamiltonian operator, \(\Psi\) is the wavefunction, and \(E\) is the energy eigenvalue. The Hamiltonian for a many-electron system is composed of the kinetic energy operators for the electrons and nuclei, the electron-nucleus potential energy, the electron-electron repulsion energy, and the nucleus-nucleus repulsion energy:

Solving this equation exactly for systems with more than one electron is not possible due to the electron-electron interaction term, and approximation is needed.

The Hartree-Fock Approximation

The Hartree-Fock method simplifies the many-electron problem by making some key approximations:

1. The Born-Oppenheimer approximation: It is based on the notion that nuclear motion is insignificant in comparison to electron mobility. As a result, we can regard the nuclei as fixed point charges and separate the electronic and nuclear wavefunctions.

2. The independent-particle approximation: According to HF, each electron travels on its own within an average potential created by the other electrons. By breaking down the many-body Schrödinger equation into a collection of one-electron equations, this approximation makes the issue easier to understand.

With these approximations, we can write the HF wavefunction as an antisymmetrized product of \(N\) single-electron wavefunctions, or spin-orbitals, called the Slater determinant:

Here, \(\phi_i(\mathbf{r}_j)\) represents the spin-orbital of the \(i\)-th electron at position \(\mathbf{r}_j\). The antisymmetrization ensures that the wavefunction obeys the Pauli exclusion principle.

The Hartree-Fock Equations

To find the optimal set of spin-orbitals, we need to minimize the energy expectation value, subject to the constraint that the spin-orbitals remain orthonormal. Using the method of Lagrange multipliers, we can derive the Hartree-Fock equations:

where \(\hat{f}_i\) is the Fock operator for electron \(i\), \(\phi_i(\mathbf{r})\) is the spin-orbital, and \(\epsilon_i\) is the orbital energy. The Fock operator consists of the one-electron Hamiltonian \(\hat{h}_i\) and the two-electron interaction terms:

where

and

The \(\hat{J}_j\) and \(\hat{K}_j\) operators are the Coulomb and exchange operators, respectively, representing the electron-electron interactions in the Hartree-Fock framework.

To solve the Hartree-Fock equations, we use an iterative process called the self-consistent field (SCF) method. Here is a step-by-step outline of the SCF procedure:

1. Choose an initial set of spin-orbitals and calculate the molecular integrals, which include the overlap matrix, the kinetic energy, the nuclear attraction, and the two-electron integrals.

2. Build the Fock matrix using the molecular integrals and the current spin-orbitals.

3. Solve the Hartree-Fock equations to obtain new spin-orbitals and orbital energies.

4. Check for convergence. If the energy and spin-orbitals have not converged, return to step 2 with the updated spin-orbitals. If converged, proceed to the next step.

5. Calculate the total Hartree-Fock energy using the converged spin-orbitals and orbital energies:

1. If needed, perform a population analysis to obtain charge and spin density distributions.

Molecular Integrals: Overlap, Kinetic Energy, and Nuclear Attraction

To solve the Hartree-Fock equations, we need to compute several molecular integrals, including the overlap integral, the kinetic energy integral, and the nuclear attraction integral.

The overlap integral measures the extent to which two atomic orbitals overlap:

The kinetic energy integral represents the kinetic energy of the electrons:

The nuclear attraction integral quantifies the attractive potential between the electrons and the nuclei:

By calculating these integrals and constructing the respective matrices, we can proceed with the SCF procedure to obtain the Hartree-Fock energy and wavefunction of the system.

Two-Electron Integrals

In addition to the overlap, kinetic energy, and nuclear attraction integrals, the Hartree-Fock method requires the evaluation of two-electron integrals, which account for the electron-electron repulsion. The two-electron integrals are given by:

Exploring with a bit of code

We will now set up a simple code for solving the Hartree-Fock equations:

# Parameters
R = 1.4  # Bond length in atomic units
Z = 1.0  # Atomic number
alpha = [3.42525091, 0.62391373, 0.16885540]  # STO-3G exponents

# Overlap integral
S(a, b, R) =  (pi / (a + b))^1.5 * exp(-a * b * R^2 / (a + b))

# Kinetic energy integral
T(a, b, R) = a * b * (3 - 2 * a * b * R^2) / (a + b)^2 * (pi / (a + b))^1.5 * exp(-a * b * R^2 / (a + b))

# Nuclear attraction integral
V(a, b, R, Z) = -2 * pi * Z / (a + b) * exp(-a * b * R^2 / (a + b))

# Two-electron integral
two_electron_integral(a, b, R) = (pi^(2.5) * (a * b)^(0.75) * exp(-a * b * R^2 / (a + b))) / (a + b)^(2.5)

# Calculate the core Hamiltonian
H = zeros(2, 2)
for i in 1:2, j in 1:2
    H[i, j] = T(alpha[i], alpha[j], R) + V(alpha[i], alpha[j], R, Z)
end

# Calculate the overlap matrix
S_matrix = zeros(2, 2)
for i in 1:2, j in 1:2
    S_matrix[i, j] = S(alpha[i], alpha[j], R)
end

function run_HF()
    # Calculate the two-electron integrals
    two_e = two_electron_integral(alpha[1], alpha[2], R)

    # Diagonalize the overlap matrix
    S_eigen = eigen(S_matrix)
    X = S_eigen.vectors * inv(Diagonal(sqrt.((S_eigen.values).+1))) * S_eigen.vectors'

    # Perform the Self-Consistent Field (SCF) procedure
    F = H
    C = X' * F * X
    F_prime = X * C * X'
    P = zeros(2, 2)
    E_old = 0.0
    tolerance = 1e-6
    max_iterations = 100
    for iteration in 1:max_iterations
        for i in 1:2, j in 1:2
            P[i, j] = C[i, 1] * C[j, 1]
        end
        G = zeros(2, 2)
        for i in 1:2, j in 1:2
            G[i, j] = P[i, 1] * P[1, j] * (two_e - 0.5 * two_e)
        end

        F = H + G
        C = X' * F * X
        F_prime = X * C * X'

        # Diagonalize the Fock matrix
        F_eigen = eigen(F_prime)
        C = X * F_eigen.vectors  # Update the coefficient matrix

        # Calculate the energy
        E = 0.0
        for i in 1:2, j in 1:2
            E += P[i, j] * (H[i, j] + F[i, j])
        end
        E += Z^2 / R

        # Check for convergence
        if abs(E - E_old) < tolerance
            break
        else
            E_old = E
        end
        if iteration == max_iterations
            println("SCF did not converge")
        end
        println("SCF converged after $$(iteration) iterations")
        println("Hartree-Fock energy: $$E atomic units")
    end
send

run_HF()

And we get the following output:

SCF converged after 1 iterations
Hartree-Fock energy: 0.15692369722853372 atomic units
SCF converged after 2 iterations
Hartree-Fock energy: -0.7505891906476986 atomic units
SCF converged after 3 iterations
Hartree-Fock energy: -0.7608913979870314 atomic units
SCF converged after 4 iterations
Hartree-Fock energy: -0.75979282480282 atomic units
SCF converged after 5 iterations
Hartree-Fock energy: -0.7599135205357536 atomic units
SCF converged after 6 iterations
Hartree-Fock energy: -0.7599003062831767 atomic units
SCF converged after 7 iterations
Hartree-Fock energy: -0.7599017535802811 atomic units

Conclusion

The Hartree-Fock method provides an accessible entry point into the world of electronic structure calculations. By making key approximations, it transforms the many-electron Schrödinger equation into a set of tractable one-electron equations. While it has limitations, the HF method lays the foundation for more sophisticated approaches that incorporate electron correlation effects such as ab-initio methods.