Harris functional

Kohn–Sham equations are the one-electron equations that must be solved in a self-consistent fashion in order to find the ground state density of a system of interacting electrons:

${\displaystyle \left({\frac {-\hbar ^{2}}{2m}}\nabla ^{2}+v_{\rm {H}}[n]+v_{\rm {xc}}[n]+v_{\rm {ext}}(r)\right)\phi _{j}(r)=\epsilon _{j}\phi _{j}(r).}$

The density, ${\displaystyle n,}$ is given by that of the Slater determinant formed by the spin-orbitals of the occupied states:

${\displaystyle n(r)=\sum _{j}f_{j}\vert \phi _{j}(r)\vert ^{2},}$

where the coefficients ${\displaystyle f_{j}}$ are the occupation numbers given by the Fermi–Dirac distribution at the temperature of the system with the restriction ${\textstyle \sum _{j}f_{j}=N}$, where ${\displaystyle N}$ is the total number of electrons. In the equation above, ${\displaystyle v_{\rm {H}}[n]}$ is the Hartree potential and ${\displaystyle v_{\rm {xc}}[n]}$ is the exchange–correlation potential, which are expressed in terms of the electronic density. Formally, one must solve these equations self-consistently, for which the usual strategy is to pick an initial guess for the density, ${\displaystyle n_{0}(r)}$, substitute in the Kohn–Sham equation, extract a new density ${\displaystyle n_{1}(r)}$ and iterate the process until convergence is obtained. When the final self-consistent density ${\displaystyle n(r)}$ is reached, the energy of the system is expressed as:

${\displaystyle E[n]=\sum _{j\in {\text{occupied}}}\epsilon _{j}-{\tfrac {1}{2}}\int v_{\rm {H}}[n]n(r)\,\mathrm {d} r-\int v_{\rm {xc}}[n]n(r)\,\mathrm {d} r+E_{\rm {xc}}[n]}$.

Assume that we have an approximate electron density ${\displaystyle n_{0}(r)}$, which is different from the exact electron density ${\displaystyle n(r)}$. We construct exchange-correlation potential ${\displaystyle v_{\rm {xc}}(r)}$ and the Hartree potential ${\displaystyle v_{\rm {H}}(r)}$ based on the approximate electron density ${\displaystyle n_{0}(r)}$. Kohn–Sham equations are then solved with the XC and Hartree potentials and eigenvalues are then obtained; that is, we perform one single iteration of the self-consistency calculation. The sum of eigenvalues is often called the band structure energy:

${\displaystyle E_{\rm {band}}=\sum _{i}\epsilon _{i},}$

where ${\displaystyle i}$ loops over all occupied Kohn–Sham orbitals. The Harris energy functional is defined as

${\displaystyle E_{\rm {Harris}}[n_{0}]=\sum _{i}\epsilon _{i}-\int \mathrm {d} r^{3}v_{\rm {xc}}[n_{0}](r)n_{0}(r)-{\tfrac {1}{2}}\int \mathrm {d} r^{3}v_{\rm {H}}[n_{0}](r)n_{0}(r)+E_{\rm {xc}}[n_{0}]}$

It was discovered by Harris that the difference between the Harris energy ${\displaystyle E_{\rm {Harris}}}$ and the exact total energy is to the second order of the error of the approximate electron density, i.e., ${\displaystyle O((\rho -\rho _{0})^{2})}$. Therefore, for many systems the accuracy of Harris energy functional may be sufficient. The Harris functional was originally developed for such calculations rather than self-consistent convergence, although it can be applied in a self-consistent manner in which the density is changed. Many density-functional tight-binding methods, such as CP2K, DFTB+, Fireball,[2] and Hotbit, are built based on the Harris energy functional. In these methods, one often does not perform self-consistent Kohn–Sham DFT calculations and the total energy is estimated using the Harris energy functional, although a version of the Harris functional where one does perform self-consistency calculations has been used.[3] These codes are often much faster than conventional Kohn–Sham DFT codes that solve Kohn–Sham DFT in a self-consistent manner.

While the Kohn–Sham DFT energy is a variational functional (never lower than the ground state energy), the Harris DFT energy was originally believed to be anti-variational (never higher than the ground state energy).[4] This was, however, conclusively demonstrated to be incorrect.[5][6]