Fractional spins and static correlation error in density functional theory
FFractional Spins and Static Correlation Error in Density Functional Theory
Aron J. Cohen, Paula Mori-S´anchez and Weitao Yang
Department of Chemistry, Duke University, Durham, North Carolina 27708 (Dated: October 3, 2018)Electronic states with fractional spins arise in systems with large static correlation (stronglycorrelated systems). Such fractional-spin states are shown to be ensembles of degenerate groundstates with normal spins. It is proven here that the energy of the exact functional for fractional-spinstates is a constant, equal to the energy of the comprising degenerate pure spin states. Dramaticdeviations from this exact constancy condition exist with all approximate functionals, leading tolarge static correlation errors for strongly correlated systems, such as chemical bond dissociationand band structure of Mott insulators. This is demonstrated with numerical calculations for severalmolecular systems. Approximating the constancy behavior for fractional spins should be a majoraim in functional constructions and should open the frontier for DFT to describe strongly correlatedsystems. The key results are also shown to apply in reduced density-matrix functional theory.
PACS numbers: 31.10.+z,31.15.E-,31.15.eg
Density functional theory (DFT) [1, 2] is a rigorous ap-proach for describing the ground state of any electronicsystem. The success or failure of DFT is based on thequality of the density functional approximation (DFA).One of the dramatic failures is in strongly correlated sys-tems, characterized by the presence of degeneracy or neardegeneracy [3], having large static correlation. The sim-plest example is the dissociation of H molecule [4, 5] forwhich commonly used DFAs over-estimate the energy bymore than 50 kcal/mol. Closedly related are the bandstructure of Mott insulators [6] which are described asmetallic by known DFAs and problems in describing su-perconducting cuprates [7].The improvement of the DFA is, therefore, a majorgoal that critically depends on the underlying theoreti-cal construction. One of the most useful developmentsis the extension of DFT to fractional charges developedby Perdew et. al. [8] in a grand canonical ensemble,which was also established later in a pure state formula-tion [9]. For a system with fractional charges, the exactenergy is a straight line interpolating the energies of theinteger electron systems. The violation of this exact con-dition leads to two types of errors [10], the delocalizationerror (DE) of most functionals like LDA, GGA and hy-brids [11, 12, 13, 14], and the localization error (LE) ofthe Hartree-Fock functional. DE captures the tendencyof commonly used functionals to bias toward a delocal-ized description of electrons with widespread implicationsfrom molecular reactions[15] to the band-gap of solids.Addressing this error resulted in the construction of theMCY3 and rCAMB3LYP functionals [16], that correctmany of the errors of previous functionals. In particularthey correctly predict the ionization energy and electronaffinity from their single-electron eigenvalues, and hencethe energy gaps in molecules [17].In this Letter, we make an extension of DFT to frac-tional spin systems and prove that the exact energy func-tional of a fractional spin state is that of the comprisingdegenerate normal spin states. We show that states withfractional spins arise in systems with large static correla- tion (strongly correlated systems) and that the dramaticdeviation from the proven exact condition accounts forlarge static correlation errors. We also introduce a quan-titative measure for static correlation error (SCE).Our starting point is the exact result for an ensem-ble of degenerate densities derived by Yang, Zhang andAyers (YZA) [9]: For a N -electron system in the exter-nal potential v ( r ) that has g -fold degenerate orthogonalground state wavefunctions (Φ i , i = 1 , , . . . , g ) with cor-responding densities ( ρ i , i = 1 , , . . . , g ) and ground stateenergy E v ( N ), the ensemble density is ρ = (cid:80) gi =1 C i ρ i , where 0 ≤ C i ≤ (cid:80) gi =1 C i = 1 . The exact energyfunctional satisfies the following equation E v (cid:34) g (cid:88) i =1 C i ρ i (cid:35) = E v [ ρ i ] = E v ( N ) , (1)if E v ( N ) ≤ ( E v ( N + 1) + E v ( N − / . Note that in thederivation of Eq. (1) only pure states were used and theensemble densities appear in the limit of large separationof fragments [9].We now examine the application of Eq. (1) to frac-tional spin systems. Consider a N -electron atom ormolecule that is a doublet, with total spin S = , forexample, the H atom with N = 1. It has two degener-ate spin states labeled with the z -component of the spin m s = and m s = − . Now we construct ρ fs ( S, γ ), anensemble density with fractional spins as ρ fs (cid:0) , γ (cid:1) = (cid:18)
12 + γ (cid:19) ρ ( , ) + (cid:18) − γ (cid:19) ρ ( , − ) , (2)where ρ ( S, m s ) is the ground state density with m s and γ ( − ≤ γ ≤ ) is the net z -component of the spin in thefractional-spin state. ρ fs ( S, γ ) represents many fractionalspin states. In particular, γ = 0 represents a state thathas half a spin-up electron and half a spin-down electronoccupying the same spatial orbital, its total spin den-sity being equal to zero everywhere. Applying the YZArelation (Eq. (1)) leads to E v (cid:2) ρ fs (cid:0) , γ (cid:1)(cid:3) = E v (cid:2) ρ ( , ) (cid:3) = E v (cid:2) ρ ( , − ) (cid:3) , (3) a r X i v : . [ c ond - m a t . o t h e r] M a y -150-100-50 0 50 100 150 200 1 2 3 4 5 6 7 8 9 D E kc a l / m o l R/AngstromH -0.25 0 0.25 0.5 g H-0.5 HFB3LYPLDA
Figure 1: Binding curve of H calculated with spin-restrictedKS and fractional spins of H atom calculated with spin-unrestricted KS (multiplied by 2). γ = 0 is a H atom withhalf an α electron and half a β electron, which is the disso-cation limit of H . All calculations are self-consistent using acc-pVQZ basis set. which shows that for the exact density functional, all thefractional-spin densities ρ fs (cid:0) , γ (cid:1) have the same degen-erate energy. However, all known DFAs fail dramaticallyand give much too high an energy for E v (cid:2) ρ fs (cid:0) , γ (cid:1)(cid:3) ,as illustrated in the right-hand side of Fig. 1. In theself-consistent spin-unrestricted Kohn-Sham (KS) calcu-lations, we use the spin densities ρ σ fs ( S, γ ) ( σ = α, β ),which are represented by a noninteracting system withfractional occupation numbers, ρ σ fs ( S, γ ) =
HOMO (cid:88) i n σi | φ σi | , (4)where only the highest occupied molecular orbital(HOMO) for each spin is fractionally occupied, with n α HOMO = + γ and n β HOMO = − γ. At γ = 0, theground state is spin-unpolarized with ρ α fs ( ,
0) = ρ β fs ( , E v [ ρ fs ] = min ˜ ρ fs E v [˜ ρ fs ] , (5)where the domain of variation for ˜ ρ fs is all the ensembledensities constructed from Eq. (2) but with arbitraryspin densities ρ ( , ) and ρ ( , − ). We have also as-sumed that such ensemble densities can be representedby a noninteracting system with fractional occupation.Details of the proof for Eq. (5) can be found in Ref. [18].This formalism is particularly interesting because thefractional spin density ρ σ fs ( ,
0) describes the dissociationlimit of a single chemical bond. For example, at the disso-ciation limit of the H molecule, a singlet system ( S = 0) is obtained, which consists of two fractional spin H atomsseparated by a large distance. This system can be prop-erly described by multi-configurational wave functionmethods. However in DFT, spin-restricted KS calcula-tions having the correct spin state ( S = 0), give muchtoo high an energy, with DFAs. The over-estimation inthe energy for the dissociation of H matches exactly theover-estimation for the H atom with fractional spin den-sity ρ σ fs ( , molecule from the ground state unre-stricted atoms (with integer spins i.e. one alpha electronand zero beta electrons or vice versa, corresponding to ρ σ fs ( , ± )) and the right-hand side shows the differencein energy of the H atom with fractional spins, ρ σ fs ( , γ ),from the energy of the same ground state unrestrictedatom (multiplied by two for direct comparsion with thebinding curve). The energy should be constant with thechange in γ but all the energy functionals have a verylarge error, ranging from 30 kcal/mol to 170 kcal/molfor the midpoint, γ = 0. HF has the largest error andLDA has the smallest error, but both functionals over-estimate the energy for fractional spin states. B3LYP,as expected, has a behavior inbetween LDA and HF.Other functionals also suffer from large errors [18], withGGA functionals performing roughly the same as LDAand other hybrid functionals somewhere inbetween LDAand HF (this also includes coulomb attenuated function-als with reduced DE). This suggests that the calculationof strongly correlated systems, where this error is impor-tant, will qualitatively fail if any of the above functionalsare used. There are many attempts in the literature tocircumvent this error, for example breaking the spin sym-metry, which gives reasonable energies but wrong spindensities.Our discussion for the single bond dissociation can beextended to multiple bond dissociation. Using the nota-tion in Eq. (2), for a system with total spin S , we canconstruct the fractional spin density ρ fs ( S, γ ) from thetwo degenerate spin states with maximum | m s | = S ; ρ fs ( S, γ ) = (cid:18)
12 + γ S (cid:19) ρ ( S, S ) + (cid:18) − γ S (cid:19) ρ ( S, − S ) , (6)where − S ≤ γ ≤ S . Applying the YZA relation leads to E v [ ρ fs ( S, γ )] = E v [ ρ ( S, S ))] = E v [ ρ ( S, − S )] . (7)As in the case of S = where the fractional spin state ρ fs ( ,
0) describes the dissociation limit of a single chem-ical bond, the fractional spin state ρ fs ( S,
0) describes thedissociation limit of a multiple chemical bond. This isdemonstrated in Fig. 2 for the dissociation of a doublebond, C , a triple bond, N , and a sextuple bond, Cr ,into two S = 1, S = and S = 3 fractional spin atoms re-spectively. For the molecules we perform spin-restrictedKS calculations and show the binding curve with respect -200-100 0 100 200 300 400 D E kc a l / m o l C C HFB3LYPLDA-300-200-100 0 100 200 300 400 500 600 700 D E kc a l / m o l N N-200 0 200 400 600 800 1000 1200 1400 2 3 4 5 6 7 8 9 D E kc a l / m o l R/AngstromCr -0.25 0 0.25 0.5 g /2SCr-0.5 Figure 2: The same as Fig. 1 but for top, C (double bond),middle, N (triple bond) and bottom, Cr (sextuple bond). to the spin-unrestricted ground-state atoms calculatedwith no fractional spin ( m s = S ). On the right handside of Fig. 2 we show spin-unrestricted calculations onthe atoms with fractional spins, also relative to the nor-mal ( m s = S ) spin-unrestricted atoms. The expressionfor the density, ρ fs ( S, γ ), is exactly the same as Eq. (4)but with fractional occupation, n α HOMO = + γ S and n βHOMO = − γ S , for the top 2 S multiple HOMOs (e.g.the N atom, ρ fs ( ,
0) has half an α electron and half a β electron in the top three 2 p orbitals). To compare tomolecular dissociation the two densities mixed in Eq. (7)must have the same symmetry. All DFAs violate the con-stancy condition, Eq. (7). The over-estimation in the en-ergy for molecular dissociation matches exactly the over-estimation for the dissociating atoms with fractional-spindensity ρ σ fs ( S, E v [ ρ fs ( S, − E v [ ρ ( S, S ))] . (8)Static correlation can be described with the use of a fewSlater determinants for small molecules. However, forlarge and bulk systems, this becomes impractical. It isnow clear that SCE is an inconsistency in the commonlyused DFAs.The errors are massive and increase with the numberof bonds. It is also very significant to see that the er-ror at the dissociation limit can already dominate thebehavior close to the bonding region, making the limit-ing behavior analysis of E [ ρ fs [ S,
0] broadly relevant. ForCr , SCE make HF and B3LYP fail to describe boundmolecules and LDA has a very small range of bonding.Note that these cases are not only challenges for DFTbut also for single reference wave-function methods. Thecases considered here are homonuclear diatomics but thesame arguments apply to the dissociation of heteronu-clear diatomics and more complicated molecules.For the fractional-spin states ρ fs ( S, γ ), we have only ex-plored the consequence of the two-state ensemble whichleads to an understanding of static correlation. Thereare, however, more general fractional-spin states: ρ fs ( S, { C m s } ) = S (cid:88) m s = − S C m s ρ ( S, m s ) , (9)where 0 ≤ C m s ≤ (cid:80) Sm s = − S C m s = 1 . Based on theYZA relation the fractional spin constancy relation is E v [ ρ fs ( S, { C m s } )] = E v [ ρ ( S, m s ))] , (10)which will also have important consequence for moleculesand solids. What may hinder the exploration of Eq. (9)is the difficulty with which DFT deals with ρ ( S, m s ) for | m s | < S . Usually only the state ρ ( S, S ) is calculated, asit can be constructed easily from a KS determinant. Itis difficult, in general, to construct a KS determinant forother states ρ ( S, m s ) with | m s | < S .It is also possible to have fractional-spin states arisingfrom an ensemble of states which are degenerate because D E kc a l / m o l d LDAB3LYPHF
Figure 3: B atom with fractional alpha spins, δ = corre-sponds to the spherical B atom. of other symmetries (e.g. spatial): ρ fs ( S, { C i,m s } ) = g (cid:88) i =1 S (cid:88) m s = − S C i,m s ρ i ( S, m s ) , (11)for a g -fold degenerate system. The exact energy func-tional gives constant energy for all { C i,m s } . For examplea spherical atom density given by (cid:80) Li = − L L +1 ρ i ( S, S )has alpha fractional spin occupation of the spatially de-generate states. If we consider the case of the B atom,which is three-fold degenerate, the lowest energy statepredicted by DFAs is given by a non-spherical den-sity. We now examine the energy of ρ fs ( S, { C δ } ) = δ [ ρ ( S, S ) + ρ ( S, S )] + (1 − δ ) ρ ( S, S ) such that at ρ fs ( S, { C } ) corresponds to the normal ground statenon-spherical atom and ρ fs ( S, { C } ) corresponds to thespherical B atom. The energy of the fractional spin statesrelative to the ground state of the B atom is plotted inFig. 3, and shows again the violation of the constancyrelation for DFAs. The error of HF is similar to the casefor ρ fs ( S, γ ) but with pure DFT the error is much smaller(only ≈ − E [ ρ fs ] to offer a new challenge for E xc [ ρ ]and open new frontiers of DFT for strongly correlatedsystems.This work has been supported by the National ScienceFoundation. Appendix A: FRACTIONAL SPIN VARIATIONALPRINCIPLE AND ITS GENERALIZATION
Our starting point is the exact result for an ensem-ble of degenerate densities derived by Yang, Zhang andAyers (YZA) [9]: For a N -electron system in the exter-nal potential v ( r ) that has g -fold degenerate orthogonalground state wavefunctions (Φ i , i = 1 , , . . . , g ) with cor-responding densities ( ρ i , i = 1 , , . . . , g ) and ground stateenergy E v ( N ), the ensemble density is ρ = g (cid:88) i =1 C i ρ i , (A1)where 0 ≤ C i ≤ (cid:80) gi =1 C i = 1 . The exact energyfunctional satisfies the following equation E v (cid:34) g (cid:88) i =1 C i ρ i (cid:35) = E v [ ρ i ] = E v ( N ) , (A2)if E v ( N ) ≤ ( E v ( N + 1) + E v ( N − / . In deriving Eq. (A2), YZA used only pure statesand the three requirements on the density functional:(1) correct for each degenerate orthogonal ground state( ρ i , i = 1 , , . . . , g ), (2) translationally invariant and (3)size-consistent. The ensemble densities thus appear inthe limit of large separation of fragments [9]. In thissense, we see that with density functional theory (DFT)and density matrix functional theory, we are forced todefine functionals of ensemble densities.To construct the functional of ensemble densities of thetype of Eq. (A1), we consider the following trial ensembledensity matrix ˜Γ = g (cid:88) i =1 C i ˜Γ i , (A3)where 0 ≤ C i ≤ (cid:80) gi =1 C i = 1 . The density matrix˜Γ i is for a pure state corresponding to i th degeneratestate (in terms of spin, symmetry), but not neccessarilythe ground state. Then we can use the constrained search[19, 20] and define the functional for the trial ensembledensity ˜ ρ as F [˜ ρ ] = min ˜Γ → ˜ ρ Tr (cid:16) ˜Γ( ˆT + ˆV ee ) (cid:17) . (A4)The total energy functional is then E v [˜ ρ ] = F [˜ ρ ] + (cid:90) ˜ ρ ( r ) v ( r ) d r . (A5)The minimum of E v [˜ ρ ] is the ground state energy E v ( N ) = min ˜ ρ E v [˜ ρ ] , (A6)independent of the mixing coefficients { C i } , because thetrial ensemble density matrix, Eq. (A3), cannot havean energy lower than the ground state energy. This isthe general variational principle used in the text for theparticular fractional spin systems.While it may appear as a direct consequence of en-semble DFT, the variational principle of Eq. (A6), un-like a general ensemble DFT theory, connects directly tonormal pure-state DFT calculations without ensembles.The key is the following: the particular trial ensembledensity of Eq. (A1), which consists of densities from or-thogonal degenerate ground states, arises directly as thedissociation limit of normal pure systems [9]. Such purestates are calculated with normal pure-state DFT withone KS determinant, making the YZA analysis and ourvariational principle of Eq. (A6) directly relevant.In carrying out corresponding self-consistent KS cal-culations, we have also assumed that such trial ensembledensities ˜ ρ can be represented by a noninteracting sys-tems with fractional occupation. This parallels the devel-opment in the self-consistent KS calculations of fractionalcharge systems [15] where the fractional charge ensembleis also represented by a noninteracting system with frac-tional occupation.While the validity of such a KS representation has notbeen mathematically established here, just as in the caseof the KS representation for normal pure state densities,the most important justification is that such fractionaloccupation KS calculations reveal important features inthe functionals for fractional spin ensembles as reportedin this work and fractional charge ensembles as reportedearlier [15]. These important failures of density function-als are observed in normal DFT calculations without us-ing fractional charge nor spin. When they are analyzed interms of the fractional occupation KS calculations, they become transparent as a clear violation of the fundamen-tal equalities: linearity in the fractional charge case andconstancy in the fractional spin case. Appendix B: RESULTS FOR THE H ATOM
The performance of several different functionals isshown in Fig. 1 for the fractional spin H atom. Thesmallest error is seen for LDA and GGA functionals[11, 12, 13] with hybrid functionals B3LYP [14], MCY2[21], B05 [5] and PBE0[22] all having very similar perfor-mance despite their varied forms. Functionals with im-proved behaviour on delocalization error, such as MCY3and rCAMB3LYP [16], which have a much straighter line D E kc a l / m o l g -0.5 HFrCAMB3LYPMCY3CAMB3LYPPBE0B05MCY2B3LYPPBEBLYPLDA Figure 4: Spin-unrestricted KS calculations, using several dif-ferent functionals, on the H atom with fractional spins (mul-tiplied by 2). γ = 0 is a H atom with half an α electron andhalf a β electron, which is the dissocation limit of H . Allcalculations are self-consistent using a cc-pVQZ basis set. for fractional charges, have a very poor performance forthis fractional spin problem. Note also that functionalswhich are exact for the integer spin H atom such as HF,MCY2 and B05 have large errors for fractional spin. [1] W. Kohn and L. Sham, Phys. Rev. A , 1133 (1965).[2] R. G. Parr and W. Yang, Density-Functional Theoryof Atoms and Molecules (Oxford University Press, NewYork, 1989).[3] A. Savin, in
Recent Developments and Applications ofModern Density Functional Theory , edited by J. M. Sem-inario (Elsevier, Amsterdam, 1996), p. 327.[4] E. J. Baerends, Phys. Rev. Lett. , 133004 (2001).[5] A. D. Becke, J. Chem. Phys. , 064101 (2005).[6] P. Fazekas, Lecture notes on electron correlation andmagnetism in Series in Modern Condensed Matter Physics, vol. 5 (World Scientific, 1999).[7] J. K. Perry, J. Phys. Chem. A , 2438 (2000).[8] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz Jr.,Phys. Rev. Lett. , 1691 (1982).[9] W. Yang, Y. Zhang, and P. W. Ayers, Phys. Rev. Lett. , 5172 (2000).[10] P. Mori-S´anchez, A. J. Cohen, and W. Yang, Phys. Rev.Lett. , 146401 (2008).[11] A. D. Becke, Phys. Rev. A , 3098 (1988).[12] C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B , 785(1988). [13] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996).[14] A. D. Becke, J. Chem. Phys. , 5648 (1993).[15] P. Mori-S´anchez, A. J. Cohen, and W. Yang, J. Chem.Phys. , 201102 (2006).[16] A. J. Cohen, P. Mori-S´anchez, and W. Yang, J. Chem.Phys. , 191109 (2007).[17] A. J. Cohen, P. Mori-S´anchez, and W. Yang, Phys. Rev.B , 6062 (1979).[20] S. M. Valone, J. Chem. Phys. , 4653 (1980).[21] P. Mori-Sanchez, A. J. Cohen, and W. Yang, Journal ofChemical Physics , 091102 (2006).[22] C. Adamo and V. Barone, J. Chem. Phys.110