Importance of complex orbitals in calculating the self-interaction-corrected ground state of atoms
aa r X i v : . [ phy s i c s . a t m - c l u s ] A ug Importance of complex orbitals in calculating the self-interaction corrected groundstate of atoms
Simon Kl¨upfel, ∗ Peter Kl¨upfel, and Hannes J´onsson
Science Institute and Faculty of Science, VR-III, University of Iceland, Reykjav´ık, Iceland (Dated: July 27, 2018)The ground state of atoms from H to Ar was calculated using a self-interaction correction tolocal and gradient dependent density functionals. The correction can significantly improve thetotal energy and makes the orbital energies consistent with ionization energies. However, whenthe calculation is restricted to real orbitals, application of the self-interaction correction can givesignificantly higher total energy and worse results, as illustrated by the case of the Perdew-Burke-Ernzerhof gradient dependent functional. This illustrates the importance of using complex orbitalsfor systems described by orbital density dependent energy functionals.
PACS numbers: 31.15.xr, 31.15.E-
Density functional theory (DFT) of electronic systemshas become a widely used tool in calculations for solids,liquids, and molecules [1]. The most commonly used ap-proximation to the exact energy functional for extendedsystems is the so-called generalized gradient approxima-tion (GGA), where the functional only depends on the to-tal spin density and its gradient, an improvement on thelocal density approximation (LDA), where gradients arenot included. The kinetic energy is commonly evaluatedby introducing orthonormal orbitals, ϕ i ( r ), consistentwith a given total electron density, ρ ( r ) = P i | ϕ i ( r ) | . Inmany cases, DFT significantly improves the total energyover the Hartree-Fock (HF) method and gives acceptableaccuracy with smaller computational effort. However, anumber of shortcomings are also known: The bond en-ergy tends to be too large while the activation energyfor atomic rearrangements tends to be underestimated[2]. There is also a tendency to over delocalize spin den-sity, sometimes making localized electronic defects unsta-ble with respect to delocalization [3]. For finite systems,an ionization energy (IE) can be determined as the en-ergy difference between the ground states of charged andneutral species. These values often agree well with ex-periment, but ionization from deeper energy levels andionization from solids cannot be estimated this way. Un-like in HF theory, the energy associated with the orbitals(single particle eigenvalues) is in practice neither a goodnor a theoretically justified estimate of ionization energy.Even the exact DFT functional would give an estimateof only the first ionization energy.A leading source of error is the spurious self-interactionintroduced when the Hartree energy E H is evaluated fromthe total electron density ρ ( r ): E H [ ρ ] = 12 Z d r d r ′ ρ ( r ) ρ ( r ′ ) | r − r ′ | (1)If the orbital densities ρ i ( r ) = | ϕ i ( r ) | represent single ∗ [email protected] particle distributions, a more accurate expression is E ODDH [ ρ N ] = 12 X i = j Z d r d r ′ ρ i ( r ) ρ j ( r ′ ) | r − r ′ | . (2)Here, ρ N denotes the set of orbital densities { ρ , . . . , ρ N } corresponding to the set of occupied orbitals, { ϕ , . . . , ϕ N } , denoted by ϕ N . This expression forthe energy is explicitly orbital-density dependent(ODD). The difference between the two expressions isthe diagonal terms ( i = j ), which can be interpretedas the interaction of the electron in each orbital withitself. In Hartree-Fock calculations, where the Hartreeenergy is evaluated as in Eq. (1), the exchange energyincludes equally large self-interaction with opposite sign,so the self-interaction cancels out exactly. In LDA andGGA [collectively referred to as Kohn-Sham (KS) here],the exchange-correlation energy, E xc , is approximateand the cancellation is incomplete. Perdew and Zunger[4] proposed an orbital-based self-interaction correction(SIC) E SIC (cid:2) ρ N (cid:3) = E KS [ ρ ] − N X i =1 E SI [ ρ i ] (3)using E SI [ ρ i ] = E H [ ρ i ]+ E xc [ ρ i ]. Other estimates of SICcan be formulated [5, 6], but here we take Eq. (3) to bethe definition. Originally, SIC was proposed for LDA, butit can in principle be applied to other functionals. Theseare examples of a more general class of functionals, ODDfunctionals, where the Hartree energy is evaluated fromEq. (2).Variational optimization of orbital dependent function-als is typically carried out by adding the orthonormalityconstraints multiplied by Lagrange multipliers, λ ji , tothe energy functional to give S (cid:2) ρ N (cid:3) = E SIC (cid:2) ρ N (cid:3) − N X i,j =1 λ ji ( h ϕ i | ϕ j i − δ ij ) (4)and finding a stationary point with respect to variationof each orbital | ϕ i i and its complex conjugate. This givestwo sets of equations for the optimal orbitals [7–9], b H i | ϕ i i = N X i =1 λ ji | ϕ j i and λ = λ † , (5)with b H i | ϕ i i = δE SIC /δϕ i and λ ji = h ϕ j | b H i | ϕ i i . TheODD functional form leads to orbital-dependent Hamil-tonians, b H i . In the case of SIC, E SI is the orbital densitydependent part of the energy while E KS only depends onthe total spin density, so that b H i = b H KS [ ρ ] + ˆ v [ ρ i ] (6)where ˆ v [ ρ i ] = − δE SI [ ρ i ] /δρ i = − (ˆ v H [ ρ i ] + ˆ v x c [ ρ i ]). Atthe minimum, the constraint matrix λ is Hermitian andcan be diagonalized to give orbital energies ε i and corre-sponding eigenfunctions, the canonical orbitals, ψ N . Interms of these, condition (5) can be written as b H | ψ i i = ε i | ψ i i , ε i δ ij = ( W λ W † ) ij (7)with | ψ i i = P Nk =1 W ki | ϕ k i . The non-local operator b H isdefined in terms of the N optimal orbitals ϕ N and theircorresponding Hamiltonians b H i . In this way, the single-particle equations (5) can be decoupled to give traditionaleigenvalue equations [9]. The calculation is carried outusing two sets of orbitals, ϕ N and ψ N , while keepingtrack of the transformation matrix, W , relating the twosets. In HF and KS-DFT, the energy is independent ofthe unitary transformation and the optimal orbitals aretypically chosen to be the same as the canonical orbitals,i.e., W = .At intermediate steps of the variational optimization ofODD functionals, λ is in general not Hermitian, but canbe made Hermitian by finding the unitary transformationthat zeros the matrix κ = λ − λ † [7], where κ ij = Z d r ϕ ∗ i ( r ) ϕ j ( r ) ( v [ ρ i ]( r ) − v [ ρ j ]( r )) = 0 . (8)With a Hermitian λ at each iteration, Eq. (7) can beused during the minimization of the energy in a schemeanalogous to what is done in KS-DFT and HF.The results presented here were obtained with the pro-gram quantice [10] using the Cartesian representationof the Gaussian-type correlation-consistent polarized va-lence quadruple zeta (cc- pVQZ)[11] basis set and nu-merical grids of 75 radial shells of 302 points [12]. Theconvergence criterion in the optimization was a squaredresidual norm below 10 − eV . For LDA, Slater exchangeand the Perdew-Wang parameterization of correlation isused (SPW92)[13].Figure 1 compares the energy of the atoms H to Ar,calculated using various functionals, with accurate refer-ence values [14]. While the inclusion of gradients in theGGA type PBE functional [15] improves on the resultsobtained with LDA, the energy is still significantly toohigh and the error per electron tends to increase withthe size of the atom. SIC applied to LDA reduces the (b) PBEPBE0PBE-SICPBEPBE0PBE-SIC-0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 H He Li Be B C N O F Ne Na Mg Al Si P S Cl Ar ( E D FT - E r e f ) / N e l ( e V ) (a) LDALDA-SICLDALDA-SIC-1.5-1.0-0.5 0.0 0.5 1.0 1.5 2.0 2.5 ( E D FT - E r e f ) / N e l ( e V ) FIG. 1. (Color online) Energy of atoms H to Ar using LDA,LDA-SIC, PBE, PBE-SIC, and PBE0 functionals comparedwith accurate, non-relativistic estimates [14], normalized tothe number of electrons. Results obtained with complex or-bitals (solid symbols) and real orbitals (open symbols) are alsocompared. Vertical lines indicate the s , p , and p electronconfigurations. The gray shading at ± magnitude of the error but gives a strong overcorrection.The overcorrection also increases with the atomic num-ber. When SIC is applied to the gradient-dependent PBEfunctional, the error is reduced to ∼ ±
5% with experiment for theatoms H to Ar, both for PBE and PBE-SIC. However,in extended systems, subject to periodic boundaries, thecharged system can not be calculated rigorously, so the IEhas to be extracted from ground-state properties. In Fig.2, the energy eigenvalue of the highest occupied canoni-cal orbital is compared with the calculated IE using thefunctionals PBE, PBE0, and PBE-SIC. As is well known,for the commonly used LDA and GGA functionals, theeigenvalues obtained from KS-DFT do not give good es- - ε H O M O ( e V ) (b) exp.PBE PBE0 PBE-SICPBE PBE0 PBE-SIC-50-40-30-20-1001020 (- ε H O M O - ∆ E ) / ∆ E ( % ) (a) FIG. 2. (Color online) (a) Comparison of the eigenvalue of thehighest occupied orbital ε HOMO with ∆ E = E cation − E neutral .(b) Comparison of PBE-SIC eigenvalues with experimentalionization energy [19]. Open and solid symbols are used asin Fig. 1. PBE-SIC eigenvalues agree much better than PBEand PBE0 with both the calculated energy difference and ex-perimental data. timates of ionization energy. The IE estimates from PBEand PBE0 eigenvalues are too low, with errors of ∼ s and p type and the total den-sity is spherical. The optimal orbitals, however, dif-fer significantly in shape. The real orbitals are alignedin the traditional tetrahedral sp configuration and canbe constructed from the canonical orbitals as ϕ r = ( s + p x + p y + p z ) followed by three consecutive fourfoldimproper rotations about the z axis, ˆ S ( z ). The com-plex optimal orbitals are also, despite their uncommontetragonal configuration, sp -hybrid orbitals. Such a setcan be constructed from ϕ c = ( s + p x + p y + i p z ) andapplication of the rotations ˆ S ( z ). The shape of a realand complex orbital is compared in Fig. 4. The densityof the real orbital has axial symmetry and a nodal sur-face. The complex orbital has lower symmetry and lacksthe nodal surface since orthogonality is achieved by thephase of the complex expansion coefficients.The large increase in total energy that occurs when (a) (c)(b) FIG. 3. (Color online) Orbital densities of a Ne atom obtainedfrom a PBE-SIC calculation. (a) Electron density isosurfacesof the complex canonical orbitals show clear correspondencewith s and p orbitals. The real canonical orbitals have asimilar shape. Electron density isosurfaces from (b) real and(c) complex optimal valence orbitals shown from the top andside views for one spin component. The gray spheres showthe total density.FIG. 4. (Color online) Isosurface and contour plots of anoptimal valence orbital of neon obtained with PBE-SIC using(a) real and (b) complex expansion coefficients. The contourplots show the orbital density in three planes through thenucleus. The complex orbital has lower symmetry and nonodal surface. the SIC calculation is restricted to real orbitals was alsoobserved for other GGA-SIC functionals as well as forexchange-only SIC calculations. This effect can be ex-plained by the fundamentally different structure of thefunctional as compared to KS or HF. There, the energyis invariant with respect to unitary transformations of theorbitals so the full flexibility of complex expansion coef-ficients is not needed. The SIC energy, however, dependsexplicitly on the orbital densities. Nodal surfaces, whichare required for orthogonality of real orbitals, impose astrong constraint on the shape of the orbital densities(and their gradient, resulting in a stronger effect on theenergy for SIC applied to GGA functionals [6]). Thiscan be illustrated by a simple example: Given a basis setconsisting of two Cartesian p -type orbitals, ϕ c ( r ) = N xe − β r , ϕ c ( r ) = N ye − β r , (9a)a second set can be constructed as the complex sphericalrepresentation, ϕ s ( r ) = N √ re i φ e − β r , ϕ s ( r ) = N √ re − i φ e − β r , (9b) FIG. 5. Orbital densities obtained from the Cartesian (upperpanel) and spherical basis (lower panel) using real expansioncoefficients. The Cartesian basis only gives localized orbitalsbut allows for all orientations. The spherical basis makes atransition from delocalized to localized orbitals possible butis restricted to certain orientations of the localized orbitals. where r = p x + y . The two sets of orthogonal orbitalsaccessible by real linear combinations are defined by asingle parameter α : e ϕ x = cos( α ) ϕ x + sin( α ) ϕ x (10a) e ϕ x = − sin( α ) ϕ x + cos( α ) ϕ x , x ∈ { c, s } (10b)The orbital densities are illustrated in Fig. 5 for severalvalues of α . For the Cartesian basis (9a), the variationof α results in a rotation of the orbitals about the z axis,while for the spherical basis a transition from delocalizedto localized orbitals occurs. The SIC calculated from theCartesian basis is independent of α , but for the spherical basis the energy of the localized and delocalized orbitalswill differ. While the spherical basis (9b) can producedelocalized orbitals, it is incapable of describing localizedorbital densities pointing in the “diagonal” orientation ifreal linear combination coefficients are used. Only whencomplex coefficients are used, can both basis sets give thefull range of possibilities.While errors in the total energy of atoms can, to alarge extent, cancel out when the energy differences oftwo atomic configurations are calculated, as, for example,in calculations of bond energy, we find that the calculatedbinding energy in diatomic molecules is significantly af-fected both by the inclusion of SIC and then also by arestriction to real orbitals. The binding energy of N ispredicted by PBE to be ∼ , PBE predicts a binding energy that is1.02 eV too large. Here, PBE-SIC overcorrects and givesa value that is 0.54 eV too small, while a calculationrestricted to real orbitals gives an even larger overcorrec-tion, a binding energy that is 1.05 eV too small.The calculations presented here for atoms demonstratethat PBE-SIC gives substantial improvement in the totalenergy and physically meaningful orbital energies. TheSIC applied to KS functionals represents only a smallset of possible orbital density dependent functionals. Itis likely that other functionals of ODD form allow foran even more accurate modeling of the electronic groundstate providing a more flexible and computationally effi-cient alternative to hybrid functionals. The derivation ofa functional that makes optimal use of the more generalODD form appears to be a promising prospect. But, it isclear that an assessment of the quality of any ODD func-tional requires a minimization in the variational space ofcomplex orbitals. [1] W. Kohn, Rev. Mod. Phys. , 1253 (1999).[2] P. Nachtigall et al. , J. Chem. Phys. , 148 (1996).[3] G. Pacchioni et al. , Phys. Rev. B , 054102 (2000).[4] J.P. Perdew and A. Zunger, Phys. Rev. B , 5048(1981).[5] U. Lundin and O. Eriksson, Int. J. Quantum Chem. ,247 (2001); C. Legrand, E. Suraud and P.-G. Reinhard, J. Phys. B , 1115 (2002).[6] O.A. Vydrov and G.E. Scuseria, J. Chem. Phys. ,191101 (2006); O.A. Vydrov et al. , ibid. , 094108(2006).[7] M.R. Pederson, R. Heaton and C.C. Lin, J. Chem. Phys. , 1972 (1984); , 2688 (1985).[8] S. Goedecker, C.J. Umrigar, Phys. Rev. A , 1765(1997); A. Svane, Phys. Rev. B , 4275 (1996).[9] J. Messud et al. , Phys. Rev. Lett. , 096404 (2008);Ann. Phys. (N.Y.) , 955 (2009).[10] http://theochem.org/quantice/. [11] D.E. Woon and J.T. Dunning,
J. Chem. Phys. , 2975(1994).[12] M.E. Mura and P.K. Knowles,
J. Chem. Phys. , 9848(1996); V.I. Lebedev and D.N. Laikov,
Dokl. Math. ,477 (1999).[13] J. C. Slater, Phys. Rev. 81, 385 (1951); J. P. Perdew andY. Wang, Phys. Rev. B 45, 13244 (1992).[14] S.J. Chakravorty and E.R. Davidson, J. Phys. Chem . , 6167 (1996).[15] J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996).[16] O.A. Vydrov and G.E. Scuseria, J. Chem. Phys. ,8187 (2004).[17] A.D. Becke,
J. Chem. Phys. , 1372 (1993); , 5648(1993).[18] C. Adamo and V. Barone, J. Chem. Phys. 110, 6158(1999).[19] NIST, http://nist.gov/pml/data/. [20] D.A. Shirley et al. , Phys. Rev. B15