The Atom Machine

CHAPTER THREE The master equation · paper eq. (1)

The sentence that runs the universe

There is one equation that decides which wave functions nature permits. Erwin Schrödinger wrote it in 1926, and the nablaDFT paper puts it on its very first page. In its time-independent form — the version for a molecule sitting still in its lowest state — it is just three symbols:

Equation 1 · the Schrödinger equation $$ \hat{H}\,\Psi = E\,\Psi $$

You already met this shape at the end of Chapter 2: an operator acts on a wave and gives back the same wave times a number. Let's meet all three pieces, because once this single line is yours, you own the spine of the whole paper.

Ψ — the wave function

The blurry map of all the electrons, from Chapter 2. It's the unknown we're solving for — like the $x$ in a school equation, except instead of one number it's an entire function spread across space. For a molecule with many electrons, $\Psi$ depends on the positions of all of them at once. Hold that thought; it's where the trouble starts.

E — the energy

A plain number: the total energy of the molecule in that state. Energy is the universe's accounting system. Like a ball on a hillside, a molecule always rolls toward the lowest energy it can reach — and that lowest-energy state, the ground state, is almost always the one we want, because it's the one molecules actually sit in at rest.

Ĥ — the Hamiltonian, built from scratch

The star. The hat means "operator," and $\hat H$ is the special operator that measures total energy. Rather than hand it to you as a black box, let's build it, because it's made of exactly two ideas — and you already know both.

Piece one: motion energy (kinetic). A moving thing has energy of motion. From Chapter 2, momentum is "how fast the wave wiggles," via the operator with a derivative in it. Kinetic energy in physics is $\tfrac{p^2}{2m}$ (momentum squared over twice the mass). Squaring the momentum operator means applying that derivative twice — a second derivative. In 3-D, the "second derivative in every direction" has a symbol of its own, the one that names this whole paper:

The kinetic-energy operator $$ \hat{T} \;=\; -\,\frac{\hbar^2}{2m}\,\nabla^2 $$

Decode it — and meet nabla

$\nabla^2$ — the upside-down triangle squared, read "nabla squared" or "the Laplacian." It just means "add up the second derivatives in the $x$, $y$, and $z$ directions." It measures how sharply the wave curves.
A wave that curves sharply (wiggles fast) → big $\nabla^2$ → lots of kinetic energy. Smooth, lazy wave → little. Exactly the Chapter-2 picture.
$m$ is the electron's mass, $\hbar$ is our friend h-bar. The minus sign and the constants are just there to make the energy come out right.
This is where the paper gets its name. The symbol $\nabla$ ("nabla") sits at the heart of the energy operator, so the dataset is christened nablaDFT, often written $\nabla$DFT — a quiet nod to the calculus living inside $\hat H$.

Piece two: position energy (potential). Charges that attract or repel store energy depending on how far apart they are — that's Coulomb's law from Chapter 1. We bundle all of that into a potential-energy operator $\hat V$. Put the two pieces together and you have the whole machine:

The Hamiltonian = motion + position $$ \hat{H} \;=\; \hat{T} \;+\; \hat{V} \;=\; -\frac{\hbar^2}{2m}\nabla^2 \;+\; \hat{V} $$

The real molecular Hamiltonian, written out in full

For an actual molecule — many electrons, many nuclei — $\hat V$ is not one term but a sum of every Coulomb interaction in the system. Here is the complete thing. It looks fearsome; we'll dismantle it term by term, and you'll find there are no surprises in it:

The full molecular Hamiltonian $$ \hat{H} = \underbrace{-\sum_i \tfrac{\hbar^2}{2m}\nabla_i^2}_{\text{electrons moving}} \;-\; \underbrace{\sum_{i,A}\tfrac{k\,Z_A e^2}{r_{iA}}}_{\substack{\text{electrons pulled}\\\text{to nuclei}}} \;+\; \underbrace{\sum_{i

Four terms, read left to right

Term 1 — electrons moving. One kinetic piece for each electron $i$. Lots of $\nabla^2$'s. This is the only term that raises energy through motion.
Term 2 — the good attraction. Every electron $i$ pulled toward every nucleus $A$. It's negative (energy-lowering) because opposite charges attract. $Z_A$ is the nucleus's charge, $r_{iA}$ the distance — pure Coulomb's law from Chapter 1.
Term 3 — the troublemaker. Every electron pushing every other electron, since like charges repel. The \(iThis is the term that ruins everything, and the next box explains why.
Term 4 — nuclei repelling nuclei. The positive nuclei push apart too. Once we freeze the nuclei (next section), this becomes a mere constant we can add at the end.

Why term 3 is a catastrophe Look at term 3: each electron's behaviour depends on the distance $r_{ij}$ to every other electron. So you cannot solve for one electron without already knowing where all the others are — but their positions depend on this one in turn. Every electron is tangled with every other, all at once, like a knot that tightens the moment you pull any strand. This single coupling term is why $\hat H\Psi=E\Psi$ has no pen-and-paper solution for any molecule bigger than hydrogen. The entire rest of this book — DFT, the matrices, the neural networks — is a campaign against term 3.

A first mercy: freeze the nuclei (Born–Oppenheimer)

Here's the first great simplification, and it's allowed by a simple observation. The lightest nucleus (a single proton) is about 1,836 times heavier than an electron. Heavier means slower — vastly slower. From an electron's frantic point of view, the nuclei are practically nailed in place, like boulders to a swarm of flies.

So we cheat, productively: treat the nuclei as frozen, solve for the lightning-fast electrons in that fixed scaffold, and only afterward worry about how the nuclei themselves might shift. This is the Born–Oppenheimer approximation, and it's the unspoken floor beneath essentially all of computational chemistry — including every calculation in nablaDFT.

What Born–Oppenheimer buys us Freezing the nuclei does two things. It kills the nuclei's own motion energy, and it turns term 4 (nucleus–nucleus repulsion) into a fixed number we just add at the end. What's left is "the electrons' problem": find the electron wave function for a given arrangement of nuclei. That arrangement is exactly what the paper calls a "conformation" — a frozen 3-D pose of the molecule. nablaDFT computes the electron answer for over five million such poses.

Reading the whole sentence — and what "eigen" means

Now put equation 1 back together. $\hat H\Psi = E\Psi$ says: "when the energy-machine acts on the right wave, it returns that very same wave, just scaled by a number — the energy." That's the pattern from the momentum example in Chapter 2.

This pattern is so important it has a name. $\Psi$ is an eigenfunction of $\hat H$, and $E$ is its eigenvalue. (German eigen = "own, characteristic.") Most waves, fed into $\hat H$, come out scrambled into a different shape. Only special "characteristic" waves come back unchanged-but-scaled. Those, and only those, are the wave functions nature allows — and each one comes with its own energy.

Eigen-problems are easiest to feel with a small matrix instead of an operator, so let's do one entirely by hand.

Worked example · finding eigenvalues of a 2×2 matrix by hand

Take the matrix $\mathbf{A} = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$. We want the special vectors $\mathbf v$ that satisfy $\mathbf A \mathbf v = \lambda \mathbf v$ — "matrix in, same vector out, times a number $\lambda$." Same shape as Schrödinger, just smaller.

The trick: $\mathbf A \mathbf v = \lambda \mathbf v$ rearranges to $(\mathbf A - \lambda \mathbf I)\mathbf v = 0$. For a non-zero $\mathbf v$ to exist, the matrix $(\mathbf A - \lambda \mathbf I)$ must be "squashable" — its determinant must be zero.

Write it out: $\det\!\begin{bmatrix} 2-\lambda & 1 \\ 1 & 2-\lambda \end{bmatrix} = (2-\lambda)^2 - 1 = 0.$

Expand: $(2-\lambda)^2 = 1$, so $2-\lambda = \pm 1$, giving $\lambda = 1$ or $\lambda = 3$. Those are the two eigenvalues — the only two scaling-numbers this matrix allows.

Find the matching vectors: $\lambda=3$ gives $\mathbf v = (1,1)$ — point along the diagonal and the matrix just triples you. $\lambda=1$ gives $\mathbf v = (1,-1)$ — point along the anti-diagonal and it leaves you alone.

That's the whole game. In a real molecule the matrix is enormous and the eigenvalues are the allowed energies — but the method is precisely this. nablaDFT's hardest task is producing that matrix; once you have it, finding the energies is just this, scaled up.

Build that 2×2 yourself below. Drag the matrix entries and watch its eigenvectors — the special directions that don't get knocked off course — snap to new angles.

Try it: hunt for eigenvectors

A grey vector is being transformed by the matrix into a coloured one. Spin the grey vector with the slider. When the coloured arrow lines up perfectly with the grey one (just longer or shorter, not rotated), you've found an eigenvector — and the stretch factor is its eigenvalue.

Matrix a,b / c,d

Rotate test vector 20°

Just how impossible is the real thing?

Our 2×2 had two numbers. A real molecule's wave function lives in a space so large it defies picturing. Here's the catastrophe in arithmetic.

Worked example · the curse of dimensionality

To store a wave on a computer, chop space into a grid. Say a modest $100$ points along each direction for a single electron in 3-D — that's already $100^3 = 1{,}000{,}000$ numbers. Manageable.

But $\Psi$ for $N$ electrons depends on all of their positions together. The grid sizes don't add — they multiply. For $N$ electrons you need roughly $(100^3)^N = 10^{6N}$ numbers.

For a tiny molecule with just $N=10$ electrons: $10^{60}$ numbers. That is more storage than there are atoms in the observable universe (about $10^{80}$, and we'd blow past it by $N=14$).

Conclusion: you cannot store the exact wave function of even a small molecule, let alone solve for it. The honest exact answer is forever out of reach. Every method that follows is a brilliant, principled approximation — a way to get a good-enough answer without ever writing down the impossible object.

So where does this leave us? We have the perfect equation ($\hat H\Psi=E\Psi$) and a perfect understanding of what's in it — and it's unsolvable for anything real. That tension is the engine of this entire field. The next chapter is the first great escape: stop trying to store the smooth, infinite wave, and turn it into a small grid of numbers a computer can actually hold.

CHAPTER THREE ½ The rules underneath · a rigorous interlude

The laws of energy

Before we start turning the Schrödinger equation into something a computer can chew on, we need to slow down and meet the laws that energy obeys. These are not decorations — they are the load-bearing walls. The variational principle in particular is the single idea that makes Hartree–Fock, DFT, and almost every method in both papers actually work. This chapter is more demanding than the ones around it; everything here is derived, not asserted. Take it slowly, and the rest of the book becomes inevitable rather than mysterious.

First, five small tools (so nothing later feels like a leap)

The proofs ahead lean on a handful of bits of notation and a few facts that are genuinely easy once spelled out — but baffling if dropped on you mid-derivation. So let's collect them here, slowly and with no proof of its own required. Come back to this box any time a later step looks like it skipped something; the missing rung is almost certainly here.

Tool 1 · the bracket ⟨ ⟩ is just "multiply and add up" You met the ket $|\psi\rangle$ in Chapter 4 — it's "a wave." The bracket $\langle\phi|\psi\rangle$ (a "bra" $\langle\phi|$ closing a "ket") means one specific, concrete thing: multiply the two functions together at every point and add up the result over all space — an integral, $\langle\phi|\psi\rangle = \int \phi^*(\vec r)\,\psi(\vec r)\,d\vec r$. The little star $^*$ means "flip the sign of the imaginary part" (complex conjugate); for the real functions we mostly use, it does nothing. So a bracket is a single number measuring "how much these two waves overlap." That's all. The overlap integral $S_{ij}$ from Chapter 4 was exactly a bracket, $\langle\phi_i|\phi_j\rangle$.

Tool 2 · ⟨ψ|Ĥ|ψ⟩ is the average energy, also just a number Put an operator in the middle, $\langle\psi|\hat H|\psi\rangle$, and the recipe is: first let $\hat H$ act on $|\psi\rangle$ (giving a new function), then take the bracket of $\langle\psi|$ with that result. The output is again a single number — the average energy you'd measure for a molecule in state $\psi$. We write it $\langle\hat H\rangle$ for short. When $\psi$ happens to be an eigenstate, $\hat H|\psi\rangle = E|\psi\rangle$, so $\langle\psi|\hat H|\psi\rangle = E\langle\psi|\psi\rangle = E$ — the average energy is just the energy, exactly as you'd hope.

Tool 3 · any wave is a sum of the "pure note" waves (completeness) A musical chord can be broken into pure notes. In the same way, the special eigenstates $\{\psi_n\}$ of $\hat H$ — the allowed "pure note" waves from Chapter 3 — are complete: any wave at all can be written as a weighted sum of them, $\psi = \sum_n c_n\psi_n$, with numbers $c_n$ saying how much of each note is present. This is the same "build it from building blocks" idea as the basis sets of Chapter 4, but now the blocks are the true eigenstates. We never need to know what they actually are — only that this sum is always possible. That single guarantee is the seed of the variational proof below.

Tool 4 · the pure notes don't overlap (orthonormality & the δ symbol) The pure-note waves are orthonormal: each one is normalized ($\langle\psi_n|\psi_n\rangle = 1$), and any two different ones don't overlap at all ($\langle\psi_m|\psi_n\rangle = 0$ when $m\ne n$). Both facts at once are written with a shorthand symbol, the Kronecker delta: $\langle\psi_m|\psi_n\rangle = \delta_{mn}$, where $\delta_{mn}$ simply means "1 if $m=n$, 0 otherwise." Its one job in a sum is to act like a switch: $\sum_n \delta_{mn}(\text{stuff}_n)$ keeps only the $n=m$ term and deletes all the rest — which is how a double sum collapses to a single sum in the proofs.

Tool 5 · the trace, and one identity it obeys The trace of a square matrix, $\operatorname{Tr}(\mathbf M)$, is the plain sum of its diagonal entries — nothing fancier. We'll need exactly one fact about it: you can cyclically rotate a product inside a trace without changing the answer, $\operatorname{Tr}(\mathbf A\mathbf B) = \operatorname{Tr}(\mathbf B\mathbf A)$. (It follows from writing out the diagonal sums, but you can take it on faith.) This is the one lever used much later, in Chapter 6½, to show the density matrix counts the electrons.

That's the whole toolkit. None of it is hard; it's just vocabulary plus two switch-like facts (the δ symbol and the trace identity). With these in hand, every derivation in this chapter and the rigorous ones later will read as a sequence of small, checkable steps — no leaps.

First, a cleaner system of units

Physicists working on atoms got tired of writing $\hbar$, the electron mass $m_e$, the electron charge $e$, and the constant $4\pi\varepsilon_0$ in every single equation. So they simply defined all four to equal 1. This is the system of atomic units (a.u.), and it is what both papers — and all of quantum chemistry — actually use.

Atomic units, and why they're natural

Set $\hbar = m_e = e = 4\pi\varepsilon_0 = 1$. Every formula loses its clutter.
The unit of length becomes the Bohr radius $a_0 \approx 0.529$ Å — about the size of a hydrogen atom. Lengths are now measured in "how many hydrogen atoms across."
The unit of energy becomes the Hartree, $E_{\text h} \approx 27.2$ eV $\approx 627.5$ kcal/mol. It is, by construction, twice the binding energy of a hydrogen atom.
The kinetic operator $-\frac{\hbar^2}{2m}\nabla^2$ collapses to a clean $-\frac{1}{2}\nabla^2$, and Coulomb's law $\frac{1}{4\pi\varepsilon_0}\frac{q_1q_2}{r}$ becomes just $\frac{q_1q_2}{r}$. We'll use these tidy forms from here on.

Keep the conversion $1\,E_{\text h} \approx 627.5$ kcal/mol in mind: the "chemical accuracy" target from Chapter 8 (1 kcal/mol) is therefore about $1/627.5 \approx 1.6 \times 10^{-3}$ Hartree. Every error bar in both papers is quoted in these units.

Law 1 — energy is conserved (and why)

So far we've used the time-independent Schrödinger equation, $\hat H\psi = E\psi$. It is a snapshot. The full law includes time, and reads:

The time-dependent Schrödinger equation $$ i\hbar\,\frac{\partial \Psi}{\partial t} \;=\; \hat{H}\,\Psi $$

Let's prove that energy doesn't drift over time. The expected energy of any state is $\langle \hat H\rangle = \langle\Psi|\hat H|\Psi\rangle$. We show its time derivative is zero.

Derivation · conservation of energy

Differentiate using the product rule: $\dfrac{d\langle\hat H\rangle}{dt} = \left\langle\dfrac{\partial\Psi}{\partial t}\Big|\hat H\Big|\Psi\right\rangle + \left\langle\Psi\Big|\hat H\Big|\dfrac{\partial\Psi}{\partial t}\right\rangle$ (the operator $\hat H$ itself has no explicit time dependence for an isolated molecule).

From the time-dependent equation, $\dfrac{\partial\Psi}{\partial t} = \dfrac{1}{i\hbar}\hat H\Psi$. Substitute it into both slots.

The first term becomes $\dfrac{-1}{i\hbar}\langle\hat H\Psi|\hat H|\Psi\rangle$, and the second becomes $\dfrac{1}{i\hbar}\langle\Psi|\hat H\,\hat H|\Psi\rangle$. (Why the minus sign on the first? When a factor moves out of the left slot of a bracket it gets conjugated — Tool 1 — and the conjugate of $\tfrac{1}{i\hbar}$ is $\tfrac{-1}{i\hbar}$, since conjugating flips the sign of $i$.)

Because $\hat H$ is self-adjoint — a Hermitian operator, meaning you can move it from one slot of a bracket to the other without changing the value, $\langle\hat H\phi|\psi\rangle = \langle\phi|\hat H\psi\rangle$ — we have $\langle\hat H\Psi|\hat H\Psi\rangle = \langle\Psi|\hat H\hat H\Psi\rangle$. The two terms are now identical except for their opposite signs.

They cancel: $\dfrac{d\langle\hat H\rangle}{dt} = 0$.

energy is constant in time

The deeper reason — Noether's idea That cancellation is not a coincidence; it's a shadow of one of the most beautiful theorems in physics. Emmy Noether proved in 1918 that every continuous symmetry of a physical system corresponds to a conserved quantity. Energy conservation is precisely the consequence of time-translation symmetry: the laws governing the molecule are the same today as tomorrow, and that sameness is energy conservation. (Momentum conservation comes from space-translation symmetry; angular-momentum conservation from rotational symmetry — the same rotational symmetry that will dominate Chapter 6½.)

Law 2 — the variational principle (the engine of everything)

Here is the most important single statement in this entire book. It is what lets us find approximate answers without ever solving the impossible equation exactly. Take any normalized trial wave function $\psi$ you like — a guess, however crude. Compute its expected energy. The claim is that this number can never dip below the true ground-state energy:

The variational principle

For any normalizable trial state $\psi$, the energy expectation value is an upper bound on the true ground-state energy $E_0$: $$ E[\psi] \;=\; \frac{\langle\psi|\hat H|\psi\rangle}{\langle\psi|\psi\rangle} \;\ge\; E_0. $$ Equality holds if and only if $\psi$ is the true ground state.

The proof is short and it is the reason the principle is trustworthy. The trick is to expand the unknown trial state in the (also unknown, but guaranteed to exist) complete set of true eigenstates.

Derivation · why no guess can beat the ground state

By Tool 3, the true eigenstates $\{\psi_n\}$ form a complete orthonormal set, with $\hat H\psi_n = E_n\psi_n$ and energies ordered $E_0 \le E_1 \le E_2 \le \cdots$. So any trial state, including ours, expands as $\psi = \sum_n c_n \psi_n$ — even though we'll never need to know what the $\psi_n$ actually are.

Normalization plus orthonormality (Tool 4) gives $\langle\psi|\psi\rangle = \sum_n |c_n|^2 = 1$. (Reason: $\langle\psi|\psi\rangle = \sum_{m,n} c_m^* c_n\langle\psi_m|\psi_n\rangle = \sum_{m,n} c_m^* c_n\,\delta_{mn} = \sum_n |c_n|^2$, the $\delta$ switch keeping only $m=n$.)

Now the energy. Substitute the expansion and use that $\hat H$ acts on each note separately ($\hat H\psi_n = E_n\psi_n$): $\langle\psi|\hat H|\psi\rangle = \sum_{m,n} c_m^* c_n \langle\psi_m|\hat H|\psi_n\rangle = \sum_{m,n} c_m^* c_n E_n \langle\psi_m|\psi_n\rangle$.

The same $\delta$ switch (Tool 4) collapses the double sum to a single one: $\langle\psi|\hat H|\psi\rangle = \sum_n |c_n|^2 E_n$.

Now the key step. Every $E_n \ge E_0$, so replacing each $E_n$ by the smallest, $E_0$, can only decrease the sum: $\sum_n |c_n|^2 E_n \ge \sum_n |c_n|^2 E_0 = E_0 \sum_n |c_n|^2 = E_0$.

Therefore $\langle\psi|\hat H|\psi\rangle \ge E_0$. Equality requires all the weight on $c_0$, i.e. $\psi = \psi_0$.

any trial energy sits above the truth

Why this changes everything The variational principle converts an unsolvable differential equation into a minimization problem. Instead of solving $\hat H\psi = E\psi$ head-on, we write down a trial wave function with adjustable knobs, compute its energy, and turn the knobs to push the energy as low as it will go. Whatever we reach is a rigorous upper bound on the truth — and the lower we get, the better our wave function must be. Hartree–Fock does exactly this. The Roothaan equation of Chapter 4 is what you get when you apply this principle to a basis-set expansion. And the second Hohenberg–Kohn theorem (Chapter 5) is this same principle, rewritten for the density. Hold this thought; it is the thread connecting the whole book.

Worked example · a variational estimate by hand

Consider a particle in the potential $V(x) = \frac{1}{2}x^2$ (a harmonic well, in atomic units with the mass set to 1). We'll guess the ground state without solving anything.

Take a Gaussian trial function with one adjustable width parameter $\alpha$: $\psi_\alpha(x) = e^{-\alpha x^2/2}$. Wider or narrower is set by $\alpha$.

Now compute its average energy, $E(\alpha) = \langle\hat T\rangle + \langle\hat V\rangle$, one piece at a time. We need two standard Gaussian integral facts (the same bell-shaped integrals from Chapter 4): for $\psi_\alpha = e^{-\alpha x^2/2}$, the normalization integral is $\langle\psi|\psi\rangle = \int e^{-\alpha x^2}dx = \sqrt{\pi/\alpha}$, and $\int x^2 e^{-\alpha x^2}dx = \tfrac{1}{2\alpha}\sqrt{\pi/\alpha}$.

Potential piece. $\langle\hat V\rangle = \dfrac{\int \tfrac12 x^2\,e^{-\alpha x^2}dx}{\int e^{-\alpha x^2}dx} = \tfrac12\cdot\dfrac{\tfrac{1}{2\alpha}\sqrt{\pi/\alpha}}{\sqrt{\pi/\alpha}} = \dfrac{1}{4\alpha}$ — the $\sqrt{\pi/\alpha}$'s cancel, leaving just $1/4\alpha$.

Kinetic piece. Applying $-\tfrac12\tfrac{d^2}{dx^2}$ to $e^{-\alpha x^2/2}$ and taking the same ratio gives $\langle\hat T\rangle = \dfrac{\alpha}{4}$ (one factor of $\alpha$ appears because each derivative pulls one down, exactly as in the virial scaling argument).

Add them: $E(\alpha) = \dfrac{\alpha}{4} + \dfrac{1}{4\alpha}$. A tug-of-war — the kinetic term wants $\alpha$ large (a narrow, sharply-curved wave), the potential term wants $\alpha$ small (a wide wave sitting low in the well).

Minimize: $\dfrac{dE}{d\alpha} = \dfrac14 - \dfrac{1}{4\alpha^2} = 0$, giving $\alpha = 1$.

Then $E(1) = \frac14 + \frac14 = \frac12$. The exact ground-state energy of this well is also $\frac12$ — our guess happened to contain the true answer, so the bound is tight. In a real molecule it won't be exact, but it is always an honest ceiling. That guarantee is what we lean on.

Law 3 — forces come for free (Hellmann–Feynman)

Chapter 7 claimed that once you know the energy, the forces follow by differentiating, $\vec F = -\nabla E$. That is true, but there's a subtlety worth making rigorous: the wave function itself changes when you move a nucleus, so naively you'd expect extra terms. The Hellmann–Feynman theorem says those extra terms vanish.

The Hellmann–Feynman theorem

If $\psi_\lambda$ is a normalized eigenstate of a Hamiltonian $\hat H_\lambda$ depending on a parameter $\lambda$, then the energy's derivative only sees the operator's explicit dependence: $$ \frac{dE}{d\lambda} \;=\; \left\langle \psi_\lambda \left| \frac{\partial \hat H_\lambda}{\partial \lambda} \right| \psi_\lambda \right\rangle. $$

Derivation · the wave-function terms cancel

Start from $E = \langle\psi_\lambda|\hat H_\lambda|\psi_\lambda\rangle$ with $\langle\psi_\lambda|\psi_\lambda\rangle = 1$. This is a product of three $\lambda$-dependent things — the bra, the operator, and the ket — so the product rule gives three terms, one for differentiating each in turn: $\dfrac{dE}{d\lambda} = \big\langle\partial_\lambda\psi|\hat H|\psi\big\rangle + \big\langle\psi|\partial_\lambda\hat H|\psi\big\rangle + \big\langle\psi|\hat H|\partial_\lambda\psi\big\rangle$, writing $\partial_\lambda$ for "$\tfrac{d}{d\lambda}$ of."

But that bracket is the derivative of the normalization: $\dfrac{d}{d\lambda}\langle\psi|\psi\rangle = \langle\partial_\lambda\psi|\psi\rangle + \langle\psi|\partial_\lambda\psi\rangle = \dfrac{d}{d\lambda}(1) = 0$.

So the first and third pieces vanish together, leaving only the middle one.

only the operator's own change survives

Now specialize $\lambda$ to the position of a nucleus. The only part of $\hat H$ that depends on a nuclear coordinate is the electron–nucleus attraction and nucleus–nucleus repulsion — both ordinary Coulomb terms. So the force on a nucleus is just the classical electrostatic force exerted on it by the electron cloud and the other nuclei. This is the rigorous foundation under "forces are the gradient of the energy" — the principle that drives the relaxation trajectories of Chapter 9. A neural network that predicts energy well, and is smooth, automatically predicts these forces.

Law 4 — kinetic and potential are locked together (the virial theorem)

One last law, and a strikingly specific one. For any system bound purely by Coulomb forces — which is to say, every atom and molecule — the average kinetic and potential energies are not independent. They are chained in a fixed ratio.

The virial theorem (for Coulomb systems)

At any true stationary state of a molecule, $$ \langle \hat T\rangle = -\tfrac{1}{2}\langle \hat V\rangle, \qquad\text{equivalently}\qquad E = -\langle\hat T\rangle = \tfrac{1}{2}\langle\hat V\rangle. $$

We can derive it cleanly from the variational principle itself, using a scaling trick — which also shows how deeply the variational idea runs.

Derivation · the scaling argument

Take the true ground state $\psi(\vec r)$ and squeeze or stretch it by a scale factor $\eta$: define $\psi_\eta(\vec r) = \eta^{3N/2}\psi(\eta\vec r)$ (the prefactor keeps it normalized for $N$ electrons in 3-D).

Now see how each energy piece responds to the squeeze. Kinetic energy carries two derivatives ($-\tfrac12\nabla^2$), and each derivative with respect to $\vec r$ brings down one factor of $\eta$ when $\vec r\to\eta\vec r$ — so kinetic energy picks up $\eta\times\eta=\eta^2$: $\langle\hat T\rangle_\eta = \eta^2\langle\hat T\rangle$. Coulomb potential energy goes as $1/r$, and $1/(\eta r)$ is just $\eta$ times $1/r$: $\langle\hat V\rangle_\eta = \eta\langle\hat V\rangle$. (Both statements use that the normalization prefactor in step 1 keeps the integrals honest; the upshot is the clean powers $\eta^2$ and $\eta^1$.)

The scaled energy is $E(\eta) = \eta^2\langle\hat T\rangle + \eta\langle\hat V\rangle$. Since the true ground state is the variational minimum, $\eta = 1$ must minimize $E(\eta)$.

Set the derivative to zero at $\eta = 1$: $\dfrac{dE}{d\eta}\Big|_{\eta=1} = 2\langle\hat T\rangle + \langle\hat V\rangle = 0$.

Rearrange: $\langle\hat T\rangle = -\tfrac12\langle\hat V\rangle$. And since $E = \langle\hat T\rangle + \langle\hat V\rangle$, substituting gives $E = -\langle\hat T\rangle$.

kinetic and potential are fixed in a 1:−2 ratio

A sanity check on hydrogen, and why this is useful The hydrogen atom's ground-state energy is $E = -\tfrac12$ Hartree. The virial theorem then predicts, with no further calculation, that its average kinetic energy is $\langle\hat T\rangle = +\tfrac12$ and its average potential energy is $\langle\hat V\rangle = -1$ Hartree — and direct calculation confirms both. Beyond elegance, the virial theorem is a practical lie-detector: when a DFT calculation finishes, codes often check the virial ratio $-\langle\hat V\rangle / \langle\hat T\rangle$ against 2. A large deviation is a red flag that the basis set was too small or the calculation didn't converge. It is one of the quiet quality controls behind the "ground truth" energies in both datasets.

Four laws: energy is conserved, no trial state beats the ground state, forces are electrostatic, and kinetic and potential energies are chained together. With these in hand, the methods of the next chapters stop looking like clever hacks and start looking like the only reasonable things to do. Onward — to turning the wave into numbers, now knowing exactly why minimizing energy is the right move.

CHAPTER FOUR The great trick · paper eqs. (2)–(7)

Turning a wave into a grid of numbers

We ended Chapter 3 stuck: the perfect equation is unsolvable, and we can't even store its answer. This chapter is the escape, and it rests on one of the most powerful ideas in all of mathematics. You half-know it already — it's how mixing paint works.

Why you can build any function from a fixed set of pieces

Think about ordinary arrows (vectors) on a flat page. Pick two reference arrows — one pointing right, one pointing up. Any arrow on the page can be built as "so much right plus so much up." The arrow pointing northeast at some length is just, say, $3$ rights $+ 2$ ups. Two numbers, $(3, 2)$, capture the whole arrow. You've turned a geometric object into a short list of numbers.

Here is the leap: functions behave just like arrows. A wave function is a kind of "arrow" in an enormous space where each "direction" is a simple, standard shape. Pick a fixed set of those standard shapes — the paper calls them basis functions, written $\lvert\phi_i\rangle$ — and any wave function can be built as a weighted sum of them, exactly like the arrow. The weights are a short list of numbers. That list is something a computer can hold.

The same idea, twice. Left: a geometric arrow becomes the two numbers $(3,2)$. Right: a complicated wave (purple) becomes a list of "how much of each basis shape." Both swap an awkward continuous object for a tidy list. This is what makes chemistry computable.

The paper writes this build-from-pieces idea as equation 2 — a single electron's orbital $\lvert\psi_m\rangle$ assembled from basis pieces $\lvert\phi_i\rangle$:

Equation 2 · mixing building-block waves $$ \lvert \psi_m \rangle \;=\; \sum_{i} c^{\,i}_{m}\,\lvert \phi_i \rangle $$

Decode it slowly

$\lvert \psi_m \rangle$ — the real orbital we want (the $m$-th one). The angled brackets $\lvert\;\rangle$, called a "ket," are just physicist packaging meaning "a quantum wave." Read "psi-m."
$\lvert \phi_i \rangle$ — basis piece number $i$: a fixed, known shape from our box of standard pieces.
$c^{\,i}_{m}$ — the mixing amounts (the "3" and "2" from the arrow). These coefficients are the entire prize. Find them and you've found the orbital.
$\sum_i$ — "add up over all the pieces." The big Greek S (sigma) means sum.

What the basis pieces actually are: Gaussians

In real calculations — and in nablaDFT — the basis pieces are Gaussian functions: little bell-shaped bumps, $e^{-\alpha r^2}$, each centred on an atom. Why bells? Because there's a mathematical miracle: when you multiply two Gaussians, or integrate them, you get another Gaussian. Those integrals (which we're about to need by the millions) become quick and exact instead of nightmarish. Chemistry runs on Gaussians for this one practical reason.

Decoding "def2-SVP" — the basis the paper uses nablaDFT computes everything in a specific basis named def2-SVP. It looks like a licence plate; it's actually a recipe. def2 is the family (designed in Karlsruhe, Germany). SV = "split valence" — use a couple of Gaussians for each important outer electron so the orbital can flex in size. P = "polarization" — toss in a few extra higher-shape (p- and d-like) bells so orbitals can lean and distort when atoms bond. Bigger basis = more pieces = more accurate and more expensive. def2-SVP is a sensible middle: good enough to trust, small enough to run a million times.

Now derive the matrix equation — actually derive it

Here's the payoff, and we're going to earn it rather than be handed it. Start from Schrödinger and substitute our basis expansion. Watch the smooth equation collapse into a grid of numbers in four moves.

Derivation · from Schrödinger to a matrix

Move 1. Start with $\hat H \psi = E\psi$, and write $\psi$ as its expansion: $\hat H \sum_j c_j \lvert\phi_j\rangle = E \sum_j c_j \lvert\phi_j\rangle.$

Move 2. The operator $\hat H$ acts on each fixed piece, and the constants $c_j$ slide out: $\sum_j c_j\, \hat H\lvert\phi_j\rangle = E \sum_j c_j \lvert\phi_j\rangle.$ This is still one equation about whole functions — not yet numbers.

Move 3. To extract numbers, "test" the equation against each basis piece $\langle\phi_i\rvert$ in turn (this means multiply by $\phi_i$ and integrate — a way of asking "how much does each side overlap with piece $i$?"): $\sum_j c_j \langle\phi_i\rvert \hat H\lvert\phi_j\rangle = E \sum_j c_j \langle\phi_i\vert\phi_j\rangle.$

Move 4. Now name the two kinds of integral that appeared. Call $\langle\phi_i\rvert \hat H\lvert\phi_j\rangle = H_{ij}$ and $\langle\phi_i\vert\phi_j\rangle = S_{ij}$. Each is just a number — one per pair $(i,j)$. The equation becomes $\sum_j H_{ij}\,c_j = E \sum_j S_{ij}\,c_j$, which is exactly a matrix equation. We're done — the smooth wave equation is now arithmetic on grids.

That last line, packed into matrix shorthand, is precisely how the paper writes equation 3 — known as the Roothaan equation:

Equation 3 · the Roothaan equation (Schrödinger, as a matrix) $$ \mathbf{F}\,\mathbf{c}_\sigma \;=\; \varepsilon_\sigma\,\mathbf{S}\,\mathbf{c}_\sigma $$

Same eigen-idea, matrix costume

$\hat H$ the operator → $\mathbf{F}$, the Fock matrix (a grid; also called the Hamiltonian matrix $\mathbf H$). The numbers $H_{ij}$ from the derivation are its entries.
$\Psi$ the wave → $\mathbf{c}_\sigma$, the plain list of mixing numbers (a vector). The $\sigma$ is the spin label you met in Chapter 2.
$E$ the energy → $\varepsilon_\sigma$ (small epsilon), this orbital's energy.
$\mathbf{S}$, the overlap matrix, is the genuinely new character — the grid of $S_{ij}$'s. The next section is entirely about why it has to be there.

And the two grids are filled by exactly the integrals our derivation produced (the paper's equations 4 and 5):

Equations 4 & 5 · how to fill the two grids $$ H_{ij} = \langle \phi_i \lvert \hat{H} \rvert \phi_j \rangle \qquad\qquad S_{ij} = \langle \phi_i \vert \phi_j \rangle $$

Where this equation really comes from: minimizing energy

We derived the Roothaan equation by substituting the basis expansion and "testing" against each basis function. That works, but it hides why this is the right equation. The honest origin is the variational principle from Chapter 3½: the best coefficients are the ones that make the energy as low as possible. Let's show that minimizing the energy forces exactly this matrix equation into existence.

Derivation · the variational principle gives Roothaan–Hall

Write a single orbital as $\psi = \sum_i c_i \phi_i$. Its energy (Chapter 3½'s variational ratio) is $E = \dfrac{\langle\psi|\hat H|\psi\rangle}{\langle\psi|\psi\rangle} = \dfrac{\sum_{ij} c_i^* c_j H_{ij}}{\sum_{ij} c_i^* c_j S_{ij}}$, using exactly the $H_{ij}$ and $S_{ij}$ integrals above.

The variational principle says the true orbital minimizes $E$. So set the derivative with respect to each coefficient to zero: $\dfrac{\partial E}{\partial c_k^*} = 0$ for every $k$.

Differentiate the quotient. Write the top as $P = \sum_{ij}c_i^* c_j H_{ij}$ and the bottom as $N = \sum_{ij}c_i^* c_j S_{ij}$, so $E = P/N$. The quotient rule gives $\dfrac{\partial E}{\partial c_k^*} = \dfrac{1}{N}\dfrac{\partial P}{\partial c_k^*} - \dfrac{P}{N^2}\dfrac{\partial N}{\partial c_k^*}$. Now $\partial P/\partial c_k^* = \sum_j H_{kj}c_j$ and $\partial N/\partial c_k^* = \sum_j S_{kj}c_j$ (only the $c_k^*$ term in each sum survives the derivative). Substituting, and using $P/N = E$, this is $\dfrac{\partial E}{\partial c_k^*} = \dfrac{1}{N}\left[\sum_j H_{kj}c_j - E\sum_j S_{kj}c_j\right]$. Set it to zero.

Since $N \ne 0$, the bracket must vanish for every $k$: $\sum_j H_{kj}c_j = E\sum_j S_{kj}c_j$.

That is the Roothaan equation, row by row: $\mathbf{H}\mathbf{c} = E\,\mathbf{S}\mathbf{c}$. (With electron–electron interaction folded in, $\mathbf H$ becomes the Fock matrix $\mathbf F$, giving $\mathbf F\mathbf c = \varepsilon\mathbf S\mathbf c$ exactly.)

minimizing energy is the matrix eigenproblem

The same idea, every time Notice what just happened: "solve the Schrödinger equation in a basis" and "minimize the energy over the coefficients" turned out to be the same equation. This is the variational principle doing its job. Every method that follows — Hartree–Fock, and the Kohn–Sham equations of DFT in the next chapter — is built by writing down an energy and minimizing it. If you understand this derivation, you understand the skeleton of all of quantum chemistry.

Why the overlap matrix S exists (and isn't just the identity)

With the perfect arrow basis — right and up — the two reference arrows were at a clean right angle, totally independent. Overlap between them was zero, and the "overlap grid" would just be the boring identity (1's on the diagonal, 0's elsewhere). But our Gaussian bells, parked on neighbouring atoms, do spill into each other's space. They are not independent. $\mathbf S$ is the bookkeeping that records exactly how much each pair leans on its neighbour. Let's compute one entry.

Worked example · two overlapping bells

Two basis bells, $\phi_1$ on atom A and $\phi_2$ on atom B. The overlap is $S_{12} = \int \phi_1(\vec r)\,\phi_2(\vec r)\,dV$ — "how much do they share the same space?"

If the atoms are far apart: where $\phi_1$ is big, $\phi_2$ is essentially zero, so their product is zero everywhere. $S_{12} \approx 0$. The bells ignore each other.

If the atoms sit right on top of each other: the bells coincide, the product is just $\phi^2$, and (since each bell is normalized) the integral is $1$. $S_{12} = 1$.

In a real bond, somewhere between: a partial overlap like $S_{12} = 0.4$. That number literally measures how strongly the two atoms' clouds mingle — which is to say, how strong the bond is.

So $\mathbf S$ is never the boring identity for real molecules. It's a live map of which atoms are touching. Forgetting it would be solving the wrong equation.

Why all this matters so much A computer cannot solve a smooth wave equation. But finding the eigenvalues and eigenvectors of a grid of numbers is something it does millions of times a second. By turning $\hat H\Psi=E\Psi$ into $\mathbf F\mathbf c = \varepsilon\mathbf S\mathbf c$, we converted the impossible into the merely expensive. This translation is the bedrock under every single method in the paper.

The catch nobody warns you about: F depends on its own answer

There's a snake in the grass. Remember term 3 from Chapter 3 — electrons pushing each other? To build the Fock matrix $\mathbf F$, you need to know the electron cloud's shape, because each electron feels the average push of all the others. But the cloud's shape is the answer — it comes from the very coefficients $\mathbf c$ you're trying to find. You need the answer to build the equation whose answer you need. A perfect chicken-and-egg.

The escape is wonderfully simple: guess, then improve, then repeat until nothing changes. This loop is called the Self-Consistent Field (SCF) method, and it's the actual engine cranking inside every DFT calculation in the paper.

Guess, refine, repeat. Each loop, the cloud you get out is fed back in to build a better $\mathbf F$. After a handful of rounds the cloud stops changing — it has become self-consistent, agreeing with the field it creates. That converged cloud is the answer. This loop, run for every one of five million molecular poses, is the expensive heart of what nablaDFT computed.

The orbital, rebuilt in space — and the density matrix

Once the loop converges, you have the mixing numbers $C$, and you can rebuild the actual orbital at any point in space. The paper writes this (equation 6) using a shorthand called Einstein summation, where a repeated index secretly means "sum over it":

Equation 6 · rebuild the orbital in space $$ \psi^{\sigma}_{i}(\vec{r}_1) \;=\; C^{\sigma}_{\mu i}\,\phi_{\mu}(\vec{r}_1) $$

It says the same thing as equation 2 — "the orbital is the basis pieces $\phi_\mu$ weighted by the mixing numbers $C$" — now read at an actual location $\vec r_1$. (The repeated $\mu$ means "sum over all pieces.") From these weights we assemble the single most useful object in the whole method, the density matrix (equation 7):

Equation 7 · the density matrix $$ D^{\sigma}_{ij} \;=\; C^{\sigma}_{ik}\,C^{\sigma}_{jk} $$

Worked example · building a density matrix by hand

Suppose a molecule has one filled orbital, and after the SCF loop its mixing numbers came out as $\mathbf c = (0.8,\; 0.6)$ over two basis pieces. (Check: $0.8^2 + 0.6^2 = 1$ — properly normalized.)

Equation 7 with a single filled orbital is just $D_{ij} = c_i\,c_j$ — multiply the entries pairwise. So:

$D_{11} = 0.8\times0.8 = 0.64$, $\quad D_{22} = 0.6\times0.6 = 0.36$, $\quad D_{12} = D_{21} = 0.8\times0.6 = 0.48.$

The density matrix is $\mathbf D = \begin{bmatrix} 0.64 & 0.48 \\ 0.48 & 0.36 \end{bmatrix}.$ The diagonal (0.64, 0.36) says how much electron sits on each piece; the off-diagonal (0.48) measures the shared, bonding cloud between them.

Notice the diagonal sums to $1.00$ — exactly one electron's worth, conserved. The density matrix is the cloud, in number form.

Why we love the density matrix D

It's the electron density in compact form — literally how thickly the cloud is spread, which (unlike the wave function) is something you can measure in a lab.
It's built straight from the mixing numbers $C$ (the repeated $k$ sums over all filled orbitals).
Almost every property you'd want — energy, bond strengths, charges — can be squeezed out of $\mathbf D$. It is the bridge to the next chapter, where it becomes the single key ingredient in DFT's energy recipe.

CHAPTER FIVE The clever shortcut · paper eqs. (8)–(9)

DFT: stop chasing every electron

Term 3 from Chapter 3 — electrons pushing electrons — is what makes the wave function explode into that $10^{60}$-number monster. Density Functional Theory, or DFT, is a Nobel-winning dodge around it, and it's the method the entire nablaDFT paper is built on. The idea is so audacious it was barely believed at first.

The astonishing claim at the heart of DFT

Instead of tracking the full tangled wave function $\Psi$ — which depends on the positions of all the electrons at once — DFT asks a far humbler question:

"How thick is the electron cloud at each point?"

That thickness is the electron density, written $\rho(\vec r)$ ("rho"). It's just one number at each point in ordinary 3-D space — no matter how many electrons there are. Compare that to $\Psi$, which needs a separate coordinate for every electron. The density is dramatically simpler. The only question is: could something so simple possibly contain enough information? The shocking answer is yes.

The two Hohenberg–Kohn theorems (1964) — the bedrock DFT rests on two proven theorems by Hohenberg and Kohn. Theorem 1: the ground-state density $\rho(\vec r)$ uniquely determines everything about the system — the same density can't come from two different molecules. So the density secretly contains all the information the giant wave function did. Theorem 2: there exists an energy "functional" $E[\rho]$ that is lowest exactly at the true density — so you can find the right density by hunting for the one that minimizes the energy, the same "roll downhill" principle from Chapter 3. Together: the simple cloud is enough, and you can find it by minimizing energy. That's the whole licence for DFT.

Proving Theorem 1 — the density really does fix everything

It is worth seeing why Theorem 1 is true, because the proof is one of the most elegant short arguments in physics — a clean proof by contradiction — and because it shows the theorem isn't a hopeful approximation but an exact fact. First, the precise statement.

Hohenberg–Kohn Theorem I

The ground-state electron density $\rho(\vec r)$ determines the external potential $v(\vec r)$ (the nuclear arrangement) uniquely, up to an irrelevant additive constant. Since $v$ fixes the entire Hamiltonian, the density therefore determines the wave function and every property of the system.

The setup: two molecules might have different nuclear arrangements — different external potentials $v$ and $v'$, hence different Hamiltonians $\hat H$ and $\hat H'$ and different ground states $\psi$ and $\psi'$. We assume, for contradiction, that they nonetheless produce the same density $\rho$, and derive an impossibility.

Derivation · reductio ad absurdum

Suppose $v \ne v'$ (differing by more than a constant) but both give the same ground-state density $\rho$. Their Hamiltonians differ only in the potential: $\hat H' = \hat H + (\hat V' - \hat V)$.

The first term is just $E_0'$ (the ground-state energy of $\hat H'$). The second is where the density enters: the potential difference $\hat V'-\hat V$ is an ordinary function of position, and the average of any such function over a state is that function weighted by the density and integrated — $\langle\psi'|\hat V - \hat V'|\psi'\rangle = \int \big(v(\vec r) - v'(\vec r)\big)\rho(\vec r)\,d\vec r$, because $\rho$ is by definition "how much electron sits at each point." So $E_0 < E_0' + \int (v - v')\rho\,d\vec r$.

Now repeat the argument the other way around — use $\psi$ as a trial state for $\hat H'$: $E_0' < E_0 + \int (v' - v)\rho\,d\vec r$. (Same density $\rho$, by assumption, so the same integral appears with the opposite sign.)

Add the two inequalities. The density integrals are equal and opposite, so they cancel, leaving $E_0 + E_0' < E_0' + E_0$.

That says a number is strictly less than itself — impossible. The assumption $v \ne v'$ must be false.

one density cannot come from two different molecules

What this buys, exactly The chain is: $\rho \Rightarrow v \Rightarrow \hat H \Rightarrow \psi \Rightarrow$ everything. A single scalar function of three variables — the cloud thickness $\rho(x,y,z)$ — secretly encodes the same complete information as the monstrous $3N$-variable wave function from Chapter 3. This is the licence to throw away the wave function entirely and work only with the density. Theorem 2 then adds: among all candidate densities, the true one is the one that minimizes the energy functional — which is, once again, the variational principle of Chapter 3½, now rephrased in the language of $\rho$ rather than $\psi$.

The trade that makes chemistry computable. A swarm of mutually-reacting electrons (left) becomes one smooth cloud (right) — and the theorems guarantee no information is lost. This is why DFT handles thousand-electron molecules where exact methods choke at a dozen.

Kohn and Sham's brilliant cheat

There's still a problem. The theorems promise the magic functional $E[\rho]$ exists, but they don't tell you its formula — and the kinetic energy of a cloud of interacting electrons is itself fiendish to write down. Kohn and Sham's 1965 fix is the trick that made DFT actually usable, and it's beautifully sneaky.

Invent a pretend molecule. Imagine a fake system of electrons that don't interact at all — no term 3, no tangle — but which are rigged to have the exact same density as the real, messy molecule. Because the fake electrons are independent, each one obeys its own simple, solvable equation. We get nearly all the energy easily from the fake system, then sweep every difference between fake and real into a single correction term. That correction has a name you'll see in the equations: exchange–correlation.

Why this is genius, not cheating The hard part of the real problem — the electron tangle — is quarantined into one term, $E_{\text{xc}}$. Everything else becomes the easy, independent-electron problem we already know how to solve with the matrix machinery of Chapter 4. We've taken an impossible problem and isolated all of its difficulty into a single box. The catch, of course, is that nobody knows the exact contents of that box — which is the next section.

The Kohn–Sham equations, written down properly

Let's make the "pretend molecule" rigorous. We split the exact energy functional into pieces we can compute exactly and one piece we can't:

The Kohn–Sham energy decomposition $$ E[\rho] \;=\; \underbrace{T_s[\rho]}_{\substack{\text{kinetic energy of}\\\text{fake free electrons}}} \;+\; \underbrace{\int v_{\text{ext}}(\vec r)\,\rho(\vec r)\,d\vec r}_{\text{pull from nuclei}} \;+\; \underbrace{\tfrac12\!\int\!\!\int \frac{\rho(\vec r)\rho(\vec r\,')}{|\vec r - \vec r\,'|}\,d\vec r\,d\vec r\,'}_{\text{Hartree (cloud on itself)}} \;+\; \underbrace{E_{\text{xc}}[\rho]}_{\substack{\text{everything}\\\text{we got wrong}}} $$

The trick in one line

$T_s[\rho]$ is the kinetic energy of the non-interacting fake electrons — computable exactly from their orbitals. It is not the true kinetic energy, but it's most of it.
The external and Hartree terms are simple, exact integrals.
The definition of $E_{\text{xc}}$ is precisely "the exact energy minus those three pieces." It absorbs (a) the difference between true and fake kinetic energy, and (b) all the non-classical electron–electron effects. Small, but not zero — and not known exactly.

Now apply the variational principle one more time — minimize this $E[\rho]$ subject to keeping $N$ electrons. Because $\rho$ is built from the fake orbitals $\phi_i$ (with $\rho = \sum_i |\phi_i|^2$), minimizing over the orbitals yields a set of single-particle equations that look exactly like a Schrödinger equation for one electron moving in an effective potential. (The machinery is the same "set the derivative to zero" move you saw give Roothaan–Hall in Chapter 4 — here each derivative is a functional derivative, the continuous cousin defined just below, but the spirit is identical.) The result:

The Kohn–Sham equations $$ \left[\,-\tfrac{1}{2}\nabla^2 \;+\; v_{\text{eff}}(\vec r)\,\right]\phi_i(\vec r) \;=\; \varepsilon_i\,\phi_i(\vec r), \qquad v_{\text{eff}} = v_{\text{ext}} + v_{\text{H}} + v_{\text{xc}} $$

Reading the effective potential

$v_{\text{ext}}$ — the bare pull of the nuclei (ordinary Coulomb).
$v_{\text{H}}(\vec r) = \int \dfrac{\rho(\vec r\,')}{|\vec r - \vec r\,'|}d\vec r\,'$ — the Hartree potential, the average electrostatic push from the whole cloud.
$v_{\text{xc}}(\vec r) = \dfrac{\delta E_{\text{xc}}[\rho]}{\delta \rho(\vec r)}$ — the exchange–correlation potential, defined as the functional derivative of $E_{\text{xc}}$ with respect to the density. (A functional derivative answers: "if I add a tiny bit of density right here at $\vec r$, how much does the energy change?" — the continuous analogue of a partial derivative.)

The catch, made precise — and the self-consistency loop returns Look carefully: $v_{\text{eff}}$ depends on $\rho$ (through $v_{\text{H}}$ and $v_{\text{xc}}$), but $\rho$ is built from the orbitals $\phi_i$ that solve the equation — which need $v_{\text{eff}}$ to be written down in the first place. This is precisely the chicken-and-egg of Chapter 4, now in DFT clothing. The cure is identical: the Self-Consistent Field loop. Guess $\rho$, build $v_{\text{eff}}$, solve for orbitals, rebuild $\rho$, repeat until it stops changing. Every DFT calculation in both papers is this loop, and the only genuinely unknown ingredient in it is $v_{\text{xc}}$.

The DFT energy machine, rebuilt in three parts

With the Kohn–Sham trick, the energy grid $\mathbf F$ from Chapter 4 (now properly called the Kohn–Sham matrix) splits into three understandable pieces. The paper writes equation 8:

Equation 8 · the Kohn–Sham matrix, in three parts $$ F^{\sigma}_{ij} \;=\; H^{\text{c}\,\sigma}_{ij} \;+\; J^{\sigma}_{ij} \;+\; V^{\text{xc}}_{ij} $$

The three ingredients

$H^{\text{c}}_{ij}$ — the core part: each electron's own motion energy plus its attraction to the nuclei (terms 1 and 2 from Chapter 3's full Hamiltonian). The easy, well-understood stuff.
$J^{\sigma}_{ij}$ — the Coulomb / Hartree part: the average electric push of the whole electron cloud on itself. Written as an integral it's $J \sim \int\int \frac{\rho(\vec r)\,\rho(\vec r')}{|\vec r - \vec r'|}\,dV\,dV'$ — "every patch of cloud pushing every other patch," straight from Coulomb's law.
$V^{\text{xc}}_{ij}$ — the exchange–correlation part: the magic correction box from the Kohn–Sham trick, holding all the subtle quantum effects the first two parts miss. The quality of your whole answer lives or dies on this one term.

A subtle bug, and what "exchange" fixes The Coulomb term $J$ has a hidden flaw: it lets the cloud push against itself, including the absurdity of a single electron repelling its own smeared-out charge — a "self-interaction" that shouldn't exist. Part of the exchange–correlation term's job is to cancel exactly this error. "Exchange" also encodes a deep quantum rule (from the Pauli principle of Chapter 2) that same-spin electrons instinctively avoid one another. "Correlation" captures the gentler way all electrons dodge each other in real time. Bundling these effects is why the term is hard — and why getting it right is a whole science.

The exchange–correlation zoo, and where the paper lives

Since nobody knows the exact $E_{\text{xc}}$, scientists have built a ladder of ever-better approximations — physicists half-jokingly call it "Jacob's Ladder," climbing from earth toward the heaven of the exact answer. Each rung feeds the functional more information about the density:

Climbing toward the exact answer. Higher rungs are more accurate but more expensive. nablaDFT uses a hybrid functional, ωB97X-D — a high, costly rung that mixes in a dose of "exact exchange." The -D tacked on the end stands for a dispersion correction: a patch for the gentle stickiness (van der Waals forces) that even hybrids otherwise miss. This is a deliberately accurate, expensive choice — part of why the dataset weighs 100 terabytes.

What the bottom rung actually is, derived

The lowest rung, the Local Density Approximation (LDA), is worth deriving because its exchange energy can be pinned down exactly — and the form is surprisingly clean. The idea: at each point in the molecule, pretend the electron gas has the uniform density found right there, and use the known exchange energy of a uniform gas. The remarkable part is that we can find the density-dependence by pure scaling, without any hard integral.

Derivation · the 4/3 power in LDA exchange

Exchange energy comes from the electron–electron interaction $\sim 1/r$, which has units of (1 / length). Density $\rho$ has units of (1 / length$^3$), so a characteristic length scale at density $\rho$ is $\ell \sim \rho^{-1/3}$.

The exchange energy per unit volume must be built from $\rho$ alone and carry units of (energy / volume) = (1/length · 1/length$^3$) = (1 / length$^4$) in atomic units. So it scales as $\ell^{-4} \sim (\rho^{-1/3})^{-4} = \rho^{4/3}$.

Therefore the exchange energy density is proportional to $\rho^{4/3}$, and integrating over space gives the total: $E_x^{\text{LDA}}[\rho] = -C_x \displaystyle\int \rho(\vec r)^{4/3}\,d\vec r$.

The constant works out (from the full uniform-gas calculation) to $C_x = \tfrac{3}{4}\left(\tfrac{3}{\pi}\right)^{1/3}$. Dimensional analysis gave us the entire functional form; only the prefactor needed the detailed integral.

exchange goes as the 4/3 power of density

Every higher rung of the ladder refines this. GGA adds a dependence on the density's gradient $\nabla\rho$ (how fast the cloud is changing, not just how thick it is). Meta-GGA adds the kinetic-energy density $\tau$. And hybrids mix in a fraction of "exact exchange" computed the expensive Hartree–Fock way. Each rung feeds the functional strictly more information about the local shape of the cloud.

The bug that the functional must fix: self-interaction

Here's a concrete flaw worth seeing explicitly, because it explains why exchange is non-negotiable. The Hartree term $J[\rho]$ — the cloud's electrostatic push on itself — is summed over the whole density. But consider a hydrogen atom, which has exactly one electron. That single electron has no other electron to repel. Yet:

Derivation · self-interaction error in one-electron hydrogen

For one electron, the density is just $\rho = |\phi|^2$ from its single orbital. The Hartree energy is $J[\rho] = \tfrac12\displaystyle\int\!\!\int \dfrac{\rho(\vec r)\rho(\vec r\,')}{|\vec r - \vec r\,'|}\,d\vec r\,d\vec r\,'$.

This integral is manifestly positive — it's the electrostatic self-repulsion of a positive charge distribution. For hydrogen's ground state it evaluates to $J = \tfrac{5}{16}$ Hartree, a sizeable spurious energy.

But physically there is nothing for the lone electron to repel. The true electron–electron energy is exactly zero. So $J$ has introduced a fake $+\tfrac{5}{16}$ Hartree of "an electron repelling itself."

The exchange energy must therefore exactly cancel this for a one-electron system: $E_x = -J$. A correct functional has to undo its own Hartree over-counting.

exchange exists, in part, to delete self-repulsion

Anatomy of ωB97X-D — decoding the paper's exact choice Now the dataset's functional name is fully readable. The trouble with hybrids is choosing how much exact exchange to mix in — too little and self-interaction leaks back, too much and other errors grow. Range separation is the clever fix, and it's what the $\omega$ stands for. Split the Coulomb interaction into a short-range and a long-range part using the error function: $$ \frac{1}{r} = \underbrace{\frac{1 - \operatorname{erf}(\omega r)}{r}}_{\text{short range}} + \underbrace{\frac{\operatorname{erf}(\omega r)}{r}}_{\text{long range}} $$ Use cheap density-functional exchange up close and exact (Hartree–Fock) exchange far away — where self-interaction error does the most damage. The parameter $\omega$ sets the crossover distance. B97 is the underlying functional family (Becke 1997); X marks the exchange variant; and -D bolts on an explicit dispersion correction of the form $-\sum_{AB} C_6^{AB}/R_{AB}^6$ — the attractive van der Waals tail between atom pairs $A,B$ that no local functional captures on its own. Every one of these pieces is a deliberate accuracy choice, and together they are why a single ∇²DFT conformation can take CPU-hours and the dataset fills hundreds of terabytes.

The grand total: the energy equation

Finally, the molecule's complete energy. The paper writes it as equation 9. It looks dense, but every piece is now an old friend — we've spaced it out so you can see the seams:

Equation 9 · total energy of the molecule $$ E_{\text{total}} \;=\; \underbrace{D^{T}_{ij}\big(T_{ij}+V_{ij}\big)}_{\text{motion + nuclear pull}} \;+\; \underbrace{\tfrac{1}{2}\,D^{T}_{ij}\,D^{T}_{\lambda\beta}\,(ij\,\vert\,\lambda\beta)}_{\text{cloud pushing on itself}} \;+\; \underbrace{E_{\text{xc}}\big[\rho_\alpha,\rho_\beta\big]}_{\text{quantum weirdness}} $$

Read the three lumps left to right

First lump — the density matrix $D$ (from Chapter 4!) multiplied by $T$ (kinetic / motion energy) plus $V$ (the nuclei's pull). "Energy of moving electrons and their attraction to nuclei." The little $T$ superscript just means the matrix is transposed — a bookkeeping flip so the grids multiply correctly.
Second lump — the density multiplied by itself, with $(ij\,\vert\,\lambda\beta)$ measuring how two patches of cloud repel. The $\tfrac12$ stops us double-counting each pushing-pair. This is the Coulomb $J$ term in full.
Third lump — $E_{\text{xc}}[\rho_\alpha,\rho_\beta]$, the exchange–correlation energy, fed the spin-up density $\rho_\alpha$ and spin-down density $\rho_\beta$ (there's that spin bookkeeping again). The square brackets are the tell: this is a functional — a machine that eats an entire function (the density) and returns one number. That is literally where "Density Functional Theory" gets its name.

Notice the star of Chapter 4 — the density matrix $D$ — appears in every single term. Get the density, and the total energy falls right out. That's the engine. And the whole point of nablaDFT is to make this engine run without the crushing cost of the SCF loop — by teaching a machine to skip straight to the answer.

CHAPTER SIX The object we want to predict · paper eq. (10)

The matrix of everything

Everything so far has been climbing toward one object. The grand prize of the whole calculation — the thing nablaDFT ultimately wants a computer to guess — is the Hamiltonian matrix $\mathbf{H}$. Hold it, and the energy, the orbitals, the density, and a dozen other properties all unlock. Let's look at how it's organised, why its structure matters, and what its hidden eigenvalues tell you about the real world.

It's a grid of blocks, one per pair of atoms

The paper draws $\mathbf H$ (equation 10) as a big grid that is itself carved into smaller rectangles called blocks:

Equation 10 · the Hamiltonian in block form $$ \mathbf{H} = \begin{bmatrix} \mathbf{H}_{11} & \cdots & \mathbf{H}_{1j} & \cdots & \mathbf{H}_{1n}\\ \vdots & \ddots & \vdots & & \vdots \\ \mathbf{H}_{i1} & \cdots & \mathbf{H}_{ij} & \cdots & \mathbf{H}_{in}\\ \vdots & & \vdots & \ddots & \vdots \\ \mathbf{H}_{n1} & \cdots & \mathbf{H}_{nj} & \cdots & \mathbf{H}_{nn} \end{bmatrix} $$

How to read this grid of grids

There are $n$ atoms. The grid is $n$ blocks across and $n$ blocks down.
Each block $\mathbf H_{ij}$ describes the relationship between atom $i$ and atom $j$ — how their electrons feel one another.
Diagonal blocks ($\mathbf H_{11}, \mathbf H_{22}, \dots$) describe each atom with itself; off-diagonal blocks describe pairs.
The paper notes a block $\mathbf H_{ij}$ lives in $\mathbb R^{\,n_{\text{ao},i}\times n_{\text{ao},j}}$ — fancy notation meaning its size depends on how many basis pieces ($n_{\text{ao}}$, "atomic orbitals") each atom brings. Big bromine brings many; small hydrogen brings few. So the rectangles are different sizes.

The molecule's whole quantum identity, as a grid. Gold block = atom 1 with itself (biggest, because the gold atom carries the most basis pieces). Grey blocks describe pairs. To predict the matrix, predict each block — and each block is a question about a pair of atoms.

Why the paper stresses the block structure The paper says this block layout "underlies the construction of interaction modules in the models." Plainly: because each block is about a pair of atoms, the neural networks in Chapter 7 are built to think in pairs too. They learn a description of each atom and its neighbourhood, then combine pairs of those descriptions to fill in each block. The shape of the answer dictates the shape of the machine.

A free gift: the matrix is symmetric

Here's a small mercy hiding in the structure. The influence of atom $i$ on atom $j$ is the same as the influence of $j$ on $i$ — so $H_{ij} = H_{ji}$. The matrix is a perfect mirror image across its diagonal (mathematicians say Hermitian, or for real numbers, symmetric). That means you only ever need to work out half of it — the bottom triangle — and the top half comes free. A useful saving, and the neural networks exploit it.

The payoff: eigenvalues are orbital energies

Now for why this matrix is worth so much. Remember Chapter 3: solving $\mathbf F\mathbf c = \varepsilon\mathbf S\mathbf c$ means finding the eigenvalues of the matrix. Each eigenvalue $\varepsilon$ is the energy of one orbital. So the matrix secretly contains a whole ladder of energy levels — and electrons fill that ladder from the bottom up, two per rung (one spin-up, one spin-down, per Chapter 2's Pauli rule).

Two rungs on that ladder have special names and enormous importance:

HOMO and LUMO — the two rungs that matter most

HOMO — the Highest Occupied Molecular Orbital: the top filled rung, where the most loosely-held electrons sit. These are the electrons a molecule offers up in a reaction.
LUMO — the Lowest Unoccupied Molecular Orbital: the lowest empty rung, the first place an incoming electron would land.
The gap between them — the HOMO–LUMO gap — is one of the most predictive numbers in chemistry. A small gap means the molecule is reactive and easily excited; a large gap means it's stable and inert.

That gap even sets a molecule's colour. A photon of light with just the right energy can kick an electron from HOMO to LUMO; the molecule absorbs that colour and reflects the rest. The size of the gap picks which colour gets absorbed. The redness of a tomato and the green of a leaf are HOMO–LUMO gaps you can see. Let's compute one.

Worked example · reading energies out of a tiny matrix

Suppose a toy molecule's Hamiltonian (in some energy unit, after accounting for overlap) came out as $\mathbf H = \begin{bmatrix} -10 & -2 \\ -2 & -10 \end{bmatrix}.$ We want its orbital energies — the eigenvalues.

Use the method from Chapter 3: set $\det(\mathbf H - \varepsilon\mathbf I) = 0$, giving $(-10-\varepsilon)^2 - (-2)^2 = 0.$

So $(-10-\varepsilon)^2 = 4$, meaning $-10-\varepsilon = \pm 2$, giving $\varepsilon = -12$ or $\varepsilon = -8.$

Two rungs: a low one at $-12$ (the bonding orbital) and a higher one at $-8$ (the antibonding orbital). With two electrons, both drop into the lower rung — so this molecule is bonded and stable, its energy lowered by forming the bond.

If these were the HOMO ($-12$, filled) and LUMO ($-8$, empty), the HOMO–LUMO gap is $-8-(-12) = 4$ energy units. That single number — which we pulled straight out of a 2×2 grid — predicts how reactive the molecule is and what colour of light it absorbs. This is why the matrix is the grand prize.

Fill the energy ladder yourself below. Add electrons, watch them stack two-per-rung, and see the HOMO–LUMO gap — and the molecule's predicted colour — change.

Try it: fill the orbital ladder

Electrons fill from the bottom, two per rung (↑↓). The gap between the highest filled rung (HOMO) and lowest empty one (LUMO) sets reactivity and colour. Change the number of electrons and the spacing of the levels.

Electrons 6

Gap size 3

The catch that makes AI worth it — and the wall

The paper adds a crucial warning: the entries of this matrix are not smooth. Nudge an atom a hair, and the numbers can jump in jagged, hard-to-predict ways. That bumpiness is exactly what flexible machine-learning models handle well and rigid formulas handle badly.

And there's the cost. Solving for these eigenvalues by the honest SCF route scales roughly as $N^3$ — triple the atoms and the work grows about 27-fold. Build the demo below to feel that wall, and to see why a fast neural-network guess is worth so much.

Try it: how big does the problem get?

Every atom you add doesn't just add a row — it adds a whole row and column of blocks, and the solve-cost grows like $N^3$. Drag to add atoms and watch both explode.

Atoms in molecule 3

That $N^3$ wall is the whole motivation for the next chapter. If a machine could look at a molecule and simply predict this matrix — skipping the SCF loop entirely — chemistry would speed up by orders of magnitude. So: can it?

CHAPTER SIX ½ Inside the grand prize · a rigorous interlude

Anatomy of the Hamiltonian matrix

Chapter 6 showed that the Hamiltonian matrix is the object both papers ultimately want to predict, and that it's organised in blocks, one per pair of atoms. That was the sketch. This chapter is the dissection: where each number in the matrix actually comes from, why the matrix has the symmetries it has, why most of it is nearly zero — and crucially, why that near-emptiness is the entire reason a neural network can predict it at all. We'll also finish the eigenvalue story properly, with the orthogonalization step real codes use. Like Chapter 3½, this is a rigorous detour; the reward is that the machines of Chapter 7 will make complete structural sense.

Each basis function has a shape and a direction

Recall from Chapter 4 that we expand orbitals in basis functions $\phi$. We can now be precise about what those functions look like. Each one factorizes into a radial part (how it fades with distance from its atom) times an angular part (its shape and orientation):

The structure of a basis function $$ \phi_{n\ell m}(\vec r) \;=\; \underbrace{R_{n\ell}(r)}_{\text{radial: how far}} \;\times\; \underbrace{Y_{\ell m}(\theta,\varphi)}_{\text{angular: what shape}} $$

The angular momentum label $\ell$

$Y_{\ell m}$ are the spherical harmonics — the same ones that will star in Chapter 7's PhiSNet. The number $\ell = 0, 1, 2, \dots$ is the angular momentum, and it sets the shape you met in Chapter 2:
$\ell = 0$: the s shape — a sphere. There is $2\ell+1 = 1$ of them.
$\ell = 1$: the p shapes — dumbbells. There are $2\ell+1 = 3$ (pointing along $x$, $y$, $z$).
$\ell = 2$: the d shapes — cloverleaves. There are $2\ell+1 = 5$.
A complete set of $2\ell+1$ functions at a given $\ell$ is called a shell. The key fact: a shell is exactly what you need so the set transforms cleanly into itself when you rotate the molecule. Rotating a single $p_x$ lobe turns it into a mix of $p_x, p_y, p_z$ — never out of the shell.

This is why the matrix blocks have the sizes they do. The block $\mathbf H_{ij}$ coupling atoms $i$ and $j$ is a rectangle whose dimensions are (number of basis functions on $i$) × (number on $j$). Let's count for real atoms in the def2-SVP basis the papers use.

Worked example · the real size of a water molecule's matrix

In def2-SVP, hydrogen gets 2 s-shells + 1 p-shell = $2(1) + 1(3) = 5$ basis functions. Oxygen gets 3 s-shells + 2 p-shells + 1 d-shell = $3(1) + 2(3) + 1(5) = 14$ basis functions.

Water is H$_2$O: two hydrogens and one oxygen. Total basis size $= 14 + 5 + 5 = 24$ functions.

So water's full Hamiltonian matrix is $24 \times 24 = 576$ numbers. The oxygen–oxygen diagonal block is $14\times14$; each O–H block is $14\times5$; the H–H block is $5\times5$.

The matrix is carved into these atom-pair rectangles of different sizes — exactly the ragged block structure drawn in Chapter 6, now with real dimensions. A big drug molecule with 40 heavy atoms reaches matrices of many hundreds of rows, which is why storing and predicting them is a serious undertaking.

Why the matrix is symmetric — and what that guarantees

Chapter 6 mentioned the matrix is symmetric as a "free gift." Here's the rigorous reason. The Hamiltonian operator $\hat H$ is self-adjoint (Hermitian) — a physical requirement, because it guarantees energies come out as real numbers rather than imaginary ones. Self-adjointness of the operator forces a symmetry on its matrix:

Derivation · Hermiticity of the matrix

By definition, $H_{ij} = \langle\phi_i|\hat H|\phi_j\rangle$ and $H_{ji} = \langle\phi_j|\hat H|\phi_i\rangle$.

Self-adjointness means $\langle\phi_i|\hat H|\phi_j\rangle = \langle\hat H\phi_i|\phi_j\rangle = \overline{\langle\phi_j|\hat H|\phi_i\rangle}$ (the overbar is complex conjugation).

So $H_{ij} = \overline{H_{ji}}$: the matrix equals its own conjugate-transpose, $\mathbf H^\dagger = \mathbf H$. It is Hermitian.

With the real-valued Gaussian basis used in practice, the conjugation does nothing, so simply $H_{ij} = H_{ji}$: the matrix is real and symmetric, a mirror image across its diagonal.

only the lower triangle is independent

Two payoffs from one symmetry First, the practical one: you only need to compute and store half the matrix — the diagonal and the lower triangle — which is exactly the saving the matrix-predicting networks exploit. Second, the deep one: a real symmetric matrix is guaranteed by the spectral theorem to have a complete set of real eigenvalues and orthogonal eigenvectors. Those real eigenvalues are the orbital energies. The physics (real energies) and the linear algebra (real spectrum) are the same statement viewed from two sides.

Why most of the matrix is almost zero — and why that matters most of all

Here is the single most important structural fact for understanding why machine learning works on this problem. The off-diagonal blocks of the matrix — the couplings between distant atoms — are not just small. They die off exponentially with distance. Let's see why, from the Gaussian basis.

Derivation · exponential decay of distant couplings

A basis function on atom $A$ is a Gaussian $\phi_A \sim e^{-\alpha|\vec r - \vec R_A|^2}$; one on atom $B$ is $\phi_B \sim e^{-\beta|\vec r - \vec R_B|^2}$. The overlap and Hamiltonian elements between them involve the product $\phi_A\,\phi_B$.

The product of two Gaussians is itself a Gaussian (the "Gaussian product theorem"), centred between the atoms — and crucially, multiplied by a constant prefactor that depends on the separation: $\phi_A\phi_B \propto \exp\!\left(-\dfrac{\alpha\beta}{\alpha+\beta}|\vec R_A - \vec R_B|^2\right)\times(\text{Gaussian in }\vec r)$.

That prefactor is the whole story. Writing $R_{AB} = |\vec R_A - \vec R_B|$ for the atom separation, the matrix element scales as $H_{AB} \sim \exp\!\left(-\dfrac{\alpha\beta}{\alpha+\beta}R_{AB}^2\right)$.

This falls off as a Gaussian in the distance: double the separation and the coupling drops not by a factor but by an exponential of the squared distance. Atoms a few bond-lengths apart have essentially zero direct matrix element.

distant atom-pairs barely couple

"Nearsightedness of electronic matter" — and the licence for machine learning Walter Kohn (the same Kohn of Kohn–Sham) elevated this observation into a principle he called the nearsightedness of electronic matter: what happens to the electrons at one point depends only on the environment within a finite range around it; the far-away parts of a molecule are invisible to it. The Hamiltonian matrix is sparse precisely because of this nearsightedness — the meaningful numbers cluster near the diagonal and between bonded neighbours, while distant blocks fade to nothing. This is the deep reason the neural networks of Chapter 7 can work. A message-passing network only ever lets each atom talk to its near neighbours — and that is exactly the right inductive bias, because the true matrix is local. A model built on locality is not a crude approximation here; it matches the physics. If electronic matter were not nearsighted, no local network could ever predict these matrices.

Finishing the eigenvalue story: the generalized problem

Chapter 6 read orbital energies off the matrix as eigenvalues. But there's a wrinkle we can now resolve rigorously. Because the basis isn't orthonormal — the overlap matrix $\mathbf S$ isn't the identity (Chapter 4) — the equation to solve isn't the textbook $\mathbf H\mathbf c = \varepsilon\mathbf c$ but the generalized eigenproblem $\mathbf H\mathbf c = \varepsilon\mathbf S\mathbf c$. Here is how real codes turn it back into a standard one.

Derivation · Löwdin orthogonalization

First, what $\mathbf S^{-1/2}$ even means. A symmetric matrix can be rotated into a frame where it's diagonal — just a list of (positive) numbers down the diagonal. In that frame, "inverse square root" is the obvious thing: replace each diagonal number $s$ by $1/\sqrt{s}$. Rotate back, and you have a genuine matrix $\mathbf S^{-1/2}$ with the defining property $\mathbf S^{-1/2}\mathbf S^{-1/2}\mathbf S = \mathbf I$. (This works because $\mathbf S$, the overlap matrix, is symmetric and positive-definite — all its diagonal numbers are positive, so the square roots are real.)

Start from $\mathbf H\mathbf c = \varepsilon\mathbf S\mathbf c$. Left-multiply by $\mathbf S^{-1/2}$ and insert $\mathbf S^{-1/2}\mathbf S^{1/2}$ before $\mathbf c$: $\mathbf S^{-1/2}\mathbf H\mathbf S^{-1/2}\,(\mathbf S^{1/2}\mathbf c) = \varepsilon\,\mathbf S^{-1/2}\mathbf S\mathbf S^{-1/2}\,(\mathbf S^{1/2}\mathbf c)$.

The right side simplifies, because $\mathbf S^{-1/2}\mathbf S\,\mathbf S^{-1/2} = \mathbf I$. Define a transformed matrix $\tilde{\mathbf H} = \mathbf S^{-1/2}\mathbf H\mathbf S^{-1/2}$ and transformed vector $\tilde{\mathbf c} = \mathbf S^{1/2}\mathbf c$.

The equation becomes $\tilde{\mathbf H}\,\tilde{\mathbf c} = \varepsilon\,\tilde{\mathbf c}$ — an ordinary symmetric eigenproblem. Solve it with standard tools, then recover the real coefficients via $\mathbf c = \mathbf S^{-1/2}\tilde{\mathbf c}$.

the overlap is absorbed; orbital energies fall out

The eigenvalues $\varepsilon_i$ of $\tilde{\mathbf H}$ are the orbital energies. Order them, fill two electrons into each from the bottom (the aufbau principle, enforced by Pauli from Chapter 2), and the highest filled and lowest empty levels are the HOMO and LUMO whose gap you computed by hand in Chapter 6. This orthogonalization-then-diagonalization is the step that costs $O(N^3)$ — the matrix operations grow as the cube of the basis size — which is the precise origin of the "$N^3$ wall" the interactive demo let you feel.

Reading chemistry straight off the matrix

Finally, the payoff that makes the matrix worth predicting: once you have it (and the density matrix $\mathbf D$ built from the occupied eigenvectors, Chapter 4), real chemical quantities come out by simple operations. Two rigorous examples.

Derivation · the density matrix counts the electrons

The density matrix is $\mathbf D = \sum_{i \in \text{occ}} \mathbf c_i \mathbf c_i^\dagger$, summing over occupied orbitals (Chapter 4's equation 7). A key identity follows from how it pairs with the overlap.

Consider the trace of the product $\mathbf D\mathbf S$. Using the cyclic identity from Tool 5 (you can rotate factors inside a trace), $\operatorname{Tr}(\mathbf D\mathbf S) = \sum_{i\in\text{occ}} \operatorname{Tr}(\mathbf c_i\mathbf c_i^\dagger \mathbf S) = \sum_{i\in\text{occ}} \mathbf c_i^\dagger \mathbf S\,\mathbf c_i$ — the rotation turns each outer-product-times-$\mathbf S$ into a plain number.

But $\mathbf c_i^\dagger\mathbf S\,\mathbf c_i = \langle\psi_i|\psi_i\rangle = 1$ — each occupied orbital is normalized (that's what $\mathbf S$ measures, from Chapter 4).

So $\operatorname{Tr}(\mathbf D\mathbf S) = \sum_{i\in\text{occ}} 1 = N_{\text{occ}}$ — exactly the number of occupied orbitals, hence (×2 for spin) the number of electrons. The matrix knows how many electrons it holds.

Tr(DS) = electron count

More chemistry hiding in D and S

Atomic charges. Splitting $\operatorname{Tr}(\mathbf D\mathbf S)$ by which atom each basis function sits on gives the electrons "belonging" to each atom — the Mulliken population. Subtract from the nuclear charge and you have the partial charge on every atom, the thing that drives how molecules attract and dock.
Bond orders. The off-diagonal blocks of $\mathbf D$ between two atoms measure their shared electron density — large means a strong bond, near-zero means none. This is a direct, quantitative read-out of the molecular structure.
Idempotency. The density matrix satisfies $\mathbf D\mathbf S\mathbf D = \mathbf D$ — a mathematical fingerprint that it describes a valid set of fully-occupied orbitals. Codes check this as another quality control.

Why all this anatomy matters for the papers Step back and see what the structure bought us. The Hamiltonian matrix is symmetric (so predict half), built from angular-momentum shells (so it transforms predictably under rotation — the hook for PhiSNet's equivariance), and sparse with exponential decay (so a local, neighbour-based network is the physically correct architecture). Every design choice in the Chapter 7 models is a response to a structural fact proven in this chapter. When ∇²DFT reports that PhiSNet predicts these matrices well but struggles to generalize across diverse molecules, it's measuring how well a network has learned this structure — symmetry, locality, and all — from data alone.

CHAPTER SEVEN The guessing machines · paper eqs. (11)–(15)

Teaching a computer to skip the hard sum

Solving DFT honestly means grinding the SCF loop from Chapter 4 until it converges — hours for a medium molecule. nablaDFT's plan: pay that cost once for a giant pile of molecules, then train a neural network to learn the pattern from input (the atoms and their positions) to output (the energy or the matrix). Afterward, a brand-new molecule's answer takes a heartbeat. This chapter meets the machines, and the genuine physics built into each one.

What a neural network really is, and how it learns

A neural network is a long chain of simple steps. At each step it takes numbers, multiplies and adds them in adjustable ways (the adjustable numbers are called weights), and passes the result on. None of the physics is programmed in. Instead, the network learns by example.

Learning works like this. You show the network a molecule whose true energy you already computed with DFT. The network makes a guess. You measure how wrong it was with a loss function — typically the squared error:

The loss · how wrong the guess is $$ \mathcal{L} \;=\; \big(E_{\text{guess}} - E_{\text{true}}\big)^2 $$

Then comes the clever part, gradient descent. Think of the loss as a landscape, with the network's thousands of weights as the coordinates and the height as how-wrong-you-are. You want the lowest valley. So you feel the slope under your feet — the gradient, our friend $\nabla$ again — and take a small step downhill. Repeat millions of times, across millions of molecules, and the network slowly descends into a valley where its guesses are good.

Gradient descent · step downhill on the loss $$ w_{\text{new}} \;=\; w_{\text{old}} \;-\; \eta\,\frac{\partial \mathcal{L}}{\partial w} $$

Decode the learning step

$w$ — one of the network's adjustable weights (it has millions; this rule applies to each).
$\frac{\partial \mathcal L}{\partial w}$ — the slope of the loss with respect to that weight: "if I nudge this knob, does the error go up or down, and how fast?"
The minus sign — step against the slope, i.e. downhill, toward less error.
$\eta$ ("eta") — the learning rate, the size of each step. Too big and you overshoot the valley; too small and training crawls.

Training is rolling downhill. Each step, the network feels the slope of its error and nudges its weights to reduce it. After enough steps over enough examples, it settles where its predictions match the DFT truth. The slope it feels is a gradient — the same $\nabla$ that named the paper.

Every model in the paper shares the same four-stage skeleton; they differ in how they build the third stage, the interaction, where atoms share information.

The shared blueprint. Atoms in as numbers, then rounds of "interaction" where each atom updates itself based on its neighbours, then a readout. The models are variations on the interaction step.

The big idea underneath: message passing on a graph

Picture the molecule as a graph — atoms are dots, and nearby atoms are joined by lines. The networks learn by playing a kind of telephone game on this graph: every atom sends a little "message" (a bundle of numbers describing itself) to its neighbours, each atom gathers the messages it receives and updates its own description, and you repeat for a few rounds. After round one, each atom knows about its immediate neighbours; after round two, about neighbours-of-neighbours; and so on, until information has rippled across the whole molecule. This is called message passing, and it's the shared soul of all four models.

Model 1 — SchNet: let distance do the talking

SchNet's insight: how much one atom influences another should depend on how far apart they are. It writes the update as equation 11:

Equation 11 · SchNet's distance-aware update $$ x^{\,l+1}_{i} \;=\; \big(X^{l} * W^{l}\big)_{i} \;=\; \sum_{j} x^{l}_{j} \circ W^{l}\!\big(r_i - r_j\big) $$

Decode SchNet

$x^{l}_{i}$ — the bundle of numbers describing atom $i$ at round $l$. The $l+1$ on the left is the updated bundle.
$\sum_j$ — "add up the messages from every other atom $j$." Atom $i$ listens to all its neighbours.
$W^{l}(r_i - r_j)$ — a learned filter that depends only on the distance between atoms $i$ and $j$. Close neighbour → strong influence; far → weak. This is the heart of it.
$\circ$ — multiply element by element (a way of blending two lists of numbers).

How a distance becomes numbers: the RBF trick

A subtlety: you can't just feed the raw distance (say, "1.4 ångström") into the network — a single number is too blunt, and the network learns smoother patterns if distance is spread across many gentle features. So SchNet expands each distance into a whole bank of overlapping bell-shaped detectors, called radial basis functions (RBFs):

RBF expansion · one distance → a bank of features $$ e_k(d) \;=\; \exp\!\left(-\,\frac{(d - \mu_k)^2}{2\sigma^2}\right) $$

Decode the RBF

$d$ — the distance between two atoms.
$\mu_k$ — the centre of detector number $k$. You place a row of them at, say, 0.5 Å, 1.0 Å, 1.5 Å, and so on.
Each detector "lights up" most when the distance $d$ is near its own centre $\mu_k$ and fades smoothly away otherwise — that's the bell shape (a Gaussian, like our basis functions in Chapter 4!).
The point: a distance of 1.4 Å becomes a smooth pattern of which detectors fired and how strongly. Smooth inputs → smooth, learnable filters → and crucially, the result can be cleanly differentiated (which we'll need for forces).

Why SchNet's answers don't change when you spin the molecule

Because the filter depends only on distance — never on which way the molecule points — SchNet's predictions are unchanged if you rotate the whole molecule. This property is called invariance, and it's exactly what real physics demands: a molecule's energy doesn't care which way it's facing. Let's make "invariance" precise, because the next models hinge on a sharper version of it.

Invariance vs equivariance — the crucial distinction Let $R$ be a rotation (a matrix that spins coordinates). Apply it to every atom's position.

Invariant quantity: doesn't change at all. Energy is invariant — $E(R\cdot\text{molecule}) = E(\text{molecule})$. Spin the molecule, same energy. A single number, so "doesn't change" is the only sensible behaviour.
Equivariant quantity: changes in step with the rotation. A force is an arrow; spin the molecule and the force arrow must spin the same way — $\vec F(R\cdot\text{molecule}) = R\cdot\vec F(\text{molecule})$. It doesn't stay fixed; it co-rotates.

SchNet handles invariant energies. But the Hamiltonian matrix is equivariant — spin the molecule and the matrix must transform along with it, in a precise way. A model that gets this wrong gives different answers for the same molecule at different tilts. That's the problem the last model solves.

The bend that makes learning possible

Between rounds, SchNet needs to "bend" its numbers so it can learn curved relationships — a perfectly straight machine can only ever learn straight lines. It uses a smooth bending function, the shifted softplus (equation 12):

Equation 12 · the shifted-softplus bend $$ \operatorname{ssp}(x) \;=\; \ln\!\big(0.5\,e^{x} + 0.5\big) $$

What this little function does

It's perfectly smooth — no sharp corners. This matters intensely because the paper computes forces by differentiating the energy, and you cannot cleanly differentiate across a corner. (More on forces just below.)
For big positive $x$ it rises almost like a straight line; for big negative $x$ it flattens toward zero. That gentle kink is what lets a stack of them mould to any curve.
The $0.5$'s and the $\ln$ are tuned so that $\operatorname{ssp}(0)=0$ — a tidy zero point. That's the "shifted" part.

The bend that lets learning happen. Flat-left, straight-right, smoothly curved between. Stack thousands of these and the network can shape itself to almost any pattern in the data.

Forces close the loop with Chapter 3

Why does smoothness matter so much? Because a force is how energy changes as you move an atom — and "how something changes as you move" is a derivative. In physics, force is the negative gradient of energy:

Force is the slope of energy $$ \vec{F}_i \;=\; -\,\nabla_i E $$

So once the network can predict energy $E$, it gets the forces for free by differentiating — but only if every step in the network is smooth enough to differentiate. That's the whole reason for the gentle shifted-softplus instead of a function with corners. The $\nabla$ that named the paper shows up one last time, right here, turning predicted energies into the forces that push atoms around.

Model 2 — DimeNet: distances aren't enough, you need angles

DimeNet spotted a gap in SchNet. Distances alone can make two genuinely different molecules look identical to the machine. The fix: also pay attention to angles — the bend formed by triples of atoms. DimeNet passes messages that carry this directional information (equation 13):

Equation 13 · DimeNet's directional message $$ m^{(l+1)}_{ji} \;=\; f_{\text{update}}\!\Big( m^{(l)}_{ji},\; \sum_{k\,\in\,N_j\setminus\{i\}} f_{\text{int}}\big( m^{(l)}_{kj},\, e^{(ji)}_{\text{RBF}},\, a^{(kj,ji)}_{\text{SBF}} \big) \Big) $$

Decode DimeNet (less scary than it looks)

$m_{ji}$ — a message from atom $j$ to atom $i$. DimeNet sends messages along atom pairs, not just atoms.
$f_{\text{update}}$, $f_{\text{int}}$ — two small learned mini-networks ("update" and "interact"). Adjustable blenders; don't fret their insides.
$\sum_{k \in N_j\setminus\{i\}}$ — "sum over every neighbour $k$ of atom $j$, except $i$ itself." These extra neighbours are what create the angles.
$e^{(ji)}_{\text{RBF}}$ — the distance from $j$ to $i$, in the RBF bell-bank form you just met.
$a^{(kj,ji)}_{\text{SBF}}$ — the angle made by the $k\!-\!j\!-\!i$ trio, encoded with "spherical" basis functions (SBF). This is the new ingredient SchNet lacked.

After all messages have flowed, DimeNet gathers them into a final description of each atom (equation 14):

Equation 14 · gather the messages into an atom $$ h_i \;=\; \sum_{j\,\in\,N_i} m_{ji} $$

Plainly: "atom $i$'s final description $h_i$ is the sum of all messages $m_{ji}$ its neighbours sent it." From $h_i$ the network reads the energy. The paper reports DimeNet (and its faster sequel DimeNet++) clearly beats SchNet — angles really do carry information distances miss.

Why angles matter. Both have the same two bond lengths from the centre atom — only the angle differs, making them different molecules with different energies. SchNet (distances only) could confuse them; DimeNet can't.

Models 3 & 4 — SchNOrb and PhiSNet: predicting the whole matrix

SchNet and DimeNet predict a single number (the energy). But the real grand prize from Chapter 6 was the entire Hamiltonian matrix. Two models chase it. SchNOrb extends SchNet to output the matrix blocks. PhiSNet goes further and bakes in the rotation rule from the invariance/equivariance box above.

Recall the subtlety: spin a molecule and its energy is unchanged (invariant — easy), but its Hamiltonian matrix must co-rotate in a precise, mathematically-required way (equivariant — hard). A model ignorant of this could give different matrices for the same molecule at different tilts — plainly wrong. PhiSNet builds the rule in from the ground up by storing its numbers in a special rotation-aware form, which the paper describes by its shape (equation 15):

Equation 15 · PhiSNet's rotation-aware features $$ \mathbb{R}^{\,F\times (L+1)^2} $$

Spherical harmonics, and where (L+1)² comes from

The clever $(L+1)^2$ deserves unpacking, and it ties straight back to the orbital shapes from Chapter 2. The features are organised by spherical harmonics — the natural set of "wave shapes on the surface of a sphere," labelled by a level $\ell = 0, 1, 2, \dots$. You've already seen the first few:

Counting the spherical harmonics

Level $\ell=0$: the sphere (the s shape from Chapter 2). There is $1$ of these.
Level $\ell=1$: the dumbbells (the p shapes). There are $3$ of these — one for each direction $x, y, z$.
Level $\ell=2$: the cloverleaves (the d shapes). There are $5$ of these.
The pattern: level $\ell$ has $2\ell+1$ shapes. Add them all up to level $L$: $1 + 3 + 5 + \dots + (2L+1) = (L+1)^2$. That's exactly the count in the formula. (The sum of the first $L+1$ odd numbers is a perfect square — a fact the ancient Greeks knew.)

So what is $\mathbb{R}^{\,F\times (L+1)^2}$?

$\mathbb R^{\dots}$ just means "a block of real numbers of this size" — it describes the container the model stores its features in.
$F$ — the number of feature channels: how many separate descriptive numbers the model keeps per atom (like having many sliders per atom).
$(L+1)^2$ — the spherical-harmonic slots we just counted. And here's the magic: spherical harmonics transform in a known, exact way when you rotate them. By storing features in these slots, PhiSNet guarantees its predicted matrix co-rotates correctly with the molecule. The physics is built into the plumbing, not hoped for. The paper calls this being SE(3)-equivariant.

The pattern across all four models Each model bakes in more real physics than the last. SchNet respects rotation for energies. DimeNet adds angles. PhiSNet guarantees the whole matrix rotates correctly. The more physical truth you build into the machine's bones, the less it must learn from scratch — and the better it does. That's the quiet thesis threaded through the entire paper, and it's the lesson the field keeps relearning: the best machine-learning is the kind that knows some physics before it ever sees the data.

CHAPTER EIGHT What they actually found

The verdict

So nablaDFT built a colossal library and trained these machines on it. What did they learn? The honest answer has three parts — and some of them did not go the way you'd hope. That candour is exactly what makes it good science.

First, how do you even score a guess? Meet the MAE

Before any result makes sense, you need a yardstick for "how wrong." The paper uses the Mean Absolute Error (MAE): for each test molecule, take the size of the gap between the guess and the truth, then average over all of them.

Mean Absolute Error · the scorecard $$ \text{MAE} \;=\; \frac{1}{N}\sum_{i=1}^{N}\big|\,E^{\text{pred}}_i - E^{\text{true}}_i\,\big| $$

Worked example · computing an MAE

Say a model predicts energies for 3 molecules. Truth: $-100,\ -200,\ -150$. Guesses: $-98,\ -205,\ -148$.

The absolute errors (sizes of the gaps, signs dropped): $|{-98}-({-100})| = 2$, then $|{-205}-({-200})| = 5$, then $|{-148}-({-150})| = 2$.

Average them: $(2 + 5 + 2)/3 = 3$. So the MAE is $3$ energy units. Smaller MAE = better model. Every number in the paper's results tables is one of these averages — now you can read them.

What counts as "good"? The chemical-accuracy bar

Chemists have a famous target called chemical accuracy: get within about $1$ kilocalorie per mole of the truth, and your prediction is trustworthy enough to do real chemistry with. In the energy units these models use (Hartree, written $E_{\text h}$), that bar is roughly:

Chemical accuracy, converted $$ 1\ \tfrac{\text{kcal}}{\text{mol}} \;\approx\; 0.0016\ E_{\text h} \;=\; 1.6\times 10^{-3}\ E_{\text h} $$

So the game is to push the MAE below roughly $1.6\times10^{-3}$ Hartree. Keep that number in mind as the finish line while we look at how the models actually did.

Second finding: the library is enormous — and more data helps

The dataset holds 1,004,918 molecules and 5,340,152 conformations (a conformation being one frozen 3-D pose — the Born–Oppenheimer snapshots from Chapter 3). Every one was computed with full DFT at the ωB97X-D/def2-SVP level using the Psi4 software. The whole thing fills around 100 terabytes.

A bigger library to learn from. And crucially, nablaDFT stores the full Hamiltonian matrix — the grand-prize object from Chapter 6 — for every conformation, which the others mostly don't.

When they trained the matrix-predictor PhiSNet on larger and larger slices, its error fell steadily — on one test, its Hamiltonian MAE dropped from 7.4 down to 2.9 (in units of $10^{-3}$ Hartree) as the training set grew. The hopeful reading: these models haven't hit their ceiling — feed them more and they keep improving.

One telling exception proved the point. The plain linear-regression baseline — the simplest possible model — barely improved at all with more data (it sat near 4.6, in units of $10^{-2}$ Hartree, no matter what). That's the signature of a model too simple to use what the extra data offers. More food only helps if you're able to digest it.

Third finding, the sobering one: this is much harder than it looked

Here's the result the paper is most candid about. On the famous older benchmarks (called QM9 and MD17), these models are spectacular — DimeNet++ reaches an energy MAE around 0.00023 Hartree on QM9, comfortably past the chemical-accuracy finish line. But on nablaDFT's far more varied, drug-like molecules, the same kind of model managed only about 3.2 ×10⁻² Hartree on the hardest test — more than a hundred times worse, and nowhere near the bar.

Why the collapse? Because the older benchmarks ask an easier question. Understanding which question is the key to the whole paper — and it's best seen as the difference between kinds of exams.

The four exams — how the paper tests generalization The paper splits its data several ways to test how well a model truly generalizes versus merely memorizes:

Structure split — train and test on different conformations (poses) of the same molecules. Like being tested on problems you practised, with the numbers changed. Easiest.
Scaffold split — test on molecules with different chemical backbones than any seen in training. Like being tested on a topic you never studied. Much harder, and far more honest about real-world use, where you'll meet molecules unlike your training set.
Conformation split — focuses on getting the energy ranking of a single molecule's poses right.

The older benchmarks are mostly the "easy exam." nablaDFT insists on the hard ones — and on the hard exams, today's celebrated models stumble. The victory laps were run on too gentle a track.

The headline result The paper's central message is a caution: today's best models look brilliant when tested on molecules similar to their training set, but their accuracy falls apart when asked to handle a genuinely diverse range of chemistry, or to generalize to unfamiliar molecular scaffolds. The field had been grading itself on too easy a curve.

Two smaller findings sharpen the picture. Some deep models, trained on only small slices of data, did worse than the dumb linear baseline — flexibility is useless without enough examples to pin it down. And one matrix-model, SchNOrb, wouldn't even train successfully on the larger sets, a reminder that these machines remain finicky and far from solved.

Why this is a gift, not a disappointment

It would be easy to read all this as failure. It is the opposite. The paper's real contribution is an honest, hard yardstick — plus the giant, carefully-computed dataset needed to attack it. Earlier tests flattered the models; nablaDFT shows the field exactly where the real work remains. (The 100-terabyte size is itself a consequence of that honesty: a high, expensive hybrid functional and a real basis set for five million poses is what rigor costs.)

A good benchmark isn't there to crown winners. It's there to show you, honestly, how far you still have to go.

The whole journey, in one breath Chemistry is electrons obeying Coulomb's law (Ch. 1). Electrons are waves described by $\Psi$, and $|\Psi|^2$ is where they're found (Ch. 2). The Schrödinger equation $\hat H\Psi=E\Psi$ decides which waves are allowed — but its electron-tangle term makes it unsolvable, and Born–Oppenheimer freezes the nuclei to help (Ch. 3). Those waves obey hard energy laws — conservation, and above all the variational principle that says any guess sits above the truth, which turns "solve the equation" into "minimize the energy" (Ch. 3½). So we turn waves into a grid of numbers via a basis, and minimizing energy gives $\mathbf F\mathbf c=\varepsilon\mathbf S\mathbf c$, cranked by the self-consistent loop (Ch. 4). DFT's density shortcut — licensed by a clean contradiction proof of Hohenberg–Kohn, made practical by the Kohn–Sham equations, with all the difficulty quarantined in an exchange–correlation term we could even partly derive — gets the energy from the simple cloud (Ch. 5). It all lands in the block Hamiltonian matrix, whose eigenvalues are orbital energies and whose HOMO–LUMO gap sets reactivity and colour (Ch. 6) — a matrix that is symmetric, built from angular-momentum shells, and sparse with exponential decay, which is the very reason a local network can predict it (Ch. 6½). That matrix is slow to compute — $N^3$ — so we train neural networks that grow smarter as they bake in more physics, from SchNet's distances to DimeNet's angles to PhiSNet's rotation-aware spherical harmonics (Ch. 7). And the honest verdict: the guessing genuinely works, more data helps, but real chemistry is far harder than the easy benchmarks let on (Ch. 8).

CHAPTER NINE The sequel · ∇²DFT (2024)

What happened next: nabla-squared

Two years after the first paper, the same team came back with a sequel — and gave it a sequel's name. If the original was nablaDFT ($\nabla$DFT), this one is ∇²DFT, "nabla-squared DFT." The little square is a chemist's in-joke: you met $\nabla^2$ back in Chapter 3 as the operator for an electron's motion energy. Here it just means "more" — a bigger, richer, more honest second pass at the same dream. Everything you learned still holds; this chapter is what the team added once they had the first benchmark to learn from.

The same problem, twice as much of it

The first move was simply scale. ∇²DFT roughly doubles the molecules and triples the conformations of the original. The numbers, straight from the paper:

Bigger in every direction. 1,936,929 molecules and 12,676,264 conformations — and once the new relaxation trajectories are counted in, 15,716,667 conformations in all, weighing about 220 terabytes. Crucially, the chemistry is unchanged: this is the same problem from Chapters 1–8, just with far more examples to learn from.

From one answer to a whole movie: relaxation trajectories

The headline new idea isn't size — it's a new kind of data, and it solves a problem we quietly skipped. All through this book, we computed the energy of a molecule frozen in one pose (a "conformation," from Chapter 3). But molecules are floppy, and a randomly-generated pose is usually an awkward, high-energy one. What chemists actually want is the relaxed pose — the molecule settled into its most comfortable, lowest-energy shape, like a dropped string slithering into a heap.

Finding that comfortable shape is called geometry optimization (or relaxation), and it works exactly like the gradient descent from Chapter 7 — but now you're rolling the physical molecule downhill in energy instead of rolling a network's weights downhill in error. At each step you compute the forces (remember, $\vec F = -\nabla E$), nudge every atom a little in the direction it's being pushed, and repeat until nothing wants to move. The sequence of poses along the way — from awkward start to relaxed finish — is a trajectory.

A trajectory is the whole downhill walk. ∇²DFT is the first dataset to provide these relaxation movies for a large number of drug-like molecules — about 60,000 optimizations across roughly 17,000 molecules, adding up to around 3 million individual geometries. Why bother storing the whole walk instead of just the destination? Because that's what teaches a neural network to do the walking.

Why trajectories are the point Here's the real prize. Every DFT relaxation needs dozens of expensive energy-and-force calculations — one at each step of the walk. If a neural network could predict those forces instead, it could relax a molecule in seconds rather than CPU-hours. But to train such a network, you need examples of the walk itself, not just its endpoint. ∇²DFT provides the walks — making it, in the paper's word, a universal dataset: not just "what is this molecule's energy?" but "how does this molecule find its shape?"

A richer haul from every calculation

The first dataset mainly stored the Hamiltonian and overlap matrices. ∇²DFT keeps those but saves far more from each DFT run — 17 molecular properties plus the complete wavefunction object (the full output of the Psi4 solver). That last part matters: with the entire wavefunction in hand, you can compute almost any property you want later, without rerunning anything. Several of these quantities are old friends from earlier chapters:

A few of the 17 properties — you already know them

Total energy, formation energy — the molecule's energy budget (Chapters 3 & 5), including the formation energy, which measures how much energy was released building it from loose atoms.
Forces — $\vec F = -\nabla E$ on every atom (Chapter 7), the thing that drives relaxation.
HOMO, LUMO, and the HOMO–LUMO gap — the reactivity-and-colour numbers you computed by hand and played with on the ladder in Chapter 6.
Dipole moment — how lopsided the molecule's charge is, which sets how it dissolves and docks.
The exchange–correlation, one-electron, two-electron, and nuclear-repulsion energies — the individual lumps of equation 9 from Chapter 5, now reported separately.

Three tasks, and a bigger stable of machines

The original paper tested two things: predicting the Hamiltonian, and predicting energy. ∇²DFT builds a proper benchmark with three tasks and a tidy, extendable framework holding ten models — including several that postdate the first paper. The three tasks:

The new third task is the star. Tasks ① and ② extend the first paper (though ① now predicts the full Fock matrix, not just the easier "core" part). Task ③ — using a network to relax a molecule into its lowest-energy shape — is entirely new, and it's what the trajectories were collected to support.

The ten models, sorted by what they do

Hamiltonian predictors: SchNOrb and PhiSNet (both from Chapter 7) plus a newcomer, QHNet, another equivariant matrix-predictor.
Neural network potentials (energy + forces): the old guard — linear regression, SchNet, DimeNet++ — joined by a wave of newer, more powerful architectures: PaiNN, GemNet-OC, Graphormer3D, EquiformerV2, and eSCN. Several of these are equivariant in the deep sense from Chapter 7's spherical-harmonics box.
And to study how much data each model actually needs, every model is trained on a ladder of dataset sizes — tiny, small, medium, large — across 12 predefined splits of the structure / scaffold / conformation kind you met in Chapter 8's "four exams."

What they found — the honest verdict, sharpened

The results rhyme with the first paper's, and deepen its central warning.

The three big findings

More data reliably helps. Across every task, models improved as the training set grew — underscoring exactly why building a 15-million-conformation dataset was worth the effort.
The generalization gap is real and large. On the matrix task, the best model (PhiSNet) scored an MAE around $0.2$–$0.4\times10^{-3}$ Hartree on ∇²DFT — yet the same kind of model reaches $1.8\times10^{-5}$ on the easier, less-diverse MD17 set. A diverse, drug-like test set is still dramatically harder than the gentle benchmarks the field once celebrated — the very point the first paper made, now confirmed at larger scale.
Neural networks can relax molecules — and beat the classics. On the new optimization task, NNPs clearly outperformed traditional fast methods (like the MMFF force field and semi-empirical xTB), and got better still when fine-tuned on the relaxation trajectories. A genuine, useful win.

There were honest stumbles too, faithfully reported: PhiSNet, the best matrix-predictor, failed to finish training on the largest data splits even after 1,920 hours of GPU time — a reminder that these machines remain finicky and costly. And the surprise from the first paper held: modelling the full electronic structure (as SchNOrb does) helped it predict plain energies about as well as the specialised energy models, hinting that teaching a network real quantum structure makes it better at everything.

Same dream, sharper tools, and the same honest answer: getting closer, not yet there.

Where the two papers leave us The first paper (Chapters 1–8) built the ladder: from a single electron's wave all the way up to neural networks guessing the Hamiltonian, and a hard benchmark proving how far there was to go. The sequel, ∇²DFT, doubled the data, added the wavefunction and the all-important relaxation trajectories, broadened the contest to ten models and three tasks — and showed that machines can now not only describe a molecule but help it find its shape. The chase that has run for a hundred years is still on. But the tools are getting very, very good — and the honest scorecards that papers like these provide are exactly how the field knows which way is downhill.

The end · for now

Find the
wave.

You just climbed from "what is an atom" to fifteen equations out of a real research paper — through Coulomb's law, the wave function, the Schrödinger equation, the basis-set translation, the self-consistent field, density functional theory, the block Hamiltonian, gradient descent, and four neural networks built to guess it all — and then on to its sequel, with relaxation trajectories and a stable of ten models. You met the same chase that has occupied physicists for a hundred years, and you didn't skip the maths. Not bad for one little book.

Based on nablaDFT: Large-Scale Conformational Energy and Hamiltonian Prediction benchmark and dataset, Khrabrov et al., Phys. Chem. Chem. Phys., 2022, 24, 25853–25863; and its sequel ∇²DFT: A Universal Quantum Chemistry Dataset of Drug-Like Molecules and a Benchmark for Neural Network Potentials, Khrabrov et al., NeurIPS 2024 (arXiv:2406.14347). Equations and figures (dataset sizes, error values, the level of theory) are drawn from those papers. All illustrations original.

How to read this book

Atoms, and the rule of attraction

First, how unbelievably empty an atom is

The force that holds it all together

The periodic table is just an electron filing cabinet

An electron is a wave

Why would anyone believe such a thing?

Meet psi: the wave function

The wave has to add up to one electron

Try it: an electron in a box

In a real atom, the waves are called orbitals

One more property: spin

Operators: machines that ask questions of a wave