Spin Contamination in UHF

In unrestricted Hartree-Fock (UHF), the single determinant wavefunction is not an eigenfunction of the $\hat{S}^{2}$ operator. This leads to spin contamination — the wavefunction mixes in higher-spin states.

The $⟨ S^{2} ⟩$ Expectation Value

For a UHF determinant with $N_{α}$ alpha and $N_{β}$ beta electrons, the expectation value of the total spin operator is:

$⟨ S^{2} ⟩ = S (S + 1) + N_{β} - i = 1 \sum N_{α} j = 1 \sum N_{β} ∣ ⟨ ϕ_{i}^{α} ∣ ϕ_{j}^{β} ⟩ ∣^{2}$

where:

$S = \frac{N _{α} - N _{β}}{2}$ is the expected spin quantum number
$ϕ_{i}^{α}$ and $ϕ_{j}^{β}$ are the occupied alpha and beta molecular orbitals
The overlap integral is $⟨ ϕ_{i}^{α} ∣ ϕ_{j}^{β} ⟩ = \sum_{μν} C_{μ i}^{α} S_{μν} C_{ν j}^{β}$

Ideal vs. Actual Values

System	Multiplicity	$S$	Ideal $⟨ S^{2} ⟩ = S (S + 1)$
Doublet (1 unpaired e⁻)	2	½	0.75
Triplet (2 unpaired e⁻)	3	1	2.00
Quartet (3 unpaired e⁻)	4	3/2	3.75

The actual $⟨ S^{2} ⟩$ from a UHF calculation equals the ideal value only when the alpha and beta occupied orbitals are perfectly orthogonal ( $⟨ ϕ_{i}^{α} ∣ ϕ_{j}^{β} ⟩ = 0$ for all $i, j$ ). Any non-zero overlap between occupied alpha and beta orbitals increases $⟨ S^{2} ⟩$ above the ideal value.

Quantifying Spin Contamination

The degree of spin contamination is measured by the difference:

$Δ ⟨ S^{2} ⟩ = ⟨ S^{2} ⟩_{UHF} - ⟨ S^{2} ⟩_{ideal}$

$Δ ⟨ S^{2} ⟩ \approx 0$ : Pure spin state, minimal contamination
$0 < Δ ⟨ S^{2} ⟩ ≲ 0.1$ : Moderate contamination, usually acceptable
$Δ ⟨ S^{2} ⟩ > 0.1$ : Significant contamination; results may be unreliable

In OpenQuantum, a warning is printed when $Δ ⟨ S^{2} ⟩ > 0.01$ :

Spin Information:
  N(alpha): 5
  N(beta):  4
  <S²>:     0.850000 (ideal: 0.750000)
  WARNING: Significant spin contamination detected (0.1000)

Physical Interpretation

Spin contamination arises because the UHF determinant allows alpha and beta electrons to occupy different spatial orbitals. This flexibility lowers the energy but mixes in higher-spin configurations:

For a doublet ( $S = 1/2$ ), the contaminated wavefunction contains quartet ( $S = 3/2$ ) and higher-spin character
For a triplet ( $S = 1$ ), it contains quintet ( $S = 2$ ) character
The energy lowering is artificial — it comes from spin-symmetry breaking, not physical correlation

Consequences

Energies: UHF energies are variationally lower than RHF, but part of this lowering is due to spin contamination rather than physical correlation
Properties: Properties like dipole moments and gradients may be affected
Post-HF methods: MP2 and CCSD(T) built on a contaminated UHF reference inherit spin contamination
Geometry optimization: Contaminated forces can lead to incorrect geometries

Mitigation Strategies

Strategy	Description
ROHF	Restricted open-shell HF — forces alpha/beta orbitals to be identical in the open shell; spin-pure but often higher energy
Spin projection	Approximate projection of $⟨ S^{2} ⟩$ (e.g., Yamaguchi formula)
Spin-flip methods	Use spin-flip TDDFT or spin-flip CC to target spin-pure states
Use RHF for closed-shell	Avoid UHF entirely when not needed

Implementation in OpenQuantum

The computation is in crates/scf/src/fock.rs::compute_s_squared():

#![allow(unused)]
fn main() {
pub fn compute_s_squared(
    coefficients_alpha: &DMatrix<f64>,
    coefficients_beta: &DMatrix<f64>,
    overlap: &DMatrix<f64>,
    n_alpha: usize,
    n_beta: usize,
) -> (f64, f64) {
    // Ideal S²: S(S+1) where S = (N_α - N_β)/2
    let s = (n_alpha as f64 - n_beta as f64) / 2.0;
    let s_squared_ideal = s * (s + 1.0);

    // Overlap of occupied orbitals: S_αβ = C_occ^α^T · S · C_occ^β
    let c_alpha_occ = coefficients_alpha.columns(0, n_alpha);
    let c_beta_occ = coefficients_beta.columns(0, n_beta);
    let s_alpha_beta = c_alpha_occ.transpose() * overlap * c_beta_occ;

    // Sum of squared overlaps
    let overlap_sum: f64 = s_alpha_beta.iter().map(|x| x * x).sum();

    // ⟨S²⟩ = S_ideal + N_β - Σ|S_αβ|²
    let s_squared = s_squared_ideal + n_beta as f64 - overlap_sum;

    (s_squared, s_squared_ideal)
}
}

The result is returned in UhfResult and formatted in format_uhf_energy_output() with a warning for significant contamination.