The Kohn-Sham equations are what turned Density Functional Theory from an elegant existence proof into the computational workhorse of modern materials science. They provide a practical recipe for computing the exact ground-state density of an interacting electron system by solving a set of deceptively simple single-particle equations. This post unpacks how that trick works and what its outputs really mean.
Why We Need the Kohn-Sham Trick
The Hohenberg-Kohn theorems guarantee that the ground-state energy is a functional of the electron density \(n(\mathbf{r})\), and that minimizing this functional yields the true ground state. But they say nothing about how to write the functional — in particular, the kinetic energy of interacting electrons as a functional of density is unknown and notoriously difficult to approximate.
Early orbital-free approaches like Thomas-Fermi approximated the kinetic energy directly as a functional of \(n(\mathbf{r})\). The results were catastrophic for chemistry: Thomas-Fermi theory cannot even bind molecules, because the kinetic-energy approximation is too crude. The kinetic energy is simply too large and too sensitive to get wrong.
Kohn and Sham’s 1965 insight (Phys. Rev. 140, A1133) was to recover almost all of the kinetic energy exactly by reintroducing orbitals — but for a cleverly chosen fictitious system.
The Kohn-Sham Ansatz: A Non-Interacting Reference System
The core idea is to replace the real system of \(N\) interacting electrons with a fictitious system of \(N\) non-interacting electrons that is constructed to have exactly the same ground-state density \(n(\mathbf{r})\).
This is a bold assumption — it presumes that for any physical interacting density, there exists a non-interacting system reproducing it (the non-interacting v-representability assumption). For the vast majority of real systems this holds, and it buys us something invaluable: a non-interacting system’s wavefunction is a single Slater determinant of single-particle orbitals, and its kinetic energy can be computed exactly from those orbitals.
flowchart LR A[Interacting electrons] -->|same density n r| B[Non-interacting KS system] B --> C[Orbitals phi_i] C --> D[T_s exact from orbitals] C --> E[n r = sum of phi_i squared]
The density is built from the occupied orbitals \(\phi_i(\mathbf{r})\):
$$ n(\mathbf{r}) = \sum_{i \in \mathrm{occ}} \lvert\phi_i(\mathbf{r})\rvert^2 $$and the non-interacting kinetic energy is:
$$ T_s = \sum_{i \in \mathrm{occ}} \int \phi_i^*(\mathbf{r})\,\Bigl(-\tfrac{1}{2}\nabla^2\Bigr)\,\phi_i(\mathbf{r})\, d\mathbf{r} $$(atomic units, where \(\hbar = m_e = e = 1\)). This \(T_s\) is not the true interacting kinetic energy, but it captures the overwhelming majority of it. The small remainder is absorbed elsewhere.
Decomposing the Total Energy
Kohn and Sham partitioned the total energy functional as:
$$ E[n] = T_s[n] + E_{\mathrm{ext}}[n] + E_H[n] + E_{\mathrm{xc}}[n] $$Each term has a clear physical meaning:
- \(T_s[n]\) — kinetic energy of the non-interacting reference electrons.
- \(E_{\mathrm{ext}}[n]\) — interaction of the electrons with the external potential (nuclei), \(\int v_{\mathrm{ext}}(\mathbf{r})\, n(\mathbf{r})\, d\mathbf{r}\).
- \(E_H[n]\) — the classical Hartree Coulomb repulsion of the density with itself.
- \(E_{\mathrm{xc}}[n]\) — the exchange-correlation energy, defined as everything left over.
The genius is in the definition of \(E_{\mathrm{xc}}\). It is constructed to contain exactly two corrections: the difference between the true and non-interacting kinetic energies, and the non-classical part of the electron-electron interaction (exchange plus correlation):
$$ E_{\mathrm{xc}}[n] = \bigl(T[n] - T_s[n]\bigr) + \bigl(V_{\mathrm{ee}}[n] - E_H[n]\bigr) $$Because \(T - T_s\) is small, \(E_{\mathrm{xc}}\) is a small fraction of the total energy — which is precisely why approximating it works so well in practice.
Deriving the Kohn-Sham Equations
Minimizing \(E[n]\) subject to keeping the orbitals normalized, and requiring the density to match the non-interacting form, leads to a set of single-particle Schrödinger-like equations:
$$ \Bigl(-\tfrac{1}{2}\nabla^2 + v_{\mathrm{eff}}(\mathbf{r})\Bigr) \phi_i(\mathbf{r}) = \varepsilon_i\, \phi_i(\mathbf{r}) $$These are the Kohn-Sham equations. They look exactly like independent-particle equations, but the potential \(v_{\mathrm{eff}}(\mathbf{r})\) cleverly encodes all the many-body physics.
The Effective (Kohn-Sham) Potential
The effective potential is the sum of three contributions:
$$ v_{\mathrm{eff}}(\mathbf{r}) = v_{\mathrm{ext}}(\mathbf{r}) + v_H(\mathbf{r}) + v_{\mathrm{xc}}(\mathbf{r}) $$- \(v_{\mathrm{ext}}(\mathbf{r})\) — the external potential from the nuclei.
- \(v_H(\mathbf{r})\) — the Hartree potential, \(\int n(\mathbf{r}') / \lvert\mathbf{r} - \mathbf{r}'\rvert\, d\mathbf{r}'\), the electrostatic field of the electron cloud.
- \(v_{\mathrm{xc}}(\mathbf{r})\) — the exchange-correlation potential, defined as the functional derivative \(\delta E_{\mathrm{xc}} / \delta n(\mathbf{r})\).
The exchange-correlation potential is where all approximation enters. Choosing its form — LDA, GGA, meta-GGA, or hybrid — is the central modeling decision in any calculation, a topic we explore in our post on exchange-correlation functionals.
The Self-Consistent Field Loop
Here is the catch that makes the Kohn-Sham equations nonlinear: \(v_{\mathrm{eff}}(\mathbf{r})\) depends on \(n(\mathbf{r})\) (through \(v_H\) and \(v_{\mathrm{xc}}\)), but \(n(\mathbf{r})\) is built from the orbitals \(\phi_i\), which are themselves the solutions of the equations. The potential depends on its own answer.
The resolution is iteration to self-consistency:
flowchart TD
A[Guess initial density n r] --> B[Construct v_eff from density]
B --> C[Solve KS equations for phi_i and eps_i]
C --> D[Rebuild density from occupied orbitals]
D --> E[Mix new density with previous]
E --> F{Density or energy converged?}
F -->|No| B
F -->|Yes| G[Converged ground state]- Guess an initial density \(n(\mathbf{r})\) — often a superposition of atomic densities.
- Construct \(v_{\mathrm{eff}}(\mathbf{r})\) from the current density.
- Solve the Kohn-Sham equations to obtain new orbitals \(\phi_i\) and eigenvalues \(\varepsilon_i\).
- Rebuild the density from the occupied orbitals.
- Mix the new density with the previous one (linear, Broyden, or Pulay/DIIS mixing) to aid convergence.
- Check convergence: if the density or total energy is stable to within a threshold, stop; otherwise return to step 2.
Why Density Mixing Matters
Naively feeding the output density straight back as the next input almost always diverges or oscillates, especially in metals, due to charge sloshing. Mixing schemes damp these oscillations. Common thresholds and controls include:
| Parameter | Typical value | Role |
|---|---|---|
| Energy convergence | \(10^{-6}\) Ry | Stops SCF when total energy is stable |
| Density convergence | \(10^{-8}\) to \(10^{-10}\) | Tighter criterion on charge density |
| Mixing factor (\(\beta\)) | 0.1 to 0.7 | Fraction of new density blended in |
| Smearing (metals) | 0.01 to 0.02 Ry | Broadens occupations near the Fermi level |
In Quantum ESPRESSO these map to conv_thr, mixing_beta, and degauss in the input file. Poor SCF convergence is most often a symptom of too-aggressive mixing, insufficient smearing in a metal, or a bad starting geometry.
What Do the Kohn-Sham Orbitals and Eigenvalues Mean?
This is the most misunderstood aspect of the theory, and it deserves a careful statement.
The orbitals \(\phi_i\) are auxiliary mathematical objects. They belong to the fictitious non-interacting system, not to real electrons. Only two things about them are formally guaranteed:
- Their squared sum reproduces the exact ground-state density.
- Their total kinetic energy is \(T_s\).
The eigenvalues \(\varepsilon_i\) are not, in general, electron removal or addition energies. There is one rigorous exception: Janak’s theorem relates the highest occupied eigenvalue to minus the ionization potential for the exact functional. The others carry no strict physical meaning.
That said, in practice Kohn-Sham eigenvalues are routinely — and successfully — interpreted as approximate quasiparticle energies. Band structures plotted from \(\varepsilon_i(\mathbf{k})\) correctly predict whether a material is a metal or semiconductor and usually reproduce band shapes well.
The Band-Gap Problem
The most notorious consequence is that the Kohn-Sham gap systematically underestimates the true fundamental band gap. Two effects combine: the approximate functionals are imperfect, and the exact exchange-correlation potential has a derivative discontinuity that Kohn-Sham eigenvalue differences miss entirely. This is why silicon’s gap comes out near 0.6 eV with standard GGA versus its experimental 1.17 eV. Remedies include hybrid functionals and many-body GW corrections.
Practical Cost and Scaling
Solving the Kohn-Sham equations dominates the runtime of a DFT calculation. The formal bottleneck is orthogonalizing and diagonalizing orbitals, which scales as \(N^3\) with system size. Each SCF iteration also requires evaluating the Hartree potential, typically done efficiently with fast Fourier transforms in plane-wave codes like Quantum ESPRESSO.
For large supercells the number of SCF iterations and the cost per iteration both grow, which is exactly where GPU acceleration delivers the largest wins.
From Orbitals to Observable Properties
Once the SCF loop converges, the Kohn-Sham orbitals and density unlock a broad set of properties. The total energy is the foundation for everything downstream. Its derivatives with respect to atomic positions give the forces (via the Hellmann-Feynman theorem), which drive geometry relaxation and molecular dynamics; derivatives with respect to the cell give the stress tensor. Plotting the eigenvalues along a path in reciprocal space produces the band structure, while summing states within energy windows yields the density of states.
Second derivatives, computed with density functional perturbation theory (DFPT), give phonon frequencies, dielectric constants, and Born effective charges without the need for large finite-difference supercells. This is one of the most powerful features of the Kohn-Sham framework: a single converged ground state seeds an entire cascade of predictions.
Common Convergence Pitfalls
Even experienced practitioners hit SCF trouble. The most frequent causes are:
- Metallic systems without smearing — sharp occupation cutoffs at the Fermi level cause oscillations; a small Gaussian or Marzari-Vanderbilt smearing fixes this.
- Too-large a mixing factor — aggressive mixing amplifies charge sloshing in large cells; reduce \(\beta\).
- A poor initial guess — magnetic systems in particular need a sensible starting spin configuration to avoid falling into the wrong local minimum.
- Insufficient basis or k-points — an under-converged setup can appear to converge to the wrong answer.
Diagnosing these early saves enormous compute time, especially on large supercells where each iteration is expensive.
Run it on Simatra
Simatra runs these Kohn-Sham SCF calculations on GPU-accelerated clusters using the GPU-Opt-V2 instance type, achieving up to 5x faster convergence and supporting supercells up to ~2,000 atoms. Pick from two computational engines — the open-source Quantum ESPRESSO or Simatra’s native KRONOS engine, a C++20, GPL-3.0 code with CUDA, HIP, and Metal backends. Spin up your first self-consistent calculation today with a free trial and $100 in credits.