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

The density functional theory (DFT) approximations that are the most accurate for the calculation of band gap of bulk materials are hybrid functionals like HSE06, the MBJ potential, and the GLLB-SC potential. More recently, generalized gradient approximations (GGA), like HLE16, or meta-GGAs, like (m)TASK, have proven to be also quite accurate for the band gap. Here, the focus is on 2D materials and the goal is to provide a broad overview of the performance of DFT functionals by considering a large test set of 298 2D systems. The present work is an extension of our recent studies [Rauch et al., Phys. Rev. B 101, 245163 (2020) and Patra et al., J. Phys. Chem. C 125, 11206 (2021)]. Due to the lack of experimental results for the band gap of 2D systems, $G_{0}W_{0}$ results were taken as reference. It is shown that the GLLB-SC potential and mTASK functional provide the band gaps that are the closest to $G_{0}W_{0}$. Following closely, the local MBJ potential has a pretty good accuracy that is similar to the accuracy of the more expensive hybrid functional HSE06.


I. INTRODUCTION
The calculation of the band gap with density functional theory 1,2 (DFT) is computationally efficient, but 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 band gap of bulk solids are the modified Becke-Johnson (MBJ) potential 3 and the GLLB-SC potential. 4,5 The absolute errors with respect to 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][9][10][11][12] Such errors are much lower than the average errors obtained with the standard PBE functional 13 that are in the range 1−2 eV (around 50% for the relative error). Proposed recently, the HLE16 14 generalized gradient approximation (GGA), and the HLE17, 15 MGGAC, 16 TASK, 17 and modified TASK 18 (mTASK) meta-GGAs (MGGA) are examples of other fast semilocal DFT functionals that can be quite accurate for the band gap of solids as well (see also . Thus, there are fast semilocal xc functionals that are useful alternatives to the much more expensive hybrid functionals and quasiparticle GW methods, 24 albeit the latter are more accurate if performed self-consistently and with vertex corrections. 25,26 Actually, it should be mentioned that dielectric-dependent hybrid functionals are also more accurate than the hybrids using a fixed amount of exact exchange. [27][28][29][30] The majority of benchmarks of xc functionals for the band gap are done using test sets composed of bulk solids. Furthera) Electronic mail: tran@theochem.tuwien.ac.at more, only bulk solids are typically used when free parameters in a 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 band gaps can be found in Refs. 18, 31-42. We note that an important difference between 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 band gap and the gap obtained from optical experiments should in principle not be done.
Alternatively, fundamental band gaps calculated from the GW quasiparticle method can be used as reference, since GW is viewed as the state-of-the-art for band structure calculations. In Refs. 35, 38, 39, and 42 for instance, GW band gaps 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 band gap of 2D materials. It 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 benchmark tests of theoretical methods are especially important for 2D materials nowadays. 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. 47 The 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.

II. METHODS AND COMPUTATIONAL DETAILS
We describe below the tested xc functionals. The Perdew-Burke-Ernzerhof (PBE) functional 13 is the standard GGA functional for solid-state total-energy calculations, however it is not accurate for band gaps. The GGA EV93PW91, which consists of the exchange from Engel and Vosko (EV93) 48 and Perdew-Wang (PW91) correlation 49 slightly improves over PBE for the band gap. 8,12,50 Also considered is the high local exchange (HLE16) functional, 14 which is one of the most accurate GGAs for the band gap of bulk solids. 8,11 We mention that the Armiento-Kümmel (AK13) functional 19 is another GGA that is good for the band gap 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, AK13 will not be considered in the present work.
The MGGA functionals, which represent the next level of approximation, 51 can be more accurate. The ones that are tested here are r 2 SCAN, 52 which is a numerically more stable version of the well-known strongly constrained and appropriately normed (SCAN), 53 HLE17, 15 TASK 17 and its modified version mTASK, 18 and MGGAC. 16 SCAN and r 2 SCAN are general purpose functionals that were constructed by satisfying as many constraints as possible. Results for the band gap of bulk solids show that SCAN gives average errors similar to EV93PW91, i.e., the improvement with respect to PBE is moderate. 12 However, for strongly correlated systems SCAN is better than EV93PW91. 8,54 Like HLE16, HLE17 was modelled specifically for band gaps and molecular excitation energies, and therefore yields quite accurate band gaps with errors very similar to HLE16. 12 TASK was constructed nonempirically and also leads to accurate band gaps of bulk solids as shown in Refs. 12 and 17. However, TASK does not seem to be a general purpose functional like SCAN/r 2 SCAN, since the lattice constants are very inaccurate. 55 The same can probably be said about mTASK, which differs from TASK by the value of two parameters that were modified to make the enhancement factor more nonlocal and thus to increase the band gap. 18 MGGAC contains parameters, some of them were determined using mathematical constraints (e.g., uniform electron gas limit or the tight Lieb-Oxford bound), while 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 5 and local MBJ (LMBJ) 56 that were modelled at the level of the xc potential and are not derivative of an energy functional. GLLB-SC is parameter-free and has been shown to be very accurate for the band gap of bulk solids, 5,9,10,[57][58][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 band gap is calculated. While for the other functionals the band gap is calculated just as the difference between the conduction band minimum (CBM) and valence band maximum (VBM), for GLLB-SC it is calculated by adding a derivative discontinuity to the CBM−VBM difference. 5 GLLB-SC band gaps have been calculated for 2D materials in the computational 2D materials database (C2DB) database, 32,35,60 and it was shown that the agreement with the G 0 W 0 (one-shot GW 61 ) band gaps is excellent. 35 A similar conclusion was drawn in our recent work where selected 2D systems were considered. 42 The LMBJ potential is an adaptation of the MBJ potential for systems with vacuum (molecules, thin films, surfaces) and interfaces. 56 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 case of supercells including vacuum and for interfaces. Instead, LMBJ uses a local average of |∇ρ| /ρ [that is slightly modified, see Eq. (4)] so that it can be used for any kind of system. The LMBJ potential is given by where v BR x is the Becke-Roussel potential, 62 τ is the kineticenergy density, and c is given by In Eq. (2)ḡ where  63, while σ is the smearing parameter that determines the size of the region over which the average of g is done. Finally, ρ th is the threshold density, which corresponds to a Wigner-Seitz radius r th A technical aspect of the LMBJ potential and its parents BJ 64 and MBJ 3 as well as AK13 19 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 systemdependent 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 band gap. 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 −9 e/bohr 3 seems appropriate for LMBJ.
The LMBJ potential was tested in Ref. 39 on the 2D materials of the C2DB database; in that work it was shown that LMBJ is almost as accurate as HSE06 in reproducing the G 0 W 0 band gaps. In the present work, two sets of LMBJ band gaps will be shown and discussed. One set was obtained with the parameters from Rauch et al. listed above, while the other one was obtained with β = 0.6, which leads to results which agree better with G 0 W 0 .
The potential LB94 of van Leeuwen and Baerends, 67 which is also not a functional derivative, will be tested as well. Note that the correlation in LB94 is LDA. 68 While hybrid functionals are in general much more accurate than the standard GGA PBE for the band gap 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 which are dielectric-dependent) for the band gap of solids.
The test set of 2D materials considered in the present work is the same as the one used by Rauch et al. 39 It consists of 298 nonmagnetic systems taken from the C2DB database. 35 For all of them, G 0 W 0 band gaps calculated with the GPAW code 72 are available and will be used as reference values. For 2D materials, experimental band gaps with excitonic effect subtracted are scarce in the literature. For a few systems, comparison of the G 0 W 0 band gap with the one inferred from experiment shows that the agreement is pretty good. 34,35 However, it may well be that some of the G 0 W 0 values that we use as reference are not as accurate for various reasons, in particular because of the dependency on the input orbitals and eigenvalues (PBE was used as the reference ground-state functional. 35 ) Nevertheless, due to the large number of systems, we believe that possible inaccuracies in the G 0 W 0 data should be small for the statistics and conclusions.
The SCF calculations with all xc functionals were done with the WIEN2k code, 73,74 which is based on the augmentedplane-wave plus local orbitals method. 75,76 Spin-orbit coupling was included for all systems. Parameters like the basisset size and number of k-points in the Brillouin zone were chosen to be large enough so that the band gap is converged to within a few 0.01 eV. The self-consistent implementation of MGGA functionals in WIEN2k is very recent 55 and uses the subroutines from the library of exchange-correlation functionals Libxc. 77,78

III. RESULTS AND DISCUSSION
The functionals listed in Sec. II were used to calculate the band gap of 298 of the 2D materials in the C2DB database. As mentioned, the 298 systems were selected by Rauch et al. 39 and G 0 W 0 results are used as reference. Table I shows Table I are also shown graphically on Fig. 1.
The errors obtained with 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 PBE. Thus, the correct asymptotic behavior −1/r of the LB94 potential far from nuclei does not seem to be as useful for 2D materials as it is for molecules. 79 Turning to the MGGAs, r 2 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 as well as TASK are more accurate than the expensive hybrid functional HSE06. However, the DFT band gaps that agree best with the G 0 W 0 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. 35 for another subset of about 250 materials of the C2DB database.
As previously observed 39 the performance of LMBJ with the original parameters (β = 0.5) and HSE06 are pretty similar, but of course the advantage of the LMBJ potential is to be computationally more efficient and to scale better with 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 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 MGGAC and LMBJ(β = 0.6). Concerning the linear fit, the functionals with the slope (a) that is the 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 GLLB-SC along with LB94 lead 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 band gaps, more particularly FeCl 2 , FeBr 2 , and FeI 2 , that are strongly  I. Summary statistics for the error in the DFT band gaps with respect to G 0 W 0 reference values for the set of 298 2D materials. M(P)E, MA(P)E, and S(P)D denote the mean (percentage) error, mean absolute (percentage) error, and standard (percentage) deviation, respectively. a and b are the coefficients of the linear fit y = ax + b (shown in Fig. 2), r is the Pearson correlation coefficient, and the last column is for the number of materials wrongly described as metallic. The type of approximation of the functionals is indicated in parenthesis (GLLB-SC depends on the eigenvalues ε i and the exchange derivative discontinuity ∆ x is added to CBM − VBM). The units of the ME, MAE, SD, and b are eV. 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 0 W 0 band gap for the three FeX 2 systems, in particular since also HSE06 shows large deviation, which is somehow surprising. The correlation on the Fe atom may be quite strong, such that PBE is not the appropriate functional to generate the orbitals and eigenvalues for G 0 W 0 . 80 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 band gaps, for CaF 2 and SrF 2 , are largely overestimated with GLLB-SC. Concerning the number of false metals (shown in Table I), we note that there are 7 such cases with LB94, as for instance AsIn (1.69 eV with G 0 W 0 ) or PtTe 2 (1.24 eV with G 0 W 0 ). HLE16 leads to three false metals, HLE17 and GLLB-SC 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 band gap. As discussed in previous works, 8,81 HLE16 and HLE17 have rather oscillatory potentials and are therefore unpredictable.
By considering all results discussed so far, we can conclude that the best approximation for the calculation of the fundamental band gap of 2D materials is GLLB-SC. It is in fact the best (or nearly the best) approximation for all 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, and actually at least as accurate as the hybrid HSE06. However, as noted above, the SPD of LMBJ(β = 0.6) is very large, which is due to a large spread of the errors for the materials with small band gaps.
In one of our previous works on 2D materials, 42 plots of the xc potentials were shown to provide an explanation to some of the observed trends. For instance, for XS 2 , XSe 2 , and XTe 2 , the PBE and LMBJ(β = 0.5) band gaps 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 2 , which, as the other 2D systems with Cr, seems to be problematic for all semilocal functionals. The G 0 W 0 and HSE06 band gaps 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 band gap (0.19 eV), while a value around 0.45 eV is obtained by PBE, EV93PW91, r 2 SCAN, (m)TASK, and GLLB-SC. Figure 3 shows the VBM and CBM in CrO 2 , where we can see that both are of Cr-3d z 2 character (the VBM has also 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 band gap is difficult to achieve with a multiplicative potential (a GGA or LMBJ). Actually, for CrO 2 even the non-multiplicative MGGAs and GLLB-SC can not increase the band gap with respect to PBE. Other systems with a d−d gap that is underestimated with all semilocal functionals are the other isoelectronic Cr materials (e.g., CrS 2 or CrSTe), but also for instance TiCl 2 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 r th s in Eqs.
(2)-(4) would be tedious. Therefore, only one or two parameters were simultaneously varied, which however should be sufficient to tell us if 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 β = 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 band gap for the vast majority of materials. Changing the value of σ and/or r th s leads to a deterioration of the results, which is rather expected since the original values were already optimized for 2D materials by Rauch et al. 39 β = 0.6 (with α, σ , and r th s 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 bẽ 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 done.
To finish this section, we mention that LMBJ has also been tested by some of us 82 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.

IV. SUMMARY
In this work we tested a variety of xc functionals for the calculation of the band gap of 2D materials. The test set comprises 298 2D materials for which G 0 W 0 band gaps are available and were used as reference. The tested xc functionals are the most accurate currently available for band gaps. 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 that 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 83 ). Furthermore, GLLB-SC depends on the eigenvalues and the Fermi energy, such that some kind of sizeconsistency is not satisfied. All these problems do not occur with MGGAs like TASK which consist of a energy functional. However, we should also mention that the lattice parameters of bulk solids are very inaccurate (strongly overestimated) with TASK. 55 This is probably the unavoidable price to pay when a GGA or MGGA energy functional is constructed with the aim to provide very accurate band gaps.
The results presented in this work represent the most comprehensive study about the performance of 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 0 W 0 results that we have used here as reference.

SUPPLEMENTARY MATERIAL
See the supplementary material for the band gap of all materials calculated with all xc methods. ACKNOWLEDGMENTS J. Doumont and P. Blaha acknowledge support from the Austrian Science Fund (FWF) through project W1243 (Solids4Fun). L. Kalantari acknowledges support from the TU-D doctoral college (TU Wien). S. Botti and P. Borlido are supported by the European Commission in the framework of the H2020 FET Open project SiLAS (GA no. 735008). M. A. L. Marques and S. Botti acknowledge partial support from the DFG though the projects TRR 227 (project B09), SFB-1375 (project A02), and BO 4280/8-1, respectively. S. Botti and T. Rauch acknowledge funding from the Volkswagen Stiftung (Momentum) through the project "dandelion". S. Jana acknowledges funding from NISER, Bhubaneswar, India.

DATA AVAILABILITY
The data that supports the findings of this study are available within the article [and its supplementary material].