One of the well-known limitations of Kohn–Sham density functional theory is the tendency to strongly underestimate bandgaps. Meta-generalized gradient approximations (mGGAs), which include the kinetic energy density in the functional form, have been shown to significantly alleviate this deficiency. In this study, we explore the mechanisms responsible for this improvement from the angle of the underlying local densities. We find that the highest occupied and lowest unoccupied states are distinct in the space of the underlying descriptors. The gap opening is compared to a simple scaling of the local density approximation, and two mechanisms responsible for opening the mGGA gaps are identified. First of all, the relatively large negative derivative of the functional form with respect to reduced kinetic energy tends to elevate the lowest unoccupied state. Second, the curvature of functional, which ensures that it is bounded, tends to lower the highest occupied state. Remarkably, these two mechanisms are found to be transferable over a large and diverse database of compounds.

## I. INTRODUCTION

Kohn–Sham density functional theory (KS-DFT)^{1} has become the workhorse of electronic structure calculations due to its favorable trade-off between accuracy and computational cost. However, the accuracy of KS-DFT depends on approximate exchange–correlation (xc) functionals that have a number of known problems. Probably the best-known weakness of KS-DFT is the systematic underestimation of bandgaps, which results in narrow-band-gap semiconductors being predicted to be metals and the gap of insulators being underestimated by approximately a factor of two.^{2,3} The problems of bandgap prediction are usually attributed to the derivative discontinuity, which is absent in KS-DFT but must be a part of the exact functional.^{4,5} This is corroborated by the fact that functionals that belong to the generalized KS (gKS) scheme^{6} usually result in more realistic gaps than those obtained with KS-DFT.^{3} Even meta-generalized gradient approximations (mGGAs), which include a dependency on the kinetic energy density in the functional form and represent probably the simplest gKS functionals, result in substantial improvements in bandgap predictions over KS functionals.^{3}

Approximate xc functionals will typically involve a trade-off between the accuracy that can be obtained for different properties. Recently, we systematically explored these trade-offs by optimizing 25 different mGGA exchange functionals to reproduce lattice constants, cohesive energies, and bandgaps to different levels of accuracy.^{7} While a trade-off in the accuracy was observed for the ground-state properties and the bandgap, these functionals still managed to strike a balance.^{7} It was thus possible to optimize mGGAs that perform well for the ground-state properties while simultaneously systematically improving bandgap predictions. Additionally, it was possible to construct a mGGA, which yields bandgap predictions with an accuracy similar to hybrid functionals. This corroborates the results obtained for the TASK functional, which was designed with a focus on the non-locality^{8} and also performs very well for bandgap predictions.

Interestingly, the optimized functional forms were also similar to the designed forms.^{7} Specifically, we found that functionals where a large weight was put on minimizing the bandgap error tended to have an enhancement factor with a large negative gradient with respect to the reduced kinetic energy. This data driven result is in agreement with the design idea of the TASK^{8} functional. This provides a first idea of how mGGAs are able to improve the bandgap prediction. However, a more direct understanding in terms of the underlying densities is somewhat lacking. This approach has proven useful in understanding how different xc functionals impact ground state properties.^{9–11} It is worth noting that the overall usefulness of a negative gradient of the functional with respect to the reduced kinetic energy was first recognized^{10} in an analysis of an early trained functional.^{12}

^{13–15}However, the gKS implementation of the mGGAs relies on a non-multiplicative xc potential operator, which makes it non-trivial to directly analyze the xc contribution to the bandgap. In this paper, we aim at a more direct understanding of the influence of the functional form on the bandgap. Take the definition of the fundamental bandgap as the difference between its ionization energy (

*I*) and electron affinity (

*A*),

*E*(

*N*) corresponds to the total ground-state energy of the system with

*N*electrons. This expression can be evaluated directly using the densities obtained by removing one electron from highest occupied (HO) state or adding one to the lowest unoccupied (LU) state, respectively.

^{16}Compared to simple (semi-)local KS density functionals, advanced functionals should be better at discriminating regions typically occupied by the HO or LU states. The mGGA energy functional relies only on the local kinetic and electron densities, and Eq. (1) thus makes it possible to obtain a direct local understanding of the factors influencing the mGGA gaps.

In this work, we aim to explore the mGGAs bandgaps from the angle of orbital shapes and the underlying local densities. We use Eq. (1) to directly evaluate the difference of the ionization energies and electron affinities. The qualitatively different behavior of some materials with the local density approximation (LDA) exchange functional suggests that the exact shape of the HO and LU orbitals, which can strongly vary depending on the specific materials, plays an important role in the bandgap prediction. We investigate how the removal and addition of an electron to the HO and LU states influences the underlying densities used in mGGA functionals. We identify two mechanisms directly related to the functional form and its derivatives and show how these mechanisms are surprisingly transferable over a wide range of materials.

## II. METHODOLOGY

### A. mGGA gaps

*n*,

*F*is the specific functional form. The

*p*and

*t*parameters are the reduced gradient of the density,

*τ*

^{KS}= (1/2)∑|∇

*ψ*

_{i}|

^{2}is the KS kinetic energy density and denominator is the kinetic energy density of the homogeneous electron gas.

^{16}We thus define

*N*

_{k}is the number of

*k*points used for the specific material and the

*m*, 0, and

*p*subscripts correspond to the one-electron-removed, the neutral, and the one-added-electron cases. We note here that

*n*

_{HO}and

*τ*

_{HO}defined like this contain the negative sign. The gap is then given as

*E*

_{x}will give the corresponding exchange gap,

*G*

_{x}. The corresponding normalized value

*p*

_{HO}can be evaluated from

*n*

_{0}+

*n*

_{HO}and

*t*

_{HO}from

*n*

_{0}+

*n*

_{HO}and

*τ*

_{0}+

*τ*

_{HO}and similar for the LU.

In this paper, we will focus on the KTBMgap functional, but all conclusions are quite general to mGGAs of the form given by Eq. (2). The KTBMgap functional was part of a group of functionals constructed to explore trade-offs in accuracy that can be obtained with a functional form given by Eq. (2). Illustrations of some of the functional forms are given in the inset in Fig. 1. The functional forms for functionals optimized toward a low mean absolute relative error (MARE) of the cohesive energy (KTBM5), the lattice constant (KTBM9), and the bandgap are shown (KTBM23). Furthermore, a general purpose functional (KTBM13) and a functional (KTBMgap) optimized toward the mean absolute error (MAE) of the bandgap is illustrated. The functional form of KTBM13 agrees qualitatively with a designed general-purpose mGGA, such as SCAN,^{19} and the KTBMgap form with TASK, which was designed with a focus on the non-locality.^{8} As illustrated in Table I, the results obtained for the KTBM13 and KTBMgap also agree well with the corresponding designed functionals. The main plot of Fig. 1 furthermore illustrates the derivative with respect to *α* = *t* − 5/3*p*. Comparing KTBM23 and KTBMgap, it is clear that the priority of the MAE over the MARE, which in effect prioritizes the large gaps, leads to a much stronger dependence on *α*. Furthermore, *F* must be bounded from both below and above for any “reasonable” functional. This means that the derivative tends to zero for larger *α* values.

. | MAE(a_{0}) (Å)
. | MAE(E_{coh}) (eV/atom)
. | MAE(E_{g}) (eV)
. |
---|---|---|---|

LDA | 0.07 | 0.77 | 1.18 |

SCAN | 0.03 | 0.19 | 0.62 |

TASK | 0.23 | 1.34 | 0.47 |

KTBM13 | 0.03 | 0.19 | 0.62 |

KTBMgap | 0.06 | 0.57 | 0.46 |

. | MAE(a_{0}) (Å)
. | MAE(E_{coh}) (eV/atom)
. | MAE(E_{g}) (eV)
. |
---|---|---|---|

LDA | 0.07 | 0.77 | 1.18 |

SCAN | 0.03 | 0.19 | 0.62 |

TASK | 0.23 | 1.34 | 0.47 |

KTBM13 | 0.03 | 0.19 | 0.62 |

KTBMgap | 0.06 | 0.57 | 0.46 |

### B. Database

Our database contains 288 materials and is a subset of the database used in Borlido *et al.*’s large-scale study of xc functionals for bandgaps.^{3} While self-consistent mGGA calculations are now becoming standard,^{20–23} they are obviously not practical when optimizing a functional. We therefore performed our evaluation of Eq. (1) non-self-consistently using the densities obtained with the RPBE functional^{24} as it has been shown to result in the lowest non-self-consistency error for mGGA functionals.^{16} However, since RPBE predicts no gap for some compounds in the original database, which prevents the correct identification of the HO and LU states, we removed any material with an RPBE bandgap of less than 0.1 eV. Furthermore, we also removed any compound where the difference between the self-consistent bandgap and the non-self-consistent gap calculations resulted in a difference greater than 0.25 eV for the TASK functional.

## III. RESULTS AND DISCUSSION

### A. The LDA gap

The LDA functional results in 1.18 eV mean absolute error (MAE), Table I, and 50% mean absolute relative error (MARE) for bandgaps on our 288 material database with all gaps being underestimated. It is, in principle, relatively straightforward to obtain better average bandgaps even with the very limited LDA exchange functional. Simply scaling up the exchange energy would thus lead to better average gaps. Scaling the LDA exchange energy density, Eq. (3), up by 21%, corresponding to a constant *F* = 1.21 in Eq. (2), results in the MAE for the bandgap dropping to 0.67 eV and the MARE to 35%. Further increasing *F* to a constant 1.28 results in a MAE of 0.62 eV and a MARE of 37%. Similar ideas of simply scaling functionals to obtain better gaps have already been employed on multiple levels of density functional approximations.^{25–27}

Two compounds found on the opposite sides of the dashed line in Fig. 2 are isostructural Si and Ge, which we will analyze in more detail. The $\epsilon \u0304xLDA$ values are given in the first column of Table II. The corresponding LDA exchange gaps, Eq. (10), have opposite signs, $GxLDA=2.87$ eV for Si and $GxLDA=\u22122.85$ eV for Ge. In Fig. 3(a), the cumulative HO and LU densities of Si and Ge are shown. The curves are obtained by integrating |*n*_{HO}| and *n*_{LU} in radial concentric spheres in the Voronoi cell of the diamond structure of the two compounds. The curves confirm how the Si HO-state is more localized than the Si LU-state as is expected for the majority of the database. Ge, on the other hand, has a more localized LU state in accordance with its negative $GxLDA$. Also shown in Fig. 3(a) are approximate region attributed, through a simple derivative criteria, to the HO and LU states. These regions can be compared to the radial distribution functions of the atomic states of Si and Ge shown in Fig. 3(b). The comparison shows how the anomalously localized LU-state in germanium can be attributed to the well-known contribution from 3*d*-states.

KTBMgap . | TASK . | KTBMgap . | TASK . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|---|

$\epsilon \u0304xLDA$ (eV) . | $\Delta Gxn$ (eV) . | F^{eff}
. | $\Delta Gxn$ (eV) . | F^{eff}
. | $\epsilon \u0304xpt$ (eV) . | $\Delta Gxpt$ (eV) . | dF^{eff}
. | $\Delta Gxpt$ (eV) . | dF^{eff}
. | ||

Si | HO | 11.87 | 0.77 | 1.01 | 0.72 | 0.99 | −9.48 | 0.19 | −0.002 | −0.28 | 0.016 |

LU | −9.00 | 0.93 | 0.91 | 21.53 | 0.008 | −0.006 | |||||

Ge | HO | 13.11 | 0.80 | 0.97 | 0.89 | 0.95 | −18.73 | 0.33 | −0.039 | 0.03 | −0.035 |

LU | −15.96 | 0.92 | 0.90 | 32.85 | −0.012 | −0.019 | |||||

NaF | HO | 21.70 | 2.32 | 1.05 | 2.20 | 1.02 | −26.67 | 2.15 | −0.051 | 2.65 | −0.034 |

LU | −9.46 | 0.87 | 0.82 | 70.76 | 0.011 | 0.025 |

KTBMgap . | TASK . | KTBMgap . | TASK . | ||||||||
---|---|---|---|---|---|---|---|---|---|---|---|

$\epsilon \u0304xLDA$ (eV) . | $\Delta Gxn$ (eV) . | F^{eff}
. | $\Delta Gxn$ (eV) . | F^{eff}
. | $\epsilon \u0304xpt$ (eV) . | $\Delta Gxpt$ (eV) . | dF^{eff}
. | $\Delta Gxpt$ (eV) . | dF^{eff}
. | ||

Si | HO | 11.87 | 0.77 | 1.01 | 0.72 | 0.99 | −9.48 | 0.19 | −0.002 | −0.28 | 0.016 |

LU | −9.00 | 0.93 | 0.91 | 21.53 | 0.008 | −0.006 | |||||

Ge | HO | 13.11 | 0.80 | 0.97 | 0.89 | 0.95 | −18.73 | 0.33 | −0.039 | 0.03 | −0.035 |

LU | −15.96 | 0.92 | 0.90 | 32.85 | −0.012 | −0.019 | |||||

NaF | HO | 21.70 | 2.32 | 1.05 | 2.20 | 1.02 | −26.67 | 2.15 | −0.051 | 2.65 | −0.034 |

LU | −9.46 | 0.87 | 0.82 | 70.76 | 0.011 | 0.025 |

### B. Origin of the mGGA gap

Clearly, the mGGAs have more flexibility than can be obtained by a simple constant scaling of a purely local functional. As expected, the mGGAs also reduce the gap MAEs, Table I, more than the “enhanced LDAs” discussed in Sec. III A. Importantly, the mGGAs increase the predicted gap for all compounds in the database.

*F*(

*p*

_{0}+

*p*

_{HO},

*t*

_{0}+

*t*

_{HO}) as

*F*(

*p*

_{0},

*t*

_{0}) in Eq. (13) is a distribution, not a simple constant. The mGGAs are able to exploit this to systematically improve the bandgap predictions. While the distribution

*F*(

*p*

_{0},

*t*

_{0}) is the same for the HO and LU in Eq. (13), it is weighted different because $\Delta ex,HOLDA$ is different from $\Delta ex,LULDA$. Figure 5 illustrates $\Delta exLDA$’s role as a weighting factor and reveals how mGGAs functionals are able to have different enhancement factors for the HO and LU electrons. Visually the LU state gives the impression that it is more spread out in

*p*−

*t*space than the HO state. Interestingly, this seems to be the case for both Si and Ge despite the densities of especially the LU states being different for the two compounds. Since the electrons are spread out in relatively large regions in

*p*−

*t*space, we define the weighted averages,

*p*–

*t*values and the similarity of Si and Ge in

*p*–

*t*space. Based on these average

*p*and

*t*values, we can make a rough approximation of the mGGA exchange gap as $Gxn\u2248\epsilon \u0304x,HOLDAF(p\u0304HO,t\u0304HO)+\epsilon \u0304x,LULDAF(p\u0304LU,t\u0304LU)$. As $\epsilon \u0304xLDA$ is negative, a $F(p\u0304LU,t\u0304LU)<F(p\u0304HO,t\u0304HO)$ will increase the gap. The decreasing

*F*with

*t*thus causes an increase in the gap. This increase is, to a small degree, balanced by the increase in

*F*with increasing

*p*. However, as expected from the strict Lieb–Oxford bound,

^{28}the dependency of

*F*on

*p*can be seen to be small in Fig. 5.

*p*and

*t*values, with the difference being that Eq. (19) is exact. The effective enhancement factors are given for Si and Ge in Table II and confirm the analysis based on the average

*p*and

*t*values. While $FHOeff\u22481$, the LU of Eq. (19) is effectively shifted upwards compared to LDA by a $FLUeff<1$. The feature is surprisingly general. In Fig. 6, the calculated effective enhancement factors are shown for the entire database. While there is some spread, the values are centered around the average values of $FHOeff=1.01$ and $FLUeff=0.93$ close to the values found for Si and Ge, thus indicating that this is quite a general mechanism. The mechanism for opening the mGGA gap thus relies on an upward shift of the LU and is quite different from the simple LDA scaling. In the LDA scaling, both HO and LU were shifted downwards but the HO to a greater extent. Naively, based on the usual GGA naming of

*F*as an enhancement factor, one would have expected

*F*

_{0}to serve a similar role as the simple LDA scaling. However, the name “enhancement factor” is misleading for mGGAs as

*F*< 1 for large areas of

*p*–

*t*space of typical mGGAs. The gap opening comes from these areas being stronger weighted for the LU.

In the context of $GxLDA$ and Fig. 3(a), it was observed how the densities of particularly the LU orbitals are even qualitatively different between Si and Ge. It is noteworthy that despite these differences, the $FLUeff$ values, which are basically density weighted *F*(*p*_{0}, *t*_{0}), are so similar for the two compounds. This similarity is evident in both the visual inspection of Fig. 5, the average values $p\u0304LU,t\u0304LU$, and the effective enhancement factors, $FLUeff$, given in Table II. To gain an understanding of how dissimilar *n*_{LU} densities in Si and Ge result in comparable effective enhancement factors, the atomic contributions to *p*_{0} and *t*_{0} are depicted in Fig. 7. Like in Fig. 3, the colored shaded backgrounds indicate the location of the LU and HO electron. The Si LU is dominated by the outer parts of the 3*s* and 3*p* valence shell, whereas the density of the Ge-LU is much more localized with a contribution from the 3*d* semi-core. However, both the outer part of the valence, such as Si LU region (for *r* > 2.8*a*_{0}), and the inner part of the Ge LU region (1*a*_{0} < *r* < 1.2*a*_{0}) are characterized by larger *p* and *t* values than the HO regions, Fig. 7(a). Thereby, the similar fingerprint of the Si and Ge LU seen in Fig. 5 and the similar effective enhancement factor can be understood.

*p*and

*t*values and the atomic shell structure can be obtained with the electron localization function (ELF).

^{29}The ELF is defined as

*τ*

^{vW}= |∇

*n*|

^{2}/8, and is known to be an indicator of orbital localization.

^{29}It is depicted in Fig. 7(d), and one can indeed relate the atomic shells on Fig. 3(b) with the regions of high ELF in Fig. 7(b). It illustrates how the Si-LU has strong contributions from the valence

*n*= 3 shell characterized by high

*p*and

*t*values. In contrast, the Ge-LU has contributions from the outer part of semi-core shell (which is also

*n*= 3). The ELF locates the semi-core shell approximately between 0.4 and 1.2

*a*

_{0}with the outer part (1

*a*

_{0}<

*r*< 1.2

*a*

_{0}) contributing to higher

*p*and

*t*values.

*dF*

^{eff}, which will depend on functional, orbital, and material,

*dF*

^{eff}values are plotted in Fig. 6(b). Like the

*F*

^{eff}depicted in Fig. 6(a), there is a reasonable transferability. However, it is the HO with typically smaller

*p*and

*t*values that results in the largest shift. Like $\Delta exLDA$ for

*F*

^{eff}, the integrand of Eq. (21) can be viewed as a materials dependent weighting function of the functional dependent part. If the functional was linear with

*p*and

*t*, the effective derivative with respect to the HO, $dFHOeff$, would be the same as that of the LU, $dFLUeff$. However, for any “reasonable” functional,

*F*must be bounded from both below and above. This means that overall the second derivatives will tend to the opposite sign of the first derivative, and assuming no oscillatory behavior, at large

*p*and

*t*values, both partial derivatives in Eq. (22) have to tend to zero, which was also illustrated in Fig. 1.

As a result, the derivative of the derivatives tending toward zero for larger *p* and *t* values and the LU tending toward larger *p* and *t* values $dFLUeff$ will be absolutely smaller than $dFHOeff$. The average values over the 288 compounds in the database for the KTBMgap functional are −0.036 and −0.002 for the HO and LU, respectively. This is also seen in Fig. 6(b) and explains the aforementioned dominance of the HO shift on this part.

Hence, while the bandgap opening due to *F*^{eff} is related to the slope of *F*, the opening due to *dF*^{eff} is related to the curvature of *F*. The effect is weighted by the distribution of *p* and *t* values and is rather small for the small-gap compounds Si and Ge, Fig. 4. However, $\Delta Gxpt$ can be quite significant (Fig. 4). In NaF, the two contributions are similar $\Delta Gxn=2.3$ eV and $\Delta Gxpt=2.2$ eV. Table II shows that this is mainly caused by $dFHOeff$ in agreement with the average values. For the TASK functional, the LU contribution in NaF is actually the largest. However, NaF is less representative for TASK than for KTBMgap. The average values for TASK (−0.028 for HO and −0.005 for LU) are reasonably comparable to the KTBMgap values, thereby indicating that the mechanism is quite general.

## IV. CONCLUSION

To summarize, the success of mGGA functionals for bandgap prediction boils down to the HO and LU being distinct in weighted *p* and *t* space. Based on this observation, we have identified two mechanisms responsible for opening the mGGA gaps. First of all, the tendency of *F* to exhibit relatively large negative derivative with respect to *t* tends to shift the LU state upwards compared to the HO state. Second, the bending of *F* to ensure that it is bounded for large *p* and *t* values tends to move the HO downward compared to the LU. While the actual mechanisms, thus, are complex effects where the decreasing dependency of the functionals on *t* leads to an increase in the smaller bandgaps and the upper and lower limits give a contribution to the larger bandgaps, the mechanisms are found to be surprisingly transferable over a large database of compounds.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Péter Kovács**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Peter Blaha**: Formal analysis (equal); Software (supporting); Writing – review & editing (equal). **Georg K. H. Madsen**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (supporting); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*GW*theory