Density functional theory & Kohn–Sham

Hartree–Fock gave us a mean-field wavefunction but missed correlation. Recovering it by improving the wavefunction (post-HF methods) is accurate but expensive. Density functional theory (DFT) takes a completely different, and cheaper, route: work not with the \(3N\)-dimensional wavefunction but with the electron density \(\rho(\mathbf r)\) — a function of just three variables. This is the workhorse of modern quantum chemistry, and it is what mychk.scf(ref="r", xc="...").run() computes.

The density is enough: Hohenberg–Kohn

It is not obvious that the simple density \(\rho(\mathbf r)\) contains as much information as the full wavefunction. The Hohenberg–Kohn theorems say it does.

Theorem 2 (Hohenberg–Kohn)

1. The ground-state density determines everything. For a system of electrons in an external potential \(v_{\text{ext}}\) (the nuclei), \(v_{\text{ext}}\) is fixed — up to a constant — by the ground-state density \(\rho(\mathbf r)\). Hence every ground-state property, including the energy, is a functional of the density: \(E = E[\rho]\).

2. A variational principle for the density. There is a universal functional such that \(E[\rho] \ge E_0\) for any trial density, with equality at the true ground-state density \(\rho_0\).

In principle this is miraculous: the impossible \(3N\)-dimensional problem collapses to finding one 3-D function. The catch is equally severe — the exact functional \(E[\rho]\) is unknown, and its worst piece is the kinetic energy, which is very hard to write in terms of \(\rho\) alone.

The Kohn–Sham trick

Kohn and Sham found the escape that makes DFT practical. Reintroduce orbitals, but only to handle the kinetic energy: imagine a fictitious system of non-interacting electrons that has the same density as the real one. Its kinetic energy \(T_s[\rho]\) is easy (it is just the orbital kinetic energy, as in HF), and we lump everything we do not know into one term. The energy splits as

\[ E[\rho] = \underbrace{T_s[\rho]}_{\text{non-interacting kinetic}} + \underbrace{\int v_{\text{ext}}(\mathbf r)\,\rho(\mathbf r)\,d\mathbf r}_{\text{nuclear attraction}} + \underbrace{J[\rho]}_{\text{Hartree (classical Coulomb)}} + \underbrace{E_{\text{xc}}[\rho]}_{\text{exchange–correlation}} , \]

where the exchange–correlation functional \(E_{\text{xc}}[\rho]\) absorbs everything left over — exchange, correlation, and the small kinetic-energy difference between the real and non-interacting systems. It is the only unknown.

Minimizing this energy gives the Kohn–Sham equations, one-electron eigenvalue equations that look exactly like Hartree–Fock:

\[ \left[-\tfrac12\nabla^2 + v_{\text{eff}}(\mathbf r)\right]\psi_i = \varepsilon_i\,\psi_i, \qquad v_{\text{eff}} = v_{\text{ext}} + v_{\text{H}} + v_{\text{xc}}, \quad v_{\text{xc}} = \frac{\delta E_{\text{xc}}}{\delta \rho}, \]

with the density rebuilt from the occupied Kohn–Sham orbitals, \(\rho(\mathbf r) = \sum_i |\psi_i(\mathbf r)|^2\). Because \(v_{\text{eff}}\) depends on \(\rho\), these are solved by the same self-consistent-field loop as Hartree–Fock — the same run(log=...) cycle table, the same DIIS. In fact, if you set \(E_{\text{xc}}\) to the exact-exchange-only expression, you recover Hartree–Fock. DFT and HF are two points on one continuum.

Choosing a functional: Jacob’s ladder

Everything hinges on the approximation to \(E_{\text{xc}}\). These are ranked on a famous “Jacob’s ladder” of increasing sophistication — the choices behind xc=:

• LDA (local density approximation) — \(E_{\text{xc}}\) depends on \(\rho\) at each point only. The simplest rung.

• GGA (generalized gradient approximation, e.g. pbe) — also uses the density gradient \(\nabla\rho\).

• meta-GGA (e.g. scan) — adds the kinetic-energy density \(\tau\) (and sometimes \(\nabla^2\rho\)).

• Hybrids (e.g. b3lyp, pbe0) — mix in a fraction of exact (Hartree–Fock) exchange. This rung literally blends the HF of the last chapter with DFT, and is the most popular choice for molecules.

qc-rs passes the functional name straight to libxc, so any of its LDA/GGA/meta-GGA/hybrid functionals works in xc=.

The numbers, and an important caveat

Here is water again (RHF vs three functionals, all cc-pVDZ):

method	xc=	rung	energy / \(E_h\)
Hartree–Fock	—	—	−76.026772
LDA	lda	local	−75.854689
PBE	pbe	GGA	−76.333442
B3LYP	b3lyp	hybrid	−76.420369

Notice something that trips up beginners: LDA’s energy is higher than Hartree–Fock’s. That is not a bug — it is a fundamental difference from the variational picture of chapters 2–3. HF gives a rigorous upper bound to the exact energy (it uses the exact Hamiltonian), but a DFT total energy uses an approximate \(E_{\text{xc}}\), so it is not a bound at all. You therefore cannot rank functionals by their total energy — a “better” functional need not give a lower number. Functionals are judged by how well they reproduce reference data (geometries, reaction energies, experiment), not by the energy itself. (B3LYP’s lower value here does reflect the correlation it captures, but that comparison is not a rigorous ranking.)

Strengths, weaknesses, and where to go

DFT’s appeal is its cost–accuracy trade-off: for roughly the price of Hartree–Fock, a good functional folds in much of the correlation HF misses — which is why it dominates practical chemistry. Its limits follow from the same source:

• No systematic route to the exact answer. Unlike wavefunction methods (HF → MP2 → coupled cluster → …), there is no ladder you can climb to convergence; you rely on functionals validated against data.

• Known failure modes. Most functionals miss long-range dispersion (fixed by the DFT-D3/D4 corrections in Solvation & dispersion), and many suffer self-interaction error.

You have now met the two pillars of practical electronic structure: Hartree–Fock and DFT. With the theory in hand, the User guide shows how to run them well — choosing references, functionals, guesses, and convergence — starting from molecular input.