Density Functional Theory and the Kohn-Sham Equations

A carbon atom sitting inside a heme cofactor is surrounded by roughly a dozen electrons. Their motions are correlated in complicated ways. Yet everything a chemist actually needs — the bond lengths, the spin state, the barrier to O₂ binding — is written into a single scalar field: the electron density ρ(r)\rho(\mathbf{r}), a function of only three spatial coordinates. Density Functional Theory (DFT) takes this observation and turns it into a rigorous, practical theory. Instead of chasing a many-electron wavefunction Ψ(r1,r2,…,rN)\Psi(\mathbf{r}_1,\mathbf{r}_2,\dots,\mathbf{r}_N) that lives in a 3N3N-dimensional space, one works with ρ(r)\rho(\mathbf{r}), a three-dimensional object. This object contains, in principle, everything there is to know about the ground state. This is not a hand-wave. The Hohenberg–Kohn theorems establish that the ground-state density determines the external potential (up to a constant). It therefore determines the Hamiltonian. It therefore determines every ground-state property. This includes energy, dipole moment, and force constants. All of these follow through a universal functional of the density alone. The catch, discovered immediately after the theorems were proved, is that nobody knows this universal functional in closed form. Kohn and Sham's resolution is the single most consequential trick in modern computational chemistry. It replaces the intractable interacting problem with an auxiliary system of noninteracting electrons. This auxiliary system is built, by construction, to reproduce the exact density. All the hard many-body physics is then hidden in one unknown piece called the exchange-correlation functional. Why does this matter so much? Because it converts an exponentially hard problem into something manageable. Wavefunction methods scale steeply with electron number. They struggle badly past a few dozen atoms. DFT instead produces a set of effective one-electron equations. These scale far more gently, typically as N3N^3 with system size. This is the machinery behind essentially every DFT calculation run today. It spans geometry optimizations of drug-like molecules to band-structure calculations of battery cathodes. The method inherits the same self-consistent iterative structure familiar from the The Self-Consistent-Field Procedure used in Hartree–Fock theory. But the orbitals it produces are not approximations to a many-electron wavefunction. They are exact-in-principle objects of a fictitious noninteracting reference system. This fictitious system happens to share the real system's density. This article builds the theory from its two load-bearing theorems. It proceeds through the constructive Kohn–Sham mapping. It ends at the self-consistent Kohn–Sham equations actually solved on a computer. It stays agnostic about the specific approximate exchange-correlation functional. That choice, and its consequences, is the subject of its own treatment. Instead, this article nails down exactly what DFT promises. It also clarifies what DFT does not promise. Finally, it explains why the Kohn–Sham reformulation was necessary at all.

Key takeaways
  • The ground-state electron density ρ(r)\rho(\mathbf{r}) determines the external potential and hence the full Hamiltonian and all ground-state properties (Hohenberg–Kohn Theorem 1).
  • A universal energy functional E[ρ]E[\rho] exists whose minimum, over densities integrating to NN electrons, equals the true ground-state energy (Hohenberg–Kohn Theorem 2, a variational principle in ρ\rho).
  • The exact functional's kinetic energy piece is unknown for interacting electrons, so Kohn and Sham introduce a fictitious noninteracting system with the same density, whose kinetic energy is exactly computable from orbitals.
  • All the missing many-body physics — the difference between the true kinetic energy and the noninteracting one, plus all non-classical electron-electron effects — is swept into one term, the exchange-correlation functional Exc[ρ]E_{xc}[\rho].
  • Minimizing the Kohn–Sham energy expression leads to one-electron Schrödinger-like equations with an effective potential that must be solved self-consistently, exactly as in Hartree–Fock.
  • DFT is formally exact but practically approximate: everything hinges on the (unknown) form of Exc[ρ]E_{xc}[\rho], whose construction is a separate, ongoing enterprise.

Key idea: The 3N-dimensional wavefunction problem can be replaced, without approximation in principle, by a 3-dimensional density problem, at the cost of hiding all the hard many-body physics inside one unknown functional.

Where this lives today

Kohn–Sham DFT is the workhorse of computational chemistry and materials science precisely because its cost scales gently with system size while retaining a formally exact foundation. Pharmaceutical computational chemists use it routinely to compute binding conformations, tautomer stabilities, and reaction barriers for drug candidates docked into protein active sites, often paired with implicit solvation models. In materials science, plane-wave Kohn–Sham DFT calculations predict formation energies, band gaps, and defect chemistry in battery cathode materials, photovoltaic absorbers, and catalytic surfaces — the entire discipline of high-throughput computational materials screening (materials genome-style databases) rests on running Kohn–Sham calculations across thousands of candidate compounds. In biology, QM/MM schemes embed a DFT-treated active site (an enzyme's metal cofactor, say) within a classical force-field-treated protein environment, allowing mechanistic study of enzymatic catalysis at a cost ordinary correlated wavefunction methods could not touch. Spectroscopists use time-dependent extensions of Kohn–Sham theory to predict UV-Vis absorption spectra, and instrumentation groups use DFT-computed vibrational frequencies to assign infrared and Raman spectra of newly synthesized compounds. The choice of which approximate exchange-correlation functional to trust for a given problem — a decision that can change predicted barrier heights by tens of kJ/mol — is itself a major applied subfield, covered under Exchange-Correlation Functionals and weighed against wavefunction-based alternatives in Choosing a Computational Method.

:::quiz[In the Kohn-Sham scheme, what is T_s[ρ] designed to represent?]

  • The exact kinetic energy of the interacting electron system
  • The classical Coulomb repulsion energy of the density with itself
  • The exactly computable kinetic energy of a fictitious noninteracting system with the same density ✓
  • The total exchange-correlation energy

:::quiz[A converged closed-shell Kohn-Sham calculation gives orbital energies ε₁ = -1.10 and ε₂ = -0.60 hartree (each doubly occupied), J[ρ] = 0.55 hartree, E_xc = -0.30 hartree, and ∫v_xc ρ dr = -0.34 hartree. What is the total electronic energy?]

  • -3.40 hartree
  • -2.71 hartree ✓
  • -3.02 hartree
  • -0.85 hartree

Ask the tutor: How does the choice of approximate exchange-correlation functional (LDA vs. GGA vs. hybrid) change the accuracy of Kohn–Sham predictions for reaction barriers?

The physical picture

Think first about what a wavefunction actually costs to store. For a molecule with NN electrons, Ψ\Psi depends on 3N3N spatial coordinates (plus spin), so doubling the number of electrons squares — no, exponentiates — the size of the space you must describe. The electron density, by contrast, is defined as

ρ(r)=N∫∣Ψ(r,r2,…,rN)∣2 dr2⋯drN\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

and no matter how many electrons a system has, ρ(r)\rho(\mathbf{r}) is always a function of the same three variables: it tells you the probability of finding an electron (any electron, since they're indistinguishable) at point r\mathbf{r}. Its integral over all space is NN, the total electron count, and it is directly observable — X-ray diffraction literally measures it. The astonishing claim of Hohenberg–Kohn theory is that this humble three-dimensional function, in its ground-state form, contains the same information as the full wavefunction. The route from "density" to "energy" runs through the variational principle familiar from ordinary quantum mechanics, but recast entirely in terms of ρ\rho. In the interacting real system, the Hamiltonian has three pieces: kinetic energy of the electrons, electron–electron repulsion, and the interaction of electrons with an external potential v(r)v(\mathbf{r}) (usually the nuclear attraction). Hohenberg and Kohn show that the kinetic and electron-repulsion pieces, summed into a single object F[ρ]F[\rho], are universal — they take the same functional form regardless of what molecule or solid you are studying — while all system-specific information enters only through v(r)v(\mathbf{r}). This split is the conceptual heart of the theory: chemistry-specific information (which atoms, where) lives in v(r)v(\mathbf{r}), and everything else about how electrons kinetically move and repel each other is captured once and for all in F[ρ]F[\rho]. The trouble is that F[ρ]F[\rho], while universal, is not known. In particular, nobody has an exact expression for the kinetic energy of an interacting electron gas as a functional of ρ\rho alone. Early attempts (Thomas–Fermi theory) approximated the kinetic energy directly from the density using a local electron-gas formula, and the results were qualitatively poor — Thomas–Fermi theory famously predicts no chemical bonding at all. Kohn and Sham's 1965 insight sidesteps this failure by refusing to approximate the kinetic energy directly. Instead, they introduce a fictitious system of NN noninteracting electrons, moving in some effective potential vs(r)v_s(\mathbf{r}), chosen so that this fictitious system's ground-state density exactly equals the true interacting density. Because the fictitious electrons are noninteracting, their many-electron wavefunction is exactly a single Slater determinant of one-electron orbitals {ϕi}\{\phi_i\} — and the kinetic energy of a Slater determinant is exactly computable, no approximation needed:

Ts[ρ]=∑i=1N⟨ϕi∣−12∇2∣ϕi⟩T_s[\rho] = \sum_{i=1}^N \left\langle \phi_i \left| -\tfrac{1}{2}\nabla^2 \right| \phi_i \right\rangle

This TsT_s is not the true kinetic energy of the interacting system (it's smaller, generally), but it captures the overwhelming majority of it. The remainder — plus everything about how electron repulsion differs from its naive classical (Coulomb, or "Hartree") approximation — gets packaged into the exchange-correlation functional Exc[ρ]E_{xc}[\rho]. This is where all of the theory's approximation is concentrated, and it is the price paid for making the kinetic energy tractable.

Figure
From interacting electrons to a fictitious noninteracting reference with the same density
Real interacting system Hamiltonian: T + V_ee + V_ext wavefunction Ψ(r₁,…,r_N) density ρ(r) — same in both! Kohn–Sham mapping Fictitious noninteracting system effective potential v_s(r) orbitals φ₁,…,φ_N (Slater determinant) exact kinetic energy T_s[ρ] same density ρ(r) as real system All the missing physics (T − T_s, and non-classical V_ee) packed into E_xc[ρ]

The self-consistent field structure appears the moment you ask what vs(r)v_s(\mathbf{r}) must be. Kohn and Sham show that if you set

vs(r)=vext(r)+∫ρ(r′)∣r−r′∣ dr′+vxc(r)v_s(\mathbf{r}) = v_{\text{ext}}(\mathbf{r}) + \int \frac{\rho(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}\,d\mathbf{r}' + v_{xc}(\mathbf{r})

then the noninteracting system built from vsv_s has exactly the true ground-state density. But vsv_s depends on ρ\rho through the Hartree (classical Coulomb) term and through vxc=δExc/δρv_{xc} = \delta E_{xc}/\delta \rho, and ρ\rho in turn is built from the orbitals that solve the Schrödinger equation with potential vsv_s. This circularity is resolved exactly the way it is in The Self-Consistent-Field Procedure: guess a density, build vsv_s, solve for orbitals, rebuild the density, and iterate to convergence. The output orbitals are not "approximate many-electron wavefunctions" in the Hartree–Fock sense; they are exact objects of an auxiliary noninteracting problem, constructed purely as a computational device to get the correct interacting density and, from it, the correct energy.

Figure
Self-consistent Kohn–Sham loop
Guess ρ⁽⁰⁾(r) Build v_s[ρ] = v_ext + v_H + v_xc Solve KS equations for φᵢ, εᵢ New ρ(r) = Σ|φᵢ|² Converged when ρ_new ≈ ρ_old → total energy E[ρ]

The core ideas

Hohenberg–Kohn Theorem 1 (existence/uniqueness). For a system of interacting electrons in an external potential vext(r)v_{\text{ext}}(\mathbf{r}), the ground-state density ρ(r)\rho(\mathbf{r}) determines vext(r)v_{\text{ext}}(\mathbf{r}) uniquely, up to an additive constant.

ρ(r)  ⟺  vext(r)  ⟹  H^  ⟹  Ψ0  ⟹  all ground-state properties\rho(\mathbf{r}) \;\Longleftrightarrow\; v_{\text{ext}}(\mathbf{r}) \;\Longrightarrow\; \hat{H} \;\Longrightarrow\; \Psi_0 \;\Longrightarrow\; \text{all ground-state properties}

where: ρ(r)\rho(\mathbf{r}) is the ground-state electron density; vext(r)v_{\text{ext}}(\mathbf{r}) is the potential electrons feel from the nuclei (or any external field); H^\hat{H} is the full electronic Hamiltonian, fixed once vextv_{\text{ext}} and NN are fixed; Ψ0\Psi_0 is the resulting nondegenerate ground-state wavefunction. Reading: the density is not merely a byproduct of the wavefunction — it is informationally equivalent to it, for ground states. Two different potentials (differing by more than a constant) cannot produce the same ground-state density. This is why ρ\rho can serve as the fundamental variable rather than Ψ\Psi.

Hohenberg–Kohn Theorem 2 (variational principle). There exists a universal functional F[ρ]F[\rho], independent of vextv_{\text{ext}}, such that the total energy functional

E[ρ]=F[ρ]+∫vext(r) ρ(r) drE[\rho] = F[\rho] + \int v_{\text{ext}}(\mathbf{r})\,\rho(\mathbf{r})\,d\mathbf{r}

is minimized, over all NN-representable trial densities, by the true ground-state density, and the minimum value equals the true ground-state energy E0E_0. where: F[ρ]=T[ρ]+Vee[ρ]F[\rho] = T[\rho] + V_{ee}[\rho] bundles the (unknown) universal kinetic and electron–electron repulsion functionals; the second term is the only piece that depends on which molecule you're studying. Reading: this converts electronic structure into a minimization problem over a 3-D function rather than a diagonalization or energy-minimization over a 3N3N-D wavefunction. It is exact — if F[ρ]F[\rho] were known exactly, this would be an exact theory, not an approximation.

Kohn–Sham decomposition. Introduce the noninteracting kinetic energy Ts[ρ]T_s[\rho] and the classical Hartree repulsion J[ρ]J[\rho], and define everything left over as Exc[ρ]E_{xc}[\rho]:

E[ρ]=Ts[ρ]+J[ρ]+Exc[ρ]+∫vext(r)ρ(r) drE[\rho] = T_s[\rho] + J[\rho] + E_{xc}[\rho] + \int v_{\text{ext}}(\mathbf{r})\rho(\mathbf{r})\,d\mathbf{r}

where:

Ts[ρ]=∑i=1N⟨ϕi∣−12∇2∣ϕi⟩,J[ρ]=12∬ρ(r)ρ(r′)∣r−r′∣ dr dr′T_s[\rho] = \sum_{i=1}^{N}\left\langle \phi_i\left|-\tfrac{1}{2}\nabla^2\right|\phi_i\right\rangle, \qquad J[\rho] = \frac{1}{2}\iint \frac{\rho(\mathbf{r})\rho(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}\,d\mathbf{r}\,d\mathbf{r}'
Exc[ρ]≡(T[ρ]−Ts[ρ])+(Vee[ρ]−J[ρ])E_{xc}[\rho] \equiv \big(T[\rho]-T_s[\rho]\big) + \big(V_{ee}[\rho]-J[\rho]\big)

Reading: ExcE_{xc} is defined as a residual — literally "whatever is left after you've done the easy parts exactly." It contains the true kinetic-energy correction (correlation kinetic energy) and all non-classical electron–electron effects: exchange (Pauli antisymmetry) and Coulomb correlation, connecting directly to the physics discussed under Electron Correlation and the Correlation Energy.

Kohn–Sham one-electron equations. Minimizing E[ρ]E[\rho] subject to orbital orthonormality, with ρ(r)=∑i∣ϕi(r)∣2\rho(\mathbf{r}) = \sum_i |\phi_i(\mathbf{r})|^2, gives

[−12∇2+vext(r)+∫ρ(r′)∣r−r′∣ dr′+vxc(r)]ϕi(r)=εi ϕi(r)\left[-\tfrac{1}{2}\nabla^2 + v_{\text{ext}}(\mathbf{r}) + \int\frac{\rho(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}\,d\mathbf{r}' + v_{xc}(\mathbf{r})\right]\phi_i(\mathbf{r}) = \varepsilon_i\,\phi_i(\mathbf{r})

where: vxc(r)=δExc[ρ]/δρ(r)v_{xc}(\mathbf{r}) = \delta E_{xc}[\rho]/\delta\rho(\mathbf{r}) is the exchange-correlation potential, a functional derivative; εi\varepsilon_i are Kohn–Sham orbital energies (not directly physical excitation energies, though the highest occupied one has a known meaning — see below); the bracketed operator is the effective one-electron Kohn–Sham Hamiltonian h^KS\hat{h}^{KS}. Reading: this looks exactly like a one-electron Schrödinger equation, and it is solved the same iterative way as the Hartree–Fock equations, but the physics inside vxcv_{xc} is, in principle, exact — it's just that nobody has the exact formula.

Derivation AFrom the Hohenberg–Kohn variational principle to the Kohn–Sham equations

Goal: show how demanding stationarity of E[ρ]E[\rho] under the Kohn–Sham decomposition, with the constraint that ρ\rho comes from NN orthonormal orbitals, produces the effective one-electron equations above.

Step 1 — Write the energy as an orbital functional. Substitute ρ(r)=∑i=1N∣ϕi(r)∣2\rho(\mathbf{r}) = \sum_{i=1}^N |\phi_i(\mathbf{r})|^2 into the Kohn–Sham energy expression:

E[{ϕi}]=∑i⟨ϕi∣−12∇2∣ϕi⟩+∫vextρ dr+J[ρ]+Exc[ρ]E[\{\phi_i\}] = \sum_i \left\langle \phi_i\left|-\tfrac{1}{2}\nabla^2\right|\phi_i\right\rangle + \int v_{\text{ext}}\rho\,d\mathbf{r} + J[\rho] + E_{xc}[\rho]

Every term is now expressed through the orbitals; only the first term needs orbitals explicitly (the rest depend on ρ\rho, itself built from orbitals).

Step 2 — Impose orthonormality with Lagrange multipliers. We minimize EE subject to ⟨ϕi∣ϕj⟩=δij\langle \phi_i|\phi_j\rangle = \delta_{ij}. Form the Lagrangian

L=E[{ϕi}]−∑i,jεij(⟨ϕi∣ϕj⟩−δij)L = E[\{\phi_i\}] - \sum_{i,j}\varepsilon_{ij}\left(\langle\phi_i|\phi_j\rangle - \delta_{ij}\right)

and require δL/δϕi∗(r)=0\delta L/\delta \phi_i^*(\mathbf{r}) = 0 for each ii. This is the same variational machinery used to derive the Hartree–Fock equations under The Self-Consistent-Field Procedure; the only difference is what sits inside the effective potential.

Step 3 — Differentiate each energy piece. The kinetic term gives −12∇2ϕi-\tfrac12\nabla^2\phi_i directly. The external and Hartree terms give, via the chain rule through ρ\rho,

δδϕi∗(r)∫vextρ dr=vext(r)ϕi(r),δJ[ρ]δϕi∗(r)=(∫ρ(r′)∣r−r′∣dr′)ϕi(r)\frac{\delta}{\delta\phi_i^*(\mathbf{r})}\int v_{\text{ext}}\rho\,d\mathbf{r} = v_{\text{ext}}(\mathbf{r})\phi_i(\mathbf{r}), \qquad \frac{\delta J[\rho]}{\delta\phi_i^*(\mathbf{r})} = \left(\int\frac{\rho(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}d\mathbf{r}'\right)\phi_i(\mathbf{r})

and by definition of the functional derivative,

δExc[ρ]δϕi∗(r)=δExcδρ(r)⋅δρ(r)δϕi∗(r)=vxc(r) ϕi(r)\frac{\delta E_{xc}[\rho]}{\delta \phi_i^*(\mathbf{r})} = \frac{\delta E_{xc}}{\delta \rho(\mathbf{r})}\cdot\frac{\delta\rho(\mathbf{r})}{\delta\phi_i^*(\mathbf{r})} = v_{xc}(\mathbf{r})\,\phi_i(\mathbf{r})

Step 4 — Collect terms and diagonalize the multipliers. Setting δL/δϕi∗=0\delta L/\delta\phi_i^* = 0 gives

[−12∇2+vext(r)+vH(r)+vxc(r)]ϕi(r)=∑jεijϕj(r)\left[-\tfrac12\nabla^2 + v_{\text{ext}}(\mathbf{r}) + v_H(\mathbf{r}) + v_{xc}(\mathbf{r})\right]\phi_i(\mathbf{r}) = \sum_j \varepsilon_{ij}\phi_j(\mathbf{r})

Because the effective Hamiltonian h^KS\hat h^{KS} is Hermitian, the Lagrange multiplier matrix εij\varepsilon_{ij} can be unitarily diagonalized (exactly as with Fock matrix off-diagonal multipliers), giving canonical orbitals that satisfy the single-index eigenvalue equation quoted above, with eigenvalues εi\varepsilon_i.

Step 5 — Interpret the self-consistency requirement. Because vHv_H and vxcv_{xc} both depend on ρ\rho, and ρ\rho depends on the ϕi\phi_i that solve this very equation, the equations must be solved iteratively: guess ρ\rho, build h^KS\hat h^{KS}, diagonalize to get new orbitals and a new ρ\rho, and repeat until the density (and hence the total energy) stops changing.

Go deeper · Derivation B · 1 minThe energy expression in terms of orbital energies, and the double-counting correction
Derivation BThe energy expression in terms of orbital energies, and the double-counting correction

Goal: show that the total Kohn–Sham energy is not simply the sum of occupied orbital energies, and derive the correction term explicitly — a frequent point of confusion analogous to the same issue in Hartree–Fock theory.

Step 1 — Sum the orbital energies. Multiply the Kohn–Sham equation for orbital ii on the left by ϕi∗\phi_i^*, integrate, and sum over occupied orbitals:

∑iεi=∑i⟨ϕi∣−12∇2∣ϕi⟩+∫vextρ dr+∫vH(r)ρ(r) dr+∫vxc(r)ρ(r) dr\sum_i \varepsilon_i = \sum_i \left\langle \phi_i\left|-\tfrac12\nabla^2\right|\phi_i\right\rangle + \int v_{\text{ext}}\rho\,d\mathbf{r} + \int v_H(\mathbf{r})\rho(\mathbf{r})\,d\mathbf{r} + \int v_{xc}(\mathbf{r})\rho(\mathbf{r})\,d\mathbf{r}

Step 2 — Compare to the true total energy. The actual Kohn–Sham total energy is

E=∑i⟨ϕi∣−12∇2∣ϕi⟩+∫vextρ dr+J[ρ]+Exc[ρ]E = \sum_i\left\langle \phi_i\left|-\tfrac12\nabla^2\right|\phi_i\right\rangle + \int v_{\text{ext}}\rho\,d\mathbf{r} + J[\rho] + E_{xc}[\rho]

Note that ∫vHρ dr=2J[ρ]\int v_H \rho\, d\mathbf{r} = 2J[\rho], because vH(r)=∫ρ(r′)/∣r−r′∣ dr′v_H(\mathbf{r}) = \int \rho(\mathbf{r}')/|\mathbf{r}-\mathbf{r}'|\,d\mathbf{r}' integrates against ρ\rho to give twice the Hartree energy (the same factor-of-two "double counting" familiar from the Hartree–Fock Coulomb-energy bookkeeping).

Step 3 — Subtract to isolate the correction. Comparing the orbital-energy sum to EE term by term:

E=∑iεi−J[ρ]+Exc[ρ]−∫vxc(r)ρ(r) drE = \sum_i \varepsilon_i - J[\rho] + E_{xc}[\rho] - \int v_{xc}(\mathbf{r})\rho(\mathbf{r})\,d\mathbf{r}

Step 4 — State the physical meaning. The correction −J[ρ]-J[\rho] removes the double-counted Hartree repulsion (each electron pair's repulsion was counted once in εi\varepsilon_i and once in εj\varepsilon_j), and Exc[ρ]−∫vxcρ drE_{xc}[\rho]-\int v_{xc}\rho\,d\mathbf{r} replaces the (incorrectly double-counted) exchange-correlation potential energy with the correct exchange-correlation energy functional value. This is structurally identical to the Hartree–Fock double-counting correction, generalized to include ExcE_{xc}.

Worked examples

Example 1 — Foundational: verifying the Hartree double-counting factor for a simple density. Consider a toy one-dimensional "density" used purely to fix the bookkeeping numerically (not a real atom, but numbers chosen off any special symmetry so the algebra is checkable): suppose we model ρ(r)\rho(\mathbf{r}) as a normalized two-Gaussian mixture contributing charges q1=1.3q_1 = 1.3 and q2=0.7q_2 = 0.7 (summing to N=2N=2) concentrated effectively at two points separated so that their mutual Coulomb integral evaluates to g12=0.42g_{12} = 0.42 hartree, and each has a self-repulsion integral of g11=0.90g_{11}=0.90 and g22=0.55g_{22}=0.55 hartree (in some model units). The classical Hartree energy is

J[ρ]=12(q12g11+q22g22+2q1q2g12)J[\rho] = \tfrac12\left(q_1^2 g_{11} + q_2^2 g_{22} + 2q_1q_2 g_{12}\right)

Plugging in:

q12g11=(1.3)2(0.90)=1.521,q22g22=(0.7)2(0.55)=0.2695q_1^2 g_{11} = (1.3)^2(0.90) = 1.521, \qquad q_2^2 g_{22} = (0.7)^2(0.55) = 0.2695
2q1q2g12=2(1.3)(0.7)(0.42)=0.76442q_1 q_2 g_{12} = 2(1.3)(0.7)(0.42) = 0.7644
J[ρ]=12(1.521+0.2695+0.7644)=12(2.5549)=1.27745 hartreeJ[\rho] = \tfrac12(1.521 + 0.2695 + 0.7644) = \tfrac12(2.5549) = 1.27745 \text{ hartree}

Now check the claimed identity ∫vHρ dr=2J[ρ]\int v_H \rho\,d\mathbf{r} = 2J[\rho]. The Hartree potential felt by charge 1 is vH(1)=q1g11+q2g12=1.3(0.90)+0.7(0.42)=1.17+0.294=1.464v_H^{(1)} = q_1 g_{11} + q_2 g_{12} = 1.3(0.90) + 0.7(0.42) = 1.17 + 0.294 = 1.464, and similarly vH(2)=q1g12+q2g22=1.3(0.42)+0.7(0.55)=0.546+0.385=0.931v_H^{(2)} = q_1 g_{12} + q_2 g_{22} = 1.3(0.42)+0.7(0.55) = 0.546+0.385=0.931. Then

∫vHρ dr=q1vH(1)+q2vH(2)=1.3(1.464)+0.7(0.931)=1.9032+0.6517=2.5549\int v_H\rho\,d\mathbf{r} = q_1 v_H^{(1)} + q_2 v_H^{(2)} = 1.3(1.464) + 0.7(0.931) = 1.9032 + 0.6517 = 2.5549

which is indeed exactly 2J[ρ]=2(1.27745)=2.55492J[\rho] = 2(1.27745) = 2.5549, confirming the factor-of-two relation used in the double-counting derivation, with genuinely unequal, non-round input numbers. Limiting case: if the two "electrons" merge into a single point charge of magnitude NN (or if q1=q2=1q_1=q_2=1 with g11=g22=g12=gg_{11}=g_{22}=g_{12}=g), the formula collapses to J=12(g+g+2g)=2gJ = \tfrac12(g+g+2g) = 2g and ∫vHρ=2(2g)=4g=2J\int v_H\rho = 2(2g)=4g = 2J — the identity survives the symmetric limit as it must, but the asymmetric numbers above are what actually stress-test the formula.

Example 2 — Applied: estimating the HOMO-LUMO-adjacent quantity for a diatomic using Koopmans-like reasoning in Kohn–Sham theory. Unlike Hartree–Fock, Kohn–Sham orbital energies do not generally approximate ionization energies via Koopmans' theorem — except for one exact result: the highest occupied Kohn–Sham orbital energy equals minus the exact ionization energy, provided the exact exchange-correlation functional is used,

εHOMO=−I\varepsilon_{\text{HOMO}} = -I

Suppose a calculation on a diatomic radical gives a HOMO Kohn–Sham eigenvalue of εHOMO=−0.3187\varepsilon_{\text{HOMO}} = -0.3187 hartree using some approximate functional. Converting to eV (1 hartree =27.2114=27.2114 eV):

Ipredicted=0.3187×27.2114=8.672 eVI_{\text{predicted}} = 0.3187 \times 27.2114 = 8.672 \text{ eV}

If the experimental ionization energy is Iexp=9.845I_{\text{exp}} = 9.845 eV, the discrepancy is

Δ=9.845−8.672=1.173 eV,relative error=1.1739.845×100%=11.9%\Delta = 9.845 - 8.672 = 1.173\text{ eV}, \qquad \text{relative error} = \frac{1.173}{9.845}\times100\% = 11.9\%

This sizable error is typical of common approximate functionals (which suffer from self-interaction error), even though the underlying theorem is exact for the true ExcE_{xc}. This is precisely why orbital energies from approximate Kohn–Sham calculations are treated as qualitative indicators rather than quantitative ionization energies — a caveat carried forward whenever Exchange-Correlation Functionals are compared against each other. Limiting case: as the approximate functional approaches the exact (unknown) functional, Δ→0\Delta \to 0 by the exactness of the εHOMO=−I\varepsilon_{\text{HOMO}}=-I theorem; the 11.9% error quantifies how far a typical practical functional sits from that exact limit for this system.

Example 3 — Exam/competition level: energy bookkeeping from orbital energies for a four-electron closed-shell system. A closed-shell four-electron system (two doubly occupied Kohn–Sham orbitals) is converged with orbital energies ε1=−1.204\varepsilon_1 = -1.204 hartree and ε2=−0.512\varepsilon_2 = -0.512 hartree. The converged Hartree energy is J[ρ]=0.687J[\rho] = 0.687 hartree, the exchange-correlation energy is Exc[ρ]=−0.398E_{xc}[\rho] = -0.398 hartree, and the exchange-correlation potential energy (the expectation value ∫vxcρ dr\int v_{xc}\rho\,d\mathbf{r}) is −0.451-0.451 hartree. Find the total electronic energy. Using the double-counting correction derived above,

E=∑iniεi−J[ρ]+Exc[ρ]−∫vxcρ drE = \sum_i n_i\varepsilon_i - J[\rho] + E_{xc}[\rho] - \int v_{xc}\rho\,d\mathbf{r}

where ni=2n_i=2 for each doubly-occupied orbital. First sum the orbital contribution:

∑iniεi=2(−1.204)+2(−0.512)=−2.408−1.024=−3.432 hartree\sum_i n_i \varepsilon_i = 2(-1.204) + 2(-0.512) = -2.408 - 1.024 = -3.432 \text{ hartree}

Now apply the correction:

E=−3.432−0.687+(−0.398)−(−0.451)E = -3.432 - 0.687 + (-0.398) - (-0.451)
E=−3.432−0.687−0.398+0.451E = -3.432 - 0.687 - 0.398 + 0.451
E=−4.066 hartreeE = -4.066 \text{ hartree}

So the total electronic energy is −4.066-4.066 hartree, noticeably different from the naive orbital-energy sum of −3.432-3.432 hartree — a difference of 0.6340.634 hartree (≈17.3\approx 17.3 eV), illustrating why summing Kohn–Sham eigenvalues alone (as one might naively do by analogy with a simple particle-in-a-box picture) badly misrepresents the total energy, exactly as it would in Hartree–Fock theory. Limiting case: if ExcE_{xc} and its potential term happened to satisfy Exc=∫vxcρ drE_{xc} = \int v_{xc}\rho\,d\mathbf{r} (true only for potentials linear in ρ\rho, not realistic ExcE_{xc} functionals), those two terms would cancel and EE would reduce to ∑niεi−J[ρ]\sum n_i\varepsilon_i - J[\rho], the pure Hartree-theory double-counting formula — showing the ExcE_{xc} correction is an additive refinement on top of the classical bookkeeping, not a replacement for it.

Limiting cases, scope & common traps

SituationWhat happensWhy
Exact Exc[ρ]E_{xc}[\rho] usedKohn–Sham theory gives the exact ground-state energy and densityHohenberg–Kohn theorems are exact; all approximation is isolated in ExcE_{xc}
N→1N\to 1 (one electron)Ts[ρ]T_s[\rho] becomes exact, and the exact ExcE_{xc} should cancel the Hartree self-repulsion J[ρ]J[\rho] entirelyA single electron cannot repel itself; approximate functionals that fail to do this exhibit self-interaction error
Noninteracting limit (Vee→0V_{ee}\to 0)Kohn–Sham reduces exactly to ordinary one-electron quantum mechanics; Exc→0E_{xc}\to 0The fictitious system is the real system when there's no interaction to hide
Near-degenerate ground statesHohenberg–Kohn theorem 1's uniqueness argument (originally for nondegenerate states) requires care; extensions (Levy constrained search) handle degeneracyEnsemble/generalized formulations restore rigor
Charge-transfer and Rydberg excited states (via response-theory extensions)Common approximate functionals perform poorlySelf-interaction error and incorrect asymptotic potential decay
Comparing orbital energies to spectroscopic energiesOnly εHOMO=−I\varepsilon_{\text{HOMO}} = -I is exact (with exact ExcE_{xc}); other εi\varepsilon_i have no rigorous spectroscopic meaning in ground-state DFTKohn–Sham orbitals are auxiliary constructs, not approximate quasiparticle wavefunctions

The chemist behind it

Origins · 1 minWalter Kohn and how this was discovered
WK
Walter Kohn
1923–2016 · Refugee physicist turned architect of density-functional theory, University of California San Diego and Santa Barbara

Walter Kohn was born in Vienna and, as a Jewish teenager, escaped Nazi-annexed Austria via the Kindertransport in 1939 while his parents perished in the Holocaust. He resettled in Canada, was briefly interned as an "enemy alien" during the Second World War, and went on to earn a PhD in physics at Harvard under Julian Schwinger. Kohn's early career centered on solid-state physics — pseudopotential theory, and the physics of electrons in semiconductors and metals — before he turned, with postdoctoral collaborator Pierre Hohenberg and later graduate student Lu Jeu Sham, to the foundational question of whether electron density alone could serve as the basic variable of quantum many-body theory. The 1964 Hohenberg–Kohn paper and the 1965 Kohn–Sham paper together transformed a formally elegant but practically inert existence theorem into a genuinely computable scheme, and Kohn spent much of his later career, at UC San Diego and then UC Santa Barbara, championing DFT's application across physics, chemistry, and materials science. He shared the 1998 Nobel Prize in Chemistry with John Pople, a striking recognition that a physicist's density-based reformulation of quantum mechanics had become indispensable to chemistry.

History
- 1927 — Thomas and (independently) Fermi propose an early density-based model of atomic electronic structure, with the kinetic energy approximated directly from a local electron-gas formula.
- 1964 — Pierre Hohenberg and Walter Kohn prove the two foundational theorems establishing the electron density as a legitimate fundamental variable for the ground state.
- 1965 — Walter Kohn and Lu Jeu Sham introduce the fictitious noninteracting reference system, yielding the practical, self-consistent Kohn–Sham equations.
- 1970s–1980s — Local density and generalized-gradient exchange-correlation functionals are developed, making DFT viable for real molecules and solids.
- 1998 — Walter Kohn shares the Nobel Prize in Chemistry with John Pople for the development of density-functional methods in quantum chemistry.
Is DFT exact, or is it inherently an approximate method?

The Hohenberg–Kohn theorems and the Kohn–Sham reformulation are exact statements: if you had the true universal functional F[ρ]F[\rho], or equivalently the true exchange-correlation functional Exc[ρ]E_{xc}[\rho], Kohn–Sham theory would reproduce the exact ground-state energy and density with no approximation at all. In that sense DFT is a formally exact theory, on the same conceptual footing as the Schrödinger equation itself. The approximation enters only because nobody has ever written down the exact Exc[ρ]E_{xc}[\rho] in closed form; every practical DFT calculation uses one of many approximate functionals, and it is those approximations — not the underlying theory — that introduce error. This distinction matters enormously when judging DFT's failures: they are failures of a specific approximate functional, not of the density-functional framework.

Why can't we just use the electron density directly with an approximate kinetic-energy functional, skipping Kohn-Sham orbitals entirely?

This is exactly what Thomas–Fermi theory (and its descendants, orbital-free DFT) attempts, and historically it failed badly for chemistry: the Thomas–Fermi kinetic-energy functional, built from a local electron-gas approximation, is not accurate enough to reproduce shell structure or chemical bonding, and famously predicts that molecules do not bind at all. The kinetic energy of even a noninteracting electron system, expressed purely as a functional of ρ\rho, has no known accurate closed form that avoids using orbitals. Kohn and Sham's move — computing Ts[ρ]T_s[\rho] exactly from a Slater determinant of auxiliary orbitals rather than approximating it directly from ρ\rho — is precisely what fixes this failure and makes DFT chemically accurate.

Do Kohn-Sham orbitals mean anything chemically, like frontier molecular orbitals in Hartree-Fock theory?

They are widely used that way in practice — HOMO/LUMO pictures, orbital symmetry arguments, and qualitative bonding discussions routinely borrow Kohn–Sham orbitals — and empirically they often resemble Hartree–Fock orbitals closely for many systems. But formally, Kohn–Sham orbitals are defined only as the machinery that reproduces the correct total density of the real interacting system; they are solutions of an auxiliary noninteracting problem, not approximations to pieces of the true many-electron wavefunction. Their qualitative use is generally reliable for rationalizing trends, but treating individual orbital energies as spectroscopic quantities (beyond the exact HOMO-ionization-energy theorem) is not formally justified.

How does the Kohn-Sham self-consistent field procedure differ numerically from the Hartree-Fock SCF procedure?

The iterative structure is essentially identical: guess a density or set of orbitals, build an effective one-electron Hamiltonian, diagonalize it to get new orbitals, rebuild the density, and repeat until convergence — the same self-consistent loop used throughout The Self-Consistent-Field Procedure. The difference lies entirely in what sits inside the effective potential. Hartree–Fock includes an exact, nonlocal exchange term derived from antisymmetrizing a single Slater determinant, with no correlation whatsoever. Kohn–Sham instead includes a (generally local or semi-local) exchange-correlation potential vxc(r)v_{xc}(\mathbf{r}) derived from an approximate density functional, which folds in both exchange and correlation effects simultaneously, at typically lower computational cost per iteration than exact Hartree–Fock exchange.

Why does the Hohenberg-Kohn framework matter if we still need approximate functionals to compute anything?

Because it tells you, with certainty, that a three-variable function contains all the ground-state information that a 3N3N-variable wavefunction does — a claim that is not obvious a priori and that licenses the entire enterprise of density-based methods. Without the Hohenberg–Kohn theorems, there would be no guarantee that minimizing some energy expression over densities gives anything physically meaningful; the theory could easily have been an ill-posed dead end, the way naive density-only guesses often are. The theorems supply the rigorous foundation; the ongoing search for better exchange-correlation functionals (see Exchange-Correlation Functionals) is the practical engineering built on that foundation.

Does the Kohn-Sham approach handle electron correlation exactly, in principle?

Yes, in the sense that the exact Exc[ρ]E_{xc}[\rho] would include all correlation effects — both the "dynamical" correlation missing from Hartree–Fock and any static/near-degeneracy correlation — exactly. The entire content of Electron Correlation and the Correlation Energy is, formally, hidden inside Exc[ρ]E_{xc}[\rho] along with exchange effects. In practice, most widely used approximate functionals handle dynamical correlation reasonably well but struggle with strongly correlated, near-degenerate situations (some transition-metal complexes, bond-breaking along a potential energy surface), which is why wavefunction-based multireference methods remain preferred for those cases, a comparison made explicit under Choosing a Computational Method.

What does N-representability mean, and why does it matter for the variational principle in Hohenberg-Kohn theory?

A trial density is called NN-representable if it can, in principle, arise from some valid NN-electron antisymmetric wavefunction (equivalently, if it is non-negative everywhere and integrates to NN). The Hohenberg–Kohn variational theorem states that minimizing E[ρ]E[\rho] over all NN-representable densities yields the true ground-state energy, with equality only at the true ground-state density. This matters because the minimization is only guaranteed to behave correctly within this restricted class of densities; feeding the functional a density that could never come from any physical wavefunction voids the guarantee, though in practice Kohn–Sham densities (built from an actual Slater determinant) are automatically NN-representable by construction.

Before this

Where this leads