3 Density functional theory

3.1 Why first principles?

Classical potentials have fundamental limitations: the accuracy is bounded by the functional form you choose, and someone has to develop and validate a potential for your specific system.

First-principles methods address this. Instead of parameterising interactions, we calculate \(E(r)\) from quantum mechanics—solving (approximately) for the electrons in your system. The physics is built in, not fitted.

3.2 From wavefunctions to density

In principle, the energy of a collection of electrons and nuclei comes from solving the Schrödinger equation. The ground-state wavefunction \(\psi(r_1, r_2, \ldots, r_n)\) gives everything — the energy, the electron density, all ground-state properties. The problem is dimensionality. For \(N\) electrons, the wavefunction is a function of \(3N\) spatial coordinates. The complexity grows exponentially with the number of electrons — intractable for real materials.

The electron density \(\rho(\mathbf{r})\) offers an alternative. Unlike the wavefunction, the density is always three-dimensional — ten electrons or ten thousand, it is still just \(\rho(x, y, z)\). The first Hohenberg-Kohn theorem (1964) established that for ground states, the energy can be written as a functional of the density: \(E = E[\rho]\). A functional maps a function to a number — here \(\rho(\mathbf{r})\) is a function of position, and \(E\) is a single number. The theorem does not specify what this functional is, only that it exists. The second Hohenberg-Kohn theorem establishes that this functional is variational: the density that minimises \(E[\rho]\) is the true ground-state density.

Together, these theorems imply that the exact ground-state energy could be found by minimising over three-dimensional densities rather than \(3N\)-dimensional wavefunctions — if we knew \(E[\rho]\). We can write:

\[E[\rho] = T[\rho] + V_\mathrm{ee}[\rho] + \int v(\mathbf{r})\rho(\mathbf{r})\,d\mathbf{r}\]

where \(T[\rho]\) is the kinetic energy, \(V_\mathrm{ee}[\rho]\) is the electron-electron repulsion, and the final term is the electron-nuclear attraction. The last term is straightforward — we know \(v\) from the nuclear positions. The first two are not: we do not know how to write them as explicit functionals of the density.

3.3 The Kohn-Sham approach

Kohn and Sham (1965) addressed this by separating each unknown term into a part that can be calculated plus a remainder.

For the electron-electron interaction, the classical Coulomb energy of a charge distribution is:

\[J[\rho] = \frac{1}{2}\int\int \frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}\, \mathrm{d}\mathbf{r}\, \mathrm{d}\mathbf{r}'\]

This can be evaluated directly from the density. What it misses are exchange (electrons with parallel spin avoid each other, a consequence of the Pauli principle) and correlation (electrons avoid each other due to their mutual repulsion).

For the kinetic energy, there is a system for which this is straightforward: non-interacting electrons. For such a system, the ground state is a Slater determinant of one-electron orbitals \(\phi_i\), and the kinetic energy is the sum of one-electron contributions:

\[T_s = -\frac{1}{2}\sum_i \langle\phi_i|\nabla^2|\phi_i\rangle\]

The true kinetic energy \(T[\rho]\) differs from \(T_s[\rho]\) because interacting electrons move differently from non-interacting ones — but the difference is small. Both systems have the same density, so electrons are confined to the same regions of space, and kinetic energy is largely determined by spatial confinement.

The remainders from both separations are bundled into a single correction, the exchange-correlation functional:

\[E_\mathrm{xc}[\rho] = (T[\rho] - T_s[\rho]) + (V_\mathrm{ee}[\rho] - J[\rho])\]

The total energy becomes:

\[E[\rho] = T_s[\rho] + J[\rho] + \int v(\mathbf{r})\rho(\mathbf{r})\,d\mathbf{r} + E_\mathrm{xc}[\rho]\]

The first three terms can be calculated; all the unknowns are in \(E_\mathrm{xc}\).

To evaluate \(T_s\) requires finding the orbitals of the non-interacting system. These satisfy one-electron Schrödinger equations in an effective potential \(v_\mathrm{eff} = v + v_J + v_\mathrm{xc}\), where \(v\) is the nuclear potential, \(v_J\) is the classical Coulomb potential from the electron density, and \(v_\mathrm{xc}\) is the functional derivative of \(E_\mathrm{xc}\) with respect to \(\rho\). The density is constructed from the occupied orbitals:

\[\rho(\mathbf{r}) = \sum_i |\phi_i(\mathbf{r})|^2\]

Since \(v_J\) and \(v_\mathrm{xc}\) depend on \(\rho\), which depends on the orbitals, which depend on \(v_\mathrm{eff}\), the equations must be solved iteratively until self-consistent.

Nothing in this development is an approximation — it is exact repackaging into terms that can be evaluated plus a correction. If \(E_\mathrm{xc}[\rho]\) were known exactly, Kohn-Sham DFT would give the exact ground-state energy. All the approximation in practical DFT lives in one place: \(E_\mathrm{xc}\).

3.4 Exchange-correlation functionals

The simplest approximation is the local density approximation (LDA): at each point, use the exchange-correlation energy of a uniform electron gas with that density. This is surprisingly effective, but tends to overbind—predicting bonds that are too strong and too short.

Generalised gradient approximations (GGAs) improve on this by including how the density varies—not just \(\rho\) but also \(\nabla \rho\). PBE is the most common GGA in materials science.

Meta-GGAs go further, including the kinetic energy density as well. r²SCAN is increasingly the reference level for materials science calculations.

Hybrid functionals mix in a fraction of exact exchange from Hartree-Fock theory. These often give better band gaps and reaction energies, but are significantly more expensive—you’re reintroducing the costly exact exchange integrals.

3.5 Strengths and limitations

DFT is first-principles and transferable—the same functional works across different chemistries without system-specific parameterisation. It’s accurate for structures and relative energies, and scales as roughly N³, making it practical for systems of hundreds of atoms.

The limitations are significant. The computational cost means DFT is practical for hundreds of atoms, not thousands, and for static calculations or short dynamics rather than the long timescales often needed to study diffusion or phase transitions. For those problems, you need classical potentials or machine-learned interatomic potentials.

More fundamentally, DFT offers no systematic path to the exact answer. In wavefunction methods, you can improve accuracy by including more excitations or larger basis sets—the path is clear even if expensive. In DFT, accuracy depends on choosing a good functional, and a better functional requires new physical insight, not just more computer time. This also means different functionals give different answers, and it’s not always obvious which to trust.

Finally, approximate functionals suffer from self-interaction error: an electron spuriously repels itself because \(J[\rho]\) includes the interaction of each electron with the total density, including its own contribution. This causes particular problems for localised states and transition-metal chemistry, where the error doesn’t cancel out.