The Many-Body Problem

Learning Objectives

Sections 4.1, 4.2 and 4.3 each told the same lie, well: that an electron in a crystal can be described as if it were alone. The free-electron gas, the nearly-free-electron model, and tight binding all solve a single-particle Schrödinger equation and hope the other 10²³ electrons average out. Sometimes they do. Mostly they don't, and the gap between “single-particle intuition” and “real material” is the central difficulty of the rest of this book.

This short section faces the lie head-on. We write down the full many-electron Hamiltonian, ask how big its wavefunction would have to be on a computer, and convince ourselves that the brute-force approach is hopeless — not difficult, not slow, physically impossible. We then walk the same hierarchy every electronic-structure code follows: Hartree ignores Pauli, Hartree–Fock repairs Pauli with a Slater determinant, and density-functional theory (Section 4.5) bypasses the wavefunction altogether by working with the density $\rho(\mathbf{r})$.

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

1. Write down the non-relativistic many-electron Hamiltonian in a fixed external potential and identify which term is responsible for everything difficult.
2. Quantify the curse of dimensionality — show that storing $\psi(\mathbf{r}_1, \dots, \mathbf{r}_N)$ on a grid scales as $M^{3N}$, exponential in the electron count.
3. State the Pauli antisymmetry principle and recognise why a Slater determinant is the simplest wavefunction that respects it.
4. Explain the difference between the Hartree and Hartree–Fock approximations, and connect the exchange integral $K$ to the cancellation of self-interaction.
5. Identify exchange and correlation as the two missing pieces of mean-field theory and read them off the energy decomposition of a real atom.
6. Map every approximation in this section to a corresponding tag (ALGO, LHFCALC, NELECT) in a real VASP INCAR, and understand why VASP solves Kohn–Sham equations rather than the bare many-body problem.

Where Sections 4.1–4.3 Left Off

In Section 4.1 we filled a Fermi sphere of non-interacting plane waves and got the heat capacity of sodium right within a factor close to 1. In 4.2 we let a periodic potential lift the zone-boundary degeneracies and conjured up band gaps for the price of a 2 × 2 diagonalisation. In 4.3 we abandoned the plane-wave basis and replaced it with atomic orbitals, recovering the same bands by hopping integrals.

Three different starting points, one stunningly common feature: each electron lived in its own private universe. The Hamiltonian in every case had the form

$\hat{H}_{\text{1e}} = -\frac{\hbar^2}{2 m_e}\nabla^2 + V_{\text{eff}}(\mathbf{r})$

with some effective potential $V_{\text{eff}}$ that pretended to encode the rest of the universe. None of those models has any way to know that two electrons can't occupy the same state, that they repel each other through the bare Coulomb interaction $1/|\mathbf{r}_i - \mathbf{r}_j|$, or that the wavefunction must change sign when two electrons swap places. Time to put all that back.

Modeling Assumptions

Before writing the full Hamiltonian, pin down the contract. We are not trying to solve all of quantum electrodynamics. We are solving the standard electronic-structure problem used by VASP and most solid-state DFT codes:

• Born–Oppenheimer nuclei: nuclei are fixed at positions $\mathbf{R}_I$ while electrons relax around them.
• Non-relativistic electrons: spin-orbit coupling, scalar relativistic corrections, and magnetic terms can be added later, but they are not part of the baseline Hamiltonian here.
• Coulomb interactions only: electrons repel through$1/|\mathbf{r}_i-\mathbf{r}_j|$, and nuclei attract through their external potential.
• Atomic units when convenient: later equations set$\hbar = m_e = e = 4\pi\varepsilon_0 = 1$, which removes constants without changing the physics.

The Real Hamiltonian — All The Terms We Have Been Hiding

For a non-relativistic, Born–Oppenheimer-frozen crystal of $N$ electrons in the field of fixed nuclei at positions $\mathbf{R}_I$ with charges $Z_I$, the electronic Hamiltonian is:

$\hat{H} = -\sum_{i=1}^{N} \frac{\hbar^2 \nabla_i^2}{2 m_e} - \sum_{i=1}^{N}\sum_{I} \frac{Z_I e^2}{4\pi\varepsilon_0 |\mathbf{r}_i - \mathbf{R}_I|} + \sum_{i<j} \frac{e^2}{4\pi\varepsilon_0 |\mathbf{r}_i - \mathbf{r}_j|}$

Read in plain English, three pieces:

• Kinetic energy $\hat{T} = -\sum_i \hbar^2 \nabla_i^2 / 2 m_e$ — separable: each electron contributes independently.
• External potential $\hat{V}_{\text{ext}} = \sum_i v(\mathbf{r}_i)$ where $v(\mathbf{r}) = -\sum_I Z_I/|\mathbf{r} - \mathbf{R}_I|$ (in atomic units). Also separable. This is the only term that remembers we are in a particular crystal and not a free electron gas.
• Electron–electron interaction $\hat{V}_{\text{ee}} = \sum_{i<j} 1/|\mathbf{r}_i - \mathbf{r}_j|$ — not separable. Every electron knows about every other electron. This single term is the entire many-body problem.

The wavefunction we are looking for has $3N$ spatial arguments plus $N$ spin labels: $\psi(\mathbf{r}_1\sigma_1, \dots, \mathbf{r}_N \sigma_N)$. It must satisfy the eigenvalue equation $\hat{H}\psi = E\psi$ and one non-negotiable additional rule we will derive in a moment: antisymmetry under exchange of any two electron labels.

The Wall: The Curse of Dimensionality

The Schrödinger equation has been around since 1926. Why hasn't someone simply solved it on a big computer? Let us count.

Suppose we represent each spatial coordinate on a grid of $M$ points per axis — say $M = 100$, hardly extravagant. For a single electron in 3D the wavefunction has $M^3 = 10^6$ amplitudes — a few megabytes. Comfortable.

Now add a second electron. The joint wavefunction lives on the product grid, so it has $M^{3 \cdot 2} = 10^{12}$ amplitudes: eight terabytes. Still painful but feasible.

Three electrons: $10^{18}$ amplitudes — eight exabytes, more than the world's combined hard-drive capacity. Four electrons: $10^{24}$ — about Avogadro's number of bytes. Ten electrons: more amplitudes than there are atoms in the human body. Eighteen electrons (a single Ar atom, or one CdSe formula unit): more than the number of atoms in the observable universe.

The hard limit: direct discretisation of the many-body wavefunction is not a question of slow computers. There is not enough matter in the visible universe to even store the wavefunction of a single noble-gas atom at modest resolution. Brute force is not just expensive — it is physically excluded.

Interactive: How Big Is the Many-Body Wavefunction?

The red bar below is the storage required for $\psi(\mathbf{r}_1, \dots, \mathbf{r}_N)$ on an $M$-per-axis 3D grid. The green bar is the storage required for the one-body density $\rho(\mathbf{r})$. Move the sliders. Watch the green bar barely twitch while the red one punches through the data on Earth, the atoms in your body, and the atoms in the Earth, in that order.

The Pauli Antisymmetry Principle

Electrons are fermions: spin-½ particles whose many-body wavefunction must be totally antisymmetric under the exchange of any two particles:

$\psi(\dots, \mathbf{x}_i, \dots, \mathbf{x}_j, \dots) = -\,\psi(\dots, \mathbf{x}_j, \dots, \mathbf{x}_i, \dots)$

where $\mathbf{x}_i = (\mathbf{r}_i, \sigma_i)$ bundles the spatial position and spin of electron $i$. This is not a calculational trick — it is an empirical fact about the universe (and a theorem in relativistic quantum field theory: the spin-statistics theorem). Observable consequences of the minus sign include the entire periodic table, the rigidity of solids, and the existence of white dwarfs.

From the equation above follows the famous corollary:

$\psi(\dots, \mathbf{x}, \dots, \mathbf{x}, \dots) = -\,\psi(\dots, \mathbf{x}, \dots, \mathbf{x}, \dots) \;\;\Longrightarrow\;\; \psi = 0$

Two electrons with the same spin cannot occupy the same point in space — the wavefunction is forced to vanish there. This is the geometric content of the Pauli exclusion principle, and it will show up below as a literal zero on the diagonal of a 2D plot.

Slater Determinants — Antisymmetry on a Page

How do we build a wavefunction that automatically antisymmetrises itself? In 1929 John Slater wrote down the cleanest answer: take a determinant.

Given $N$ orthonormal single-particle spin-orbitals $\chi_1, \chi_2, \dots, \chi_N$, define

$\psi_{\text{SD}}(\mathbf{x}_1, \dots, \mathbf{x}_N) = \frac{1}{\sqrt{N!}}\,\det \begin{pmatrix} \chi_1(\mathbf{x}_1) & \chi_2(\mathbf{x}_1) & \cdots & \chi_N(\mathbf{x}_1) \\ \chi_1(\mathbf{x}_2) & \chi_2(\mathbf{x}_2) & \cdots & \chi_N(\mathbf{x}_2) \\ \vdots & \vdots & \ddots & \vdots \\ \chi_1(\mathbf{x}_N) & \chi_2(\mathbf{x}_N) & \cdots & \chi_N(\mathbf{x}_N) \end{pmatrix}$

Determinants are antisymmetric under row swaps. Swapping electrons $i$ and $j$ corresponds to swapping rows $i$ and $j$ of the matrix above, which flips the sign of $\det$. Antisymmetry done. Determinants vanish when two rows are equal. Putting two electrons at the same $\mathbf{x}$ makes two rows equal, so $\psi_{\text{SD}} = 0$. Pauli enforced for free.

For the simplest case, two electrons in two spin-orbitals $\chi_a$, $\chi_b$:

$\psi_{\text{SD}}(\mathbf{x}_1, \mathbf{x}_2) = \frac{1}{\sqrt{2}}\big[\chi_a(\mathbf{x}_1) \chi_b(\mathbf{x}_2) - \chi_b(\mathbf{x}_1) \chi_a(\mathbf{x}_2)\big]$

Compare this to the naïve Hartree product $\chi_a(\mathbf{x}_1)\chi_b(\mathbf{x}_2)$: the second term is the antisymmetrising correction. The minus sign is small ink and enormous physics.

Interactive: Antisymmetry Under a Particle Swap

Below is the two-electron wavefunction $\psi(x_1, x_2)$ for two orbitals of a 1D infinite well: $\phi_a(x) = \sqrt{2}\sin(\pi x)$ and $\phi_b(x) = \sqrt{2}\sin(2\pi x)$. Pick the wavefunction style with the buttons at the top right. Drag the two probe sliders to read out $\psi(x_1, x_2)$ and its swapped twin $\psi(x_2, x_1)$. Look in particular at the yellow dashed diagonal — the Pauli line $x_1 = x_2$.

• Hartree product — the swap ratio is some generic non-±1 number; the diagonal is alive and well. Pauli is violated.
• Symmetric combination — the swap ratio is $+1$ identically; the diagonal is the peak of $|\psi|^2$. This is what bosons do, not electrons.
• Slater determinant — the swap ratio is $-1$ identically and the diagonal is machine-zero. Pauli enforced geometrically.

The Hartree Approximation: Mean Field Without Pauli

Even after enforcing antisymmetry, the determinantal wavefunction contains $N!$ terms — still hopeless for anything beyond a few atoms. The Hartree (1928) approximation cuts the Gordian knot by ignoring antisymmetry altogether and positing a simple product:

$\psi_{\text{H}}(\mathbf{r}_1, \dots, \mathbf{r}_N) = \phi_1(\mathbf{r}_1)\,\phi_2(\mathbf{r}_2)\cdots \phi_N(\mathbf{r}_N)$

Each electron then sees an average potential generated by the charge density of all the others:

$v_H^{(i)}(\mathbf{r}) = \sum_{j \neq i} \int \frac{|\phi_j(\mathbf{r}')|^2}{|\mathbf{r} - \mathbf{r}'|}\,d\mathbf{r}'$

and satisfies a self-consistent Schrödinger equation

$\left[-\tfrac{1}{2}\nabla^2 + v_{\text{ext}}(\mathbf{r}) + v_H^{(i)}(\mathbf{r})\right]\phi_i(\mathbf{r}) = \varepsilon_i\,\phi_i(\mathbf{r}).$

This is intuitive and computable: solve for the orbitals, build the density, recompute $v_H$, iterate until self-consistent. But it has two visible disasters. It violates Pauli — nothing forces the product to change sign under exchange. And it has a self-interaction error: in the formula $v_H^{(i)} = \sum_{j \neq i} \int |\phi_j|^2/r\,d\mathbf{r}'$ we excluded $j = i$ by hand, but in practice we usually compute $v_H = \sum_j \int |\phi_j|^2/r\,d\mathbf{r}'$ and tolerate the spurious self-repulsion. For a single electron this means it repels itself, lifting its own energy.

The Hartree–Fock Approximation: Mean Field Done Right

Hartree–Fock (1930) is the natural fix. Replace the simple product with a single Slater determinant of spin-orbitals, $\psi_{\text{HF}} = |\chi_1 \chi_2 \cdots \chi_N\rangle$, and minimise $\langle\psi|\hat{H}|\psi\rangle$ over all such determinants. The variational equations that come out look almost identical to Hartree, but with one extra non-local term:

$\hat{f}\,\chi_i(\mathbf{x}) = \Big[\hat{h} + \hat{J} - \hat{K}\Big]\chi_i(\mathbf{x}) = \varepsilon_i\,\chi_i(\mathbf{x})$

Here $\hat{J}$ is the same Hartree (Coulomb) operator as before, and $\hat{K}$ is the new exchange operator:

$\hat{K}\chi_i(\mathbf{x}) = \sum_j \chi_j(\mathbf{x}) \int \frac{\chi_j^*(\mathbf{x}')\chi_i(\mathbf{x}')}{|\mathbf{r} - \mathbf{r}'|}\,d\mathbf{x}'$

Notice that $\hat{K}$ is non-local: its action on $\chi_i$ at the point $\mathbf{x}$ depends on the value of $\chi_i$ at every other point $\mathbf{x}'$. That non-locality is the whole point: it is precisely what subtracts the spurious self-interaction of $\hat{J}$. For the special case $j = i$, the $\hat{J}$-self-interaction and the $\hat{K}$-self-exchange cancel exactly. Hartree–Fock is self-interaction free.

The energy expression splits into the Coulomb integral $J_{ij}$ and the exchange integral $K_{ij}$:

$E_{\text{HF}} = \sum_i h_{ii} + \tfrac{1}{2}\sum_{ij}\big(J_{ij} - K_{ij}\big)$

with

$J_{ij} = \iint \frac{|\chi_i(\mathbf{x}_1)|^2\,|\chi_j(\mathbf{x}_2)|^2}{|\mathbf{r}_1 - \mathbf{r}_2|}\,d\mathbf{x}_1\,d\mathbf{x}_2,$

$K_{ij} = \iint \frac{\chi_i^*(\mathbf{x}_1)\chi_j(\mathbf{x}_1)\,\chi_j^*(\mathbf{x}_2)\chi_i(\mathbf{x}_2)}{|\mathbf{r}_1 - \mathbf{r}_2|}\,d\mathbf{x}_1\,d\mathbf{x}_2.$

$K_{ij}$ is non-zero only when $\chi_i$ and $\chi_j$ share the same spin: opposite spins are orthogonal in spin space and the spin part of the integrand integrates to zero. So exchange is a same-spin effect.

The intuitive content of exchange: two electrons with the same spin actively avoid each other in space, even before any Coulomb repulsion is taken into account, because antisymmetrising their joint wavefunction forces it to vanish at coincidence. The hole that opens up around each electron is called the exchange hole; $K_{ii}$ is its electrostatic energy.

Exchange vs Correlation — What HF Gets and What It Misses

Hartree–Fock is exact for any system of non-interacting fermions in a fixed external potential, and an extremely good approximation for a single closed-shell atom. But it is built on a single Slater determinant, and that turns out not to be flexible enough.

The exact ground state can be written as a (typically infinite) linear combination of Slater determinants. The energy difference between the best single-determinant solution and the exact solution is, by definition, the correlation energy:

$E_c \equiv E_{\text{exact}} - E_{\text{HF}}$

$E_c$ is always negative (HF is variational, so it can only lie above the true energy). For light atoms it is small in absolute terms — about $-0.04$ Ha for helium, less than 2% of the total energy — but it is responsible for chemistry: bond strengths, reaction barriers, dispersion forces, and the difference between a fragile metal and a strong one. A “1% effect” that is bigger than every chemical bond energy on Earth.

Interactive: Anatomy of the Helium Atom Energy

Helium is the smallest non-trivial atom: two electrons, opposite spins, both in the 1s orbital. Its total energy splits cleanly into kinetic, electron–nucleus, classical Coulomb, exchange, and correlation pieces. Toggle the buttons to see what each approximation captures.

Three lessons jump out. One: the kinetic and electron-nucleus pieces are huge ($\sim$ 7 Ha), but they nearly cancel. Two: Hartree-Fock captures $\sim 99\%$ of the total energy in Helium — not because correlation is small but because kinetic and external dominate. Three: the leftover correlation energy $E_c \approx -0.04$ Ha is bigger than the first ionisation energy of caesium. Chemistry lives in the last percent.

Code Walk-Through: Two Electrons, Three Wavefunctions

We can build all three wavefunctions explicitly on a 1D grid and watch the abstract concepts turn into numbers. The script below constructs (1) a Hartree product and (2) a Slater determinant for two electrons in two box orbitals; verifies antisymmetry and the Pauli-line vanishing; computes the Coulomb ($J$) and exchange ($K$) integrals; and finally prints the storage required for the full $\psi$ on a 3D grid as $N$ grows. Click any line on the right to walk through it.

shape (200, 200) — first electron carries x₁ along rows, second carries x₂ along columns.

φ_a(x₁) φ_b(x₂)

φ_b(x₁) φ_a(x₂)

shape (200, 200), antisymmetric matrix: psi_S.T == −psi_S to within rounding.

transpose of psi_S — encodes ψ(x₂, x₁)

elementwise negation

Pairwise differences, shape (200, 200): D[i, j] = x_i − x_j

shape (200, 200): largest entry ≈ 1000 on the diagonal, smallest ≈ 1.0 on the corners.

|φ_a(x₁)|² · 1/|x₁−x₂| · |φ_b(x₂)|² — the double-integral kernel.

φ_a(x₁) φ_b(x₁) (1/|x₁−x₂|) φ_a(x₂) φ_b(x₂)

1import numpy as np
2
3# Two electrons in a 1D box of length L = 1, on a grid of M points.
4M = 200
5x = np.linspace(0, 1, M)
6dx = x[1] - x[0]
7
8# Single-particle orbitals: ground and first excited state of the box
9phi_a = np.sqrt(2) * np.sin(np.pi * x)
10phi_b = np.sqrt(2) * np.sin(2 * np.pi * x)
11
12# 1) Hartree product (no antisymmetry):  psi_H(x1, x2) = phi_a(x1) phi_b(x2)
13psi_H = np.outer(phi_a, phi_b)
14
15# 2) Slater determinant (antisymmetrized 2-electron wavefunction)
16psi_S = (np.outer(phi_a, phi_b) - np.outer(phi_b, phi_a)) / np.sqrt(2)
17
18# 3) Verify antisymmetry under x1 <-> x2:  psi_S(x2, x1) == -psi_S(x1, x2)
19print("antisymmetric? ", np.allclose(psi_S.T, -psi_S))
20
21# 4) On the Pauli line (x1 = x2) the Slater determinant must vanish identically
22print(f"max |psi_S(x, x)| = {np.abs(np.diag(psi_S)).max():.2e}")
23
24# 5) Coulomb (J) and exchange (K) integrals with softened 1/|x1 - x2|
25inv_r = 1.0 / (np.abs(x[:, None] - x[None, :]) + 1e-3)
26J = ((phi_a[:, None]**2) * inv_r * (phi_b[None, :]**2)).sum() * dx**2
27K = (phi_a[:, None] * phi_b[:, None] * inv_r
28     * phi_a[None, :] * phi_b[None, :]).sum() * dx**2
29print(f"J = {J:.4f}   K = {K:.4f}   HF e-e energy = J - K = {J - K:.4f}")
30
31# 6) Storage of full psi(r1,...,rN) on an M-per-axis 3D grid vs density rho(r)
32for N in [1, 4, 10, 18]:
33    full = M ** (3 * N) * 8
34    rho  = M ** 3 * 8
35    print(f"N={N:2d}: psi ~ {full:.1e} B,  rho ~ {rho:.1e} B")

Run it and the output reads:

1antisymmetric?  True
2max |psi_S(x, x)| = 1.11e-16
3J = 4.5023   K = 1.6177   HF e-e energy = J - K = 2.8846
4N= 1: psi ~ 6.4e+07 B,  rho ~ 6.4e+07 B
5N= 4: psi ~ 3.3e+28 B,  rho ~ 6.4e+07 B
6N=10: psi ~ 8.6e+69 B,  rho ~ 6.4e+07 B
7N=18: psi ~ 5.5e+125 B, rho ~ 6.4e+07 B

The first three lines are little theorems we have verified numerically: the Slater determinant is antisymmetric to machine precision, vanishes on the Pauli line to machine precision, and the exchange integral $K$ shaves a healthy chunk off the bare Coulomb $J$. The last four lines are the wall: by $N = 18$ the wavefunction would need 10¹²⁵ bytes — about $10^{45}$ times the number of atoms in the observable universe — while the density still fits in 64 MB.

The Road to DFT — Why a Density Is Enough

Hartree–Fock's problem is that a single Slater determinant misses correlation. Going beyond a single determinant (configuration interaction, coupled cluster) recovers correlation but at a cost that scales like $N^7$ or worse and runs out of steam at about a dozen heavy atoms. We need a different idea.

That idea is staring at us in the green bar of the curse-of-dimensionality plot. The one-body density

$\rho(\mathbf{r}) = N \int |\psi(\mathbf{r}, \mathbf{r}_2, \dots, \mathbf{r}_N)|^2\,d\mathbf{r}_2 \cdots d\mathbf{r}_N$

is a single 3D function, no matter how many electrons we have. It contains far less raw information than $\psi$. Hohenberg and Kohn proved in 1964 that it is nevertheless sufficient: in a fixed external potential, the ground-state density $\rho_0(\mathbf{r})$ uniquely determines every observable, including the wavefunction itself. The explosion in the red bar is irrelevant; the green bar carries all the physics. That theorem is the subject of Section 4.5.

One-line preview of DFT: instead of solving a Schrödinger equation for $\psi(\mathbf{r}_1, \dots, \mathbf{r}_N)$, solve a self-consistent set of one-electron equations for orbitals $\phi_i(\mathbf{r})$ whose density $\rho = \sum_i |\phi_i|^2$ equals the true many-body density. The exchange-correlation hole is folded into a single functional $E_{xc}[\rho]$.

VASP Connection: Why Real Codes Skip the Many-Body Wavefunction

VASP — like every production electronic-structure code — does not store $\psi(\mathbf{r}_1, \dots, \mathbf{r}_N)$. It cannot. What it stores is a set of single-particle orbitals $\phi_i(\mathbf{r})$ (Kohn–Sham orbitals, Section 4.6) and the resulting density $\rho(\mathbf{r})$. Three knobs in the INCAR directly correspond to choices we made above.

1. NELECT — how many electrons there are

The integer $N$ from the entire discussion. VASP usually reads it from the POTCAR (one valence count per species), but you can override it for charged cells.

1# INCAR — electron count for a 4-formula-unit CdSe supercell
2NELECT = 72        # 18 × 4  electrons in the cell
3ISPIN  = 2         # allow up- and down-spin orbitals to differ

2. ALGO and LHFCALC — which Hamiltonian to solve

The same ALGO tag that selects the iterative diagonaliser also picks the level of theory. The default is Kohn–Sham DFT. Setting LHFCALC = .TRUE. turns on Hartree–Fock exchange — the exact $\hat{K}$ we wrote above — either as part of a hybrid functional (PBE0, HSE) or as pure HF.

1# Pure Hartree-Fock (no correlation): exact exchange, no E_c
2LHFCALC = .TRUE.
3AEXX    = 1.0          # 100% exact exchange
4ALDAC   = 0.0          # 0% LDA correlation
5AGGAC   = 0.0          # 0% PBE correlation
6ALGO    = ALL          # robust algorithm for HF
7
8# HSE06 hybrid (24 of exact exchange, screened) — usually the right call
9LHFCALC = .TRUE.
10HFSCREEN = 0.2
11AEXX    = 0.25
12ALGO    = Damped

3. Why pure many-body methods do not appear in VASP

VASP supports many-body perturbation theory via GW (ALGO = GW0) and the Bethe–Salpeter equation, but it does not store the full configuration-interaction wavefunction — the curse of dimensionality bars the door. Practical codes climb only as high up Jacob's ladder as needed, and the rungs are: LDA → GGA → meta-GGA → hybrid → GW → CI/CC. VASP lives roughly on the hybrid rung; quantum-chemistry codes live higher. The choice is dictated by the size of the system and the property you need.

4. Reading total-energy lines

After every SCF cycle, VASP prints decompositions like

1free  energy   TOTEN  =      -240.812345 eV
2  energy without entropy =  -240.812345  energy(sigma->0) =  -240.812345
3
4  alpha Z        PSCENC =      ...
5  Ewald energy  TEWEN   =      ...
6  -1/2 Hartree  DENC    =      ...     <-- this is -E_H/2 (avoids double-counting)
7  -exchange  EXHF       =      ...     <-- exact exchange when LHFCALC = .TRUE.
8  -V(xc)+E(xc) XCENC    =      ...     <-- the DFT exchange-correlation contribution

Every name on the right of the equals sign is something we have defined in this section: $E_H$ is the Hartree (Coulomb) energy, $E_x$ the exchange, and the DFT exchange-correlation lump $E_{xc}$ is the functional approximation we will dissect in Section 4.7.

Summary

1. The full electronic Hamiltonian has three pieces — kinetic, external, and electron-electron Coulomb. Only the e-e term is non-separable, and it alone constitutes the many-body problem.
2. Direct discretisation of the many-body wavefunction scales as $M^{3N}$. By a dozen electrons it exceeds every storage budget in the universe; this is a hard physical wall, not an engineering inconvenience.
3. Electrons are fermions: their many-body wavefunction must be antisymmetric under exchange. The simplest wavefunction that enforces this is a Slater determinant. Antisymmetry forces $\psi$ to vanish whenever two electrons coincide — the geometric content of the Pauli principle.
4. Hartree theory is a product wavefunction with self-consistent mean-field potentials. It captures classical Coulomb but violates Pauli and suffers self-interaction error.
5. Hartree-Fock replaces the product with a single Slater determinant and adds the non-local exchange operator $\hat{K}$. Exchange exactly cancels self-interaction and acts as a same-spin repulsion.
6. The energy missing from the best Slater-determinant solution is the correlation energy $E_c = E_{\text{exact}} - E_{\text{HF}}$. It is small in absolute terms but dominates chemistry.
7. The way out is the one-body density $\rho(\mathbf{r})$: a single 3D function that, by Hohenberg–Kohn, contains all ground-state information. Section 4.5 turns that statement into a usable theorem; Section 4.6 builds the Kohn–Sham equations VASP actually solves.

We have replaced one impossible problem (10¹²⁵-byte wavefunction) with a manageable one (a few-MB density plus a clever functional). The price will be that the functional itself is approximate, and the rest of this chapter is essentially a long argument about how to design and use it.