This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

THE JOURNAL OF CHEMICAL PHYSICS 140, 18A513 (2014)

Self-interaction corrections in density functional theory Takao Tsuneda1,a) and Kimihiko Hirao2 1

Fuel Cell Nanomaterials Center, University of Yamanashi, Kofu 400-0021, Japan Computational Chemistry Unit, RIKEN Advanced Institute for Computational Science, Kobe, Hyogo 650-0047, Japan 2

(Received 17 December 2013; accepted 14 February 2014; published online 4 March 2014) Self-interaction corrections for Kohn-Sham density functional theory are reviewed for their physical meanings, formulations, and applications. The self-interaction corrections get rid of the selfinteraction error, which is the sum of the Coulomb and exchange self-interactions that remains because of the use of an approximate exchange functional. The most frequently used self-interaction correction is the Perdew-Zunger correction. However, this correction leads to instabilities in the electronic state calculations of molecules. To avoid these instabilities, several self-interaction corrections have been developed on the basis of the characteristic behaviors of self-interacting electrons, which have no two-electron interactions. These include the von Weizsäcker kinetic energy and long-range (far-from-nucleus) asymptotic correction. Applications of self-interaction corrections have shown that the self-interaction error has a serious effect on the states of core electrons, but it has a smaller than expected effect on valence electrons. This finding is supported by the fact that the distribution of self-interacting electrons indicates that they are near atomic nuclei rather than in chemical bonds. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4866996] I. SELF-INTERACTIONS OF ELECTRONS A. Historical background

The self-interaction correction (SIC) is one of the most frequently used corrections in density functional theory (DFT).1 DFT was developed as a way to calculate the electronic states of solids, and it uses pure functionals of electron density, ρ, and the gradient of density, ∇ρ. Consequently, DFT has been predominantly used in solid-state calculations. Despite this, various problems have been reported attributed to pure functionals, especially in the band calculations of solids. One of the main problems is underestimation of the semi-conductor band gaps, which are usually calculated as orbital energy gaps. This underestimation problem has been attributed to the self-interaction error (SIE) of pure exchange functionals. Based on the Pauli exclusion principle, SIE is the interactions of electrons with themselves that essentially should not exist. More precisely, the SIE indicates the sum of the Coulomb and exchange self-interactions, which remains due to the use of exchange functionals as a substitute for the Hartree-Fock exchange integral in the exchange part of the Kohn-Sham equation, E SIE =

n

(Jii + Exc [ρi ]) ,

(1)

i

where Jii is the Coulomb self-interaction of the ith orbital electron, Exc is the exchange-correlation energy functional, ρ i is the electron density of the ith orbital, and n is the number of electrons. The ratio of SIE in exchange functionals is so large

a) Electronic mail: [email protected]

0021-9606/2014/140(18)/18A513/13/$30.00

that it has been attributed to be the cause of the underestimation of the band gaps. The derivation of SIC can be traced back to the early days of quantum mechanics. In 1926, Hartree suggested a method for calculating many-electron systems, i.e., the Hartree method, in which the only electron-electron interactions are Coulomb interactions.2 This method implicitly excludes Coulomb self-interactions, Jii in Eq. (1), from the calculations. This SIC is then incorporated by cancelling the Coulomb and exchange self-interactions in the HartreeFock method, which was developed in 1930.3 Fermi and Amaldi developed the SIC for DFT. They suggested an SIC for Thomas-Fermi method,4, 5 which is the first DFT, as a correction for the Coulomb interactions.6 The wellknown form of this correction is the one for the Coulomb potential, n 1 ˆ Jj [ρ], (2) VJFA [ρ] = 2 1 − n j where Jˆj is the Coulomb operator for the interaction with the jth orbital electron. This correction assumes homogeneously distributed electrons like those in a uniform electron gas and it removes the SIE by eliminating the Coulomb potential of one electron. As an SIC for early DFT, the corrections based on this concept were applied to the Thomas-Fermi method,7 the local density approximation (LDA) exchange potential,8 and the Kohn-Sham method.9 Although this potential correction is based on a rather rough assumption, it is still used to determine an exchange-correlation potential in the constrained search method10 and in one type of the optimized effective potential (OEP) method.11 The most frequently used SIC is the Perdew-Zunger (PZ) correction,12 which simply removes the SIE from the total

140, 18A513-1

© 2014 AIP Publishing LLC

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

18A513-2

T. Tsuneda and K. Hirao

J. Chem. Phys. 140, 18A513 (2014)

energy density to identify localized electronic groups.20 The ELF is expressed as21

electronic energy, E PZ = E KS −

n

(Jii + Exc [ρi ]) .

(3)

ELF =

i

The SIE is also eliminated from the exchange potential, n δEx ˆ δEx Ji + − [ρi (r)] . VxPZ (r) = δρ δρ i

(4)

Similar SICs were suggested for the LDA exchange functional.13–16 This correction is simple and does not require many computations, and it has been used in many calculations including band calculations of solids. The PZ-SIC is explained in Sec. II A. B. Self-interacting electrons

The density matrix of self-interacting electrons, which have no two-electron interaction, can be easily derived in a way that cancels out the self-interactions17 stemming from Eq. (3). For σ -spin self-interacting electrons, the density matrix is represented as PσSI (r1 , r2 ) = ρσ1/2 (r1 )ρσ1/2 (r2 ),

(5)

where ρ σ is σ -spin electron density. Using this density matrix, the kinetic energy density, τ in the kinetic energy T ≡ d3 r τ , becomes the von Weizsäcker one, τ W , such as17 1 2 SI τ = τσ = − ∇ Pσ (r1 , r2 ) → τW 2 σ σ =

τσW =

σ

1 |∇ρσ |2 σ

8

ρσ

.

(6)

On the basis of this density matrix, the exchange energy density, xσ in the exchange energy Ex ≡ σ d3 r xσ , has an asymptotic behavior at long ranges (far from the nucleus)18 such that |P SI (r, r )|2 r→∞ 1 ρσ (r) . (7) −−−→ − xσ (r) = − d 3 r σ 2 |r − r | 2r Furthermore, the parallel-spin correlation energy density is assumed to be zero for self-interacting electrons, because the pair density matrix, which corresponds to one diagonal term of the parallel-spin second-order reduced density matrix, P2σ σ , for the self-interaction density matrix becomes zero to eliminate the parallel-spin correlation energy density, cσ σ in the parallel-spin correlation energy Ecσ σ ≡ d3 r cσ σ , as19

1 ρσ (r1 )ρσ (r2 )−|PσSI (r1 , r2 )|2 = 0 2 = ⇒ cσ σ = 0. (8)

Dσ /

σ

τσLDA

Dσ = τσ − τσW ,

2 ,

(9)

(10)

and τσLDA is the LDA kinetic energy density for σ -spin electrons, τσLDA =

3 (6π 2 )2/3 ρσ5/3 . 10

(11)

Since the von Weizsäcker kinetic energy density, τσW in Eq. (6), is positive and does not exceed the total kinetic energy density, τ σ , the Dσ in Eq. (10) is restricted to 0 ≤ Dσ ≤ 1 and therefore the ELF is also in the range of 0 ≤ ELF ≤ 1. This ELF is useful for elucidating chemical bonds and electron pairs,21–23 because the ELF is much smaller than 1 near chemical bonds in cases where the electron-electron interactions are far from the self-interactions. The long-range (far-from-nucleus) asymptotic behavior, which determines the exchange potentials, Vx , in the regions far from nuclei, has been incorporated in exchange functionals (see Sec. II C). For example, the Becke 1988 (B88) exchange functional24 gives the asymptotic behavior of the potential as 1 r→∞ Vx −−−→ − , r

(12)

for Slater-type wavefunctions. This asymptotic behavior is also used to explicitly correct the forms of exchange functionals in the regions far from atomic nuclei. Corrections of this sort are surveyed in Sec. II C. C. Self-interaction-free conditions

There are conditions for determining if the functionals are self-interaction-free. Zhang and Yang25 proposed the selfinteraction-free conditions for one-electron systems, which contain only the exchange self-interaction in two-electron interactions, for determining whether or not the kinetic, exchange, and correlation energies are self-interaction-free:

SI P2σ σ (r1 , r2 ) =

Similarly to the PZ-SIC in Eq. (3), these relations are also used to remove the SIE and to develop self-interaction-free functionals. The von Weizsäcker kinetic energy density in Eq. (6) can be used as an index to sort out self-interacting electrons. Section II B reviews the corrections based on the sorted selfinteracting electrons. Here, let us introduce the electron localization function (ELF), which uses the von Weizsäcker kinetic

σ

where Dσ is the difference between the total σ -spin kinetic energy density, τ σ , and von Weizsäcker one, τσW , in Eq. (6),

r1 =r2

1+

1

T [qρ1 ] = O(q),

(13)

Ex [qρ1 ] = q 2 Ex [ρ1 ],

(14)

Ec [qρ1 ] = 0,

(15)

and

where q is a fractional number in 0 ≤ q ≤ 1 and ρ 1 is the electron density of one-electron systems. These conditions have been used to show that the underestimation of the reaction barriers in conventional Kohn-Sham calculations is due to the SIE. Note that most exchange functionals have been developed as corrections to the LDA exchange functional. Since

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

18A513-3

T. Tsuneda and K. Hirao

the LDA exchange functional,26 3 3 1/3 ExLDA = − d 3 rρ 4/3 (r), 4 π

J. Chem. Phys. 140, 18A513 (2014)

(16)

4/3

is proportional to ρσ , the condition of Eq. (14) for this functional can be written as Ex [qρ1 ] = q 4/3 Ex [ρ1 ] ≤ q 2 Ex [ρ1 ].

(17)

This indicates that the exchange functional energies are overestimated around transition states, which are supposed to have a fractional q, and cause the calculated reaction barrier heights to be too low. In Ref. 25, this is confirmed for the dissociation of H+ 2 , which contains only the self-interaction for the electron-electron interaction. However, as will be explained in Sec. III C, it has been shown that the SIE is not the only cause of reaction barriers being underestimated in general chemical reactions. Self-interaction-free conditions have also been derived for multiple electrons. Mori-Sánchez et al.27 used the total energies for fractional occupations, the polarizabilities of polymers, and the dissociation potentials of molecules. The energy linearity theorem for fractional occupations28, 29 establishes that the total energies vary linearly as a fractional occupation number, i.e., p p q −p = E(n + 1) + E n+ E(n). (18) q q q

Eq. (18) and the nondivergent polarizabilities of long-chain polyenes are interpreted as the self-interaction-free conditions for the exchange-correlation integral kernel. In addition, the dissociation potentials of diatomic molecules can be used to check the extent that functionals are self-interactionfree, because the SIE is supposed to dominate them when the ionization potential and electron affinity of the atoms are close to each other. It was found that LC functionals give better potentials for the potential of CO− molecule than those of other functionals.27 In fact, LC functionals come close to satisfying the self-interaction-free conditions for multiple electrons, which are used as the self-interactionfree conditions for the exchange-correlation integral kernel. Note, however, that LC functionals are not exactly selfinteraction-free since their short-range parts usually contain SIE.

D. Self-interaction regions

(20)

Self-interaction regions, where electrons actually have no two-electron interaction, can be identified using the relation in Eq. (6).38 That is, those regions where the von Weizsäcker kinetic energy density is close to the total kinetic energy density, are self-interaction regions. Note that the von Weizsäcker kinetic energy density is the minimum of the total one, i.e., τσtotal ≥ τσW . Therefore, the ratio of the von Weizsäcker kinetic energy density to the total is a way to specify the selfinteraction regions in molecules. Figure 1 illustrates this ratio for the formaldehyde molecule.19 As shown in this figure, the self-interaction regions are localized only in the regions near the atomic nuclei and around the hydrogen atoms opposite the chemical bonds. That is, the self-interacting electrons are present only in the core orbitals and the nonbonding regions of valence 1s orbitals. Other regions are classified into those around bonds, where the ratio is far less than 1, and closelysituated regions, where the ratio is not so far from 1. In regions where the ratio is far less than 1, the electrons are assumed to be in motion like the free electrons in transition metals. Therefore, these regions are interpreted as being free-electron regions.38 The exchange-correlation potential in these regions can be well described by a density functional. Since the electron distributions vary slowly in these regions, the molecular

As mentioned in Sec. III A, the long-range corrected (LC) functional32 yields a curve close to being a straight line for the fractional occupation dependence of the total energies, while other functionals give concave curves.33 Pure and hybrid generalized-gradient-approximation (GGA) functionals have also been reported to overestimate the polarizabilities and hyperpolarizabilities of long-chain polymers. Moreover, although the PZ-SIC improves the polarizabilities and hyperpolarizabilities of H2 chains,34–36 it does not resolve the problem. Rather, for the polarizabilities and hyperpolarizabilities of long-chain polymers, a long-range correction dramatically improves the estimations.37 Note that the accuracy of LC functionals comes from their small SIE in the exchange-correlation integral kernel, δ 2 Exc /δρ 2 , rather than in the exchange-correlation energy or the exchangecorrelation potential.1, 33 Therefore, the energy linearity in

FIG. 1. Ratio of the von Weizsäcker kinetic energy density to the total one, τW /τtotal in Eq. (6), in formaldehyde molecule with three regions of electrons.

Combined with the Janak theorem,30 ∂E = i , ∂ni

(19)

where ni and i are the occupation number and orbital energy of the ith orbital electron, this theorem establishes that occupied and unoccupied orbital energies are, respectively, the corresponding negatives of the ionization potentials and electron affinities. Therefore, it is taken as the necessary condition for giving correct orbital energies that satisfy Eq. (18). This relation is also the basis of the orbital energy invariance theorem that the orbital energies of the outermost orbitals stay constant for the fractional occupation variation, i.e.,31 n+1 (n + 1) = n+1 (n).

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

18A513-4

T. Tsuneda and K. Hirao

J. Chem. Phys. 140, 18A513 (2014)

structures can even be described by pure GGA functionals. In these regions, the kinetic, exchange, and correlation energies have a physical connection.38 The remaining intermediate regions, where the ratio is not so far from 1, as long-range interaction regions, in which widely separated electrons interact with each other. Since these regions separate from the core orbitals (Fig. 1), this interpretation is consistent with the calculated results showing that the long-range correction hardly affects core orbitals (Sec. III A). It is also supported by analyses showing that in order to describe chemical bonds accurately, the long-range exchange and correlation interactions are required in addition to the self-interaction correction.39 Therefore, Fig. 1 shows that the SIC is needed in addition to the corrections for the long-range exchange and correlation interactions, to reproduce correct electronic states in the core and far-from-nucleus regions. II. SELF-INTERACTION CORRECTIONS A. Perdew-Zunger correction

As mentioned in Sec. I A, the PZ-SIC12 is the most frequently used SIC for removing the one-electron SIE from the exchange energy and potential functionals. This method simply eliminates the Coulomb and exchange self-interactions from the functionals: n (Jii + Exc [ρi ]) (21) E PZ = E KS − i

and VxPZ (r)

n δEx ˆ δEx Ji + − [ρi (r)] . = δρ δρ i

(22)

Note that this correction has a problem in that the Kohn-Sham equation is not invariant when making a unitary transformation of the occupied orbitals. By contrast, in the Hartree-Fock equation, the variations of the Coulomb self-interaction energy and its potential, when making a unitary transformation of occupied orbitals, cancel out with those of the exchange self-interaction; these self-interaction energies are not compensated in the Kohn-Sham equation even after the correction. Therefore, the effect of the SIC depends on the difference in occupied orbitals before and after the unitary transformation. To remove this difference, the initially proposed method is to localize the orbitals before performing the SIC.40 Note, however, that there are various orbital localization methods, and the effect of the SIC inevitably depends on them. Another presently accepted method is the explicit minimization of the SIC-DFT energy for the orbital variance.41, 42 Although this method provides unique solutions for the SIC-DFT, it is hard to be solved due to the non-unitary invariance. The PZ-SIC has been used in various chemical calculations. In particular, in their exploration of the effect of the correction, Vydrov and Scuseria43 showed that the SIE of valence orbitals does not depend on whether LDA, GGA, or meta-GGA functionals are used, and the correction adversely affects the reaction enthalpies. They concluded that the straightforward removal of the SIE in the PZ formalism is not a good strategy, because it upsets the balance of er-

ror cancellation.43 Many studies have reported instabilities in the PZ-SIC-DFT especially when it comes to calculating chemical reactions and properties of molecules.42, 44, 45 This means that the straightforward elimination of the SIE in the PZ-SIC is not appropriate in chemical calculations requiring well-balanced electron-electron interactions. Combining PZ-SIC with the OEP method46 may be an efficient way to overcome this instability. Tong and Chu47 proposed this combination method using the Krieger-Li-Iafrate (KLI) approximation.48 In this method, the total energy and exchange potential of the OEP-KLI equation are replaced with Eqs. (21) and (22). This method was found to improve the orbital energies (Sec. III A), ionization potentials, electron affinities,47 and excitation energies.49 Note, however, since the OEP method has the serious problems of poor selfinteraction correction (SCF) convergence and long computational time, orbital energy calculations have been restricted to atoms and small molecules. B. von Weizsäcker kinetic energy density-based corrections

Next, let us focus on the SIC using the von Weizsäcker kinetic energy density. Although this correction is not so popular, it has been implicitly incorporated in several exchange functionals. This is because the term attributed to this SIC naturally appears as the next term after the electron density term in the expansion of the exchange hole. Becke50 proposed that the exchange hole, hxσ in the exchange energy, ρσ (r1 )hxσ (r1 , r2 ) 1 Ex = − , (23) d 3 r1 d 3 r2 2 σ r12 can be expanded as 2 1 2 ∇ ρσ − 4Dσ r12 + · · · , (24) 6 where r12 = |r2 − r1 |. In Eq. (24), Dσ in Eq. (10) is the SIC term based on the von Weizsäcker kinetic energy density, and Eq. (6) means that it disappears for self-interacting electrons. Since the term with the density Laplacian, ∇ 2 ρ σ , integrates to zero, the exchange energy in Eq. (23) becomes the selfinteraction energy. On the basis of the exchange hole expansion, Becke and Roussel51 developed an exchange functional to recover the exact exchange energy for the hydrogen atom. This exchange functional naturally has the long-range (farfrom-nucleus) asymptotic behavior of Eq. (7).52 The von Weizsäcker kinetic energy density has also been explicitly used in several SICs. Tsuneda et al.19 proposed the regional SIC (RSIC), which corrects only for the exchange functionals in the self-interaction regions, where the ratio of the von Weizsäcker kinetic energy density to the total one, hxσ (r1 , r2 ) = ρσ +

tσ = τσW /τσtotal ,

(25)

53

is close to one. Jaramillo et al. suggested a similar SIC functional called the local hybrid (LH) functional using this tσ . The RSIC cuts out the self-interaction regions by using a region-separation function, fRS , and replaces the exchange DF DF DF and Vxσ = δxσ /δρσ , in the self-interaction functionals, xσ regions with the exchange energy density and potential of the

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

18A513-5

T. Tsuneda and K. Hirao

J. Chem. Phys. 140, 18A513 (2014)

SI SI exact exchange self-interactions, xσ and Vxσ , i.e., RSIC DF SI xσ = [1 − fRS (t)]xσ + fRS (t)xσ

(26)

RSIC DF SI Vxσ = [1 − fRS (t)]Vxσ + fRS (t)Vxσ .

(27)

and

19

The initial scheme used a spherically averaged σ -spin exchange energy density of the 1s orbital as the exact exchange energy density, ρσ 1s (r) = − [1 − (1 + αr)e−2αr ], (28) xσ 2r where α = |∇ρ σ |/(2ρ σ ); subsequently, a modified RSIC (mRSIC)54 was developed that used the exchange energy densities of hydrogen-like atomic orbitals, e.g., ρσ 1s (r) = − [1 − (1 + Zr)e−2Zr ], (29) xσ 2r

3 1 ρσ 1 2s 1 − 1 + Zr + Z 2 r 2 + Z 3 r 3 e−Zr , (r) = − xσ 2r 4 4 8 (30)

FIG. 2. Dependence of the region-separation function on the ratio of the von Weizsäcker kinetic energy density, tσ . The a in Eqs. (33) and (34) is 0.95.

the exact one, Pσexact , by using the density matrix similarity metric, σ , LH fRS = [1 − σ (r)]k tσl ,

1 3 exact SS σ (r) = 1 − d r12 Pσ (r, r + r12 )Pσ (r, r + r12 ) , ρ(r)

and

3 1 ρσ 1 2p xσ 1 − 1 + Zr + Z 2 r 2 + Z 3 r 3 e−Zr , (r) = − 2r 4 4 24 (31)

where Z is the charge of the nucleus to which the target electron belongs. However, these RSICs do not give accurate chemical properties. The pseudospectral RSIC (PSRSIC) uses the pseudospectral exchange energy density, χ ∗ (r )χμ (r ) 1 xPS = − Pμν Pλκ χν∗ (r)χλ (r) d 3 r κ , |r − r| 4 μνλκ

(37) where (k, l) = (1, 0) or (1.28, 2.78) for the LDA metric, , and (k, l) = (2, 0) or (0.88, 3.0) for the Perdew-BurkeLDA σ 58, 59 The spherically symmetErnzerhof (PBE) metric, PBE σ . SS ric density matrix, Pσ , is derived using the Iikura-TsunedaYanai-Hirao60 ansatz, i.e., 9π 1/2 1/3 kσ = ρσ , (38) Kσ 4/3 where Kσ is defined in Ex ≡ −(1/2) d 3 rρσ Kσ , as PσSS (r, r + r12 ) =

(32) where Pμν is the density matrix based on atomic orbitals and χ μ is the μth atomic orbital, and it can be used with the LC-DFT.55, 56 Note that the LH functional also uses the pseudospectral exchange energy density.53 The region-separation function, fRS , is the main difference between the PSRSIC and LH functionals, i.e.,

5(tσ − a) 1 RSIC 1 + erf , (33) fRS = 2 1−a 0 (tσ < a) PSRSIC fRS (34) = 1 (tσ ≥ a), and LH fRS = btσ ,

(35)

where a is the number of tσ s needed to cut out the selfinteraction regions, a = 0.95 in RSIC and a = −0.0204 Z + 1.0728, and b is the coefficient of tσ (b = 1 in LH(2003),53 b = 0.48 in LH(2007),57 and b = 0.706 in LH(2012)36 ). There are other types of LH functionals in which the exchange functionals are corrected for the regions where the corresponding spherically symmetric density matrices, PσSS , are dissimilar to

(36)

3j1 (kσ r12 ) ρσ (r). kσ r12

(39)

Figure 2 displays the region-separation functions in terms of tσ only for the functions of tσ . As shown in the figure, the RSIC and PSRSIC cut out only the self-interaction regions in contrast to the LH functionals mixing monotonically increasing pseudospectral exchange energy density. That is, the RSIC and PSRSIC functionals replace the exchange functional with the Hartree-Fock exchange energy only in the selfinteraction regions. On the other hand, the LH functionals mix the Hartree-Fock exchange integral in proportion to the similarity of the wavefunction to that of self-interacting electrons. C. Long-range (far-from-nucleus) asymptotic corrections

The long-range (far-from-nucleus) asymptotic behavior has also been used as an SIC for exchange functionals. As shown in Eq. (7), the asymptotic behavior of the exact exchange energy density, r→∞

xσ −−−→ −

ρσ , 2r

(40)

is derived from the density matrix of self-interacting electrons. To reproduce Eq. (40), the B88 exchange energy

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

18A513-6

T. Tsuneda and K. Hirao

J. Chem. Phys. 140, 18A513 (2014)

density24 is derived as 1 4/3 2ζ xσ2 3 1/3 B88 + xσ = − ρσ 3 , 2 4π 1 + 6ζ xσ sinh−1 xσ (41) 4/3 |∇ρσ |/ρσ

where xσ = and ζ = 0.0042. In Eq. (41), the second term in the parenthesis is used to reproduce Eq. (40), 2ζ xσ2 1 1 1/3 r→∞ (42) −−−→ − , − ρσ −1 2 2r 1 + 6ζ xσ sinh xσ for the Slater-type orbital, e.g., ρσ1s = π −1 exp(−2r), by using the relation,

r→∞ sinh−1 const.r b exp(ar) −−−→ ar, (43) where a and b are arbitrary real numbers. This relation has been used in the developments of many potentials such as the van Leeuwen-Baerends 1994 (LB94) potential mentioned later.61 Note that this feature is often not inherited by functionals derived from the B88 exchange functional such as the Perdew-Wang 1991 exchange functional.62, 63 There are SICs that explicitly use the long-range (farfrom-nucleus) asymptotic behavior to correct the exchange potentials. Surprisingly, the history of this asymptotic correction dates from 1955. Latter64 suggested an asymptotic correction using Eq. (12), i.e., 1 Latter LDA Vxσ (r) = max Vxσ (r), − , (44) r LDA is the LDA exchange potential.26 The LB94 powhere Vxσ tential mentioned above was developed by following the B88 exchange functional form61 xσ2 3 1/3 LB94 1/3 Vxσ = ρσ +β −2 , 4π 1 + 3βxσ sinh−1 xσ

(45) where β is a fitting parameter. However, these corrections have been found to give incorrect values for the electronic states of molecules. By using the asymptotic form of the exact exchange potential,65 nHOMO,σ r→∞ + IPHOMO,σ + HOMO,σ , (46) Vxσ (r) −−−→ − r where nHOMO is the occupation number of HOMO, Tozer and Handy66 and Casida et al.67 introduced a shift considering the r → ∞ limit of the exchange potential in the potential, 1 TH CCS Vxσ (r) = max Vxσ (r), − + IPHOMO,σ + HOMO,σ , r (47) where IPHOMO, σ and HOMO, σ are, respectively, the ionization potential and orbital energy corresponding to the σ -spin HOMO. Taking advantage of the LB94 and shifted potentials, Casida and Salahub68 proposed an asymptotically corrected LDA (AC-LDA) potential, LDA AC−LDA Vxσ (r) = max Vxσ (r) + IPHOMO,σ LB94 + HOMO,σ , Vxσ (r) . (48)

As mentioned in Sec. III B, the asymptotic corrections using the shift were found to improve the Rydberg excitation energies, which had been underestimated in the timedependent Kohn-Sham (TDKS) calculations. Note that the IPσ + HOMO, σ shift in Eq. (46) is introduced due to the fact that most Kohn-Sham calculations give poor HOMO energies. That is, this shift is an artificial term that is not required if the functional gives accurate orbital energies. As detailed in Sec. III A, LC functionals simultaneously provide accurate occupied and unoccupied orbital energies without requiring a self-interaction correction or artificial shift.33 Moreover, it has been shown that the TDKS method using LC functionals gives accurate excitation energies for Rydberg and charge transfer excitations.69 This indicates that the accurate Rydberg excitation energies of the potentials in Eqs. (47) and (48) are a result of the shift conducing the effect of the long-range correction.

D. Self-interaction errors in correlation functionals

SIE is often included in correlation functionals. The SIE in correlation functionals can be revealed by calculating the correlation energies of one-electron systems, which essentially must be zero as shown in Eq. (15). For example, LDA correlation functionals63, 70 generally contain parallel-spin electron correlation terms, which erroneously give correlation energies for one-electron systems. Since many GGA correlation functionals contain LDA correlation functionals,63, 71 these functionals also include SIE with no self-interaction correction. The Lee-Yang-Parr (LYP) correlation functional72 produces parallel-spin electron correlations, which give nonzero energies for one-electron systems in the form. However, a prefactor is multiplied with this functional in order for it to yield no electron correlation for one-electron systems. The one-parameter progressive (OP) correlation functional73 has no SIE, at least explicitly, because it contains only oppositespin electron correlations. Note, however, that the OP functional has an exchange functional term, so it is possible for it to have SIE for multiple electrons implicitly, although the errors are negligible in comparison with those in other correlation functionals. As mentioned above, despite the possibility of there being SIE in correlation functionals, the selfinteraction correction has rarely been used because the errors are assumed to be much smaller than those in exchange functionals. III. APPLICATIONS OF SELF-INTERACTION CORRECTED DFT A. Orbital energy calculations

SIE has been assumed to cause errors in band gap calculations of solids. Since the band gaps can be naturally viewed as the excitation energies of solids, they should be calculated with an excited state calculation method such as the TDKS method. However, in order to avoid the huge computational times required by the TDKS calculations, band gaps have usually been calculated under the assumption that they are close in value to the orbital energy gaps, despite there being no clear evidence that is so. Perdew74 and Görling and Levy75

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

18A513-7

T. Tsuneda and K. Hirao

J. Chem. Phys. 140, 18A513 (2014)

suggested that the discontinuity responsible for the inaccurate orbital energies, xc = {IP − EA} − {n+1 (n) − n (n)},

(49)

can be attributed to the error in the exchange part for the HOMO and LUMO, ∗ (r)φn+1 (r ) − φn∗ (r)φn (r )]Vxnl (r, r ) x = d 3 rd 3 r [φn+1 −

d 3 r[ρn+1 (r) − ρn (r)]Vx (r),

(50)

where Vxnl is the nonlocal exchange potential, Vxnl (r) =

n i

d 3 r φi∗ (r)

1 φi (r ). |r − r |

(51)

That is, the poor results for the band gaps have been explained as the difference in the nonlocal exchange potential and exchange potential functionals. The PZ-SIC tends to enlarge the exchange energies of occupied orbitals, and it has been reported to improve the accuracy of the band gaps estimated by the LDA functional.76 However, recent calculations using atom-centered basis sets for core orbitals have shown that the band gaps of semiconductors can be well reproduced even by pure functionals without the self-interaction correction77 and that this correction instead leads to overestimation of the band gap. It has also been found that the band gaps of insulators are overestimated by using the PZ-SIC.78 In fact, the lack of long-range exchange interactions is generally seen as a cause of underestimation of the band gaps of insulators.77 Next, let us consider the effect of the PZ-SIC on the valence orbital energies of molecules. There are several studies that have calculated orbital energies using the PZ-SIC KohnSham method.79–81 By calculating 44 HOMO energies and 32 LUMO energies, Vydrov and Scuseria79 showed that the

mean absolute errors of the calculated HOMO and LUMO energies are 1.67 and 1.61 eV for the PZ-SIC-PBE functional and 1.70 and 1.68 eV for the PZ-SIC-Tao-PerdewStaroverov-Scuseria meta-GGA functional. Vydrov, Scuseria, and Perdew82 also performed fractional occupation calculations in order to clarify the reproducibility of the orbital energies. Figure 3 shows the calculated total electronic energies and orbital energies of the PZ-SIC-PBE functional with respect to the fractional occupation number of the outermost orbitals of the carbon atom in order to check whether this functional obeys the energy linearity theorem in Eq. (18). In this figure, the functional using this self-interaction correction obeys neither the energy linearity theorem (Eq. (18)) nor the outermost orbital energy invariance theorem (Eq. (20)). Therefore, the PZ-SIC Kohn-Sham method does not give accurate valence orbital energies. Note, however, that this method affects the GGA exchange-correlation functionals, which violate these theorems, in the direction of improvement, although the overcorrection worsens the behavior. This result clearly suggests that the SIE is not the main cause for the poor quality of the valence orbital energies, although it may be partially responsible. Several studies using other SICs besides the PZ-SIC have calculated valence molecular orbital energies and valence band gaps. The LDA+U method83 is one of the most frequently used corrections taking SIC into account in band calculations. In this method, the difference of the Coulombexchange interactions from their averaged value is added to the total electronic energy of the LDA functional. Since this correction intends to incorporate spin-fluctuation effects, it is performed only for d-d and f-f orbital interactions. This method has been reported to give accurate band gaps.84 However, since this method is adjusted to band calculations, no orbital energy calculation of molecules has been reported for this method, as far as we know. The atomic SIC (ASIC)85 and variational pseudo-SIC (VPSIC)86 were developed for

FIG. 3. Calculated total electronic energies (left) and orbital energies (right) of the carbon atom with respect to the fractional occupation number. The PerdewZunger (PZ) self-interaction corrected (SIC) PBE functional leads to inaccuracies in the Kohn-Sham energies and orbital energies especially near integer occupation numbers. Reprinted with permission from O. A. Vydrov, G. E. Scuseria, and J. P. Perdew, J. Chem. Phys. 126, 154109 (2007). Copyright 2007 American Institute of Physics.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

18A513-8

T. Tsuneda and K. Hirao

J. Chem. Phys. 140, 18A513 (2014)

TABLE I. Mean absolute deviations (MADs) and mean absolute percent deviations (MA%Ds) of orbital energies from the corresponding ionization potentials (IPs) for HOMO and 1s orbitals or from the corresponding electron affinities (EAs) for LUMO in eV for each method.

HF+OP

BOP

B3LYP

LC-BOP

LC-PR -BOP

1.56 12.65 0.28

4.59 37.52 1.83

3.37 27.18 1.41

0.84 6.38 0.14

0.47 3.57 0.17

1.36 6.86 0.25

7.30 39.04 1.60

5.50 29.57 1.34

2.69 4.65 0.29

0.88 4.48 0.29

Core 1s orbitals of typical molecules (13 second-row and 5 third-row atoms) MAD 18.76 16.59 1s + IP1s (Second-row atoms) MA%D 4.13 3.79 MAD 30.29 29.99 1s + IP1s (Third-row atoms) MA%D 1.30 1.29

26.77 5.97 69.11 2.96

17.83 3.98 49.38 2.11

21.54 4.73 63.90 2.74

1.94 0.40 8.39 0.36

Deviation

HF

Valence orbitals of typical molecules (20 molecules) HOMO + IPHOMO MAD 1.80 MA%D 15.98 MAD 0.27 LUMO + EALUMO Valence orbitals of H and rare gas atoms (4 atoms) MAD 1.63 HOMO + IPHOMO MA%D 8.65 MAD 0.04 LUMO + EALUMO

the LDA exchange functional as ways to improve the calculations of molecular HOMO energies and band gaps, and they use an SIC potential with a nonlocal pseudopotentiallike projector,87 which is an approximation of the PZ-SIC potential. They also use atom-centered basis sets for the core orbitals. These SICs provide the accurate orbital energies for the HOMOs and even LUMOs. However, the LUMO energies are calculated as the HOMO energies of the anion radicals, and the studies in question did not compare the HOMO-LUMO gaps with IP+EA in these studies.85, 86 The AC-LDA potential in Eq. (48) has also been used to make the orbital energy calculations.88 In this correction, the calculated HOMO energies of the OEP-KLI method are used to evaluate the IPσ + HOMO, σ shift. A comparison of the occupied and unoccupied orbital energies found that the AC-LDA potential gives fairly close values to the orbital energies of the OEP-KLI potentials for many molecules.88 However, poor results were obtained for the unoccupied orbitals of several molecules, and the orbital energies of the OEP-KLI method themselves were reported to be not so accurate.1 Apparently, although the SIC positively affects the orbital energies, it is not a clear-cut solution for getting correct orbital energies. The lack of long-range exchange interactions has been suggested to be the main cause for the poor valence orbital energies.33 Table I summarizes the mean absolute (percent) errors in the calculated orbital energies of typical molecules from the corresponding negative ionization potentials and electron affinities.56 A comparison of the errors of pure BOP (B88 exchange + OP correlation) and LC-BOP functionals clearly shows that the long-range correction simultaneously improves the HOMO and LUMO energies of these typical molecules. Note that no Kohn-Sham method has been shown to provide accurate occupied and unoccupied orbitals simultaneously, even if highly sophisticated exchange-correlation potentials are used.89, 90 The accurate orbital energies of the LC-BOP functional have been verified by fractional occupation calculations. Figures 4 and 5 illustrate the total electronic energies and outermost orbital energies of the CO2 molecule

in terms of the fractional occupation variations. Note that the calculated total electronic energies are plotted from the energies on a straight line connecting the calculated total electronic energies of the cation and the neutral for the HOMO and of the neutral and the anion for the LUMO. According to the energy linearity theorem (Eq. (18)) and orbital energy invariance theorem (Eq. (20)), the total electronic energy should linearly vary and the outermost orbital energy should stay constant with respect to the fractional occupation number. The figures clearly show that the LC-BOP functional gives approximately constant HOMO and LUMO energies (though these HOMO energies slightly increase). Moreover, LC-BOP gives total electronic energies much closer to a straight line than those of other functionals. These results indicate that the lack of long-range exchange interactions is the main cause for the poor orbital energies in the Kohn-Sham calculations. However, Table I also shows that even the LC-BOP functional is unable to reproduce the HOMO energies of hydrogen and rare gas atoms and the core 1s orbital energies of typical molecules. What the SIC significantly improves is the core orbital energies. Tu et al.91 showed that the PZ-SIC improves the core orbital energies, but they did not give any results for valence orbitals. Figure 1 shows that self-interacting electrons are mainly found in core regions. Since the core regions separate from the long-range interaction regions, the above finding is supported by the fact the long-range correction hardly affects the calculated core orbital energies.54, 56 Recently, Nakata et al.55 combined the PSRSIC of Eqs. (32) and (34) with the LC functional to get what they called the LC-PR functional, and used it to perform the orbital energy calculations. Table I lists the errors in the calculated orbital energies of the LC-PR functional. Not only the core orbital energies of typical molecules but the HOMO energies of hydrogen and rare gas atoms are improved. Note that the PSRSIC also improves the accuracy of the HOMO energies for typical molecules, while maintaining the accuracy of the LUMO energies. This indicates that this LC-PR functional comprehensively

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

18A513-9

T. Tsuneda and K. Hirao

J. Chem. Phys. 140, 18A513 (2014)

FIG. 4. Calculated total electronic energies (a) and outermost orbital energies (b) of the carbon dioxide molecule in terms of the occupation number variation, n, of the HOMOs from the cation (n = −1) to the neutral (n = 0). The calculated total electronic energies are plotted from the energies on a straight line connecting the calculated total electronic energies of the cation and the neutral.

reproduces the energies of core and valence occupied and valence unoccupied orbitals. Figures 4 and 5 also show that this method maintains or slightly improves the physical validity of LC functionals in valence orbital energy calculations. However, this method has room for improvement when it comes to the core orbital energies. Figure 6 illustrates the the total electronic energies and outermost orbital energies of CO2 molecule in terms of the fractional occupation variation of the core 1s orbital. As the figure shows, the total electronic energy of the LC-PR functional changes linearly for small occupation number variation, but yields a concave curve for large occupation variation. Since the PSRSIC corrects only for the short-range part of the LC functional, it is supposed that the self-interaction regions of the short-range exchange interactions narrow as the occupation number of the 1s orbitals de-

creases. This suggests that there is something lacking as far as obtaining the correct core orbital energies goes. B. Excited state calculations

The most frequently performed excited state calculations that use the SIC are band calculations. However, as described in Sec. III A, the band gaps are usually evaluated as orbital energy gaps in the Kohn-Sham calculations. This section focuses on the SIC in the TDKS method for excited state calculations. Long-range (far-from-nucleus) asymptotic corrections in Sec. II C have been used to improve TDKS calculations of Rydberg excitation energies. In particular, it has been reported that an asymptotic correction with the r → ∞ limit of the

FIG. 5. Calculated total electronic energies (a) and outermost orbital energies (b) of the carbon dioxide molecule in terms of the occupation number variation, n, of the LUMOs from the neutral (n = 0) to the anion (n = 1). The calculated total electronic energies are plotted from the energies on a straight line connecting the calculated total electronic energies of the neutral and the anion.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

18A513-10

T. Tsuneda and K. Hirao

J. Chem. Phys. 140, 18A513 (2014)

FIG. 6. Calculated total electronic energies (a) and outermost orbital energies (b) of the carbon dioxide molecule in terms of the occupation number variation, n, of the core 1s orbitals from the cation (n = −1) to the neutral (n = 0). The calculated total electronic energies are plotted from the energies on a straight line connecting the calculated total electronic energies of the cation and the neutral.

exchange potential in Eqs. (47) and (48) enhances the accuracy of Rydberg excitation energies and maintains the accuracy of the valence excitation energies.66–68 However, as mentioned in Sec. II C, the asymptotic correction requires the IPσ + HOMO, σ shift in Eq. (46). Since the LC functionals provide accurate HOMO energies, this shift conduces the effect of the long-range correction. The long-range correction, actually, improves charge transfer energies besides the Rydberg excitation energies in the TDKS calculations.69 Moreover, incorporation of the PSRSIC in Eq. (32), which gives the correct −1/r behavior, hardly affects the Rydberg excitation energies of the LC-TDKS method.55 This seems to indicate that the asymptotic correction is not so important for getting accurate Rydberg excitations. However, since the LC functional uses the B88 exchange functional giving the correct −1/r behavior for the Slater-type wavefunction and the long-range correction might cover the asymptotic behavior in the far-fromnucleus regions, it is still an open question whether the −1/r behavior is significant to Rydberg excitations. There are several studies that have applied the PZSIC to the real-time TDKS method of calculating excitation energies.49, 92 By applying the PZ-SIC to the realtime TDKS method using the OEP-KLI-LDA potential, Hofmann et al.92 calculated the excitation spectra of SiH4 and 4-(N,N-dimethylamino)benzonitrile (DMABN) molecules. They found that the OEP-KLI-SIC improves the excitation energies of SiH4 , while slightly raising the main peak energies of DMABN. They also applied this method to long-range charge transfer (CT) excitation energy calculations of a dipeptide molecule and found that it improved these CT energies as well. However, they used the orbital product of the HOMO and LUMO of this molecule to represent the CT excitations and did not clarify the effect of the SIC on the long-range CT excitations. That the long-range correction overcomes the underestimation of problem for the long-range CT energies was proven in Refs. 69 and 93. Orbital energy calculations have shown that the SIC mainly affects core excitation energies and the long-range

correction hardly affects the core excitation energies. Nakata et al.54 calculated the core excitation energies by using the TDKS method with the mRSIC (Eqs. (29)–(31)). It was found that the mRSIC improves the core excitation energies, which are significantly underestimated even by LC functionals. These authors also found that the PSRSIC dramatically improved the TDKS core excitation energy while maintaining the accuracy of the KS core ionization potential calculations.55 The accuracy of the core orbital energy calculations56 comes from the SIC KS method giving the correct orbital energies. This indicates that the SIC must be incorporated in the TDKS calculations in order for them to yield correct core excitation energies.

C. Chemical reaction calculations

As mentioned in Sec. I C, the underestimated chemical reaction barriers have been attributed to the SIE. A few studies have calculated the reaction barriers with SIC-DFT. For example, it was suggested that the PZ-SIC slightly improves the accuracies of the reaction barriers calculated by the LDA and GGA functionals.94 The RSIC was also found to improve the reaction barriers of many reactions.19 Note, however, that SICs have been reported to worsen a significant fraction of reaction barriers.19, 94 Despite that, there are a few SICs that do give accurate reaction barriers, e.g., the LC-LDA exchange + short-range LDA correlation functional.36 These results suggest that the SIE is not the only cause of the underestimated barrier heights. Gräfenstein et al.39 investigated the effect of the SIE in bond dissociations on the basis of the dissociation behavior of radical cations with three-electron bonds. In their study, the SIE in bond dissociations are estimated in terms of the bond distance R as

E

SIE

= (1 − αHF )

1 1 SI −C J + , 2 4R

(52)

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

18A513-11

T. Tsuneda and K. Hirao

J. Chem. Phys. 140, 18A513 (2014)

where α HF is the mixing rate of the Hartree-Fock exchange integral in the exchange functional, C ≈ 2−1/3 , and JSI is the Coulomb self-interaction. A comparison with the results of highly sophisticated coupled-cluster singles, doubles and perturbative triples (CCSD(T)) calculations revealed that not only the SIE but also the approximate description of the interelectronic exchange contributes to the incorrect description of the bond dissociations. It was also found that the correct description of the dissociating three-electron bonds requires long-range corrections, which are neglected in standard DFTs including SIC-DFTs. The authors concluded that accurate three-electron bonds could not be obtained by incorporating SIC and that what was needed was the exchange-correlation functionals containing the exact exchange and long-range correlation effects. The necessity of long-range exchange-correlation effects in chemical reaction calculations is also pointed to in the calculations of LC-DFT and LC-DFT in combination with van der Waals correlation interactions,95–100 which is called the LC+vdW method.1, 32, 101 In the calculations of Database/3 (72 reactions) and other benchmark reactions (78 reactions),102–104 the long-range correction is found to halve the mean absolute errors of the calculated reaction barrier heights and enthalpies simultaneously.96 A remarkable improvement is also made in the calculations of condensation reactions.97 LC+vdW calculations using the local response dispersion functional105 indicate that the longrange correction significantly improves the isodesmic reaction enthalpies,98 isomerization reaction enthalpies,99 and DielsAlder reaction barriers and enthalpies.100 Note, however, that even the LC+vdW method may be unable to give accurate reaction enthalpies for the ring-openings of small molecules.99 Since the core orbitals obviously affect the reactions in these ring-opening reactions, the errors are supposed to be due to the lack of SIC. D. Chemical property calculations

Finally, let us discuss a number of interesting chemical property calculations that have used the SIC. For the optical properties, the SIC has been applied mainly to the polarizability and hyperpolarizability calculations of H2 chains. As mentioned in Sec. I C, the PZ-SIC has been reported to improve the overestimated polarizabilities and hyperpolarizabilities.34–36 However, the applications are restricted to the H2 chains, in which the self-interactions are dominant in the electron-electron interactions. The PZSIC is also applied to the magnetic property calculations. Patchkovskii et al.106 suggested that the PZ-SIC for the LDA functional in the OEP-KLI method significantly improves the isotropic nuclear magnetic resonance chemical shifts for N, O, F, while it worsens those for C and H. They also found that the PZ-SIC improves the spin-spin couplings for small Fermi contact terms, although it worsens them for large Fermi contact terms. Ruiz et al.107 discussed the effect of the PZSIC on the exchange coupling constants, J, of the Heisenberg Hamiltonian, Hˆ = −J Sˆ1 Sˆ2 ,

(53)

where Sˆn is the operator giving the spin of the nth paramagnetic center. As a result, they found that the inclusion of the PZ-SIC cancels the nondynamical electron correlation contributions, which are implicitly involved in exchange functionals, so as to improve the accuracy of exchange coupling constants. These studies support the conclusion that the PZ-SIC is not a clear-cut solution for improving the poor results of Kohn-Sham calculations. IV. CONCLUSIONS

The self-interaction correction removes the selfinteraction error which is the sum of the Coulomb and exchange self-interactions that should be cancelled out but remains as a result of the use of an approximate exchange functional. For self-interacting electrons, the electronic motion states are made known: self-interacting electrons have the von Weizsäcker kinetic energy, move in the −1/r exchange potential far from nuclei, and undergo no electron correlation with electrons of the same spin. Exchange functionals usually contain the SIE in the exchange energies, Ex , and potentials, Vx = δEx /δρ. However, these functionals also contain the SIE in the exchange integral kernels, fx = δ 2 Ex /δρ 2 , which is often taken as the SIE of multiple electrons. Note that although the SIE in fx affects the reproducibility of orbital energies, usual SICs are unable to reduce it in valence orbitals. The regions where self-interacting electrons are present, which are called self-interaction regions, are near the atomic nuclei, around the hydrogen atoms opposite the chemical bonds, and far from nuclei. So far, the most frequently used SIC is the PerdewZunger SIC that directly removes the SIE. However, this SIC has an unstable nature in the electronic state calculations of molecules. Several types of SICs have been developed in an effort to get rid of this instability: The combination of the PZ-SIC with the OEP-KLI method is a good strategy, but it requires substantial computational time. There is a type based on the von Weizsäcker kinetic energy density, e.g., an exchange functional which implicitly has a form that gives zero energy for self-interacting electrons and SICs that modify exchange functionals when the kinetic energy is nearly equal to the von Weizsäcker one. The other type includes a long-range (far-from-nucleus) asymptotic correction using the −1/r potential. SICs have been applied to various calculations, especially those of orbital energies, excitation energies, and chemical reactions. In particular, the PZ-SIC has long been used in the calculations of band energies, which are essentially excited state energies. However, it has recently been reported in solid state physics that the PZ-SIC overestimates the valence band gaps. SICs show promise in core orbital energy calculations, and a long-range (far-from-nucleus) asymptotic correction has been used to improve the accuracy of Rydberg excitation energies. SICs are of the greatest benefit in core excitation energy calculations. On the other hand, in chemical reaction calculations, SIE is considered to be not the only cause of the underestimated barrier heights. Nevertheless, detailed analyses have shown that the lack of long-range exchangecorrelation interactions are more significant than the SIE.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

18A513-12

T. Tsuneda and K. Hirao

Combining SIC with a correction for long-range exchange-correlation interactions will make it possible to perform precise DFT calculations on a wide variety of electronic states. ACKNOWLEDGMENTS

This research was supported by the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) (Grant Nos. 23225001 and 24350005). 1 T.

Tsuneda, Density Functional Theory in Quantum Chemistry (Springer, Tokyo, 2014). 2 D. R. Hartree, Math. Proc. Cambridge Philos. Soc. 24, 89 (1928). 3 V. Fock, Z. Phys. 61, 126 (1930). 4 L. H. Thomas, Math. Proc. Cambridge Philos. Soc. 23, 542 (1927). 5 E. Fermi, Z. Phys. 48, 73 (1928). 6 E. Fermi and E. Amaldi, Accad. Ital. Rome 6, 119 (1934). 7 C. A. Coulson and C. S. Sharma, Proc. Phys. Soc. 79, 920 (1962). 8 R. D. Cowan, Phys. Rev. 163, 54 (1967). 9 G. W. Bryant and G. D. Mahan, Phys. Rev. B 17, 1744 (1978). 10 Q. Zhao, R. C. Morrison, and R. G. Parr, Phys. Rev. A 50, 2138 (1994). 11 N. I. Gidopoulos and N. N. Lathiotakis, J. Chem. Phys. 136, 224109 (2012). 12 J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). 13 I. Lindgren, Int. J. Quantum Chem. 5, 411 (1971). 14 M. S. Gopinathan, Phys. Rev. A 15, 2135 (1977). 15 H. Stoll, C. M. E. Pavlidou, and H. Preuss, Theor. Chim. Acta 149, 143 (1978). 16 J. P. Perdew, Chem. Phys. Lett. 64, 127 (1979). 17 R. M. Dreizler and E. K. U. Gross, Density-Functional Theory. An Approach to the Quantum Many-Body Problem (Springer, Berlin, 1990). 18 M. Levy, J. P. Perdew, and V. Sahni, Phys. Rev. A 30, 2745 (1984). 19 T. Tsuneda, M. Kamiya, and K. Hirao, J. Comput. Chem. 24, 1592 (2003). 20 A. D. Becke and K. E. Edgecombe, J. Chem. Phys. 92, 5397 (1990). 21 A. Savin, A. D. Becke, J. Flad, R. Nesper, H. Preuss, and H. G. von Schnering, Angew. Chem., Int. Ed. Engl. 30, 409 (1991). 22 B. Silvi and A. Savin, Nature 371, 683 (1994). 23 A. Savin, R. Nesper, S. Wengert, and T. F. Fässler, Angew. Chem., Int. Ed. Engl. 36, 1808 (1997). 24 A. D. Becke, Phys. Rev. A 38, 3098 (1988). 25 Y. Zhang and W. Yang, Phys. Rev. Lett. 80, 890 (1998). 26 P. A. M. Dirac, Math. Proc. Cambridge Philos. Soc. 26, 376 (1930). 27 P. Mori-Sánchez, A. Cohen, and W. Yang, J. Chem. Phys. 125, 201102 (2006). 28 J. P. Perdew, R. G. Parr, M. Levy, and J. L. J. Balduz, Phys. Rev. Lett. 49, 1691 (1982). 29 W. Yang, Y. Zhang, and P. W. Ayers, Phys. Rev. Lett. 84, 5172 (2000). 30 J. F. Janak, Phys. Rev. B 18, 7165 (1978). 31 L. J. Sham and M. Schlüter, Phys. Rev. B 32, 3883 (1985). 32 T. Tsuneda and K. Hirao, “Long-range correction for density functional theory,” WIREs Comput. Mol. Sci. doi:10.1002/wcms.1178 (2013). 33 T. Tsuneda, J.-W. Song, S. Suzuki, and K. Hirao, J. Chem. Phys. 133, 174101 (2010). 34 A. Ruzsinszky, J. P. Perdew, G. I. Czonka, G. E. Scuseria, and O. A. Vydrov, Phys. Rev. A 77, 060502 (2008). 35 C. D. Pemmaraju, S. Sanvito, and K. Burke, Phys. Rev. B 77, 121204 (2008). 36 A. V. Arbuznikov and M. Kaupp, J. Chem. Phys. 136, 014111 (2012). 37 M. Kamiya, H. Sekino, T. Tsuneda, and K. Hirao, J. Chem. Phys. 122, 234111 (2005). 38 T. Tsuneda, M. Kamiya, N. Morinaga, and K. Hirao, J. Chem. Phys. 114, 6505 (2001). 39 J. Gräfenstein, E. Kraka, and D. Cremer, Phys. Chem. Chem. Phys. 6, 1096 (2004). 40 B. G. Johnson, C. A. Gonzales, P. M. W. Gill, and J. A. Pople, Chem. Phys. Lett. 221, 100 (1994). 41 M. R. Pederson, R. A. Heaton, and C. C. Lin, J. Chem. Phys. 80, 1972 (1984). 42 S. Goedecker and C. J. Umrigar, Phys. Rev. A 55, 1765 (1997). 43 O. A. Vydrov and G. E. Scuseria, J. Chem. Phys. 121, 8187 (2004).

J. Chem. Phys. 140, 18A513 (2014) 44 V.

Polo, E. Kraka, and D. Cremer, Mol. Phys. 100, 1771 (2002). Vieira and K. Capelle, J. Chem. Theory Comput. 6, 3319 (2010). 46 J. D. Talman and W. F. Shadwick, Phys. Rev. A 14, 36 (1976). 47 X.-M. Tong and S.-I. Chu, Phys. Rev. A 55, 3406 (1997). 48 J. B. Krieger, Y. Li, and G. J. Iafrate, Phys. Rev. A 45, 101 (1992). 49 D. Hofmann and S. Kümmel, J. Chem. Phys. 137, 064117 (2012). 50 A. D. Becke, Int. J. Quantum Chem. 23, 1915 (1983). 51 A. D. Becke and M. R. Roussel, Phys. Rev. A 39, 3761 (1989). 52 R. Neumann and N. C. Handy, Chem. Phys. Lett. 246, 381 (1995). 53 J. Jaramillo, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118, 1068 (2003). 54 A. Nakata, T. Tsuneda, and K. Hirao, J. Comput. Chem. 30, 2583 (2009). 55 A. Nakata, T. Tsuneda, and K. Hirao, J. Phys. Chem. A 114, 8521 (2010). 56 A. Nakata and T. Tsuneda, J. Chem. Phys. 139, 064102 (2013). 57 H. Bahmann, A. V. Rodenberg, A. V. Arbuznikov, and M. Kaupp, J. Chem. Phys. 126, 011103 (2007). 58 B. G. Janesko and G. E. Scuseria, J. Chem. Phys. 127, 164117 (2007). 59 B. G. Janesko and G. E. Scuseria, J. Chem. Phys. 128, 084111 (2008). 60 H. Iikura, T. Tsuneda, T. Yanai, and K. Hirao, J. Chem. Phys. 115, 3540 (2001). 61 R. van Leeuwen and E. J. Baerends, Phys. Rev. A 49, 2421 (1994). 62 J. P. Perdew, in Electronic Structure of Solids, edited by P. Ziesche and H. Eschrigh (Akademie Verlag, Berlin, 1991). 63 J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). 64 R. Latter, Phys. Rev. 99, 510 (1955). 65 C. O. Almbladh and U. von Barth, Phys. Rev. B 31, 3231 (1985). 66 D. J. Tozer and N. C. Handy, J. Chem. Phys. 109, 10180 (1998). 67 M. E. Casida, K. C. Casida, and D. R. Salahub, Int. J. Quantum Chem. 70, 933 (1998). 68 M. E. Casida and D. R. Salahub, J. Chem. Phys. 113, 8918 (2000). 69 Y. Tawada, T. Tsuneda, S. Yanagisawa, T. Yanai, and K. Hirao, J. Chem. Phys. 120, 8425 (2004). 70 S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980). 71 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 72 C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988). 73 T. Tsuneda, T. Suzumura, and K. Hirao, J. Chem. Phys. 110, 10664 (1999). 74 J. P. Perdew, Density Functional Methods in Physics, NATO Advanced Science Institutes Series Vol. 123, edited by R. M. Dreizler and J. d. Providencia (Plenum Press, New York, 1985). 75 A. Görling and M. Levy, Phys. Rev. A 52, 4493 (1995). 76 A. Svane and O. Gunnarsson, Phys. Rev. Lett. 65, 1148 (1990). 77 I. G. Gerber, J. G. Angyan, M. Marsman, and G. Kresse, J. Chem. Phys. 127, 054101 (2007). 78 M. Arai and T. Fujiwara, Phys. Rev. B 51, 1477 (1995). 79 O. A. Vydrov and G. E. Scuseria, J. Chem. Phys. 122, 184107 (2005). 80 S. Klüpfel, P. Klüpfel, and H. Jónsson, Phys. Rev. A 84, 050501 (2011). 81 S. Klüpfel, P. M. Dinh, P.-G. Reinhard, and E. Suraud, Phys. Rev. A 88, 052501 (2013). 82 O. A. Vydrov, G. E. Scuseria, and J. P. Perdew, J. Chem. Phys. 126, 154109 (2007). 83 A. I. Lichtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467(R) (1995). 84 V. I. Anisimov, F. Aryasetiawan, and A. I. Lichtenstein, J. Phys.: Condens. Matter 9, 767 (1997). 85 C. D. Pemmaraju, T. Archer, D. Sánchez-Portal, and S. Sanvito, Phys. Rev. B 75, 045101 (2007). 86 A. Filippetti, C. D. Pemmaraju, S. Sanvito, P. Delugas, D. Puggioni, and V. Fiorentini, Phys. Rev. B 84, 195127 (2011). 87 A. Filippetti and N. A. Spaldin, Phys. Rev. B 67, 125109 (2003). 88 S. Hamel, M. E. Casida, and D. R. Salahub, J. Chem. Phys. 116, 8276 (2002). 89 I. V. Schweigert and R. J. Bartlett, J. Chem. Phys. 129, 124109 (2008). 90 A. M. Teale, F. De Proft, and D. J. Tozer, J. Chem. Phys. 129, 044110 (2008). 91 G. Tu, V. Carravetta, O. Vahtras, and H. Ågren, J. Chem. Phys. 127, 174110 (2007). 92 D. Hofmann, T. Körzdörfer, and S. Kümmel, Phys. Rev. Lett. 108, 146401 (2012). 93 M. Chiba, T. Tsuneda, and K. Hirao, J. Chem. Phys. 124, 144106 (2006). 94 S. Patchkovskii and T. Ziegler, J. Chem. Phys. 116, 7806 (2002). 95 M. Kamiya, T. Tsuneda, and K. Hirao, J. Chem. Phys. 117, 6010 (2002). 96 J.-W. Song, S. Tokura, T. Sato, M. A. Watson, and K. Hirao, J. Chem. Phys. 127, 154109 (2007). 45 D.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05

18A513-13 97 R.

T. Tsuneda and K. Hirao

K. Singh, T. Tsuneda, and K. Hirao, Theor. Chem. Acc. 130, 153 (2011). 98 J.-W. Song, T. Tsuneda, T. Sato, and K. Hirao, Org. Lett. 12, 1440 (2010). 99 J.-W. Song, T. Tsuneda, T. Sato, and K. Hirao, Theor. Chem. Acc. 130, 851 (2011). 100 R. K. Singh and T. Tsuneda, J. Comput. Chem. 34, 379 (2013). 101 T. Tsuneda, and T. Taketsugu, in π -Stacked Polymers and Molecules: Theory, Synthesis, and Properties, edited by T. Nakano (Springer, Tokyo, 2013).

J. Chem. Phys. 140, 18A513 (2014) 102 B.

J. Lynch and D. G. Truhlar, J. Phys. Chem. A 107, 3898 (2003). Zhao, N. González-Garcia, and D. G. Truhlar, J. Phys. Chem. A 109, 2012 (2005). 104 J. Z. Pu and D. G. Truhlar, J. Phys. Chem. A 109, 773 (2005). 105 T. Sato and H. Nakai, J. Chem. Phys. 131, 224104 (2009). 106 S. Patchkovskii, J. Autschbach, and T. Ziegler, J. Chem. Phys. 115, 26 (2001). 107 E. Ruiz, S. Alvarez, J. Cano, and V. Polo, J. Chem. Phys. 123, 164110 (2005). 103 Y.

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 131.111.164.128 On: Fri, 14 Nov 2014 22:39:05