More-Realistic Band Gaps from Meta-Generalized Gradient Approximations: Only in a Generalized Kohn-Sham Scheme

More-Realistic Band Gaps from Meta-Generalized Gradient Approximations: Only in a Generalized Kohn-Sham Scheme


Unlike the local density approximation (LDA) and the generalized gradient approximation (GGA), calculations with meta-generalized gradient approximations (meta-GGA) are usually done according to the generalized Kohn-Sham (gKS) formalism. The exchange-correlation potential of the gKS equation is non-multiplicative, which prevents systematic comparison of meta-GGA bandstructures to those of the LDA and the GGA. We implement the optimized effective potential (OEP) of the meta-GGA for periodic systems, which allows us to carry out meta-GGA calculations in the same KS manner as for the LDA and the GGA. We apply the OEP to several meta-GGAs, including the new SCAN functional [Phys. Rev. Lett. 115, 036402 (2015)]. We find that the KS gaps and KS band structures of meta-GGAs are close to those of GGAs. They are smaller than the more realistic gKS gaps of meta-GGAs, but probably close to the less-realistic gaps in the band structure of the exact KS potential, as can be seen by comparing with the gaps of the EXX+RPA OEP potential. The well-known grid sensitivity of meta-GGAs is much more severe in OEP calculations.

71.15 Mb, 31.15 es

I Introduction

Semiconductor devices play an important role in modern technologies, and the rapid development of electronic structure theory methods has made computational design of such devices possible. The band gap and the band structure are undoubtly the most important properties of semiconductors, since these are the properties that distinguish semiconductors from other periodic systemsYu and Cardona (2010). Computational evaluation of the band gap and the band structure is thus a topic of active research.

The fundamental band gap is a ground-state property, and it is defined as , where I is the ionization energy, and A is the electron affinity. I and A are ground-state energy differences. is also an excited-state property since it is the unbound limit of the exciton series. is very difficult to calculate for periodic systems, since there is no systematic way of adding/removing one electron to/from the solid in a periodic calculation, and the bulk limit can only be approached by the calculation of very big clusters. Many-body methods such as the GW methodHedin (1965) calculate and the quasiparticle band structure accurately, but the computational cost is high.

The density-functional theory (DFT)Hohenberg and Kohn (1964); Kohn and Sham (1965); Fiolhais et al. (2003) is a formally exact electronic structure method for the ground-state energy and electron density with an excellent balance of accuracy and computational efficiency, which is achieved by mapping the real interacting system to a fictitious Kohn-Sham (KS) system of non-interacting electrons with a multiplicative effective exchange-correlation (xc) potential (the functional derivative of the exchange-correlation energy with respect to the density). The exact Kohn-Sham potential yields the exact density but not the exact quasi-particle band structure and gap. Though the exact energy functional of the DFT is unknown, there exists a plethora of approximations, which has been ordered into the ‘Jacob’s ladder’Perdew and Schmidt (2001) hierarchy. The first and the second rungs of the Jacob’s ladder are the local density approximation (LDA) and the generalized gradient approximation (GGA), and they severely underestimate the fundamental gap in periodic systems. For periodic systems, KS DFT cannot calculate from its definition, and one commonly approximates with the KS gap , where and are the KS orbital energies of the highest occupied orbital and of the lowest unoccupied orbital, respectively. However, is not equal to even with the exact functional, due to the derivative discontinuity (DD)Perdew et al. (1982). The band gap problem has been an obstacle in the application of DFT to periodic systems.

The generalized Kohn-Sham (gKS)Seidl et al. (1996) scheme is a different formulation of the DFT, which allows a non-multiplicative but still Hermitian xc potential operator. The gKS gap can be a better approximation to than is the KS gapCasida (1999). The third rung of the Jacob’s ladder, the meta-generalized gradient approximation (meta-GGA), is commonly implemented in the gKS scheme according to the method of Neumann, Nobes and Handy (also denoted as gKS in this paper)Neumann et al. (1996). The gKS meta-GGA band gap of periodic systems improvesZhao and Truhlar (2009) over the KS GGA gaps as expected. In this work we find that, with the recently proposed strongly constrained and appropriately normed (SCAN) functionalSun et al. (2015, ), the gKS meta-GGA gap corrects about 20%50% of the difference between the experimental fundamental gap and the GGA KS gap.

Due to the restriction in the functional form of GGA, a GGA cannot perform well for finite systems and periodic systems at the same time.Perdew et al. (2008) On the other hand, the functional form of meta-GGA can satisfy more exact conditions and has a wider range of applicability than the GGA form. The SCAN functional is a non-empirical functional that satisfies all known exact constraints appropriate to a semilocal functional, and is expected to perform well for systems with very different kinds of bonds. The computational accuracies of the SCAN functional for many properties improve over those of the GGAs, with only marginal increase of computational cost.Sun et al. (2015) We find that SCAN also improves band gaps, but the comparsion is between meta-GGA gKS gaps and GGA KS gaps, which are not the same quantity. It is unclear whether the KS gap itself is improved, or just the gKS gap is improved. One needs to do meta-GGA calculations in the KS scheme to be able to have a systematic comparsion of the KS band gaps between meta-GGA and GGA.

In the gKS formalism for an orbital functional such as a meta-GGA or hybrid functional, we find the “optimal variational potential”, a non-multiplicative self-adjoint operator that minimizes the energy functional with respect to the orbitals. This potential is a differential operator for a meta-GGA and an integral operator for a hybrid functional. In the KS formalism for an orbital functional, we find the optimized effective potential (OEP)Sharp and Horton (1953); Talman and Shadwick (1976), the multiplicative xc potential that minimizes the energy. The OEP meta-GGA yields the KS gap of the meta-GGA, which can be compared directly with the GGA KS gap. With the OEP, the meta-GGA xc potential also becomes comparable with that of the GGA. The OEP meta-GGA has been studied only in finite systems previously.Arbuznikov and Kaupp (2003); Eich and Hellgren (2014)

In this work, we provide the first study of the OEP meta-GGA in periodic systems. We find that the meta-GGA KS gap is not significantly improved over the GGA KS gaps, and the band structure of the OEP meta-GGA is very close to that of the GGA (and presumably to that of the exact KS potential). The xc potential of the meta-GGA has more details than that of the GGA, but the change is not as big as the change between the LDA and the GGA. Though the gKS meta-GGA is known to be sensitive to the real-space grid used in the calculation, we find that the OEP meta-GGA has worse grid sensitivity. The reason for the grid sensitivity is discussed in this work.

Within the gKS implementation of the meta-GGA form, it is possible to fit to energy gaps of solidsTran and Blaha (2009), and the result can be useful for the prediction of gapsKosobutsky and Gordienko (2015). But the fact that a property can be fitted in DFT is not evidence that it should be, or that other and more appropriate properties will not deteriorate as a result.

Ii OEP meta-GGA in periodic systems

ii.1 Theory

The meta-GGA xc energy functional has the form


where is the xc energy density, and and are the spin density and the kinetic energy density of spin respectively. and are


where is the volume of the first Brillouin zone (BZ), is the band index, is the KS orbital normalized in one unit cell, and is the Heaviside step function, with being the Fermi energy, and being the KS orbital energy. The non-interacting kinetic energy density of the Kohn-Sham orbitals is, like the orbitals themselves, a functional of the electron density. The OEP of the meta-GGA is defined by the functional derivative of the xc energy with respect to the density, which is


Eq. (4) can be partitioned into a GGA-like part and a -dependent part:


is a multiplicative potential, but does not have a closed form. can be written as


On one hand, the gKS meta-GGA potential operator is obtained by multiplying both side of Eq. (6) by and integrating over :


On the other hand, inserting


where is the KS potential, into Eq. (6) yields the OEP integral equation for :


Though a direct solution for the OEP is possibleIvanov et al. (1999); Eich and Hellgren (2014), it is both computationally expensive and numerically unstable. The common practice is to approximate the solution of Eq. (9). Ref. Arbuznikov and Kaupp, 2003 employed the local-Hartree-Fock (LHF)Della Sala and Goerling (2001) approximation to derive approximated of meta-GGA for finite systems. Here we do the derivation with the effective local potential (ELP)Staroverov et al. (2006); Izmaylov et al. (2007) approximation. The ELP is equivalent to the LHF and the common-energy-denominator approximation (CEDA)Gritsenko and Baerends (2001); Grüning et al. (2002), and the Krieger-Li-Iafrate (KLI)Krieger et al. (1992) approximation is their simplification.

The ELP minimizes the matrix norm of the commutator , where , and is the spin density matrix. becomes valid for both the KS and the gKS system when . is minimized when , which yields the ELP approximation of :


The ELP approximates the OEP in a least-squares senseIzmaylov et al. (2007).

Ref. Izmaylov et al., 2007 shows that the ELP approximation is equivalent to the CEDA, which reduces to the KLI approximation by setting in Eq. (10)Gritsenko and Baerends (2001). If only the first term of Eq. (10) is kept, one obtains the so-called Slater approximationSlater (1951).

ii.2 Implementation in the BAND code

The most efficient methods for calculating periodic systems are pseudopotential methods with planewave basis sets, with the widely used PAW method as an extensionBlöchl (1994); Kresse and Joubert (1999), but the OEP is hard to implement in such codes. The sum over bands in Eq. (10) also includes the core bands, which are not directly available in pseudopotential codes. Furthermore, the OEP is inherently incompatible with the PAW formalism. The charge density in the PAW formalism is composed of the pseudo-charge-density, the on-site all-electron density, and the on-site pseudo-charge-densityBlöchl (1994); Kresse and Joubert (1999), and the xc energy and potential are partitioned into three parts accordingly, with each part only depending on the corresponding part of the density. However, the OEP of meta-GGA is hard to partition this way. Due to these obstacles, we implement the OEP meta-GGA in the BAND codete Velde and Baerends (1991); Wiesenekker and Baerends (1991); Franchini et al. (2013, 2014); Philipsen et al. ().

BAND is an all-electron DFT code for periodic systems. It uses atom-centered numerical functions as the basis set. All the real-space integrations are done numerically on a Becke fuzzy-cell gridBecke (1988). The BZ is discretized with the commonly used Monkhorst-Pack gridMonkhorst and Pack (1976), so the -integrations of Eq. (10) become sums over -points in the irreducible Brillouin zone (IBZ):


where is the contribution to the integral from the orbital , is the occupation number, is the weight of the -point times the -space integration weight in BANDWiesenekker and Baerends (1991), is the total number of symmetry operations, is the space-group symmetry operator, and is the rotational part of .

The matrix elements in Eq. (10) vanish for , since has the perodicity of the unit cell. Eq. (10) is implemented in the BAND code as


where and are


and denotes the Hessian matrix of . The integrals are done as discrete sums on the real-space grid. The depend on , which is unknown beforehand. They can be determined by solving linear equationsArbuznikov and Kaupp (2003), but this is costly for perodic systems. Instead, we determine the integrals and iteratively by the following steps:

  1. Set the integrals to 0.

  2. Calculate with Eq. (12) using the integrals of the last iteration.

  3. Calculate the integrals using the new .

  4. Repeat step 2 and 3 if , where is the iteration count, and is the error tolerance.

is only determined up to an additive constant , since yields the same band structure. has to be fixed for determining the convergence. The integrals contain , so their values include . We fix by requiring that the term of containing the integrals averages to 0, i.e. the constant is determined by


In practice, the iterative procedure converges after 200 iterations on average for .

Ref. Grüning et al., 2002 finds that the total energies and orbital energies obtained from CEDA and KLI are very close for atoms and molecules. We find that the same is true for periodic systems. The ELP is equivalent to the CEDAIzmaylov et al. (2007), and its total energy is always lower than that of the KLI, but the total energy differences and the KS gap differences of the materials in Table 1 between that of the ELP and that of the KLI are at most 1 meV. Since the improvement of the ELP over the KLI is insignificant, we only report the KLI results in the following.

Iii Results

We calculate the KS and gKS gaps of 20 semiconductors and insulators with the SCAN, meta-GGA made simple 2 (MS2)Sun et al. (2012, 2013), meta-GGA made very simple (MVS)Sun et al. (2014) and TPSSTao et al. (2003) meta-GGAs and the PBEPerdew et al. (1996) GGA, and the results are collected in Table 1. Ge, CdO, and InN are not listed since their calculated gaps vanish. All the calculations are done with the TZ2P basis setvan Lenthe and Baerends (2003). A Monkhorst-Pack -gridMonkhorst and Pack (1976) is used for most of the materials except InN, CdS and CdSe, for which a -grid is used. We find that the OEP meta-GGA is sensitive to the real-space grid, and large grids are used to ensure convergence. The details of the grid problem are discussed in Section IV. Though GGA and gKS meta-GGA converge properly with smaller grids, we use the same grid as the OEP meta-GGA, so that the results are comparable.

The scalar relativistic effect is included in the calculations by the ZORAPhilipsen et al. (1997) method. We ignore the spin-orbit coupling since the effect is small (0.03 eV for InP). The relativistic effect is implicitly included through the pseudopotential in planewave codes, but it needs to be explicitly included in BAND. We find that the scalar relativistic effect has a big impact on the KS and gKS gaps, since it shifts down the orbital energies of -type bands, such as the lowest conduction band of GaAsBachelet and Christensen (1985). The change in the KS and gKS gaps due to the scalar relativistic effect can be as big as 0.7 eV (GaAs).

Figs. 1 and 2 compare the OEP of SCAN meta-GGA and the LDA and the PBE xc potentials for the Ne atom and bulk Si, respectively. The exact xc potential of the Ne atomRyabinkin and Staroverov (2012) is plotted as a reference. The exact xc potential of Ne has a bump between the two shells of the charge density. The of LDA does not have this feature. The ’s of PBE and SCAN both have this bump, and they are roughly in the same position as that of the exact , but these ’s are shallower than the exact . The ’s of SCAN also have a few small bumps that are not in the exact . In the asymptotic region, the exact decays as , and the ’s of all the semi-local functionals decay exponentially. Fig. 1 shows that the of SCAN has the same decay as that of LDA and PBE. For periodic systems, there is no asymptotic region, and the approximated ’s of bulk Si in Fig. 2 can be expected to be closer to the exact than those in Fig. 1. In Fig. 2, the of SCAN are similar to the of PBE, and they only differ in small details. Similar to Fig. 1, the ’s of SCAN for bulk Si also have small bumps. Though the meta-GGA is a higher rung functional on the Jacob’s ladder than the GGA, the improvement in is small going from GGA to meta-GGA, unlike going from LDA to GGA. The differences between the KS meta-GGA gaps and the PBE gaps in Table 1 are small as a consequence.

Figure 1: (color online) Comparison of ’s of the Ne atom for SCAN, LDA, and GGA. The exact is from Ref. Ryabinkin and Staroverov, 2012.
Figure 2: (color online) Comparison of ’s of bulk Si along the Si-Si bond for SCAN, LDA, and GGA. The Si atoms are located at and Å. The vertical dashed line is a numerical artifact and does not affect the band structure.
Si 1.17 0.60 0.71 1.11 0.97 0.78 0.19 1.20 0.80 0.41 1.04 0.72 0.32 0.80
InP 1.42 0.50 0.72 1.52 1.06 0.77 0.29 1.14 0.81 0.34 1.99 1.09 0.89 0.86
GaAs 1.52 0.30 0.53 1.41 0.8 0.45 0.34 0.94 0.48 0.46 2.151 1.26 N/A 0.68
BAs 1.602 1.21 1.26 1.71 1.51 1.32 0.19 1.63 1.33 0.30 1.56 1.16 0.40 1.26
CdSe 1.73 0.44 0.71 1.66 1.10 0.76 0.33 1.06 N/A3 N/A 2.14 0.50 1.64 0.92
BP 2.10 1.36 1.43 1.79 1.74 1.52 0.22 1.94 1.49 0.46 1.64 1.42 0.22 1.48
GaP 2.35 1.53 1.69 2.09 1.94 1.72 0.21 1.97 1.73 0.24 2.23 1.74 0.49 1.74
CdS 2.48 0.96 1.23 2.27 1.62 1.20 0.42 1.60 1.27 0.33 2.39 1.50 0.88 1.43
-GaN 3.17 1.70 1.69 2.97 2.03 1.84 0.20 1.69 1.70 -0.01 2.50 1.89 0.62 1.60
ZnS 3.72 1.87 2.12 3.32 2.63 2.16 0.47 2.52 2.08 0.45 3.35 N/A4 N/A 2.30
C5 5.50 4.14 4.17 4.94 4.58 4.26 0.32 4.79 4.16 0.62 4.15 4.06 0.09 4.20
BN 6.20 4.42 4.53 5.39 5.04 4.73 0.30 5.01 4.56 0.45 5.14 4.58 0.56 4.54
CaO 6.93 3.62 3.75 5.30 4.29 3.78 0.50 4.13 N/A6 N/A 4.56 3.51 1.05 3.89
MgO 7.90 4.70 4.74 6.46 5.62 4.80 0.82 5.20 5.47 -0.27 6.05 4.58 1.47 4.86
NaCl 8.97 4.70 5.08 6.42 5.86 5.25 0.61 5.85 6.51 -0.66 6.61 4.68 1.93 5.43
LiF 13.6 8.84 9.04 11.4 9.97 9.11 0.86 9.507 9.84 N/A 10.64 8.77 1.86 9.19
solid Ar 14.3 8.44 8.92 10.33 9.91 8.89 1.02 9.95 9.36 0.58 10.98 9.19 1.80 9.56
MAE 2.08 1.90 0.88 1.41 1.84 1.28 1.63 1.13 1.89 1.76
MARE 0.46 0.40 0.13 0.27 0.38 0.26 0.35 0.19 0.36 0.36
Table 1: Calculated KS/gKS gaps and the derivative discontinuities (). All energies are in eV. The experimental gaps are from Ref. Madelung, 2004. The BAs GW gap is from Ref. Surh et al., 1991. The mean absolute errors (MAE) and mean absolute relative error (MARE) for the band gaps are listed. HSEHeyd et al. (2003, 2005, 2006) gaps calculated with the VASPKresse and Hafner (1993, 1994); Kresse and Furthmüller (1996a, b); Kresse and Joubert (1999) code are also listed for reference. The PBE and SCAN gaps of Ge, CdO, and InN are vanishing, so they are not listed. Their experimental gaps are 0.74 eV, 0.84 eV, and 1.95 eV, respectively. Their HSE gaps are 0.65 eV, 0.94 eV, and 0.77 eV, respectively. The KS(KLI) TPSS results are not listed as most calculations fail to converge.
Figure 3: The absolute errors of PBE, KS(KLI) SCAN, and gKS SCAN (comparing with the experimental gap) in the band gap. The EXX+RPA OEP gaps for Si, LiF and Ar, from Ref. Grüning et al., 2006, can be considered very close to the gaps for the corresponding exact KS potential.

The band structures of Si and GaAs calculated with PBE and SCAN are plotted in Figs. 4 and 5. The KS(KLI) SCAN band structure is very close to the PBE band structure, due to the corresponding being similar to the PBE . The gKS SCAN band structure has the same overall shape as that of the PBE and the KS(KLI), and the main difference is in the band gap.

Though the gKS meta-GGA band gaps improve over the PBE gaps in general, it is disappointing that gKS meta-GGA gaps for Ge, InN, and CdO still vanish. However, it is possible for meta-GGAs to open the gap for gapless materials in GGA. gKS SCAN hasKitchaev et al. (2016) a 0.4 eV gap for , which is gapless in GGA, and the value is close to the experimental value 0.3 eV. The M06L metaGGA was reported to open the gap of Ge at 0.14 eV.Zhao and Truhlar (2009); Peverati and Truhlar (2012)

The improvement of the band gap occurs since, unlike the KS gap, the gKS gap is an approximation to the fundamental gap of the meta-GGA. A Janak-typeJanak (1978) theorem has been proven for the OEPCasida (1999), and it states that the gKS gap approximately equals the fundamental gap for the same functional, assuming fixed orbitals. This assumption does not apply to finite systems, but it is true for periodic systems, since the charge density and the orbitals of a periodic system undergo only an infinitesimal change when the number of electrons changes by one.

GGA band gaps should be compared with the OEP meta-GGA band gaps for a fair comparison between approximated functionals, since the OEP meta-GGA band gap is the KS gap. The SCAN functional is the only functional that satisfies all the known exact conditions, but the KS(KLI) SCAN gaps do not have significant improvements over the PBE gaps. This is probably due to the fact that the GGA and SCAN OEP gaps closely approximate the exact KS gap, which underestimates the fundamental gap. This has been illustrated in Fig. 3, where the errors of the EXX+RPA(OEP) KS gapsGrüning et al. (2006) are also plotted. EXX+RPA (exact exchange plus random phase approximation for correlation) is a high-level (fifth rung) method, and its OEP gaps are expected to be very close to those of the corresponding exact KS potential. Fig. 3 shows that both PBE and KS SCAN gaps are already good approximations to the exact KS gap.

Some of the gKS band gaps of MS2 and TPSS are smaller than the corresponding OEP band gaps. We do not find this behavior in other functionals. Many of the KS(KLI) TPSS calculations fail to converge. This is probably a numerical issue in the calculation of , due to the complicated functional form of TPSS.

Figure 4: (color online) The band structure of Si calculated with PBE, gKS SCAN, and KS(KLI) SCAN.
Figure 5: (color online) The band structure of GaAs calculated with PBE, gKS SCAN, and KS(KLI) SCAN.

The energy functional of the exact DFT has derivative discontinuities at integer electron numbersPerdew et al. (1982), where . The exact KS potential jumps up by the positive constant as the electron number crosses the value that makes the solid electrically neutral. LDA and GGA miss much or all of the derivative discontinuity due to the convexity of the approximated energy functionalMori-Sánchez et al. (2008), and they underestimate the gap as a consequence. In periodic systems, such discontinuities are only visible at band gaps. Though the derivative discontinuity has been thoroughly studied for atoms and molecules, it is difficult to obtain for periodic systems, since and cannot be calculated directly.

With the OEP meta-GGA provided in this work, we are able to estimate the of solids. can be approximated by the derivative gap , which is equal to the gKS meta-GGA gapCasida (1999); Mori-Sánchez et al. (2008); Yang et al. (2012). is the OEP meta-GGA gap. The results are shown in Table 1.

Comparison of the ’s of meta-GGAs with the exact is impossible for periodic systems, since no exact KS potential is available. However, Ref. Godby et al., 1986 provides an OEP of the GW method of bulk Si, and its can be seen as a good approximation to the exact . The of bulk Si in Ref. Godby et al., 1986 is 0.58 eV, which is larger than that of all the meta-GGAs in Table 1. This is not unexpected because the meta-GGA’s tested in this work are not exact for in a solid.

Iv Real-space grid dependency

Material solid Ar LiF NaCl MgO CaO BN C8 ZnS -GaN CdS GaP BP CdSe BAs GaAs InP Si
69 169 82 74 87 74 65 153 209 96 91 74 96 87 169 177 137
0.06 0.09 0.22 0.59 0.54 0.27 0.17 0.0008 0.03 0.17 0.05 0.03 0.14 0.01 0.03 0.06 0.10
537 705 638 571 672 571 504 638 873 739 705 571 739 672 672 739 571
-0.003 -0.01 -0.005 -0.002 -0.01 -0.02 -0.02 -0.006 -0.03 -0.0007 -0.002 -0.02 -0.05 -0.02 -0.01 -0.001 -0.02
Table 2: The total energy differences between KS(KLI) SCAN and SCAN with PBE orbitals of different grids. is the number of radial grid points. . All energies are in eV.

The gKS meta-GGA requires a larger real-space integration grid than GGAJohnson et al. (2009), and this is caused by the sharp variation in regions far away from nuclei of quantities containing , such as used in TPSS, and used in SCAN, MS2 and MVS, where is the von Wiezsäcker kinetic energy density, and is the kinetic energy density of the uniform electron gas. The Becke fuzzy-cell grid is used as the real-space integration grid in most of the DFT codes. It is constructed by combining atom-centered spherical grids. The radial part of the spherical grid is dense near the nuclei, and is sparse away from the nuclei. The integration weights for the grid points in the sparse region would be larger than those in the dense region. Therefore a function with sharp features in regions away from the nuclei cannot be properly represented on the grid, and the error of integrations involving this function would be large. One needs larger grids in gKS meta-GGA calculations than those used in GGA calculations.

We find that the OEP meta-GGA is more sensitive to the real-space integration grid than the gKS meta-GGA. In the gKS case, the sensitivity to the grid shows up in the potential energy surfaceJohnson et al. (2004, 2009) as spurious oscillations, but the sensitivity is not obvious in a single-point calculation. In the OEP case, a single-point calculation is enough to demonstrate the grid sensitivity by comparing the total energy. The total energy of the OEP meta-GGA is variationally minimized with respect to the charge density, but the total energy of the gKS meta-GGA is variationally minimized with respect to the orbitals. Since the gKS has bigger variational degrees of freedom, the gKS total energy should always be lower than the OEP total energy, and both should be lower than the meta-GGA total energy evaluated with non-variational orbitals. However, though the built-in integration grids in the BAND code are good enough for gKS calculations, they are not sufficient for an OEP calculation: the OEP total energy calculated with these grids is always higher than the total energy evaluated with PBE orbitals. For example, the best built-in grid in BAND for the Ne atom has 120 radial points, but at least 266 radial points are required for the SCAN OEP total energy to be lower than the SCAN total energy evaluated with PBE orbitals, and 504 radial points are required for the SCAN OEP total energy to converge with respect to the grid (error0.01 meV). The total energy is not sensitive to the number of angular grid points. Table 2 lists the total energy differences using different grids.

The OEP of the Ne atom, dimer, and bulk Si (along the Si-Si bond) are plotted in Figs. 6, 7, and 8. The OEP develops unnatural peaks close to nuclei when too few grid points are used. These peaks have a strong effect on the KS band gap and the total energy. For bulk Si with 134 points (dashed line in Fig. 8), the system becomes gapless.

The grid-dependence of the OEP meta-GGA is caused by in Eq. (12). Figs. 6, 7, and 8 show that has very sharp oscillations both close to and away from the nucleus. The dots in these plots show the grid point locations when the calculation is done with a built-in grid in BAND. These grid points are sufficient for an gKS meta-GGA calculation since they are dense enough to properly describe the oscillations in , but they are insufficient for an OEP meta-GGA calculation since the peaks in are much narrower than those of .

Figure 6: (color online) The KS(ELP) SCAN (a.u.), the radial component of (a.u.), and of the Ne atom evaluated with two different grids. The solid lines are results obtained by using 504 radial grid points, and the dots and dashed lines are results obtained by using 65 radial grid points.
Figure 7: (color online) The KS(KLI) SCAN (a.u.), the component of (a.u.), and of the Ar dimer evaluated with two different grids, plotted along the Ar-Ar axis. The Ar atoms are located at Å, so that comparison with Ref. Johnson et al., 2009 is possible. Since the system is symmetric, only the part is plotted. The solid lines are results obtained by using 806 radial grid points, and the dots and dashed lines are results obtained by using 134 radial grid points.
Figure 8: (color online) The KS(KLI) SCAN (a.u.), the component of (a.u.), and of bulk Si evaluated with two different grids, plotted along the Si-Si bond. The Si atoms are located at and Å. The solid lines are results obtained by using 571 radial grid points, and the dots and dashed lines are results obtained by using 134 radial grid points.

Fig. 9 shows of the Ne atom for SCAN, TPSS, MS2, and MVS functionals. Though of all functionals have oscillations, the sharpness of the peaks is different, and consequently the numbers of grid points required for convergence are also different: SCAN has very sharp peaks, and it requires 352 points to converge the total energy within 1% error; The peaks in MVS are broader, and it only requires 252 points to converge with the same error criterion.

The oscillations in for SCAN, MS2 and MVS are centered at or close to in Fig. 9. All these functionals use to incorporate the kinetic energy density into the functional, and their energy densities all have the form , where and are constructed for and , and is an interpolation function that decreases monotonically with from 1 at to 0 at to negative values for . The oscillation in then implies a peak in . corresponds to the uniform electron gas limit, and the of SCAN and MS2 are constructed to have there to recover the gradient expansion approximationSun et al. (2012). Therefore has a flat- or linearly-topped peak at , which explains the oscillation. For MVS, vanishes at , and its oscillation in occurs there. By choosing other functional forms for , it is possible to get rid of the oscillations and the grid sensitivity. TPSS does not have this feature since it uses instead of .

Figure 9: (color online) The radial component of of the Ne atom of different functionals. ’s of different functionals are similar, so only that of SCAN is plotted.

Even though a small grid cannot properly represent the oscillations in , Fig. 10 shows that the of the Slater approximation (which contains explicitly) evaluated on a small grid is already similar to the KLI evaluated on a big grid. The KLI approximation is supposed to be an improvement over the Slater approximation, but it introduces big errors in when using a small grid. The grid sensitivity actually enters the OEP indirectly through the integrals of Eq. (14). Since the Slater term contains , the integrals cannot be done accurately with a small grid.

The gKS meta-GGA potential operator in Eq. (7) also contains , but the grid sensitivity of gKS is much lower than that of OEP. There is no contradiction since does not have to be evaluated directly in gKS meta-GGA. Using integration by parts, the matrix elements becomes


and is less sensitive to the grid than . The grid sensitivity of the OEP is also expected to show up in the time-dependent density-functional theory (TDDFT)Runge and Gross (1984) with a meta-GGA xc kernel.

Figure 10: (color online) KS(KLI) SCAN of bulk Si along the Si-Si bond. The Si atoms are located at and Å.

Aside from the grid problem, we find that the OEP meta-GGA in general needs more self-consistent-field (SCF) cycles to reach convergence than the gKS meta-GGA. For small gap materials, the OEPs of some meta-GGA functionals do not converge, while the corresponding gKS calculations converge normally. We tested SCAN, MS2, MVS, and TPSS functionals, and SCAN has the least convergence problem.

V Conclusion

In conclusion, we implemented the ELP and the KLI approximations for the OEP meta-GGA of periodic systems, with which we study the meta-GGA band gaps of 20 semiconductors and insulators. Comparing with the GGA band gaps, the new SCAN meta-GGA in a generalized Kohn-Sham scheme is found to improve the band gaps over the GGA band gaps. The non-empirical SCAN meta-GGA outperforms the non-empirical PBE GGA, not only for diversely-bonded materials near equilibriumSun et al. (2015), but also for the band gaps of solids. The improvement is achieved without using the expensive exact exchange, as in fourth-rung hybrid functionals. For materials wrongly predicted to be gapless in GGA, the resultKitchaev et al. (2016) of with SCAN shows that meta-GGAs can open the gap.

Consider the ratio of the calculated to the experimental energy gap. For the 17 solids of Table 1, this ratio varies from 0.35 to 0.79 with an average of 0.60 for the PBE GGA, from 0.53 to 0.94 with an average of 0.73 for the SCAN meta-GGA, and from 0.72 to 1.07 with an average of 0.89 for the HSE hybrid functional. The ratio improves uniformly from PBE to SCAN to HSE. While the hybrid functional is more accurate for the gap than SCAN is, it is also more empirical and more computationally expensive.

For periodic systems, the OEP or ungeneralized Kohn-Sham ’s of meta-GGAs are close to the GGA ’s, and they only differ in small details. Consequently the OEP meta-GGA gaps (like the OEP EXX+RPA gaps where available) are not improved significantly over the GGA gaps, and the band structures of OEP meta-GGAs are similar to those of GGAs. We think it is likely that the band gap and band structure in the exact Kohn-Sham potential are close to those of GGA and OEP meta-GGA (and of the OEP hybrid) in periodic systems. Aside from the band gap, the gKS meta-GGA band structures are also similar to the GGA band structures.

Due to the sharp features in , the OEP meta-GGA is sensitive to the real-space grid used in computation, more so than the gKS meta-GGA. Different meta-GGAs have different requirements for the minimal grid, and in general SCAN needs the biggest grid to converge of the meta-GGAs tested in this work. TDDFT with meta-GGA xc kernels is also expected to show this grid sensitivity. It is possible to avoid the grid sensitivity in the design of the functional, although this may interfere with the satisfaction of some exact conditions.


ZY, JS and JPP are supported by the National Science Foundation(Grant No. DMR-1305135). HP, JS, and JPP are also supported by the Center for the Computational Design of Functional Layered Materials, an Energy Frontier Research Center (EFRC) funded by the U.S.  Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0012575. The EFRC grant contributed to the computational aspects of this work. We thank Kieron Burke for suggestions on the manuscript.


  1. converged to the wrong ground state
  2. GW gap
  3. not converged
  4. not converged
  5. diamond
  6. not converged
  7. converged to the wrong ground state
  8. diamond


  1. P. Y. Yu and M. Cardona, Fundamentals of semiconductors: physics and materials properties (Springer, Berlin, 2010).
  2. L. Hedin, Phys. Rev. 139, A796 (1965).
  3. P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  4. W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  5. C. Fiolhais, F. Nogueira,  and M. Marques, eds., A primer in density functional theory, Lecture notes in physics (Springer, Berlin, 2003).
  6. J. P. Perdew and K. Schmidt, in Density functional theory and its application to materials, edited by V. van Doren, C. van Alsenoy,  and P. Geerlings (AIP Press, New York, 2001).
  7. J. P. Perdew, R. G. Parr, M. Levy,  and J. L. J. Balduz, Phys. Rev. Lett. 49, 1691 (1982).
  8. A. Seidl, A. Görling, G. P. Vogl, J. A. Majewski,  and M. Levy, Phys. Rev. B 53, 3764 (1996).
  9. M. E. Casida, Phys. Rev. B 59, 4694 (1999).
  10. R. Neumann, R. H. Nobes,  and N. C. Handy, Mol. Phys. 87, 1 (1996).
  11. Y. Zhao and D. G. Truhlar, J. Chem. Phys. 130, 074103 (2009).
  12. J. Sun, A. Ruzsinszky,  and J. P. Perdew, Phys. Rev. Lett. 115, 036402 (2015).
  13. J. Sun, R. C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky, H. Peng, Z. Yang, A. Paul, U. Waghmare, X. Wu, M. L. Klein,  and J. P. Perdew, Nature Chem., to appear.
  14. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou,  and K. Burke, Phys. Rev. Lett. 100, 136406 (2008).
  15. R. T. Sharp and G. K. Horton, Phys. Rev. 90, 317 (1953).
  16. J. D. Talman and W. F. Shadwick, Phys. Rev. A 14, 36 (1976).
  17. A. V. Arbuznikov and M. Kaupp, Chem. Phys. Lett. 381, 495 (2003).
  18. F. G. Eich and M. Hellgren, J. Chem. Phys. 141, 224107 (2014).
  19. F. Tran and P. Blaha, Phys. Rev. Lett. 102, 226401 (2009).
  20. A. V. Kosobutsky and A. B. Gordienko, Phys. Solid State 57, 1972 (2015).
  21. S. Ivanov, S. Hirata,  and R. J. Bartlett, Phys. Rev. Lett. 83, 5455 (1999).
  22. F. Della Sala and A. Goerling, J. Chem. Phys. 115, 5718 (2001).
  23. V. N. Staroverov, G. E. Scuseria,  and E. R. Davidson, J. Chem. Phys. 125, 081104 (2006).
  24. A. F. Izmaylov, V. N. Staroverov, G. E. Scuseria, E. R. Davidson, G. Stoltz,  and E. Cancès, J. Chem. Phys. 126, 084107 (2007).
  25. O. V. Gritsenko and E. J. Baerends, Phys. Rev. A 64, 042506 (2001).
  26. M. Grüning, O. V. Gritsenko,  and E. J. Baerends, J. Chem. Phys. 116, 6435 (2002).
  27. J. B. Krieger, Y. Li,  and G. J. Iafrate, Phys. Rev. A 45, 101 (1992).
  28. J. C. Slater, Phys. Rev. 81, 385 (1951).
  29. P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  30. G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  31. G. te Velde and E. J. Baerends, Phys. Rev. B 44, 7888 (1991).
  32. G. Wiesenekker and E. J. Baerends, J. Phys.: Cond. Matt. 3, 6721 (1991).
  33. M. Franchini, P. H. T. Philipsen,  and L. Visscher, J. Comput. Chem. 34, 1818 (2013).
  34. M. Franchini, P. H. T. Philipsen, E. van Lenthe,  and L. Visscher, J. Chem. Theory Comput. 10, 1994 (2014).
  35. P. H. T. Philipsen, G. te Velde, E. J. Baerends, J. A. Berger, P. L. de Boeij, M. Franchini, J. A. Groeneveld, E. S. Kadantsev, R. Klooster, F. Kootstra, P. Romaniello, D. G. Skachkov, J. G. Snijders, C. J. O. Verzijl, G. Wiesenekker,  and T. Ziegler, “Band2014 (a modified version is used),” SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands,
  36. A. D. Becke, J. Chem. Phys. 88, 2547 (1988).
  37. H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  38. J. Sun, B. Xiao,  and A. Ruzsinszky, J. Chem. Phys. 137, 051101 (2012).
  39. J. Sun, R. Haunschild, B. Xiao, I. W. Bulik, G. E. Scuseria,  and J. P. Perdew, J. Chem. Phys. 138, 044113 (2013).
  40. J. Sun, J. P. Perdew,  and A. Ruzsinszky, Proc. Nat. Acad. Sci. 112, 685 (2014).
  41. J. Tao, J. P. Perdew, V. N. Staroverov,  and G. E. Scuseria, Phys. Rev. Lett. 91, 146401 (2003).
  42. J. P. Perdew, K. Burke,  and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  43. E. van Lenthe and E. J. Baerends, J. Comput. Chem. 24, 1142 (2003).
  44. P. H. T. Philipsen, E. van Lenthe, J. G. Snijders,  and E. J. Baerends, Phys. Rev. B 56, 13556 (1997).
  45. G. B. Bachelet and N. E. Christensen, Phys. Rev. B 31, 879 (1985).
  46. I. G. Ryabinkin and V. N. Staroverov, J. Chem. Phys. 137, 164113 (2012).
  47. O. Madelung, Semiconductors: Data Handbook (Springer, 2004).
  48. M. P. Surh, S. G. Louie,  and M. L. Cohen, Phys. Rev. B 43, 9126 (1991).
  49. J. Heyd, G. E. Scuseria,  and M. Ernzerhof, J. Chem. Phys. 118, 8207 (2003).
  50. J. Heyd, J. E. Peralta, G. E. Scuseria,  and R. L. Martin, J. Chem. Phys. 123, 174101 (2005).
  51. J. Heyd, G. E. Scuseria,  and M. Ernzerhof, J. Chem. Phys. 124, 219906 (2006).
  52. G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  53. G. Kresse and J. Hafner, Phys. Rev. B 49, 14251 (1994).
  54. G. Kresse and J. Furthmüller, Comput. Mat. Sci. 6, 15 (1996a).
  55. G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
  56. M. Grüning, A. Marini,  and A. Rubio, J. Chem. Phys. 124, 154108 (2006).
  57. D. A. Kitchaev, H. Peng, Y. Liu, J. Sun, J. P. Perdew,  and G. Ceder, Phys. Rev. B 93, 045132 (2016).
  58. R. Peverati and D. G. Truhlar, J. Chem. Phys. 136, 134704 (2012).
  59. J. F. Janak, Phys. Rev. B 18, 7165 (1978).
  60. P. Mori-Sánchez, A. J. Cohen,  and W. Yang, Phys. Rev. Lett. 100, 146401 (2008).
  61. W. Yang, A. J. Cohen,  and P. Mori-Sánchez, J. Chem. Phys. 136, 204111 (2012).
  62. R. W. Godby, M. Schlüter,  and L. J. Sham, Phys. Rev. Lett. 56, 2415 (1986).
  63. E. R. Johnson, A. D. Becke, C. D. Sherrill,  and G. A. DiLabio, J. Chem. Phys. 131, 034111 (2009).
  64. E. R. Johnson, R. A. Wolkow,  and G. A. DiLabio, Chem. Phys. Lett. 394, 334 (2004).
  65. E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description