Density-Functional Theory Workflow

A density-functional calculation combines a physical approximation for the exchange-correlation energy with numerical representations and a nonlinear solve. A converged Kohn–Sham density removes the iteration error for that setup; it does not remove functional error, basis error, or finite-size error. This page focuses on the computational workflow introduced by the Kohn–Sham construction.

Helpful background. Hartree–Fock computation introduces orbital matrices, and mean-field self-consistency explains nonlinear iteration.

Kohn–Sham equations in a representation

For a spin-unpolarized functional of the density and its gradient (LDA or GGA) in atomic units, solve

$$\left[-\frac12\nabla^2+v_{\rm ext}(\mathbf r) +v_H[n](\mathbf r)+v_{\rm xc}[n](\mathbf r)\right]\phi_i =\varepsilon_i\phi_i,$$

with $n(\mathbf r)=\sum_i f_i|\phi_i(\mathbf r)|^2$ and $\sum_i f_i=N_e$. Each spatial orbital has $0\leq f_i\leq2$. The Hartree potential integrates the Coulomb density, and $v_{\rm xc}=\delta E_{\rm xc}/\delta n$. For a GGA this functional derivative includes spatial derivatives of the density-gradient dependence. Kinetic-energy-density-dependent meta-GGAs and hybrid functionals generally require additional orbital-dependent operators in a generalized Kohn–Sham treatment; they are not covered by the displayed local-potential equation.

In a localized basis, integrate the effective potential to construct a matrix and solve a generalized eigenproblem. In a plane-wave representation, converge both the orbital cutoff and density or potential grid, and match pseudopotentials to the functional and valence treatment. These choices define the represented problem.

The total energy is not the sum of occupied eigenvalues. For the displayed local-potential equation,

$$E=\sum_i f_i\varepsilon_i-E_H[n] +E_{\rm xc}[n]-\int n v_{\rm xc}\,d^3r+E_{\rm nuc}.$$

The subtractions remove contributions counted in the one-particle eigenvalues. Hybrid functionals require a consistent exchange double-counting expression.

A controlled sequence of convergence tests

First converge the density and orbital residuals at a fixed representation. Track electron count, spin polarization, the chosen energy or free energy, and the density change. Then refine numerical quadrature, basis or plane-wave cutoff, reciprocal-space sampling, and cell size independently as applicable. Forces and stress frequently require stricter settings than total energies.

For metals or small gaps, fractional occupations and smearing improve iteration. Record the smearing distribution and width and distinguish an electronic free energy from an extrapolated zero-temperature energy. A smaller smearing width can require denser reciprocal sampling. Mixing parameters influence convergence, but a converged solution should not depend on them when the same stable branch is reached.

For molecular grids, rotate the molecule without changing internal coordinates and check whether energies and forces remain unchanged within tolerance. Excessive orientation dependence indicates insufficient quadrature. For periodic systems, compare symmetry-equivalent reciprocal points and consistent reciprocal weights.

One electron and a diagnostic cancellation

For a normalized one-electron density, the Hartree functional contains an unphysical self-repulsion. The exact exchange-correlation functional cancels it:

$$E_{\rm xc}[n]=-E_H[n]$$

for the exact one-electron ground-state problem, with compatible spin treatment. Common approximate functionals need not satisfy this cancellation. Increasing a grid cannot remove that approximation error. A hydrogen-atom test therefore separates numerical convergence from a physical limitation of the chosen functional.

An independent check computes $N_e=\int n(\mathbf r)d^3r$ from the real-space grid and from orbital occupations or $\operatorname{Tr}(PS)$. Agreement tests density assembly and quadrature but does not certify an accurate energy functional.

Interpreting the result

Kohn–Sham orbitals are auxiliary orbitals that reproduce a density within the construction. Their eigenvalue differences are not automatically neutral excitation energies or quasiparticle gaps. Geometry relaxation adds another optimization layer; forces must include basis-motion terms when the basis depends on nuclear positions. Compare competing spin and structural minima when the physical question allows them.

A useful calculation record states geometry, charge, spin, functional, basis or pseudopotential, grid, reciprocal sampling, occupation scheme, and convergence thresholds. Quantum-chemistry benchmarks explains how to make a numerical comparison that fixes these choices.

Exercises

Why can two fully converged calculations using different exchange-correlation functionals disagree even after their basis and grid limits are established?

Solution

They solve different approximate energy functionals. Numerical convergence removes representation and iteration errors within each model, not the difference between the models and the exact functional. Their spread is useful sensitivity information but is not automatically a calibrated uncertainty interval.

References

• Kohn, W., and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Physical Review 140, A1133–A1138 (1965), doi:10.1103/PhysRev.140.A1133 — auxiliary-orbital construction and energy accounting.
• Psi4 developers, “DFT: Density Functional Theory” documentation, accessed 13 October 2026 — functional classes, quadrature and numerical controls.