The H2 binding curve

Quantum Chemistry

Everything is in place. h2_scf(R) takes a bond length and returns a Hartree-Fock total energy. Calling it across a range of $R$ traces out the potential energy curve, which is the only thing this whole series has been working toward. This page runs that scan, locates the minimum, and looks at where Hartree-Fock gets the answer right and where it spectacularly fails.

The code (new this page)

Two short routines: one that computes the curve at a list of $R$ values, one that finds the minimum by scipy.optimize.minimize_scalar and reports the binding energy.

def h2_binding_curve(R_values):
    """Run h2_scf at each R and return arrays of (R, E_total, E_elec)."""
    R_values = np.asarray(R_values)
    E_total  = np.zeros_like(R_values)
    E_elec   = np.zeros_like(R_values)
    for i, R in enumerate(R_values):
        Et, _, _ = h2_scf(R)
        E_total[i] = Et
        S, T, V = h2_one_electron_matrices(R)
        ERI = h2_eri_tensor(R)
        Ee, _, _, _ = scf(T + V, S, ERI, n_occ=1, verbose=False)
        E_elec[i] = Ee
    return R_values, E_total, E_elec


def find_equilibrium():
    """Locate the minimum of E_total(R) and report (R_eq, E_eq, D_e)."""
    res = minimize_scalar(lambda R: h2_scf(R)[0], bracket=(1.0, 1.5, 2.0))
    R_eq, E_eq = res.x, res.fun

    # binding energy: D_e = 2 E(H atom, STO-3G) - E(H2 minimum)
    E_H = hydrogen_energy_sto3g()
    D_e = 2 * E_H - E_eq
    return R_eq, E_eq, D_e

Run it

H2 binding curve (STO-3G, RHF):

  R / bohr       E_elec        1/R         E_total
   0.500       -2.403326   +2.000000    -0.403326
   0.700       -2.265702   +1.428571    -0.837130
   1.000       -2.065999   +1.000000    -1.065999
   1.200       -1.943667   +0.833333    -1.110334
   1.300       -1.886102   +0.769231    -1.116871
   1.346       -1.861330   +0.743176    -1.117506   <- minimum
   1.400       -1.831000   +0.714286    -1.116714
   1.500       -1.778363   +0.666667    -1.111696
   1.700       -1.680241   +0.588235    -1.092006
   2.000       -1.549171   +0.500000    -1.049171
   2.500       -1.365794   +0.400000    -0.965794
   3.000       -1.218608   +0.333333    -0.885275
   4.000       -1.011082   +0.250000    -0.761082
   6.000       -0.811743   +0.166667    -0.645077
   8.000       -0.735039   +0.125000    -0.610039

Equilibrium:    R_eq = 1.345919 bohr  (0.7123 angstrom)
                E_eq = -1.117506 Hartree
                D_e  = +0.184342 Hartree  (5.016 eV)

Reference:      experimental R_eq    = 1.401 bohr (0.7414 A)
                experimental D_e     = 0.1745 Hartree (4.747 eV)
                exact-basis HF D_e   = 0.1336 Hartree (3.636 eV)

The shape of the curve

Three regimes, easy to read off the table.

Short bond, $R ≲ 1$ bohr. The nuclear-nuclear repulsion $1/ R$ dominates and the energy shoots up. The electronic part keeps getting more negative as the electrons feel both nuclei more strongly, but it can't keep up with the $1/ R$ divergence. The wall on the left of the curve is the Pauli + nuclear-repulsion barrier that stops the molecule from collapsing.

Equilibrium, $R_{e q} = 1.346$ bohr. The energy minimum sits at $- 1.1175$ Hartree. Compared to two non-interacting STO-3G hydrogen atoms ( $2 E_{H} = - 0.9332$ ), the molecule is bound by $D_{e} = 0.184$ Hartree — about 5 eV. The experimental bond length is $1.401$ bohr; STO-3G is a touch too short by 4 percent. The experimental dissociation energy is $0.1745$ Hartree, which means STO-3G modestly over-binds when measured against its own atomic reference, and under-binds when measured against the exact $- 1$ Hartree per atom (the STO-3G H-atom energy itself is 33 mHa too high — the basis is missing the cusp).

Dissociation, $R \to \infty$ . The energy at $R = 8$ bohr is $- 0.610$ Hartree. The correct answer is two free hydrogen atoms, $2 E_{H} = - 0.933$ in STO-3G — a difference of $323$ millihartree. Restricted Hartree-Fock does not dissociate H₂ correctly. This is the most famous failure of single-determinant methods, and it's the moment to look at it carefully because every quantum-chemistry method ever invented since 1960 has had to address it.

Why RHF can't dissociate H₂

The RHF wavefunction is a single Slater determinant with both electrons in the bonding molecular orbital. At equilibrium, that orbital is roughly $σ_{g} \propto ϕ_{A} + ϕ_{B}$ — a symmetric combination of the two atomic 1s functions. As $R \to \infty$ , the bonding MO becomes $σ_{g} \to (ϕ_{A} + ϕ_{B}) / 2$ , and the two-electron wavefunction is

ψ (r_{1}, r_{2}) \propto [ϕ_{A} (r_{1}) + ϕ_{B} (r_{1})] [ϕ_{A} (r_{2}) + ϕ_{B} (r_{2})] \cdot (singlet spin) .

Expand the product:

ψ \propto ionic: H^{-} + H^{+} ϕ_{A} (1) ϕ_{A} (2) + ϕ_{B} (1) ϕ_{B} (2) + covalent: H + H ϕ_{A} (1) ϕ_{B} (2) + ϕ_{B} (1) ϕ_{A} (2) .

Half ionic, half covalent. The covalent half is the right answer: one electron on each atom, energy $2 E_{H}$ . The ionic half is wrong: at infinite separation, putting both electrons on the same proton costs $E (H^{-}) + E (H^{+}) \approx - 0.5 + 0 = - 0.5$ Hartree, which is half a Hartree higher than the covalent answer. The 50/50 mixture sits halfway between, around $- 0.75$ Hartree, which is exactly the order-of-magnitude error we're seeing in the table at $R = 8$ .

The fix requires letting the wavefunction be a sum of two Slater determinants — one bonding-bonding, one antibonding-antibonding — with the relative coefficient varying from 0 at equilibrium to $- 1$ at infinity, killing the ionic contribution. That's two-configuration self-consistent field (TCSCF), or more generally any of full CI, CASSCF, MR-CI. Single-determinant HF can't do it. UHF (unrestricted Hartree-Fock) can recover the correct dissociation limit by spin-symmetry-breaking — letting the alpha and beta orbitals localise on different atoms — at the cost of producing a wavefunction that is no longer a spin eigenfunction.

For our purposes — a single-determinant series ending at H₂ — the right take is: HF is excellent near equilibrium, where the system is well-described by a single closed-shell determinant; HF is structurally wrong at dissociation, and the wrongness is large enough to be visible in even this small calculation. Knowing the failure mode is almost as important as getting the equilibrium answer right.

What we actually built

Eleven steps ago we started with kinetic energy of a 1s Slater orbital — a one-line function that looked like a parlor trick. The progression: variational principle, Gaussian basis set, Slater determinants, the Hartree-Fock equations, Roothaan-Hall matrix form, SCF iteration, multi-centre integrals over Gaussians. No single step was hard. The cumulative effect is a from-scratch quantum chemistry program that reproduces Szabo & Ostlund's canonical worked example for H₂ and locates a binding minimum that's within 4 percent of experiment, in about 200 lines of Python.

The same code, with multi-angular-momentum primitive integrals instead of just $s$ -type, would do hydrocarbons. With DIIS instead of plain iteration, it would converge transition metals. With density fitting and Schwarz screening, it would scale to dozens of atoms. The architecture doesn't change. What we have is the kernel that every electronic-structure code in chemistry is built on. The next page lists what we left out and where to look for it.

Is the over-binding (relative to STO-3G atoms) physical or a basis-set artifact?

Largely a basis-set artifact. STO-3G is bad at the H atom (the contraction can't reproduce the cusp) and gets less bad in H₂ at equilibrium (the second atom's basis functions improve the description of the first atom's electron, by a phenomenon called basis-set superposition error). The result is that the molecule looks artificially stable relative to its constituent atoms in the same basis. Larger basis sets reduce this effect; the counterpoise correction tries to remove it explicitly.

Why is the equilibrium bond length 4% short?

Two compounding effects, both basis-related. (1) STO-3G under-describes the bonding region between the atoms — its tails are wrong — and the SCF compensates by pulling the atoms closer. (2) Without correlation, HF tends to slightly over-bind closed-shell single bonds at equilibrium (it's missing dispersion, but for short bonds the correlation hole contributes too). Bigger basis sets and post-HF corrections push $R_{e q}$ back out toward the experimental value.

How does UHF recover the dissociation limit?

UHF allows the alpha and beta spin orbitals to be different functions, rather than forcing them to share a spatial part. At long $R$ , the energy minimum is found at a symmetry-broken solution: alpha electron localises on atom A, beta electron localises on atom B (or vice versa). At that configuration the wavefunction is essentially $ϕ_{A} (1) α ϕ_{B} (2) β$ , which has no ionic component, so the energy goes correctly to $2 E_{H}$ . The price is that this UHF wavefunction has $⟨ S^{2} ⟩ = 1$ instead of $0$ for a singlet — it's a spin contamination, half singlet and half triplet. Spin-projected UHF and multireference methods clean this up. The takeaway: HF only works when the closed-shell determinant is a good description, and that condition fails at dissociation.