5.2.5. Helium Atom Motivation#
For the sake of these notes we will consider a two electron Helium Hamiltonian of the form
\( \hat{H} = -\frac{1}{2}\nabla_1^2 - \frac{2}{r_1} - \frac{1}{2}\nabla_2^2 - \frac{2}{r_2} + \frac{1}{r_{12}} \),
where \(r_1\) is the distance of electron \(1\) from the nucleus, and \(r_{12}\) is the distance between the two electrons. The goal, as always, is to solve the Schrödinger equation with this Hamiltonian. This will yield the energies and wavefunctions of the system.
As we have seen, the Schrödinger equation with this Hamiltonian cannot be solved analytically. Instead, we have previously utilized both perturbation theory and the variational approach to solve for the approximate energies and wavefunctions. The Hartree-Fock approach is based heavily on the variational approach so will focus on that.
For the variational approach, we posit a normalized trial wavefunction with variational parameters. We previously used
\(\psi_t(r_1,r_2) = \psi_{1s}(r_1)\psi_{1s}(r_2) = \frac{\zeta^3}{\pi}e^{-\zeta(r_1+r_2)}\),
where \(\psi_{1s}(r) = \sqrt{\frac{\zeta^3}{\pi}}e^{-\zeta r}\) are normalized Hydrogen 1s-like spatial orbitals. The first step in the variational approach is to solve for the expectation value of the Hamiltonian in this trial wavefunction as a function of \(\zeta\), the variational parameter. We showed that
\(\langle E \rangle_t = \langle \psi_t | \hat{H} | \psi_t \rangle = \zeta^2 - \frac{27}{8}\zeta\)
We then minimize \(\langle E \rangle_t\) with respect to the variational parameter
\(\frac{d \langle E \rangle_t}{d\zeta}|_{\zeta_{min}} = 0\)
\(\Rightarrow \zeta_{min} = \frac{27}{16}\)
\(\Rightarrow E_{min} = -2.8477\) Hartree
In this procedure, we somewhat ignored the spin of the electron. This turns out to be fine for this specifc system, but moving to systems with larger number of electrons we have to be careful. Note that postulate 6 of quantum mechanics states that a wavefunction must be antisymmetric with respect to exhange of electrons (fermions). We account for this with the HF framework using slater determinants as described below.
5.2.6. Hartree-Fock#
The Hartree-Fock procedure is to effectively treat a multi-electron Schrödinger equation as a sum of independent electron Schrödinger equations with the electron-electron coupling term treated as an effective potential (in a mean field). This leads to a coupled set of equations that are solved iteretively in a self-consistent field (SCF).
For the Helium atom, we write the effective one-electron Hamiltonians as
Where \(V_1^{eff}(r_i)\) is the effective electron-electron repulsion potential. Technically, it is the potential felt by electron \(1\) at position \(r_1\) due to the density of electron 2 and is defined as
5.2.6.1. Antisymmetric Wavefunctions#
We need to ensure that our wavefunctions our antisymmetric with respect to exchange of electrons (postulate 6 of QM). This will be achieved using slater determinants of spin-orbits. For the minimal basis helium atom we will have two spin orbits
\(\chi_1(1) = \psi_1(1) \alpha(1)\)
\(\chi_2(1) = \psi_1(1) \beta(1)\)
where \(\psi_1\) is the spatial function and \(\alpha\) and \(\beta\) are the spin functions. The \(1\) in parantheses, \((1)\), denotes that these are functions of the coordinates of electron \(1\). Note that the spin coordinates and spacial coordinates are completely independent and that the spin functions are orthonormal: \(\langle \alpha | \alpha \rangle = \langle \beta | \beta \rangle = 1\) and \(\langle \alpha | \beta \rangle = \langle \beta | \alpha \rangle = 0\). We construct an antisymmetric wavefunction using Slater determinants:
\(|\Psi(1,2)\rangle = \frac{1}{\sqrt{2}} \begin{vmatrix}\chi_1(1) & \chi_1(2)\\ \chi_2(1) & \chi_2(2) \end{vmatrix} = \frac{1}{\sqrt{2}} (\chi_1(1)\chi_2(2) - \chi_2(1) \chi_1(2)) \)
This is generalized to \(N\) particles.
The HF method can be seen as the variational method for multielectron systems in which you consider the trial wavefunction as a separable, single Slater determinant.
5.2.6.2. General HF Equations#
The \(N\) electron Born-Oppenheimer Hamiltonian can be written as
5.2.6.3. Coulomb and Exhange Integrals#
We now want to determine the expectation value of the Hamiltonian in the antisymmetric trial wavefunction. Let’s start by rewriting the Hamiltonian as
where \(h(i) = -\frac{1}{2}\nabla_i^2 - \frac{2}{r_i}\) is the one electron operator for electron \(i\). Using this form of Hamiltonian, we can write the expectation value
where the last equality holds because there is no difference between the two electrons and two spatial orbitals in this problem.
We can now evaluate these terms separately.
The two electron integral is more interesting.
where the final equality holds due to commutative multiplication in this case. The first term in the last line is called a Coulomb integral, \(J\), and the second term is called an Exchange integral, \(K\).
Note that the exhange integral is zero for Helium, \(\langle\chi_1(1)\chi_2(2) |\frac{1}{r_{12}} | \chi_2(1)\chi_1(2) \rangle = 0\), because of the orthogonality of the spin functions.
Ultimately, for the Helium atom, this yields (for the He atom)
Differentiating this equation with respect to the orbitals yields the Fock equations:
where \(\hat f\) is the the Fock operator given by (for He):
and \(\varepsilon_i\) is the energy or orbital \(i\).
5.2.6.4. Self-consistent Field (SCF)#
To solve the Fock equations for orbital energies we must have the fock operator, which implicitly depends on the orbitals. Thus, this is a coupled problem that can be solved by iteration. We start by proposing the orbitals, create the Fock operator based on this guess (mean-field), and then solve the Fock equations for orbitals and orbital energies. The new orbitals are then used to generate a new, updated Fock operator and this procedure is iterated until convergence. This is called the self-consistent field method.
5.2.6.5. Roothan-Hall equations#
The above is all very general. If the spatial orbitals are expressed as a linear combination of functions, and the variational parameters are the linear coefficients, then the Hartree-Fock procedure becomes the Hartree-Fock Roothan-Hall procedure. This can be readily solved in matrix equations. We will go into these in more detail later (though I will use the results below).
5.2.7. Worked Example for the Helium Atom#
Consider the Helium atom where the overall spatial trial wavefunction is
where each function \(\psi(r)\) is a linear combination of 1s Slater-type orbitals with variational coefficients \(c_1\) and \(c_2\):
where
Note that we are not treating \(\zeta_1\) and \(\zeta_2\) as variational parameters rather they must be given. In this example we will use \(\zeta_1 = 1.45363\) and \(\zeta_2 = 2.91093\).
The variationally optimal \(c\) coefficients can be solved by determining the solution vector, \(c\), to the matrix equation
where \(\mathbf{S}\) is the basis overlap matrix and \(\mathbf{F}\) is the Fock matrix defined by
where \(\mathbf{H}\) is the core 1-electron Hamiltonian matrix and \(\mathbf{G}\) is the two electron effective potential with elements given by
and \(\mathbf{P}\) is the density matrix with elements defined by
\(G\) depends on the coefficients \(c\) so we have a coupled set of equations. The general procedure for solving these is:
Compute \(\mathbf{S}\) and \(\mathbf{H}\) - these don’t change over the course of an SCF calculation
Compute two electron integrals - these don’t change over the course of an SCF calculation
Guess coefficient vector
Compute \(\mathbf{G}\) and then \(\mathbf{F} = \mathbf{H} + \mathbf{G}\)
Diagonalize \(\mathbf{S}^{-1}\mathbf{F}\)
Update coefficient matrix and repeat steps 3-5 until convergence
5.2.7.1. Analytic integrals for same-center 1s STOs#
Let \(a,b,c,d\) index 1s STOs with exponents \(\zeta_a,\zeta_b,\zeta_c,\zeta_d\), and define
Overlap
Kinetic
Nuclear attraction
Two-electron repulsion
import numpy as np
def S_sto_1s(a, b):
"""Overlap <a|b> for normalized 1s STOs (same center)."""
return 8.0 * (a*b)**1.5 / (a + b)**3
def T_sto_1s(a, b):
"""Kinetic <a| -1/2 ∇^2 |b> for normalized 1s STOs (same center)."""
return 4.0 * (a*b)**2.5 / (a + b)**3
def V_sto_1s(a, b, Z):
"""Nuclear attraction <a| -Z/r |b> for normalized 1s STOs (same center)."""
return -4.0 * Z * (a*b)**1.5 / (a + b)**2
def H_core_sto_1s(a, b, Z):
return T_sto_1s(a, b) + V_sto_1s(a, b, Z)
def ERI_sto_1s(a, b, c, d):
"""
(ab|cd) for normalized same-center 1s STOs.
Corrected so that (ζζ|ζζ) = 5ζ/8.
"""
p = a + b
q = c + d
pref = 32.0 * (a*b*c*d)**1.5
num = p*p + 3.0*p*q + q*q
den = (p*p) * (q*q) * (p+q)**3
return pref * num / den
def build_integrals(zetas, Z):
"""Build S, H, and ERI tensor for same-center 1s STO basis."""
M = len(zetas)
S = np.zeros((M, M))
H = np.zeros((M, M))
eri = np.zeros((M, M, M, M))
for i in range(M):
for j in range(M):
S[i, j] = S_sto_1s(zetas[i], zetas[j])
H[i, j] = H_core_sto_1s(zetas[i], zetas[j], Z)
for i in range(M):
for j in range(M):
for k in range(M):
for l in range(M):
eri[i, j, k, l] = ERI_sto_1s(zetas[i], zetas[j], zetas[k], zetas[l])
return S, H, eri
def sym_orthogonalizer(S):
"""Return X = S^{-1/2} using symmetric orthogonalization."""
w, U = np.linalg.eigh(S)
if np.any(w <= 1e-12):
raise ValueError("S is singular/ill-conditioned.")
return U @ np.diag(w**-0.5) @ U.T
def build_fock_from_density(H, eri, P):
"""
Build Fock matrix from density P (closed shell RHF):
F_mn = H_mn + sum_ls P_ls [ (mn|ls) - 1/2 (ml|ns) ]
"""
M = H.shape[0]
F = H.copy()
for m in range(M):
for n in range(M):
g = 0.0
for l in range(M):
for s in range(M):
g += P[l, s] * (eri[m, n, l, s] - 0.5 * eri[m, l, n, s])
F[m, n] += g
return F
def rhf_energy(H, F, P):
"""Total RHF electronic energy: E = 1/2 * sum P (H + F)."""
return 0.5 * np.sum(P * (H + F))
import numpy as np
def normalize_c(c, S):
c = np.asarray(c, dtype=float).reshape(-1, 1)
n2 = (c.T @ S @ c).item()
if n2 <= 0:
raise ValueError("Non-positive norm under S.")
return c / np.sqrt(n2)
def enforce_positive_phase(c):
"""
Fix arbitrary sign of eigenvector: make the largest-magnitude coefficient positive.
"""
c = np.asarray(c).reshape(-1, 1)
idx = int(np.argmax(np.abs(c[:, 0])))
if c[idx, 0] < 0:
c = -c
return c
def rhf_scf_helium(zetas, Z=2.0, c_guess=None, max_iter=100, convE=1e-10, convP=1e-8, verbose=True):
S, H, eri = build_integrals(zetas, Z)
X = sym_orthogonalizer(S)
# ---- Initial guess for occupied MO coefficients ----
if c_guess is not None:
C_occ = enforce_positive_phase(normalize_c(c_guess, S)) # shape (M,1)
else:
# default: core Hamiltonian guess
Fp = X.T @ H @ X
eps, Cp = np.linalg.eigh(Fp)
C = X @ Cp
C_occ = enforce_positive_phase(C[:, [0]])
# initial density
P = 2.0 * (C_occ @ C_occ.T)
E_old = None
for it in range(1, max_iter + 1):
F = build_fock_from_density(H, eri, P)
if verbose and it==1:
E = rhf_energy(H, F, P)
c_print = C_occ.ravel()
print(f"iter {0:2d}: E = {E:+.12f} c = {c_print}")
Fp = X.T @ F @ X
eps, Cp = np.linalg.eigh(Fp)
C = X @ Cp
C_occ_new = enforce_positive_phase(C[:, [0]])
#new density
P_new = 2.0 * (C_occ_new @ C_occ_new.T)
# consistent energy with the *new* density
F_new = build_fock_from_density(H, eri, P_new)
E = rhf_energy(H, F_new, P_new)
dP = np.linalg.norm(P_new - P)
dE = np.inf if E_old is None else abs(E - E_old)
if verbose:
c_print = C_occ_new.ravel()
print(f"iter {it:2d}: E = {E:+.12f} c = {c_print} ||dP|| = {dP:.3e}")
if E_old is not None and dE < convE and dP < convP:
return E, eps[0], C_occ_new.ravel(), S, H, eri
E_old = E
P = P_new
raise RuntimeError("SCF did not converge.")
# --- Run with your exponents ---
zetas = np.array([1.45363, 2.91093])
E, eps0, c, S, H, eri = rhf_scf_helium(zetas, Z=2.0, c_guess=[0.5, 0.5], verbose=True)
print("\nFinal")
print("-----")
print(f"E_RHF = {E:+.12f} Hartree")
print(f"epsilon_1 = {eps0:+.12f} Hartree")
print(f"c = {c}")
print(f"c^T S c = {c @ (S @ c):.12f}")
iter 0: E = -2.605576861805 c = [0.52163718 0.52163718]
iter 1: E = -2.861227739088 c = [0.83091555 0.19508357] ||dP|| = 1.008e+00
iter 2: E = -2.861671533529 c = [0.84315634 0.18139176] ||dP|| = 4.956e-02
iter 3: E = -2.861672595174 c = [0.84375381 0.18072184] ||dP|| = 2.443e-03
iter 4: E = -2.861672597745 c = [0.8437832 0.18068888] ||dP|| = 1.203e-04
iter 5: E = -2.861672597751 c = [0.84378465 0.18068725] ||dP|| = 5.920e-06
iter 6: E = -2.861672597751 c = [0.84378472 0.18068717] ||dP|| = 2.914e-07
iter 7: E = -2.861672597751 c = [0.84378472 0.18068717] ||dP|| = 1.434e-08
iter 8: E = -2.861672597751 c = [0.84378472 0.18068717] ||dP|| = 7.060e-10
Final
-----
E_RHF = -2.861672597751 Hartree
epsilon_1 = -0.917935386525 Hartree
c = [0.84378472 0.18068717]
c^T S c = 1.000000000000