Validation of Pseudopotential Calculations for the Electronic Band Gap of Solids

Nowadays pseudopotential (PP) density functional theory calculations constitute the standard approach to tackle solid-state electronic problems. These rely on distributed PP tables that were built from all-electron atomic calculations using few popular semilocal exchange-correlation functionals, while PPs based on more modern functionals, such as meta-generalized gradient approximation and hybrid functionals, or for many-body methods, such as GW, are often not available. Because of this, employing PPs created with inconsistent exchange-correlation functionals has become a common practice. Our aim is to quantify systematically the error in the determination of the electronic band gap when cross-functional PP calculations are performed. To this end, we compare band gaps obtained with norm-conserving PPs or the projector-augmented wave method with all-electron calculations for a large data set of 473 solids. We focus, in particular, on density functionals that were designed specifically for band gap calculations. On average, the absolute error is about 0.1 eV, yielding absolute relative errors in the 5–10% range. Considering that typical errors stemming from the choice of the functional are usually larger, we conclude that the effect of choosing an inconsistent PP is rather harmless for most applications. However, we find specific cases where absolute errors can be larger than 1 eV or others where relative errors can amount to a large fraction of the band gap.


Introduction
Since its origin more than 50 years ago, density functional theory 1,2 (DFT) has become the standard approach to tackle the electronic structure of solids.A workable approach to DFT is attained via the Kohn-Sham formulation, 2 leading to equations that can be solved efficiently with modern computational resources.[5][6] From a purely theoretical point of view, the xc functional is the only approximation in DFT.However, in practice, the Kohn-Sham equations -a system of coupled, non-linear, partial differential equations -are solved numerically, introducing further approximations.
Several different approaches have been developed and are commonly used in the physics and chemistry communities.A basic distinction exists between all-electron and pseudopotential (or effective-core potential) methods.In the former, all electrons are explicitly included in the calculation, and the electron-nuclear attraction is described by the standard Coulomb potential.The latter method is based on the distinction between core and valence electrons.
By replacing the effect of the nucleus and the core electrons by an effective pseudopotential, efficient plane-wave 7 or real-space 8,9 techniques can be used to solve the Kohn-Sham equations.As a middle way, the projector-augmented wave method (PAW) 10,11 was developed, combining the advantages of pseudopotentials with a reconstruction of the all-electron wavefunction.
In materials science DFT calculations are often performed with effective-core methods. 12,13It is generally believed that well-tuned pseudopotentials can yield very precise calculations of many properties of solids at a much lower cost than all-electron methods.
However, in order to guarantee this, careful optimization is necessary to ensure the quality and transferability of the pseudopotential for the situation at hand.Several well-tested tables are available to the community, but all these pseudopotentials are built based on two approximations to the xc energy functional: the local density approximation (LDA) [14][15][16][17] ) and the generalized-gradient approximation (GGA), and more specifically the Perdew-Burke-Ernzerhof (PBE) 18 functional, which is the de facto standard in the physics community.This is true for norm-conserving 19 and ultra-soft 20 potentials, as well as for the PAW method. 106][27][28][29] As a result of this situation, the great majority of calculations performed with improved functionals makes use of an inconsistent pseudopotential, built from LDA or GGA atomic calculations.This problem is relevant for both total-energy and band-structure calculations.For the latter, which are the focus of the present work, the most accurate calculations use either hybrid functionals (such as the Heyd-Scuseria-Ernzernhof 2006 functional 28,29 ), GW methods 30,31 or specialized semilocal functionals, such as the modified Becke-Johnson potential 32 (mBJLDA).For example, despite the publication of several works on the generation of pseudopotentials for Hartree-Fock and hybrid functionals, [33][34][35][36][37][38] the lack of readily available tables means that most of these calculations are done with PBE ones.Furthermore, some codes do not give the user the possibility of changing the pseudopotentials.The situation is even more complicated for the mBJLDA potential, as it is defined for periodic systems and can not be applied to atoms, unless specific schemes are used. 39,40Note that a similar problem exists also for many-body GW 30,31 or LDA+U 41,42 calculations, that are often performed with LDA or PBE pseudopotentials.This common practice creates an inconsistency, which invariably introduces some uncontrollable errors in the calculations.Even if this is well-known, relatively little attention has been given to quantify the effect of cross-functional pseudopotential calculations. 36,39,43,44th this in mind, we decided to study the error coming from using an inconsistent pseudopotential in the calculation of band gaps.This was done by comparing the results obtained with three codes: WIEN2K, 45,46 an all-electron code that uses the augmented plane wave plus local orbitals (APW+lo) method, 47 ABINIT 48,49 with norm-conserving pseudopotentials, and VASP 11,50 with PAW setups.Note that our comparison goes beyond the ∆-test, 13 as we compare band gaps, that also include the influence of unoccupied bands, and we take into account errors coming from using inconsistent pseudopotentials.

Dataset
We performed calculations for all materials contained in the dataset developed in Ref., 51 that counts 473 non-magnetic semiconductors.This dataset covers the majority of the periodic table, and includes materials with a wide range of band gaps.All calculations were performed at the experimental geometry (see Ref. 51 for more details on the materials dataset).

Functionals
Besides the standard LDA and PBE functionals, our choice of xc functionals (or potentials) was based on two criteria: (i) their availability in all three software packages we used (see Sec. for more information on the codes) and (ii) their relevance for the calculation of band gaps.Our final choice includes the Perdew-Wang (PW92) 16 and the Perdew-Zunger (PZ81) 15 parametrization of the LDA correlation, the local Slater potential (SLOC) 52 ), which is a simple modification of the LDA exchange to approximate the Slater potential.We consider also several GGA functionals: PBE, 18 revised PBE (RPBE), 53 Engel-Vosko (EV93) 54 combined with Perdew-Wang (PW91) correlation 55 and high-local exchange (HLE16). 56Finally we used the meta-GGA mBJLDA. 32The PW92, PZ81 and PBE functionals are generalpurpose approximations with widespread use in the solid-state community for total-energy calculations.The mBJLDA potential is known to yield excellent band gaps, which is also the case of the simpler SLOC and HLE16, although they are slightly inferior to mBJLDA. 51,57te that meta-GGA energy functionals are not supported self-consistently in the reference APW+lo WIEN2k code, and therefore we did not include in our analysis the recent meta-GGAs HLE17 58 and TASK, 59 which also perform very well for band gaps.Hybrid functionals are available in the three considered codes, however they lead to calculations which are several orders of magnitude more expensive, in particular if parameters for highly converged calculations are used.Therefore, since our test set is very large, we refrained from using hybrid functionals in the present work.For the same reason we did not consider GW methods.

Codes
The all-electron calculations were done with the WIEN2k package, 45,46 which uses the APW+lo basis set. 47The calculations were done with sufficiently large parameters (e.g., basis-set size) to ensure convergence of band gaps within ∼ 0.05 eV.The WIEN2k results will serve as reference for the comparison with the other codes.
The PAW calculations were performed using a custom version of the Vienna ab initio simulation package (VASP; version 5.4) 11,50 that is interfaced to LIBXC . 22,23A plane wave cutoff of 520 eV was used for all species along with the same k-grids as in Ref. 51 All meta-GGA calculations were performed accounting for non-spherical contributions of the density gradient inside the augmentation spheres.VASP calculations are in general restricted to the PAW sets included in the distribution, and we were therefore able to use PBE and LDA PAW pseudopotentials.
The ABINIT package 48 was used to test norm-conserving pseudopotential calculations.
For LDA and PBE pseudopotentials we resorted to the Pseudo Dojo distribution 60 (version 0.4, with stringent accuracy).Although this set covers most of the periodic table, some elements like thorium are absent and had unfortunately to be left out of the calculations.
Since ABINIT does not currently support pseudopotentials with non-linear core corrections for meta-GGAs, mBJLDA calculations with this code were not performed.For SLOC and HLE16 we generated a set of pseudopotentials using the ONCVPSP package 61 (version 3.3.1).As a starting point we used the input files from the stringent set of the Pseudo Dojo distribution, changing them to include the desired functionals.Whenever a particular input leads to unsuccessful/difficult calculations, we exchanged it for that of the standard set.Further small changes to the local part of the pseudopotential were performed in order to avoid spurious effects (e.g.ghost states) that were detected by the post-processing tools of ONCVPSP.We note that although care was taken in this process, we did not perform further tests (∆-test, GBVR test, etc. [62][63][64] ) or optimizations.Therefore, and although these pseudopotentials yield generally consistent results for band gaps, care should be taken for a more general use.The pseudopotentials are given as Supplementary Information, and the whole set can also be downloaded from. 65te that the three codes, as well as ONCVPSP, are linked to the library of xc functionals LIBXC, 22,23 which allows to access several hundreds of functionals, including the ones considered here.
All calculations were performed neglecting spin-orbit coupling.This term is expected to contribute on average about 0.1 eV to the band gap.This amount is considerably smaller than the typical average error of the xc functionals.Anyway, as all calculations were performed consistently without this term, its exclusion does not affect our comparison between the codes, that is the main purpose of the present work.

Statistics
For the analysis of the results we compare band gaps calculated with norm-conserving pseudopotentials and PAW methods to all-electron values.Our analysis is restricted to those materials that were not determined by WIEN2k to have a theoretical band gap smaller than 0.01 eV, despite being measured to be semiconductors.From the original 473 entries, WIEN2k predicts between 13 and 40 of such materials, depending on the functional.
Whenever presenting the results we use the notation <calculation xc>@<pseudo xc>, and if the code used is ambiguous, we precede this string with its name.For example, ABINIT:PBE@LDA would make reference to the set of values computed with ABINIT, using the PBE calculations with LDA pseudopotentials.
The statistical analysis is based on the determination of the mean absolute error, MAE = i=1 (y i − y i,exp )/(n y i,exp ); and the maximum absolute error in the calculation of band gaps with respect to experimental values.The complete set of results is presented in Tables SI-SVII of the Supplementary Information, while Table 1 shows a summary of the most important statistical quantities.
After a preliminary analysis of the results, it became apparent that relative quantities (such as the MAPE) were being extremely affected by materials with very small band gaps.This is easy to understand as small errors lead to a rather large relative error for systems with band gaps in the range 0.1-0.2eV, skewing significantly the statistical averages.Therefore we opted to consider in Table 1 percentage quantities (MAPE and MPE) for the subset of systems with band gaps larger than 0.25 eV.Absolute quantities were still computed for the entire dataset.

Results
We start our analysis by looking at the results computed with the generally available LDA and PBE datasets.As visible in Table 1, band-gap calculations performed with LDA and PBE xc functionals on top of the corresponding pseudopotentials are in excellent agreement with all-electron calculations.Not only are the MAE and ME for these calculations very small (in absolute value smaller than 0.03 eV), but also the dispersion of results is quite localized.This is visually represented in the error histograms of Fig. 1, and also in the corresponding standard deviation σ (smaller than 0.06 eV).The maximum absolute error in these conditions is around 0.2-0.3eV.However, one can see from the distribution of errors that absolute values larger than 0.1 eV are rare.These conclusions are valid for both pseudopotential and PAW calculations, although the MPE and MAPE of the latter are consistently larger.This may be explained by the fact that the Pseudo Dojo sets 60 are much more recent than the PAW sets available in VASP, and that they were systematically optimized.In any case, these errors are certainly acceptable for the large majority of applications.This confirms the generally accepted idea that effective-core methods (either norm-conserving or PAW approach) are reliable for the calculations of band gaps.
Moving to the cross-functional calculations, we have to distinguish two different cases.
Rather good results are still found for RPBE@PBE and EV93@PBE, which is probably due to the fact that the functional form of RPBE and EV93 does not differ too much from PBE, so that the PBE pseudopotentials are still accurate for RPBE and EV93.The increase in the mean errors is the largest for EV93@PBE with ABINIT, 0.06 eV for the MAE and 4.51% for the MAPE.
The situation is different for functionals that give accurate band gaps (HLE16, SLOC, and mBJLDA).Indeed, we can see a clear drop in accuracy, as seen by the increment in   the various mean errors.For instance, the MAPE lies between 4.5% and 6.5% with both ABINIT and VASP, which is two or three times larger than when LDA/PBE calculations are done with a consistent pseudopotential.Note that these are still relatively small values, clearly smaller than the average errors with respect to experiment, which are of the order of 30%. 51The loss in performance is more strikingly seen in terms of the dispersion of the errors (see Fig. 1).In particular, the standard deviation σ reaches values in the range 0.14 − 0.26 eV, which represents a threefold increase.This increase in dispersion leads inevitably to an increase in the maximum error, that can reach 1 eV or more as shown in Tables 2 and 3.For example, ABINIT:SLOC@LDA and ABINIT:HLE16@PBE yield errors of 1.33 and 1.10 eV for LaF 3 , respectively, corresponding to very large relative differences of ∼ 50% and ∼ 30% with respect to WIEN2k values.VASP-PAW setups do not perform better.Maximum errors of −1.17 eV with VASP:HLE16@PBE (for NaCl) and −0.76 eV with VASP:SLOC@LDA (for BiF 3 ) are obtained.The situation is particularly worrying for mBJLDA@LDA, where we find a maximum deviation of −3.87 eV (for Ne) with respect to WIEN2k, while three other error values are above 1 eV.
The solution to avoid such inaccuracies when using a new functional is rather clear and straightforward, namely a pseudopotential has to be generated, optimized and tested specifically for the new functional.This is what we have done for HLE16 and SLOC.The results are shown in Table 4, where we list the MAE, standard deviation σ, and MAPE for all functional combinations of LDA, PBE, HLE16, and SLOC.As expected, using a LDA potential in a PBE calculation or vice-versa has relatively small effect on the overall quality of the band gaps.However, using an inconsistent potential in a SLOC or HLE16 calculation has a much larger effect (by a factor of at least three), confirming what was already discussed above.Since, grossly speaking, the difference between LDA and PBE results is smaller than the one between PBE/LDA and HLE16/SLOC, it is not entirely unexpected that an inconsistent pseudopotential for these functionals yields worse results.
Performing consistent calculations brings back the pseudopotential error to the normal range exchange potential via a density-dependent parameter c: where In Eq. ( 1), n is the electron density, τ is the kinetic-energy density and v BR x is the Becke-Roussel potential. 67α = −0.012and β = 1.023Bohr 1/2 are parameters which were fitted specifically for band gaps.Originally, mBJLDA was implemented in an all-electron code, meaning that the quantities in Eqs. ( 1) and ( 2) are defined with respect to the total density n of the system, comprising core and valence electrons.However, in the case of a pseudopotential code, such as ABINIT, one has to work with a pseudodensity n that differs significantly from the total density in the regions around the atoms (i.e., it is much smoother).As such, it is doubtful that all-electron and pseudopotential calculations of Eq. ( 2) should give the same result.In fact, at the very least, the coefficients α and β should be reoptimized for pseudodensities (see, e.g., Ref. 68 for such a procedure with ABINIT).
There is also some incongruousness within the VASP implementation of mBJLDA.Even if theoretically the PAW method allows to recover the true density of the system, VASP performs an additive separation of c.This is certainly better than ignoring the contribution of the core density to c, but it is a numerical approximation.In practice this leads to an underestimation of c (for most cases, see Fig. S1) and therefore of the band gap (since the band gaps increases in a monotonous way with c).In order to illustrate more concretely the effect of an inaccurate value of c on the band gap, we show in Table 5 the mBJLDA band gap obtained with VASP and WIEN2k for some of the materials where the discrepancy between the two codes is the largest (from 1 to 4 eV).However, if the VASP calculations are done by fixing the value of c to the one obtained from the WIEN2k calculation, then the band gap is clearly much closer to the WIEN2k value.For instance, for Ne the VASP error gets reduced from about 4 eV to less than 1 eV.This shows that a large portion of the error with VASP is due to an inaccurate calculation of c.
The definition of the c parameter also creates some complications for the generation of a mBJLDA pseudopotential.As seen in Eq. ( 2), this parameter becomes ill-defined for (semi-)finite systems.This problem can however be bypassed.For example, Bartók and Yates 39 have recently proposed using a constant value for c during the pseudopotential generation, obtaining good agreement with all-electron calculations.Another possible solution is the use of the localized version of the mBJLDA potential recently proposed in Ref.

Conclusions
In summary, we performed a series of band gap calculations on a test set of 473 materials using all-electron, norm-conserving pseudopotential, and PAW methods.Our goal was to estimate the error on the band gap when standard LDA/PBE pseudopotentials or PAW setups are used inconsistently to perform electronic structure calculations using other density functionals.From our results we concluded the following.
1.As expected, consistent pseudopotential calculations are perfectly suited for the evaluation of band gaps.The errors with respect to the all-electron reference results (MAE ∼ 0.02 eV) are considerably smaller, by one or two orders of magnitude, than the errors (with respect to experiment) due to the approximation to the xc functional.
2. Using an inconsistent pseudopotential, i.e., one that was generated and optimized for another xc functional, increases the mean absolute errors by a factor of three or more, so that the MAE, for instance, can reach values around 0.1 eV.However, this is still smaller than the error due to other theoretical approximations, such as the choice of the xc functional, therefore the pseudopotential approach is still justified for several applications, especially when one is interested in average quantities.
3. Nevertheless, the error in few specific (and unpredictable) cases can be quite large, with band gaps sometimes differing from the all-electron results by several eV.Therefore, when one requires precise numerical estimations of band gaps, the use of consistent pseudopotentials or all-electron calculations is strongly recommended.
Of course, these conclusions can be relevant not only for semi-local functionals: the use of inconsistent pseudopotentials with hybrid functionals, LDA+U , or many-body approaches like GW should also be investigated more thoroughly.
The test set employed here did not contain magnetic materials.Effective-core methods are known to perform sometimes badly for this type of systems (see for example Ref. 69 ).
Preliminary calculations indicate that using the wrong pseudopotential in this situation can lead to errors as bad or worse than the ones observed here, but future tests are necessary to be able to properly quantify this effect.
As a final word, band gaps are only one of many quantities of interest for the solidstate community.Large scale studies of the effect of cross-functional calculations on other properties such as formation energies, geometries, absorption spectra, etc., are important and should be encouraged.Only with access to such quantitative data we can make informed decisions on the choice of methods for the calculation of electronic properties of materials.
• Spreadsheet containing the dataset and the calculated gaps.
• Tables SI-SVII: summary of the relevant statistical quantities for the different error distributions described in the main text.
• Figure SI: comparison of c values from mBJ obtained using WIEN2k and VASP.

n
i=1 |y i − y i,exp |/n; the mean error, ME = n i=1 (y i − y i,exp )/n; the standard deviation of the errors, σ = n i=1 (y i − y i,exp − ME) 2 /n; the median error (MnE); the interquartile range (IQR); the median of the absolute deviations from the median (MADM); the mean absolute percentage error, MAPE = 100 × n i=1 |y i − y i,exp |/(n y i,exp ); the mean percentage error MPE = 100 × n

Figure 1 :
Figure 1: Histograms of errors (with respect to the all-electron results) for band gaps computed with VASP (left) and ABINIT (right).Boxes have a width of 0.02 eV.

Table 1 :
M(A)E (in eV), standard deviation (in eV), M(A)PE (in %, for band gaps larger than 0.25 eV), and maximum absolute error (in eV) with respect to the all-electron results obtained with the considered functionals when a standard (LDA or PBE) pseudopotential is used.

Table 2 :
List of the ten materials with highest absolute error (with respect to the allelectron values) in the band gap (in eV) obtained with the ABINIT code for the SLOC and HLE16 functionals.In parentheses is indicated the all-electron band gap of the corresponding material (in eV).La 2 Ti 3 O 10 0.60(1.17)

Table 3 :
List of the ten materials with highest absolute error (with respect to the allelectron values) in the band gap (in eV) obtained with the VASP code for the SLOC, HLE16 and mBJLDA functionals.In parentheses is indicated the all-electron band gap of the corresponding material (in eV).

Table 4 :
MAE (in eV), standard deviation σ (in eV) and MAPE (in %) with respect to the all-electron results for the band gaps obtained using the ABINIT code for all crossfunctionals pseudopotential (PP) combinations.
66elow 0.05 eV for the MAE and σ, and at ∼ 2% for the MAPE).Unfortunately, generating consistent pseudopotentials is far from obvious for more sophisticated methods like hybrid functionals, many-body theories, or even the meta-GGA mBJLDA potential.As the mBJLDA is currently the best performing semilocal functional for the prediction of band gaps, 51,57 a more detailed discussion is in order.The exchange component of mBJLDA (the correlation component consists of LDA) is essentially a rescaling of the Becke-Johnson66

Table 5 :
mBJLDA band gaps (in eV) obtained with VASP and WIEN2k.The second set of VASP results were obtained with the parameter c in Eq. (1) fixed to the value obtained from the WIEN2k calculation.

TABLE SI .
Statistical data regarding the dispersion of results for LDA calculations with respect

TABLE SI .
(cont.) Statistical data regarding the dispersion of results for LDA calculations with respect to WIEN2k reference values.

TABLE SII .
(cont.) Statistical data regarding the dispersion of results for PBE calculations with respect to WIEN2k reference values.

TABLE SIII .
Statistical data regarding the dispersion of results for HLE16 calculations with respect to WIEN2k reference values.

TABLE SIII .
(cont.) Statistical data regarding the dispersion of results for HLE16 calculations with respect to WIEN2k reference values.

TABLE SIV .
Statistical data regarding the dispersion of results for SLOC calculations with respect

TABLE SIV .
(cont.) Statistical data regarding the dispersion of results for SLOC calculations with respect to WIEN2k reference values.

TABLE SV .
Statistical data regarding the dispersion of results for EV93 calculations with respectto WIEN2k reference values.

TABLE SVI .
Statistical data regarding the dispersion of results for RPBE calculations with respectto WIEN2k reference values.

TABLE SVII .
Statistical data regarding the dispersion of results for mBJLDA calculations with respect to WIEN2k reference values.Comparison of computed values of c (see Eq.(??)) for VASP and WIEN2k.