# Restricted Hartree–Fock in C++, from scratch

A closed-shell (RHF) Hartree–Fock code with **no external dependencies** — the
molecular integrals, the linear algebra, the generalized eigensolver, and the
SCF loop with DIIS are all written from scratch. No libint, no LAPACK, no Eigen.

Validated on H₂/STO-3G at R = 1.401 bohr against Szabo & Ostlund:

```
total energy:      -1.1166856295 Hartree     (Szabo: -1.1167)
orbital energies:  -0.578, +0.670 Ha         (Szabo: -0.578, +0.670)
```

## What's inside

| file | purpose |
|---|---|
| `src/boys.cpp`      | Boys function F_n(t): series for small t, asymptotic + downward recursion elsewhere |
| `src/cart_int.cpp`  | Obara–Saika recurrences for Gaussian overlap/kinetic/nuclear-attraction/ERI over Cartesian primitives |
| `src/integrals.cpp` | assembles S, T, V, and the two-electron integral tensor over the contracted basis |
| `src/basis.cpp`     | STO-3G contractions, H₂ dimer builder |
| `src/matrix.cpp`    | dense matrix type, Cholesky, Jacobi eigensolver, generalized symmetric eigenproblem FC = SCε |
| `src/hf.cpp`        | Fock build, RHF energy, density, DIIS, the SCF loop |
| `tests/`            | numeric checks: known 2×2 eigenvalues, CᵀSC = I, H₂ core-Hamiltonian eigenvalues |

## Build & run

```
mkdir build && cd build
cmake -DCMAKE_BUILD_TYPE=Release .. && make
./h2_sto3g        # the validation run above
./qchem_tests     # unit checks
```

## The bug that H₂ couldn't see

Found while preparing this code for hosting (2026-07-06). `build_fock` had

```cpp
jk += P(l,s) * (2.0*(uv|ls) - (ul|sv));   // wrong
```

where the closed-shell Fock for a density carrying the factor 2
(P = 2·C_occ·C_occᵀ) is

```cpp
jk += P(l,s) * ((uv|ls) - 0.5*(ul|sv));   // right
```

— the two-electron part of F was **doubled**. Yet H₂/STO-3G converged to the
correct energy anyway: with two identical atoms, symmetry pins the eigenvectors
to [1,1] and [1,−1] no matter how the two-electron term is scaled, so the
density — and therefore the energy, which is recomputed from the density — came
out right. Only the orbital energies betrayed it: +0.097/+1.815 Ha instead of
−0.578/+0.670. The diagnosis fell out of the arithmetic: the printed σg energy
equaled ε_core + 2·G exactly, where the correct value is ε_core + G
(ε_core = −1.2525 Ha, which the unit tests independently pin down).

Same lesson as the [Skyrme–HF solver](/code/skyrme-hf-vb)'s two bugs, which were
invisible on ¹⁶O: **the most symmetric test case is the weakest test**.
H₂/STO-3G is the ¹⁶O of quantum chemistry.

## References

- A. Szabo, N. S. Ostlund, *Modern Quantum Chemistry* (Dover, 1996) — H₂/STO-3G benchmark, RHF equations.
- S. Obara, A. Saika, *J. Chem. Phys.* **84**, 3963 (1986) — integral recurrences.
- P. Pulay, *Chem. Phys. Lett.* **73**, 393 (1980) — DIIS.