Accurate and efficient band-gap predictions for metal halide perovskites at finite temperature

We develop a computationally efficient scheme to accurately determine finite-temperature band gaps. We here focus on materials belonging to the class ABX3 (A = Rb, Cs; B = Ge, Sn, Pb; and X = F, Cl, Br, I), which includes halide perovskites. First, an initial estimate of the band gap is provided for the ideal crystalline structure through the use of a range-separated hybrid functional, in which the parameters are determined nonempirically from the electron density and the high-frequency dielectric constant. Next, we consider two kinds of band-gap corrections to account for spin-orbit coupling and thermal vibrations including zero-point motions. In particular, the latter effect is accounted for through the special displacement method, which consists in using a single distorted configuration obtained from the vibrational frequencies and eigenmodes, thereby avoiding lengthy molecular dynamics. The sequential consideration of both corrections systematically improves the band gaps, reaching a mean absolute error of 0.17 eV with respect to experimental values. The computational efficiency of our scheme stems from the fact that only a single calculation at the hybrid-functional level is required and that it is sufficient to evaluate the corrections at the semilocal level of theory. Our scheme is particularly convenient for large-size systems and for the screening of large databases of materials.


Abstract
We develop a computationally efficient scheme to accurately determine finite-temperature band gaps for metal halide perovskites belonging to the class ABX3 (A = Rb, Cs; B = Ge, Sn, Pb; and X = F, Cl, Br, I). First, an initial estimate of the band gap is provided for the ideal crystalline structure through the use of a range-separated hybrid functional, in which the parameters are determined nonempirically from the electron density and the high-frequency dielectric constant. Next, we consider two kinds of band-gap corrections to account for spinorbit coupling and thermal vibrations including zero-point motions. In particular, the latter effect is accounted for through the special displacement method, which consists in using a single distorted configuration obtained from the vibrational frequencies and eigenmodes, thereby avoiding lengthy molecular dynamics. The sequential consideration of both corrections systematically improves the band gaps, reaching a mean absolute error of 0.17 eV with respect to experimental values. The computational efficiency of our scheme stems from the fact that only a single calculation at the hybrid-functional level is required and that it is sufficient to evaluate the corrections at the semilocal level of theory. Our scheme is thus convenient for the screening of large databases of metal halide perovskites, including large-size systems.

Introduction
Metal halide perovskite solar cells have shown remarkable progress in power conversion efficiency, which has been boosted up to 25% within only a few years [1]. Approximately 2000 perovskites [2] can be synthesized but their suitability for photoelectric devices strongly depends on their electronic band gap. In particular, this property varies with the specific phase in which these materials occur, going from the highly symmetric cubic phase at high temperatures to structures of lower symmetry at lower temperatures. Given the large variety of available perovskites, it is important to develop a computational tool to screen large databases in search of an optimal material. Such a tool should provide band gaps at finite temperature in a computationally efficient way without compromising on accuracy. In the case of cubic halide perovskites, which have small unit cells, Wiktor et al. have shown that it is possible to achieve an accurate description of the experimental band gaps by combining, on the one hand, high-level many-body perturbation theory for the electronic structure and, on the other hand, ab initio molecular dynamics to account for the finite temperature, the two schemes implying extensive computational resources of comparable order. [3] Next, Bischoff et al. further demonstrated that the accuracy by which the electronic structure is described could be preserved through the use of computationally more advantageous nonempirical hybrid functionals. [4] However, this development needs to be accompanied with a similar reduction of computational cost for treating the effects of finite temperature in order to deploy such methodology to a large set of perovskite materials in an efficient way.
Density functional theory [5,6] based on semilocal approximations for the exchange and correlation energy represents a powerful tool for predicting ground state properties in materials sciences. However, it is well-known that such semilocal approximations generally underestimate the electronic band gap [7]. In contrast with this general notion, several previous studies focusing on perovskites have reported that the semilocal functional proposed by Perdew, Burke, and Ernzerhof (PBE) gives band gaps in good agreement with experimental data. [8][9][10] This agreement should be considered accidental because of the neglect of various effects [3], such as spin-orbit coupling (SOC), nuclear quantum effects (NQEs), and thermal vibrations. Many-body perturbation theory in the GW approximation [11] is currently recognized as the most accurate scheme for band-gap calculations, especially in its selfconsistent quasiparticle formulation including vertex corrections in the screening [12][13][14][15]. However, such calculations are computationally highly demanding and cannot be envisaged as an efficient tool for screening large databases.
Hybrid functionals like the global PBE0 [16,17] or the range-separated HSE [18,19] are obtained by admixing a fraction of Fock exchange to the semilocal exchange potential and can serve as valuable alternatives to GW calculations. These hybrid functionals are defined through fixed mixing parameters, and thereby generally fail in providing accurate band gaps for a large class of materials. [20][21][22] It has been shown that the higher accuracy can be achieved through material-specific dielectric-dependent hybrid (DDH) functionals. [14,[23][24][25][26][27][28] Alkauskas et al. [23] and Marques et al. [24] linked the incorporated amount of Fock exchange to the inverse high-frequency dielectric constant 1/ǫ ∞ . Next, range-separated DDH functionals were developed in which 1/ǫ ∞ determines the fraction of Fock exchange in the long range (LR), while retaining a suitable fraction of Fock exchange in the short range (SR). [25,29,30] Very recently, further developments led to the suggestion of two kinds of nonempirical range-separated DDH functionals. [26,27] In the first one (denoted DD-RSH-CAM) the parameters are determined through the first-principles calculation of the dielectric screening function, [26] while in the second one (denoted DSH) the parameters are obtained through combining metallic and dielectric screening.
[27] These two advanced DDH schemes have been shown to yield band gaps of comparable accuracy to state-of-the-art GW calculations for a large variety of materials, [26,27] but at a significantly lower computational cost.
To achieve high accuracy at finite temperature, it is further necessary to account for electron-phonon interactions, which have been found to significantly affect the band gap through the zero-point motions and the thermal vibrations. [3] Widely used supercell methods to estimate the electronphonon interaction are based on the adiabatic approximation and include statistical sampling through Monte-Carlo [31,32] and molecular dynamics methods. [3,33,34] Within the context of the adiabatic approximation, one also finds the so-called special displacement method (SDM), in which only one single optimal configuration of the atomic positions is sufficient to evaluate the band-gap renormalization due to NQEs and thermal vibrations. [35,36] For polar materials, Poncé et al. pointed out that the adiabatic approximation leads to a divergence at small momentum transfer, [37] but Zacharias and Giustino remarked that the observed divergency is not intrinsic to the adiabatic approximation. [36] In the SDM, the divergence is found not to build up, [35,36,38,39] because the periodic boundary conditions effectively short-circuit the long-range electric field associated with longitudinal optical phonons. [35] Recent reports show that it is necessary to resort to the nonadiabatic formulation in the evaluation of the zero-point phonon renormalization in the case of infrared-active materials. [37,[40][41][42] Nonadiabatic effects are usually included through a generalized Fröhlich model, [40,42] which in its one-band formulation can also be used for an estimation of the correction. [39,41,43] Furthermore, the challenge of calculating band gaps of metal halide perovskites is exacerbated by other issues associated with thermal vibrations, such as anharmonic effects [44], strong dependence of spin-orbit coupling [3] and phonon renormalization [38,41,45] on the adopted electronic-structure theory, and cross-coupling between spin-orbit coupling and phonon renormalization. [3,46] In this work, we achieve accurate band-gap predictions for metal halide perovskites at finite temperature in a highly efficient way by combining dielectric-dependent hybrid functionals for the electronic structure with the special displacement method for the nuclear quantum effects and the thermal vibrations. To highlight the advantage of using dielectric-dependent hybrid functionals, our investigation also includes standard hybrid functionals, such as PBE0 and HSE06. Our best band-gap predictions are obtained by applying the following procedure. We start by calculating initial band-gap estimates with the dielectric dependent hybrid functional for the pristine crystalline structure. Next, two corrections are applied. The first correction results from spin-orbit coupling effects. The second correction accounts for the nuclear quantum effects and the thermal vibrations. In this case, the special displacement method is applied and nonadiabatic effects are estimated through the one-band Fröhlich model. The band gap of the supercell structure with the special displacements is extracted through a procedure similar to the Tauc analysis of experimental spectra. [47] We demonstrate that both corrections can be determined at the semilocal level of theory, which does not lead to any significant loss of accuracy with respect to hybrid-functional results of higher level and thus only implies a moderate computational effort. The sequential inclusion of the two corrections leads to band gaps that progressively improve with respect to available experimental references, thereby confirming the validity of our approach.
This manuscript is organized as follows. In Sec. 2, we first describe the model structures and the calculation of the dielectric constants. We then report our band-gap predictions focusing on the corrections due to spin-orbit coupling and to electron-phonon interactions. More discussions are drawn in Sec. 3. In Sec. 4, we briefly outline the main theoretical methodologies used in this work, i.e. the construction of the dielectric-dependent hybrid functional and the special displacement method. We also presented the calculation details.

Structural models
Our general motivation is to predict band gaps at room temperature for metal halide pervoskites. When various structures of a given material exist at room temperature, we take under consideration the available competitive structures to demonstrate the reliability of our approach upon structural variation. For benchmark purposes, we also consider materials at higher temperatures when a detailed experimental characterization is available.
The materials under consideration in this work all have the structural formula ABX 3 , where A and B are metal atoms and X is a halide atom. These materials come in different structures identified by the prefix β, γ, δ, R, and M, which indicate that their phase is tetragonal, orthorhombic with cornersharing octahedra, orthorhombic with edge-sharing octahedra, rhombohedral, and monoclinic, respectively. Within this group, all the structures with cornersharing BX 6 octahedra (R, γ, and β) are metal halide perovskites, while the other structures (δ and M) do not formally belong to this class. Representative atomic structures covering all space groups considered are illustrated in the Supplementary Fig. 1.
In this study, we do not consider cubic phases because they are generally unstable at 300 K, the targeted temperature for photovoltaic materials. The instability results from soft phonon modes, [46,48] which can be related to phase transitions. [48] This indicates that an harmonic description of the phonons would be inappropriate. Indeed, cubic structures are generally subject to important anharmonic effects, [44,49] which cause the excess free energy surface to change as a function of temperature. [49] The cubic phases at high temperature correspond to average structures, which differ noticeably from the local atomic scale description. [3,49] For instance, the average B-X-B bond angle corresponds to 180 • for the cubic structure, but in molecular dynamics simulations the distribution is peaked at a lower angle (∼168 • , Ref. 3). These structural changes depend on temperature and have noticeable effects on the band gap. Furthermore, the cross coupling between SOC and thermal effect has been found to be sizable in cubic systems. [3,46] However, our test calculations on noncubic δ-CsPbI 3 and R-CsPbF 3 show that the error associated to such cross-coupling effects is ∼0.1 eV (see the Supplementary Table 1). For these reasons, cubic structures would require a different theoretical treatment than the one proposed in this work.
Here, we make use of the experimentally determined structures given in the Supplementary Table 2. These structures are used in the calculations of the high-frequency dielectric constant and the band gap. Since these structures are taken from experiment, they already include the thermal expansion effect corresponding to the temperature at which their structure has been characterized (cf. discussion in Ref. 3).

Dielectric constants
The high-frequency dielectric constants used in the DSH functional are calculated through finite electric fields [50]. For fixed atomic positions, we calculate the variation of the polarization ∆P as a function of the electric field (cf. Supplementary Fig. 2 for δ-CsPbI 3 ). More specifically, ∆P = P ′ − P 0 , where P ′ and P 0 are the polarizations with and without the applied electric field. The dependence is close to linear up to a critical value, above which the calculation no longer converges [50]. For all materials, we use in this work a value of 0.001 a.u. for the electric field and determine the high-frequency dielectric constant through where V and E ǫ are the volume of the supercell and the electric field, respectively. The dielectric constants ǫ ∞ calculated at the PBE level are reported in Table 1. The values of the dielectric constants used in this work correspond to the average over the three principal directions of the dielectric tensor. In these perovskite materials, the dielectric constant in a given direction remains generally close to this average value (Supplementary Table 3). For instance, the material showing the largest spread is δ-RbGeI 3 , for which we found 4.93, 5.81 and 4.75, with an average of 5.16. Our calculated values show good agreement with previous PBE+U results from Ref. 51, which have been integrated as benchmarks in the Materials Project [52]. However, there are only few benchmarks for materials belonging to M and γ phases. Therefore, to further validate our calculation of dielectric constants, we take γ-RbGeI 3 and M-CsSnCl 3 as representatives of theses phases and calculate ǫ ∞ with the higher-level DSH functional that we derived. We obtain 5.09 and 3.15, which show differences of at most 7% compared to the plain PBE results of 5.10 and 3.38, respectively. This error is noticeably smaller than the typical mean absolute percentage error of 15% found for PBE results compared with experiment.
[26] Furthermore, the effect of this error on the band gap is to a large extent attenuated because ǫ ∞ enters through its inverse into the hybrid functional.
[26] It has also been shown that the DSH functional with PBE values for ǫ ∞ is generally more accurate for band-gap predictions than a DSH functional with self-consistently determined dielectric constants [27,53]. Hence, we take ǫ ∞ values calculated at the PBE level to derive the parameter α LR of the DSH functionals used in this work.
The inverse screening lengths µ of these materials vary between 0.67 and 0.82 bohr −1 , with an average value of 0.73 bohr −1 (cf.  Table 4 the µ values for a selection of materials as obtained with three methods, namely Thomas-Fermi, DD-RSH-CAM, and DSH. We observe that DSH and DD-RSH-CAM give close values of µ for most materials, which result in similar band-gap estimates. At variance, the Thomas-Fermi values of µ generally differ more significantly, particularly for wide band-gap materials. These deviations generally lead to inaccuracies as far as the calculated band gaps are concerned. As an example, the reader is referred to the case of LiF discussed in the Supplementary Table 4.

Band gaps
The band structures of the class of materials considered in this work, both perovskites and non-perovskites, have been heavily investigated in the literature through a variety of electronic structure methods going from semilocal to GW methods. [54][55][56][57][58] The calculated band gaps typically vary by as much as a factor of two depending on the method. In our study, we first determine the fundamental band gap of the ternary ABX 3 perovskites with PBE, PBE0(0.25), HSE06, and DSH (Table 2), including neither spin-orbit effects nor coupling to phonons. As reference, we also list in Table 2 band gaps calculated through high-level QSGW under the same assumptions. [54,59,60] As expected, the functional PBE strongly underestimates the band gaps showing the worst mean absolute error (MAE) of 1.30 eV and the worst mean absolute relative error (MARE) of 55%. The hybrid functional HSE06 opens the band gap with respect to the PBE value by 0.50 to 1.05 eV, and the hybrid functional PBE0(0.25) further opens the band gap by about 0.7 eV, in accord with general considerations [61]. We remark that the MAEs with respect to the accurate QSGW values progressively decrease when adopting in sequence PBE, HSE06, PBE0(0.25), and DSH (cf. Table 2 and Fig. 1). In particular, DSH yields a MAE of 0.22 eV and a MARE of 9%. To further support the quality of the DSH description, we perform for β-CsSnI 3 a comparison between the DSH band structure calculated here and the GW band structure obtained by Huang and Lambrecht. [54] As shown in the Supplementary Fig. 3, the band gap occurring at the Z point is well reproduced and the overall trends of the band structure are very similar. This level of agreement between hybrid functionals and GW matches that found for other materials in the literature. [62] We next analyze the general band-gap properties of the materials considered comparing variations of the phase. As shown in Table 2, the R, γ, and β phases all have direct band gaps, while the δ phases show an indirect band gap. For the monoclinic phases, the band gaps can be either direct or indirect. However, due to the flatness of the bands, [54] the difference between the direct and the indirect bands gaps is generally not large. For instance, in the case of M-CsSnF 3 , we find an indirect band gap of 6.59 eV to be compared with the smallest direct band gap of 6.94 eV. Additionally, we remark that the δ phases generally have larger band gaps than the corresponding β or γ phases. This can be readily understood by their band structures. To illustrate this property, we calculate the band structures of β-CsSnI 3 , γ-CsSnI 3 , and δ-CsSnI 3 with the PBE functional (cf. Supplementary Fig. 4). It is found that the valence bands of δ-CsSnI 3 are much flatter bands than those of β-CsSnI 3 and γ-CsSnI 3 , which can be related to the distance between nearest neighbor halogen atoms. This provides an explanation for the systematically wider band gaps found for the δ phase compared to the other phases.
Considering chemical variations within the same symmetry class, we first notice that the calculated band gap E bare for all classes decreases for variations of the halide atom going from Cl to I while keeping the composition in terms of metal atoms unmodified. For instance, E This property can be understood invoking the electronegativity, which reduces from Cl to I. Since the character of the valence band is dominated by the halide atoms, [54,63,64] such a downward move in the periodic table results in an upward shift of the valence band and consequently in a narrower band gap. Second, we consider band-gap variations when keeping the A and X composition unmodified and varying the B atom. In this case, the band gap of Pb-based compounds is always larger than the corresponding Sn-based counterpart. This effect can also be attributed to a reduction of electronegativity going from Pb (2.33) to Sn (1.96). In the absence of any change associated with the halide atom, the smaller electronegativity of Sn leads to an upward shift for both the VBM and CBM. [64] However, Tao et al. [64] found that the replacement of Pb with Sn results in a larger shift for the VBM than for the CBM because of their different ratios of s and p character. Therefore, smaller band gaps are found for Sn-based compounds. Since the electronegativity of Ge (2.01) falls in between those of Pb and Sn, we infer the following trend: . Third, when changing the A atom, we don't have a uniform rule for the band-gap variation. Indeed, we find E . This effect can be related to the fact that the A-atom states lie far away from the VBM and CBM, and thus do not influence their energy levels directly. [54,63,64] In the following, the calculated band gaps are benchmarked with respect to nine materials for which a careful experimental characterization of the band gap is available. This experimental set of ABX 3 materials comprises eight perovskites and one non-perovskite. The B site is occupied by Ge, Sn, and Pb atoms, while the X site is occupied by Cl, Br, and I atoms. The A site is solely occupied by Cs atoms, but this site is known not to affect to band gap in a significant way. [54,63,64] The band gaps covered by this experimental set of materials go from 1.25 to 4.46 eV, which is the relevant range for photovoltaic materials. Overall, these materials offer a diversified set to benchmark our method.
Many of the materials considered contain heavy elements, which are expected to give significant SOC effects. To estimate SOC corrections, we consider three different functionals: PBE, PBE0, and DSH. The band-gap variations are given in Table 3. We observe similar variations for the three functionals studied. The largest differences between PBE and DSH are found for R-CsPbF 3 and γ-CsPbF 3 and amount to only 0.10 and 0.09 eV. In a study comprising 103 materials but not including perovskites, Huhn et al. also found that the SOC corrections calculated with PBE and HSE06 are overall very similar, with most differences falling below 60 meV and the most notable one being 189 meV. [65] However, we remark that Wiktor et al. found that the SOC corrections on the band gap are underestimated by as much as 0.3 eV for Pb-based cubic phases, [3] indicating that SOC effects need careful evaluation for high-symmetry structures containing Pb atoms. Hence, we conclude that the PBE functional is sufficiently accurate to estimate SOC corrections for all considered materials in this work, none of which shows the cubic symmetry.
The calculated SOC corrections can be found in Fig. 2(a) and Table 4. These corrections always lead to a band-gap reduction with respect to E bare . The largest corrections correspond to Pb-based compounds and reach values around −1 eV for γ-CsPbCl 3 and γ-CsPbBr 3 . The largest correction for a Snbased compound amounts to ∼ −0.4 eV and is found for γ-CsSnI 3 . We remark that the SOC corrections do not only depend on composition but also on the underlying atomic structure. Indeed, in the case of CsPbBr 3 , we find −0.37 and −1.02 eV for the δ and γ structure, respectively. This indicates that SOC effects may be large and highly nontrivial.
We then obtain the band gaps E SOC by adding the SOC corrections to E bare (cf. Table 4). For the subset of materials for which experimental estimates are available, the accuracy of the corrected band gaps improves significantly, with the MAE reducing from 0.57 to 0.26 eV. Correspondingly, the MARE drops from 27% to 11%. These errors clearly illustrate the importance of considering SOC effects for achieving accurate band gaps.
In this section, we focus on the coupling to phonons, which affect the band gap through ZPR and thermal vibrations. For the ZPR and the vibrational properties, we use experimental lattice constants at the temperatures given in the Supplementary Table 2. In this way, thermal expansion effects are already accounted for. For room temperature, we take T = 300 K in this work. All the modes of the structures considered in this work show positive frequencies.
To obtain the band gap renormalization caused by zero-point motion and thermal vibrations, we use the special displacement method [35] described in Sec. 4.2. The supercell sizes adopted in this work contain between 160 and 625 atoms (cf. Supplementary Table 5) and have been chosen in order to achieve accuracies better than 0.05 eV. We checked this explicitly for all selected supercells by performing the same calculation with smaller supercells, leading to differences smaller than 0.05 eV for the band-gap predictions.
The special displacement method does not account for the long-range Fröhlich coupling resulting from the nonadiabatic treatment. [36,37,[39][40][41][42][43] To obtain an estimate for this contribution, we use the one-band Fröhlich model (see Supplementary Note 1). [43] For a series of perovskites belonging to the class of materials under investigation in our work, we obtain estimates ranging between −0.05 and −0.13 eV, with a mean value of −0.09 eV (see Supplementary Table 6). As will be seen in the following, the error resulting from the neglect of this effect is smaller than the final MAE claimed in our work (0.17 eV). Hence, it can be concluded that the long-range Fröhlich coupling coming from the nonadiabatic treatment of the electron-phonon interaction is not susceptible to affect in a significant way the accuracy of our results for metal halide perovskite materials.
The SDM method provides supercells with distorted atomic structures, which we use to extract the band gap. More specifically, these disordered supercells simulating the effect of finite temperature cannot be characterised with the same Brillouin zone as the ideal material, because of the lack of translational symmetry. To conform with the practice in the experimental determination of optical band gaps at finite temperature, [47] we therefore rely on extrapolations of the wings in the electronic density of states (EDOS) to obtain the VBM and CBM. For simplicity, we here perform a linear extrapolation of the band wings. [34,66] In order to perform reliable extrapolations of the band wings, the calculated eigenvalues in the vicinity of the band edges should be sufficiently dense.
This condition is not guaranteed for all materials when using the cp2k code, which solely uses the Γ point. In some cases, the states in vicinity of the band edge are sparsely distributed and do not connect with the band wing upon reasonable Gaussian smearing. It is then necessary to use denser k-point meshes or alternatively larger supercells. We illustrate this issue for an ideal structure of γ-CsSnI 3 in Fig. 3(a). First, we apply the code vasp [67,68] to the primitive unit cell and establish a converged benchmark by increasing the k-point sampling of the Brillouin zone (cf. Supplementary Fig. 5). Then, we compare the converged reference shown in Fig. 3(a) with a cp2k calculation obtained with a 6×4×6 supercell and find good agreement for the band gap within less than 0.005 eV. Hence, the band gap can be obtained accurately provided a sufficiently large supercell is used. In the Supplementary Table 5, we give suitable supercell sizes to obtain converged band gaps for the materials considered in this work. Once a sufficiently large supercell is ensured, the EDOS is determined for the ideal structure and for two displaced structures accounting for ZPR and ZPR+T , as illustrated for M-RbSnF 3 in Fig. 3(b). The linear extrapolation of the band wings then allows us to extract the VBM [ Fig. 3(c)] and the CBM [ Fig. 3(d)]. Such an analysis is performed for each material (cf. Supplementary Figs. 6 and 7).
The use of hybrid functionals for the EDOS with a dense k-point mesh or with a large supercell would undermine the efficiency of our band-gap determinations. There are reports in the literature with ZPR calculations showing a dependence on the adopted electronic-structure method. [38,41,45] Variations are particularly important for silicon [45] and diamond [38] and might amount up to ∼0.25 eV. For other materials like MgO and LiF, the variation between functionals is much smaller (∼0.05 eV). [41] Since it is difficult to develop a rationale for understanding the size of this dependence, we explicitly investigate whether thermal corrections for perovskite materials can be accurately evaluated with the semilocal PBE functional. For this, we select M-RbSnF 3 and M-CsSnF 3 , which show noticeable differences between the band gaps calculated with PBE0(0.25) and DSH (cf. Table 2). We check the effect of the functional by extracting the ZPR+T band-gap corrections from the EDOS obtained with PBE, PBE0(0.25) and PBE0(0.40) (see Table 5 and Supplementary Fig. 8). The functional PBE0(0.40) can be taken as representative of DSH for M-RbSnF 3 and M-CsSnF 3 , since the bare band gaps calculated with these two functionals differ by less than 0.08 eV. Additionally, we use PBE and PBE0(0.25) to test R-CsGeCl 3 , which shows the largest band gap renormalization due to ZPR+T in our work, i.e. of 0.4 eV (cf. Table 6). In all cases studied, we find that the corrections ∆E ZPR+T obtained with PBE coincide with those obtained with the hybrid functionals within at most 0.03 eV. Thus, this justifies the use of the PBE functional for estimating the effect of the band-gap corrections due to ZPR and thermal vibrations in the case of perovskite materials.
The band gap renormalizations due to zero-point motion and thermal vibrations are given in Table 6 and are illustrated in Fig. 2(b). In particular, we also give the corrections at T = 0 K corresponding to the ZPR. The compounds containing the light F atoms show the most prominent ZPR effects, with the largest one found for M-RbSnF 3 and amounting to −0.20 eV. The next group of compounds showing sizable ZPRs are those containing Cl atoms, among which R-CsGeCl 3 shows the most significant value of −0.14 eV. Increasing the temperature from 0 to 300 K produces further band-gap renormalization. Together with the ZPR, this leads to ∆E ZPR+T corrections ranging between −0.40 and 0.23 eV (Table 6), which cannot be neglected when aiming at high accuracies. We account for the combined effect of SOC and phonons by determining the two band-gap renormalizations separately, as these effects have been found to give a negligible cross coupling (see Supplementary Table 1). Adding ∆E ZPR+T to E SOC , we then obtain the band gap estimates E theory , which further improve the comparison with experiment leading to a MAE of 0.17 eV and a MARE of 6% (cf. Table 4).
The band-gap renormalization due to thermal vibrations can be compared with previous work in the literature. The band-gap variation of γ-CsPbBr 3 between 0 and 300 K has been investigated through constant-volume molecular dynamics simulations [69] yielding a shift of 0.01 eV in good agreement with our estimate of ∆E T = 0.02 eV (cf. Table 6). Moreover, for large-size nanocrystals of γ-CsPbI 3 , temperature-dependent luminescence spectra show that the observed band gap remains invariant with temperature. [70] This behavior stems from the combined effect of lattice expansion and phonon renormalization. To establish a connection with our work, we estimate the correction due to lattice expansion using the experimental thermal lattice-expansion coefficient of 0.39 × 10 −4 K −1 . [71] We find a band-gap increase of 0.07 eV over a thermal range of 300 K, which almost perfectly compensates the corresponding band-gap reduction of 0.06 eV due to thermal vibrations (cf. ∆E T in Table  6). Hence, our correction is fully consistent with the null result observed in the experiment. [70] Our method relies on the harmonic approximation for describing the vibrational properties of the investigated materials. Patrick et al. performed a careful study including anharmonic effects in both cubic and noncubic CsSnI 3 . [44] To evaluate possible contributions coming from anharmonic effects, we applied our harmonic scheme to the noncubic β-CsSnI 3 and γ-CsSnI 3 at 380 and 300 K, respectively, finding band-gap renormalizations of 0.23 and 0.11 eV due to ZPR + T . These values are very close to the respective results of 0.16 and 0.11 eV found by Patrick et al. including anharmonic effects. This suggests that anharmonic effects only lead to small corrections for noncubic perovskites. However, we emphasize that anharmonic effects can be larger in the cubic phases because the structure then undergoes noticeable distortion moving away from the highly symmetric structure, which only holds on average. [49] The correction ∆E ZPR+T can either increase or reduce the band gap. Most materials in our selection exhibit the conventional Varshni effect, i.e. their band gap decreases with increasing temperature (red shift). [72] The ZPR and the band gap reduction due to thermal vibrations then contribute in the same direction with corrections being as significant as −0.40 eV in the case of R-CsGeCl 3 . At variance, some materials, most particularly the Sn-based γ and β phases, show a band gap opening with increasing temperature (blue shift). This unusual effect is referred to as anomalous band gap shift and has experimentally been observed for a variety of hybrid organic-inorganic lead halide perovskites. [73] In the case of octahedral structures, the origin of the band-gap renormalizations can be understood in terms of B-X bond-length and B-X-B bond-angle variations (see Supplementary Table 7). Indeed, it is well-known that reducing the B-X bond lengths or the B-X-B bond angles generally leads to lower band gaps due to enhanced orbital coupling. [74][75][76] Practically all octahedral structures in our selection show reduced bond angles and increased bond lengths when going from 0 to 300 K (cf. Supplementary Table 7). These trends oppositely impact the band gap with a combined effect that can either lead to band-gap opening or to band-gap closure. More specifically, for the rhombohedral structures, the effect due to the large bond-angle reductions (5 • -10 • ) apparently dominates, leading to overall band-gap closing. The larger the bond-angle reduction, the lower the band gap. For γ phases, the bondangle reductions are generally smaller (1 • -4 • ) and their effect is similar to the effect due to the bond-length increases, resulting together in small bandgap variations, either red-shifted or blue-shifted. The structural variations of the δ phases are generally less significant, with bond-angle reductions and bond-length increases amounting to at most 1 • and 0.02Å, respectively. Nevertheless, their combined effect still yields quite sizable band-gap reductions, with the largest being 0.33 eV in the case of δ-CsPbBr 3 .

Discussion
In Fig. 4, we show how the calculated corrections progressively improve the theoretical band gaps. We here consider all the materials for which experimental band gap values are available. The corrections due to SOC and ZPR+T are added in sequence to the bare band gaps calculated with the DSH hybrid functional. Moreover, our best estimate for the fundamental band gaps E theory are compared with optical band gaps from experiment. We hereby neglect excitonic effects, which are estimated to be on the order of 10 to 100 meV for this class of materials. [77][78][79][80] For instance, photoluminescence experiments on γ-CsSnI 3 suggest a binding energy of 18 meV. [78] Calculations based on the Bethe-Salpeter equation also give very small exciton binding energies for CsGeX 3 (X = Cl, Br, I), on the order of at most ∼1 meV. [80] The residual deviations in Fig. 4 between our theoretical estimates and the experimental values should be assigned to the intrinsic limitations of the adopted theory. The fact that the errors are not systematic reinforces this point of view. In fact, neglected effects such as exciton binding energies, polar corrections due to LO phonon coupling, anharmonic vibrational effects, errors resulting from the use of PBE for SOC and phonon renormalizations, and cross coupling between SOC and temperature effects, have all been demonstrated to affect the results to a lower extent than the residual error. Residual deviations, such as encountered for R-CsGeBr 3 , γ-CsPbCI 3 , and M-CsSnCl 3 should therefore not be attributed to one particular effect. It is likely that the largest residual errors ranging up to 0.49 eV for an individual perovskite results from limitations of the adopted electronic-structure method. It should be remarked that even state-of-the-art GW calculations show non-systematic behavior with individual errors ranging up to 0.44 eV. [15] As seen in Fig. 4, the agreement with experiment systematically improves when applying corrections due to ∆E SOC and ∆E ZPR+T in sequence. On average, ∆E SOC and ∆E ZPR+T amount to −0.36 and −0.15 eV, respectively, corresponding to an average global correction E full of 0.51 eV (cf. Fig. 2). The MAEs of E bare , E SOC , and E theory are 0.57, 0.26, and 0.17 eV with respect to experiment, with respective MAREs of 27%, 11%, and 6% (see also Table 4). The finally achieved accuracy of 0.17 eV (6%) informs us that the applied corrections are essential to achieve this level of agreement with experiment, since they are of the same size or larger on average. Furthermore, the achieved accuracy is comparable to that of state-of-the-art GW calculations for extended sets of materials. [12,13,15,81] To support the use of the material-specific functional DSH, we carry out a comparison starting from standard functionals, such as PBE, HSE06, and PBE0(0.25). As discussed, our band-gap corrections can all accurately be calculated using the PBE functional. The use of a different functional then only affects the calculation of the bare band gap. The detailed results of these calculations are given in the Supplementary Tables 8, 9, and 10 for PBE, HSE06, and PBE0(0.25), respectively, and can be directly compared with the DSH results given in Table 4. To assess the performance of all the investigated functionals, we have summarized in Table 7 their respective MAEs and MAREs for the calculated band gaps with respect to the experimental values. We notice that for the PBE and HSE06 functionals the agreement with experimental band gaps steadily deteriorates upon including the corrections. In particular, the E bare from HSE06 already shows relatively low MAE and MARE of 0.36 eV and 14%, respectively. However, this result should be considered accidental and due to error cancellation. The results in Table 7 also indicate that the performance of PBE0(0.25) for the materials investigated is rather good (MAE of 0.28 eV; MARE of 11%) showing systematic improvement upon the inclusion of the corrections. As seen in Table 7, DSH outperforms all the other functionals with a MAE of 0.17 eV and a MARE of 6%. We stress that the best agreement between theory and experiment is obtained when all the various contributing physical effects are properly estimated.
In hindsight, it is of interest to question whether the approach developed is indeed superior to simple band-gap predictions made on the basis of a trivial PBE calculations. To address this issue, we develop a simple model based on the PBE band gap and test it on the nine materials for which experimental band gaps are available. From Supplementary Fig. 9, one infers that the relation between PBE and experimental band gaps is approximately linear. To obtain a model prediction, we then fit the data with a linear function: E model gap = 1.39 · E PBE gap + 0.47 eV. Such a model prediction yields a MAE of 0.30 eV with a maximum error as high as 0.92 eV. These errors correspond to the data fitted and are susceptible to increase when a larger set of materials is taken under consideration. Nevertheless, the MAE value and the maximum error of this simple model are already noticeably worse than the MAE of 0.17 eV and maximal error of 0.49 eV found with the scheme introduced in this work. This provides further support to the validity of our scheme.
We remark that the present methodology has here successfully been applied to the large class of noncubic metal halide perovskite materials. To validate this methodology to this class of materials, we verified that the consideration of anharmonicity, cross-coupling between SOC and phonon effects, and nonadiabatic effects in phonon renormalization due to long-range polar interactions, all lead to negligible contributions to the band gap. Thus, it should be understood that the present methodology can reach the high level of accuracy found in our work only to extent that these conditions hold.
To summarize, we have investigated the fundamental band gaps for a set of inorganic halide perovskites through the use of various functionals. To ensure a meaningful comparison with experiment, we have included band-gap corrections due to spin-orbit coupling and to thermal vibrations including zero-point motions. In particular, we use the dielectric-dependent hybrid functional DSH with parameters fixed through the dielectric function. Among the functionals considered, the functional DSH stands out providing a high level of accuracy compared to experimental band gaps (MAE = 0.17 eV; MARE = 6%). The achieved accuracy relies on the consideration of both corrections, which lead to a step by step improvement of the agreement with experiment. Our final accuracy is comparable to the state of the art as far as band-gap estimates are concerned.
Additionally, our scheme is designed to be computationally efficient. The band gap calculations with the hybrid functional are performed only once for a single crystalline structure for each material. We demonstrate that both band-gap corrections due to spin-orbit coupling and to thermal vibrations including zero-point motions can accurately be determined with the computationally more convenient semilocal PBE functional. In particular, the effects due to zero-point motions and thermal vibrations are estimated through the special displacement method, which gives the band-gap correction through a one-shot calculation, thereby speeding up the effort with respect to the more lengthy molecular dynamics simulations. Consequently, our scheme is not only highly accurate but also provides band gaps in a computationally efficient way. These features make our approach thus convenient for material screening procedures [82][83][84] in which large databases of metal halide perovskites are taken under consideration.

Range-separated dielectric-dependent hybrid functional
We use a range-separated hybrid functional formulation in which the parameters are determined from the dielectric response of the system. In such a scheme, the Coulomb potential is partitioned into a short-range (SR) and a long-range (LR) component through an error function [18, 25-27, 53, 85]: where µ is the range-separation parameter. Next, the SR and LR components appearing in the nonlocal exchange potential are separately decomposed into nonlocal Fock and semilocal PBE [86] exchange terms through admixtures defined by α SR and α LR : Numerous widely used density functionals can be recovered from Eq. , also belong to this class of functionals and can be found by setting α SR = 1 and α LR = 1/ǫ ∞ . However, these two DDH schemes differ in the way the parameter µ is set. In DD-RSH-CAM, µ is derived from a fit to the dielectric function calculated through first-principles linear response [4,26], while in DSH this parameter is approximated through an empirical expression that can be evaluated analytically [27]. As we show below, the two proposed values of µ are nevertheless rather close (see Supplementary Table 4). Therefore, for speeding up the computational effort, we adopt in this work the dielectricdependent functional denoted DSH in Ref. 27 with the empirical expression for µ proposed therein. More specifically, the empirical expression for µ proposed in ref. 27 stems from the nonlocal screened Coulomb potential derived by Shimazaki and Asai [87]: Here,k 2 TF = k 2 TF (ǫ ∞ − 1) −1 + 1 /α and k TF = 2 6 (3n/π) is the Thomas-Fermi screening length, which depends on the valence charge density n. Following refs 25, 27, we include the d electrons in the valence density when the valence band is dominated by d orbitals, namely for Ge, Sn, Pb, Br, and I.
The coefficient α = 1.563 is considered to be independent of material. [88,89] Equation (4) can be well approximated by the error-function expression [87], provided one sets µ = 2k TF /3. In Eqs. (4) and (5), the first term considers only the SR contribution that corresponds to metallic screening, whereas the second term considers the full-range (FR) of the Coulomb interaction representing dielectric screening.
[27] Therefore, the nonlocal Fock exchange found by Shimazaki and Asai can be expressed as: To use this expression as nonlocal exchange in a hybrid functional formulation, one needs to add a compensating semilocal exchange in the long range.
[27] Thus, one obtains This expression is formally the same as the one obtained by Chen et al. following an alternative derivation path [26]. This expression is also consistent with Eq. (3), through which it can be obtained by setting α SR = 1 and α LR = 1/ǫ ∞ . Hence, we use in this work the exchange potential given in Eq. (7) with µ = 2k TF /3.

Special displacement method
Based on the theory introduced by Williams [90] and Lax [91] (WL) and the use of the harmonic approximation, the temperature-dependent band gap can be expressed as [35]: where the product runs over all modes ν, and Q is used to indicate collectively the configuration defined by the normal coordinates Q ν . This expression can be readily understood as the thermal average of E Q gap with weights ν exp −Q 2 ν / 2σ 2 ν,T / √ 2πσ ν,T , in which σ ν,T is a spatial Gaussian broadening given by with l ν = /(2M p ω ν ) and n ν,T = {exp [ ω ν / (k B T )]−1} −1 for the zero-point vibrational amplitude and the Bose-Einstein distribution, respectively. M p and ω ν denote the mass of the proton and the frequency of the νth normal mode. The broadening σ ν,T accounts for both the zero-point and thermal effects, and is denoted as σ ZPR+T . For T → 0 and n ν,T → 0, the broadening only results from the zero-point amplitude through The average in Eq. (8) for a temperature T can be obtained by considering a single distorted atomic configuration, which is created by displacing the atoms by an amount of ∆τ kα from the equilibrium structure [35]: where ν runs over all the non-translational modes, α indicates a Cartesian direction x, y or z, M k is the mass of the kth nucleus, and e kα,ν is the vibrational eigenmode of the νth normal mode. We note that depending on whether we use σ ν,T from Eq. (9) or from Eq. (10), we either describe the ZPR and thermal effects together or just the ZPR effects alone. From the study of Zacharias and Giustino [35], we infer that the band gap calculated for the distorted structure defined by the displacements in Eq. (11) reproduces the thermally averaged band gap defined in Eq. (8) within an accuracy of about 50 meV when the adopted supercell contains 150 atoms or more. [35] We consider this level of accuracy sufficient for the purpose of our work. The special displacement method allows one to obtain band-gap renormalizations due to phonon-coupling within the context of the adiabatic approximation. [35,36] The problematic limit at small momentum transfers in the summation over the Brillouin zone [37] is effectively handled in our supercell calculations by considering larger and larger supercells. We systematically considered supercells of varying size (containing a number atoms ranging between 160 to 625) and obtained constant results within 0.05 eV. In practice, the frequencies remain always finite and the ω 0 of the recipe of Zacharias and Giustino [36] could be taken to correspond to the lowest frequency found for the series of considered supercell sizes. For any lower ω 0 , the results would remain unchanged and the obtained values correspond to a finite result. However, the adiabatic approximation does not account for the long-range polar coupling to the longitudinal optical phonons, which results from a nonadiabatic treatment. [43] In this work, this effect has been estimated through the one-band Fröhlich model (see Sec. 2.3 and Supplementary Note 1).

Calculation details
The DSH band gaps of the ideal crystalline systems are determined through the implementation of DD-RSH-CAM functionals [26] available in the Quantumespresso software suites. [92] We use fully-relativistic pseudopotentials generated by the optimized norm-conserving Vanderbilt pseudopotential scheme [93] to account for spin-orbit coupling. All the calculations are carried out with the stringent set of the Pseudo Dojo (vailable at http://www.pseudodojo.org) to ensure band gaps converged within 0.1 eV. A kinetic plane-wave cutoff of 100 Ry is used for all materials together with sufficiently dense k-point grids.
In the band-gap calculations, we use unit cells with atomic structures and lattice parameters taken from experiment (cf. Supplementary Table 2). The high-frequency dielectric constant ǫ ∞ is calculated at the PBE level through the application of finite electric fields [50,94] and is used to fix the parameter α LR of the DSH functional.
For the application of the special displacement method, [35,36] we use the cp2k suite of codes, [95] which comprises an efficient tool for determining the vibrational frequencies and modes through finite differences. The core-valence interactions are described through Goedecker-Teter-Hutter (GTH) pseudopotentials. [96,97] We use double-zeta basis sets of MOLOPT quality. [98] The plane-wave energy cutoff for the electron density is set to 600 Ry to ensure the convergence of the total energy. The vibrational frequencies ω ν and eigenmodes e kα,ν required for the distorted structure specified in Eq. (11) are obtained through finite displacements of 0.01Å from fully relaxed atomic positions. We use Γ-point sampling with the supercells defined in the Supplementary Table 5. As an alternative scheme for obtaining vibrational properties, the non-diagonal approach proposed by Lloyd-Williams et al. [99] could be used to further reduce the involved computational cost.
The band gaps are determined from linear extrapolations of the EDOS obtained for the supercells specified in the Supplementary Table 5. The cp2k results are systematically obtained at the PBE level of theory, with the exception of the test results in Table 5, which required the use of hybrid functionals. In this case, the auxiliary density matrix method (ADMM) is employed to speed up the calculations. [100] Our efficient band-gap calculations contain several parts that contribute to the overall computational cost. The first contribution results from the determination of the dielectric constant, but the cost of this calculation is not significant since it is carried out at the PBE level. The first important contribution comes from the band-gap calculation for the ideal crystalline system with the DSH functional. Since the cost of this calculation is easy to evaluate, we use it here to set the unit. The calculation of the SOC correction is carried out at the PBE level and is negligible on the scale of our cost unit. The second important contribution comes from the determination of the vibrational properties at the PBE level. The corresponding cost of this part depends on the number atoms considered, but it is comparable to the cost unit for systems containing 270 atoms. For some materials, the EDOS needs to be determined with rather large supercells, but still at the PBE level. For the largest case considered (2880 atoms), this part gives a cost that amounts to half a cost unit and is thus quite negligible in most cases. All together, the dominant costs for the band-gap determination result from the vibrational properties and from the DSH hybrid-functional calculation in comparable amounts.
Data availability statement. The structure used in the paper can be found in the materials cloud: https://archive.materialscloud.org/record/2022.35.
Code availability statement. The relevant codes in this study are available from the corresponding authors upon the reasonable request.           (Table 2), ∆E SOC is the correction due to SOC (Table 3), and ∆E ZPR+T the correction accounting for ZPR and thermal vibrations (    [59,60]. The percentages correspond to the mean absolute relative errors (MAREs). γ-C sS nB r3 γ-C sP bI3 R -C sG eC l3 R -C sG eI3 γ-C sS nI3 γ-C sP bB r3 M -C sS nC l3