Hartree-Fock Method

The Hartree-Fock (HF) method is the foundation of modern electronic structure theory. It provides a mean-field approximation to the many-electron Schrödinger equation by approximating the wavefunction as a Slater determinant of one-electron orbitals.

This notebook covers:

  1. Basic Theory - The Hartree-Fock approximation

  2. Initial Guess - Methods for starting the SCF procedure

  3. Restricted Hartree-Fock (RHF) - For closed-shell systems

  4. Unrestricted Hartree-Fock (UHF) - For open-shell systems

  5. Restricted Open-shell Hartree-Fock (ROHF) - For open-shell systems with spatial symmetry

  6. Practical Examples - Using MiniHF’s HFSolver

1. Hartree-Fock Theory

The Many-Electron Problem

The time-independent Schrödinger equation for a molecule with \(N\) electrons is:

\[ \hat{H} \Psi = E \Psi \]

where the Hamiltonian operator is:

\[ \hat{H} = -\frac{1}{2} \sum_{i=1}^N \nabla_i^2 - \sum_{i=1}^N \sum_{\alpha=1}^M \frac{Z_{\alpha}}{r_{i\alpha}} + \sum_{i=1}^N \sum_{j>i}^N \frac{1}{r_{ij}} \]

The Hartree-Fock Approximation

The HF method approximates the many-electron wavefunction as a Slater determinant of spin-orbitals:

\[\begin{split} \Psi_{HF} = \frac{1}{\sqrt{N!}} \begin{vmatrix} \phi_1(\mathbf{r}_1, \sigma_1) & \phi_2(\mathbf{r}_1, \sigma_1) & \cdots & \phi_N(\mathbf{r}_1, \sigma_1) \\ \phi_1(\mathbf{r}_2, \sigma_2) & \phi_2(\mathbf{r}_2, \sigma_2) & \cdots & \phi_N(\mathbf{r}_2, \sigma_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_1(\mathbf{r}_N, \sigma_N) & \phi_2(\mathbf{r}_N, \sigma_N) & \cdots & \phi_N(\mathbf{r}_N, \sigma_N) \end{vmatrix} \end{split}\]

The Fock Matrix

Minimizing the energy leads to the Hartree-Fock equations:

\[ \hat{F} \phi_i = \epsilon_i \phi_i \]

where \(\hat{F}\) is the Fock operator:

\[ \hat{F}(\mathbf{r}) = \hat{h}(\mathbf{r}) + \sum_{j=1}^N \left[ \hat{J}_j(\mathbf{r}) - \hat{K}_j(\mathbf{r}) \right] \]

with:

  • \(\hat{h}(\mathbf{r})\): One-electron Hamiltonian (kinetic + nuclear attraction)

  • \(\hat{J}_j(\mathbf{r})\): Coulomb operator (classical electron repulsion)

  • \(\hat{K}_j(\mathbf{r})\): Exchange operator (quantum mechanical exchange)

Matrix Formulation

In a basis set \(\{\chi_{\mu}\}\), the HF equations become:

\[ \sum_{\nu} F_{\mu\nu} C_{\nu i} = \epsilon_i \sum_{\nu} S_{\mu\nu} C_{\nu i} \]

where:

  • \(F_{\mu\nu}\): Fock matrix element

  • \(S_{\mu\nu}\): Overlap matrix element

  • \(C_{\nu i}\): Molecular orbital coefficient

  • \(\epsilon_i\): Orbital energy

2. Initial Guess

The Hartree-Fock equations are solved iteratively through the Self-Consistent Field (SCF) procedure, which requires an initial guess for the molecular orbital coefficients. The quality of this initial guess matters significantly for three reasons:

  • Convergence: A good guess reduces the number of SCF iterations

  • Stability: A poor guess may lead to convergence to a local minimum

  • Performance: Some guesses work better for specific systems (e.g., transition metals)

Naive Initial Guess Methods

These methods are computationally inexpensive and easy to implement, but often yield poor convergence or converge to local minima. Thus, this project implements only the naive initial guess approaches (Core Hamiltonian, Hückel, and Extended Hückel) as described below.

2.1 Core Hamiltonian Guess

The simplest ab initio guess. It neglects all electron–electron repulsion and uses only one‑electron terms.

Procedure:

  1. Compute the core Hamiltonian matrix
    \(\mathbf{H}_{\text{core}} = \mathbf{T} + \mathbf{V}_{\text{ne}}\), where
    \(\mathbf{T}\) = kinetic energy, \(\mathbf{V}_{\text{ne}}\) = nuclear attraction.

  2. Solve the generalized eigenvalue problem
    \(\mathbf{H}_{\text{core}} \mathbf{C} = \mathbf{S} \mathbf{C} \epsilon\).

  3. Use the resulting coefficients \(\mathbf{C}\) as the initial MO guess.

Note: This guess is very crude. It contains no electron‑electron repulsion, so it often fails for metals, small gap systems, or when SCF is near instabilities.

2.2 Simple Hückel Guess (Non‑iterative)

A semi‑empirical π‑electron method. Not commonly used as a standalone HF initial guess today, but historically relevant.

Form (for planar π systems):

  • Set \(\mathbf{H}_{ii} = \alpha\), \(\mathbf{H}_{ij} = \beta\) (if bonded)

  • Solve \(\mathbf{H}^{\text{Hückel}} \mathbf{C} = \mathbf{S}^{\text{Hückel}} \mathbf{C} \epsilon\)
    (often \(\mathbf{S} = \mathbf{I}\) is assumed)

Limitation: Only π electrons, empirical parameters, not suitable for general HF initial guesses.

2.3 Extended Hückel Guess (EHT)

Uses empirical atom‑specific parameters (\(H_{ii}\) = valence orbital ionization energy) and a distance‑dependent off‑diagonal term (e.g., Wolfsberg–Helmholz formula).

Features:

  • Includes all valence electrons.

  • Constructs a modified \(\mathbf{H}^{\text{EHT}}\), then solves \(\mathbf{H}^{\text{EHT}} \mathbf{C} = \mathbf{S} \mathbf{C} \epsilon\).

  • Provides a much better initial guess than the core Hamiltonian, especially for organic molecules and transition metal complexes.

Advantage over core guess: It mimics some effective potential via empirical data, reducing SCF oscillations. Many modern codes (e.g., ORCA, XTB) use variants of EHT‑like guesses as default.

📋 Advanced Initial Guess
  • SAD (Superposition of Atomic Densities): Constructs the initial density matrix by superimposing spherically averaged atomic densities. Significantly better than the core Hamiltonian guess, especially for large basis sets and large molecules. Default guess in Q-Chem.
  • Harris: Diagonalizes the Harris functional to obtain the initial guess. Default option for HF and DFT calculations in Gaussian.
  • MINAO (Minimum Basis Atomic Orbitals): Builds the initial guess using a minimum basis set of atomic orbitals. Often serves as the underlying construction for other advanced guesses (e.g., SAD). Widely used in ORCA and PySCF.

3. Restricted Hartree-Fock (RHF)

RHF is used for closed-shell systems where all electrons are paired (even number of electrons, singlet spin state).

Key Features of RHF

  1. Spin Restriction: \(\alpha\) and \(\beta\) spin orbitals share the same spatial part:

    \[ \phi_i^{\alpha}(\mathbf{r}, \sigma) = \psi_i(\mathbf{r}) \alpha(\sigma) \]
    \[ \phi_i^{\beta}(\mathbf{r}, \sigma) = \psi_i(\mathbf{r}) \beta(\sigma) \]

  2. Density Matrix: A single density matrix describes both spins:

    \[ P_{\mu\nu} = 2 \sum_{i=1}^{N/2} C_{\mu i} C_{\nu i}^* \]

  3. Fock Matrix: Simplified due to spin restriction:

    \[ F_{\mu\nu} = H_{\mu\nu}^{core} + \sum_{\lambda\sigma} P_{\lambda\sigma} \left[ (\mu\nu|\lambda\sigma) - \frac{1}{2} (\mu\lambda|\nu\sigma) \right] \]

RHF SCF Algorithm

  1. Initialize: Compute integrals and form initial guess for \(\mathbf{C}\)

  2. Build Density: \(P_{\mu\nu} = 2 \sum_{i=1}^{N/2} C_{\mu i} C_{\nu i}^*\)

  3. Build Fock Matrix: \(\mathbf{F} = \mathbf{H}_{core} + \mathbf{J} - \frac{1}{2}\mathbf{K}\)

  4. Solve Roothaan Equations: \(\mathbf{F} \mathbf{C} = \mathbf{S} \mathbf{C} \mathbf{\epsilon}\)

  5. Compute New Density: \(\mathbf{P}^{new}\) from new \(\mathbf{C}\)

  6. Check Convergence: Compare \(\mathbf{P}^{new}\) with \(\mathbf{P}^{old}\)

  7. Iterate: Repeat until convergence

Convergence Criteria

  • Density convergence: \(||\mathbf{P}^{new} - \mathbf{P}^{old}|| < 10^{-8}\)

  • Energy convergence: \(|E^{new} - E^{old}| < 10^{-10}\) Hartree

  • Maximum iterations: Typically 50-200

4. Unrestricted Hartree-Fock (UHF)

UHF is used for open-shell systems where electrons are not all paired (odd number of electrons or triplet states).

Key Features of UHF

  1. Spin Unrestriction: \(\alpha\) and \(\beta\) spin orbitals have different spatial parts:

    \[ \phi_i^{\alpha}(\mathbf{r}, \sigma) = \psi_i^{\alpha}(\mathbf{r}) \alpha(\sigma) \]
    \[ \phi_i^{\beta}(\mathbf{r}, \sigma) = \psi_i^{\beta}(\mathbf{r}) \beta(\sigma) \]

  2. Separate Density Matrices: One for each spin:

    \[ P_{\mu\nu}^{\alpha} = \sum_{i=1}^{N_{\alpha}} C_{\mu i}^{\alpha} C_{\nu i}^{\alpha *} \]
    \[ P_{\mu\nu}^{\beta} = \sum_{i=1}^{N_{\beta}} C_{\mu i}^{\beta} C_{\nu i}^{\beta *} \]

  3. Separate Fock Matrices: One for each spin:

    \[ F_{\mu\nu}^{\alpha} = H_{\mu\nu}^{core} + \sum_{\lambda\sigma} (P_{\lambda\sigma}^{\alpha} + P_{\lambda\sigma}^{\beta}) (\mu\nu|\lambda\sigma) - \sum_{\lambda\sigma} P_{\lambda\sigma}^{\alpha} (\mu\lambda|\nu\sigma) \]
    \[ F_{\mu\nu}^{\beta} = H_{\mu\nu}^{core} + \sum_{\lambda\sigma} (P_{\lambda\sigma}^{\alpha} + P_{\lambda\sigma}^{\beta}) (\mu\nu|\lambda\sigma) - \sum_{\lambda\sigma} P_{\lambda\sigma}^{\beta} (\mu\lambda|\nu\sigma) \]

UHF vs RHF

Property

RHF

UHF

Spin Restriction

\(\alpha=\beta\)

\(\alpha\neq\beta\)

Density Matrices

1

2

Fock Matrices

1

2

Applicable System

Closed-Shell

Open-Shell

Spin Contamination

None

May occur

Spin Contamination

UHF wavefunctions can be contaminated with higher spin states:

\[ \langle S^2 \rangle = S(S+1) + \delta S \]

where \(\delta S > 0\) is the spin contamination.

5. Restricted Open-shell Hartree-Fock (ROHF)

ROHF is a compromise between RHF and UHF for open-shell systems.

Key Features of ROHF

  1. Partial Spin Restriction: Core orbitals are doubly occupied with same spatial function, virtual orbitals are unoccupied. Only the open-shell orbitals have different \(\alpha\) and \(\beta\) occupation.

  2. Avoids Spin Contamination: Maintains spatial symmetry while handling open-shell character.

  3. More Complex Equations: Requires separate treatment of core, open-shell, and virtual orbitals.

ROHF Energy Expression

\[ E_{ROHF} = \sum_{i}^{core} 2 \langle \psi_i | \hat{h} | \psi_i \rangle + \sum_{i}^{open} \langle \psi_i^{\alpha} | \hat{h} | \psi_i^{\alpha} \rangle + \sum_{i}^{open} \langle \psi_i^{\beta} | \hat{h} | \psi_i^{\beta} \rangle \]

plus Coulomb and exchange terms.

When to Use ROHF

  • Open-shell systems where spin contamination is a concern

  • Systems with high symmetry

  • Transition metal complexes with unpaired electrons

6. Generalized Hartree-Fock (GHF)

GHF is the most general form of Hartree-Fock theory.

Key Features of GHF

  1. No Spin Restrictions: Orbitals can be arbitrary linear combinations of \(\alpha\) and \(\beta\) spins:

    \[ \phi_i(\mathbf{r}, \sigma) = \psi_i^{\alpha}(\mathbf{r}) \alpha(\sigma) + \psi_i^{\beta}(\mathbf{r}) \beta(\sigma) \]

  2. Complex Orbitals: GHF orbitals are generally complex-valued.

  3. Most Flexible: Can describe any spin state.

GHF Density Matrix

The density matrix becomes a 2x2 matrix in spin space:

\[\begin{split} \mathbf{P} = \begin{pmatrix} \mathbf{P}^{\alpha\alpha} & \mathbf{P}^{\alpha\beta} \\ \mathbf{P}^{\beta\alpha} & \mathbf{P}^{\beta\beta} \end{pmatrix} \end{split}\]

When to Use GHF

  • Systems with strong spin-orbit coupling

  • Systems where spatial symmetry is broken

  • Theoretical studies of spin dynamics

Drawbacks of GHF

  • Computationally more expensive

  • Complex-valued orbitals

  • Less commonly implemented

7. Practical Examples with MiniHF

Let’s demonstrate the Hartree-Fock solver in MiniHF with practical examples.

# Import required modules
from minihf.molecule import Molecule
from minihf.basis import BasisSet
from minihf.hamiltonian import Hamiltonian
from minihf.scf.hf_solver import HFSolver
import numpy as np

# Example 1: RHF calculation for water (closed-shell)
print("\n--- RHF Calculation: Water Molecule ---")

# Create water molecule
h2o_atoms = ['O 0.0, 0.0, 0.0', 'H 0.0, -0.757, 0.587', 'H 0.0, 0.757, 0.587']
h2o = Molecule.build(atoms=h2o_atoms, charge=0, multiplicity=1)

# Build basis set
h2o_basis = BasisSet.build(h2o, 'sto-3g.nwchem')

# Build Hamiltonian
h2o_ham = Hamiltonian.build(h2o, h2o_basis)

# Create HF solver and run RHF
hf_solver = HFSolver(h2o_ham)
rhf_results = hf_solver.rhf()

# Print results
print(f"Molecule: Water (H2O)")
print(f"Basis set: STO-3G")
print(f"Number of basis functions: {len(h2o_basis)}")
print(f"\nTotal energy: {rhf_results['energy_total']:.6f} Hartree")
print(f"Electronic energy: {rhf_results['energy_elec']:.6f} Hartree")
print(f"Nuclear repulsion: {rhf_results['energy_nuc']:.6f} Hartree")

# Print orbital energies
print("\nMolecular orbital energies (eV):")
for i, energy in enumerate(rhf_results['mo_energy']):
    print(f"  MO {i+1}: {energy * 27.2114:>8.4f} eV")

# HOMO-LUMO gap
homo_idx = 4  # Water has 10 electrons = 5 occupied MOs (0-indexed: 0-4)
homo_energy = rhf_results['mo_energy'][homo_idx] * 27.2114
lumo_energy = rhf_results['mo_energy'][homo_idx + 1] * 27.2114
print(f"\nHOMO energy: {homo_energy:.4f} eV")
print(f"LUMO energy: {lumo_energy:.4f} eV")
print(f"HOMO-LUMO gap: {lumo_energy - homo_energy:.4f} eV")

--- RHF Calculation: Water Molecule ---

Molecule: Water (H2O)
Basis set: STO-3G
Number of basis functions: 7

Total energy: -73.218632 Hartree
Electronic energy: -90.581924 Hartree
Nuclear repulsion: 17.363292 Hartree

Molecular orbital energies (eV):
  MO 1: -556.8357 eV
  MO 2: -45.9198 eV
  MO 3: -25.6964 eV
  MO 4: -15.9000 eV
  MO 5: -14.8850 eV
  MO 6:  26.9861 eV
  MO 7:  40.8066 eV

HOMO energy: -14.8850 eV
LUMO energy: 26.9861 eV
HOMO-LUMO gap: 41.8711 eV

# Example 2: UHF calculation for methyl radical (open-shell)
print("\n--- UHF Calculation: Methyl Radical ---")

# Create methyl radical (CH3)
ch3_atoms = ['C 0.0, 0.0, 0.0', 'H 0.629, 0.629, 0.629', 'H -0.629, -0.629, 0.629', 'H -0.629, 0.629, -0.629']
ch3 = Molecule.build(atoms=ch3_atoms, charge=0, multiplicity=2)  # Doublet state

# Build basis set
ch3_basis = BasisSet.build(ch3, 'sto-3g.nwchem')

# Build Hamiltonian
ch3_ham = Hamiltonian.build(ch3, ch3_basis)

# Create HF solver and run UHF
hf_solver = HFSolver(ch3_ham)
uhf_results = hf_solver.uhf()

# Print results
print(f"Molecule: Methyl Radical (CH3)")
print(f"Basis set: STO-3G")
print(f"Multiplicity: {ch3.multiplicity} (doublet)")
print(f"Number of basis functions: {len(ch3_basis)}")
print(f"\nTotal energy: {uhf_results['energy_total']:.6f} Hartree")
print(f"Electronic energy: {uhf_results['energy_elec']:.6f} Hartree")
print(f"Nuclear repulsion: {uhf_results['energy_nuc']:.6f} Hartree")

# Compare alpha and beta orbital energies
print("\nOrbital energies (eV):")
print(f"{'MO':<5} {'Alpha':<12} {'Beta':<12} {'Difference':<12}")
print("-" * 45)
for i, (e_alpha, e_beta) in enumerate(zip(uhf_results['mo_energy'][0], uhf_results['mo_energy'][1])):
    diff = (e_alpha - e_beta) * 27.2114
    print(f"{i+1:<5} {e_alpha * 27.2114:<12.4f} {e_beta * 27.2114:<12.4f} {diff:<12.4f}")

# Check spin contamination
if 's_squared' in uhf_results:
    s_squared = uhf_results['s_squared']
    expected = (ch3.multiplicity - 1) * ch3.multiplicity / 4  # S(S+1) for doublet = 0.75
    print(f"\n<S^2> = {s_squared:.6f}")
    print(f"Expected = {expected:.6f}")
    print(f"Spin contamination = {s_squared - expected:.6f}")

--- UHF Calculation: Methyl Radical ---

Molecule: Methyl Radical (CH3)
Basis set: STO-3G
Multiplicity: 2 (doublet)
Number of basis functions: 8

Total energy: -37.146975 Hartree
Electronic energy: -55.355186 Hartree
Nuclear repulsion: 18.208211 Hartree

Orbital energies (eV):
MO    Alpha        Beta         Difference  
---------------------------------------------
1     -311.4240    -310.6645    -0.7594     
2     -34.2493     -32.0569     -2.1924     
3     -22.5441     -21.5784     -0.9657     
4     -22.5441     -21.5784     -0.9657     
5     -16.0028     3.3759       -19.3788    
6     39.3965      40.4194      -1.0229     
7     39.3965      40.4194      -1.0229     
8     51.8817      53.5871      -1.7053

8. Summary

In this notebook, we’ve covered the main variants of the Hartree-Fock method:

Method Comparison

Method

Description

Best For

RHF

Restricted Hartree-Fock

Closed-shell systems (even electrons, singlet)

UHF

Unrestricted Hartree-Fock

Open-shell systems (odd electrons, triplets)

ROHF

Restricted Open-shell HF

Open-shell systems with symmetry constraints

GHF

Generalized Hartree-Fock

Systems with strong spin-orbit coupling

Key Takeaways

  1. SCF Procedure: HF equations are solved iteratively until self-consistency

  2. Initial Guess: Important for convergence; core Hamiltonian guess is common

  3. RHF: Most efficient for closed-shell systems

  4. UHF: Handles open-shell but may have spin contamination

  5. ROHF: Balances symmetry and open-shell treatment

  6. GHF: Most general but computationally expensive