Density functionals with spin-density accuracy for open shells

Density functional theory (DFT)^1,2 and spin-density functional theory (SDFT)^3,4 are two reformulations of the many-electron problem. Their computational advantage is replacing the solution of Schrödinger’s equation for the multi-dimensional ground state of the physical electronic system of interest by the solution of a much less demanding equation that yields the ground-state (total) density ρ (DFT) or the spin density (ρ^↑, ρ^↓) (SDFT) of the real system.⁵ 

Both the theories can, in principle, be applied equally well to study electronic systems when the external magnetic field, B, vanishes (B = 0). Then, obviously, the predicted values for any observable quantity must be the same between the two, formally exact, theories. Specifically, the ground-state (total) density and energy of the electronic system must be the same in both the theories.

Even though for B = 0, both exact theories can equally well describe an electronic system, this changes when density functional approximations (DFAs) are introduced. The flexibility of two basic variables in SDFT is thought to offer a better modeling of the exchange-correlation (xc) energy, E_xc[ρ^↑, ρ^↓], compared with just a single variable in DFT, E_xc[ρ]. As Parr and Yang write,⁶ in local and semi-local density functional approximations “the exchange-correlation energy of the electrons is approximated locally by the results for the homogeneous spin-compensated electron gas. Such a procedure is not appropriate for systems with unpaired electrons, such as open-shell molecules. A better description for such systems will be obtained through the use of the exchange-correlation energy of the homogeneous spin-polarized electron gas.” As a result, open-shell systems are, in practice, usually treated in SDFT rather than DFT.

The first aim of this Communication is to revisit the above widely held view, namely, that approximations in SDFT are inherently superior to those of DFT. As it will transpire, this view is a misconception, perhaps abetted by a qualitative error in the standard treatment of open shells with Kohn–Sham (KS)-DFT approximations. We shall explain the source of this error and show that correcting for it yields KS-DFT results that are (almost) as accurate as the corresponding KS-SDFT results.

Although this result goes against current wisdom, its corollary is positive for the DFT literature, in general. Take, for example, the xc derivative discontinuity, the shift of the KS-DFT xc potential when the number of electrons crosses an integer number of particles from N to N + 1. Obviously, either N or N + 1 is an odd number, and in the analysis of the derivative discontinuity, it is assumed that it is possible to calculate the xc potential for an odd number of electrons accurately, without qualitative errors. With the present work, this assumption is justified, even for approximate xc potentials.⁷ 

In (S)DFT, the value of any observable quantity becomes a functional of the (spin) density. So far, little is known about the functional dependence of most observable quantities, except the total energy and identically the total density (DFT) or the spin density (SDFT). For open-shell systems and for B = 0, the virtual KS systems of DFT and SDFT are not the same (despite sharing the same total density), and the single-particle KS equations of exact SDFT do not reduce to those of exact DFT. In DFT, the (unknown) density functional for the spin density is required to obtain the physical spin density.

The second aim of this Communication is to investigate the link between exact DFT and SDFT when B = 0. In general, a pair of external potentials is necessary for the formulation of SDFT (conjugate to the pair of spin densities), whereas when B = 0, there is only one external potential. So, either SDFT becomes singular in this limit or the KS-SDFT equations must reduce to a set of (yet unknown) single-particle equations of DFT. In Sec. IV of this Communication, we show that when B = 0, the KS-SDFT equations reduce to the equations of generalized KS-DFT (GKS-DFT), establishing a link between DFT and SDFT. In addition, the spin density of the GKS-DFT system emerges as the unknown density functional for the physical (real) spin density.

For simplicity, we shall discuss the special case of SDFT for a collinear magnetization density m = −μ_B(ρ^↑ − ρ^↓). We shall also restrict to densities that are pure-state v-representable (both for the interacting and for the KS system⁸), with ground states having an integer number of spin-up and spin-down electrons.⁹ Finally, we consider that the KS states are Slater determinants.

We return now to the first aim of our Communication. For B = 0, the main error in modeling open-shell electronic systems as unpolarized arises in the exchange energy functional, E_x[ρ] (DFT) vs E_x[ρ^↑, ρ^↓] (SDFT). When a local or semi-local approximate exchange energy expression is employed, depending explicitly only on the density ρ and its derivatives, it is essentially assumed that the spin-up and spin-down densities are equal to each other and to one-half of the total density. For open-shell systems, this amounts to mixing partly the spin-up ρ^↑ with the spin-down ρ^↓ densities in the exchange energy expression,

This mixing of the spin densities in the exchange energy leads to a spurious error that we call the ghost-exchange energy error, in analogy to the “ghost interaction” error of Ref. 10. We define it by the difference of the exchange energies with mixed spin densities from the reference exchange energy, where the spin densities are separate,

In DFT, the ghost-exchange error is not caused by the approximate or exact expression for the exchange energy but by the mixing of two spin densities when an open-shell system is treated as if it were closed-shell, considering half its electrons as spin-up and half as spin-down.

In the Secs. II and III, we show that open-shell systems in zero magnetic field can be modeled within DFT, without any mixing of spin-up and spin-down KS orbitals and avoiding the ghost-exchange error. Then, we explore the link between exact DFT and SDFT in the limit B = 0 in Sec. IV. Finally, in Sec. V, we discuss the famous spin-symmetry dilemma¹¹ in the context of our work.

Before proceeding, we note that there are various extensions of KS-DFT targeted for open-shell systems. For example, we mention the restricted open-shell KS (ROKS) method,^12,13 KS-DFT with complex, spin-restricted KS orbitals,¹⁴ and the recent progress in ensemble DFT for degenerate/excited states.^15–20 Methods with multi-determinant KS reference states, or using ensemble DFT, are, indeed, intrinsically advantageous for open-shell systems. Our approach is restricted to single-determinant references states and cannot be expected to compete in accuracy. Nevertheless, we maintain the practical advantage of KS-DFT, by employing its standard and popular approximations, while our aim has been to understand better the basic theory and to reveal the missing link between KS-DFT and KS-SDFT.

To continue, we consider the exact exchange energy in KS-DFT. It is given by the Fock expression in terms of the spin-orbitals of the KS Slater determinant,

where ρ^σ(r, x), σ = ↑, ↓, is the spin-σ, one-body reduced density matrix of KS orbitals.

The exact exchange energy functional separates in two disjoint terms,

in which the subsets of spin-up and spin-down KS orbitals do not mix. Hence, in order to avoid cross-exchange effects between opposite spin-electrons also in approximations, the approximate exchange energy density functionals must also satisfy equality (4).

Since modeling the approximate exchange energy in terms of the total density and its derivatives, as in the local density approximation (LDA) and in semi-local DFAs, always violates equality (4), we have to model the exchange energy using the DFT-KS spin density.²¹ This modeling is still within DFT (not SDFT) since in DFT, the KS spin density (the spin density of the KS determinant) is an implicit functional of the total density.

Equality (4) is satisfied by the local spin-density approximation (LSDA) and as far as we know all spin-dependent density functionals (DFs) for the exchange energy, approximate or the exact one. Therefore, we shall use expressions for the exchange energy that depend on the total density indirectly, or implicitly, via the KS spin density [spin-density functionals (SDFs)],

where $ExSDF$ must be the sum of two disjoint terms (4).

The acronyms iDF and SDF denote that the exchange energy on the lhs is an implicit density functional and that the exchange energy expression on the rhs depends on the KS spin density.

Equation (5) not only (re)defines for any DFA a new, ghost-exchange-error free, x energy functional, but it is also a relation satisfied by the exact x energy functionals of DFT and SDFT. In other words, for exact exchange: $Ex[ρ]=ExiDF[ρ]$⁠.

In order to keep the treatment of xc uniform, we decompose the correlation energy DF in an implicit DF, following (5), plus a remainder,

where ΔE_c[ρ] ≠ 0, in general, for open-shell systems. In the following, we omit ΔE_c (which makes the rest of the treatment approximate) and write $Ec[ρ]≈EciDF[ρ]$⁠. In (semi)-local DFAs, this redefinition of the correlation energy as an implicit DF yields a small accuracy gain, though this is typically a minor contribution compared to the ghost-exchange error. The whole exchange-correlation (xc) density functional is, thus, written as

Hence, the total energy DF is now given by

where U[ρ] is the Hartree energy and v_en(r) is the external (electron–nuclear) potential of the interacting system.

As the xc energy functional is an implicit DF, the xc potential, given by the functional derivative,

must be determined using the optimized effective potential (OEP) method,^22,23

with

where χ^σ(r, r′) is the spin-σ density–density response function.²⁴ 

The derivation of Eqs. (8)–(11) constitutes our first important result. These equations determine the ghost-exchange-error-corrected xc energy and potential for open-shell systems in (semi)-local DFAs. From (11), we obtain the DFT-KS xc potential v_xc, which is used in the usual way during the self-consistent KS calculation. It is given as the weighted sum of the spin-dependent xc potentials $vxcσ$⁠, with weighting factors being the spin-dependent response functions χ^σ.

For closed-shell systems, the solution of (11) reduces to the familiar KS xc functional derivative. For fully spin-polarized systems (e.g. spin-up, ρ^↑ = ρ, ρ^↓ = 0), the spin-down response function vanishes, and the solution of (11) reduces to the spin-up xc potential, $vxc[ρ]=vxc↑[ρ,0]$⁠. Even for partially spin-polarized, open-shell systems, the spin-up (i.e., the majority spin) xc potential is still an accurate approximation to the full solution of (11). First, this is because the latter must interpolate smoothly between the two limits just described, where $vxc=vxc↑$⁠. Second, because the weighting of the spin-σ response function in (11) depends roughly on the inverse of the HOMO–LUMO gap, Δ^σ, of the spin-σ KS orbitals, and since Δ^↑ < Δ^↓, the spin-up response function dominates. Another meaningful approximation for the solution is the weighted average, $vxc≈(Δ↓vxc↑+Δ↑vxc↓)/(Δ↑+Δ↓)$⁠. In the supplementary material, we compare these two approximations with the full solution of the OEP equation (11).

We have implemented Eqs. (5)–(11) in the Gaussian basis set code HIPPO.²⁵ The solution of finite-basis OEP equations hides non-analytic behavior,²⁶ which leads to oscillations of the potential near the nuclei.^27,28 We overcome these issues routinely following the method in Refs. 7 and 26 (please see Refs. 7 and 29 for details of the computational implementation). Unless otherwise stated, all results use cc-pVTZ bases for the KS orbitals and uncontracted cc-pVDZ auxiliary bases^30,31 for the expansion of the “electron repulsion” or “screening” density.^7,32,33

Ground-state energies calculated with our iDF KS method compare favorably with those from SDFT-KS. This is demonstrated in Table I for some atoms and molecules at their equilibrium geometries, using the L(S)DA functional. We emphasize again that we seek only the ground-state energy in our approach, and thus, where a state is referred to as “doublet,” this simply means

and likewise for the triplet and so on. We also note the effect of the ghost-exchange (and correlation) error, which is particularly pronounced for triplet states; we shall now explore some examples in which this error significantly affects the results.

In Fig. 1, we have calculated the energy dissociation curve for the H₂ molecule with the L(S)DA functional. As is well-known in the literature³⁴ and can be seen in Fig. 1, the standard restricted solution yields a qualitatively incorrect dissociation curve; however, in our new method, once the bond is stretched enough, the triplet state becomes lower in energy than the singlet, and the energy tends to the correct limit. This transition does not occur when the ghost-exchange energy is present: in this case, the triplet energy is higher than the singlet at all bond distances.

We see a similar picture emerge for the stretched OH radical with the Perdew–Burke–Ernzehof (PBE)³⁵ functional in Fig. 2. In this case, the energy dissociation curve is again qualitatively inaccurate for the solution contaminated with ghost-exchange errors; in our method that removes the ghost-exchange error, the energy correctly becomes equal to the energy of the two atoms in the infinitely stretched limit. In both the examples, the transition region is interesting: the unrestricted solution yields a smooth dissociation curve, whereas in our single-determinant method, the transition is abrupt. Of course, our method does not correct the ubiquitous inability of typical semi-local functionals to capture static correlation effects.³⁶ Naturally, using a spin-restricted formalism mandates constraints absent in an unrestricted scheme: for example, the HOMO and LUMO KS energy levels in a “doublet” state are by definition degenerate, though this is not typically the case for the corresponding energy levels in an unrestricted scheme.

We now turn our attention to the limit of SDFT when B = 0 and compare with DFT. We start with the universal internal energy density functional,

The minimizing state Ψ_ρ depends on the total density ρ. When ρ is the ground-state density of an interacting electronic system bound by a local potential v (i.e., when ρ is pure-state interacting v-representable), Ψ_ρ is the corresponding ground state. We remind the reader that we restricted our study to densities that are pure-state v-representable, both for the interacting and the KS systems.⁸ The numbers of spin-up/down electrons N^↑, N^↓ are good quantum numbers for Hamiltonians with collinear B, including B = 0. It follows that the ground state Ψ_ρ can always be selected to have integer N^↑, N^↓. In DFT-KS theory, the KS reference state Φ_ρ is obtained from a similar minimization to (15) for the non-interacting kinetic energy density functional, $Ts[ρ]=minΦ→ρ⟨Φ|T̂|Φ⟩$⁠. Under the same assumptions, the KS ground state Φ_ρ can also be chosen to have integer N^↑, N^↓ and be a Slater determinant.

To proceed, we follow Levy³⁷ and separate the minimization in (15) into two separate minimizations with the same minimum,

The inner minimization, performed first, is over all states Ψ with spin density (ρ^↑, ρ^↓). The outer minimization is over all spin densities that add up to the total density ρ. The first minimization inside the brackets on the rhs defines SDFT’s internal energy functional; we conclude that the DFT internal energy functional can be obtained from the SDFT functional with an extra optimization over spin densities with the same ρ (see also Appendix B in Ref. 11),

The two minimizations in (16) have the same overall minimum as the minimization in (15). Therefore, the minimizing spin density, $(ρρ↑,ρρ↓)$⁠, is equal to the spin density of Ψ_ρ in (15) and, hence, equal to the spin density of the real interacting system, even though the minimization (17) is at fixed total density (not fixed spin density), i.e., within DFT, not SDFT.

The minimization in (17) is worth investigating further. We invoke SDFT’s KS system with spin density (ρ^↑, ρ^↓) to expand the internal energy functional in the usual way, F[ρ^↑, ρ^↓] = T_s[ρ^↑, ρ^↓] + E_xc[ρ^↑, ρ^↓] + U[ρ^↑ + ρ^↓]. Dropping the superscript SDF from the functionals of SDFT, we obtain

Since every SDFT-KS system has a different spin density, the minimization effectively searches over all SDFT-KS systems with common total density ρ and returns that with the correct spin density $(ρρ↑,ρρ↓)$⁠. The KS Slater determinant state is $Φρρ↑,ρρ↓$⁠. The minimizing SDFT-KS system depends only on the total density ρ. The corresponding SDFT-KS potential (functional of the spin density) at the specific spin density $(ρρ↑,ρρ↓)$ is also an implicit functional of ρ. In the supplementary material, we show how this unrestricted KS potential emerges directly from the minimization in (18).

The SDFT-KS system that minimizes (18) defines a new non-interacting system in DFT (when B = 0) represented by the state $Φρρ↑,ρρ↓$⁠. We call it the generalized KS (GKS) system, as its derivation is analogous to the well-known GKS scheme in the literature.³⁸ However, unlike typical GKS schemes, which employ a nonlocal potential term, the potential in the present DFT-GKS scheme is local and spin-dependent.

We conclude that the SDFT-KS equations reduce to the DFT-GKS equations for B = 0. The DFT-GKS Slater determinant state $Φρρ↑,ρρ↓$ gives not only the true total density but also the true (i.e., physical or observable) spin density. The elusive exact density functional for the spin density is the spin density of the DFT-GKS system.

One of the oft-cited issues with SDFT in the absence of an external magnetic field is the spin-symmetry dilemma, so-called because SDFT approximate results yield either accurate total energies but with a poor (broken symmetry) prediction for the spin density or an accurate prediction for the spin density with poor total energies.^11,39,40 The dilemma lies in that the spin density is the key quantity and hence SDFT (exact) results are supposed to yield both the exact total energy and the exact spin density. For example, H₂ should dissociate into two hydrogen atoms with zero net magnetization (difference of spin densities): indeed, the lowest energy solution of local and semi-local approximations in SDFT gets the energy right, but it wrongly dissociates the molecule into two atoms of opposite spin (broken symmetry solution). There is also a self-consistent solution where the spin densities are correctly equal to each other, but the total energy is too high.

In their famous paper,¹¹ Perdew, Savin, and Burke “escape” the spin-symmetry dilemma, by postulating a “nearly exact alternative theory,” where the spin densities obtained from (semi)-local functionals “are not physical spin densities, but are only intermediate objects such as the Kohn–Sham orbitals.” Instead, the on-top electron pair density P(r, r), and the total density, become the two fundamental variables of the new theory. As an aside, we note that P(r, r′) is itself a useful quantity and has been used to develop improved LSD approximations for multi-determinant reference states.⁴¹ We share the principle of this Communication, namely, that the spin density is no longer the basic quantity in the absence of a magnetic field, but in our approach, we do not consider any alternative theory: our (ghost-exchange-error-free) theory yields the correct ground-state energy and the correct total density in the dissociation limit, and thus, the spin-symmetry dilemma is weakened as a fundamental problem, since the KS spin density is not expected equal to the observable spin density.

However, we also propose that the DFT-GKS spin density is the exact density functional of the observable spin density; knowledge of the (exact) SDFT xc functional E_xc[ρ^↑, ρ^↓] is necessary to obtain from the density the (exact) DFT-GKS spin density. The value of the latter is the same as the SDFT-KS spin density,⁴² which, in local and semi-local approximations, demonstrates the spin-symmetry dilemma. We believe the explanation is that the approximate DFT-GKS spin density is not a sufficiently accurate density functional for the real spin density. According to Cohen, Mori-Sánchez, and Yang^34,36 on the fractional spin error, the crux of the problem is the degeneracy of different Slater determinants with the same energy and total density but different spin densities. This kind of strong correlation cannot be captured with local or semi-local functionals.

We have addressed two related conceptual or fundamental theoretical questions in DFT and spin-DFT. The first is the challenge how to apply the KS equations of DFT for open shells, avoiding a serious qualitative error, which we called the “ghost-exchange error.” Associated with this problem is the widespread belief that spin-DFT approximations are inherently more flexible and, hence, can model more accurately the xc energy for open shells than the corresponding approximations of DFT. In the literature, this common view seemed plausible, almost obvious. We demonstrate that actually, density functionals can be as accurate as spin-density functionals, provided the ghost-exchange error is corrected.

For zero external magnetic field (B = 0), both DFT and spin-DFT can obviously be applied to study open-shell electronic systems, and both exact theories should give identical results for any observable quantity. Then, intuitively, one would expect in the limit B = 0 the single-particle KS equations of spin-DFT to reduce to a set of single-particle equations of DFT. We demonstrate that this intuition is correct and that the SDFT equations reduce to the generalized KS equations of DFT.

See the supplementary material for the comparison of the full solution of the OEP equation (11) with two simple approximations and an alternative derivation of the DFT-GKS system.

The authors thank N. Lathiotakis for help with the code HIPPO. The authors acknowledge enlightening discussions with K. Burke, E. K. U. Gross, and N. Lathiotakis. N.G. acknowledges financial support from The Leverhulme Trust, through a Research Project (Grant No. RPG-2016-005). B.P. acknowledges financial support from the EPSRC (Grant No. EP/R513039/1).

During the final stage of the project, T.C. was funded by the Center for Advanced Systems Understanding (CASUS), which is financed by Germany’s Federal Ministry of Education and Research (BMBF) and by the Saxon Ministry for Science, Culture and Tourism (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament.

The authors state that they have no conflict of interest to disclose.

The data that support the findings of this study are openly available in Rodare at http://doi.org/10.14278/rodare.1404, Ref. 43.