Skip to content

Kohn–Sham Density Functional Theory

This chapter begins with the interacting electronic problem at fixed nuclear positions and then summarizes its spin-dependent Kohn–Sham formulation. No spatial discretization or periodicity is assumed. Hartree atomic units are used unless stated otherwise.

Born–Oppenheimer Electronic Problem

Consider nuclei with positions \(\{\boldsymbol R_i\}\), charges \(\{Z_i\}\), and masses \(\{M_i\}\). Within the Born–Oppenheimer approximation, the nuclear positions are fixed parameters of the electronic Hamiltonian,

\[ \hat H_{\mathrm e}(\{\boldsymbol R_i\}) = \hat T_{\mathrm e} +\hat V_{\mathrm{ext}}^{(N_{\mathrm e})} +\hat V_{\mathrm{ee}}. \]

Here \(i\) labels nuclei. For a nonrelativistic all-electron system without applied fields, let \(p\) and \(q\) label electrons and let \(\boldsymbol r_p\) be the position of electron \(p\). The kinetic operator and the local electron–nucleus potential are

\[ \hat T_{\mathrm e} = -\frac{1}{2}\sum_p \nabla_p^2, \qquad v_{\mathrm{ext}}(\boldsymbol r) = -\sum_i \frac{Z_i} {\lvert\boldsymbol r-\boldsymbol R_i\rvert}, \]

with the corresponding many-electron external operator

\[ \hat V_{\mathrm{ext}}^{(N_{\mathrm e})} = \sum_p v_{\mathrm{ext}}(\boldsymbol r_p). \]

The electron–electron interaction is

\[ \hat V_{\mathrm{ee}} = \frac{1}{2} \sum_{p\ne q} \frac{1} {\lvert\boldsymbol r_p-\boldsymbol r_q\rvert}. \]

The clamped-nuclei electronic ground-state energy \(\mathcal E_0\) defines the Born–Oppenheimer energy

\[ E_{\mathrm{BO}}(\{\boldsymbol R_i\}) = \mathcal E_0(\{\boldsymbol R_i\}) +E_{\mathrm{NN}}(\{\boldsymbol R_i\}), \]

where \(E_{\mathrm{NN}}\) is the nucleus–nucleus interaction. For a finite set of point nuclei,

\[ E_{\mathrm{NN}} = \frac{1}{2} \sum_{i\ne j} \frac{ Z_i Z_{j} }{ \lvert \boldsymbol R_i-\boldsymbol R_{j} \rvert }. \]

Relativistic corrections and effective-core treatments modify the one-electron part of the Hamiltonian and may introduce nonlocal and spin-dependent operators. They do not change the separation between the fixed nuclear geometry and the electronic problem.

Kohn–Sham Spinor Space

The Kohn–Sham construction replaces the interacting ground-state problem by an auxiliary noninteracting spinor system constrained to reproduce the same local spin-density matrix. The Kohn–Sham spinors construct this density and its one-particle density operator; they are not additional basic density variables. The one-electron Hilbert space is

\[ \mathcal H_1 = L^2(\mathcal V)\otimes\mathbb C^2, \]

where \(\mathcal V\) is the spatial domain and \(\mathbb C^2\) is spin space. A Kohn–Sham state labeled by \(a\) is represented by a two-component spinor,

\[ \psi_a(\boldsymbol r) = \begin{pmatrix} \psi_{a\uparrow}(\boldsymbol r)\\ \psi_{a\downarrow}(\boldsymbol r) \end{pmatrix}. \]

The spatial domain and its boundary conditions are part of the physical problem. Isolated systems use boundary conditions appropriate to finite systems. Translation symmetry and Bloch boundary conditions are introduced separately in Periodic Systems.

A spin-matrix quantity \(\mathbf X\) admits the Pauli decomposition

\[ \mathbf X = X^0\sigma_0 +\boldsymbol X\cdot\boldsymbol\sigma, \]

where \(\sigma_0\) is the spin identity and \(\boldsymbol\sigma=(\sigma_x,\sigma_y,\sigma_z)\) collects the Pauli matrices. The components are

\[ X^0 = \frac{1}{2}\operatorname{tr}_s\mathbf X, \qquad X^\alpha = \frac{1}{2} \operatorname{tr}_s(\sigma_\alpha\mathbf X), \qquad \alpha=x,y,z. \]

Here \(\operatorname{tr}_s\) denotes the trace over spin indices. This decomposition is algebraic; the scalar and vector parts acquire their physical meaning from the quantity being decomposed.

Density Operator and Local Density

Let \(\{\lvert\psi_a\rangle\}\) be orthonormal Kohn–Sham spinors and let \(f_a\) be their occupations. In the spinor formulation, \(0\le f_a\le1\). The one-particle density operator is

\[ \hat\rho = \sum_a f_a \lvert\psi_a\rangle \langle\psi_a\rvert. \]

For an infinite periodic system, the Kohn–Sham states also carry a \(\boldsymbol k\) label, and the density is obtained by a normalized Brillouin-zone average over the occupied Bloch states. Electron numbers and other extensive quantities are reported per cell because their totals for the infinite crystal are not finite. The detailed convention is defined in Periodic Systems.

For a finite system, the real-space kernel has spin components

\[ \rho_{ss'}(\boldsymbol r,\boldsymbol r') = \langle\boldsymbol r,s \lvert\hat\rho\rvert \boldsymbol r',s'\rangle = \sum_a f_a \psi_{as}(\boldsymbol r) \psi_{as'}^*(\boldsymbol r'). \]

Hermiticity of \(\hat\rho\) implies

\[ \rho_{ss'}(\boldsymbol r,\boldsymbol r') = \rho_{s's}^*(\boldsymbol r',\boldsymbol r). \]

The local spin-density matrix \(\mathbf n(\boldsymbol r)\) is formed from the diagonal kernel:

\[ n_{ss'}(\boldsymbol r) \equiv \rho_{ss'}(\boldsymbol r,\boldsymbol r), \qquad \mathbf n(\boldsymbol r) = \left( n_{ss'}(\boldsymbol r) \right)_{ss'}. \]

Its Pauli decomposition is

\[ \mathbf n(\boldsymbol r) = \frac{1}{2} \left[ n(\boldsymbol r)\sigma_0 +\boldsymbol m(\boldsymbol r) \cdot\boldsymbol\sigma \right], \]

where the scalar electron density and Pauli spin-polarization density are

\[ n(\boldsymbol r) = \operatorname{tr}_s\mathbf n(\boldsymbol r), \qquad m^\alpha(\boldsymbol r) = \operatorname{tr}_s \left[ \sigma_\alpha\mathbf n(\boldsymbol r) \right]. \]

Conversion of \(\boldsymbol m\) to a magnetic-moment density introduces the corresponding magnetic-moment factor and sign convention.

For a finite system, the electron number is

\[ N_{\mathrm e} = \operatorname{tr}\hat\rho = \int_{\mathcal V} \mathrm d\boldsymbol r\, \operatorname{tr}_s \mathbf n(\boldsymbol r). \]

Here \(\operatorname{tr}\) without a subscript denotes the trace over the full one-electron Hilbert space. In an infinite periodic problem, the corresponding quantity per cell is a normalized Brillouin-zone average. For a one-body operator \(\hat X\),

\[ \langle\hat X\rangle = \operatorname{tr}(\hat\rho\hat X). \]

The notation distinguishes three forms of the same one-particle object: \(\hat\rho\) is the density operator, \(\rho_{ss'}(\boldsymbol r,\boldsymbol r')\) are the components of its real-space kernel, and \(\mathbf D\) is its density matrix after a basis has been chosen. In a nonorthogonal basis, \(\mathbf D\) must not be identified with the covariant matrix elements of \(\hat\rho\).

Energy Functional

The Kohn–Sham energy functional is written as

\[ E_{\mathrm{KS}}[\hat\rho] = T_{\mathrm s}[\hat\rho] +E_{\mathrm{ext}}[\hat\rho] +E_{\mathrm H}[n] +E_{\mathrm{xc}}[n,\boldsymbol m] +E_{\mathrm{NN}}. \]

Here \(T_{\mathrm s}\) is the kinetic energy of the noninteracting Kohn–Sham system, \(E_{\mathrm{ext}}\) is its interaction with the external one-electron operator, \(E_{\mathrm H}\) is the classical Hartree energy, and \(E_{\mathrm{xc}}\) contains the nonclassical electron–electron contribution and the difference between the interacting and noninteracting kinetic energies. \(E_{\mathrm{NN}}\) is the nuclear interaction at the fixed nuclear geometry.

The noninteracting kinetic energy and external-potential contribution are

\[ T_{\mathrm s}[\hat\rho] = \operatorname{tr} \left( \hat\rho \left[-\frac{1}{2}\nabla^2\sigma_0\right] \right), \qquad E_{\mathrm{ext}}[\hat\rho] = \operatorname{tr} \left( \hat\rho\hat V_{\mathrm{ext}} \right). \]

For the local spin-independent external potential introduced above,

\[ E_{\mathrm{ext}} = \int \mathrm d\boldsymbol r\, v_{\mathrm{ext}}(\boldsymbol r) n(\boldsymbol r). \]

For a finite system with isolated Coulomb boundary conditions, the Hartree energy is the classical Coulomb energy of the electron density,

\[ E_{\mathrm H}[n] = \frac{1}{2} \iint \mathrm d\boldsymbol r\, \mathrm d\boldsymbol r'\, \frac{ n(\boldsymbol r)n(\boldsymbol r') }{ \lvert\boldsymbol r-\boldsymbol r'\rvert }. \]

Its functional derivative with respect to \(n\) defines the Hartree potential,

\[ v_{\mathrm H}(\boldsymbol r) = \frac{\delta E_{\mathrm H}} {\delta n(\boldsymbol r)} = \int \mathrm d\boldsymbol r'\, \frac{ n(\boldsymbol r') }{ \lvert\boldsymbol r-\boldsymbol r'\rvert }. \]

For a periodic Coulomb system, the electron–electron, electron–ion, and ion–ion terms must instead be defined with one consistent periodic electrostatic convention, including a common treatment of the \(\boldsymbol G=0\) component. The bare contributions need not be finite separately even when the unit cell is neutral.

The scalar and spin-vector components of the exchange–correlation potential are likewise defined as functional derivatives,

\[ v_{\mathrm{xc}}^0(\boldsymbol r) \equiv \frac{\delta E_{\mathrm{xc}}} {\delta n(\boldsymbol r)}, \qquad v_{\mathrm{xc}}^\alpha(\boldsymbol r) \equiv \frac{\delta E_{\mathrm{xc}}} {\delta m^\alpha(\boldsymbol r)}, \qquad \alpha=x,y,z. \]

Equivalently, for arbitrary infinitesimal variations of \(n\) and \(\boldsymbol m\), the first variation of the functional is

\[ \delta E_{\mathrm{xc}} = \int \mathrm d\boldsymbol r\, \left[ v_{\mathrm{xc}}^0(\boldsymbol r) \delta n(\boldsymbol r) + \boldsymbol v_{\mathrm{xc}}(\boldsymbol r) \cdot \delta\boldsymbol m(\boldsymbol r) \right], \]

where \(\delta\) denotes a functional variation, not the difference between two SCF iterations. The \(2\times2\) multiplicative potential acting on spinors is therefore

\[ \mathbf v_{\mathrm{xc}}(\boldsymbol r) = v_{\mathrm{xc}}^0(\boldsymbol r)\sigma_0 +\boldsymbol v_{\mathrm{xc}}(\boldsymbol r) \cdot\boldsymbol\sigma. \]

With the Pauli convention above, the same definition can be written without components as

\[ \delta E_{\mathrm{xc}} = \int \mathrm d\boldsymbol r\, \operatorname{tr}_s \left[ \mathbf v_{\mathrm{xc}}(\boldsymbol r) \delta\mathbf n(\boldsymbol r) \right]. \]

This notation covers spin-unpolarized, collinear, and noncollinear spin-density functionals. Orbital-dependent functionals may lead instead to a generalized Kohn–Sham Hamiltonian with additional nonlocal terms.

Kohn–Sham Hamiltonian

Variation of the energy with respect to the Kohn–Sham spinors, subject to their orthonormality, gives the effective one-electron Hamiltonian

\[ \hat H_{\mathrm{KS}}[\mathbf n] = -\frac{1}{2}\nabla^2\sigma_0 +\hat V_{\mathrm{ext}} +v_{\mathrm H}[n](\boldsymbol r)\sigma_0 +\mathbf v_{\mathrm{xc}}[\mathbf n](\boldsymbol r). \]

The external operator may include nonlocal and spin-dependent terms, including effective core operators and spin–orbit coupling. The Kohn–Sham equations are

\[ \hat H_{\mathrm{KS}}[\mathbf n] \lvert\psi_a\rangle = \varepsilon_a \lvert\psi_a\rangle, \qquad \langle\psi_a\vert\psi_b\rangle = \delta_{ab}. \]

The solutions and their occupations reconstruct \(\hat\rho\) as defined above. Its diagonal real-space kernel gives the updated \(\mathbf n(\boldsymbol r)\) and closes the self-consistency condition

\[ \mathbf n \longmapsto \hat H_{\mathrm{KS}}[\mathbf n] \longmapsto \hat\rho \longmapsto \mathbf n. \]

At self-consistency, \(\hat\rho\) commutes with \(\hat H_{\mathrm{KS}}\), so they may be diagonalized simultaneously. If occupations are assigned by a scalar energy-dependent rule, with a common occupation throughout each degenerate eigenspace, then

\[ \hat\rho = f(\hat H_{\mathrm{KS}}-\mu), \]

for an occupation function \(f\) and chemical potential \(\mu\).

Occupations

For a finite system in the spinor formulation,

\[ \sum_a f_a=N_{\mathrm e}. \]

With cell-normalized Bloch states, the periodic electron count per cell is

\[ N_{\mathrm e} = \frac{1}{\Omega_{\mathrm{BZ}}} \int_{\mathrm{BZ}} \mathrm d\boldsymbol k\, \sum_a f_{a\boldsymbol k}. \]

The normalization and Brillouin-zone convention are detailed in Periodic Systems.

At zero temperature, a ground state with integer occupations has \(\hat\rho^2=\hat\rho\). Degeneracies at the Fermi level, ensemble treatments, finite electronic temperature, and numerical smearing may produce fractional occupations. For a physical electronic temperature,

\[ f_a = \frac{1}{ 1+\exp[(\varepsilon_a-\mu)/(k_{\mathrm B}T)] }. \]

The stationary thermodynamic potential at finite temperature includes the corresponding electronic entropy. Numerical smearing schemes need not represent a physical temperature and define their own energy corrections.

A spinless nonmagnetic formulation may combine two degenerate spin states into one orbital with occupation up to two. Such a factor is a reduced formulation convention; it is not present in the spinor formulation above.

Spin Structure

The general object is the full \(2\times2\) spin-density matrix.

For a nonmagnetic density, \(\boldsymbol m(\boldsymbol r)=0\). This does not by itself exclude spin–orbit coupling: a time-reversal-symmetric spinful Hamiltonian may have a nontrivial spin structure while its ground-state magnetization vanishes.

In a collinear calculation, a global spin axis can be chosen such that \(\mathbf n\) and \(\hat H_{\mathrm{KS}}\) are block diagonal in spin. The two spin channels may then be represented separately. In the noncollinear case, the off-diagonal spin components are retained and the Kohn–Sham equations do not separate into independent spin channels.

Spin–orbit coupling is a term in the Hamiltonian rather than a magnetic state. It couples spin space to the spatial degrees of freedom and generally removes spin as a conserved quantum number.

Total Energy and Forces

At zero electronic temperature, for a multiplicative Kohn–Sham exchange–correlation potential, define the occupation-weighted eigenvalue sum by

\[ E_{\mathrm{band}} = \begin{cases} \displaystyle \sum_a f_a\varepsilon_a, & \text{finite system},\\ \displaystyle \frac{1}{\Omega_{\mathrm{BZ}}} \int_{\mathrm{BZ}} \mathrm d\boldsymbol k\, \sum_a f_{a\boldsymbol k}\varepsilon_{a\boldsymbol k}, & \text{periodic system, per cell}. \end{cases}. \]

The self-consistent total energy is then

\[ E_{\mathrm{KS}} = E_{\mathrm{band}} -E_{\mathrm H}[n] +E_{\mathrm{xc}}[n,\boldsymbol m] -\int \mathrm d\boldsymbol r\, \operatorname{tr}_s \left[ \mathbf n(\boldsymbol r) \mathbf v_{\mathrm{xc}}(\boldsymbol r) \right] +E_{\mathrm{NN}}. \]

For a periodic system, the resulting total energy is per unit cell. The Coulomb terms in this decomposition use the common periodic electrostatic convention stated above and must not be interpreted as separately convergent bare infinite-crystal sums.

The eigenvalue sum is therefore not the total energy. Generalized Kohn–Sham functionals require their formulation-specific double-counting corrections. A norm-conserving effective-core calculation also uses dataset-consistent ionic and reference terms, while PAW uses its complete smooth-plus-one-center energy.

Forces are derivatives of the stationary Born–Oppenheimer energy,

\[ \boldsymbol F_i = -\frac{\mathrm d E_{\mathrm{BO}}} {\mathrm d\boldsymbol R_i}. \]

References