Bandgap of two-dimensional materials: Thorough assessment of modern exchange–correlation functionals

The calculation of the bandgap with density-functional theory^1,2 (DFT) is computationally efficient and can also be quite accurate, provided that a suitable approximation for the exchange-correlation (xc) effects is chosen. Particularly interesting are the semilocal xc functionals, which are the fastest to evaluate. The most used semilocal xc functionals for the bandgap of bulk solids are the modified Becke–Johnson (MBJ) potential³ and the GLLB-SC potential.^4,5 The absolute errors with respect to the experiment that are obtained with the hybrid functionals (e.g., HSE06^6,7), MBJ, and GLLB-SC are, on average, in the range 0.5–0.8 eV (15%–40% for the absolute relative error) depending on the test set.^8–12 Such errors are much lower than the average errors obtained with the standard Perdew–Burke–Ernzerhof (PBE) functional¹³ that are in the range 1–2 eV (around 50% for the relative error). Proposed recently, the HLE16¹⁴ generalized gradient approximation (GGA), HLE17,¹⁵ MGGAC,¹⁶ TASK,¹⁷ and modified TASK¹⁸ (mTASK) meta-GGAs (MGGA) are examples of other fast semilocal DFT functionals that can be quite accurate for the bandgap of solids as well (see also Refs. 19–23).

Thus, there are fast semilocal xc functionals that are useful alternatives to the much more expensive hybrid functionals and quasiparticle GW methods,²⁴ albeit the latter are more accurate if performed self-consistently and with vertex corrections.^25,26 Actually, it should be mentioned that the dielectric-dependent hybrid functionals are also more accurate than the hybrids using a fixed amount of exact exchange.^27–30 

The majority of the benchmarks of xc functionals for the bandgap are performed using the test sets composed of bulk solids. Furthermore, only bulk solids are typically used when free parameters in an xc functional are optimized.^3,20 Comparatively, there are much less similar studies on atomically thin films, often referred to as two-dimensional (2D) systems. Such works on 2D systems related to the assessment of xc functionals for bandgaps can be found in Refs. 18 and 31–42. We note that an important difference between the bulk solids and 2D systems concerns the magnitude of the excitonic effect (present in optical measurements), which is an effect beyond the Kohn–Sham (KS) and GW methods. It is usually small for most bulk solids (except wide-gap insulators), of the order of tens of meV, but can reach several eV for systems with reduced dimensionality (see, e.g., Ref. 43). Therefore, a direct comparison between the theoretical quasiparticle bandgap and the gap obtained from optical experiments should, in principle, not be made.

Alternatively, fundamental bandgaps calculated using the GW quasiparticle method can be used as reference since GW is viewed as the state-of-the-art in band structure calculations. In Refs. 35, 38, 39, and 42, for instance, the GW bandgaps were used for the testing of various xc DFT functionals on 2D materials (see Ref. 44 for bulk solids). The goal of the present work is to provide a thorough assessment of DFT approximations for the bandgap of 2D materials. This is a follow-up study of our recent works^39,42 devoted to a more systematic comparison of xc functionals. In particular, some of the most recent MGGA functionals will be considered. Studies reporting the benchmark tests of theoretical methods are especially important nowadays for 2D materials.^{32,35,39,45,46} Indeed, these systems can show exceptional properties that are not found in bulk materials and can be very relevant for technological applications.⁴⁷ 

This paper is organized as follows. In Sec. II, a brief description of the tested functionals and test set is given. Section III discusses the results, and Sec. IV gives a summary of the conclusions.

We describe here the tested xc functionals. The PBE functional¹³ is the standard GGA functional for solid-state total-energy calculations; however, it is not accurate for bandgaps. The GGA EV93PW91, which consists of the exchange from Engel and Vosko (EV93)⁴⁸ and Perdew–Wang (PW91) correlation,⁴⁹ slightly improves over PBE for the bandgap.^8,12,50 Also considered is the high-local exchange (HLE16) functional,¹⁴ which is one of the most accurate GGAs for the bandgap of bulk solids.^8,11 We mention that the Armiento–Kümmel (AK13) functional¹⁹ is another GGA that is good for the bandgap of bulk solids;^8,12 however, it leads to numerical problems in 2D materials due to the presence of vacuum. For instance, no self-consistent field (SCF) convergence could be achieved for many of the systems studied here. Therefore, the AK13 will not be considered in the present work.

The MGGA functionals, which represent the next level of approximation,⁵¹ can be more accurate. The ones that are tested here are r²SCAN,⁵² which is a numerically more stable version of the well-known strongly constrained and appropriately normed (SCAN),⁵³ HLE17,¹⁵ TASK¹⁷ and its modified version mTASK,¹⁸ and MGGAC.¹⁶ SCAN and r²SCAN are general purpose functionals that were constructed by satisfying as many constraints as possible. The results for the bandgap of bulk solids show that SCAN gives average errors similar to EV93PW91, i.e., the improvement with respect to PBE is moderate.¹² However, for strongly correlated systems, SCAN is better than EV93PW91.^8,54 Like HLE16, HLE17 was modeled specifically for bandgaps and molecular excitation energies and therefore yields quite accurate bandgaps with errors very similar to HLE16.¹² TASK was constructed nonempirically and also leads to accurate bandgaps of bulk solids, as shown in Refs. 12 and 17. However, TASK does not seem to be a general purpose functional like SCAN/r²SCAN since the lattice constants are very inaccurate.⁵⁵ The same can be probably said about mTASK, which differs from TASK in the value of two parameters that were modified to make the enhancement factor more nonlocal and thus to increase the bandgap.¹⁸ MGGAC contains parameters; some of them were determined using mathematical constraints (e.g., uniform electron gas limit or the tight Lieb–Oxford bound), while the others were fitted to the exchange energy of noble gas atoms or lattice constants of bulk solids. Within the general purpose functionals, MGGAC provides results close to HSE06.

All the xc approximations listed so far are based on an energy functional. This is not the case for GLLB-SC⁵ and local MBJ (LMBJ)⁵⁶ that were modeled at the level of the xc potential and are not derivatives of an energy functional. GLLB-SC is parameter-free and has been shown to be very accurate for the bandgap of bulk solids^{5,9,10,57–59} although it does not perform as well as the MBJ potential, as shown in Refs. 9 and 10. Compared to all other DFT approximations considered here, GLLB-SC differs in the way the bandgap is calculated. While for the other functionals, the bandgap is calculated just as the difference between the conduction band minimum (CBM) and the valence band maximum (VBM), for GLLB-SC, it is calculated by adding a derivative discontinuity to the CBM–VBM difference.⁵ GLLB-SC bandgaps have been calculated for 2D materials in the Computational 2D Materials Database (C2DB),^32,35,60 and it was shown that the agreement with the G₀W₀ (one-shot GW⁶¹) bandgaps is excellent.³⁵ A similar conclusion was drawn in our recent work where selected 2D systems were considered.⁴² 

The LMBJ potential is an adaptation of the MBJ potential for systems with vacuum (molecules, thin films, and surfaces) and interfaces.⁵⁶ MBJ depends on the average of $∇ρ/ρ$ over the unit cell (ρ being the electron density), a quantity that is meaningful in periodic bulk solids, but not in the case of supercells including vacuum and for interfaces. Instead, LMBJ uses a local average of $∇ρ/ρ$ [slightly modified, see Eq. (4)] so that it can be used for any kind of system. The LMBJ potential is given by

where $vxBR$ is the Becke–Roussel potential,⁶² τ is the kinetic-energy density, and c is given by

In Eq. (2),

where

The values of the parameters in Eqs. (2)–(4) chosen by Rauch et al. (see the erratum of Ref. 39) are α = 0.488, β = 0.5 bohr, σ = 3.78 bohr, and ρ_th = 6.96 × 10⁻⁴ e/bohr³. α and β are originally from Ref. 63, while σ is the smearing parameter that determines the size of the region over which the average of g is calculated. Finally, ρ_th is the threshold density, which corresponds to the Wigner–Seitz radius $rsth=(4/3)πρth−1/3=7$ bohr.

A technical aspect of the LMBJ potential and its parents BJ⁶⁴ and MBJ³ as well as AK13¹⁹ should be mentioned. As discussed in Refs. 64–66, these potentials do not tend to zero in the region far from the nuclei, but rather to a system-dependent constant. Since it is customary to set any xc potential to zero or to the LDA value, where ρ is very low (below some chosen density threshold), to avoid numerical instabilities, one has to be careful not to do it in a region of space where the CBM extends. Otherwise, the CBM will be artificially shifted to a wrong energy, which would also lead to a wrong bandgap. Thus, one has to use a density threshold that is small enough to avoid this problem, but also not too small to avoid numerical instabilities. A value of 10⁻⁹ e/bohr³ seems appropriate for LMBJ.

The LMBJ potential was tested in Ref. 39 on the 2D materials of the C2DB; in that work, it was shown that LMBJ is almost as accurate as HSE06 in reproducing the G₀W₀ bandgaps. In the present work, two sets of LMBJ bandgaps will be shown and discussed. One set was obtained with the parameters from the work of Rauch et al., while the other one was obtained with β = 0.6, which leads to results that agree better with G₀W₀.

The potential LB94 of van Leeuwen and Baerends,⁶⁷ which is also not a functional derivative, will be tested as well. Note that the correlation in LB94 is LDA.⁶⁸ 

While the hybrid functionals are, in general, much more accurate than the standard GGA PBE for the bandgap of solids (see Refs. 69–71 for recent works), we will consider only the well-known HSE06 since the focus of the present work is on the fast semilocal functionals. However, this should have no impact on the conclusion of our work since previous benchmark studies^11,12,71 have shown that HSE06 is already among the very best hybrid functionals (excluding those that are dielectric-dependent) for the bandgap of solids.

The test set of 2D materials considered in the present work is the same as the one used by Rauch et al.³⁹ It consists of 298 nonmagnetic systems taken from the C2DB.³⁵ For all of them, the G₀W₀ bandgaps calculated with the GPAW code⁷² are available and will be used as reference values. For 2D materials, experimental bandgaps with excitonic effect subtracted are scarce in the literature. For a few systems, a comparison of the G₀W₀ bandgap with the one inferred from experiment shows that the agreement is pretty good.^34,35 However, it is not excluded that some of the G₀W₀ values that we use as reference are not accurate for various reasons, in particular because of the dependency on the input orbitals and the eigenvalues (PBE was used as the reference ground-state functional³⁵). Nevertheless, due to the large number of systems, we believe that possible inaccuracies in the G₀W₀ data should be small for the statistics and conclusions.

SCF calculations with all xc functionals were performed with the WIEN2k code,^73,74 which is based on the augmented-plane-wave plus local orbital method.^75,76 Spin–orbit coupling was included for all the systems. Parameters such as the basis-set size and number of k-points in the Brillouin zone were chosen to be large enough so that the bandgap is converged to within a few 0.01 eV. The self-consistent implementation of MGGA functionals in WIEN2k is very recent⁵⁵ and uses the subroutines from the library of exchange–correlation functionals—Libxc.^77,78

The functionals listed in Sec. II are used to calculate the bandgap of 298 of the 2D materials in the C2DB. As mentioned, the 298 systems were selected by Rauch et al.³⁹ and the G₀W₀ results are used as reference. Table I shows various statistical quantities: the mean (percentage) error [M(P)E], mean absolute (percentage) error [MA(P)E], and standard (percentage) deviation [S(P)D]. The coefficients a and b of a linear fit as well as the Pearson correlation coefficient r are also listed in Table I. For comparison purposes, the results from Ref. 35 obtained with the screened hybrid functional HSE06^6,7 using the GPAW code⁷² are shown, as well. The numerical values of the bandgaps for all the materials and xc functionals can be found in the supplementary material. The values in Table I are also shown graphically in Fig. 1.

The errors obtained with the PBE are large since the MAE is 1.50 eV and the MAPE is 51%. The two other GGAs improve over PBE. The improvement is moderate with EV93PW91 but more visible with HLE16. The latter leads to a MAE close to 1 eV and a MAPE of 39%. The LB94 potential leads to even worse results than the PBE. Thus, the correct asymptotic behavior −1/r of the LB94 potential far from the nuclei does not seem to be as useful for 2D materials as it is for molecules.⁷⁹ 

Turning to the MGGAs, r²SCAN and HLE17 perform roughly as the GGAs EV93PW91 and HLE16 with a MAE slightly above 1 eV and a MAPE just below 40%. The three other MGGAs are more accurate, and this is particularly the case of mTASK, which leads to a MAE of 0.52 eV and a MAPE of 21%. Therefore, mTASK and TASK are more accurate than the expensive hybrid functional HSE06. However, the DFT bandgaps that agree best with the G₀W₀ reference values are those obtained with GLLB-SC. The MAE and MAPE with GLLB-SC are as low as 0.42 eV and 21%, respectively. In passing, we note that our MAE obtained with GLLB-SC agrees very well with the MAE of 0.38 eV reported by Haastrup et al.³⁵ for another subset of about 250 materials of the C2DB.

As previously observed,³⁹ the performance of LMBJ with the original parameters (β = 0.5) and HSE06 is pretty similar, but of course, the advantage of the LMBJ potential is to be computationally more efficient and to scale better with the system size. The MAE of LMBJ with β = 0.5 is 0.78 eV, while the MAPE is 35%. However, if the parameter β in Eqs. (2) and (4) is increased to 0.6, then the MAE and MAPE are reduced. This is especially the case for the MAE, which is now 0.50 eV, similar to mTASK.

The individual results are shown graphically for most of the DFT functionals in the panels of Fig. 2. This is a convenient way to get an idea of the S(P)D and the linear fit coefficients. From Table I, we can see that the SD is the lowest (in the range 0.50–0.56 eV) for mTASK, GLLB-SC, TASK, and LMBJ (β = 0.5) and the largest (close to 1 eV) for LB94, PBE, and EV93PW91. In terms of SPD, the lowest value is 16%, obtained with EV93PW91 (one of the worst functionals for the SD). The largest values, around 150%, are obtained with the MGGAC and LMBJ (β = 0.6). Concerning the linear fit, the functionals with the slope (a) that is closest to 1 are GLLB-SC and LMBJ (β = 0.6), which is also visible from Fig. 2. PBE leads to the worst value of a (0.57). For the offset b, it is noteworthy that the GLLB-SC along with the LB94 leads to the values that differ the most from zero. Regarding the correlation coefficient r, most functionals have a value of 0.97 or 0.98, while LMBJ (β = 0.6) leads to 0.96 and LB94 to 0.95. Thus, except for LB94, which is the worst functional for most quantities, the correlation coefficient does not really seem to be a useful quantity.

Concerning the outliers visible in Fig. 2, we note the following. There are a couple of materials with small bandgaps, more particularly FeCl₂, FeBr₂, and FeI₂, that are strongly overestimated with HSE06 and LMBJ, and to a lesser degree with MGGAC, GLLB-SC, and mTASK. Nevertheless, it is legitimate to question the accuracy of the G₀W₀ bandgap for the three FeX₂ systems, in particular, since HSE06 shows large deviations, which is somehow surprising. The correlation in the Fe atom may be quite strong such that PBE is not the appropriate functional to generate the orbitals and eigenvalues for G₀W₀.⁸⁰ Furthermore, it should be noted that the real ground state of those transition-metal halides, as well as of some other systems, is ferromagnetic.

The two largest bandgaps, for CaF₂ and SrF₂, are largely overestimated with GLLB-SC. Concerning the number of false metals (shown in Table I), we note that there are seven such cases with LB94, as for instance AsIn (1.69 eV with G₀W₀) or PtTe₂ (1.24 eV with G₀W₀). HLE16 leads to three false metals, and HLE17 and GLLB-SC lead to one such case, while no false metals are obtained with the other functionals. In this respect, LB94 is the worst functional, which is expected since it leads to the strongest underestimation in the bandgap. As discussed in previous works,^8,81 HLE16 and HLE17 have rather oscillatory potentials and are therefore unpredictable.

By considering all the results discussed so far, we can conclude that the best approximation for the calculation of the fundamental bandgap of 2D materials is GLLB-SC. It is, in fact, the best (or nearly the best) approximation for all the statistical quantities, except the SPD and coefficient b of the linear fit. mTASK, which is the second best approximation, is overall rather close to GLLB-SC. LMBJ (β = 0.6) and TASK are also pretty accurate, actually at least as accurate as the hybrid HSE06. However, as discussed above, the SPD of LMBJ (β = 0.6) is very large, which is due to the large spread of the errors for the materials with small bandgaps.

In one of our previous works on 2D materials,⁴² plots of the xc potentials were shown to provide an explanation to some of the observed trends. For instance, for XS₂, XSe₂, and XTe₂, the PBE and LMBJ (β = 0.5) bandgaps are basically the same when X = Mo or W, while it is not at all the case when X = Zr or Hf. We could understand these results by comparing the curves of the PBE and LMBJ potentials. Here, we consider CrO₂, which, as the other 2D systems with Cr, seems to be problematic for all the semilocal functionals. The G₀W₀ and HSE06 bandgaps are 1.23 and 1.14 eV, respectively; however, the semilocal methods give values that are in the range 0.19–0.46 eV, which is much smaller. More specifically, LMBJ (β = 0.6) leads to the smallest bandgap (0.19 eV), while a value around 0.45 eV is obtained by PBE, EV93PW91, r²SCAN, (m)TASK, and GLLB-SC. Figure 3 shows the VBM and CBM in CrO₂, where we can see that both are of Cr-$3dz2$ character (the VBM has some O-2p_z component). In such a case where the VBM and CBM are located in the same position of space and are, in addition, of the same character, an opening of the bandgap is difficult to achieve with a multiplicative potential (a GGA or LMBJ). Actually, for CrO₂, even the non-multiplicative MGGAs and GLLB-SC cannot increase the bandgap with respect to PBE. Other systems with a d–d gap that is underestimated with all the semilocal functionals are the other isoelectronic Cr materials (e.g., CrS₂ or CrSTe), TiCl₂, and those of the same family.

For this work, we also investigated if an improvement of the LMBJ potential as originally proposed in Ref. 39 was possible. A full scan of the four-dimensional space of the parameters α, β, σ, and $rsth$ in Eqs. (2)–(4) would be tedious. Therefore, only one or two parameters were simultaneously varied, which should be sufficient to tell us whether a clear improvement of the LMBJ accuracy is possible or not. Without going too much into detail, we observed the following. The MAE can be reduced by increasing either α or β, and this, as exemplified above using β = 0.6, improves the results compared to the original value of β = 0.5. However, the MAPE was barely reduced (from 35% to 32%), while there was a clear increase in the SPD. The latter quantity is nearly doubled. Up to some point, increasing β further would continue to lower the MAE and rise the MAPE, SD, and SPD. We also observed that a larger β leads to a larger bandgap for the vast majority of the materials. Changing the value of σ and/or $rsth$ leads to a deterioration of the results, which is rather expected since the original values were already optimized for the 2D materials by Rauch et al.³⁹ β = 0.6 (with α, σ, and $rsth$ unmodified) is a choice, among others, that leads to rather well balanced errors overall if one considers the four mean errors and coefficient a as the most important quantities. As discussed above, LMBJ is quite satisfying overall and does not lag far behind GLLB-SC and mTASK.

We mention that in the search of an alternative to Eq. (3) for systems with vacuum, a possibility could be

where

In Eqs. (5) and (6), the contribution to the integral comes only from the region of space where the electron density ρ has a non-negligible value. Although potentially interesting, Eq. (5) has not led to improved results. However, a more careful exploration needs to be performed.

To finish this section, we mention that LMBJ has also been tested by some of us⁸² for the ionization potential (IP) of molecules and the electronic properties of the system consisting of a F6-TCNNQ molecule adsorbed on a hydrogenated Si(111) surface. While LMBJ is accurate for the IP of molecules, it is not so for the charge transfer between the molecule and the surface.

In this work, we tested a variety of xc functionals for the calculation of the bandgap of 2D materials. The test set comprises 298 2D materials for which G₀W₀ bandgaps are available and were used as reference. The tested xc functionals are the most accurate among those that are currently available for bandgaps. The results show that the two most accurate are the GLLB-SC potential and the mTASK functional. The LMBJ (β = 0.6) potential and TASK functional can also be considered as accurate and follow quite closely GLLB-SC and mTASK.

At this point, it is also important to remind some technical aspects. GLLB-SC and LMBJ have no associated energy functional, which is inconvenient (e.g., no geometry optimization is possible and the zero-force and zero-torque conditions are not satisfied⁸³). Furthermore, GLLB-SC depends on the eigenvalues and the Fermi energy such that some kind of size-consistency is not satisfied. All these problems do not occur with MGGAs such as TASK, which consist of an energy functional. However, we should also mention that the lattice parameters of bulk solids are very inaccurate (strongly overestimated) with TASK.⁵⁵ This is probably the unavoidable price to pay when a GGA or MGGA energy functional is constructed with the aim of providing very accurate bandgaps.

The results presented in this work represent the most comprehensive study about the performance of the semilocal xc functionals for 2D materials. They can serve as a guide for applications and future development of xc functionals. However, we believe that it will be difficult to achieve a better accuracy than GLLB-SC or mTASK, at least against the G₀W₀ results that we have used here as reference.

See the supplementary material for the bandgap of all materials calculated with all xc methods.

J. Doumont and P. Blaha acknowledge support from the Austrian Science Fund (FWF) through Project No. W1243 (Solids4Fun). L. Kalantari acknowledges support from the TU-D doctoral college (TU Wien). S. Botti and P. Borlido were supported by the European Commission in the framework of the H2020 FET Open project SiLAS (Grant Agreement No. 735008). M. A. L. Marques and S. Botti acknowledge partial support from the DFG through the projects TRR 227 (Project No. B09), SFB-1375 (Project No. A02), and BO 4280/8-1. S. Botti and T. Rauch acknowledge funding from the Volkswagen Stiftung (Momentum) through the project “dandelion.” S. Jana acknowledges funding from NISER, Bhubaneswar, India.

The data that support the findings of this study are available within this article and its supplementary material.