Band gaps of crystalline solids from Wannier-localization based optimal tuning of a screened range-separated hybrid functional
Dahvyd Wing, Guy Ohad, Jonah B. Haber, Marina R. Filip, Stephen E. Gant, Jeffrey B. Neaton, Leeor Kronik
BBand gaps of crystalline solids from Wannier-localization based optimal tuning of ascreened range-separated hybrid functional
Dahvyd Wing, Guy Ohad, Jonah B. Haber,
2, 3
Marina R. Filip, Stephen E. Gant,
2, 3
Jeffrey B. Neaton,
2, 3, 5 and Leeor Kronik Department of Materials and Interfaces, Weizmann Institute of Science, Rehovoth 76100, Israel Department of Physics, University of California, Berkeley, Berkeley, California 94720, USA Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA Department of Physics, University of Oxford, Oxford OX1 3PJ, United Kingdom Kavli Energy NanoSciences Institute at Berkeley, Berkeley, CA, 94720, USA (Dated: December 8, 2020)Accurate prediction of fundamental band gaps of crystalline solid state systems entirely withindensity functional theory is a long standing challenge. Here, we present a simple and inexpensivemethod that achieves this by means of non-empirical optimal tuning of the parameters of a screenedrange-separated hybrid functional. The tuning involves the enforcement of an ansatz that generalizesthe ionization potential theorem to the removal of an electron in an occupied state described by alocalized Wannier function in a modestly sized supercell calculation. The method is benchmarkedagainst experiment for a set of systems ranging from narrow band gap semiconductors to large bandgap insulators, spanning a range of fundamental band gaps from 0.2 to 14.2 eV and is found to yieldquantitative accuracy across the board, with a mean absolute error of ∼ ∼ The fundamental band gap of a semiconductor or in-sulator, defined as the difference between the ionizationpotential and the electron affinity of the material, is anessential material property. However, predicting it fromfirst principles using density functional theory (DFT) hasproven to be challenging [1, 2]. The Kohn-Sham (KS)lowest unoccupied - highest occupied eigenvalue gap can-not be equated with the fundamental gap even if the ex-act (and generally unknown) exchange-correlation func-tional is used [3, 4]. This is because the KS potential fea-tures a discontinuity as the number of electrons crossesinteger values [5], which results in a different reference po-tential for electron removal and addition, and thereforetypically causes KS eigenvalue differences to underesti-mate the fundamental band gap by as much as 50% [6–8]. In some cases, the discrepancy can be larger and evenlead to the spurious prediction of a metallic state [9, 10].For finite systems one can calculate accurate fundamen-tal band gaps from total energy differences between thecation, neutral, and anion systems [11]. However, forsolid state systems, the subject of this work, delocaliza-tion of the KS orbitals causes total energy differences toconverge to KS eigenvalue differences rather than to thefundamental band gap [12–16].Many DFT-based strategies for obtaining the funda-mental band gaps of solids by going beyond the KSscheme have been proposed over the years, e.g., Refs.16–40. However, two outstanding issues remain. One isachieving a level of accuracy that is on par with that ofexperiment (approximately 0.1 eV) for a wide range ofmaterials, from narrow band gap semiconductors to wideband gap insulators [41]. The other is doing so non-empirically within a formally exact framework, such thatother material properties, from defect energetics [42] tooptical absorption [43, 44], can be predicted on the samefooting. Here, we show that both these issues can be re- solved simultaneously. This is obtained by non-empirical,system specific, optimal parameter tuning of a screenedrange-separated hybrid (SRSH) functional, based on en-forcement of an ansatz generalizing the ionization poten-tial theorem to the removal of an electron from an orbitalcorresponding to a Wannier function.Our starting point is the screened range-separated hy-brid (SRSH) functional [45], which mixes a fraction ofexact exchange and semi-local exchange, as in a standardhybrid functional [46, 47], but with a generally differentfraction used in the short-range and in the long-range.This is accomplished by partitioning the exchange partof the Coulomb interaction using the identity [48]1 r = α erfc( γr ) r | {z } xx, SR + (1 − α ) erfc( γr ) r | {z } KSx, SR +1 (cid:15) ∞ erf( γr ) r | {z } xx, LR + (cid:18) − (cid:15) ∞ (cid:19) erf( γr ) r | {z } KSx, LR , (1)where the first and third terms are treated with exactexchange (xx) and the second and fourth terms are re-placed with (semi)local KS exchange (KSx). α is thenthe fraction of short-range (SR) exact exchange, (cid:15) ∞ thefraction of long-range (LR) exact exchange, and γ therange-separation parameter that determines the transi-tion from short to long range with increasing r , where r is the interelectron coordinate. The SRSH exchange-correlation functional is then given by E SRSHxc ( α, γ, (cid:15) ∞ ) = αE SR ,γ xx + (1 − α ) E SR ,γ KSx +1 (cid:15) ∞ E LR ,γ xx + (cid:18) − (cid:15) ∞ (cid:19) E LR ,γ KSx + E KSc , (2) a r X i v : . [ c ond - m a t . m t r l - s c i ] D ec where KSc denotes (semi)local KS correlation.The SRSH functional has several advantageous fea-tures. First, being a (range-separated) hybrid functional,it is a special case of the rigorous theoretical frameworkprovided by generalized Kohn-Sham (GKS) theory [49–52]. Therefore, self-consistent exchange-correlation po-tentials and kernels can be obtained in a straightforwardmanner by taking appropriate derivatives. Second, bysetting (cid:15) ∞ to the orientationally-averaged high-frequency(ion-clamped) dielectric constant, the functional pos-sesses the correct average long-range dielectric screening[23, 45]. This is known to be an important criterion forproducing accurate band gaps [23, 28, 29, 34, 37, 45, 51,53–56], optical absorption spectra [54, 57–61], and de-fect energy levels [53, 62]. Third, the short-range exactexchange fraction, α , can be chosen to balance exchangeand correlation effects and mitigate self-interaction errorsin the short range [63–65]. Here, we use semilocal ex-change and correlation components based on the Perdew-Burke-Ernzerhof (PBE) functional and adopt the defaultvalue of α = 0 .
25, as in the hybrid (PBE0) [47, 66] andshort-range hybrid (Heyd-Scuseria-Ernzerhof, HSE) [63]functionals.With values for α and (cid:15) ∞ defined uniquely, the onlyparameter left undetermined is the range-separation pa-rameter, γ . For finite systems, γ can be optimally tuned ,i.e., determined from first principles based on the satis-faction of physical criteria [51, 67]. This has typicallybeen achieved by selecting γ to satisfy the IP theorem,sometimes also known as the DFT version of Koopmans’theorem or as the generalized Koopmans’ theorem [5, 68–70]. It states that for the exact (G)KS functional (cid:15) ho = E ( N ) − E ( N − ≡ − I, (3)where (cid:15) ho is the highest occupied eigenvalue, E ( N ) isthe ground-state total energy of the system with N elec-trons, E ( N −
1) is the ground-state total energy of thesystem with one electron removed, and I is the ioniza-tion potential. For finite systems, deviation from the IPtheorem is equivalent to a missing derivative discontinu-ity, ∆ xc , in the approximate exchange-correlation func-tional (as compared to the exact functional)[71], suchthat (cid:15) ho + I ≈ ∆ xc [7]. Within the GKS framework, thepresence of a Fock-like operator “absorbs” some of thederivative discontinuity, such that each parameterizationof the Fock-like operator in a hybrid functional impliesa different exact remainder functional [49, 72]. Choos-ing γ to satisfy the IP theorem is therefore tantamountto selecting a functional form with a negligible deriva-tive discontinuity, meaning that the exact GKS eigen-value gap will equal the fundamental gap for that system[51, 67, 73].While the above method, known as the optimally-tuned range-separated hybrid (OT-RSH) functional ap-proach, has been very successful for determining the fun-damental gap in molecules (see, e.g., Refs. 51, 67, 74–78),it is not generally helpful for determining the band gapin solids. Due to the delocalized nature of the orbital corresponding to (cid:15) ho in solid state systems, the IP theo-rem is trivially satisfied for all parameterizations of theFock operator, regardless of whether the correspondingexact remainder functional has a derivative discontinu-ity or not [13–16]. Therefore, optimal tuning can not beapplied without further modification [79] and predictivepower is lost.The a priori selection of γ in the solid state has re-mained an open question and an active area of research[29, 34, 35, 37, 80]. Because it is the delocalization oforbitals that prevents the use of optimal tuning in thesolid state, the next logical step is to remedy the situa-tion by creating localized orbitals. Miceli et al. [35] haveachieved this by introducing a point defect in a supercellof a bulk solid and enforcing the IP theorem for the or-bital localized around the defect. While this approach iswell-justified physically and indeed significantly improvesthe predicted band gaps, it still suffers from sensitivityto the type of defect used for tuning [35, 36]. Further-more, ideally we would like to predict the properties ofthe pristine crystalline material without changing it.An alternate route for solving the problem of optimaltuning in the solid state is to rely on a different scheme tocreate orbital localization. Indeed, several recent strate-gies for band gap estimation have relied on localized or-bitals for obtaining correction terms that compensate forthe missing derivative discontinuity [21, 30, 31, 33, 35–38, 81]. Here, we exploit Wannier functions, which area localized orthonormal basis set obtained via a uni-tary transformation of a set of Bloch wavefunctions [82].While many different sets of Wannier functions can beproduced from a set of Bloch functions, maximally lo-calized Wannier functions [82] often match intuition forchemical bonds in solids and are “natural” localized or-bitals of the bulk system. However, Wannier functionsare not eigenfunctions of the (G)KS Hamiltonian [82];moreover, a system with N − N -electronground state with the charge density of a Wannier func-tion removed, is not the ground state of the N − ansatz that piecewise linearity of the total energy (whichis equivalent to the IP theorem [71]) be used in this casetoo, and this ansatz has been used to generate correctionterms for semilocal functionals which improved the accu-racy of the computed band gaps [30, 31, 38]. Here, weadopt this ansatz , but instead of using it to generate acorrection term, we use it to select the range-separationparameter, γ , in the SRSH functional of Eq. (1). To doso, we seek a value of γ that satisfies ∆ I γ = 0, where∆ I γ = E γ constr [ φ ]( N − − E γ ( N ) + h φ | ˆ H γ SRSH | φ i . (4)Here E γ constr [ φ ]( N −
1) is the total energy of a systemwith N − φ , is not occupied; ˆ H γ SRSH is theHamiltonian of the SRSH functional with N -electrons. h φ | ˆ H SRSH | φ i is the energy of the Wannier function,which we calculate via h φ | ˆ H SRSH | φ i = ∞ X i |h φ | ψ i i| (cid:15) i , (5)where ψ i and (cid:15) i are the GKS eigenfunctions and eigen-values of the N -electron system, and the sum is over allorbitals.To calculate the total energy of the N − E constr [ φ ]( N − f φ , that is˜ E constr [ φ ]( N −
1) = min ψ E [ { ψ } ]( N − λ N-1 X i |h ψ i | φ i| − f φ ! , (6)where { ψ } and E [ { ψ } ]( N −
1) are the eigenfunctions andtotal energy of the system with N − λ isa Lagrange multiplier. Taking the functional derivativeof Eq. (6) then yields the constrained GKS equation:ˆ H SRSH | ψ i i + λ | φ i h φ | ψ i i = (cid:15) i | ψ i i , (7)which we solve self-consistently for the N − λ determines the occupation of the Wannier func-tion, f φ , and, in this case we enforce f φ < × − bysetting λ = 15 Ry. We find that, in general, by set-ting λ > ∼ N − E constr [ φ ]( N −
1) using the Makov-Paynemonopole image charge correction for a charged system[86–88]: E constr [ φ ]( N −
1) = ˜ E constr [ φ ]( N −
1) + α mad q (cid:15) ∞ L , (8)where α mad = 2 .
837 is the Madelung constant for a sim-ple cubic cell; q is the charge of the system, e.g. q = 1;and L is the length of the supercell.Use of the above ingredients allows us to investigatea wide range of crystalline solids, from narrow bandgap semiconductors to wide band gap insulators, usingthe following four step procedure. In step 1, we cal-culate the orientationally-averaged high-frequency (ion-clamped) dielectric constant, (cid:15) ∞ , using a primitive unitcell and use this to set the fraction of long-range exactexchange in the SRSH functional [29, 57, 58]. In step2, we calculate the Wannier functions of a supercell ofan N -electron system. In step 3, we calculate ∆ I us-ing the SRSH functional with a particular value of γ . This entails calculating E ( N ) and E constr [ φ ]( N −
1) us-ing the Wannier function and the supercell from step 2.We iterate step 3 for different range-separation param-eters until we find the range-separation parameter thatyields | ∆ I | γ < .
02 eV. Here and throughout we referto the functional obtained from this choice of parametersas the Wannier-localized, optimally-tuned SRSH (WOT-SRSH) functional. In step 4, we calculate the fundamen-tal band gap for the material with this functional, usinga primitive unit cell. Importantly, no empirical fitting isintroduced in any step.For convenience we calculate (cid:15) ∞ in step 1 with theVienna ab initio simulation package ( VASP ) [89], a planewave code, using PBE-based projector-augmented waves(PAWs) for treating core electrons [90]. We use HSE tocalculate the dielectric constant for semiconductors [91],including spin-orbit coupling effects where necessary, andPBE0 to calculate the dieletric constant for insulators(i.e. GaN and materials with larger band gaps) [28]. Formore details, see the supplementary materials (SM) [92].We calculate steps 2 and 3 using supercells for threereasons. First, the larger the supercell/ k -point grid themore accurately the calculation approximates the bulksystem. Second, large k -grid calculations approach thebulk limit in which an infinitesimal fraction of an elec-tron is removed from each unit cell. This causes the IPtheorem to be trivially satisfied. We use Γ-point-onlysupercell calculations because they do not have this dif-ficulty. Third, the remaining image charge effects thatare not taken into account by the image charge correc-tion rapidly diminish with increasing supercell size. Weuse a 2 × × × × I calculations to ∼ .
02 eV and ∼ .
07 eV respectively (see SM for details[92]). When calculating the image charge correction viaEq. (8), we use the orientationally-averaged (cid:15) ∞ calcu-lated in step 1. For Wurtzite materials, we use L = √ Ωwhere Ω is the volume of the supercell. For these weaklyanisotropic systems, this produces an image charge cor-rection that is nearly equivalent to more exact methodsof incorporating anisotropy [93].We use the
Quantum Espresso [94] plane-wave codeto calculate steps 2-4. We use optimized norm-conserving(NC) Vanderbilt pseudopotentials [95] obtained from theonline repository, pseudo-dojo [96] (see SM for com-plete computational details [92]). Methods which includeexact exchange are known to be sensitive to the numberof semicore states [97]. We find that for Ge, Ga, In, Asand Sb it is important to include one complete shell ofsemicore states as valence electrons. Maximally localizedWannier functions are generated using the
Wannier90 software package [98]. We find that ∆ I is sensitive to thecharacter of the Wannier function used and that it is im-portant to choose the Wannier function that is composedof the top-most valence bands. However, we find that ∆ I is insensitive to the functional used to generate the Wan-nier function. We therefore use the PBE functional togenerate the Wannier function, except for narrow gapsemiconductors where PBE produces a spurious metallicstate, in which case we use PBE0 instead. For step 3,when α is close to (cid:15) ∞ , varying γ does not greatly changethe amount of exact exchange in the functional and thesearch for an optimal γ may fail. In this case, the valueof α is slightly increased beyond the default value of 0.25and the process of looking for the optimal γ is repeated.Our ∆ I calculations do not include spin-orbit couplingeffects. However, in step 4 we include spin-orbit couplingeffects where necessary. When calculating indirect bandgaps in step 4, we either use Wannier interpolation orsample the required specific k -point.To demonstrate how the above approach works in prac-tice, in Fig. 1(a) we plot ∆ I as a function of γ for AlP.Clearly, ∆ I varies monotonically with γ , such that thereexists a γ for which ∆ I = 0. In Fig. 1(b), we plot theband gap as a function of ∆ I and again find a monotonicdependence. The same behavior has been observed forall materials investigated in this work (see Table I for thelattice parameters and WOT-SRSH functional parame-ters used in this study).A summary of the fundamental band gaps predictedusing the WOT-SRSH approach, compared to experi-mental room-temperature fundamental band gaps [99],is given in Table I and in Fig. 2. Our calculations donot include electron-phonon coupling. Since zero pointrenormalization (ZPR) typically closes the band gap withrespect to the fixed ion approximation [109], we add theZPR energy to the experimental band gaps to form a“reference band gap” that electronic structure theoryshould ideally match. Additionally, we partially accountfor finite temperature effects by using room temperatureexperimental lattice parameters. It is readily observedthat excellent agreement between predicted and referenceband gaps is found throughout. Importantly, the meanabsolute error (MAE) is a satisfyingly small ∼ ∼ ± . I is ∼ ∼ ansatz for localized orbitals, in this case maximally-localized Wannier functions. Practically, it only requiresmodest supercell calculations with a hybrid functional.It may therefore serve as a useful means not only for theprediction of band gaps, but also for the non-empiricalprediction of other properties for which hybrid function-als are useful, such as optical absorption spectra anddefect energetics. It also lends itself as a natural andwell-motivated starting point for many-body perturba-tion theory calculations within the GW approximation. ACKNOWLEDGMENTS
This work was supported via a US-Israel National Sci-ence Foundation - Binational Science Foundation (NSF-BSF) grant, DMR-1708892, and by the Israel Ministryof Defense. Computational resources were provided bythe Texas Advanced Computing Center (TACC) and theNational Energy Research Scientific Computing Center,DOE Office of Science User Facilities supported by theOffice of Science of the U.S. Department of Energy underContract No. DE-AC02-05CH11231. [1] S. K¨ummel and L. Kronik, Rev. Mod. Phys. , 3(2008), see section IV.A and references therein.[2] G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. , 601 (2002).[3] J. P. Perdew and M. Levy, Phys. Rev. Lett. , 1884(1983).[4] L. J. Sham and M. Schl¨uter, Phys. Rev. Lett. , 1888(1983).[5] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz,Phys. Rev. Lett. , 1691 (1982).[6] R. W. Godby, M. Schl¨uter, and L. J. Sham, Phys. Rev.Lett. , 2415 (1986).[7] M. J. Allen and D. J. Tozer, Mol. Phys. , 433 (2002).[8] G. K.-L. Chan, J. Chem. Phys. , 4710 (1999).[9] S. Massidda, A. Continenza, A. J. Freeman, T. M.de Pascale, F. Meloni, and M. Serra, Phys. Rev. B , 12079 (1990).[10] N. A. Lima, L. N. Oliveira, and K. Capelle, EPL ,601 (2002).[11] D. J. Tozer and F. De Proft, J. Phys. Chem. A ,8923 (2005).[12] R. W. Godby and I. D. White, Phys. Rev. Lett. ,3161 (1998).[13] P. Mori-S´anchez, A. J. Cohen, and W. Yang, Phys. Rev.Lett. , 146401 (2008).[14] E. Kraisler and L. Kronik, J. Chem. Phys. , 18A540(2014).[15] V. Vlˇcek, H. R. Eisenberg, G. Steinle-Neumann, L. Kro-nik, and R. Baer, J. Chem. Phys. , 034107 (2015).[16] A. G¨orling, Phys. Rev. B , 245120 (2015).[17] D. M. Bylander and L. Kleinman, Phys. Rev. B , 7868(1990). γ (Bohr − ) ∆ I ( e V ) (a) ∆ I (eV) E S R S H g ( e V ) (b) FIG. 1. (a) Deviation from the IP theorem, ∆ I , for the SRSH functional, as a function of the range separation parameter, γ ,for AlP. (b) The fundamental band gap of AlP, calculated by SRSH, as a function of ∆ I .TABLE I. Parameters of WOT-SRSH calculations used in this work: experimental room-temperature lattice parameters ( a , c for Wurtzite structure); fraction of short-range exact exchange, α ; range-separation parameter, γ ; the calculated orientationally-averaged high-frequency dielectric constant; predicted fundamental band gap (* = including spin-orbit coupling effects); ref-erence band gap (sum of the experimental room temperature fundamental band gap [99] and the zero point renormalizationenergy). Also given are the mean absolute error (MAE) and the mean signed error (MSE), E WOT-SRSH g − E ref. g . a lat (˚A) α γ (Bohr − ) (cid:15) ∞ E WOT-SRSH g (eV) E ref. g (eV) ( E expt. g , ZPR) (eV)InSb 6.48 a a , 0.02 b )InAs 6.06 a a , 0.02 b )Ge 5.66 c c , 0.05 b )GaSb 6.10 a a , 0.03 b )Si 5.43 c c , 0.06 b )InP 5.87 a a , 0.05 b )GaAs 5.65 a a , 0.05 b )AlSb 6.14 a a , 0.04 b )AlAs 5.66 a a , 0.04 b )GaP 5.45 a a , 0.08 b )AlP 5.47 a a , 0.02 b )GaN 3.19, 5.19 c a , 0.17 b )C 3.57 c d , 0.38 e )AlN 3.11, 4.98 c a , 0.38 e )MgO 4.22 c f , 0.53 g )LiF 4.03 h i , 1.15 g )MAE (eV) 0.08MSE (eV) 0.01 a Ref. [100] b Ref. [101] Deduced by comparing the experimental 4K band gap to an extrapolated band gap c Ref. [102] d Ref. [103] e Ref. [104] DFPT calculation f Ref. [105] measured at 80K g Ref. [34, 106] GW calculation including both Migdal-Fan and Debye-Waller contributions. h Ref. [107] i Ref. [108] measured at 200-250K E ref .g (eV) E W O T − S R S H g ( e V ) G a N C A l N M g O L i F I n S b I n A s G e G a S b S i I n P G a A s A l S b A l A s G a P A l P FIG. 2. Fundamental band gaps predicted by WOT-SRSH, compared to the reference band gaps (fundamental experimentalband gaps plus zero-point renormalization energy). The straight line indicates perfect agreement. Inset: zoom-in on the 0 to3 eV region.[18] C. B. Geller, W. Wolf, S. Picozzi, A. Continenza,R. Asahi, W. Mannstadt, A. J. Freeman, and E. Wim-mer, Appl. Phys. Lett. , 368 (2001).[19] J. Heyd, J. E. Peralta, G. E. Scuseria, and R. L. Martin,J. Chem. Phys. , 174101 (2005).[20] M. Cococcioni and S. de Gironcoli, Phys. Rev. B ,035105 (2005).[21] V. I. Anisimov and A. V. Kozhevnikov, Phys. Rev. B , 075125 (2005).[22] L. G. Ferreira, M. Marques, and L. K. Teles, Phys. Rev.B , 125116 (2008).[23] T. Shimazaki and Y. Asai, Chem. Phys. Lett. , 91(2008).[24] Y. Zhao and D. G. Truhlar, J. Chem. Phys. , 074103(2009).[25] F. Tran and P. Blaha, Phys. Rev. Lett. , 226401(2009).[26] M. K. Y. Chan and G. Ceder, Phys. Rev. Lett. ,196403 (2010).[27] M. A. L. Marques, J. Vidal, M. J. T. Oliveira, L. Rein-ing, and S. Botti, Phys. Rev. B , 035119 (2011). [28] J. H. Skone, M. Govoni, and G. Galli, Phys. Rev. B ,195112 (2014).[29] J. H. Skone, M. Govoni, and G. Galli, Phys. Rev. B ,235106 (2016).[30] J. Ma and L.-W. Wang, Sci. Rep. , 1 (2016).[31] M. Weng, S. Li, J. Ma, J. Zheng, F. Pan, and L.-W.Wang, Appl. Phys. Lett. , 054101 (2017).[32] P. Verma and D. G. Truhlar, J. Phys. Chem. C ,7144 (2017).[33] N. L. Nguyen, N. Colonna, A. Ferretti, and N. Marzari,Phys. Rev. X , 021051 (2018).[34] W. Chen, G. Miceli, G.-M. Rignanese, andA. Pasquarello, Phys. Rev. Mater. , 073803 (2018).[35] G. Miceli, W. Chen, I. Reshetnyak, and A. Pasquarello,Phys. Rev. B , 121112 (2018).[36] T. Bischoff, I. Reshetnyak, and A. Pasquarello, Phys.Rev. B , 201114 (2019).[37] T. Bischoff, J. Wiktor, W. Chen, and A. Pasquarello,Phys. Rev. Mater. , 123802 (2019).[38] M. Weng, F. Pan, and L.-W. Wang, NPJ Comput.Mater. , 1 (2020). [39] L. A. Cipriano, G. Di Liberto, S. Tosoni, and G. Pac-chioni, J. Chem. Theory Comput. , 3786 (2020).[40] N. Tancogne-Dejean and A. Rubio, Phys. Rev. B ,155117 (2020).[41] P. Borlido, T. Aull, A. W. Huran, F. Tran, M. A. L.Marques, and S. Botti, J. Chem. Theory Comput. ,5069 (2019).[42] C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer,G. Kresse, A. Janotti, and C. G. Van de Walle, Rev.Mod. Phys. , 253 (2014).[43] N. T. Maitra, J. Chem. Phys. , 220901 (2016).[44] Y.-M. Byun, J. Sun, and C. A. Ullrich, Electron. Struct. , 023002 (2020).[45] S. Refaely-Abramson, S. Sharifzadeh, M. Jain, R. Baer,J. B. Neaton, and L. Kronik, Phys. Rev. B , 081204(2013).[46] A. D. Becke, J. Chem. Phys. , 1372 (1993).[47] J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem.Phys. , 9982 (1996).[48] This is equivalent to the more compact representationin previous papers [57], where the xx and KSx terms arecombined and α + β = 1 /(cid:15) .[49] A. Seidl, A. G¨orling, P. Vogl, J. Majewski, and M. Levy,Phys. Rev. B , 3764 (1996).[50] A. G¨orling and M. Levy, J. Chem. Phys. , 2675(1997).[51] L. Kronik, T. Stein, S. Refaely-Abramson, and R. Baer,J. Chem. Theory Comp. , 1515 (2012).[52] R. Baer and L. Kronik, Eur. Phys. J. B , 1 (2018).[53] D. Wing, J. Strand, T. Durrant, A. L. Shluger, andL. Kronik, Phys. Rev. Mater. , 083808 (2020).[54] A. K. Manna, S. Refaely-Abramson, A. M. Reilly,A. Tkatchenko, J. B. Neaton, and L. Kronik, J. Chem.Theory Comput. , 2919 (2018).[55] L. Kronik and S. K¨ummel, Adv. Mater. , 1706560(2018).[56] A. Tal, P. Liu, G. Kresse, and A. Pasquarello, Phys.Rev. Res. , 032019 (2020).[57] S. Refaely-Abramson, M. Jain, S. Sharifzadeh, J. B.Neaton, and L. Kronik, Phys. Rev. B , 081204 (2015).[58] D. Wing, J. B. Haber, R. Noff, B. Barker, D. A. Egger,A. Ramasubramaniam, S. G. Louie, J. B. Neaton, andL. Kronik, Phys. Rev. Mater. , 064603 (2019).[59] J. Sun, J. Yang, and C. A. Ullrich, Phys. Rev. Res. ,013091 (2020).[60] J. Sun and C. A. Ullrich, Phys. Rev. Mater. , 095402(2020).[61] D. Wing, J. B. Neaton, and L. Kronik, Adv. TheorySimul. n/a , 2000220 (2020).[62] M. Gerosa, C. E. Bottani, L. Caramella, G. Onida,C. Di Valentin, and G. Pacchioni, J. Chem. Phys. ,134702 (2015).[63] J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem.Phys. , 8207 (2003); J. Heyd, ibid . , 219906(2006).[64] D. L¨uftner, S. Refaely-Abramson, M. Pachler, R. Resel,M. G. Ramsey, L. Kronik, and P. Puschnig, Phys. Rev.B , 075204 (2014).[65] D. A. Egger, S. Weissman, S. Refaely-Abramson,S. Sharifzadeh, M. Dauth, R. Baer, S. K¨ummel, J. B.Neaton, E. Zojer, and L. Kronik, J. Chem. TheoryComp. , 1934 (2014).[66] C. Adamo and V. Barone, J. Chem. Phys. , 6158(1999). [67] T. Stein, H. Eisenberg, L. Kronik, and R. Baer, Phys.Rev. Lett. , 266802 (2010).[68] C.-O. Almbladh and U. von Barth, Phys. Rev. B ,3231 (1985).[69] J. P. Perdew and M. Levy, Phys. Rev. B , 16021(1997).[70] M. Levy, J. P. Perdew, and V. Sahni, Phys. Rev. A ,2745 (1984).[71] T. Stein, J. Autschbach, N. Govind, L. Kronik, andR. Baer, J. Phys. Chem. Lett. , 3740 (2012).[72] R. Garrick, A. Natan, T. Gould, and L. Kronik, Phys.Rev. X , 021040 (2020).[73] J. P. Perdew, W. Yang, K. Burke, Z. Yang, E. K. U.Gross, M. Scheffler, G. E. Scuseria, T. M. Henderson,I. Y. Zhang, A. Ruzsinszky, H. Peng, J. Sun, E. Trushin,and A. G¨orling, Proc. Nat. Acad. Sci. , 2801 (2017).[74] S. Refaely-Abramson, R. Baer, and L. Kronik, Phys.Rev. B , 075144 (2011).[75] J. Autschbach and M. Srebro, Acc. Chem. Res. , 2592(2014).[76] H. Phillips, Z. Zheng, E. Geva, and B. D. Dunietz, Org.Electron. , 1509 (2014).[77] M. E. Foster, J. D. Azoulay, B. M. Wong, and M. D.Allendorf, Chem. Sci. , 2081 (2014).[78] T. K¨orzd¨orfer and J.-L. Br´edas, Acc. Chem. Res. ,3284 (2014).[79] Except in the special case of a molecular solid, whereit could be “inherited” from the underlying molecule[45, 54].[80] I. C. Gerber, J. G. ´Angy´an, M. Marsman, and G. Kresse,J. Chem. Phys. , 054101 (2007).[81] C. Li, X. Zheng, N. Q. Su, and W. Yang, Natl. Sci. Rev. , 203 (2017).[82] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, andD. Vanderbilt, Rev. Mod. Phys. , 1419 (2012).[83] P. H. Dederichs, S. Bl¨ugel, R. Zeller, and H. Akai, Phys.Rev. Lett. , 2512 (1984).[84] K. Nakamura, R. Arita, Y. Yoshimoto, andS. Tsuneyuki, Phys. Rev. B , 235113 (2006).[85] M. S. Hybertsen, M. Schl¨uter, and N. E. Christensen,Phys. Rev. B , 9028 (1989).[86] G. Makov and M. C. Payne, Phys. Rev. B , 4014(1995).[87] M. Leslie and N. J. Gillan, J. Phys. C: Solid State Phys. , 973 (1985).[88] H.-P. Komsa, T. T. Rantala, and A. Pasquarello, Phys.Rev. B , 045112 (2012).[89] G. Kresse and J. Furthm¨uller, Phys. Rev. B , 11169(1996).[90] G. Kresse and D. Joubert, Phys. Rev. B , 1758 (1999).[91] J. Paier, M. Marsman, and G. Kresse, Phys. Rev. B ,121201 (2008).[92] See Supplemental Material at [URL will be inserted bypublisher] for more computational details and a discus-sion of the convergence of ∆ I calculations with respectto supercell size.[93] R. Rurali and X. Cartoix´a , Nano Lett. , 975 (2009).[94] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau,M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni,D. Ceresoli, M. Cococcioni, N. Colonna, I. Carn-imeo, A. D. Corso, S. de Gironcoli, P. Delugas,R. A. D. Jr, A. Ferretti, A. Floris, G. Fratesi,G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj,E. K¨u¸c¨ukbenli, M. Lazzeri, M. Marsili, N. Marzari,F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-laRoza, L. Paulatto, S. Ponc´e, D. Rocca, R. Sabatini,B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov,I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu,and S. Baroni, J. Phys.: Condens. Matter , 465901(2017).[95] D. R. Hamann, Phys. Rev. B , 085117 (2013).[96] M. van Setten, M. Giantomassi, E. Bousquet, M. Ver-straete, D. Hamann, X. Gonze, and G.-M. Rignanese,Comput. Phys. Commun. , 39 (2018).[97] M. L. Tiago, S. Ismail-Beigi, and S. G. Louie, Phys.Rev. B , 125212 (2004).[98] A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza,D. Vanderbilt, and N. Marzari, Comput. Phys. Com-mun. , 2309 (2014).[99] The experimental fundamental band gaps reported areproduced by adding estimated or calculated excitonbinding energies to the optical absorption edge or byinferring the fundamental band gap position based onthe location and identification of excitonic absorptionpeaks. [100] I. Vurgaftman, J. R. Meyer, and L. R. Ram-Mohan, J.appl. phys. , 5815 (2001).[101] M. Cardona and M. L. W. Thewalt, Rev. Mod. Phys. , 1173 (2005).[102] O. Madelung, Semiconductors Data Handbook , 3rd ed.(Springer-Verlag Berlin Heidelberg, 2004).[103] C. D. Clark, P. J. Dean, P. V. Harris, and W. C. Price,Proc. Math. Phys. Eng. Sci. , 312 (1964).[104] S. Ponc´e, Y. Gillet, J. Laflamme Janssen, A. Marini,M. Verstraete, and X. Gonze, J. Chem. Phys. ,102813 (2015).[105] R. Whited, C. J. Flaten, and W. Walker, Solid StateCommun. , 1903 (1973).[106] J. P. Nery, P. B. Allen, G. Antonius, L. Reining,A. Miglio, and X. Gonze, Phys. Rev. B , 115145(2018).[107] R. C. Leckey, 2.2.2 lithium halides: Datasheet fromlandolt-b¨ornstein - group iii condensed matter · volume23a: “subvolume a” in springermaterials (1989).[108] M. Piacentini, Solid State Commun. , 697 (1975).[109] F. Giustino, S. G. Louie, and M. L. Cohen, Phys. Rev.Lett. , 265501 (2010). upplementary Material to: Band gaps of crystalline solids from Wannier-localizationbased optimal tuning of a screened range-separated hybrid functional Dahvyd Wing, Guy Ohad, Jonah B. Haber,
2, 3
Marina R. Filip, Stephen E. Gant,
2, 3
Jeffrey B. Neaton,
2, 3, 5 and Leeor Kronik Department of Materials and Interfaces, Weizmann Institute of Science, Rehovoth 76100, Israel Department of Physics, University of California, Berkeley, Berkeley, California 94720, USA Materials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA Department of Physics, University of Oxford, Oxford OX1 3PJ, United Kingdom Kavli Energy NanoSciences Institute at Berkeley, Berkeley, CA, 94720, USA (Dated: December 8, 2020)
I. COMPUTATIONAL DETAILS
In table I we show the numerical parameters used for each step of obtaining the WOT-SRSH functional, as describedin the main text. For step 1 we calculate the ion-clamped, high-frequency dielectric constant, (cid:15) ∞ , by the change inpolarization in response to a finite electric field [1, 2]. In these calculations local field effects are included for bothHartree and exchange-correlation potentials [3]. The Ga, Ge, and In PAWs include semicore d -electrons as valenceones. The total energy convergence criteria for calculations where the small electric field is applied is 10 − eV. Mostof the dielectric constants used for this study were already reported in Refs. [4, 5].For steps 2-4, we use version 0.4 pseudopotentials from pseudo-dojo [6]. For most pseudopotentials we usestringent fully relativistic pseudopotentials. We use standard scalar relativistic pseudopotentials for Si, AlP, AlN, C,MgO, and LiF. TABLE I. Plane wave energy cutoff , E cutoff , and k -grid size for calculating the dielectric constant, (cid:15) ∞ ; deviation from satisfactionof the generalized ionization potential theorem, ∆ I ; and the fundamental band gap, E WOT-SRSH g (cid:15) ∞ ∆ I E g E cutoff (eV) k -grid E cutoff (Ry) E cutoff (Ry) k -gridInSb 300 8 × × × × × × × × × × × × × × × × × ×
11 20 30 10 × × × ×
11 60 70 8 × × × ×
11 100 110 8 × × × ×
11 75 85 8 × × × ×
11 95 105 8 × × × ×
11 90 100 8 × × × ×
11 30 40 8 × × × × × ×
6C 450 9 × × × × × × × × × × × × × × × × II. CONVERGENCE OF ∆ I WITH RESPECT TO SUPERCELL SIZE
For seven materials, we calculated ∆ I for different supercell sizes, while keeping the SRSH parameters fixed foreach system, as summarized in Table II. The results for Si, AlP, InSb, MgO and LiF show that ∆ I for a 2 × × a r X i v : . [ c ond - m a t . m t r l - s c i ] D ec × ×
2. Based on the convergence of the other materials we estimate that ∆ I is converged to 0.1 eV for AlN andGaN.The Makov-Payne image charge correction has been shown to accelerate convergence of ∆ I calculations with respectto the supercell size [7]. We use the same value of the dielectric constant in the Makov-Payne correction in Eq. (4)as we use for setting the fraction of the screened long-range exact exchange in the SRSH functional. This dielectricconstant is converged with respect to the k -grid and includes spin-orbit coupling effects in order to be consistentwith the calculation of the band gap of the material in step 4. We find that using Makov-Payne corrections withdielectric constants calculated without spin-orbit effects and with a k -grid that matches the supercell used in the˜ E γ constr [ φ ]( N −
1) calculations causes ∆ I calculations to converge somewhat slowly with respect to the supercell size,while ultimately yielding the same result. TABLE II. ∆ I (eV) as a function of supercell size for selected materialsSupercell size Si AlP InSb MgO LiF1 × × × × × × × × × × × × , 155107 (2001).[2] I. Souza, J. ´I˜niguez, and D. Vanderbilt, Phys. Rev. Lett. , 117602 (2002).[3] J. E. Northrup, M. S. Hybertsen, and S. G. Louie, Phys. Rev. Lett. , 819 (1987).[4] D. Wing, J. B. Haber, R. Noff, B. Barker, D. A. Egger, A. Ramasubramaniam, S. G. Louie, J. B. Neaton, and L. Kronik,Phys. Rev. Mater. , 064603 (2019).[5] D. Wing, J. B. Neaton, and L. Kronik, Adv. Theory Simul. n/a , 2000220 (2020).[6] M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze, and G.-M. Rignanese, Comput. Phys.Commun. , 39 (2018).[7] D. Wing, J. Strand, T. Durrant, A. L. Shluger, and L. Kronik, Phys. Rev. Mater.4