Stretched or noded orbital densities and self-interaction correction in density functional theory

Kohn-Sham density functional theory^1,2 simplifies the many-electron ground-state problem of condensed-matter physics or quantum chemistry to a self-consistent one-electron form, in a way that is formally exact for the total energy and electron density. To make this theory tractable and widely useful, an approximation (typically a computationally efficient semilocal one) must be made to the exact density functional for the exchange-correlation energy, the many-body “glue” that binds atoms to form molecules and solids.

A semilocal approximation expresses this energy as a single integral over three-dimensional space of a function of various ingredients that are available from the solution of the self-consistent one-electron Schrödinger equation. The local spin density approximation (LSDA)^1,3 uses only the local electron spin densities as ingredients, the generalized gradient approximation (GGA) adds their gradients, and the meta-GGA further adds the positive spin-resolved orbital kinetic energy densities

Square brackets in Eq. (1) indicate a functional, and round brackets indicate a function. The total electron density $nr=n↑r+n↓(r)$ is the sum of up- and down-spin contributions.

The integrand of Eq. (1) can be constructed empirically by fitting to molecular data or in a more first-principles way by satisfying known exact constraints and appropriate norms. Within a semilocal form, many but not all exact constraints or mathematical properties of the exact functional can be satisfied for all possible electron densities by the careful design of the integrand. The uniform electron gas³ is an appropriate norm for the LSDA, for the Perdew-Burke-Ernzerhof (PBE) GGA,⁴ and for the strongly constrained and appropriately normed (SCAN) meta-GGA⁵ that will be considered here. Free atoms (but not molecules) provide other appropriate norms for SCAN, which further satisfies all 17 known exact constraints that a semilocal form can satisfy. All three ingredients n_σ, |∇n_σ|, τ_σ in Eq. (1) are needed to recognize one- and two-electron regions for which there are special constraints (e.g., the correlation energy vanishes for any one-electron density). However, once a meta-GGA has been constructed, it is possible if desired to approximate τ_σ in terms of n_σ, |∇n_σ|, ∇²n_σ.⁶ 

Recent successes, e.g., Refs. 7–12, suggest that SCAN may be approaching the limiting accuracy of semilocal functionals. It may be time to satisfy more exact constraints via fully nonlocal functionals. Semilocal spin-density functionals cannot be exact for all one-electron (fully spin-polarized) densities n, for which the exact spin-density functional for the exchange-correlation energy must exactly cancel the fully nonlocal Hartree electrostatic energy

For a one-electron density, all of the U[n] term in the total energy is a spurious self-interaction. The residual self-interaction error of a semilocal functional for a one-electron density, due to an imperfect cancellation of U[n] by the semilocal exchange-correlation energy, manifests in serious errors for stretched bonds that can arise in binding energy curves (e.g., for H₂⁺,^13,14 as shown in Fig. 1), charge transfers,¹⁵ energies of transition states for chemical reactions,¹⁶ and in other ways.

In 1981, Perdew and Zunger¹⁷ proposed a self-interaction correction (SIC) to any spin-density functional approximation

Here, n_iσ(r) is the ith occupied orbital density of spin σ and, thus, a one-electron density. The resulting fully nonlocal functional is exact for any collection of separated one-electron densities, and the Perdew-Zunger (PZ) correction to the exact functional vanishes.

PZ SIC was an early attempt to improve approximate functionals via constraint satisfaction and was proposed without much concern that satisfying an additional constraint or norm (exactness for all one-electron densities) might violate pre-existing constraints or norms (such as exactness for all uniform electron gases and 0% error for the exchange-correlation energies of neutral atoms in the limit of large atomic number¹⁸). A recent study¹⁹ suggests that the SIC errors for these norms are not negligible.

When PZ SIC was proposed, LSDA was almost the only available density functional approximation for exchange and correlation. LSDA is not very accurate for the compact one- and two-electron hydrogen and helium atoms, and those errors are corrected by PZ SIC. The later GGAs and meta-GGAs also largely correct those errors, but they continue to make self-interaction errors for electron transfer and for the lobed one-electron densities that will be discussed in Sec. II.

The occupied orbital densities n_iσ in Eq. (3) are not the squares of the occupied Kohn-Sham orbitals, which can be delocalized. To make SIC size-extensive, they should be the squares of unitarily equivalent occupied localized orbitals. This was discussed in Ref. 17, but implemented to minimize the SIC energy by Pederson, Heaton, and Lin,²⁰ enabling the first SIC calculations for molecules. The recently proposed Fermi-Löwdin orbitals,²¹ which are guaranteed to be localized, correspond to a particular set of real unitary transformations. Without such a restriction on the unitary transformation, delocalized orbitals that violate size-extensivity can be obtained when the underlying functional has a positive self-interaction correction from any localized orbital. This is unlikely for atoms and molecules in LSDA-SIC, but it can happen for PBE-SIC or SCAN-SIC, where the underlying functional is more accurate for compact one-electron densities.

In the absence of an external magnetic field or a spin-orbit interaction, the Kohn-Sham orbitals are usually chosen to be real (although complex Kohn-Sham orbitals may bring some benefits for orbital functionals²²). Thus, it was long assumed that the localized orbitals should also be real. Real orbitals that overlap one another must have nodes to ensure mutual orthogonality. Vydrov and Scuseria²³ showed that PZ SIC with real localized orbitals, applied to PBE, worsens the description of equilibrium (unstretched) bonds and that conclusion is supported by later studies.^24–27 It has been suspected that noded localized orbitals limit the accuracy of PBE-SIC.²⁸ In 2011, Klüpfel, Klüpfel, and Jónsson²⁹ found that complex localized orbitals with nodeless orbital densities improve the accuracy of PBE-SIC over real ones. (While the real and imaginary parts of an orbital may have nodal surfaces, these surfaces can be different, leading to nodeless orbital densities.^29,30) Complex localized orbitals have also been found²⁵ to reduce but not eliminate PZ SIC symmetry-breaking in molecules. One argument²⁹ for their better performance is that introducing complex localized orbitals increases variational freedom and produces lower and more realistic PBE-SIC total energies. References 25 and 29 found that complex orbitals yield lower energies than real ones for most systems. The problems that noded orbitals pose for semilocal functionals were also discussed in Refs. 28 and 30. Hofmann, Klüpfel, Klüpfel, and Kümmel³⁰ say that “The smoother orbital densities without nodal planes are much closer to the realm where (semi) local functionals are considered to be appropriate.”

Here, we will argue that semilocal approximations must fail for noded orbital densities for the same fundamental reason and in the same way that they must fail for stretched orbital densities. Perhaps the best-known failures of semilocal functionals occur for stretched orbital densities. In applications to ground states, the stretched orbitals are problematic for semilocal functionals, while the noded orbitals are problematic only for SIC functionals.

PZ SIC with real localized orbitals does indeed correct errors for stretched bonds, as shown in Fig. 1, but when applied to a semilocal functional, it introduces other errors arising from the nodes of the orbital densities. The density-nodality error of PZ SIC is absent from uncorrected semilocal functionals applied to ground-state total densities, which are always nodeless for real systems. The most problematic ground states for the semilocal functionals are those of stretched radicals, such as H₂⁺, He₂⁺, and transition states of chemical reactions. Large errors in the total energies from semilocal approximations would also occur for strongly stretched H₂ and N₂ if spin symmetry were not allowed to break.^14,31

It is well-known that the exact Kohn-Sham theory is exact for ground-state energies and densities but not for most excited states of many-electron systems. The argument at the end of Sec. II shows that it is also exact for the excited states of one-electron systems.

A formally exact expression for the exchange-correlation energy is^31,32

i.e., a Coulomb interaction between the electron density at a position r and the density at r′ of the exchange-correlation hole surrounding an electron at r. Around a given electron, one electron is missing from the rest of the density, leading to the sum rule^31,32

The hole density is thus typically negative, and the more negative it is close to its electron the more negative the exchange-correlation energy will be. Although we do not usually know the hole density exactly, its known exact properties can be used to understand the power and limitations of approximate density functionals^31,32 and to provide an alternative and earlier construction³³ of the PBE GGA.⁴ 

Semilocal approximations either explicitly³³ or implicitly³⁴ model the hole density as a function of r′ that is localized (on a scale set by the local spatial variation of the electron density) around the electron position r and that satisfies the sum rule of Eq. (5). This model is reasonably correct in many cases. A functional like PBE or SCAN can be accurate both for the extended densities in solids and for the compact densities of atoms or of molecules at equilibrium geometries, while LSDA is less accurate, especially for the compact densities. However, for any one-electron density n, the exact hole density is −n(r′), independent of r, making Eq. (4) reduce to

or minus the Hartree electrostatic energy of the one-electron density. Unless the one-electron density is reasonably compact, a semilocal approximation cannot mimic the full nonlocality of Eq. (6) (Fig. 1). While LSDA cannot mimic Eq. (6), PBE and SCAN can do so for compact orbital densities (as in the small-R part of Fig. 1). However, if a one-electron density is separated into two equal lobes (as in the large-R part of Fig. 1 and in Fig. 2), then the exact hole density is shared equally between the lobes, while the semilocal hole density is concentrated on the lobe nearer to the electron, making the PBE or SCAN exchange-correlation energy too negative. In the limit of large separation between the lobes, the approximate hole is entirely on the lobe nearer the electron.

To illustrate this problem, Fig. 2 shows the electron density of the stretched one-electron ion H₂⁺. The exact exchange-correlation hole density for any electron position r is −n(r′), independent of r, or Fig. 2 turned upside down. The LSDA hole density at r′ for electron position r is a function of |r′ − r| that minimizes at the value −n(r) when r′ = r and then rises to zero over a length scale roughly proportional to n(r)^−1/3. This qualitative failure of LSDA is shared by GGA and meta-GGA for stretched H₂⁺ and other stretched radicals.

From the very beginning of density-functional theory, it was known that the exact density functional for the exchange-correlation energy is nonlocal. For example, the exact exchange-correlation potential around an atom decays as −1/r as r → ∞, like a Coulomb potential. A more radical nonlocality, which does not decay to zero with increasing separation, was identified by Perdew, Parr, Levy, and Balduz,³⁵ who were motivated in part by PZ SIC. This radical nonlocality is associated with a noninteger average electron number on a separated subsystem. Stretched H₂⁺, with half an electron (and half an exact exchange-correlation hole) around each separated proton, was identified as an example by Zhang and Yang,¹³ who found wrong dissociation limits with the BLYP GGA and B3LYP hybrid functionals and made a quantitative model to explain this. It is clear from Ref. 13 that only a fully nonlocal functional such as Hartree-Fock or PZ SIC can account for the energies of such systems. No semilocal approximation can be correct for such systems, where the error of a semilocal approximation grows gradually with increasing separation as in Fig. 3. The errors that semilocal approximations make for noded orbital densities reflect the incipient division of one electron into fragments of noninteger electron number.

GGAs and meta-GGAs can accurately describe compact one-electron densities such as the 1s density in H or He⁺. For example, the total energy of the hydrogen atom in hartree is −0.479 (LSDA), −0.500 (PBE or SCAN), or −0.500 (exact). It is well known that stretched orbitals, as in the stretched one-electron molecular ion H₂⁺, are poorly described even by GGAs and meta-GGAs.¹³ It is less well known (but see Ref. 36) that noded orbitals, which are also lobed, encounter a similar self-interaction error. Figure 2 compares the densities of two different systems built up from one electron and two protons. The first is H₂⁺ in its 1σ even-parity ground state at the stretched bond length R = 4 bohrs, with its density plotted along the bond axis. The second is the He⁺ atomic ion in its 2p_z excited state, with its density plotted along the z axis in Fig. 2. These systems are similarly lobed, and for both of them, PBE or SCAN makes large and remarkably similar relative errors (Table I). (The 2p_z state of He⁺ is introduced because its density is strongly noded, in much the same way that the H₂⁺ density in Fig. 2 is strongly stretched. A noded density is the limit of a sequence of un-noded ground-state densities, as discussed at the end of this section.)

While a compact orbital density is clearly a single system with an integer (1) electron number, stretching or noding can be regarded as early steps in the formation of separate systems with a noninteger average electron number for which semilocal approximations always fail dramatically.

Figure 3 shows that the PBE and SCAN errors almost go to zero as the bond length R of H₂⁺ is reduced to zero. Similarly, the large nodality errors for the 2p_z orbital density are greatly reduced but not eliminated for typical real SIC localized orbitals of atoms and molecules (Table II). The 2p_z orbital density is strongly noded, with a nodal plane through its center, but the real Fermi-Löwdin orbitals considered in Table II are more weakly noded, with nodal surfaces nearer their edges.

Since ground-state densities are typically nodeless, can we even define the exact ground-state exchange-correlation energy of a noded density? We can do so by regarding a noded one-electron density as the limit of a sequence of nodeless one-electron ground-state densities. For example, the 2p_z orbital density in Fig. 2(b) varies like z² near the nodal plane z = 0 and can be obtained as the a² → 0 limit of a sequence of one-electron densities that behave like z² + a² near z = 0. A similar limit can be defined for any noded one-electron density, and the sequence of un-noded densities will all be v-representable ground-state densities by inversion of the one-electron Schrödinger equation for the positive square root of the density to find the potential v(r). The integral of Eq. (2) then tells us that, for a noded and for a nodeless one-electron density n, Eq. (6) is valid. In other words, the exact ground-state exchange-correlation functional for spin-polarized one-electron ground-state densities, Eq. (6), applied to a sequence of nodeless spin-polarized one-electron ground-state densities that lie within its domain of definition implies Eq. (6) even for noded spin-polarized one-electron densities. The exact exchange-correlation energy of a spin-polarized one-electron density is Eq. (6), which depends only on the amplitude and not on the phase of the one-electron wavefunction. Thus, as claimed in Ref. 17, the PZ correction on the extreme right of Eq. (3) does vanish when applied to the exact functional. Unfortunately, this does not imply that the correction will be acceptably small when applied to a semilocal functional in a situation for which that functional is accurate, as discussed in Sec. III.

From the discussion of the previous paragraph, the exact Kohn-Sham theory is exact for the ground and excited states of any one-electron density although, of course, not for all excited states of many-electron systems, and PZ SIC is exact for all ground and excited states of one-electron systems.

GGA and meta-GGA functionals of the total density often work well, but not in regions where there are stretched orbitals. Our previous example was stretched H₂⁺,¹³ but the same effect occurs for the stretched bonds between different open-shell atoms¹⁵ and in the transition states¹⁶ of chemical reactions, where again the semilocal exchange-correlation energy is too negative, making the energy barriers too low.

A standard cure is the Perdew-Zunger (PZ) self-interaction correction (SIC) of Eq. (3). To work well, the PZ correction from each occupied localized orbital should be small (i.e., $Excapprox[niσ,0]$ must nearly cancel U[n_iσ]), except when that orbital is stretched. But, as argued in Sec. II, that cancellation will be less perfect when the SIC orbital density is noded (Tables I and II) and even less perfect as the orbital density becomes more strongly noded. As the orbital density becomes more noded, its semilocal exchange-correlation energy becomes more negative and its PZ self-interaction correction of Eq. (3) becomes more positive, even when the orbital density is unstretched. Noded SIC orbitals can be avoided, even when the canonical or generalized Kohn-Sham orbitals are real, by making the unitary transformation matrix and the resulting localized SIC orbitals complex, as in Refs. 25 and 29.

Tables III–V confirm this expectation for the atomization energies of molecules in the small representative set AE6,³⁷ for the barrier heights of chemical reactions in the small representative set BH6,³⁷ and for the 55 molecular formation energies of the G2-1 set.³⁸ For PBE and SCAN, the PZ self-interaction-corrected results are more accurate when complex SIC orbitals are used. These results for PBE-SIC are in agreement with those of earlier studies,^23–27 which also noted that only the real localized orbitals are noded. But the explanation in the previous paragraph for the better results from complex orbitals is ours.

The SCAN and SCAN-SIC atomization energies of molecules (Tables III and V) illustrate the magnitude of the problem that remains to be solved. Without fitting to any bonded system, SCAN yields errors much smaller than those of LSDA or PBE. But SCAN-SIC yields errors that are bigger than those of SCAN by a factor of four or five for real SIC orbitals and three or four for complex SIC orbitals.

Tables III–V also agree with earlier studies^23–30 showing that PZ-SIC with real localized orbitals, applied to PBE (and in our study to SCAN), improves barriers but often worsens atomization energies. Use of complex orbitals solves only part of the problem by reducing but not eliminating the problematic lobed character of the orbital densities. (See Fig. 2 of Ref. 30.) The best SIC atomization energies in Table III are those of PBE-CSIC, as in Table I of Ref. 24. Perhaps significantly, PBE-CSIC also yields the most accurate SIC total energies for large-Z atoms.

PZ SIC imposes a new exact constraint (exactness for all one-electron densities) at the cost of other exact constraints or of appropriate norms such as exactness for the uniform electron gas and for large-Z atoms.¹⁹ This in turn results in errors that do not fully cancel out in the energy differences between free atoms and molecules. Needed is a generalization of PZ SIC that treats occupied localized orbitals in the PZ way only when they do not overlap with other occupied localized orbitals but (unlike earlier attempts to scale down the self-interaction correction in many-electron regions²⁸) retains the full Hartree part of the PZ SIC. Such a generalization is being attempted³⁹ to make SIC exact for all large-Z atoms and other slowly varying densities. While this seems to be achievable, restoring some of the other exact constraints satisfied by the nonempirical semilocal functionals may be harder to achieve with a fully nonlocal approximation.

PZ SIC applied to a semilocal functional like the PBE GGA or the SCAN meta-GGA removes troublesome self-interaction errors that manifest most strongly for stretched bonds. But it introduces other errors, which manifest significantly in the atomization energies of molecules. Some of these are orbital-density nodality errors that can be removed by replacing real localized orbitals by complex ones (Tables III–V). Removing the nodes of the localized orbital densities reduces but does not eliminate their problematic lobed one-electron structure. Errors remain because the PZ-corrected functional for the total electron density violates (at least for real Fermi-Löwdin orbitals)¹⁹ exact constraints or appropriate norms that are built into PBE and especially SCAN. Hopefully a generalized PZ SIC³⁹ will be able to recover some of these correct features of SCAN, thus providing a large self-interaction correction only where it is needed. SCAN, which is already self-correlation free, is a natural candidate for a generalized self-interaction correction. As a first step, we have presented here what may be the first results for the ungeneralized PZ SIC applied to SCAN. (For a recent alternative approach, see Ref. 40.)

The PZ self-interaction correction does indeed vanish¹⁷ when applied to the exact density functional for the exchange-correlation energy, but unfortunately this conclusion does not guarantee that, when applied to a good semilocal approximate functional, this correction will be acceptably small in situations where the semilocal functional is accurate. This is because the densities of occupied localized orbitals are typically more challenging to the semilocal functional than is the total electron density.

To make PZ SIC exact for any collection of separated one-electron densities, the right orbitals for the PZ SIC must span the space of occupied generalized Kohn-Sham or canonical orbitals. Using the often-delocalized generalized Kohn-Sham orbitals themselves would lead to a well-known size-extensivity problem.¹⁷ We have found here a second reason not to use canonical orbitals: They can be highly noded (e.g., the 2p_z orbital), with disastrous results for PZ SIC applied to semilocal functionals beyond LSDA (Table I). A better choice is to use real localized orbitals that are more-weakly noded, i.e., noded near the edges but not near the center, such as the real Fermi-Löwdin orbitals.²¹ Still better for the correction of semilocal functionals are complex localized orbitals with nodeless orbital densities.²⁹ Complex Fermi-Löwdin orbitals are possible⁴¹ and also presumably guarantee the size-extensivity of SIC.

A good semilocal functional like PBE or SCAN will, by our analysis of Sec. II, produce a positive PZ self-interaction correction in Eq. (3) from any highly stretched or strongly noded localized orbital. As shown in Table VI, even the total self-interaction corrections to the exchange-correlation energies for atoms, using the weakly noded real Fermi-Löwdin orbitals, tend to be positive for these functionals. Thus, the strong nodality of some Kohn-Sham orbitals in atoms and small molecules drives the energy-minimizing localized orbitals away from the highly noded Kohn-Sham orbitals and toward the more weakly noded Fermi-Löwdin orbitals. But, even so, the PBE-SIC and SCAN-SIC exchange-correlation energies with real SIC orbitals are not negative enough and are much less accurate than those of SCAN (Table VI).

PZ SIC reduces atomization energies of molecules (often by too much) and raises energy barriers to chemical reactions. We can now propose a tentative explanation for this although the explanation is restricted to the case of real localized SIC orbitals. In the ground state, real localized SIC orbitals that overlap must be noded to achieve orthogonality. As we pass from separated atoms to the molecules they form or from separated molecules to the transition states they form, the valence orbitals acquire more orbitals with which they overlap and, hence, become more noded. As the valence orbitals become more noded, their PZ self-interaction correction of Eq. (3) becomes more positive or less negative. Thus, PZ SIC raises the energy of a molecule relative to the energy of its separated atoms and raises the energy of a transition state relative to its separated molecules. PZ SIC also raises the energy of the transition state relative to its separated molecules through bond stretching.

The correlation energy as defined in the exact Kohn-Sham theory is numerically close to its quantum-chemistry definition,⁴² the deviation of the total energy from the Hartree-Fock value. In this sense, the transition states of chemical reactions are more strongly correlated than are the ground states of molecules in equilibrium since the Hartree-Fock energy barriers to chemical reactions are seriously too high. The barriers from semilocal functionals are seriously too low, suggesting that in this sense they “over-correlate” for stretched radicals, but the barriers are rather accurate when a self-interaction correction is applied. This and other results (e.g., Refs. 9 and 43) suggest that properly self-interaction-corrected semilocal functionals might usefully describe the energetics of strongly correlated systems. By properly self-interaction-corrected, we mean exact for all one-electron densities and reasonably accurate even for the fractional electron number, as in PZ SIC, but roughly as accurate as SCAN for atomization energies and equilibrium bond lengths (an accuracy that arises from the satisfaction of exact constraints and appropriate norms, with no fitting to bonded systems). The SCAN and SCAN-SIC results in Tables III, V, and VI show how far we still are from this goal. Hybrid functionals,⁴³ which mix an empirical fraction of exact exchange with a complementary fraction of semilocal exchange, go only part of the way toward this goal by reducing but not eliminating the self-interaction error.

It is not only atomization energies but also equilibrium bond lengths that are worsened by PZ SIC. For small molecules, LSDA-SIC bond lengths are shorter and less accurate than those of LSDA, and PBE-SIC bond lengths are shorter and less accurate than those of PBE, whether real²⁸ or complex⁴⁴ SIC orbitals are used. In addition, PZ SIC breaks molecular symmetries.²⁵ We plan to return to these problems in future work.

A by-product of this work is the observation that the exact Kohn-Sham theory, with the standard exact ground-state exchange-correlation energy functional, is exact for all excited states of one-electron systems, as well as for all ground-states, and that the Perdew-Zunger self-interaction correction¹⁷ (PZ SIC) perfectly imitates the exact Kohn-Sham theory in this special case of one-electron densities. For a many-electron excited state, of course, there is no exact Kohn-Sham theory, but removing the self-interaction error from an approximate functional could improve the approximate description of an excited state.

The supplementary material presents details of the ERKALE and PySCF/PyFLOSIC codes and detailed results for the AE6, BH6, and G2-1 molecular test sets.

The work of many of us (B.S., K.A.J., J.E.P., K.T., H.M., S.A., A.R., Y.Y., T.B., R.R.Z., and J.P.P.) was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC0018331 as part of the Computational Chemical Sciences Program. The work of C.S. was supported by the Center for Complex Materials from First Principles, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC0012575. The work of P.B. and K.W. was supported by the U.S. National Science Foundation under Grant No. DMR-1607868. The work of J.K. was supported by the German Research Foundation (DFG) under No. KO1924/9-1. S.L. acknowledges support from the Academy of Finland (Suomen Akatemia, Luonnontieteiden ja Tekniikan Tutkimuksen Toimikunta), Grant No. 311149. N.K.N. acknowledges support by the National Science Foundation under Grant No. DMR-1553022. Many of us acknowledge stimulating discussions with Mark R. Pederson. J.P.P. designed the work and wrote the first draft. The other authors contributed calculations, figures, tables, references, discussions, and revisions. This research includes calculations carried out on Temple University’s HPC resources and thus was supported in part by the National Science Foundation through major research instrumentation Grant No. 1625061 and by the Army Research Laboratory under Contract No. W911NF-16-2-0189.