Build Your Own Nuclear DFT

This is a build-it tutorial in the style of the famous kilo text editor walkthrough: seven chapters, each a complete program that runs, each adding one idea, ending at the production Vautherin–Brink solver — the one that reproduces the 1972 binding energies. With one upgrade on kilo: every chapter runs in your browser. Each page shows the full program, and its Open in IDE ▸ link loads that step into the Knot IDE; press Run, watch the figures.

The theory behind every formula here is on the main page. This tutorial is the other half of understanding: watching the machine get assembled, one working piece at a time.

1. Step 1 — The grid and a guess density A radial mesh, a volume integral, and a Woods–Saxon profile holding exactly 16 nucleons.
2. Step 2 — One bound state in a fixed well The radial Schrödinger equation as a tridiagonal matrix, and an O(N) route to its ground state — with the tridiagonal and Gershgorin machinery actually explained.
3. Step 3 — Shell structure appears Centrifugal barrier + spin-orbit splitting + deflation: the magic numbers assemble themselves.
4. Step 4 — From orbitals back to densities ρ, τ, and J from the occupied orbitals — and the spin-saturation cancellation, in actual arithmetic.
5. Step 5 — The force becomes fields Differentiate the Skyrme energy density and a −94 MeV well emerges that nobody drew.
6. Step 6 — Close the loop Fields → orbitals → densities → mix → repeat. Plus the half-point trick that keeps the matrix Hermitian.
7. Step 7 — The real energy — oxygen binds The functional, Coulomb, and the species split: E = −113.353 MeV, matching the Python reference to every digit.

After step 7 — the chain to ²⁰⁸Pb

Past ¹⁶O, not one term is added to the functional — the recipe from step 5 that turns densities into an energy — so the equations are finished. What heavier nuclei change is which solutions those equations must deliver, and every upgrade in skyrme_vb_chain.knot is a physical fact wearing a code disguise. The shell table (every doubly-magic filling — both proton and neutron counts sitting at closed-shell magic numbers, the main page's five nuclei — a prefix of one master list, so ⁴⁸Ca's N = 28 falls out for free) is the magic-number ordering — the spin-orbit physics of step 3, tabulated. Deflation exists because the Pauli principle marches heavier nuclei up the radial ladder into 2s, 3s, 3p states, and one lowest eigenvector no longer covers the occupation — the list of which orbitals are actually filled. And the warm-started shifts (each state's shift seeded from its own eigenvalue at the previous SCF iteration) are there because a heavy nucleus crowds its least-bound levels within an MeV of each other — the same high level density near the Fermi surface that the effective-mass story on the main page quantifies. New solver demands, no new laws.

Edit nucleus = 5.0 and the same program produces ²⁰⁸Pb at −1616.55 MeV — every rung matching the Python reference. ⁴⁰Ca runs in the browser in about four seconds; lead is the download-and-gcc case, a quarter second native.

Two bugs surfaced on the way to lead, both worth knowing about — the main page keeps the same kind of record. First: lead's 22 neutron orbitals silently overflowed a size-20 table, and the energies came out garbage — fixed-size arrays don't warn. Second: a "smarter" Rayleigh-updated shift captured wrong states mid-flight in Pb's crowded spectrum, and the SCF dutifully converged to a wrong configuration at −1005 MeV. Fixed shifts are provably safe; the reasoning is now a comment in the eigensolver.

Start with Step 1 ▸