Computational Quantum Chemistry from Scratch

Quantum Chemistry

The plan: by the end of this series you will have a Python program, written from scratch, that computes the binding-energy curve of H2 using the Hartree-Fock method. Every page in the series extends a single growing file. Every page the program runs. Every page it does a little more.

Most computational tutorials show you the finished artifact and annotate it. The structure here is borrowed from Salvatore Sanfilippo's kilo text-editor tutorial: incremental, runnable, narrated. At every step you see what changed, why, and what the program now does — not just what it eventually became. The intermediate scaffolding is the point. The final code is the souvenir.

What you should know going in: undergraduate quantum mechanics (the time-independent Schrödinger equation, the hydrogen atom, expectation values), and linear algebra (matrices, eigenvalue problems, dot products). What you don't need to know: anything about Hartree-Fock specifically. We'll build it.

The eleven steps

  1. The setup — Why we need approximation methods at all. A 10-line Python file that computes the energy of a 1s Slater-orbital trial wavefunction for the H atom.
  2. The variational principle — Any trial wavefunction gives an upper bound on the ground-state energy. Sweep over the orbital exponent and find the minimum. The program now optimizes.
  3. Atomic basis functions — Slater orbitals are the right shape but the wrong integrals. Gaussians integrate analytically, so chemists use linear combinations of Gaussians (STO-nG). The program now uses a contracted Gaussian basis.
  4. Many-electron wavefunctions — A Hartree product is wrong because it isn't antisymmetric. The Slater determinant fixes this. The program now describes a 2-electron system.
  5. The Hartree-Fock equations — Apply the variational principle to a single Slater determinant and the result is one coupled single-particle equation per orbital.
  6. Roothaan-Hall — Project the HF equations into a finite basis and they become a generalized matrix eigenvalue problem. The program now has a Fock matrix.
  7. The SCF algorithm — The Fock matrix depends on the orbitals it's trying to find. Fixed-point iteration: guess, build, diagonalize, repeat until convergence.
  8. One-electron integrals — Overlap, kinetic energy, and nuclear attraction. Closed-form for Gaussian basis functions; we implement them.
  9. Two-electron integrals — The four-centre repulsion integral. The expensive part of any HF calculation. We compute the H2 case directly.
  10. H2 from scratch — Putting it all together. The program now reproduces the textbook H2 binding-energy curve.
  11. What we left out — Correlation energy, post-HF methods, DFT, larger molecules. A signpost for what comes next.

How to read each page

Every page has the same shape: a short prose section explaining why the next change is needed, an inline excerpt showing the new code, the full updated file at the bottom, and what the program prints when you run it. Copy the bottom file into a single Python file (qc.py) and run it after every chapter. You should see the output match what's printed on the page.

Code dependencies: NumPy, SciPy. Nothing else. No quantum-chemistry packages — that would defeat the purpose. We compute the integrals by hand from closed-form expressions and assemble the rest.