Communication: Self-interaction correction with unitary invariance in density functional theory Free

Density Functional Theory (DFT)^1,2 allows for a quantum-mechanical description of electrons without direct calculation of a many-electron wavefunction. Approximations to DFT, including the local-spin-density (LSD) approximation² and gradient approximations,^3,4 are widely and successfully used to predict, understand, and optimize many physical and chemical phenomena. However, such approximations have yet to be cast in a form that is both efficient and free of the self-interaction error (SIE).⁵ For example, in neutral systems the long-range behavior of the Kohn-Sham (KS) potential does not reduce to the −1/r form expected from general considerations. As discussed in Refs. 6,7, this and other SIEs lead to a range of related issues, sometimes referred to as delocalization errors,⁸ when using DFT approximations to understand chemistry, materials, and physics. Prior to introducing a size-extensive unitarily invariant self-interaction corrected (SIC) density functional theory, we review some physical and chemical systems where this class of functionals could offer improvements.

While standard DFT expressions allow for accurate predictions of electron affinities, based upon total energy differences, the SIE always pushes the occupied-orbital energies upward in a non-systematic way. This problem prevents the calculation of dipole-bound anions, such as ethylene carbonate,^9,10 and causes errors when DFT approximations are used for predicting the location and number of defect-levels in solids,^11,12 Rydberg excitations in atoms and molecules, or chemical processes mediated by differences in atomic- or fragment-electronegativity. With respect to the latter, the SIE is responsible for dissociation of molecules and crystals into fractionally charged fragments,¹³ and complicates the understanding of transport in systems containing gold leads and molecular islands.^14,15 While constrained DFT methods^12,16,17 have been used to address problems related to charge transfer, such constraints are also tacitly assumed when atomization energies are calculated. Although past applications of SIC may not provide empirical evidence that atomization energies can be improved in molecules,^18,19 the fact that published DFT atomization energies are generally constrained to non-fractional dissociation blunts this criticism of SIC. Calculated DFT atomic energies as a function of fractional occupation show, for at least localized systems, that the energy as a function of charge state is not linear for fractional electron numbers as required by the Perdew-Levy theorem that requires piece-wise linear energies and derivative discontinuities at integer electron states.⁶ Time-dependent DFT (TDDFT)^20,21 could be improved by inclusion of SIC. However, any new SIC method must address both DFT-based issues and issues related to the explicitly complex canonical orbitals that may result from TDDFT methods. Other manifestations of the self-interaction-induced delocalization error impact predictions of magnetic properties due to qualitatively incorrect predictions of valences in some 3d- and many 4f-electron systems,^22–25 and/or the delocalization of 3d orbitals, which causes overestimates in spin-excitations.²⁶ The broken-symmetry calculations that are often used for such magnetic calculations are also needed for describing bond-breaking of condensed phases into separated atoms. This raises questions about the so-called symmetry dilemma, due to the fact that a single determinant constructed from the KS orbitals is not always an eigenstate of symmetry or spin. There are good physical reasons²⁷ to permit such symmetry breaking in density functional approaches including ours here.

In the original formulation of the problem, for a given approximation to the exchange-correlation functional,

$Excapprox[n↑,n↓]$⁠, Perdew and Zunger⁵ suggested appending the following term to the DFT functional:

In the above equation, the SIC (localized) orbitals {ϕ_iσ} are used to define orbital densities according to:

$niσ(r)=|ϕiσ(r)|2$⁠. The terms U[n_{i, σ}] and

$E_{xc}^{approx} [n_{i,\sigma },0]$

$Excapprox[ni,σ,0]$ are the exact self-Coulomb and approximate self-exchange-correlation energies, respectively. The Perdew-Zunger (PZ) paper recognized that this formulation led to a definition for the energy functional that did not transform like the density, and posited that localized orbitals similar to those proposed by Edmiston and Ruedenberg,²⁸ rather than molecular KS orbitals, might be the most appropriate set of orbitals for defining the SIC.⁵ Shortly thereafter, Lin's Wisconsin SIC group followed up on this suggestion and introduced the concept of localized and canonical orbitals in self-interaction corrected theories.^29–32 These papers showed that the orbitals used for constructing the SIC energy must satisfy O(N²) localization equations given by

with

$ViσSIC$ the partial functional derivative of Eq. (1) with respect to the orbital density n_iσ. The Jacobi-like approach³¹ to solving these equations also depended on O(N²) Jacobi updates. The localized orbitals obtained from these equations were found to be topologically similar to sp³ hybrids in atoms, alternative energy-localized orbitals in molecules,²⁸ and Wannier functions in solids.²⁹ In addition, they satisfied an explicitly local Schrödinger-like equation that was coupled together by off-diagonal Lagrange multipliers.

The one-particle density matrices,

$ρσ(r,r′)=Σαψασ*(r)ψασ(r′)$ and densities n_σ(r) = ρ_σ(r, r), arising from single determinants or products of single determinants of each spin, are invariant under unitary transformations of the occupied KS orbitals {ψ_ασ} of each spin, and the value of such a many-electron wavefunction can only change by a phase factor under such transformations. The resulting energy expectation value from such wavefunctions is also invariant to unitary transformations within the space of orbitals used for constructing wavefunctions, transforming like the one-particle density matrices. While DFT also exhibits these symmetries, in SIC-DFT there are many possible choices for a unitary transformation from KS to SIC orbitals (real or complex,³³ satisfying energy-minimizing or other conditions). The standard ways to find SIC orbitals are computationally demanding, especially for systems of many atoms, and results have raised questions about size-extensivity⁶ in SIC-DFT.³⁴ Size-extensivity means that the total energy of separated density fragments equals the sum of the energies of the fragments. While the uncorrected functionals are size-extensive, the SIC functionals fail to be so if the SIC orbitals delocalize over distant atomic sites. These problems are solved in this Communication.

Here, a modification of the original formulation of the Perdew-Zunger self-interaction correction is introduced. For notational simplicity, we consider the large number of molecular and crystalline systems that can, in principle, be described by real KS orbitals and that have an integer number of spin up and spin down electrons. This formulation leads to a size-consistent SIC spin-density functional that is invariant to unitary transformations in the occupied orbital space and leads to a long-range effective potential that scales as −1/r. Rather than allowing the SIC localized orbitals to be any unitary transformation within the occupied orbital space, we introduce a constraint that the orbitals used for constructing Eq. (1) must be explicitly dependent on a quantity that is itself unitarily invariant. The Fermi orbital (FO)^35,36 is a specific example of such a quantity. Given a trial set of KS orbitals, the FO (F_iσ) is defined at any point in space, a_iσ, according to

In other words, the FO is simply the ratio of the one-particle spin-density matrix to the square root of the spin density and is ultimately a simple transformation of the KS orbitals. It is easy to verify that Wannier functions are a sub-class of Fermi Orbitals, and a few more comments illustrate their physical and chemical nature. At r = a_iσ, the value of the absolute square of the FO is identically equal to the total spin density at r = a_iσ. Further, the FO associated with any position, a_iσ, in space is normalized to unity. For special sets of points, the FOs can be immediately orthogonal (e.g., Wannier Functions) but for general sets of points, they are not orthogonal. Second, the absolute square of the FO is minus the exchange-hole density at r around an electron at a_iσ. The exact exchange energy of a single Slater determinant is simply

$Ex=−12Σσ∫d3r∫d3a|Σαψασ*(r)ψασ(a)|2/|r−a|.$ In Ref. 37, it was noted that minus the exchange-hole density had characteristics similar to a single-orbital density, but a simple-to-use local expression for replacing the non-local Fock operator remained difficult to obtain.

Because the FO accounts for all the spin density at a given point in space it is expected to be a rather localized function. So, to incorporate self-interaction corrections for the states with spin σ, in a system with N_↑ + N_↓ electrons, the PZ-SIC formulation can be slightly constrained by invoking the following strategy for each spin:

1. For a trial set of KS orbitals {ψ_ασ} find N_σ centroids

   $\lbrace {\bf a}_{1\sigma }, {\bf a}_{2\sigma },...,{\bf a}_{N_\sigma \sigma }\rbrace$
   ${a1σ,a2σ,...,aNσσ}$ which provide a set of N_σ normalized linearly independent, but not orthogonal FO
   $\lbrace F_{1\sigma },F_{2\sigma }...F_{N_\sigma \sigma }\rbrace$
   ${F1σ,F2σ...FNσσ}$ which, from Eq. (4), will always lie in the space spanned by the KS orbitals.

2. Use Löwdin's method of symmetric orthonormali-zation^31,38 to construct a set of localized orthonormal orbitals

   $\lbrace \phi _{1\sigma },\phi _{2\sigma },...,\phi _{N_\sigma \sigma }\rbrace$
   ${ϕ1σ,ϕ2σ,...,ϕNσσ}$ and construct the SIC-DFT energy of Eq. (1) from the set of FOs. This set is a unitary transformation on the KS orbitals.

3. Minimize the SIC-DFT energy as a function of the KS orbitals and the classical centroids of the FOs.

This approach bypasses solution of the localization equations in entirety at the lesser expense of the search for quasi-transferable and chemically appealing semi-classical FOCs (Fig. 1) on which the FO and unitary transformation depend.

Compared to the localization equations,³¹ the FO-SIC (FSIC) formalism adds a constraint that prevents consideration of all unitary transformations. An advantageous example is that the FO-formalism especially discourages the use of the KS orbitals such as Bloch functions, molecular orbitals, or orbital-angular momentum states in atoms for construction of the SIC energy (Eq. (1)). This point is key to regaining size extensivity. A simple molecular case that illustrates this is the He₂ molecule. To obtain σ_g and σ_u FO, it would be necessary to find a point where the value of the nodeless σ_g state is zero and the value of the σ_u state is non-zero, which is impossible.

The PW92 LSD functional was used to calculate all FSIC energies presented here. While future FSIC-DFT approaches will undoubtedly exploit sparsity that is provided by a direct self-consistent solution of the FSIC localized orbitals, in this work we have used the set of self-consistent KS orbitals to construct the spin-density matrices and then minimized the SIC-DFT energy of Eq. (1) by a rather tedious but straightforward direct brute-force variation of the FO-centroids {a_iσ}. To account for SIC-induced relaxation of the KS orbitals we used a strategy similar to the optimized effective potential method. Once a set of FO-centroids was obtained from a specific set of KS orbitals, we recalculated the PW92 FSIC-LSD energy using SCF orbitals obtained from different mixtures of all the functionals available in NRLMOL (the generalized-gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE), LSD(PW92), KS-Exchange Only, etc.). In some cases we found that the SCF GGA(PBE) KS orbitals provided slightly lower, albeit chemically insignificant, PW92 FSIC-LSD energies than the SCF LSD(PW92) KS orbitals. This two-step brute force process was iterated until the FSIC-LSD energy was minimized. The eigenvalues of the KS orbitals (Table II) were found by diagonalizing a Hamiltonian similar to that of Ref. 31 given by

with

$Vijiσ=⟨ϕiσ|ViσSIC|ϕjσ⟩$⁠. The calculations discussed here used a modified version of NRLMOL with large Gaussian-orbital basis sets.^39,40 It is actually very easy to find starting sets of Fermi-Orbital centroids that deliver a set of FOs that span the same space as the KS orbitals.

The FSIC formalism reduces the number of unknown parameters needed for describing the unitary transformation from O(N²) to 3N and reduces the number of equations from O(N²) to 3N. Thus, our analysis is that an optimal algorithm for solving the FSIC equations can easily scale as favorably as DFT and that, accounting for sparsity, an O(N) algorithm should be possible in the many-electron limit. In contrast, solution of the localization equations scales as poorly as O(N⁶) for a small number of electrons. After accounting for sparsity, solution of the localization equations should be slower by O(M³) with M a characteristic number of close neighbors that each atom has. This scaling argument shows that FSIC is potentially much faster than SIC and will actually be once an optimal algorithm is developed.

As a first test we have performed FSIC calculations on the six lightest closed shell-atoms (He, Be, Ne, Mg, Ar, Ca). In all cases, the FO procedure provides localized orbitals that resemble the sp³ hybrids and 1s-core orbitals identified in Ref. 32.

In Table I, we compare calculated atomization energies from LSD(PW92),³ GGA(PBE),⁴ SIC-LSD,³³ FSIC-LSD, and experiment. The FSIC-LSD provides significant improvements over LSD as compared to experiment and does as well or better than GGA(PBE) in some cases. The mean absolute error (MAE) for FSIC-LSD in Table I is significantly improved over SIC-LSD. This average improvement is almost entirely due to systematically better performance of FSIC in molecules with π bonds.

For example, the strongly covalent singlet N₂ molecule,

$1σg21σu22σg22σu21πu43σg2$⁠, which dissociates into an open-shell singlet with three unpaired 2p electrons per atom, is challenging to represent continuously as a broken-symmetry single determinant. Using FSIC-LSD, we find an atomization energy (D_e) of 9.80 eV which compares well to the experimental atomization energy (9.84 eV) and to high-accuracy complete active space self-consistent field (CASSCF) results (9.85 eV) of Li Manni et al.⁴¹ In Fig. 1, the valence FOCs are shown pictorially. As in Refs. 30,33, the FO-formalism creates doubly occupied 1s FOs on both atoms,

$1s_{A_+/A_-}=\lbrace 1\sigma ^{\prime }_g \pm 1\sigma ^{\prime }_u\rbrace /\sqrt{2}$

$1sA+/A−={1σg′±1σu′}/2$⁠, lone-pair states on the exterior of the molecule,

$2sp_{A+/A_-}=\lbrace 3\sigma ^{\prime }_g \pm 2\sigma ^{\prime }_u\rbrace /\sqrt{2}$

$2spA+/A−={3σg′±2σu′}/2$⁠, and three bond-centered banana orbitals (e.g.,

$\phi _n=\frac{1}{\sqrt{3}}[2\sigma _g^{\prime }\break-\sqrt{2}\lbrace \cos (\frac{2n\pi }{3})\pi _{ux} + \sin (\frac{2n\pi }{3}) \pi _{uy}\rbrace ]$

$ϕn=13[2σg′−2{cos(2nπ3)πux+sin(2nπ3)πuy}]$⁠, with n = −1, 0, +1 ). As in Ref. 30, the primes indicate that KS molecular orbitals of the same symmetry are mixed together by a unitary transformation within each irreducible representation to minimize Eq. (1). For example, the

$\lbrace 2\sigma _g^{\prime },3\sigma _g^{\prime }\rbrace$

${2σg′,3σg′}$ are determined by a nearly diagonal unitary mixture of the {2σ_g, 3σ_g} KS eigenstates. In Table II, we compare eigenvalues as calculated from HF, an earlier SIC-LSDX calculation,³¹ and experiment.

To summarize, we present and test a simplified and systematic theoretical framework for incorporating self-interaction corrections into a density-functional approximation. Our approach achieves computational efficiency and guarantees size-extensivity, as previous approaches may not. Compared to the accuracy of GGA(PBE), the results for FSIC-LSD are encouraging and greater accuracy is expected when FSIC is applied to meta-generalized gradient approximations.⁴² The localized orbitals obtained from the FO formulation are topologically similar to those obtained from the localization equations used earlier.³² However, the FSIC does not allow consideration of orbitals that correspond to saddle points and energy maxima whereas the original version of SIC did. Relative to the original formulation, the FO-induced constraint on the unitary transformations leading to localized orbitals allows FSIC to retain the symmetries and size-extensivity that are present in density-functional theory.

M.R.P. thanks Dr. R. A. Heaton who first suggested the possibility of using a FO in the SIC-LSD functional. J.P.P. acknowledges support from NSF Grant No. DMR-1305135.