Density Functional Theory

Learning Objectives

Section 4.4 left us face-to-face with the most embarrassing fact in electronic-structure theory: the many-body wavefunction $\Psi(\mathbf{r}_1, \mathbf{r}_2, \ldots, \mathbf{r}_N)$ is impossible to store, let alone solve, for any system you actually care about. Even a humble water molecule needs a function of 30 coordinates; a CdSe quantum dot of 100 atoms involves hundreds. The universe simply does not have enough atoms to write $\Psi$ down on a grid.

And yet VASP routinely simulates such systems on a laptop. The reason is a mathematical sleight of hand so audacious that two of its inventors won a Nobel Prize for it: throw away $\Psi$ and keep only the electron density $n(\mathbf{r})$ — a humble function of three variables. This section explains why that move is allowed, what it costs, and what it does not give you.

By the end of this section you should be able to:

1. State the Hohenberg–Kohn theorems in your own words and explain why theorem 1 is so surprising.
2. Write the total energy as a functional $E[n] = T[n] + V_{\text{ext}}[n] + V_{ee}[n]$ and identify which piece is universal and which is system-specific.
3. Apply the variational principle (HK2) to a one-parameter trial density and find the minimum by hand.
4. Implement the original Thomas–Fermi DFT for the H atom in NumPy and quantify how badly it does compared with the exact answer.
5. Articulate, precisely, what is unknown about $F[n]$ and why that motivates the Kohn–Sham construction in Section 4.6.
6. Locate the electron density in a real VASP run — CHGCAR, CHG, partial densities — and read off basic properties from it.

Where Section 4.4 Left Us Stranded

For a non-relativistic crystal, the electronic Schrödinger equation is, in atomic units,

$\hat{H}\,\Psi = \left[ -\tfrac{1}{2}\sum_{i} \nabla_i^2 \;+\; \sum_{i} v_{\text{ext}}(\mathbf{r}_i) \;+\; \tfrac{1}{2}\sum_{i \neq j} \frac{1}{|\mathbf{r}_i - \mathbf{r}_j|} \right] \Psi = E\,\Psi$

Three terms: kinetic energy of electrons, attraction to nuclei, and electron–electron repulsion. Innocent enough on paper. The catastrophe is the third term: it couples every pair of electrons, which means $\Psi$ cannot be written as a product of one-electron functions. Storing$\Psi$ on a 100-point grid for a 50-electron system needs $100^{150} \approx 10^{300}$ numbers. The observable universe contains about $10^{80}$ protons.

The Wild Idea: Forget Ψ, Keep n(r)

Notice something subtle about the Hamiltonian. The kinetic operator and the electron–electron repulsion look the same for every system: a chunk of copper, a benzene ring, a Mn:CdSe quantum dot, the inside of a neutron star. Only the external potential $v_{\text{ext}}(\mathbf{r})$ — the Coulomb attraction of this particular set of nuclei at these particular positions — tells one material from another.

So the question is: do we really need the full $\Psi$, or is some smaller object enough to reconstruct everything we want?

The boldest possible answer was given by Pierre Hohenberg and Walter Kohn in 1964: the electron density is enough. Not the wavefunction. Not the density matrix. Just the scalar field $n(\mathbf{r})$ — how many electrons per cubic bohr, at every point in space.

Wavefunction, Density, and Kohn-Sham Orbitals

DFT becomes less mysterious once the three central objects are kept separate. The exact many-body wavefunction is the formal object, the density is the physical object Hohenberg and Kohn elevate, and the Kohn–Sham orbitals are the computational device VASP actually manipulates.

What Exactly Is the Electron Density?

Given any (normalised) many-body wavefunction $\Psi(\mathbf{r}_1, \ldots, \mathbf{r}_N)$, the electron density is

$n(\mathbf{r}) = N \int |\Psi(\mathbf{r}, \mathbf{r}_2, \ldots, \mathbf{r}_N)|^2 \, d\mathbf{r}_2 \cdots d\mathbf{r}_N$

Operationally: pin one electron at $\mathbf{r}$, integrate out the other $N-1$ coordinates, and multiply by $N$ for the indistinguishability of identical particles. The result satisfies two simple sum rules every working DFT user invokes daily:

• Normalisation: $\int n(\mathbf{r}) \, d\mathbf{r} = N$ (total electron count).
• Cusps at nuclei: $\left.\frac{\partial \bar{n}}{\partial r}\right|_{r \to 0} = -2 Z_\alpha\,\bar{n}(0)$ (Kato's cusp condition; the density tells you both where the nuclei are and what charge they carry).

That second bullet is already a hint: the density is not just a blob of charge. It encodes precise structural information.

Hohenberg–Kohn Theorem 1: Density Determines Everything

HK1 (1964): For a non-degenerate ground state, the external potential $v_{\text{ext}}(\mathbf{r})$ is uniquely determined — up to an irrelevant additive constant — by the ground-state density $n_0(\mathbf{r})$.

Read that twice. The density tells you the potential. Knowing the potential, you know the entire Hamiltonian (the kinetic and repulsion pieces are universal). Knowing the Hamiltonian, you can in principle solve for the wavefunction and any other observable. Therefore, in a precise mathematical sense,

$\boxed{\quad n_0(\mathbf{r}) \;\Longleftrightarrow\; v_{\text{ext}}(\mathbf{r}) \;\Longleftrightarrow\; \hat{H} \;\Longleftrightarrow\; \Psi \;\Longleftrightarrow\; \text{everything}.\quad}$

The proof is a one-page reductio ad absurdum: assume two different potentials $v_1, v_2$ share the same ground-state density, apply the Rayleigh–Ritz variational principle to each, and add the two strict inequalities. The contradiction $E_0 < E_0$ falls out immediately.

Interactive: The Density ↔ Potential Map

Drag the red ions and watch the orange external potential and the cyan electron density shift together. They are two sides of the same coin: HK1 says the cyan curve uniquely determines the orange one. Add or remove ions and the bijection still holds — the density even encodes how many nuclei there are and what charges they carry, through cusps you can spot by eye.

Reading the density backwards: peaks pinpoint nuclei, peak heights scale with $Z$, and the area under the curve equals the total electron count. A trained DFT eye reads a CHGCAR slice the way a doctor reads a chest X-ray.

Hohenberg–Kohn Theorem 2: The Variational Principle

HK2 (1964): There exists a universal energy functional $E_v[n]$ such that, for any trial density $\tilde{n}(\mathbf{r})$ integrating to $N$,

$E_v[\tilde{n}] \;\geq\; E_v[n_0] \;=\; E_{\text{ground}},$

with equality if and only if $\tilde{n}(\mathbf{r}) = n_0(\mathbf{r})$. In other words, the ground-state energy is the global minimum of an energy functional — a function whose argument is itself a function. This is the analogue of the Rayleigh–Ritz variational principle, but lifted from the space of wavefunctions (size $\infty^{3N}$) to the space of densities (size $\infty^{3}$).

HK2 is what makes DFT computable. Pick any physically admissible trial density, plug it into $E_v[n]$, and you get an upper bound on the true ground-state energy, provided the trial density is physically admissible: non-negative, normalised to the correct $N$, and representable by at least one antisymmetric $N$-electron state. Improve the density — the bound improves. Find the minimum — you have the answer.

The Universal Functional F[n]

Following Hohenberg and Kohn, decompose the total energy as

$E_v[n] \;=\; F[n] \;+\; \int v_{\text{ext}}(\mathbf{r})\,n(\mathbf{r})\,d\mathbf{r}.$

The second term is system-specific but trivial: just the electrostatic energy of the density in the field of the (known) nuclei. The first term $F[n]$ contains the kinetic energy $T[n]$ and the electron–electron interaction $V_{ee}[n]$. Crucially, $F[n]$ is the same functional for every Coulombic system: copper, benzene, CdSe quantum dot, the corona of the sun. Once we know $F$, every electronic-structure problem in the universe becomes a one-line minimisation:

$\boxed{\quad E_{\text{ground}} = \min_{n} \left\{ F[n] + \int v_{\text{ext}}\,n\, d\mathbf{r}\right\}.\quad}$

The Catch: Nobody Knows F[n]

And here is the punchline. Hohenberg and Kohn proved $F[n]$ exists. They did not write it down. Sixty years and a Nobel Prize later, the exact form of $F[n]$ is still unknown. We have very good approximations — some of them excellent for solids — but not the truth.

The structure of the unknown is this. Split $F[n]$ into pieces we know and pieces we don't:

$F[n] \;=\; \underbrace{T_s[n]}_{\text{independent-particle kinetic}} \;+\; \underbrace{\tfrac{1}{2}\!\!\iint \frac{n(\mathbf{r})\,n(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}\,d\mathbf{r}\,d\mathbf{r}'}_{\text{Hartree, } J[n]} \;+\; \underbrace{E_{xc}[n]}_{\text{exchange-correlation}}.$

The Hartree term is an explicit density functional. $T_s$ is explicit once we introduce Kohn–Sham orbitals, but it is not known as a simple closed-form function of $n(\mathbf{r})$ alone. The third piece, $E_{xc}[n]$, contains both leftovers: the difference between the true kinetic energy and $T_s$, plus the difference between the true electron–electron interaction and the classical Hartree self-repulsion. Modelling $E_{xc}$ is the entire art of modern DFT.

First Concrete Functional: Thomas–Fermi

Long before HK proved their theorems, Llewellyn Thomas and Enrico Fermi (independently, 1927) wrote down the first explicit density functional. Their approximation:

• Borrow the kinetic-energy density of a uniform electron gas (Section 4.1), evaluated at the local density $n(\mathbf{r})$.
• Drop exchange and correlation entirely.
• Treat the Coulomb interaction classically (just $J[n]$).

The kinetic functional is

$T_{\text{TF}}[n] = C_F \int n(\mathbf{r})^{5/3}\, d\mathbf{r}, \qquad C_F = \tfrac{3}{10}\,(3\pi^2)^{2/3} \approx 2.871.$

For our Thomas–Fermi-style one-electron toy model of hydrogen, the Hartree term $J[n]$ is meaningless (there is no second electron to repel against), so the energy reduces to

$E[n] = T_{\text{TF}}[n] + V_{\text{ext}}[n] = C_F \int n^{5/3}\,d\mathbf{r} - Z\int \frac{n(\mathbf{r})}{r}\, d\mathbf{r}.$

Pick a one-parameter family of trial densities — say spherical Gaussians of width $\alpha$ —

$n_\alpha(r) = \frac{1}{(\pi\alpha^2)^{3/2}}\,e^{-r^2/\alpha^2},$

and compute the energy as a function of $\alpha$. The integrals close in elementary form:

$T_{\text{TF}}(\alpha) = \frac{C_F (3/5)^{3/2}}{\pi\alpha^2}, \qquad V_{\text{ext}}(\alpha) = -\frac{2 Z}{\sqrt{\pi}\,\alpha}.$

Setting $dE/d\alpha = 0$ for $Z = 1$ gives $\alpha^\star \approx 0.753\,a_0$ with $E^\star \approx -0.749\,\text{Ha}$. Compare the exact answer for a hydrogen atom: $E_{\text{exact}} = -0.500\,\text{Ha}$. This TF-style density-only model over-binds the electron by 50%. The functional is qualitatively right (a sane minimum exists, the density peaks at the nucleus, units are bohr and hartree) and quantitatively wrong (a critical 0.25 Ha off).

Interactive: Variational Principle in Action

Slide $\alpha$ below and watch the trial Gaussian on the left widen or narrow. The right panel decomposes the energy: a violet kinetic curve (penalises sharp densities), an orange external curve (rewards densities concentrated at the nucleus), and a cyan total. The yellow star marks the variational minimum — the “best” Thomas–Fermi answer for this trial family.

Why does the kinetic curve fall as 1/α²? Heisenberg. Confining an electron to a Gaussian of width $\alpha$ sets its momentum spread to $\Delta p \sim 1/\alpha$, so the kinetic energy $p^2/(2m)$ goes as $1/\alpha^2$. Squeezing the cloud costs kinetic energy.

Why does the external curve fall as −1/α? The Coulomb potential of a single point charge looks the same on every length scale up to a single factor; integrating it against a Gaussian of width $\alpha$ gives an answer proportional to $1/\alpha$. A more compact density sees a stronger nuclear attraction.

Code Walk-Through: A Working Density-Only Model in 30 Lines

Time to make the picture concrete. The script below evaluates the Thomas–Fermi-style total energy for a Gaussian trial density at 17 widths $\alpha$, picks the smallest one, and prints the comparison with the exact H atom. This is the variational skeleton of a DFT calculation — choose an allowed density, evaluate an energy functional, improve the density — but it is not yet a Kohn–Sham self-consistent field calculation. Click any line on the right to see what is in memory at that point.

1import numpy as np
2
3# Thomas-Fermi kinetic-energy constant (Hartree atomic units)
4C_F = (3.0 / 10.0) * (3.0 * np.pi**2)**(2.0 / 3.0)
5
6# Trial electron density: spherical Gaussian, total charge N around Z=1 nucleus
7def density(r, alpha, N=1.0):
8    """n(r) = (N / (pi*alpha^2)^(3/2)) * exp(-r^2 / alpha^2)."""
9    return (N / (np.pi * alpha**2)**1.5) * np.exp(-r**2 / alpha**2)
10
11# Radial integration helper: integrand g(r) over 0..r_max, 4*pi*r^2 included
12def radial_integral(g, alpha, n_pts=4000):
13    r = np.linspace(1e-10, 30 * alpha, n_pts)
14    return np.trapezoid(g(r) * 4.0 * np.pi * r**2, r)
15
16# DFT energy components for the trial density
17def kinetic_TF(alpha, N=1.0):
18    return C_F * radial_integral(lambda r: density(r, alpha, N)**(5.0/3.0), alpha)
19
20def external_energy(alpha, Z=1.0, N=1.0):
21    return -Z * radial_integral(lambda r: density(r, alpha, N) / r, alpha)
22
23def total_energy(alpha):
24    return kinetic_TF(alpha) + external_energy(alpha)
25
26# Variational scan: which alpha minimises E[n_alpha]?
27alphas = np.linspace(0.4, 2.0, 17)
28print("  alpha     T_TF       V_ext       E_total")
29for a in alphas:
30    print(f"  {a:5.3f}   {kinetic_TF(a):7.4f}   {external_energy(a):7.4f}   {total_energy(a):7.4f}")
31
32best = alphas[np.argmin([total_energy(a) for a in alphas])]
33print(f"\nBest alpha on grid: {best:.3f}, E = {total_energy(best):.4f} Ha")
34print(f"Exact H 1s energy : -0.5000 Ha   (TF-style model over-binds)")

Run the script and you will see:

1alpha     T_TF       V_ext       E_total
2  0.400    2.6548    -2.8209    -0.1661
3  0.500    1.6993    -2.2568    -0.5575
4  0.600    1.1799    -1.8806    -0.7007
5  0.700    0.8669    -1.6120    -0.7451
6  0.800    0.6637    -1.4105    -0.7468
7  0.900    0.5244    -1.2538    -0.7294
8  1.000    0.4248    -1.1284    -0.7036
9  ...
10  2.000    0.1062    -0.5642    -0.4580
11
12Best alpha on grid: 0.800, E = -0.7468 Ha
13Exact H 1s energy : -0.5000 Ha   (TF-style model over-binds)

Why Thomas–Fermi Fails (and Why That Matters)

Thomas–Fermi's 50% error on hydrogen is not a quirk; it is the symptom of two structural problems that Kohn and Sham would fix 12 years later.

VASP Connection: Where the Density Lives

VASP is, at its core, a machine for minimising $E_v[n]$ with a particular choice of $E_{xc}[n]$ over a plane-wave basis. The density is a first-class citizen of every output file.

1. CHGCAR — the converged charge density

After every self-consistent run VASP writes the density $n(\mathbf{r})$ on a real-space grid into CHGCAR. The format is plain text: a header copied from POSCAR (lattice + atoms), then a block of floating-point numbers giving $n$ at every grid point. The grid is set by NGXF, NGYF, NGZF in the OUTCAR and is typically twice the wavefunction grid.

1# Quickly inspect a CHGCAR
2head -10 CHGCAR                # POSCAR-style header
3sed -n '/^ *[0-9]\+ [0-9]\+ [0-9]\+$/,$p' CHGCAR | head  # find the grid line, dump first values
4
5# Check normalisation (should equal NELECT)
6python -c "
7import numpy as np
8from pymatgen.io.vasp import Chgcar
9chg = Chgcar.from_file('CHGCAR')
10print('integral n =', chg.data['total'].sum() * chg.structure.volume / chg.data['total'].size)
11"

2. PARCHG — partial densities (one band, one k-point, ...)

Set LPARD = .TRUE. in INCAR and ask for a specific band/k-point/spin range; VASP prints a PARCHG file in the same format. This is how you visualise the HOMO of a molecule, the wavefunction of a defect state, or the localisation of a Mn d orbital inside a CdSe matrix — pictures we will need in Chapter 7.

3. CHG — the running density (for restarts)

CHG is the same data as CHGCAR but truncated to fewer significant figures. VASP can read it back at the start of a new run via ICHARG = 1 (read from file) or ICHARG = 11 (read and freeze for non-self-consistent post-processing such as band structures and DOS).

1# INCAR — controlling what density VASP writes and reads
2LCHARG = .TRUE.    # write CHGCAR/CHG at the end (default)
3LWAVE  = .FALSE.   # don’t write WAVECAR — saves disk
4ICHARG = 2         # build the initial density from atomic charges (default for fresh runs)
5                   # set to 1 to read CHGCAR; 11 for non-SCF DOS / band-structure runs
6LPARD  = .FALSE.   # set to .TRUE. and add EINT/IBAND/KPUSE for partial densities

4. Sanity checks for any CHGCAR

A useful habit is to treat every charge-density file as data that should pass basic conservation checks before you interpret pictures from it.

• Integrate the total density and compare it with NELECT. A mismatch usually means you parsed the grid or volume factor incorrectly.
• Plot a planar average $\bar{n}(z) = A^{-1}\int n(x,y,z)\,dx\,dy$ for slabs, interfaces, and surfaces. It reveals charge spill-out and vacuum-region artefacts faster than a 3D isosurface.
• For charged cells or defect calculations, compare CHGCAR files with the same grid and structure before subtracting them. Density differences are meaningful only when the underlying real-space grids align.

5. Reading the density at the cusps

Kato's cusp condition lets you read off nuclear positions and charges directly from the density: every nucleus is a sharp peak, and the slope of the spherically averaged density at $r \to 0$ equals $-2 Z$. In a plane-wave code this cusp is only approximate (the basis cannot resolve a true cusp), but pseudopotential-augmented PAW datasets reconstruct it via partial-wave projections. We will return to this in Section 4.8.

Summary

1. The many-body wavefunction $\Psi(\mathbf{r}_1, \ldots, \mathbf{r}_N)$ lives in $3N$ dimensions and is hopeless to store directly. The electron density $n(\mathbf{r})$ lives in just three.
2. HK1: the ground-state density uniquely determines the external potential (and therefore the entire Hamiltonian and every observable). The map is bijective — no two different potentials can give the same density.
3. HK2: there is a universal energy functional $E_v[n]$ minimised by the true ground-state density. Every trial density gives an upper bound.
4. The unknown piece is $F[n] = T[n] + V_{ee}[n]$, conventionally decomposed as $F[n] = T_s[n] + J[n] + E_{xc}[n]$. The Hartree term is explicit in the density; $T_s$ is explicit through Kohn–Sham orbitals; and $E_{xc}[n]$ absorbs the remaining kinetic and interaction errors.
5. Thomas–Fermi (1927) gave the first concrete $F[n]$. It is qualitatively right but quantitatively poor — it over-binds H by 50% and cannot bind molecules at all. The fix is to handle the kinetic energy via orbitals (Kohn–Sham, Section 4.6) and bury the remaining many-body errors in $E_{xc}$.
6. In VASP the density lives in CHGCAR, PARCHG, and CHG. Every output file in electronic-structure post-processing is a derivative of the same three-dimensional object Hohenberg and Kohn proved was sufficient.

Next section we turn the existence theorem into an algorithm. Kohn–Sham introduce a fictitious non-interacting system whose orbitals carry the kinetic energy exactly, leaving only a (small) exchange–correlation correction to approximate. That single construction is the engine of every DFT calculation VASP has ever run — including the Mn:CdSe quantum dots waiting for us in Chapter 7.