Importance of the Kinetic Energy Density for Band Gap Calculations in Solids with Density Functional Theory.

Recently, exchange-correlation potentials in density functional theory were developed with the goal of providing improved band gaps in solids. Among them, the semilocal potentials are particularly interesting for large systems since they lead to calculations that are much faster than with hybrid functionals or methods like GW. We present an exhaustive comparison of semilocal exchange-correlation potentials for band gap calculations on a large test set of solids, and particular attention is paid to the potential HLE16 proposed by Verma and Truhlar. It is shown that the most accurate potential is the modified Becke-Johnson potential, which, most noticeably, is much more accurate than all other semilocal potentials for strongly correlated systems. This can be attributed to its additional dependence on the kinetic energy density. It is also shown that the modified Becke-Johnson potential is at least as accurate as the hybrid functionals and more reliable for solids with large band gaps.

Introduction

The calculation of the electronic properties of molecules, surfaces, and bulk solids is done mostly with the Kohn–Sham (KS) scheme of density functional theory (DFT),[1,2] which is considered as a fast method especially if the exchange and correlation effects are approximated at the semilocal level, i.e., by using functionals Exc of the local density approximation (LDA), generalized gradient approximation (GGA), or meta-GGA (MGGA).[3] Of high interest is the electronic gap Δg, which is defined as I – A, where I and A are the ionization potential and electron affinity, respectively. However, the calculation of Δg with the KS-DFT orbital energies is not obvious both from the theoretical and practical points of view. KS-DFT is in principle a ground-state method, which raises the question how well-founded is the use of KS-DFT for the calculation of Δg, which is an excited-state property. Furthermore, it has been shown that the KS band gap (defined as the conduction band minimum minus the valence band maximum) calculated with the exact (but unknown) potential vxc = δExc/δρ differs from the band gap Δg by the derivative discontinuity Δxc,[4−7] which can be of the same order of magnitude as the band gap Δg itself.[8] In solid-state physics, the Green-functions-based GW methods provide a formal way of calculating the band gap Δg from the quasi-particles band structure;[9−11] however, these methods lead to very expensive calculations and can not be applied routinely to large systems. In addition, the GW method is usually applied in a perturbative way, i.e., as a first-order correction to the orbital energies obtained from a self-consistent KS-DFT calculation, such that the results may depend strongly on the exchange-correlation functional used in the KS-DFT calculation.

Therefore, it has been very common to use KS-DFT for band gap calculations despite it may not be a fully justified procedure depending on the type of potential that is used.[12] In general, the practical problem with KS-DFT is to choose the right functional Exc (or potential vxc) for the problem at hand. It is well-known that LDA and the standard GGA functionals like PBE[13] or BLYP[14,15] severely underestimate the band gap (see, e.g., refs (16) and (17)). Major improvements in band gap prediction have been obtained by using the hybrid functionals,[18−20] which lead to an accuracy similar to the GW method.[21] Therefore, these last ten years have seen a steady increase in the number of solid-state applications using hybrid functionals, in particular the one from Heyd, Scuseria, and Ernzerhof (HSE06).[22,23] Note that hybrid functionals are usually not implemented in the strict KS framework, but in the generalized KS (gKS) framework,[24] thus leading to a potential that is nonmultiplicative since a fraction of Hartree–Fock exchange is used.[25] With such nonmultiplicative potentials, (a part of) the discontinuity Δxc is included in the orbital energies.[6] However, hybrid functionals are much more expensive than semilocal functionals (albeit less than GW calculations), such that they can not be applied routinely to large systems, in particular with codes based on plane-waves basis functions. Thus, the search for a fast semilocal and reliable DFT method for electronic structure calculation, more particularly for band gaps, is of very high interest.[26−30]

The goal of this work is to present a comparison between various DFT methods for band gap prediction. A particular focus will be on the HLE16 (high local exchange) GGA proposed by Verma and Truhlar,[17] which was shown to be of similar accuracy as the hybrid functional HSE06. However, what is clearly missing in their study is a comparison with the other semilocal potentials that were previously shown to be also very accurate for band gaps in solids (see, e.g., refs (26) and (31−36) for the modified Becke–Johnson potential[26] and refs (28), (37), and (38) for the potential of Armiento and Kümmel[28]). In order to provide such a comparison, a large test set of 76 solids will be considered. For comparison purposes, results obtained with hybrid functionals will also be shown.

Methods and Computational Details

The semilocal methods that improve over the standards like PBE for the band gap can be divided into two groups according to their type of potential operator v̂xc, multiplicative (KS framework) or nonmultiplicative (gKS framework), and the focus of the present work is to make a thorough comparative assessment of multiplicative potentials for band gaps in solids. The potentials that are considered are the LDA and a modification of it, four GGAs and one MGGA. All these LDA and GGA potentials are obtained as functional derivative vxc = δExc/δρ of energy functionals. The LDA consists of the exchange[2] and correlation[39] of the homogeneous electron gas, while Sloc (acronym for local Slater potential[30]) consists of a slightly modified but strongly enhanced exchange LDA (vxSloc = −1.67ρ0.3 compared to vxLDA ≃ −0.7386ρ1/3) with no correlation added. The GGAs are the exchange-correlation PBE from Perdew et al.,[13] the exchange of Engel and Vosko[40] (EV93PW91, combined with correlation from Perdew and Wang[41] as done in ref (42)), the exchange from Armiento and Kümmel[28,37,38] (AK13, no correlation added as done in refs (28) and (37)), and the recently proposed HLE16[17] that consists of a modification of the HCTH/407 functional[43] (the exchange and correlation components are multiplied by 1.25 and 0.5, respectively). The tested MGGA is the modified Becke–Johnson exchange potential combined with LDA for correlation (mBJLDA),[26] which is based on the exchange potential of Becke and Johnson.[44]

The test set (see Table S1 of the Supporting Information) consists of most of the solids considered in refs (17), (21), and (45) and includes a large variety of solids: rare gases, sp-semiconductors, ionic insulators, and transition-metal compounds. All systems, except the antiferromagnetic Cr2O3, Fe2O3, MnO, FeO, CoO, and NiO, were treated as nonmagnetic. Antiferromagnetic 3d transition metal oxides are the typical systems with strongly correlated 3d electrons, which are particularly difficult cases for LDA and GGA methods[46] since the self-interaction error is much more important than in the case of itinerant metals like Fe or Ni with more delocalized 3d electrons. VO2 is also a difficult case for which it was shown in ref (47) that mBJLDA leads to qualitatively correct trends. Also included in the test set are Cu1+ compounds, e.g., Cu2O, which is known to be problematic for the mBJLDA potential.[32] For comparison purposes, we will also show the results obtained with the hybrid functionals B3PW91[25] and HSE06[22,23] that have been tested very recently by Crowley et al.[21] and Garza and Scuseria,[45] respectively. However, note that the HSE06 results of the present work were actually obtained with YS-PBE0[48] (YS stands for Yukawa screened), which leads to band structures that are quasi-identical to HSE06 thanks to the similarity between the error-function- and exponential-based screened Coulomb operators.[49]

The calculations of the electronic structure were done with the all-electron code WIEN2k,[60] which is based on the linearized-augmented plane-wave method[61] that allows for a very accurate solution of the KS-DFT equations. The experimental geometry (see Table S1) was used for all calculations such that the differences in the band gaps obtained with the various methods are only due to the exchange-correlation potential. This is the most reasonable choice, in particular since the energy functionals EV93PW91, AK13, and HLE16 were shown to be very inaccurate for geometry optimization,[17,38,62] while mBJLDA is only a potential with no corresponding energy functional.[63,64] Actually, it was shown in ref (17) that the HLE16 band gaps calculated at the corresponding HLE16 equilibrium geometry (severely underestimated in many cases) are much larger (in some sense artificially) than the band gaps evaluated at a more reasonable geometry. The calculation parameters like the size of the basis set or the number of k-points for Brillouin zone integrations were chosen such that the band gaps, including those calculated with the hybrid functionals, should be converged with an error bar below 0.05 eV. We note that a comparison of our HSE06 (YS-PBE0) and B3PW91 band gaps with those obtained from codes using Gaussian basis functions[21,45] shows in general a rather fair agreement; however, large discrepancies in the order of 1–2 eV were obtained in a few cases, e.g., MnO, NiO,[45] and systems with very large band gaps like LiF, NaCl, or SiO2.[21] The agreement of our results with the HSE results from refs (53) and (65), which were obtained with the projector augmented-wave method, is in general clearly better, thus pointing to some problems with the Gaussian basis sets.

Results

The results for the fundamental band gap are shown in Table and graphically in Figure a–f (Figures S1–S9 of the Supporting Information show the results individually for each method), while the summary statistics is shown in Table , where M(R)E, MA(R)E, and STD(R)E denote the mean (relative) error, the mean absolute (relative) error, and the standard deviation of the (relative) error, respectively. The largest errors are obtained with LDA and PBE, which systematically underestimate the band gap and lead to ME of about −2 eV, while in terms of relative error the underestimation is above 50%. The band gaps obtained with EV93PW91 are slightly closer to the experimental values since the MAE and MARE are 1.55 eV and 36%, which, however, can still be considered as rather large with a clear trend for underestimation. The recently proposed AK13, Sloc, and HLE16 lead to further (and substantial) reductions of the error in the band gap with MAE in the range 0.75–0.90 eV and MARE of 30% for Sloc and ∼25% for AK13 and HLE16. However, the best agreement with experiment is obtained with the MGGA mBJLDA since the MAE and MARE drop to 0.47 eV and 15%, respectively. From the ME and MRE, the tendency of all semilocal potentials is to underestimate the band gap; however, it is rather small for the AK13 and mBJLDA potentials. The STD(R)E (see Table ) is a measure of how scattered are the (relative) errors that are obtained with a potential, and the smaller STD(R)E is, the more a (relative) error is predictable and likely to be close to the M(R)E. The results show that among all semilocal potentials, the mBJLDA potential leads to the smallest values for both the STDE and STDRE. For the STDE, the other methods lead to values that are clearly larger including HLE16 with a STDE that is about twice larger (0.57 eV for mBJLDA and 1.07 eV for HLE16). It is also interesting to note that despite Sloc is a very simple potential,[30] which depends only on the electron density ρ (no dependence on derivatives of ρ), it is of the same accuracy as HLE16, which has a much more complicated analytical form that involves ρ and its two first derivatives since the energy functional is a GGA.[17]

Table 1

Calculated and Experimental[21,50−59] Fundamental Band Gaps (in eV). The Space Group Number Is Indicated in Parenthesis.

solid	LDA	PBE	EV93PW91	AK13	Sloc	HLE16	mBJLDA	HSE06	B3PW91	exptl
Ne (225)	11.42	11.58	11.21	19.98	18.30	15.06	22.33	14.27	15.09	21.48
Ar (225)	8.18	8.70	9.26	15.13	12.57	11.74	13.84	10.37	10.89	14.15
Kr (225)	6.76	7.26	7.95	12.82	10.61	10.09	10.80	8.71	9.17	11.59
Xe (225)	5.78	6.24	7.03	10.67	9.02	8.80	8.48	7.44	7.86	9.29
C (227)	4.10	4.14	4.31	4.78	5.18	4.55	4.92	5.26	5.71	5.50
Si (227)	0.47	0.57	0.91	1.60	1.70	1.38	1.15	1.17	1.59	1.17
Ge (227)	0.00	0.06	0.58	0.70	0.00	0.12	0.83	0.82	1.08	0.74
Al2O3 (167)	6.19	6.26	6.67	7.92	7.30	7.20	8.34	8.08	8.43	8.80
SiC (216)	1.32	1.36	1.52	2.18	2.83	2.06	2.25	2.23	2.69	2.42
SiO2 (α-quartz, 152)	5.70	5.93	6.47	8.17	7.58	7.53	8.70	7.77	8.12	9.65
SiO2 (β-cristobalite, 227)	5.54	5.79	5.81	9.48	9.20	8.10	10.43	7.39	8.05	8.90
BN (216)	4.35	4.47	4.82	5.67	6.18	5.38	5.80	5.76	6.20	6.36
BP (216)	1.19	1.24	1.46	2.02	2.37	1.85	1.85	1.98	2.44	2.10
BAs (216)	1.14	1.20	1.46	1.96	2.23	1.82	1.71	1.86	2.30	1.46
AlN (216)	3.24	3.33	3.58	4.44	5.25	4.41	4.88	4.55	4.97	4.90
AlN (wurtzite, 186)	4.18	4.16	4.43	5.31	5.17	4.80	5.51	5.49	5.89	6.19
AlP (216)	1.45	1.59	2.10	3.00	3.08	2.75	2.31	2.30	2.72	2.50
AlAs (216)	1.35	1.45	1.90	2.76	1.95	2.55	2.13	2.11	2.52	2.23
AlSb (216)	1.15	1.22	1.56	2.29	1.16	1.84	1.75	1.80	2.19	1.69
GaN (216)	1.67	1.66	1.95	2.37	2.57	2.39	2.85	2.85	3.20	3.28
GaN (wurtzite, 186)	1.94	1.94	2.27	2.71	2.82	2.68	3.17	3.15	3.52	3.50
GaP (216)	1.44	1.60	2.08	2.60	1.57	2.35	2.25	2.28	2.68	2.35
GaSb (216)	0.00	0.11	0.60	0.76	0.00	0.26	0.95	0.88	1.11	0.82
GaAs (216)	0.30	0.54	1.13	1.45	0.12	0.94	1.64	1.40	1.64	1.52
InN (wurtzite, 186)	0.02	0.03	0.18	0.56	0.73	0.60	0.89	0.70	1.03	0.72
InP (216)	0.46	0.68	1.30	1.81	0.63	1.29	1.62	1.43	1.72	1.42
InAs (216)	0.00	0.00	0.36	0.73	0.00	0.23	0.67	0.45	0.69	0.42
InSb (216)	0.00	0.00	0.30	0.51	0.00	0.00	0.47	0.45	0.68	0.24
SnO2 (136)	1.09	1.24	1.64	2.20	1.40	1.62	3.19	2.88	3.17	3.59
SnSe (62)	0.48	0.55	0.76	1.10	1.17	1.03	0.89	0.98	1.31	0.90
SnTe (225)	0.09	0.07	0.07	0.22	0.74	0.47	0.15	0.17	0.07	0.36
Sb2Te3 (166)	0.04	0.00	0.20	0.18	0.35	0.34	0.24	0.31	0.13	0.28
LiH (225)	2.67	3.03	3.72	6.19	4.33	4.21	5.06	4.06	4.50	4.94
LiF (225)	8.95	9.20	9.98	12.55	11.76	11.84	12.89	11.46	11.74	14.20
LiCl (225)	6.06	6.41	7.42	9.78	8.08	8.56	8.64	7.81	8.15	9.40
NaF (225)	6.05	6.41	7.17	11.10	9.36	9.44	11.46	8.57	8.89	11.50
NaCl (225)	4.79	5.21	6.18	9.82	7.16	7.73	8.45	6.61	6.96	8.50
KF (225)	5.79	6.17	6.82	10.98	9.19	8.91	10.40	8.18	8.56	10.90
KCl (225)	4.75	5.19	6.02	9.79	7.38	7.59	8.48	6.53	6.93	8.50
BeO (wurtzite, 186)	7.58	7.65	8.22	9.42	9.38	8.91	9.66	9.48	9.82	10.60
MgO (225)	4.69	4.78	5.17	6.69	6.11	5.91	7.13	6.47	6.83	7.83
MgS (216)	3.29	3.57	4.50	6.21	4.76	5.26	5.17	4.66	5.00	4.78
MgSe (225)	1.71	1.87	2.21	3.65	2.66	2.86	2.93	2.74	3.16	2.47
MgTe (216)	2.27	2.51	3.22	4.74	3.19	3.68	3.61	3.39	3.73	3.60
CaO (225)	3.49	3.67	4.07	5.03	4.36	4.63	5.35	5.26	5.59	7.00
CaF2 (225)	6.92	7.30	8.11	9.81	9.04	9.26	10.34	9.37	9.72	11.80
BaS (225)	2.00	2.20	2.79	4.03	3.43	3.56	3.27	3.11	3.50	3.88
BaSe (225)	1.82	1.97	2.50	3.68	3.17	3.25	2.86	2.79	3.17	3.58
BaTe (225)	1.48	1.61	2.09	3.23	2.72	2.79	2.28	2.31	2.70	3.08
ScN (225)	0.00	0.00	0.18	0.70	0.41	0.40	0.88	0.90	1.26	0.90
TiO2 (rutile, 136)	1.80	1.89	2.03	2.23	1.44	1.75	2.56	3.34	3.66	3.30
TiO2 (anatase, 141)	2.00	2.11	2.30	2.54	1.74	2.06	2.92	3.57	3.87	3.40
SrTiO3 (221)	1.78	1.88	2.06	2.34	1.61	1.87	2.68	3.29	3.62	3.30
VO2 (M1, 14)	0.00	0.00	0.00	0.00	0.00	0.00	0.51	1.03	1.35	0.60
Cr2O3 (167, 146)	1.20	1.63	1.88	2.35	1.46	1.53	3.68	4.42	4.54	3.40
Fe2O3 (167, 146)	0.33	0.56	0.88	1.50	1.49	1.70	2.35	3.24	3.35	2.20
MnO (225, 166)	0.74	0.86	1.23	2.56	3.71	3.26	2.94	2.85	3.19	3.90
FeO (225, 166)	0.00	0.00	0.10	0.84	0.20	0.13	1.84	2.35	2.08	2.40
CoO (225, 166)	0.00	0.00	0.24	1.35	0.43	0.34	3.13	3.48	3.57	2.50
NiO (225, 166)	0.43	0.95	1.35	2.08	1.21	1.23	4.14	4.37	4.36	4.30
Cu2O (224)	0.53	0.53	0.57	0.84	1.27	0.81	0.81	1.98	2.34	2.17
CuSCN (160)	2.01	2.23	2.46	2.78	3.64	3.24	2.79	3.60	3.99	3.94
CuCl (216)	0.33	0.47	0.98	1.95	3.78	2.95	1.69	2.37	2.63	3.40
CuBr (216)	0.20	0.36	0.89	1.76	3.38	2.69	1.56	2.15	2.40	3.07
CuI (216)	0.95	1.12	1.66	2.36	3.32	2.98	2.20	2.65	2.89	3.12
ZnO (wurtzite, 186)	0.74	0.81	1.27	2.06	3.21	2.81	2.65	2.50	2.80	3.44
ZnS (216)	1.84	2.09	2.81	3.67	3.13	3.54	3.65	3.30	3.58	3.84
ZnSe (216)	1.02	1.27	1.97	2.67	2.04	2.50	2.75	2.37	2.63	2.82
ZnTe (216)	1.04	1.27	1.84	2.33	1.54	2.08	2.42	2.25	2.51	2.39
MoS2 (194)	0.79	0.86	1.02	1.25	1.18	1.21	1.08	1.41	1.81	1.29
AgCl (225)	0.62	0.92	1.74	2.82	3.39	3.31	2.95	2.41	2.68	3.25
AgBr (225)	0.38	0.67	1.48	2.44	2.76	2.77	2.50	2.01	2.27	2.71
AgI (216)	1.05	1.34	2.02	3.02	3.15	3.14	2.77	2.48	2.78	2.91
CdS (216)	0.88	1.16	1.91	2.84	2.21	2.60	2.67	2.14	2.44	2.50
CdSe (216)	0.36	0.63	1.36	2.14	1.46	1.87	1.99	1.52	1.79	1.85
CdTe (216)	0.51	0.76	1.39	1.98	1.12	1.61	1.79	1.57	1.84	1.61

Figure 1

Calculated versus experimental fundamental band gaps for the set of 76 solids. The right panels are zooms of the left panels focusing on band gaps smaller than 5 eV.

Table 2

Summary Statistics for the Error in the Calculated Band Gap for the Set of 76 Solids

	LDA	PBE	EV93PW91	AK13	Sloc	HLE16	mBJLDA	HSE06	B3PW91
ME (eV)	–2.17	–1.99	–1.55	–0.28	–0.76	–0.82	–0.30	–0.68	–0.36
MAE (eV)	2.17	1.99	1.55	0.75	0.90	0.90	0.47	0.82	0.73
STDE (eV)	1.63	1.56	1.55	0.89	0.93	1.07	0.57	1.21	1.14
MRE (%)	–58	–53	–35	–6	–21	–20	–5	–7	6
MARE (%)	58	53	36	24	30	25	15	17	23
STDRE (%)	23	23	23	31	37	28	22	22	35

Turning now to the comparison with hybrid functionals, Tables and 2 and Figure e,f also show the results obtained with B3PW91 and HSE06. The main conclusion is the same as above, namely, the mBJLDA potential leads to the smallest MAE and MARE; however, note that, in terms of MARE, mBJLDA and the screened hybrid HSE06 perform basically the same, while for the MAE, mBJLDA leads to a value that is nearly twice smaller since the MAE with HSE06 is 0.82 eV. These observations are on par with those obtained in ref (36) with a reduced test set of solids. The performances of HSE06 and the unscreened hybrid functional B3PW91 are rather similar, with B3PW91 leading to MAE and STDE that are slightly smaller, but MARE and STDRE that are larger as noticed in ref (45). In summary, the statistics for the absolute errors are the smallest for mBJLDA, while for the relative error, mBJLDA and HSE06 are of similar accuracy.

Besides the statistics, it is also instructive to discuss the results in more detail in order to figure out what are the outliers and the eventual trends among the classes of solids. This can be done most easily by inspecting Figures S1–S9, which show the results individually for each method. The most interesting observations are the following. As already stated above, LDA, PBE, and EV93PW91 underestimate basically all compounds, and for a few of them, e.g., Ge, InAs, and FeO, the system is described as metallic. With AK13, Sloc, and HLE16, the band gap is (slightly) overestimated in many cases, but there are still many outliers for which a large underestimation of the band gap is obtained, especially with Sloc and HLE16. Furthermore, some systems have still no band gap with these potentials (one with AK13, five with Sloc, and two with HLE16). The data points for mBJLDA are overall clearly closer to the diagonal line. The three most visible outliers are Cu compounds: Cu2O, CuCl, and CuBr for which the underestimation of the band gap is 1.36, 1.71, and 1.51 eV, respectively. Note that in these compounds, the oxidation state of the copper atom is Cu1+, while it is Cu2+ in CuO, where such an underestimation of the band gap with mBJLDA was not observed.[32]

For band gaps smaller than 5 eV, the HSE06 hybrid functional leads to results that are rather close to the experimental values, however, there is a clear underestimation for large band gaps like Ne, Ar, and LiF that is similar, although less pronounced, to what is observed with LDA, PBE, and EV93PW91. As indicated in the upper panels of Figures S1–S9, the slope of the linear regression for HSE06 is 0.73, which is close to 0.80 for HLE16, while LDA, PBE, and EV93PW91 lead to slopes in the range 0.62–0.65. In this respect, mBJLDA and AK13, with slopes of 0.97, are the most balanced potentials, not showing obvious underestimation or overestimation for small or large band gaps. For band gaps smaller than 2 eV, B3PW91 leads to the largest overestimations among all methods tested in this work, but overall the results are rather good for band gaps between 2 and 10 eV. At that point, it should be mentioned that hybrid functionals with a system-dependent fraction of Hartree–Fock exchange α have been shown to improve the band gap prediction with respect to the traditional hybrids. For instance, in refs (49), (54), and (55), the use of the dielectric function to determine α revealed to be very useful to cure the underestimation for large band gaps, while the same was concluded in ref (66) by using a density-based quantity for α similar to that for mBJLDA.

As mentioned above, antiferromagnetic transition-metal oxides are difficult cases for the usual LDA/PBE methods, which lead to qualitatively wrong results; however, these systems can be treated more accurately by LDA + U or hybrid functionals (see, e.g., refs (52), (65), (67), and (68)). Among the local potentials tested here, only mBJLDA is able to provide band gaps that are in qualitative agreement with experiment. For FeO, CoO, and NiO, the underestimation with all other semilocal potentials is at least 2 eV. Similarly, only mBJLDA and the hybrids open a band gap for nonmagnetic VO2 (see ref (47) for a much more detailed discussion).

Exchange-correlation potentials vxc in Si (diamond structure) and Kr (face-centered cubic) are shown in Figure , where significant differences between the potentials can be observed. The most striking features concern HLE16 which shows extremely large oscillations in the regions between the atoms. The peaks in the HLE16 potential, which can be of positive or negative values, are obviously a direct consequence of the energy functional HLE16, which is based on the highly parametrized HCTH/407 functional. The enhancement factor of such functionals are usually less smooth than the more simple ones like PBE, leading to strong oscillations in the potential when taking the functional derivative. As already discussed in refs (69) and (70), AK13 and, to a lesser extent, EV93 show strong positive values in the interstitial regions, which may be responsible for some of the overestimations of the band gap. Sloc (as well as LDA) is the most structureless potential; nevertheless, as already mentioned above, it improves over LDA/PBE the same way as HLE16 does. Note that also the mBJLDA potential is smooth in the interstitial region.

Figure 2

Potentials vxc plotted (a) in Si from the atom at (1/8, 1/8, 1/8) (d = 0) to the middle of the unit cell and (b) in Kr from the atom at (0, 0, 0) (d = 0) to the middle of the unit cell.

In various previous works (see, e.g., refs (70) and (71)), it was shown that usually the more the magnitudes of the potential around the atoms and in the interstitial differ, the larger is the band gap. This is the case when the orbitals of the valence band maximum (VBM) and conduction band minimum (CBM) are mainly located around an atom or in the interstitial, respectively, like in sp-semiconductors and ionic insulators. This explains why the band gap is larger (with respect to LDA/PBE) with basically all semilocal potentials tested in this work, including the very simple LDA-like Sloc potential (see Figure ). However, in systems like transition metal oxides with a d–d band gap, the VBM and CBM may be located at the same position in space (around the transition metal atom), such that the simple mechanism explained above to open the band gap does not work. This is why the band gap is still by far too small with all LDA/GGA-type potentials, including AK13, Sloc, and HLE16, for Cr2O3, Fe2O3, FeO, CoO, NiO, and VO2 as noticed above. mBJLDA depends on an additional ingredient, the kinetic energy density,[72] which provides a more efficient way for the potential to act differently on orbitals localized on the same atom according to their angular distribution (e.g., t2 versus e). This feature of the BJ-based potentials, which has been discussed in great detail in our previous works,[32,69,70,73] is crucial in explaining the opening of the band gap in strongly correlated systems. However, it should be also mentioned that it does not work in all cases, like for the Cu1+ compounds (Cu2O, CuCl, and CuBr, in particular), where the Cu-sd hybridization is essential and for which hybrid functionals are more accurate. Note that, for these systems, Sloc is the only semilocal potential leading to values close to experiment.

Summary and Conclusion

In summary, several semilocal multiplicative potentials have been compared for band gap calculations. Using a test set of 76 solids of various types, it has been shown that on average mBJLDA is the most accurate potential and clearly outperforms all others including HLE16. Except for very few outliers (the Cu1+ compounds), the mBJLDA potential is equally accurate for small and large band gaps as well as for strongly correlated systems for which all other potentials completely fail. Furthermore, mBJLDA is as accurate as the hybrid functionals for semiconductors, but more reliable for large band gap insulators. Our results strongly support the use of the kinetic energy density in the design of KS potentials for band gap calculations or to approximate the exact KS potential as done in our previous works.[69,73] To finish, we should also mention that nonmultiplicative MGGA potentials, which include the derivative discontinuity Δxc in the energy spectrum, are extremely promising for band gap calculations; however, at the moment no such potential leading to satisfying results has been proposed.[12] Therefore, at the moment the mBJLDA potential represents the best alternative to the much more expensive hybrid or GW methods.