Restoration of the Derivative Discontinuity in Kohn-Sham Density Functional Theory:
An Efficient Scheme for Energy Gap Correction
Jeng-Da Chai1,2,*and Po-Ta Chen1
1Department of Physics, National Taiwan University, Taipei 10617, Taiwan
2Center for Theoretical Sciences and Center for Quantum Science and Engineering, National Taiwan University, Taipei 10617, Taiwan (Received 5 June 2012; published 15 January 2013)
From the perspective of perturbation theory, we propose a systematic procedure for the evaluation of the derivative discontinuity (DD) of the exchange-correlation energy functional in Kohn-Sham (KS) density functional theory, wherein the exact DD can in principle be obtained by summing up all the perturbation corrections to infinite order. Truncation of the perturbation series at low order yields an efficient scheme for obtaining the approximate DD. While the zeroth-order theory yields a vanishing DD, the first-order correction to the DD can be expressed as an explicit universal functional of the ground-state density and the KS lowest unoccupied molecular orbital density, allowing the direct evaluation of the DD in the standard KS method without extra computational cost. The fundamental gap can be predicted by adding the estimated DD to the KS gap. This scheme is shown to be accurate in the prediction of the fundamental gaps for a wide variety of atoms and molecules.
DOI:10.1103/PhysRevLett.110.033002 PACS numbers: 31.15.E!, 71.15.Mb, 71.10.!w
Over the past two decades, Kohn-Sham density func- tional theory (KS DFT) [1] has become one of the most powerful theoretical methods for studying the ground-state properties of electronic systems. As the exact exchange- correlation (XC) energy functional Exc½!# in KS DFT remains unknown, functionals based on the local density approximation (LDA) and generalized gradient approxi- mations have been widely used for large systems, due to their computational efficiency and reasonable accuracy.
However, owing to their qualitative failures in a number of situations [2–5], resolving these failures at a reasonable computational cost continues to be the subject of intense research interest.
The prediction of the fundamental gap Eg has been an important and challenging subject in KS DFT [6–21]. For a system of N electrons (N is an integer) in the presence of an external potential vextðrÞ, Egis defined as
Eg¼ IðNÞ ! AðNÞ; (1)
where IðNÞ ¼ EðN ! 1Þ ! EðNÞ is the vertical ionization potential and AðNÞ ¼ EðNÞ ! EðN þ 1Þ is the vertical electron affinity, with EðNÞ being the ground-state energy of the N-electron system. Therefore, Eg can be extracted from three KS calculations for the ground-state energies of the N- and (N( 1)-electron systems. However, such mul- tiple energy-difference calculations are inapplicable for the prediction of fundamental band gaps of solid-state systems [6–11,14,17,18].
By contrast, the KS gap !KS is defined as the energy difference between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) of the N-electron system [22–24],
!KS ¼ "Nþ1ðNÞ ! "NðNÞ; (2)
where "iðNÞ is the ith KS orbital energy of the N-electron system. Therefore, !KScan be obtained from only one KS calculation for the KS orbital energies of the N-electron system. Note that Egis not simply !KSbut is given by
Eg¼ !KSþ !xc; (3)
where
!xc¼ lim
#!0þ
!$Exc½!#
$!ðrÞ
"
"
"
"
"
"
"
"Nþ#!$Exc½!#
$!ðrÞ
"
"
"
"
"
"
"
"N!#
# (4) is the derivative discontinuity (DD) of Exc½!# [23–32]. As the KS gap (even with the exact functional) severely under- estimates the fundamental gap [24,25,31], the evaluation of the DD is tremendously important. Recently, the impor- tance of the DD in the excited-state [33] and time- dependent [34,35] properties has also been highlighted.
Although several schemes have been proposed for calcu- lating the DD, they can be very computationally demand- ing for large systems, due to the use of the nonlocal energy-dependent self-energy operators [6,7,10,11,18] or of the Hartree-Fock operator [8,13,21].
In this Letter, we provide a systematic procedure for the evaluation of the DD, based on perturbation theory [36].
The lowest-order estimate of the DD can be expressed as an explicit universal (i.e., system-independent) functional of the ground-state density and the KS LUMO density, allowing very efficient and accurate calculations of the DD and, via Eq. (3), the fundamental gap in the standard KS method.
For the exact KS DFT, IðNÞ ¼ !"NðNÞ [22,23,28,37–39], and therefore AðNÞ ¼ IðN þ 1Þ ¼ !"Nþ1ðN þ 1Þ. Con- sequently, Eg[see Eq. (1)] can be expressed as
Eg ¼ "Nþ1ðN þ 1Þ ! "NðNÞ; (5)
which is simply the energy difference between the HOMOs of the N- and (Nþ 1)-electron systems [26]. By subtracting Eq. (2) from (5), !xc [see Eq. (3)] can be expressed as
!xc¼ "Nþ1ðN þ 1Þ ! "Nþ1ðNÞ
¼ h ~cNþ1j ~HKSj ~cNþ1i ! hcNþ1j HKSjcNþ1i:
(6) Here, HKS)f!2m@2er2þvextðrÞþe2R!ðr0Þ
jr!r0jdr0þvxcð½!#;rÞg andciðrÞ are, respectively, the KS Hamiltonian and the ith KS orbital of the N-electron system, with vxcð½!#; rÞ being the XC potential and !ðrÞ ¼PN
i¼1jciðrÞj2 being the ground-state density. H~KS) f!2m@2er2þ vextðrÞ þ e2R !~ðr0Þ
jr!r0jdr0þ vxcð½~!#; rÞg and ~ciðrÞ are, respectively, the KS Hamiltonian and the ith KS orbital of the (Nþ 1)- electron system, with vxcð½~!#; rÞ being the XC potential and ~!ðrÞ ¼PNþ1
i¼1 j ~ciðrÞj2being the ground-state density.
Aiming to compute !xc (and hence Eg) using only one KS calculation for the N-electron system (e.g., for being applicable to solids), we express "Nþ1ðN þ 1Þ in terms of f"iðNÞ;ciðrÞg, based on perturbation theory [36].
We choose HKS as the unperturbed Hamiltonian and suppose that the unperturbed energy levels are nondegen- erate. Let % be a dimensionless parameter, ranging continuously from 0 (no perturbation) to 1 (the full pertur- bation). Consider the perturbed Hamiltonian H%given by
H%¼ HKSþ %H0%; (7) where the perturbation H0%)e2R!~%ðr0Þ
jr!r0jdr0þvxcð½~!%#;rÞ!
e2R!ðr0Þ
jr!r0jdr0!vxcð½!#;rÞ involves ~!%ðrÞ )PNþ1
i¼1 j ~c%iðrÞj2 (filling the orbitals in order of increasing energy). Here, f ~c%iðrÞg and f"%iðN þ 1Þg are, respectively, the eigenstates and eigenvalues of H%:
H%c~%iðrÞ ¼ "%iðN þ 1Þ ~c%iðrÞ: (8) Equation (8) at %¼ 1 is simply the KS equation for the (Nþ 1)-electron system, as it can be verified that f ~ciðrÞg and f"iðN þ 1Þg are, respectively, the eigenstates and eigenvalues of H%¼1. Therefore, "Nþ1ðN þ 1Þ ¼
"%N¼1þ1ðN þ 1Þ.
Writing H%0, ~c%iðrÞ, and "%iðN þ 1Þ as a power series in
%, we have
H0%¼ H0ð0Þþ %H0ð1Þþ %2H0ð2Þþ * * * ; (9)
c~%iðrÞ ¼cð0Þi ðrÞ þ %cð1Þi ðrÞ þ %2cð2Þi ðrÞ þ * * * ; (10)
"%iðN þ 1Þ ¼ "ð0Þi þ %"ð1Þi þ %2"ð2Þi þ * * * : (11) Inserting Eqs. (7) and (9)–(11), into Eq. (8) gives
ðHKSþ%H0ð0Þþ%2H0ð1Þþ***Þðcð0Þi þ%cð1Þi þ%2cð2Þi þ***Þ
¼ð"ð0Þi þ%"ð1Þi þ%2"ð2Þi þ***Þðcð0Þi þ%cð1Þi þ%2cð2Þi þ***Þ:
(12) Expanding Eq. (12) and comparing the coefficients of each power of % yield an infinite series of simultaneous equations.
To zeroth orderð%0Þ in Eq. (12), the equation is HKScð0Þi ðrÞ ¼ "ð0Þi cð0Þi ðrÞ; (13) which is simply the KS equation for the N-electron system (i.e., the unperturbed system). We then have cð0Þi ðrÞ¼ciðrÞ and "ð0Þi ¼"iðNÞ. Therefore, "Nþ1ðNþ1Þ¼
"%N¼1þ1ðNþ1Þ+"ð0ÞNþ1¼"Nþ1ðNÞ. Correspondingly, !xc¼
"Nþ1ðN þ 1Þ ! "Nþ1ðNÞ + "Nþ1ðNÞ ! "Nþ1ðNÞ ¼ 0 and Eg¼ !KSþ !xc + !KS. Therefore, to obtain a nonvan- ishing !xc, it is necessary to go beyond the zeroth-order theory.
To first order ð%1Þ in Eq. (12) (see the Supplemental Material [40]), the first-order correction to the orbital energy is
"ð1Þi ¼ hcð0Þi j H0ð0Þjcð0Þi i ¼ hcij H0%¼0jcii; (14) and the first-order correction to the orbital is
cð1Þi ðrÞ ¼X
j!i
hcð0Þj j H0ð0Þjcð0Þi i
"ð0Þi ! "ð0Þj
cð0Þj ðrÞ: (15)
Note that ~!%¼0ðrÞ¼PNþ1
i¼1 j ~c%i¼0ðrÞj2¼PNi¼1þ1jcð0Þi ðrÞj2¼ PNþ1
i¼1 jciðrÞj2¼!ðrÞþ!LðrÞ, where !LðrÞ ) jcNþ1ðrÞj2 is the KS LUMO density of the N-electron system.
Consequently, we have H0%¼0¼ e2Z !Lðr0Þ
jr ! r0jdr0þ vxcð½! þ !L#; rÞ ! vxcð½!#; rÞ:
(16) As "Nþ1ðNþ1Þ¼"%N¼1þ1ðNþ1Þ+"ð0ÞNþ1þ"ð1ÞNþ1, we have
!xc¼ "Nþ1ðN þ1Þ!"Nþ1ðNÞ + "ð1ÞNþ1
¼ hcNþ1j H0%¼0jcNþ1i
¼ e2ZZ !LðrÞ!Lðr0Þ
jr!r0j drdr0þZ
!LðrÞfvxcð½!þ!L#;rÞ
!vxcð½!#;rÞgdr (17)
and Eg¼ !KSþ !xc + !KSþ "ð1ÞNþ1. Equation (17) is a key result, showing that the DD can be approximately expressed as an explicit universal functional of !ðrÞ and
!LðrÞ and can be calculated in the standard KS method without extra computational cost. Note that it can also be derived from Eq. (6) by assuming (‘‘frozen orbital approxi- mation’’) that
c~iðrÞ +ciðrÞ; i¼ 1; 2; 3; . . . : (18) To second order (%2) in Eq. (12) (see the Supplemental Material [40]), the second-order correction to the orbital energy is
"ð2Þi ¼ hcð0Þi j H0ð0Þjcð1Þi i þ hcð0Þi j H0ð1Þjcð0Þi i
¼X
j!i
jhcjj H0%¼0jciij2
"iðNÞ ! "jðNÞ þ hcij H0ð1Þjcii; (19) where
H0ð1Þ¼@H%0
@%
"
"
"
"
"
"
"
"%¼0
¼Z ! e2
jr ! r0jþ$vxcð½~!%#; rÞ
$ ~!%ðr0Þ
#@ ~!%ðr0Þ
@%
"
"
"
"
"
"
"
"%¼0
dr0
¼Z ! e2
jr ! r0jþ fxcð½! þ !L#; r; r0Þ
#
,
! 2NXþ1
i¼1<½c-iðr0Þcð1Þi ðr0Þ#
#
dr0: (20)
Here, fxcð½!#; r; r0Þ ) $vxcð½!#; rÞ=$!ðr0Þ is the XC ker- nel, the asterisk denotes a complex conjugate, and<½* * *#
denotes the real part of½* * *#. From Eq. (19), we have
"ð2ÞNþ1¼ X
j!Nþ1
jhcjj H%0¼0jcNþ1ij2
"Nþ1ðNÞ ! "jðNÞ
þ hcNþ1j H0ð1ÞjcNþ1i: (21) Correspondingly, "Nþ1ðN þ 1Þ ¼ "%N¼1þ1ðN þ 1Þ + "ð0ÞNþ1þ
"ð1ÞNþ1þ "ð2ÞNþ1. This gives !xc¼"Nþ1ðNþ1Þ!"Nþ1ðNÞ+
"ð1ÞNþ1þ"ð2ÞNþ1and Eg¼!KSþ!xc+!KSþð"ð1ÞNþ1þ"ð2ÞNþ1Þ.
Extending the process further, the HOMO energy of the (Nþ 1)-electron system can be obtained by summing up all the perturbation corrections to infinite order, i.e.,
"Nþ1ðNþ1Þ¼"%N¼1þ1ðNþ1Þ¼P1
n¼0"ðnÞNþ1. Therefore, we can, in principle, obtain the exact !xc¼"Nþ1ðNþ1Þ!
"Nþ1ðNÞ¼P1
n¼1"ðnÞNþ1 and the exact Eg¼!KSþ!xc¼
!KSþP1
n¼1"ðnÞNþ1.
For any finite-order truncation of the above perturbation series, if two or more unperturbed states share the same energy, degenerate perturbation theory may be needed [36]. Since the concept of the perturbation to the unper- turbed Hamiltonian HKSremains valid, this scheme could be extended to estimate the !xc(and hence the Eg) for the degenerate cases based on the corresponding degenerate perturbation theory.
As mentioned previously, the DD needs to be summed to the KS gap to give the fundamental gap. While the DD given by Eq. (4) should be the same as that given by Eq. (6) for the exact functional, this property may no longer hold true for an approximate functional. For example, for a LDA or a generalized gradient approximation, while the DD
given by Eq. (4) is shown to vanish [13,21,30], we empha- size that the DD can be favorably restored by Eq. (6) and subsequently approximated by Eq. (17). Although a more accurate approximation for the DD could be pursued by higher-order perturbation theory, we adopt the DD given by Eq. (17) (i.e., first-order correction) for simplicity.
Accordingly, the fundamental gap is predicted by summing Eqs. (2) and (17) in our !KSþ !xc scheme.
Here, we examine the performance of various schemes in the prediction of the fundamental gaps for the FG115 database [20], which consists of 115 accurate reference values for the fundamental gaps of 18 atoms and 97 molecules at their experimental geometries. The funda- mental gaps are calculated by our !KSþ !xc scheme, the
!KSscheme [by Eq. (2)], and the Egscheme [by Eq. (5)], using the LDA [41] and LB94 [42] functionals and the 6-311++G(3df,3pd) basis set, with a development version of Q-CHEM3.2 [43]. The error for each entry is defined as (error¼ theoretical value ! reference value). The nota- tion used for characterizing statistical errors is as follows:
mean signed errors (MSEs), mean absolute errors (MAEs), and root-mean-square (rms) errors. Note that, for the !KS or !KSþ !xc schemes, only one KS calculation for the N-electron system is required (i.e., applicable to solids), while, for the Egscheme, which is the !KSþ !xc scheme with !xc being exactly calculated by Eq. (6) (with no further approximations), two KS calculations for the N- and (Nþ 1)-electron systems are required (i.e., inap- plicable to solids).
The calculated gaps are plotted against the reference values in Fig. 1 (for LDA) and Fig. 2 (for LB94). For both functionals, as the !KSgaps are shown to be vanish- ingly small (some of them are even negative) for the small- gap (smaller than 10 eV) systems, the DDs are essential for the accurate prediction of the fundamental gaps. In fact,
0 5 10 15 20 25 30
0 5 10 15 20 25 30
Calculated Fundamental Gap (eV)
Reference Fundamental Gap (eV)
∆KS
∆KS + ∆xc Eg
FIG. 1 (color online). Calculated versus reference fundamental gaps for the FG115 database [20]. The fundamental gaps are calculated by three schemes (see the text for details) using the LDA functional.
even for the large-gap (larger than 15 eV) systems, the DDs remain significant fractions of the fundamental gaps. As shown in TableI, the MAE associated with the !KSþ !xc
or Eg schemes is more than three times smaller than that associated with the !KSscheme [40].
Due to the use of the frozen orbital approximation [see Eq. (18)] in the evaluation of the DD, the !KSþ !xcgaps tend to be larger than the Eggaps. For rare gas atoms (e.g., He, Ne, and Ar), where this approximation becomes excel- lent, the !KSþ !xcgaps are very close to the Eggaps [40].
For LB94, the Eg gaps are in excellent agreement with the reference values, due to the correct asymptote of the LB94 potential, which is a key factor for the accurate prediction of the HOMO energies [28,37,38,42] and, via Eq. (5), the fundamental gaps. By contrast, for LDA, the Eg gaps tend to underestimate the reference values, due to the imbalanced self-interaction errors (as the LDA potential is asymptotically incorrect) in the predicted HOMO energies of the N- and (Nþ 1)-electron systems [2,4,28,37,38,42].
As the !KSþ !xcgaps tend to overestimate the Eggaps, it appears that there is a fortuitous cancellation of errors in the predicted !KSþ !xc gaps, when compared with the reference values.
In conclusion, we have provided a systematic procedure for the direct evaluation of the DD, based on perturbation theory. The lowest-order estimate of the DD is an explicit universal functional of the ground-state density and the KS LUMO density [see Eq. (17)], presenting a simple, effi- cient, and nonempirical scheme for the direct evaluation of the DD in the standard KS method. The fundamental gap can be accurately predicted by the sum of the KS gap and the estimated DD. The validity and accuracy of this scheme have been demonstrated for a wide variety of atoms and molecules, extending the applicability of KS DFT to an area long believed to be beyond its reach. To further improve the accuracy of this scheme, a more accurate
functional and a more accurate approximation for the DD (based on higher-order perturbation theory) should be adopted, although this will necessarily be somewhat more expensive. Since the concepts of the DD and the perturbation to the unperturbed Hamiltonian HKS are still valid for solid-state systems, this scheme could be extended to estimate the DD (a correction to the KS band gap) for solids, where the prediction of accurate funda- mental band gaps is very challenging for KS DFT. Work in this direction is in progress.
This work was supported by the National Science Council of Taiwan (Grant No. NSC101-2112-M-002- 017-MY3), National Taiwan University (Grants No. 99R70304, No. 101R891401, and No. 101R891403), and the National Center for Theoretical Sciences of Taiwan. We would like to thank Professor Shih-I Chu for providing us with Refs. [34,35].
*To whom all correspondence should be addressed.
[1] P. Hohenberg and W. Kohn,Phys. Rev. 136, B864 (1964);
W. Kohn and L. J. Sham,Phys. Rev. 140, A1133 (1965).
[2] S. Ku¨mmel and L. Kronik,Rev. Mod. Phys. 80, 3 (2008).
[3] A. J. Cohen, P. Mori-Sa´nchez, and W. Yang,Science 321, 792 (2008).
[4] A. J. Cohen, P. Mori-Sa´nchez, and W. Yang,Chem. Rev.
112, 289 (2011).
[5] J.-D. Chai,J. Chem. Phys. 136, 154104 (2012).
[6] M. S. Hybertsen and S. G. Louie,Phys. Rev. B 34, 5390 (1986).
[7] R. W. Godby, M. Schlu¨ter, and L. J. Sham,Phys. Rev. Lett.
56, 2415 (1986).
[8] A. Seidl, A. Go¨rling, P. Vogl, J. A. Majewski, and M. Levy, Phys. Rev. B 53, 3764 (1996).
[9] J. Heyd, G. E. Scuseria, and M. Ernzerhof,J. Chem. Phys.
118, 8207 (2003).
[10] Y. M. Niquet and X. Gonze, Phys. Rev. B 70, 245115 (2004).
[11] M. Gru¨ning, A. Marini, and A. Rubio,J. Chem. Phys. 124, 154108 (2006).
[12] P. Mori-Sa´nchez, A. J. Cohen, and W. Yang, Phys. Rev.
Lett. 100, 146401 (2008).
[13] A. J. Cohen, P. Mori-Sa´nchez, and W. Yang,Phys. Rev. B 77, 115123 (2008).
[14] F. Tran and P. Blaha,Phys. Rev. Lett. 102, 226401 (2009).
0 5 10 15 20 25 30
0 5 10 15 20 25 30
Calculated Fundamental Gap (eV)
Reference Fundamental Gap (eV)
∆KS
∆KS + ∆xc Eg
FIG. 2 (color online). Same as Fig. 1 but using the LB94 functional.
TABLE I. Statistical errors (in eV) of the 115 fundamental gaps of the FG115 database [20], calculated by three schemes (see the text for details) using the LDA and LB94 functionals.
Error
LDA LB94
!KS !KSþ !xc Eg !KS !KSþ !xc Eg
MSE !7:22 0.78 !2:33 !7:20 2.73 0.03
MAE 7.22 2.11 2.33 7.20 2.74 0.47
rms 7.44 2.45 2.57 7.48 3.13 0.63
[15] T. Stein, H. Eisenberg, L. Kronik, and R. Baer,Phys. Rev.
Lett. 105, 266802 (2010).
[16] T. Tsuneda, J.-W. Song, S. Suzuki, and K. Hirao,J. Chem.
Phys. 133, 174101 (2010).
[17] X. Zheng, A. J. Cohen, P. Mori-Sa´nchez, X. Hu, and W.
Yang,Phys. Rev. Lett. 107, 026403 (2011).
[18] M. Jain, J. R. Chelikowsky, and S. G. Louie, Phys. Rev.
Lett. 107, 216806 (2011).
[19] X. Andrade and A. Aspuru-Guzik, Phys. Rev. Lett. 107, 183002 (2011).
[20] Y.-S. Lin, C.-W. Tsai, G.-D. Li, and J.-D. Chai,J. Chem.
Phys. 136, 154109 (2012).
[21] W. Yang, A. J. Cohen, and P. Mori-Sa´nchez, J. Chem.
Phys. 136, 204111 (2012).
[22] J. F. Janak,Phys. Rev. B 18, 7165 (1978).
[23] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz Jr., Phys. Rev. Lett. 49, 1691 (1982).
[24] J. P. Perdew and M. Levy, Phys. Rev. Lett. 51, 1884 (1983).
[25] L. J. Sham and M. Schlu¨ter, Phys. Rev. Lett. 51, 1888 (1983).
[26] L. J. Sham and M. Schlu¨ter, Phys. Rev. B 32, 3883 (1985).
[27] W. Kohn,Phys. Rev. B 33, 4331 (1986).
[28] M. E. Casida,Phys. Rev. B 59, 4694 (1999).
[29] G. K.-L. Chan,J. Chem. Phys. 110, 4710 (1999).
[30] D. J. Tozer and N. C. Handy,Mol. Phys. 101, 2669 (2003).
[31] E. Sagvolden and J. P. Perdew,Phys. Rev. A 77, 012517 (2008).
[32] P. Mori-Sa´nchez, A. J. Cohen, and W. Yang, Phys. Rev.
Lett. 102, 066403 (2009).
[33] T. Stein, L. Kronik, and R. Baer,J. Am. Chem. Soc. 131, 2818 (2009).
[34] M. Lein and S. Ku¨mmel, Phys. Rev. Lett. 94, 143003 (2005); M. Mundt and S. Ku¨mmel, Phys. Rev. Lett. 95, 203004 (2005).
[35] D. A. Telnov, J. T. Heslar, and S.-I. Chu,Chem. Phys. 391, 88 (2011).
[36] D. J. Griffiths, Introduction to Quantum Mechanics (Prentice Hall, Upper Saddle River, NJ, 2004), 2nd ed.
[37] M. Levy, J. P. Perdew, and V. Sahni,Phys. Rev. A 30, 2745 (1984).
[38] C.-O. Almbladh and U. von Barth,Phys. Rev. B 31, 3231 (1985).
[39] J. P. Perdew and M. Levy,Phys. Rev. B 56, 16021 (1997).
[40] See Supplemental Material at http://link.aps.org/
supplemental/10.1103/PhysRevLett.110.033002 for the full derivation and further numerical results.
[41] P. A. M. Dirac, Proc. Cambridge Philos. Soc. 26, 376 (1930); J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
[42] R. van Leeuwen and E. J. Baerends,Phys. Rev. A 49, 2421 (1994).
[43] Y. Shao et al.,Phys. Chem. Chem. Phys. 8, 3172 (2006).
