Notebooks · posts/h2-rhf-uhf-and-full-ci-along-the-bond

RHF, UHF and full CI for H₂ in STO-3G along the bond

Adam Zettel · September 25, 2026 · written with an agent

Pull the two atoms of H₂ apart in the smallest basis there is, and compute its energy three ways. Restricted Hartree-Fock is off by $0.021$ Ha at the equilibrium bond length and heads for the wrong limit, 0.39 Ha too high. Unrestricted Hartree-Fock splits from it at $R = 2.180$ bohr and reaches two hydrogen atoms. Full CI is exact in this basis. Every energy matches PySCF to $1 0^{- 9}$ Ha.

E_{corr} = E_{FCI} - E_{RHF} = - 0.02056 Ha (R = 1.4 bohr)

The program is self-contained: the STO-3G basis, the integrals over Gaussian functions, the Hartree-Fock and full-CI energies, all in this file. The reference is PySCF 2.5.0 (posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.py), which shares no code with it: its own basis library, integrals, SCF solvers, stability analysis and full-CI solver. The Coulson-Fischer point is checked against PySCF two ways, by the onset of the RHF triplet instability in PySCF's molecular-orbital integrals (2.179695 bohr) and by bisecting on where PySCF's UHF first falls below its RHF (2.179731, with a 10⁻¹⁰ Ha threshold, which puts it slightly past the true onset).

Units are atomic: energies in hartree (1 Ha = 27.211 eV), lengths in bohr (1 bohr = 0.5292 Å), both from CODATA 2022.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

The basis

A hydrogen molecule has two electrons and two protons. Its electrons' wave function has to be built out of something, and the smallest choice is one function per atom. STO-3G makes that function from three Gaussians, $e^{- α r^{2}}$ with fixed exponents $α$ and fixed weights, fitted once to look like a hydrogen 1s orbital. Gaussians are used because every integral the calculation needs has a closed form for them.

The two functions, $χ_{1}$ on atom 1 and $χ_{2}$ on atom 2, give two molecular orbitals, fixed by symmetry: the bonding orbital $σ_{g} \propto χ_{1} + χ_{2}$ and the antibonding orbital $σ_{u} \propto χ_{1} - χ_{2}$ .

[1]lines 42–81

PI = 3.141592653589793
# STO-3G for hydrogen, as PySCF's basis library holds it (posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out prints it)
ALPHA = [3.42525091, 0.62391373, 0.1688554]
COEF = [0.15432897, 0.53532814, 0.44463454]
# one basis function on each atom
NBF = 2

#| F₀(t) = ∫₀¹ exp(−t·u²) du = ½·√(π/t)·erf(√t),
#| and its series 1 − t/3 at tiny t.
def boys0(t) {
    if t < 1.0e-12 {
        return 1.0 - t / 3.0
    }
    return 0.5 * sqrt(PI / t) * erf(sqrt(t))
}

#| The normalization of the primitive exp(−a·r²): (2a/π)^(3/4).
def prim_norm(a) {
    return pow(2.0 * a / PI, 0.75)
}

#| The overlap of two normalized s primitives with exponents a and b whose
#| centers are a distance squared r2 apart.
def s_prim(a, b, r2) {
    p = a + b
    return prim_norm(a) * prim_norm(b) * pow(PI / p, 1.5) * exp(-a * b / p * r2)
}

#| The factor that makes the contracted function ∑ᵢ COEF[i]·gᵢ normalized.
def contraction_scale() {
    n2 = 0.0
    for i to len(ALPHA) {
        for j to len(ALPHA) {
            n2 += COEF[i] * COEF[j] * s_prim(ALPHA[i], ALPHA[j], 0.0)
        }
    }
    return 1.0 / sqrt(n2)
}

D_SCALE = contraction_scale()

PI = 3.141592653589793
# STO-3G for hydrogen, as PySCF's basis library holds it (posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out prints it)
ALPHA = [3.42525091, 0.62391373, 0.1688554]
COEF = [0.15432897, 0.53532814, 0.44463454]
# one basis function on each atom
NBF = 2

#| F₀(t) = ∫₀¹ exp(−t·u²) du = ½·√(π/t)·erf(√t),
#| and its series 1 − t/3 at tiny t.
def boys0(t) {
    if t < 1.0e-12 {
        return 1.0 - t / 3.0
    }
    return 0.5 * sqrt(PI / t) * erf(sqrt(t))
}

#| The normalization of the primitive exp(−a·r²): (2a/π)^(3/4).
def prim_norm(a) {
    return pow(2.0 * a / PI, 0.75)
}

#| The overlap of two normalized s primitives with exponents a and b whose
#| centers are a distance squared r2 apart.
def s_prim(a, b, r2) {
    p = a + b
    return prim_norm(a) * prim_norm(b) * pow(PI / p, 1.5) * exp(-a * b / p * r2)
}

#| The factor that makes the contracted function ∑ᵢ COEF[i]·gᵢ normalized.
def contraction_scale() {
    n2 = 0.0
    for i to len(ALPHA) {
        for j to len(ALPHA) {
            n2 += COEF[i] * COEF[j] * s_prim(ALPHA[i], ALPHA[j], 0.0)
        }
    }
    return 1.0 / sqrt(n2)
}

D_SCALE = contraction_scale()

The integrals

With the atoms at $z = 0$ and $z = R$ , the program needs the overlap $S_{pq}$ , the one-electron energy $h_{pq}$ (kinetic energy plus attraction to both protons), and the electron-repulsion integrals $(pq ∣ r s)$ , each a sum over three primitives per function. For two s Gaussians with exponents $a$ and $b$ , the product is one Gaussian of exponent $p = a + b$ centered between them, and the standard results (Szabo and Ostlund, Modern Quantum Chemistry, Appendix A) are

S = (\frac{π}{p})^{3/2} e^{- \frac{ab}{p} R_{A B}^{2}}, T = \frac{ab}{p} (3 - \frac{2 ab}{p} R_{A B}^{2}) S, V_{C} = - \frac{2 π}{p} e^{- \frac{ab}{p} R_{A B}^{2}} F_{0} (p R_{P C}^{2}),

and for the repulsion between two such products, with exponents $p$ and $q$ and centers $P$ and $Q$ , $(ab ∣ c d) = \frac{2 π ^{5/2}}{pq p + q} e^{- \frac{ab}{p} R_{A B}^{2} - \frac{c d}{q} R_{C D}^{2}} F_{0} (\frac{pq}{p + q} R_{P Q}^{2})$ , each times the primitives' normalizations.

[2]lines 104–187

#| The AO integrals of H₂ at bond length R, written into S, h and eri; eri
#| holds (pq|rs) at index ((p·2 + q)·2 + r)·2 + s.
def ao_integrals(R, S: Arr, h: Arr, eri: Arr) {
    z = [0.0, R]
    np = len(ALPHA)
    for p to NBF {
        for q to NBF {
            S[p, q] = 0.0
            h[p, q] = 0.0
            for i to np {
                for j to np {
                    a = ALPHA[i]
                    b = ALPHA[j]
                    w = D_SCALE * D_SCALE * COEF[i] * COEF[j] * prim_norm(a) * prim_norm(b)
                    pab = a + b
                    rab2 = (z[p] - z[q]) * (z[p] - z[q])
                    kab = exp(-a * b / pab * rab2)
                    zp = (a * z[p] + b * z[q]) / pab
                    ov = pow(PI / pab, 1.5) * kab
                    S[p, q] += w * ov
                    h[p, q] += w * a * b / pab * (3.0 - 2.0 * a * b / pab * rab2) * ov
                    for c to NBF {
                        h[p, q] -= w * 2.0 * PI / pab * kab * boys0(pab * (zp - z[c]) * (zp - z[c]))
                    }
                }
            }
        }
    }
    for p to NBF {
        for q to NBF {
            for r to NBF {
                for s to NBF {
                    acc = 0.0
                    for i to np {
                        for j to np {
                            for k to np {
                                for l to np {
                                    a = ALPHA[i]
                                    b = ALPHA[j]
                                    c = ALPHA[k]
                                    d = ALPHA[l]
                                    w = COEF[i] * COEF[j] * COEF[k] * COEF[l] * prim_norm(a) * prim_norm(b) * prim_norm(c) * prim_norm(d)
                                    pp = a + b
                                    qq = c + d
                                    zp = (a * z[p] + b * z[q]) / pp
                                    zq = (c * z[r] + d * z[s]) / qq
                                    e_ab = a * b / pp * (z[p] - z[q]) * (z[p] - z[q])
                                    e_cd = c * d / qq * (z[r] - z[s]) * (z[r] - z[s])
                                    t = pp * qq / (pp + qq) * (zp - zq) * (zp - zq)
                                    acc += w * 2.0 * pow(PI, 2.5) / (pp * qq * sqrt(pp + qq)) * exp(-e_ab - e_cd) * boys0(t)
                                }
                            }
                        }
                    }
                    eri[((p * 2.0 + q) * 2.0 + r) * 2.0 + s] = pow(D_SCALE, 4.0) * acc
                }
            }
        }
    }
    return 0.0
}

# the number of AO repulsion integrals (pq|rs)
N_ERI = NBF * NBF * NBF * NBF
R_EQ = 1.4
S_eq = zeros(NBF, NBF)
h_eq = zeros(NBF, NBF)
eri_eq = zeros(N_ERI)
ao_integrals(R_EQ, S_eq, h_eq, eri_eq)
show S_eq
show h_eq
show eri_eq

test "the AO integrals at R = 1.4 are PySCF's to 1e-10" {
    # posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out, from mol.intor
    assert_near(S_eq[0, 0], 1.0, 1.0e-10)
    assert_near(S_eq[0, 1], 0.659318206134864, 1.0e-10)
    assert_near(h_eq[0, 0], -1.1204090089068204, 1.0e-10)
    assert_near(h_eq[0, 1], -0.958379964389617, 1.0e-10)
    assert_near(eri_eq[0], 0.7746059439198978, 1.0e-10)
    assert_near(eri_eq[3], 0.5696759256037501, 1.0e-10)
    assert_near(eri_eq[5], 0.2970285402769315, 1.0e-10)
    assert_near(eri_eq[1], 0.44410765803196095, 1.0e-10)
}

#| The AO integrals of H₂ at bond length R, written into S, h and eri; eri
#| holds (pq|rs) at index ((p·2 + q)·2 + r)·2 + s.
def ao_integrals(R, S: Arr, h: Arr, eri: Arr) {
    z = [0.0, R]
    np = len(ALPHA)
    for p to NBF {
        for q to NBF {
            S[p, q] = 0.0
            h[p, q] = 0.0
            for i to np {
                for j to np {
                    a = ALPHA[i]
                    b = ALPHA[j]
                    w = D_SCALE * D_SCALE * COEF[i] * COEF[j] * prim_norm(a) * prim_norm(b)
                    pab = a + b
                    rab2 = (z[p] - z[q]) * (z[p] - z[q])
                    kab = exp(-a * b / pab * rab2)
                    zp = (a * z[p] + b * z[q]) / pab
                    ov = pow(PI / pab, 1.5) * kab
                    S[p, q] += w * ov
                    h[p, q] += w * a * b / pab * (3.0 - 2.0 * a * b / pab * rab2) * ov
                    for c to NBF {
                        h[p, q] -= w * 2.0 * PI / pab * kab * boys0(pab * (zp - z[c]) * (zp - z[c]))
                    }
                }
            }
        }
    }
    for p to NBF {
        for q to NBF {
            for r to NBF {
                for s to NBF {
                    acc = 0.0
                    for i to np {
                        for j to np {
                            for k to np {
                                for l to np {
                                    a = ALPHA[i]
                                    b = ALPHA[j]
                                    c = ALPHA[k]
                                    d = ALPHA[l]
                                    w = COEF[i] * COEF[j] * COEF[k] * COEF[l] * prim_norm(a) * prim_norm(b) * prim_norm(c) * prim_norm(d)
                                    pp = a + b
                                    qq = c + d
                                    zp = (a * z[p] + b * z[q]) / pp
                                    zq = (c * z[r] + d * z[s]) / qq
                                    e_ab = a * b / pp * (z[p] - z[q]) * (z[p] - z[q])
                                    e_cd = c * d / qq * (z[r] - z[s]) * (z[r] - z[s])
                                    t = pp * qq / (pp + qq) * (zp - zq) * (zp - zq)
                                    acc += w * 2.0 * pow(PI, 2.5) / (pp * qq * sqrt(pp + qq)) * exp(-e_ab - e_cd) * boys0(t)
                                }
                            }
                        }
                    }
                    eri[((p * 2.0 + q) * 2.0 + r) * 2.0 + s] = pow(D_SCALE, 4.0) * acc
                }
            }
        }
    }
    return 0.0
}

# the number of AO repulsion integrals (pq|rs)
N_ERI = NBF * NBF * NBF * NBF
R_EQ = 1.4
S_eq = zeros(NBF, NBF)
h_eq = zeros(NBF, NBF)
eri_eq = zeros(N_ERI)
ao_integrals(R_EQ, S_eq, h_eq, eri_eq)
show S_eq
show h_eq
show eri_eq

test "the AO integrals at R = 1.4 are PySCF's to 1e-10" {
    # posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out, from mol.intor
    assert_near(S_eq[0, 0], 1.0, 1.0e-10)
    assert_near(S_eq[0, 1], 0.659318206134864, 1.0e-10)
    assert_near(h_eq[0, 0], -1.1204090089068204, 1.0e-10)
    assert_near(h_eq[0, 1], -0.958379964389617, 1.0e-10)
    assert_near(eri_eq[0], 0.7746059439198978, 1.0e-10)
    assert_near(eri_eq[3], 0.5696759256037501, 1.0e-10)
    assert_near(eri_eq[5], 0.2970285402769315, 1.0e-10)
    assert_near(eri_eq[1], 0.44410765803196095, 1.0e-10)
}

The overlap of the two atoms' functions at $R = 1.4$ bohr is $0.659$ (S_eq): they share most of their space. $(11∣11) = 0.775$ Ha is the repulsion of two electrons both in $χ_{1}$ , and it matters below.

Three approximations to the same energy

Restricted Hartree-Fock (RHF) puts both electrons, spin up and spin down, in the same orbital. For H₂ in this basis that orbital must be $σ_{g}$ , so nothing needs iterating, and the energy is $E_{RHF} = 2 h_{g g} + (g g ∣ g g) + 1/ R$ in molecular-orbital integrals.

Unrestricted Hartree-Fock (UHF) lets the spin-up electron have one orbital and the spin-down electron another. Here they can only be mixtures of $σ_{g}$ and $σ_{u}$ ; the lowest-energy pair has the form $cos θ σ_{g} \pm sin θ σ_{u}$ , which leans the spin-up electron toward one atom and the spin-down electron toward the other. The program minimizes the energy over $θ$ . At $θ = 0$ it is RHF.

Full configuration interaction (FCI) is the exact answer within the basis: the lowest eigenvalue of the Hamiltonian over every two-electron state the two orbitals allow. With the right symmetry only two states mix, $σ_{g}^{2}$ and $σ_{u}^{2}$ , so it is a 2 by 2 matrix $(E_{g g} K K E_{uu})$ with $K = (g u ∣ g u)$ .

[3]lines 216–314

#| The MO integral (ab|cd) from the AO ones, for orbital coefficient
#| vectors a, b, c and d.
def mo_eri(eri: Arr, a: Arr, b: Arr, c: Arr, d: Arr) {
    acc = 0.0
    for p to NBF {
        for q to NBF {
            for r to NBF {
                for s to NBF {
                    acc += a[p] * b[q] * c[r] * d[s] * eri[((p * 2.0 + q) * 2.0 + r) * 2.0 + s]
                }
            }
        }
    }
    return acc
}

#| The one-electron integral ⟨a|h|b⟩.
def mo_h(h: Arr, a: Arr, b: Arr) {
    acc = 0.0
    for p to NBF {
        for q to NBF {
            acc += a[p] * h[p, q] * b[q]
        }
    }
    return acc
}

#| The UHF energy with spin orbitals cos θ·σg ± sin θ·σu.
def uhf_energy(h: Arr, eri: Arr, g: Arr, u: Arr, theta, R) {
    al = cos(theta) * g + sin(theta) * u
    be = cos(theta) * g - sin(theta) * u
    return mo_h(h, al, al) + mo_h(h, be, be) + mo_eri(eri, al, al, be, be) + 1.0 / R
}

# golden-section steps for the UHF angle; 0.618^60 of pi/4 is far below
# where the energy stops changing
N_GOLDEN = 60

#| RHF, UHF and FCI at bond length R: [E_RHF, E_UHF, E_FCI, <S^2>_UHF].
def h2_energies(R) {
    S = zeros(NBF, NBF)
    h = zeros(NBF, NBF)
    eri = zeros(N_ERI)
    ao_integrals(R, S, h, eri)
    g = [1.0, 1.0] / sqrt(2.0 * (1.0 + S[0, 1]))
    u = [1.0, -1.0] / sqrt(2.0 * (1.0 - S[0, 1]))
    e_gg = 2.0 * mo_h(h, g, g) + mo_eri(eri, g, g, g, g)
    e_uu = 2.0 * mo_h(h, u, u) + mo_eri(eri, u, u, u, u)
    k = mo_eri(eri, g, u, g, u)
    e_rhf = e_gg + 1.0 / R
    e_fci = 0.5 * (e_gg + e_uu) - sqrt(0.25 * (e_uu - e_gg) * (e_uu - e_gg) + k * k) + 1.0 / R
    gr = 0.5 * (sqrt(5.0) - 1.0)
    lo = 0.0
    hi = 0.25 * PI
    for it to N_GOLDEN {
        t1 = hi - gr * (hi - lo)
        t2 = lo + gr * (hi - lo)
        if uhf_energy(h, eri, g, u, t1, R) < uhf_energy(h, eri, g, u, t2, R) {
            hi = t2
        }
        else {
            lo = t1
        }
    }
    theta = 0.5 * (lo + hi)
    e_uhf = uhf_energy(h, eri, g, u, theta, R)
    if e_uhf > e_rhf {
        e_uhf = e_rhf
        theta = 0.0
    }
    s2 = sin(2.0 * theta) * sin(2.0 * theta)
    return [e_rhf, e_uhf, e_fci, s2]
}

at_eq = h2_energies(R_EQ)
e_corr = at_eq[2] - at_eq[0]
show at_eq
show e_corr

test "at R = 1.4 RHF and FCI are PySCF's, and so is the correlation energy, -0.020562 Ha" {
    # posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out: scf.RHF and fci.FCI at conv_tol 1e-12
    assert_near(at_eq[0], -1.116714325062551, 1.0e-9)
    assert_near(at_eq[1], at_eq[0], 1.0e-12)
    assert_near(at_eq[2], -1.1372759436170439, 1.0e-9)
    assert_near(e_corr, -0.02056161855449279, 1.0e-9)
}

# the correlation energy as a share of the full CI energy, and in eV
# (the hartree in eV: CODATA 2022, physics.nist.gov/cgi-bin/cuu/Value?hrev)
HARTREE_EV = 27.211386245981
corr_percent = 100.0 * e_corr / at_eq[2]
corr_ev = -e_corr * HARTREE_EV
show corr_percent
show corr_ev

test "the correlation energy is 1.8% of the total energy and 0.56 eV" {
    assert_near(corr_percent, 1.808, 0.001)
    assert_near(corr_ev, 0.5595, 0.0001)
}

#| The MO integral (ab|cd) from the AO ones, for orbital coefficient
#| vectors a, b, c and d.
def mo_eri(eri: Arr, a: Arr, b: Arr, c: Arr, d: Arr) {
    acc = 0.0
    for p to NBF {
        for q to NBF {
            for r to NBF {
                for s to NBF {
                    acc += a[p] * b[q] * c[r] * d[s] * eri[((p * 2.0 + q) * 2.0 + r) * 2.0 + s]
                }
            }
        }
    }
    return acc
}

#| The one-electron integral ⟨a|h|b⟩.
def mo_h(h: Arr, a: Arr, b: Arr) {
    acc = 0.0
    for p to NBF {
        for q to NBF {
            acc += a[p] * h[p, q] * b[q]
        }
    }
    return acc
}

#| The UHF energy with spin orbitals cos θ·σg ± sin θ·σu.
def uhf_energy(h: Arr, eri: Arr, g: Arr, u: Arr, theta, R) {
    al = cos(theta) * g + sin(theta) * u
    be = cos(theta) * g - sin(theta) * u
    return mo_h(h, al, al) + mo_h(h, be, be) + mo_eri(eri, al, al, be, be) + 1.0 / R
}

# golden-section steps for the UHF angle; 0.618^60 of pi/4 is far below
# where the energy stops changing
N_GOLDEN = 60

#| RHF, UHF and FCI at bond length R: [E_RHF, E_UHF, E_FCI, <S^2>_UHF].
def h2_energies(R) {
    S = zeros(NBF, NBF)
    h = zeros(NBF, NBF)
    eri = zeros(N_ERI)
    ao_integrals(R, S, h, eri)
    g = [1.0, 1.0] / sqrt(2.0 * (1.0 + S[0, 1]))
    u = [1.0, -1.0] / sqrt(2.0 * (1.0 - S[0, 1]))
    e_gg = 2.0 * mo_h(h, g, g) + mo_eri(eri, g, g, g, g)
    e_uu = 2.0 * mo_h(h, u, u) + mo_eri(eri, u, u, u, u)
    k = mo_eri(eri, g, u, g, u)
    e_rhf = e_gg + 1.0 / R
    e_fci = 0.5 * (e_gg + e_uu) - sqrt(0.25 * (e_uu - e_gg) * (e_uu - e_gg) + k * k) + 1.0 / R
    gr = 0.5 * (sqrt(5.0) - 1.0)
    lo = 0.0
    hi = 0.25 * PI
    for it to N_GOLDEN {
        t1 = hi - gr * (hi - lo)
        t2 = lo + gr * (hi - lo)
        if uhf_energy(h, eri, g, u, t1, R) < uhf_energy(h, eri, g, u, t2, R) {
            hi = t2
        }
        else {
            lo = t1
        }
    }
    theta = 0.5 * (lo + hi)
    e_uhf = uhf_energy(h, eri, g, u, theta, R)
    if e_uhf > e_rhf {
        e_uhf = e_rhf
        theta = 0.0
    }
    s2 = sin(2.0 * theta) * sin(2.0 * theta)
    return [e_rhf, e_uhf, e_fci, s2]
}

at_eq = h2_energies(R_EQ)
e_corr = at_eq[2] - at_eq[0]
show at_eq
show e_corr

test "at R = 1.4 RHF and FCI are PySCF's, and so is the correlation energy, -0.020562 Ha" {
    # posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out: scf.RHF and fci.FCI at conv_tol 1e-12
    assert_near(at_eq[0], -1.116714325062551, 1.0e-9)
    assert_near(at_eq[1], at_eq[0], 1.0e-12)
    assert_near(at_eq[2], -1.1372759436170439, 1.0e-9)
    assert_near(e_corr, -0.02056161855449279, 1.0e-9)
}

# the correlation energy as a share of the full CI energy, and in eV
# (the hartree in eV: CODATA 2022, physics.nist.gov/cgi-bin/cuu/Value?hrev)
HARTREE_EV = 27.211386245981
corr_percent = 100.0 * e_corr / at_eq[2]
corr_ev = -e_corr * HARTREE_EV
show corr_percent
show corr_ev

test "the correlation energy is 1.8% of the total energy and 0.56 eV" {
    assert_near(corr_percent, 1.808, 0.001)
    assert_near(corr_ev, 0.5595, 0.0001)
}

At the equilibrium bond length the two Hartree-Fock methods agree, at $- 1.11671$ Ha, and full CI is lower by $0.02056$ Ha (e_corr): the correlation energy, the part of the electrons' repulsion a single orbital per electron cannot describe. It is 1.8% of the total energy and 0.56 eV, the same size as the energy differences chemistry cares about.

Pulling the atoms apart

Before running: at large $R$ each atom should keep one electron, and the energy should be twice that of a hydrogen atom in this basis, $2 \times (- 0.46658) = - 0.93316$ Ha. Which of the three methods gets there?

[4]lines 330–392

# the hydrogen atom in STO-3G, from PySCF's UHF of a lone H (posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out)
E_H_ATOM = -0.46658184955727533
# bond lengths along the curve (bohr)
N_CURVE = 56
R_LO = 0.6
R_HI = 8.0
R_axis = linspace(R_LO, R_HI, N_CURVE)
curve_rhf = zeros(N_CURVE)
curve_uhf = zeros(N_CURVE)
curve_fci = zeros(N_CURVE)
curve_s2 = zeros(N_CURVE)
for k to N_CURVE {
    e = h2_energies(R_axis[k])
    curve_rhf[k] = e[0]
    curve_uhf[k] = e[1]
    curve_fci[k] = e[2]
    curve_s2[k] = e[3]
}

R_CHECK = [0.8, 2.0, 2.5, 3.0, 4.0, 6.0, 10.0]
REF_RHF = [-0.9473079311160366, -1.0491709019858069, -0.9657936760289826, -0.8852750001110465, -0.761082247024236, -0.6450767495106331, -0.5959706348512762]
REF_UHF = [-0.9473079311160366, -1.0491709019858069, -0.9799475298244157, -0.9510179480528731, -0.9358423283198218, -0.9332273454117433, -0.9331637039744968]
REF_FCI = [-0.9575985320217821, -1.0884963081468997, -1.0304740010738502, -0.9851568243679045, -0.9437784716239705, -0.9334021769910897, -0.9331637120042972]
REF_S2 = [0.0, 0.0, 0.439823149071, 0.774240341562, 0.963991844797, 0.999297214498, 0.999999973454]
worst_e = 0.0
worst_s2 = 0.0
for k to len(R_CHECK) {
    e = h2_energies(R_CHECK[k])
    worst_e = max(worst_e, abs(e[0] - REF_RHF[k]))
    worst_e = max(worst_e, abs(e[1] - REF_UHF[k]))
    worst_e = max(worst_e, abs(e[2] - REF_FCI[k]))
    worst_s2 = max(worst_s2, abs(e[3] - REF_S2[k]))
}
far = h2_energies(10.0)
rhf_excess = far[0] - 2.0 * E_H_ATOM
show worst_e
show worst_s2
show rhf_excess

test "RHF, UHF and FCI match PySCF to 1e-9 Ha at seven bond lengths from 0.8 to 10 bohr" {
    assert_lt(worst_e, 1.0e-9)
}

test "the UHF spin contamination <S^2> matches PySCF's spin_square to 1e-7" {
    # the energy is flat at its minimum in theta, so theta is found only to
    # about 1e-8 (the square root of the energy's rounding), and <S^2> =
    # sin^2(2 theta) moves by up to four times that; measured worst: 2.4e-8
    assert_lt(worst_s2, 1.0e-7)
}

test "at R = 10 UHF and FCI are two hydrogen atoms, and RHF is 0.337 Ha too high" {
    assert_near(far[1], 2.0 * E_H_ATOM, 1.0e-6)
    assert_near(far[2], 2.0 * E_H_ATOM, 1.0e-6)
    # PySCF: -0.5959706348512762 - 2 (-0.46658184955727533)
    assert_near(rhf_excess, 0.33719306426327446, 1.0e-9)
}

figure("H2 in STO-3G: energy against bond length", "bond length R (bohr)", "energy (hartree)", R_LO, R_HI, -1.2, -0.5)
plot(R_axis, curve_rhf, "RHF")
plot(R_axis, curve_uhf, "UHF")
style("UHF", "dash")
plot(R_axis, curve_fci, "full CI")
hline("two H atoms: -0.9332", 2.0 * E_H_ATOM)

# the hydrogen atom in STO-3G, from PySCF's UHF of a lone H (posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out)
E_H_ATOM = -0.46658184955727533
# bond lengths along the curve (bohr)
N_CURVE = 56
R_LO = 0.6
R_HI = 8.0
R_axis = linspace(R_LO, R_HI, N_CURVE)
curve_rhf = zeros(N_CURVE)
curve_uhf = zeros(N_CURVE)
curve_fci = zeros(N_CURVE)
curve_s2 = zeros(N_CURVE)
for k to N_CURVE {
    e = h2_energies(R_axis[k])
    curve_rhf[k] = e[0]
    curve_uhf[k] = e[1]
    curve_fci[k] = e[2]
    curve_s2[k] = e[3]
}

R_CHECK = [0.8, 2.0, 2.5, 3.0, 4.0, 6.0, 10.0]
REF_RHF = [-0.9473079311160366, -1.0491709019858069, -0.9657936760289826, -0.8852750001110465, -0.761082247024236, -0.6450767495106331, -0.5959706348512762]
REF_UHF = [-0.9473079311160366, -1.0491709019858069, -0.9799475298244157, -0.9510179480528731, -0.9358423283198218, -0.9332273454117433, -0.9331637039744968]
REF_FCI = [-0.9575985320217821, -1.0884963081468997, -1.0304740010738502, -0.9851568243679045, -0.9437784716239705, -0.9334021769910897, -0.9331637120042972]
REF_S2 = [0.0, 0.0, 0.439823149071, 0.774240341562, 0.963991844797, 0.999297214498, 0.999999973454]
worst_e = 0.0
worst_s2 = 0.0
for k to len(R_CHECK) {
    e = h2_energies(R_CHECK[k])
    worst_e = max(worst_e, abs(e[0] - REF_RHF[k]))
    worst_e = max(worst_e, abs(e[1] - REF_UHF[k]))
    worst_e = max(worst_e, abs(e[2] - REF_FCI[k]))
    worst_s2 = max(worst_s2, abs(e[3] - REF_S2[k]))
}
far = h2_energies(10.0)
rhf_excess = far[0] - 2.0 * E_H_ATOM
show worst_e
show worst_s2
show rhf_excess

test "RHF, UHF and FCI match PySCF to 1e-9 Ha at seven bond lengths from 0.8 to 10 bohr" {
    assert_lt(worst_e, 1.0e-9)
}

test "the UHF spin contamination <S^2> matches PySCF's spin_square to 1e-7" {
    # the energy is flat at its minimum in theta, so theta is found only to
    # about 1e-8 (the square root of the energy's rounding), and <S^2> =
    # sin^2(2 theta) moves by up to four times that; measured worst: 2.4e-8
    assert_lt(worst_s2, 1.0e-7)
}

test "at R = 10 UHF and FCI are two hydrogen atoms, and RHF is 0.337 Ha too high" {
    assert_near(far[1], 2.0 * E_H_ATOM, 1.0e-6)
    assert_near(far[2], 2.0 * E_H_ATOM, 1.0e-6)
    # PySCF: -0.5959706348512762 - 2 (-0.46658184955727533)
    assert_near(rhf_excess, 0.33719306426327446, 1.0e-9)
}

figure("H2 in STO-3G: energy against bond length", "bond length R (bohr)", "energy (hartree)", R_LO, R_HI, -1.2, -0.5)
plot(R_axis, curve_rhf, "RHF")
plot(R_axis, curve_uhf, "UHF")
style("UHF", "dash")
plot(R_axis, curve_fci, "full CI")
hline("two H atoms: -0.9332", 2.0 * E_H_ATOM)

Only RHF misses. At $R = 10$ bohr it is $0.337$ Ha above two atoms (rhf_excess). The reason is in its one orbital. $σ_{g}^{2}$ written out in the atoms' functions is $\frac{1}{2} (χ_{1} χ_{2} + χ_{2} χ_{1}) + \frac{1}{2} (χ_{1} χ_{1} + χ_{2} χ_{2})$ : half the time the two electrons are on different atoms, and half the time they are both on one, H⁻ next to a bare proton. Once the atoms' functions stop overlapping, that ionic half costs the repulsion $(11∣11)$ of two electrons on one atom and gains the attraction $1/ R$ between H⁻ and the proton, each weighted by one half. Before looking at the next cell: what should the excess be at $R = 20$ bohr, and where does it go as $R$ grows?

[5]lines 407–432

R_FAR = 20.0
R_LIMIT = 1000.0
far20 = h2_energies(R_FAR)
excess_20 = far20[0] - 2.0 * E_H_ATOM
ionic_half = 0.5 * eri_eq[0] - 0.5 / R_FAR
excess_limit = h2_energies(R_LIMIT)[0] - 2.0 * E_H_ATOM
show excess_20
show ionic_half
show excess_limit

test "at R = 20 the RHF excess is PySCF's 0.362303, which is (11|11)/2 - 1/(2R)" {
    # posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out: E_RHF(20) - 2 E(H) = 0.3623029719599441, and
    # (11|11) = 0.7746059439198978, so (11|11)/2 - 1/40 = 0.3623029719599489
    assert_near(excess_20, 0.3623029719599441, 1.0e-9)
    assert_near(excess_20, 0.5 * 0.7746059439198978 - 0.5 / R_FAR, 1.0e-9)
}

# where the excess goes: half of the repulsion of two electrons on one atom
ionic_limit = 0.5 * eri_eq[0]
show ionic_limit

test "at R = 1000 the excess is 0.3868, on its way to (11|11)/2 = 0.3873 Ha" {
    # (11|11) = 0.7746059439198978 from PySCF, posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out
    assert_near(ionic_limit, 0.5 * 0.7746059439198978, 1.0e-10)
    assert_near(excess_limit, ionic_limit - 0.5 / R_LIMIT, 1.0e-6)
}

R_FAR = 20.0
R_LIMIT = 1000.0
far20 = h2_energies(R_FAR)
excess_20 = far20[0] - 2.0 * E_H_ATOM
ionic_half = 0.5 * eri_eq[0] - 0.5 / R_FAR
excess_limit = h2_energies(R_LIMIT)[0] - 2.0 * E_H_ATOM
show excess_20
show ionic_half
show excess_limit

test "at R = 20 the RHF excess is PySCF's 0.362303, which is (11|11)/2 - 1/(2R)" {
    # posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out: E_RHF(20) - 2 E(H) = 0.3623029719599441, and
    # (11|11) = 0.7746059439198978, so (11|11)/2 - 1/40 = 0.3623029719599489
    assert_near(excess_20, 0.3623029719599441, 1.0e-9)
    assert_near(excess_20, 0.5 * 0.7746059439198978 - 0.5 / R_FAR, 1.0e-9)
}

# where the excess goes: half of the repulsion of two electrons on one atom
ionic_limit = 0.5 * eri_eq[0]
show ionic_limit

test "at R = 1000 the excess is 0.3868, on its way to (11|11)/2 = 0.3873 Ha" {
    # (11|11) = 0.7746059439198978 from PySCF, posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out
    assert_near(ionic_limit, 0.5 * 0.7746059439198978, 1.0e-10)
    assert_near(excess_limit, ionic_limit - 0.5 / R_LIMIT, 1.0e-6)
}

At $R = 20$ bohr the excess is $0.36230$ Ha (excess_20), which is $\frac{1}{2} (11∣11) - \frac{1}{2 R}$ to nine digits, and at $R = 1000$ it is $0.3868$ (excess_limit), on its way to $\frac{1}{2} (11∣11) = 0.3873$ Ha. RHF describes the separated atoms as half atoms and half ions.

Where UHF breaks away

UHF avoids this by letting the spin-up electron sit on one atom and the spin-down electron on the other. Close to equilibrium that costs more than it gains, and UHF's lowest state is the RHF one. Past a bond length called the Coulson-Fischer point it pays, and the two curves split. The program finds it as the $R$ where the UHF energy stops curving upward in $θ$ at $θ = 0$ : the second difference $E (δ) + E (- δ) - 2 E (0)$ changes sign.

The price is the spin. The UHF state is no longer a pure singlet: its $⟨ S^{2} ⟩$ , zero for a singlet and 2 for a triplet, is $sin^{2} 2 θ$ , and it rises to 1 as the atoms separate, an even mixture of the two.

[6]lines 456–502

# the step in theta for the second difference, and the bisection steps on R
DELTA = 1.0e-3
N_BISECT = 40

#| The second difference of the UHF energy in theta at theta = 0, at bond
#| length R: positive while RHF is stable against UHF.
def uhf_curvature(R) {
    S = zeros(NBF, NBF)
    h = zeros(NBF, NBF)
    eri = zeros(N_ERI)
    ao_integrals(R, S, h, eri)
    g = [1.0, 1.0] / sqrt(2.0 * (1.0 + S[0, 1]))
    u = [1.0, -1.0] / sqrt(2.0 * (1.0 - S[0, 1]))
    e0 = uhf_energy(h, eri, g, u, 0.0, R)
    return uhf_energy(h, eri, g, u, DELTA, R) + uhf_energy(h, eri, g, u, -DELTA, R) - 2.0 * e0
}

lo = 1.5
hi = 3.5
for it to N_BISECT {
    mid = 0.5 * (lo + hi)
    if uhf_curvature(mid) > 0.0 {
        lo = mid
    }
    else {
        hi = mid
    }
}
r_cf = 0.5 * (lo + hi)
show r_cf

test "the Coulson-Fischer point is PySCF's RHF triplet instability, 2.179695 bohr, to 1e-5" {
    # posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out: brentq on e_u - e_g - (gg|uu) - (gu|gu) in PySCF's MO integrals
    assert_near(r_cf, 2.179694793840183, 1.0e-5)
}

cf_ratio = r_cf / R_EQ
show cf_ratio

test "the Coulson-Fischer point is 1.56 times R = 1.4 bohr" {
    assert_near(cf_ratio, 2.179694793840183 / 1.4, 1.0e-5)
}

figure("UHF spin contamination", "bond length R (bohr)", "$\\langle S^2 \\rangle$ of the UHF state", R_LO, R_HI, 0.0, 1.05)
plot(R_axis, curve_s2, "UHF $\\langle S^2 \\rangle$")
vline("Coulson-Fischer point", r_cf)
hline("even singlet-triplet mixture", 1.0)

# the step in theta for the second difference, and the bisection steps on R
DELTA = 1.0e-3
N_BISECT = 40

#| The second difference of the UHF energy in theta at theta = 0, at bond
#| length R: positive while RHF is stable against UHF.
def uhf_curvature(R) {
    S = zeros(NBF, NBF)
    h = zeros(NBF, NBF)
    eri = zeros(N_ERI)
    ao_integrals(R, S, h, eri)
    g = [1.0, 1.0] / sqrt(2.0 * (1.0 + S[0, 1]))
    u = [1.0, -1.0] / sqrt(2.0 * (1.0 - S[0, 1]))
    e0 = uhf_energy(h, eri, g, u, 0.0, R)
    return uhf_energy(h, eri, g, u, DELTA, R) + uhf_energy(h, eri, g, u, -DELTA, R) - 2.0 * e0
}

lo = 1.5
hi = 3.5
for it to N_BISECT {
    mid = 0.5 * (lo + hi)
    if uhf_curvature(mid) > 0.0 {
        lo = mid
    }
    else {
        hi = mid
    }
}
r_cf = 0.5 * (lo + hi)
show r_cf

test "the Coulson-Fischer point is PySCF's RHF triplet instability, 2.179695 bohr, to 1e-5" {
    # posts/refs/h2-rhf-uhf-and-full-ci-along-the-bond.out: brentq on e_u - e_g - (gg|uu) - (gu|gu) in PySCF's MO integrals
    assert_near(r_cf, 2.179694793840183, 1.0e-5)
}

cf_ratio = r_cf / R_EQ
show cf_ratio

test "the Coulson-Fischer point is 1.56 times R = 1.4 bohr" {
    assert_near(cf_ratio, 2.179694793840183 / 1.4, 1.0e-5)
}

figure("UHF spin contamination", "bond length R (bohr)", "$\\langle S^2 \\rangle$ of the UHF state", R_LO, R_HI, 0.0, 1.05)
plot(R_axis, curve_s2, "UHF $\\langle S^2 \\rangle$")
vline("Coulson-Fischer point", r_cf)
hline("even singlet-triplet mixture", 1.0)

The split comes at $R = 2.180$ bohr (r_cf), 1.56 times the 1.4 bohr bond length used above. Past it the UHF curve bends over toward two atoms, and $⟨ S^{2} ⟩$ climbs, to 0.44 at 2.5 bohr and 0.96 at 4 bohr. Full CI stays a pure singlet and goes to the right limit smoothly: it mixes in $σ_{u}^{2}$ , which cancels the ionic half of $σ_{g}^{2}$ , a little near equilibrium and half and half at large $R$ .

[7]lines 513–514

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors