Imagining the future of density functional theory

INTRODUCTION TO GROUND-STATE DENSITY FUNCTIONAL THEORY

This section is background for the later, more forward-looking sections. Those who know the fundamentals of DFT can skip to those sections. Those who don’t are reminded that this section is too short and too author-centric to fairly review all the important developments of more than 60 years of this theory.

In the early 1960s, most electronic structure calculations for atoms, molecules, and solids were made in the Hartree approximation (which ignores both exchange and correlation) or the Hartree-Fock approximation (which treats exchange exactly but ignores correlation). In the Hartree approximation, the electron-electron interaction energy is taken to be the positive classical-electrostatic Coulomb energy of the electron charge density interacting with itself, an approximation that ignores the particle nature of the electrons. The exact negative exchange energy arises from two effects: it removes the large self-interaction error of the Hartree approximation, and it accounts for the tendency of any two parallel-spin electrons to avoid one another due to the Pauli exclusion principle. The negative correlation energy further accounts for the tendency of any two electrons to avoid one another due to their mutual Coulomb repulsion. For an N-electron system, both the Hartree and the Hartree-Fock approximations self-consistently compute N orthonormal occupied spin orbitals, and find the electron density as the sum of the absolute squares of these orbitals. Physicists preferred the Hartree method, although it unrealistically underbinds the atoms in molecules and solids. Chemists preferred the Hartree-Fock method, which is better for bonding but still seriously underbinds. The inaccuracy of the Hartree-Fock approximation for atoms and molecules was known to result from the omission of interelectronic Coulomb correlation. There were also some non-rigorous but interesting density functional approximations such as the Thomas-Fermi approximation, which did not bind atoms at all, and Slater’s Xα approximation, which added a scaled local exchange correction to the Hartree approximation.

Modern density functional theory (DFT) starts with a rigorous justification for the use of the electron density n(r) to compute electronic structure. Importantly, the Hamiltonian of a system of $$ N=\int d^{3} r n(\boldsymbol{r}) $$ mutually repelling electrons in the presence of a confining static scalar potential v(r) (typically the Coulomb attraction potential from the nuclei) has universal terms depending only on N and other terms depending as well on v, and the expectation value of this non-universal term is $$ \int d^{3} r n(\boldsymbol{r}) v(\boldsymbol{r}) $$. Pierre Hohenberg and Walter Kohn in 1964 showed^[1] that the ground-state density formally determines the external potential v, and that, for a given v and N, the ground-state energy and electron density are found from the minimum principle

where Q[n] is a universal (independent of v) exact density functional, and the minimum is over all N-electron densities.

Walter Kohn and Lu Sham in 1965 showed^[2] that Q[n] can be split into a sum of a non-interacting kinetic energy that can be constructed exactly from a set of occupied orbitals, an energy $$ \int d^{3} r n(\boldsymbol{r}) v(\boldsymbol{r}) $$ of interaction between the electrons and the external potential, the Hartree electrostatic energy, and a residual exchange-correlation energy E_xc[n] that is relatively small but important for bonding, and is the only term of the energy that must, for practical reasons, be approximated. Effectively, they exactified the Hartree and Hartree-Fock approximations.

Kohn and Sham proved the existence of an exact density functional for the ground-state exchange-correlation energy. The importance of an existence theorem cannot be overstated. Knowing that something important exists, even if it is not practical to evaluate it directly, motivates the persistent search for its exact properties and for better and better approximations to it.

Kohn and Sham also proposed the simple local density approximation

where $$ \varepsilon_{x c}^{\text {unif }}(n) $$ is the exchange-correlation energy per electron of an electron gas of uniform density n. They realized that they could generalize to an external magnetic field that couples to the z component of electron spin by generalizing the scalars n and v to the vectors (n_↑, n_↓) and (v_↑, v_↓).

Walter Kohn once told the author that in 1965 he expected the local density approximation to make only a modest improvement over the Hartree approximation, because real atoms, molecules, and solids are not like the uniform electron gas for which the LSDA is exact. And there was in fact little interest in Kohn-Sham theory until the early 1970s, when numerical tests found that LDA was remarkably accurate for solids and solid surfaces. Shortly after that, Langreth and Perdew^[3], and independently Gunnarsson and Lundqvist^[4], derived the adiabatic connection or fluctuation-dissipation theorem for the exact exchange-correlation (xc) energy of any electron density (non-uniform as well as uniform), and found that LDA and LSDA (the local spin density approximation) are accurate because the xc energy is the electrostatic interaction of each electron with the xc hole in the density that the electron digs around itself, and because the LDA and LSDA xc holes are those of a uniform electron gas, with important mathematical features that are found in all electronic systems.

The work of Langreth and Perdew^[5] also found that the xc hole of the second-order gradient expansion of the energy, valid for densities that vary slowly over space, does not have those features and so this truncated expansion is less accurate than LSDA for real systems. Restoration of those features led to the first generalized gradient approximation^[6], with better-than-LSDA accuracy.

In 1979, Mel Levy^[7] reasoned that Equation (1) could be derived from the wavefunction variational principle by performing the search over normalized antisymmetric N-electron wavefunctions in two steps: first over all wavefunctions yielding a given density, defining a density functional, and then over all densities that belong to wavefunctions. This produces an exact, if computationally intractable, expression for the xc energy functional. Levy and Perdew^[8], among others, derived many mathematical properties of the exact functional that were later used to constrain approximations to the exact functional.

In 1980, Ceperley and Alder^[9] performed a quantum Monte Carlo calculation of $$ \varepsilon_{x c}^{\text {unif }}(n) $$ with high accuracy, which was fitted to analytic expressions^[10,11] including, in most cases, correct high- and low-density limits. After that, nearly all density functionals were essentially exact for uniform electron spin densities. Even in the absence of a magnetic field, spin polarization often arises in the ground state, and using it provides additional information that makes an approximation more accurate.

In part because a single numerical integral is much less expensive than a double one, the most computationally efficient density functional approximations are single integrals over 3D space:

Here τ is the positive orbital kinetic energy density, and the arguments of the function $$ \varepsilon_{x c}^{\text {unif }} $$ are actually resolved into spin-↑ and -↓ components. LSDA uses only the argument n. Generalized gradient approximations (GGAs) add ∇n, and meta-GGAs further add τ. The extra ingredients can be used to satisfy more exact constraints on E_xc[n]. For example, LSDA satisfies 8 exact constraints; the Perdew-Burke-Ernzerhof^[12] (PBE) GGA satisfies 11, and the strongly constrained and appropriately normed (SCAN^[13] and r²SCAN^[14]) meta-GGAs satisfy all 17 exact constraints that a meta-GGA can. More exact constraints make a functional overall more predictive^[15].

Although the self-interaction error is strongly reduced from the Hartree approximation to the Kohn-Sham LSDA, an important residue remains. LSDAs, GGAs, and meta-GGAs are single integrals over three-dimensional space, and so cannot be exact for all one-electron densities, in which the exact xc energy must cancel the Hartree electrostatic energy of the density (a double integral). In 1981, Perdew and Zunger^[16] found a way to make any density functional exact for all one-electron systems (at a price to be discussed in Section PROPER SELF-INTERACTION CORRECTION). This work and observations^[17] about spurious non-integer electron numbers on separated atoms in Slater Xα calculations, led in 1982 to an exact density functional theory for systems of non-integer average electron number^[18].

Because there is a strong and understood error cancellation between x and c, it is wrong to combine full exact exchange with LSDA, GGA, or meta-GGA correlation. But Axel Becke^[19] used the adiabatic connection between the Kohn-Sham non-interacting system and the real interacting system to argue that a smaller fraction of exact exchange should be mixed with the complementary larger fraction of uncorrected exchange in a hybrid functional. Becke’s discussions with the leading theoretical chemist John Pople eventually produced a phase change in which many other chemists (including eventually experimentalists) also became interested. Kohn and Pople shared the 1998 Nobel Prize in Chemistry.

The PBE GGA is in nearly all computer codes and thus highly cited by many physicists, chemists, and materials scientists who use DFT to understand and predict the properties of atoms, molecules, and solids. Since 1996, it has received 242,000 Google Scholar citations, and its annual citations are still growing. r²SCAN is overall more accurate than PBE, and the Materials Project, which computes and compiles the electronic structures of known materials, is switching from PBE to r²SCAN.

From 2020 onward, the method of constraint satisfaction used to approximate the xc energy was also used to approximate the frequency- and wavevector-dependent xc kernel for a uniform electron gas^[20]. The kernel is a key ingredient of the linear response of the density in time-dependent density functional theory. The resulting kernel was used to compute plasmon dispersion and lifetime, and to show that the symmetry-breaking static charge density wave that occurs in a uniform electron gas at low density happens because a plasmon-like density fluctuation drops to zero frequency^[21]. This work provided a numerical confirmation of P. W. Anderson’s time-centered interpretation of symmetry breaking^[22], which is discussed further in Section SYMMETRY BREAKING BY A PROPER DENSITY FUNCTIONAL APPROXIMATION TRANSFORMS STRONG TO NORMAL CORRELATION.

TOWARD A MORE ACCURATE META-GGA

Meta-GGAs employ the exact positive orbital kinetic energy density τ, which reduces to $$ \tau^{W}=\frac{|\nabla n|^{2}}{8 n} $$ in regions of one- or two-electron density and to $$ \tau^{\text {unif }}=\left(\frac{3}{10}\right)\left(3 \pi^{2}\right)^{2 / 3} n^{5 / 3} \\ $$ in regions of spin-unpolarized density varying slowly over space. This quantity is available in any Kohn-Sham calculation. SCAN^[13] and r²SCAN^[14] use the dimensionless ingredient

to interpolate between the one- and two-electron limit (α = 0) and the limit of slow spatial variation (α ≈ 1), and to extrapolate to regions of density-tail overlap ($$ alpha \gg 1 $$). The one-electron limit of these meta-GGAs is not exact, but it is very accurate for compact one-electron densities such as hydrogenic and Gaussian one-electron densities. Thus, these meta-GGAs are not fully self-interaction free. Nevertheless, self-interaction error typically reduces from LSDA to PBE to SCAN. The limit of slowly varying density is described by a gradient expansion. The interpolation is controlled by non-bonded appropriate norms such as rare-gas atoms and the strongly compressed argon dimer.

Although meta-GGAs are computationally semi-local in the sense of Equation (3), as are LSDA and GGAs, meta-GGAs have full nonlocality because the orbitals are fully nonlocal functionals of the electron density. This makes them typically very accurate in atoms, molecules, and insulating solids, where full nonlocality is needed more than it is in metals. The perfect long-range screening of the exchange hole by correlation makes the functional more short-ranged and local in a metal than in a non-metal.

As we climb the ladder of density functional approximations [Figure 1], overall accuracy increases but so does computation time or cost. However, the cost increases from LSDA to GGA to meta-GGA are rather modest compared to the increase from meta-GGA to hybrid or self-interaction correction. One reason for this is that numerical integration can be much more expensive for double integrals than for single integrals. At the same time, meta-GGAs (e.g., SCAN and its smoothed version r²SCAN) are much more accurate than GGAs for molecules and insulating solids, including quantum materials with 3d or 4f electrons. So it is not an idle question to ask how much the accuracy of SCAN and r²SCAN can be further boosted.

Figure 1. Jacob’s ladder of density functional approximations to the exchange-correlation energy of a many-electron system. This energy is approximated as an integral over three-dimensional position space of a function of one or more arguments. The argument added on each rung is listed to the left of that rung: the electron density n as a function of position r, its gradient, the positive orbital kinetic energy density τ, etc. The exact exchange information comes from the one-particle density matrix ρ (r, r') constructed from the occupied Kohn-Sham orbitals or one-electron wavefunctions, which is also the input for proper self-interaction corrections. The proper self-interaction corrections are now added to the hybrid functionals on the fourth rung. The random phase approximation or RPA, on the highest rung, employs also the unoccupied Kohn-Shan orbitals. Each of the most widely-used density-functional approximations sits on one of the rungs.

Here at Tulane, three directions are being explored: (1) Can we improve the description of magnetism in these meta-GGAs, which exaggerate magnetic moments and overstabilize magnetism in metallic Fe, Co, and Ni^[23], while LSDA and PBE do so to a much lesser degree? (2) Can we improve the accuracy of the one- and two-electron limits in these meta-GGAs, without stepping outside the meta-GGA form? (3) The LAK functional^[24] was designed much as r²SCAN was, but with a stronger dependence on the fully nonlocal ingredient τ which leads to larger and more accurate fundamental band gaps in insulators than r²SCAN yields, at the cost of less accurate lattice constants than r²SCAN yields. Is it possible for a single meta-GGA to achieve band gaps as accurate as those of LAK and lattice constants as accurate as those of r²SCAN?

These descriptions are sketchy to protect the effort of the collaborators who are working on these projects: (1) R. Maniar and T. Lebeda. (2) A. Ramasamy, J. Sun, and A. Ruzsinszky. (3) Y. Wang and A. Ruzsinszky.

PROPER SELF-INTERACTION CORRECTION

Rungs 1-3 of the Jacob’s ladder of density functional approximations are LSDA, GGA, and meta-GGA. There is a reason why hybrid functionals (often with empirical fitting) defined the original fourth rung of the ladder of density functional approximations, while self-interaction corrections were excluded: Both are fully nonlocal functionals of the occupied orbitals. But the hybrid functionals improved overall accuracy for atoms, molecules, and insulating solids, while the self-interaction correction (SIC) was more paradoxical. SIC improved some properties, like barrier heights to chemical reactions, while worsening many others, in comparison with the performance of the semi-local functional (LSDA, GGA, or meta-GGA) that SIC was intended to correct. Proper self-interaction correction, as defined below, does not have this problem, and is now added to the fourth rung of the ladder.

The Perdew-Zunger^[16] self-interaction (PZ SIC) subtracts from the uncorrected functional on an orbital-by-orbital basis the self-Hartree and self-XC energies of that orbital. The SIC orbitals must be a unitary transformation of the occupied Kohn-Sham orbitals, and must be localized to make the correction size-extensive. In 2014, the guaranteed-localized FLOSIC^[25] orbitals were found. But they did not solve the paradox of SIC.

In 2019, Santra and Perdew^[26] found that applying PZ SIC to an approximate functional loses that functional’s exactness for all uniform densities, explaining but not solving the paradox of SIC. A first step toward the solution was the locally scaled-down LSIC of Zope et al.^[27], which starts from the energy density for spin σ of the PZ SIC to LSDA, then scales it down locally by a factor of

which equals 1 in one- and two-electron density regions and 0 in uniform-density regions. This is an elegant no-parameter solution that greatly improves LSDA predictions for everything except weak (hydrogen or van der Waals) bonds, which it dramatically worsens. The problem was that LSIC was not able to distinguish between uniform-density regions and regions of density-tail overlap. The natural solution was the LSIC α of Shahi et al.^[28], which scales down the energy density for spin σ by a factor p (a, b, α_σ) of α_σ of Equation (4). Here a ≤ 2 and b > 0 are parameters fitted to a subset of the r²SCAN appropriate norms, and p is a function bounded between 0 and 1 that equals 1 when α_σ is 0 or ∞ (one- or two-electron regions), and equals 0 when α_σ is 1. LSIC α was found to strongly improve upon LSDA for all eight test sets for main-group molecules (atomization energies, barrier heights, ionization potentials, binding energy curves, molecules with strong self-interaction errors, weak bonds, Diels-Alder reaction energies, and halogen/chalcogen-bond energies). The mean absolute error of LSIC α for the weak-bond test set is 2 kcal/mol, a factor of 5 smaller than that of LSIC. Although the mean absolute error of LSDA is also 2 kcal/mol, LSDA consistently overbinds and can only be worsened by a long-range van der Waals correction, while LSIC α is more amenable to such a correction.

A proper self-interaction correction should not lose the correct limits for uniform or slowly varying densities of the approximate functional being corrected, and should “do no harm” to that functional. LSIC α seems to be a proper self-interaction correction to LSDA. The conclusion to be drawn from LSIC α and previous work seems to be that, while all exact constraints are important for the construction of accurate density functionals, the most important constraints are exactness for all uniform densities and for all one-electron densities.

Can we make a proper self-interaction correction to the much more realistic r²SCAN meta-GGA? That will be more difficult than it was for LSDA, for the following reason. The PZ SIC subtracts two terms from the approximate functional for each SIC orbital: the self-Hartree energy, which is positive, and the self-XC energy, which is negative. These two terms, after integration over space, tend to cancel. Energy densities are gauge variant, which means that the same energy can have many different energy densities, one differing from another by a term that integrates over space to zero. The standard LSDA XC energy is in the same gauge as the Hartree energy^[29], so there is a strong local cancellation as well as an integrated cancellation, and as a result both terms of the energy density of the PZ SIC correction can be scaled down by the same scaling function p (a, b, α_σ). In a GGA or meta-GGA, the energy density of the self-xc energy term is in a different gauge from that of the self-Hartree term, which complicates the construction of a proper self-interaction correction. When one is found, the result should be a reliably accurate density functional approximation for all normally correlated systems, as defined in the next section.

SYMMETRY BREAKING BY A PROPER DENSITY FUNCTIONAL APPROXIMATION TRANSFORMS STRONG TO NORMAL CORRELATION

A normally correlated wavefunction can be written as a superposition of a single Slater determinant with heavy weight and many others with small or zero weight. The one-particle density matrix for such a wavefunction has eigenvalues (natural occupation numbers between 0 and 1) that are all close to 1 or close to zero. A normal correlation energy can be zero, as in any one-electron density, but the correlation energy still plays an important role in the bonding of one atom to another. Hartree-Fock theory, which includes full exact exchange but no correlation, produces unrealistically weak bonds. The normally correlated molecules in the AE6 test set have atomization energy errors^[30] of order 6 eV for Hartree-Fock, 3 eV for LSDA, 0.6 eV for the PBE GGA, and 0.1 eV for the r²SCAN meta-GGA. In contrast, strong correlation arises when the wavefunction is a superposition in which two or more degenerate or nearly-degenerate Slater determinants are heavily weighted and strongly mixed by the interelectronic Coulomb interaction. Strong correlation energies are much more negative than normal correlation energies.

The electron gas of uniform density, whose energy is properly described by nearly all density functional approximations, is normally correlated at typical valence- and core-electron densities. While the exact density functional theory describes strong correlation, practical density functional approximations are at best reliable for normal correlation. This appears to be a serious limitation, because many interesting quantum materials, especially those with valence d or f electrons, are strongly correlated. But, at least in many cases, symmetry breaking in a density-functional calculation seems to transform strong correlation into the normal correlation that a practical density functional approximation can describe^[31,32].

The earliest and most dramatic example is found in the stretched H₂ molecule^[4]. As the bond length is stretched, the unoccupied 1s ungerade orbital approaches (but does not reach) exact degeneracy with the occupied 1s gerade orbital. The resulting strong correlation energy of about -6 eV at infinite bond length is reduced to zero correlation when the singlet spin symmetry is broken by localizing the up-spin electron on one nucleus and the down-spin electron on the other, corresponding to two independent H atoms in the limit of infinite stretch. Most density functional approximations do this well.

A subtler situation arises for the singlet ground state of the C₂ molecule at its equilibrium bond length, which is strongly correlated due to an accidental near-degeneracy. Here the extra energy due to strong correlation is only -0.6 eV. Although LSDA, PBE GGA, and SCAN meta-GGA all spontaneously lower the energy by breaking the spin symmetry, only the most reliable functional for normally correlated systems (SCAN) can predict the atomization energy within 0.1 eV by symmetry breaking^[30].

A recent perspective^[32] shows that many transition-metal compound solids that are experimentally insulators display a non-zero fundamental band gap (between unoccupied and occupied one-electron states) only from the combination of SCAN with spontaneous symmetry breaking (electronic or nuclear-positional or both), and not from PBE with or without symmetry breaking.

The Nobel Prize-winning condensed matter theorist P.W. Anderson^[22] proposed a time-centered interpretation of symmetry breaking: Although a stationary-state wavefunction is time-independent, it contains implicit time-dependent fluctuations, and a symmetry breaks when a density or spin-density fluctuation drops to zero or near-zero frequency. This interpretation has been confirmed^[20] recently for the static charge density wave that breaks the translational symmetry of a uniform electron gas at low average density.

Steps have also been taken toward a formal theory^[33] of symmetry breaking in density functional theory. The minimum principle of Equation (1) and Ref.^[7] for the ground state energy implies that the broken-symmetry energy from a proper density functional approximation (defined as one that reliably predicts the energies of normally-correlated states) is bounded from below by the exact energy for a strongly correlated state. The bound is close to the extent that the symmetry-broken density or spin-density expresses the strong correlation in the symmetry-unbroken wavefunction, and also removes the excessive one-electron degeneracies responsible for it.

Superconductivity is an effect that manifests not so much in the one-electron density as in the density of Cooper pairs. Nevertheless, it has been found that superconductivity is often indicated by characteristic patterns^[34] in the subtle breaking of crystal-lattice symmetries.

An important feature of both the Perdew-Zunger self-interaction correction (PZ-SIC) and its locally scaled-down extensions such as LSIC α is^[31] that they are exact, not only for all one-electron densities but for all collections of non-overlapping one-electron densities. This means that they are also exact^[31,35] for all electronic ground states in the broken-symmetry strongly-correlated or classical limit. This limit is one in which the kinetic energy of the electrons tends to zero (because Planck’s constant tends to zero or the electron mass tends to infinity). Then the electrons come to rest at the minima of the total external plus electron-electron potential energy, where their positions are perfectly correlated. Each electron is then in its own highly localized and non-overlapping orbital. A recent calculation^[35] shows that PZ-SIC, applied to any approximate density functional, yields the exact ground-state energy for harmonium (a model system of two real electrons trapped in a spherically symmetric external harmonic oscillator potential) in this limit, and (when applied to LSDA, PBE, or r2SCAN) an accurate ground-state energy for all values of Planck’s constant between 0 and the observed value.

The two most serious challenges to density functional theory, self-interaction error and strong correlation, might be solved together by making a proper self-interaction correction (one that leaves the corrected functional exact for all uniform densities and otherwise does no harm to its accuracy) to a good meta-GGA like r²SCAN or its successor, and then letting all symmetries break.

ARTIFICIAL INTELLIGENCE AND MACHINE LEARNING; EXTENSIONS OF DENSITY FUNCTIONAL THEORY

Advances in artificial intelligence and machine learning are increasingly part of the landscape of our daily lives and of science. They are also increasingly a part of density functional theory. They bring with them the temptation to fit large amounts of data. The author does not see a problem with that per se, but supports the viewpoint that the fitting be constrained by the known mathematical properties of the exact density functional for the exchange-correlation energy, which make a density functional approximation more widely predictive^[15]. In time, artificial intelligence may help us derive more such mathematical properties and impose them on the fitting.

This article is not a balanced review. Only a vastly larger article could review the many important and interesting developments in ground-state density functional theory. There have also been many important extensions of this theory beyond the ground state to the thermal equilibrium state^[36], time-dependent phenomena^[37], excited states^[38,39], and via generalized densities to superconductivity^[40], which will not be reviewed here. The author does want to mention that statistical mechanics shows that symmetry-broken phases are replaced during phase transitions by less symmetry-broken phases and eventually by symmetric equilibrium states as temperature increases. Zentropy^[41] is a broken-symmetry density functional theory for the statistical mechanics of phase transitions, for which artificial intelligence has been helpful^[42,43].

CONCLUSIONS

For the construction of accurate and predictive density functional approximations to the exchange-correlation energy, all exact constraints are important, but the most important constraints seem to be exactness for all uniform densities and all one-electron densities. A proper self-interaction correction (one that is exact in both limits) should describe normal correlation without symmetry breaking, and the strongly-correlated or classical limit through symmetry breaking. Proper self-interaction corrections need to be developed for functionals more advanced than the local spin density approximation, such as r2SCAN, and tested for normal correlation and for the many kinds of realistic strong correlation that never reach the strongly-correlated or classical limit.