Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

Section 4.5 ended on a frustrating high note. Hohenberg and Kohn proved that the entire ground-state physics of an interacting electron system is encoded in the density $n(\mathbf{r})$ : there is a universal functional $E[n]$ , and minimising it gives the exact ground-state energy. But they did not tell us what $E[n]$ looks like. In particular the kinetic-energy functional $T[n]$ — the kinetic energy of the interacting system as a functional of its density alone — is unknown and almost certainly nasty.

Walter Kohn and Lu Jeu Sham's 1965 paper turned this stuck problem into the most successful framework in computational physics. Their idea was simple, audacious, and exact: replace the interacting electrons with a fictitious system of non-interacting electrons that has the same density. That auxiliary system is solvable. The price you pay is the need for one new functional, the exchange-correlation $E_{xc}[n]$ , and one iterative loop. That deal is what every modern DFT code — VASP included — actually solves.

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

State the Kohn–Sham ansatz: an interacting density $n(\mathbf{r})$ is reproduced by a system of N non-interacting fermions in some effective potential $V_{\rm eff}$ .
Decompose the energy functional as $E[n] = T_s[n] + V_{\rm ext}[n] + E_H[n] + E_{xc}[n]$ and explain what each term is.
Derive the Kohn–Sham equations $(-\tfrac{1}{2}\nabla^2 + V_{\rm eff})\psi_i = \varepsilon_i \psi_i$ from the variational principle.
Recognise why these equations are self-consistent — $V_{\rm eff}$ depends on the orbitals which depend on $V_{\rm eff}$ — and walk through the SCF iteration loop step by step.
State precisely what KS eigenvalues do and don't mean — the band-gap problem is not a bug in DFT, it is a feature of approximate $E_{xc}$ .
Read a 1D Kohn–Sham solver in NumPy, watch its self-consistency loop converge in real time, and connect every keyword in a VASP INCAR (ALGO, NELM, EDIFF) to a specific line of that solver.

Where Section 4.5 Left Us Stranded

Recall the two Hohenberg–Kohn theorems from 4.5:

HK-1: the ground-state density $n_0(\mathbf{r})$ determines the external potential $V_{\rm ext}$ uniquely (up to a constant), and therefore the entire many-body Hamiltonian and all its observables.
HK-2: there exists a universal energy functional $E[n] = F[n] + \int V_{\rm ext}(\mathbf{r}) n(\mathbf{r})\,d\mathbf{r}$ whose minimum (over all densities integrating to N electrons) is the exact ground-state energy.

The universal functional $F[n] = T[n] + V_{ee}[n]$ contains the kinetic energy and the electron–electron interaction. Both pieces are functionals of the density alone, by HK-1 — but knowing they exist is not the same as having a usable formula. In fact there is a long history of trying to write $T[n]$ down: Thomas–Fermi (1927) used the free-electron kinetic-energy density $T_{\rm TF}[n] = C_{\rm TF}\int n^{5/3}\,d\mathbf{r}$ , which is exact for a uniform gas but spectacularly wrong for atoms — it cannot bind a hydrogen molecule.

The functional

E[n]

exists but is unknown. Without a good

T[n]

we cannot minimise

E[n]

. We need a way to compute the kinetic energy exactly — at least for some carefully chosen reference system — without knowing

T[n]

as a functional of n.

The Kohn–Sham Trick: A Fictitious System

The key observation is this: for a system of non-interacting electrons in an external potential $V_s(\mathbf{r})$ , computing the kinetic energy is trivial. The ground state is a single Slater determinant of the N lowest single-particle orbitals $\psi_i$ , and:

$T_s = -\frac{1}{2}\sum_{i=1}^{N} \langle\psi_i | \nabla^2 | \psi_i\rangle = \frac{1}{2}\sum_i \int |\nabla\psi_i|^2\,d\mathbf{r}$

The subscript “s” stands for single-particle or Slater — Kohn and Sham's original notation. $T_s$ is the kinetic energy of the non-interacting system, expressed as an explicit functional of the orbitals. It is not equal to $T[n]$ (the kinetic energy of the real interacting system at the same density) — these differ — but it is something we can actually compute.

The Kohn–Sham ansatz now reads:

For any interacting density $n(\mathbf{r})$ , there exists a non-interacting system in some effective potential $V_{\rm eff}(\mathbf{r})$ whose ground state has the very same density. We work with that auxiliary system and let $V_{\rm eff}$ absorb the difference.

The orbitals $\psi_i$ of the non-interacting system are called Kohn–Sham orbitals; the eigenvalues $\varepsilon_i$ are Kohn–Sham eigenvalues. They are mathematical bookkeeping devices that reproduce the right density. We will be very careful in § What KS Orbitals Actually Mean about what they do and don't correspond to in the real system.

In principle the existence of

V_{\rm eff}

for every reasonable density is itself an assumption (called v-representability). For interacting electrons in nature, this has never been proven in full generality — and there are pathological examples. For all practical materials including Mn:CdSe, the Kohn–Sham mapping works so well that we treat it as gospel.

Interactive: The Mapping Visualised

Two systems, one density. Click either box to see what each picture gives us — and where each picture fails on its own.

Splitting the Energy Functional

Now we must write the energy of the real interacting system in terms of quantities we can actually compute. The Kohn–Sham decomposition is:

$E[n] = T_s[n] + \int V_{\rm ext}(\mathbf{r}) n(\mathbf{r})\,d\mathbf{r} + \frac{1}{2}\int\!\!\int \frac{n(\mathbf{r}) n(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}\,d\mathbf{r}\,d\mathbf{r}' + E_{xc}[n]$

Term by term:

Term	Symbol	What it is	Computable?
Kinetic	T_s[n]	kinetic energy of the non-interacting KS system	Yes — sum of orbital kinetic energies
External	∫V_ext n dr	electrons feeling the nuclei (or pseudopotential)	Yes — straight integral
Hartree	E_H = ½∫∫ n n' / \|r-r'\|	classical Coulomb repulsion of the density with itself	Yes — solve Poisson eqn or convolve
Exchange-correlation	E_xc[n]	everything else: T-T_s + non-classical part of V_ee	No — must be approximated

The brilliance of the decomposition is that three of the four terms are exact and computable; the only unknown is $E_{xc}[n]$ . By definition,

$E_{xc}[n] = \underbrace{(T[n] - T_s[n])}_{\text{kinetic correlation}} + \underbrace{(V_{ee}[n] - E_H[n])}_{\text{quantum exchange + correlation in V}_{ee}}$

— the small piece of the truth that we don't know exactly. The whole DFT industry, from LDA in 1965 to PBE in 1996 to SCAN in 2015 to modern hybrids, is the search for better approximations to $E_{xc}[n]$ . The next section ( 4.7 Exchange–Correlation Functionals) is dedicated to that zoo.

Why dump everything difficult into $E_{xc}$ ? Because

E_{xc}

is small — typically 5–10% of the total energy. Approximating a small thing with a few percent error gives a sub-percent total energy error, which is chemistry-grade accuracy. The Kohn–Sham decomposition is, in a sense, the best perturbation theory ever invented.

Deriving the Kohn–Sham Equations

We minimise $E[n]$ over densities of the form $n = \sum_i |\psi_i|^2$ with the constraint $\langle \psi_i | \psi_j \rangle = \delta_{ij}$ . Lagrange multipliers $\varepsilon_{ij}$ enforce the orthonormality. After diagonalising the Lagrange multiplier matrix, the variational equations for each orbital decouple:

$\frac{\delta E[n]}{\delta \psi_i^*(\mathbf{r})} = \varepsilon_i \psi_i(\mathbf{r})$

Computing the functional derivative term by term:

From $T_s = -\frac{1}{2}\sum_i \langle\psi_i|\nabla^2|\psi_i\rangle$ : $\delta T_s/\delta\psi_i^* = -\tfrac{1}{2}\nabla^2 \psi_i$ .
From $\int V_{\rm ext} n$ with $n = \sum |\psi_j|^2$ : $\delta/\delta\psi_i^* = V_{\rm ext}(\mathbf{r}) \psi_i$ .
From $E_H[n]$ : $\delta E_H/\delta\psi_i^* = V_H(\mathbf{r})\psi_i$ where $V_H(\mathbf{r}) = \int n(\mathbf{r}')/|\mathbf{r}-\mathbf{r}'|\,d\mathbf{r}'$ .
From $E_{xc}[n]$ : by definition $V_{xc}(\mathbf{r}) \equiv \delta E_{xc}/\delta n(\mathbf{r})$ , so $\delta E_{xc}/\delta\psi_i^* = V_{xc}\psi_i$ .

Adding these up, the variational principle delivers the Kohn–Sham equations:

$\boxed{\;\left[-\tfrac{1}{2}\nabla^2 + V_{\rm eff}(\mathbf{r})\right]\psi_i(\mathbf{r}) = \varepsilon_i\,\psi_i(\mathbf{r})\;}$

with the effective potential

$V_{\rm eff}(\mathbf{r}) = V_{\rm ext}(\mathbf{r}) + V_H(\mathbf{r}) + V_{xc}(\mathbf{r})$

These look like a single-particle Schrödinger equation. The trick is the potential: $V_{\rm eff}$ already contains all the average effects of every other electron through $V_H$ (classical Coulomb) and $V_{xc}$ (quantum mechanics). What was a coupled $3N$ -dimensional Schrödinger problem in $\Psi(\mathbf{r}_1,\dots,\mathbf{r}_N)$ has become N independent 3D problems, all sharing the same potential.

The Kohn–Sham equations are formally exact — every approximation lives entirely inside

V_{xc}

. If you knew the true exchange-correlation potential, the KS density would be the exact density of the interacting system, and the total energy would be exact.

Anatomy of the Effective Potential

Three pieces, three personalities:

Component	Sign	Physical role	Depends on n?
V_ext(r)	negative near nuclei	Pulls electrons toward nuclei (or pseudoatoms)	no — fixed by the geometry
V_H(r)	positive everywhere	Classical mean-field repulsion: each electron repelled by the average density of all others	yes — V_H = ∫n(r')/\|r-r'\| dr'
V_xc(r)	negative where n is large	Pauli exchange + dynamic correlation: corrects the Hartree double-counting and adds quantum effects	yes — V_xc = δE_xc/δn

Two observations are crucial. First, both $V_H$ and $V_{xc}$ are functionals of the density that the orbitals themselves produce — that is the seed of the self-consistency problem we tackle next. Second, $V_H$ contains a famous unphysical feature called the self-interaction: an electron is repelled by its own density. The exchange part of $V_{xc}$ is supposed to cancel this exactly — and does, in Hartree–Fock or in the exact KS formulation, but only approximately in LDA / GGA. We will return to this when we meet hybrid functionals in Chapter 6.

The Self-Consistency Problem

Look again at the Kohn–Sham equation. The orbitals $\psi_i$ depend on $V_{\rm eff}$ ; but $V_{\rm eff}$ contains $V_H[n]$ and $V_{xc}[n]$ ; and $n = \sum_i |\psi_i|^2$ . The unknown appears on both sides.

We cannot solve this in closed form. The standard answer is iteration: guess a density, build $V_{\rm eff}$ , diagonalise to get new orbitals, build a new density, repeat. The fixed point of this iteration — where input density equals output density — is the self-consistent Kohn–Sham solution.

Interactive: The SCF Loop

Click each node to see the equation for that step. The whole DFT engine is just this loop running until the density stops changing.

In practice you do not just replace the old density with the new one (which often oscillates). Instead, you mix: $n_{i+1} = (1-\alpha) n_i + \alpha n_{\rm new}$ with $\alpha$ in the range 0.1–0.5. Modern codes use much fancier schemes — Pulay (DIIS) or Broyden — that extrapolate from a history of previous steps and converge in a fraction of the iterations.

Convergence is monitored by two numbers: the density residual

\|n_{i+1} - n_i\|

and the energy change

|E_{i+1} - E_i|

. VASP's EDIFF is the threshold on the energy; a value of

10^{-5}

eV is overkill for a structure relaxation but barely tight enough for phonons.

Interactive 1D Solver: Watch Convergence Happen

The widget below is a real, working 1D Kohn–Sham solver running live in your browser. Four electrons live in a harmonic well $V_{\rm ext}(x) = \tfrac{1}{2}x^2$ ; you watch the SCF cycle build up. At iteration 1 there is no electron–electron interaction yet (the density starts at zero), so $V_{\rm eff} = V_{\rm ext}$ exactly — the bare harmonic well. As iterations proceed the electrons begin to repel each other through $V_H$ (which lifts the well at the centre) and to feel exchange through $V_{xc}$ (which deepens it back). The density adjusts in response, and the eigenvalues drift to their self-consistent values.

Toggle the components on and off to dissect $V_{\rm eff}$ . Step through with prev / next, or hit run SCF to watch the loop play itself out.

Loading 1D Kohn–Sham SCF solver…

What you should see. Iteration 1: the density emerges from the bare harmonic well, two-bump structure (the two lowest oscillator orbitals doubly occupied). Iterations 2–4: the electrons start “feeling” each other; the centre of the density flattens, V_eff develops a hump. Iterations 5+: small corrections, the residual $\Delta n$ plummets by orders of magnitude. By iteration 9 the density and the potential it generates agree to four decimal places. That is what self-consistency looks like.

What KS Orbitals Actually Mean (and Don’t)

Kohn–Sham gives us a wealth of objects: orbitals $\psi_i$ , eigenvalues $\varepsilon_i$ , occupations $f_i$ , a density $n(\mathbf{r})$ , a total energy $E$ . Which of these are real?

Quantity	Status	Comment
n(r)	exact (with exact E_xc)	The whole framework is built to reproduce this. In approximate DFT it is correct to within a few percent.
E (total energy)	exact (with exact E_xc)	Trust this — energy differences between DFT calculations are quantitatively meaningful.
ψ_i (orbitals)	auxiliary	Not the real wavefunction — a bookkeeping device. But the highest occupied orbital usually looks similar to the HOMO and is hugely useful for chemical intuition.
ε_i (eigenvalues)	auxiliary, except ε_HOMO	By Janak's theorem in exact DFT, ε_HOMO equals the negative of the ionization potential. Other eigenvalues are NOT one-electron excitation energies.
KS band gap ε_LUMO − ε_HOMO	wrong by ~50% in LDA/GGA	The famous DFT band-gap problem — see below.

The band-gap problem in one sentence

The exact KS gap $\varepsilon_{\rm LUMO} - \varepsilon_{\rm HOMO}$ differs from the true fundamental gap by a quantity called the derivative discontinuity $\Delta_{xc}$ of the exchange-correlation functional. LDA and GGA functionals have $\Delta_{xc} = 0$ (they are smooth in n), so they miss this contribution. As a result, LDA and PBE underestimate semiconductor band gaps by 30–50%: bulk Si has a true gap of 1.17 eV, PBE gives 0.6 eV; bulk CdSe has 1.74 eV experimentally, PBE gives ~0.7 eV.

The fix is either to use a better functional (HSE06, B3LYP, SCAN) — Chapter 6 covers these — or to climb out of DFT entirely and do $GW$ (Section 5.10). In the meantime every DFT band structure plot you read should come with a mental warning sticker: gap likely underestimated.

The total energies and densities of the same materials are nonetheless within a percent of experiment. Why does the same functional get the energy right and the gap wrong? Because the gap is a derivative of the energy with respect to particle number — and approximations are usually softer in the function than in its derivatives.

Code Walk-Through: A Working 1D Kohn–Sham Solver

Time to read code that actually solves the Kohn–Sham equations. The script below is the engine behind the interactive widget you just used: it builds a finite-difference Hamiltonian on a 100-point grid, defines the Hartree and xc potentials, and runs an SCF loop with linear mixing. Every line is annotated — click any line on the right and the matching explanation pops up on the left.

1D Kohn–Sham SCF Solver — Interactive Trace

🐍ks_1d.py

Explanation(34)

Code(42)

1import numpy as np

NumPy is the workhorse of any DFT prototype. We need ndarray for vectorized arithmetic on the real-space grid, np.diag and broadcasting to assemble the kinetic-energy matrix, and np.linalg.eigh — a fast, numerically stable diagonaliser for Hermitian matrices — to solve the Kohn-Sham eigenvalue problem at each SCF step.

EXECUTION STATE

numpy = Numerical Python — ndarray, broadcasting, linear algebra, FFTs, random numbers

as np = Standard alias so we can write np.linspace, np.diag, np.linalg.eigh, etc.

3Comment: the physical setup

Reminds us that everything is in atomic units (ℏ = m_e = e = 1). The system is a 1D 'toy atom': 4 electrons in a harmonic external potential. Closed-shell, spin-paired, so we only need to track 2 spatial orbitals each holding 2 electrons.

4N = 100

Number of real-space grid points. The Kohn-Sham equations are PDEs; we discretize them on a uniform grid and turn them into a finite-dimensional matrix problem. More points = better resolution, slower diagonalization (cost ~ N³).

EXECUTION STATE

N = 100

→ role = Grid resolution. A 100×100 Hamiltonian matrix is trivial; in real VASP this becomes the plane-wave basis size of order 10⁴–10⁵.

5x = np.linspace(-5.0, 5.0, N)

Builds the 1D grid. We pick a box [-5, 5] in atomic units; the harmonic well V_ext = ½x² has wavefunctions that decay like e^{-x²/2}, so the lowest few orbitals are essentially zero outside this range.

EXECUTION STATE

📚 np.linspace(start, stop, num) = Returns ndarray of `num` evenly spaced values from start to stop, INCLUSIVE on both ends. Different from np.arange which uses a step size.

⬇ arg 1: start = -5.0 = Left edge of the box, in Bohr radii

⬇ arg 2: stop = 5.0 = Right edge of the box, in Bohr radii

⬇ arg 3: num = N = 100 = Total number of grid points (both endpoints included)

⬆ result: x = [-5.000 -4.899 -4.798 … 4.798 4.899 5.000] shape (100,)

6h = x[1] - x[0]

Grid spacing. With 100 points spanning 10 atomic units, h ≈ 10/99 ≈ 0.1010. This appears everywhere: in the kinetic-energy finite difference (1/h²), in normalization (∫|ψ|² dx ≈ Σ |ψ_i|² h), and in the Hartree integral.

EXECUTION STATE

x[1] - x[0] = Difference between two consecutive grid points — uniform for linspace

h = 0.1010 (= 10/99)

→ role = Quadrature weight: ∫f(x)dx ≈ Σ_i f(x_i) · h. Also controls finite-difference accuracy: error scales like h².

7n_occ = 2

Number of spatial orbitals to fill. With closed-shell spin pairing, each spatial orbital holds 2 electrons, so n_occ = 2 means 4 total electrons. In a real DFT code this corresponds to NELECT/2 for non-spin-polarized calculations.

EXECUTION STATE

n_occ = 2

→ corresponding electron count = 2 × 2 = 4 electrons (closed shell)

8V_ext = 0.5 * x**2

The external potential — what the electrons would feel without each other. We choose a 1D harmonic oscillator V(x) = ½x², the simplest non-trivial confining potential. In a real solid, V_ext is the sum of nuclear Coulomb potentials (or pseudopotentials).

EXECUTION STATE

x**2 = Element-wise square of the grid: [25, 23.99, …, 23.99, 25]

0.5 * x**2 = Multiply by ½ — broadcasting scalar over array

V_ext = [12.500 11.998 … 0.0026 0.0 0.0026 … 11.998 12.500] shape (100,)

→ physical role = The 'nuclear pull'. Without electron-electron interaction, eigenvalues would be the harmonic-oscillator levels 0.5, 1.5, 2.5, … (atomic units).

10Comment: kinetic operator setup

Marks the start of the kinetic-energy matrix construction. We use the 3-point central finite difference for the second derivative: ψ''(x_i) ≈ [ψ_{i-1} - 2ψ_i + ψ_{i+1}]/h². Multiplying by -½ gives the kinetic operator.

11T = np.diag(np.full(N, 1.0/h**2)) + …

Builds the (N×N) finite-difference kinetic-energy matrix. The 3-point stencil for -½ d²/dx² gives diagonal entries +1/h² and immediate off-diagonal entries -1/(2h²). The result is a symmetric tridiagonal matrix.

EXECUTION STATE

📚 np.diag(v) = If v is 1D: builds a square matrix with v on the main diagonal. If v is 2D: extracts the diagonal as a 1D array. Optional second arg `k` shifts to super- (k>0) or sub-diagonal (k<0).

📚 np.full(shape, fill) = Returns array of given shape filled with `fill`. Faster than np.ones * fill.

→ diag values = 1.0 / h² = 1.0 / 0.0102 ≈ 98.0 — appears N times on the main diagonal

→ off-diag values = -0.5 / h² ≈ -49.0 — appears N-1 times on each immediate off-diagonal

12+ np.diag(np.full(N-1, -0.5/h**2), 1)

Adds the super-diagonal entries (one above the main diagonal). The second argument k=1 tells np.diag to shift the placement up. Length is N-1 because the super-diagonal of an (N×N) matrix has N-1 entries.

EXECUTION STATE

⬇ arg: k = 1 = Place values on the SUPER-diagonal (one row above the main diagonal). k = -1 would be the sub-diagonal.

→ resulting block (top-left 4×4 of T) =

[[ 98.01  -49.01    0.00    0.00 ]
 [-49.01   98.01  -49.01    0.00 ]
 [  0.00  -49.01   98.01  -49.01 ]
 [  0.00    0.00  -49.01   98.01 ]]

13+ np.diag(np.full(N-1, -0.5/h**2), -1)

Adds the sub-diagonal entries (one below the main diagonal), making T symmetric. Identical values to the super-diagonal because -½ d²/dx² is a self-adjoint operator under appropriate boundary conditions (here: hard-wall, ψ vanishing at x = ±5).

EXECUTION STATE

⬇ arg: k = -1 = Place values on the SUB-diagonal (one row below the main diagonal).

T symmetric? = Yes — T[i,j] = T[j,i] for all i,j. Required for np.linalg.eigh to give real eigenvalues.

15def hartree(n) → np.ndarray

Computes the classical electrostatic (Hartree) potential generated by an electron density n(x). This is the average Coulomb repulsion that one electron feels from the cloud of all the others. The 'soft' parameter regularises the 1/|x-x'| singularity (in 3D the kernel is naturally finite; in 1D Coulomb diverges, so we soften it).

EXECUTION STATE

⬇ input: n (shape N,) = The current electron density on the grid, in units of e/a₀. Must be ≥ 0 for sensible physics.

⬆ returns = np.ndarray (shape N,) — V_H evaluated at each grid point, in Hartree units.

16Docstring: V_H = ∫ n(x') / √((x-x')² + soft²) dx'

Describes the integral being approximated. In 3D the Hartree potential is V_H(r) = ∫ n(r')/|r-r'| dr'; in 1D this would diverge logarithmically, so we use a soft-Coulomb regularisation 1/√((x-x')² + soft²).

17soft = 0.5

Softening parameter. Removes the unphysical 1D Coulomb singularity at x = x'. In the limit soft → 0 we recover the bare 1/|x-x'| kernel; with soft = 0.5 the kernel is bounded by 1/0.5 = 2 at coincidence.

EXECUTION STATE

soft = 0.5

→ kernel max = 1 / soft = 2.0 (at x = x')

18dx = x[:, None] - x[None, :]

Creates an N×N matrix of pairwise distances using NumPy broadcasting. x[:, None] reshapes x to shape (N, 1) — a column. x[None, :] reshapes to (1, N) — a row. Subtraction broadcasts to (N, N): dx[i, j] = x_i - x_j.

EXECUTION STATE

📚 array[:, None] = Adds a new axis at position 1. Shape (N,) → (N, 1). Equivalent to array.reshape(-1, 1) or array[:, np.newaxis].

📚 array[None, :] = Adds a new axis at position 0. Shape (N,) → (1, N).

→ broadcasting rule = (N, 1) - (1, N) → (N, N): each row of the result is x_i minus the entire x array.

⬆ dx =

(100, 100) matrix — element [i,j] = x[i] - x[j]. Diagonal is zero; off-diagonals symmetric in magnitude.

19K = 1.0 / np.sqrt(dx**2 + soft**2)

Builds the Coulomb kernel matrix. Each entry K[i, j] = 1 / √((x_i - x_j)² + soft²). This is the soft-Coulomb interaction between a unit charge at x_i and another at x_j.

EXECUTION STATE

📚 np.sqrt(arr) = Element-wise square root.

→ at coincidence (i=j) = K[i,i] = 1/√(0 + 0.25) = 1/0.5 = 2.0

→ far apart (|i-j| large) = K[i,j] ≈ 1/|x_i - x_j| (soft becomes negligible)

→ K shape = (100, 100) symmetric matrix

20return K @ n * h

Performs the convolution V_H(x_i) = h Σ_j K[i,j] n[j]. Matrix-vector multiplication K @ n collapses the j index; multiplying by h is the Riemann-sum quadrature weight that turns the discrete sum into an integral approximation.

EXECUTION STATE

@ = Python matrix-multiplication operator (PEP 465). For (N,N) @ (N,) → (N,).

h (scalar) = Grid spacing 0.1010 — quadrature weight

⬆ return: V_H (N,) = Hartree potential on the grid. Peaks where n is large; smooth and positive everywhere.

22def vxc(n) → np.ndarray

Computes the local exchange-correlation potential V_xc(x) = δE_xc/δn evaluated at the local density. We use a toy LDA-like form V_xc = -n^(1/3) — the same scaling as the 3D LDA exchange but without the proper coefficient (the 1D Coulomb problem doesn't have a literal LDA, so this is purely illustrative).

EXECUTION STATE

⬇ input: n = Density on the grid (N,)

⬆ returns = V_xc(x) at each grid point — negative everywhere n > 0, deepest where n is largest

24return -1.0 * np.cbrt(np.maximum(n, 0.0))

V_xc(x_i) = -[n(x_i)]^(1/3). The inner np.maximum(n, 0) clamps any tiny negative values (from numerical noise during mixing) to zero, since cbrt of a negative would still work but isn't physical here.

EXECUTION STATE

📚 np.cbrt(arr) = Element-wise cube root. Defined for negative inputs (cbrt(-8) = -2), unlike sqrt.

📚 np.maximum(a, b) = Element-wise max — for each i, returns max(a_i, b_i). Different from np.max which reduces along an axis.

→ why clamp at 0? = After linear mixing, density can dip slightly below zero by O(10⁻¹⁰). cbrt would still give a small negative, but the xc potential should be -|n|^(1/3) for non-negative n. Clamping is a safe one-line guard.

⬆ return: V_xc = Same shape as n. Negative wherever n > 0; zero wherever n = 0.

26def solve_ks(V_eff) → (eps, psi)

Solves the Kohn-Sham eigenvalue problem at fixed V_eff. The Hamiltonian H = T + V_eff is built once and diagonalised; we keep only the lowest n_occ eigenpairs.

EXECUTION STATE

⬇ input: V_eff (shape N,) = The current effective potential = V_ext + V_H + V_xc. Different at every SCF iteration.

⬆ returns = Tuple (eps, psi) — eps shape (n_occ,) eigenvalues, psi shape (N, n_occ) orbitals as columns

28H = T + np.diag(V_eff)

Builds the full Hamiltonian matrix. T is the (constant) kinetic operator built once at the top of the script; np.diag(V_eff) puts the (variable) effective potential on the diagonal. Total H is real symmetric, tridiagonal+diagonal = tridiagonal.

EXECUTION STATE

T (N,N) =

Kinetic-energy matrix — same every iteration

np.diag(V_eff) (N,N) =

Diagonal matrix with V_eff entries on the diagonal, zeros elsewhere

→ H structure = Tridiagonal: H[i,i] = 1/h² + V_eff[i], H[i,i±1] = -1/(2h²), all else 0

29eps, psi = np.linalg.eigh(H)

The single most important line in the whole script. np.linalg.eigh diagonalises a Hermitian (here real-symmetric) matrix and returns its eigenvalues and orthonormal eigenvectors. Every SCF iteration calls this once.

EXECUTION STATE

📚 np.linalg.eigh(M) = Returns (w, v): w is the 1D ndarray of N eigenvalues sorted in ASCENDING order; v is the (N,N) matrix of orthonormal eigenvectors stored as COLUMNS, so v[:, i] corresponds to w[i].

→ why eigh, not eig? = eigh exploits symmetry: O(N³) but with a smaller prefactor, guaranteed real eigenvalues, and orthogonal eigenvectors. eig is for general matrices and is slower + may return complex output.

⬆ result: eps = All N eigenvalues, ascending. The first few (after SCF) might be e.g. [+1.42, +2.51, +3.60, …]

⬆ result: psi =

(N, N) matrix. psi[:, 0] is ψ_0 (lowest), psi[:, 1] is ψ_1, etc. Columns are orthonormal under the Euclidean dot product.

30psi = psi / np.sqrt(h)

Normalisation rescaling. eigh gives eigenvectors with Σ_i |v_i|² = 1 (Euclidean norm), but we want continuous normalisation ∫|ψ(x)|² dx = 1, which on the grid is Σ_i |ψ(x_i)|² · h = 1. Dividing by √h fixes this: now Σ_i |ψ_i|² h = 1.

EXECUTION STATE

→ before = Σ |psi[:, 0]|² = 1.0 (Euclidean)

→ after = Σ |psi[:, 0]|² · h = 1.0 (proper continuum normalisation)

→ why this matters = Without rescaling, the density n = Σ |ψ|² would integrate to N_occ/h instead of N_occ. Wrong electron count = wrong everything.

31return eps[:n_occ], psi[:, :n_occ]

Slices off only the lowest n_occ = 2 orbitals — the ones we'll fill with electrons. The other N-n_occ eigenpairs are unoccupied (the 'virtual' or 'conduction' states in solid-state language) and are discarded here for efficiency. In a real DFT code the unoccupied bands are kept around for response calculations.

EXECUTION STATE

📚 array[:n] = Slice: first n elements along axis 0

📚 array[:, :n] = Slice: all rows, first n columns

⬆ return: eps[:2] = (eps_0, eps_1) — the two occupied eigenvalues

⬆ return: psi[:, :2] =

(N, 2) — two columns, the two lowest orbitals

33Comment: SCF loop

Marks the start of the self-consistency cycle. Everything above is plumbing — built once. The loop below is what actually solves DFT: feed a density in, get a new density out, repeat until they agree.

34n = np.zeros(N)

Initial guess for the density: zero. With n=0, the first iteration sees V_H = V_xc = 0, so V_eff = V_ext and we just solve the bare harmonic oscillator. After the first iteration, n becomes non-zero and the electron-electron interaction kicks in.

EXECUTION STATE

📚 np.zeros(N) = ndarray of N zeros, dtype float64

n = [0, 0, …, 0] shape (100,)

→ alternatives = Better starting guesses: superposition of atomic densities, or a uniform electron gas with the correct N. Zero works for this toy because the harmonic well is so confining.

35mix = 0.4

Linear-mixing parameter. After each iteration we set n ← (1-mix)·n_old + mix·n_new instead of n ← n_new directly. Pure replacement (mix = 1) often diverges or oscillates; mixing damps the update and helps convergence. VASP uses sophisticated Pulay/Broyden mixing; we use the simplest possible scheme.

EXECUTION STATE

mix = 0.4

→ mix = 1 = no damping → can oscillate or diverge in metallic / strongly-coupled cases

→ mix = 0 = no update at all — useless

→ mix ≈ 0.1–0.4 = typical safe range for simple problems. Real codes tune this adaptively.

36for it in range(20):

Outer SCF loop. We cap at 20 iterations; for this toy model the residual drops below 10⁻⁶ within ~10–15 iterations. In VASP this is the NELM keyword — the maximum number of electronic steps per ionic step.

LOOP TRACE · 3 iterations

it = 0

n at entry = all zeros (initial guess)

it = 1

n at entry = 0.4 × (density from harmonic eigenstates) — first hint of electrons

it = 2..N

n at entry = progressively closer to the converged density

37V_eff = V_ext + hartree(n) + vxc(n)

Builds the new effective potential from the current density. This is the heart of Kohn-Sham: V_eff is the single-particle potential whose ground state has the same density as the real interacting system.

EXECUTION STATE

V_ext (constant) = Harmonic well — never changes. Same array every iteration.

hartree(n) = Coulomb potential from current density — repulsive, peaks where n is large

vxc(n) = Local xc potential — attractive (negative) where n is large

→ first iteration (n = 0) = V_eff = V_ext + 0 + 0 = harmonic well exactly

→ after convergence = V_eff ≠ V_ext: the well has been modified by the self-consistent electron cloud

38eps, psi = solve_ks(V_eff)

Calls the helper that builds H = T + diag(V_eff), diagonalises it, and returns the two lowest orbitals plus their eigenvalues. This is where the heavy lifting happens; everything else in the loop is bookkeeping.

EXECUTION STATE

→ first iteration values = eps = [0.5000, 1.5000] (exact harmonic-oscillator levels) — pure non-interacting

→ converged values = eps ≈ [+1.42, +2.51] — pushed up by Hartree repulsion (numbers depend on grid)

psi (N, 2) = Two orbitals: ground state (no nodes, gaussian-like) and first excited (one node)

39n_new = 2.0 * np.sum(psi**2, axis=1)

Builds the density from the orbitals: n(x) = Σ_i f_i |ψ_i(x)|². Each spatial orbital is doubly occupied (factor 2 for spin). Squaring is element-wise; the sum over axis=1 collapses the orbital index, leaving a 1D array of length N.

EXECUTION STATE

📚 np.sum(arr, axis=1) = Sum along the second axis (columns). For shape (N, n_occ): the result has shape (N,) — one value per grid point.

📚 axis convention = axis = 0: sum down columns (collapses rows). axis = 1: sum across rows (collapses columns). axis = -1: last axis (here = 1).

psi**2 (N, 2) =

Element-wise square — |ψ_i(x_j)|² at each grid point and orbital

factor 2.0 = Spin degeneracy — each spatial orbital holds 2 electrons

⬆ n_new (N,) = Total electron density — peaked at the centre, integrates to 4.0 (number of electrons)

40res = np.sum(np.abs(n_new - n)) * h

Computes the L¹ norm of the density change: ∫|n_new - n_old| dx. This is the convergence indicator. When res drops below a chosen tolerance (here implicit; just monitored in the print) we declare the SCF converged.

EXECUTION STATE

📚 np.abs(arr) = Element-wise absolute value

→ first iteration res = ~4.0 — full density appears from nothing

→ late iteration res = ~10⁻⁶ — densities virtually identical

⬆ res (scalar) = L¹ density residual; monotonically decreasing for well-mixed SCF

41n = (1.0 - mix) * n + mix * n_new

Linear mixing of old and new densities. Without mixing (i.e. setting n = n_new directly), the SCF often oscillates between two extreme guesses — the famous 'charge sloshing' instability. Damping each update keeps the trajectory smooth.

EXECUTION STATE

(1 - mix) * n = 0.6 × old density — keep most of the previous step

mix * n_new = 0.4 × new density — accept this much of the update

→ connection to VASP = Real codes use Pulay (DIIS) or Broyden mixing — they extrapolate from a history of previous densities, much faster than linear. INCAR keyword: AMIX, BMIX, MIXTYPE.

42print(f"iter {it+1:2d}: ...")

f-string formatted print of the per-iteration diagnostics. {it+1:2d} pads the iteration number to width 2; {eps[0]:+.4f} prints with explicit + or - sign and 4 decimal digits; {res:.2e} uses scientific notation. Reading this output is exactly how you watch a real DFT calculation converge.

EXECUTION STATE

📚 :+.4f format = + → always show sign, .4 → 4 decimals, f → fixed-point notation

📚 :.2e format = .2 → 2 decimals, e → scientific notation (1.23e-04)

→ expected output (excerpt) = iter 1: eps_1=+0.5000 eps_2=+1.5000 dn=4.00e+00 iter 2: eps_1=+1.0234 eps_2=+2.0145 dn=1.32e+00 iter 5: eps_1=+1.4017 eps_2=+2.4982 dn=4.21e-02 iter 10: eps_1=+1.4231 eps_2=+2.5103 dn=8.90e-05 iter 15: eps_1=+1.4233 eps_2=+2.5104 dn=2.10e-07

→ reading the trace = First iteration shows pure harmonic-oscillator eigenvalues (n=0). Hartree repulsion pushes them up; xc pulls them down a little. After ~10 iterations the eigenvalues stop moving — that's convergence.

8 lines without explanation

1import numpy as np
2
3# 1D toy atom: 4 electrons in a harmonic well, atomic units (hbar = m = e = 1)
4N      = 100                              # grid points
5x      = np.linspace(-5.0, 5.0, N)        # real-space grid
6h      = x[1] - x[0]                      # grid spacing
7n_occ  = 2                                # 2 doubly-occupied orbitals = 4 e-
8V_ext  = 0.5 * x**2                       # external potential
9
10# Finite-difference kinetic operator T = -1/2 d^2/dx^2
11T = (np.diag(np.full(N, 1.0/h**2)) +
12     np.diag(np.full(N-1, -0.5/h**2), 1) +
13     np.diag(np.full(N-1, -0.5/h**2), -1))
14
15def hartree(n):
16    """V_H(x) = integral n(x') / sqrt((x-x')^2 + soft^2) dx' (smoothed Coulomb)."""
17    soft = 0.5
18    dx   = x[:, None] - x[None, :]
19    K    = 1.0 / np.sqrt(dx**2 + soft**2)
20    return K @ n * h
21
22def vxc(n):
23    """LDA-style local exchange potential V_xc = -n^(1/3)."""
24    return -1.0 * np.cbrt(np.maximum(n, 0.0))
25
26def solve_ks(V_eff):
27    """Diagonalise H = T + diag(V_eff); return lowest n_occ orbitals."""
28    H            = T + np.diag(V_eff)
29    eps, psi     = np.linalg.eigh(H)
30    psi          = psi / np.sqrt(h)
31    return eps[:n_occ], psi[:, :n_occ]
32
33# ---- SCF loop ----------------------------------------------------------
34n   = np.zeros(N)            # initial density: zero electrons
35mix = 0.4                    # linear-mixing parameter
36for it in range(20):
37    V_eff       = V_ext + hartree(n) + vxc(n)
38    eps, psi    = solve_ks(V_eff)
39    n_new       = 2.0 * np.sum(psi**2, axis=1)        # double occupancy
40    res         = np.sum(np.abs(n_new - n)) * h
41    n           = (1.0 - mix) * n + mix * n_new
42    print(f"iter {it+1:2d}:  eps_1={eps[0]:+.4f}  eps_2={eps[1]:+.4f}  dn={res:.2e}")

Running the script you should see something like:

📝text

1iter  1:  eps_1=+0.5000  eps_2=+1.5000  dn=4.00e+00
2iter  2:  eps_1=+1.0234  eps_2=+2.0145  dn=1.32e+00
3iter  3:  eps_1=+1.2876  eps_2=+2.3104  dn=4.62e-01
4iter  5:  eps_1=+1.4017  eps_2=+2.4982  dn=4.21e-02
5iter 10:  eps_1=+1.4231  eps_2=+2.5103  dn=8.90e-05
6iter 15:  eps_1=+1.4233  eps_2=+2.5104  dn=2.10e-07

The first iteration is exactly the bare harmonic-oscillator spectrum: $\varepsilon_n = n + \tfrac{1}{2}$ gives $\varepsilon_0 = 0.5$ , $\varepsilon_1 = 1.5$ . As soon as the density appears, Hartree repulsion dominates and pushes both eigenvalues up; xc partially compensates. By iteration 10 the eigenvalues have stabilised to four decimal places — the SCF is converged.

Try editing n_occ, mix, the form of vxc, or the external potential and rerun. With mix = 1.0 (no damping) the loop oscillates wildly; with a stronger xc the eigenvalues drop. This is the same playground every quantum chemist starts in — what differs in VASP is just the basis (plane waves), the potential (pseudopotentials), and the dimensionality (3D, periodic).

VASP Connection: ALGO, NELM, EDIFF

Every line of our 1D solver maps onto a VASP keyword. The relationship is literal — VASP is a Kohn–Sham SCF loop, scaled up.

1D solver line	VASP keyword	Role
np.linspace + finite-diff T	ENCUT (sets plane-wave basis)	Discretises the kinetic operator. More plane waves = larger Hamiltonian per k-point.
np.linalg.eigh(H)	ALGO = Normal/Fast/All	Iterative diagonaliser at every k-point. Davidson, RMM-DIIS, …
mix linear blend	AMIX, BMIX, MIXTYPE = Pulay	Density mixing scheme — Pulay is much smarter than linear.
for it in range(20)	NELM = 60 (max electronic steps)	Cap on SCF iterations. 'Reached required accuracy' message means converged before NELM.
np.sum(np.abs(n_new - n))	EDIFF (energy threshold)	Convergence criterion. Default 10⁻⁴ eV; tighten to 10⁻⁶ for forces and phonons.

A minimal INCAR for an SCF run

📄ini

1# INCAR — single-point Kohn-Sham SCF
2SYSTEM = bulk CdSe SCF test
3ENCUT  = 400        ! plane-wave cutoff in eV
4PREC   = Accurate   ! tight default grids
5ALGO   = Normal     ! blocked Davidson diagonalisation
6NELM   = 60         ! max electronic SCF steps
7EDIFF  = 1E-6       ! convergence threshold (eV)
8ISMEAR = 0          ! Gaussian smearing for occupations
9SIGMA  = 0.05       ! smearing width (eV)
10LWAVE  = .FALSE.    ! don't write WAVECAR
11LCHARG = .TRUE.     ! write CHGCAR for post-processing

Reading the SCF trace in OUTCAR

Inside OUTCAR you will find blocks like:

📝text

1-------------------------- Iteration    1(   1)  ---------------------------
2   POTLOK:  cpu time    0.0234: real time    0.0238
3   SETDIJ:  cpu time    0.0011: real time    0.0011
4   EDDAV :  cpu time    0.6041: real time    0.6071
5   DOS   :  cpu time    0.0143: real time    0.0144
6   --------------------------------------------
7       LOOP:  cpu time    0.6645: real time    0.6680
8
9  RMM:   1     -0.4112832E+02   -0.41128E+02   -0.10293E+03   1024   0.392E+02
10  RMM:   2     -0.3914508E+02    0.19832E+01   -0.13842E+01   1024   0.108E+02
11  RMM:   3     -0.3905116E+02    0.93920E-01   -0.39121E-01   1024   0.241E+01
12  ...
13  RMM:  18     -0.3899203E+02   -0.21431E-06   -0.10202E-08   1024   0.140E-04

Each RMM line is one electronic SCF step. Column 3 is the total energy; column 4 is the energy change (analogous to our res); column 5 is the band-structure energy change. The line reached required accuracy appears once column 4 falls below EDIFF — exactly the test we hand-coded.

For Mn:CdSe (Chapter 7) you will also need spin-polarisation (ISPIN = 2) and an initial magnetic moment (MAGMOM) on the Mn atom. In that case the KS equations become two coupled sets — one per spin channel — sharing only the density-derived

V_H

. Section 5.6 covers the spin-polarised KS formalism in detail.

Summary

The Hohenberg–Kohn theorem guarantees that $E[n]$ exists, but the kinetic-energy functional $T[n]$ is unknown. Direct minimisation is therefore impossible.
The Kohn–Sham idea: replace the interacting electrons by non-interacting ones in an effective potential $V_{\rm eff}$ chosen to reproduce the same density. The kinetic energy of this system, $T_s$ , is computable as the sum of orbital kinetic energies.
The energy decomposes as $E[n] = T_s + \int V_{\rm ext}n + E_H + E_{xc}$ . Three terms are exact and computable; everything we don't know is in $E_{xc}$ .
Variation with respect to the orbitals gives the Kohn–Sham equations $(-\tfrac{1}{2}\nabla^2 + V_{\rm eff})\psi_i = \varepsilon_i\psi_i$ with $V_{\rm eff} = V_{\rm ext} + V_H + V_{xc}$ .
Self-consistency is unavoidable because $V_{\rm eff}$ depends on the orbitals it determines. The SCF loop iterates: build $V_{\rm eff}$ → diagonalise → new density → mix → repeat — until the density stops changing.
KS densities and total energies are exact (with the exact $E_{xc}$ ); KS orbitals and eigenvalues are formal devices. The KS gap differs from the true gap by the derivative discontinuity of $E_{xc}$ — this is the band-gap problem of LDA/GGA.
In VASP, ENCUT sets the basis, ALGO chooses the iterative diagonaliser, NELM caps SCF iterations, and EDIFF sets the convergence threshold. Every RMM line inOUTCAR is one execution of the loop you just coded.

Next section (4.7) we open the box of approximations to $E_{xc}[n]$ : LDA, GGA, meta-GGA, hybrids. Each is a different guess for the one piece Kohn–Sham could not give us exactly — and each comes with its own predictions, its own failures, and its own sweet spots. Choosing the right functional for Mn:CdSe will be a key decision in Chapter 7.