Hartree-Fock for water

The input is the molecule as PySCF takes it: the atoms as text, each an element symbol and its position in angstrom, and the basis set by its name, "sto-3g", found in a table generated from the Basis Set Exchange (lib/basis/sto3g). molecule (lib/hf/mole) turns them into the three tables lib/hf computes with: the atoms (charge and position in bohr), the shells of the basis (which atom, which angular momentum, where its primitives start) and the primitives (exponent and contraction coefficient). For water in STO-3G that is 3 atoms, 5 shells and 15 primitives, giving 7 basis functions, the same functions PySCF builds for the molecule.

Everything after that is computed by lib/hf, the same code whose tests compare it against PySCF on a set of molecules and basis sets (the reference values are in lib/hfref):

1. The overlap $S$, kinetic energy $T$ and nuclear attraction $V$ matrices. lib/hf has more than one implementation of each (three for $S$ and $T$, two for $V$). Here every implementation runs, and the program stops if any two disagree by more than $10^{-12}$.
2. The electron repulsion integrals $(\mu \nu |\lambda \sigma )$, 406 of them once the eightfold permutation symmetry is used. Two independent methods compute them: McMurchie-Davidson, which expands each product of Gaussians in Hermite Gaussians, and Head-Gordon-Pople, which runs the Obara-Saika vertical recurrence and then the horizontal one. They share no recurrence. The program compares them entry by entry and stops if any integral differs by more than $10^{-12}$ relative, naming the method and the integral.
3. The self-consistent field, using only the default implementation of each step: a starting guess from a spin-averaged calculation, DIIS acceleration (which combines the Fock matrices of the last few iterations to make the commutator $FDS-SDF$ as small as possible), canonical orthogonalization, and the Coulomb and exchange matrices J and K from the packed integrals. It solves $F(D)C=SC\epsilon$ repeatedly until the change in energy and the commutator both fall below their tolerances, and $E$ is the formula above.

The energy matches PySCF's $-74.9630631297292$ within $10^{-9}$.

The figure moves the two hydrogens along their O-H bonds, keeping the angle, and reruns the whole calculation at five more bond lengths, using only the default implementations. Each energy is checked against PySCF at the same geometry. A parabola through the three points around the stored geometry puts the bottom at 0.987 angstrom. PySCF's own minimization along the same line gives 0.98704. The stored geometry has O-H at 0.958 angstrom, close to the experimental bond length, so STO-3G makes the bond about 3 percent too long.

The density, iteration by iteration

The figure below starts again from the core Hamiltonian, the density of the electrons with no repulsion between them, and takes plain Roothaan steps: build $F$ from the density, solve $FC=SC\epsilon$, fill the five lowest orbitals, repeat, with no DIIS. Each frame is the electron density $\rho (\mathbf{r})={\sum }_{\mu \nu }{D}_{\mu \nu }{\varphi }_{\mu }(\mathbf{r}){\varphi }_{\nu }(\mathbf{r})$ in the plane of the molecule, drawn as ${\log }_{10}\rho$ in bohr${}^{-3}$, and the second figure is the energy of each frame's density. The steps end on the same energy as the accelerated run above.