Theory of the special displacement method for electronic structure calculations at finite temperature

Calculations of electronic and optical properties of solids at finite temperature including electronphonon interactions and quantum zero-point renormalization have enjoyed considerable progress during the past few years. Among the emerging methodologies in this area, we recently proposed a new approach to compute optical spectra at finite temperature including phonon-assisted quantum processes via a single supercell calculation [Zacharias and Giustino, Phys. Rev. B 94, 075125 (2016)]. In the present work we considerably expand the scope of our previous theory starting from a compact reciprocal space formulation, and we demonstrate that this improved approach provides accurate temperature-dependent band structures in three-dimensional and two-dimensional materials, using a special set of atomic displacements in a single supercell calculation. We also demonstrate that our special displacement reproduces the thermal ellipsoids obtained from X-ray crystallography, and yields accurate thermal averages of the mean-square atomic displacements. At a more fundamental level, we show that the special displacement represents an exact single-point approximant of an imaginary-time Feynman’s path integral for the lattice dynamics. This enhanced version of the special displacement method enables non-perturbative, robust, and straightforward ab initio calculations of the electronic and optical properties of solids at finite temperature, and can easily be used as a post-processing step to any electronic structure code. To illustrate the capabilities of this method, we investigate the temperature-dependent band structures and atomic displacement parameters of prototypical nonpolar and polar semiconductors and of a prototypical two-dimensional semiconductor, namely Si, GaAs, and monolayer MoS2, and we obtain excellent agreement with previous calculations and experiments. Given its simplicity and numerical stability, the present development is suited for high-throughput calculations of band structures, quasiparticle corrections, optical spectra, and transport coefficients at finite temperature.


I. INTRODUCTION
The calculation of the electronic and optical properties of materials at finite temperature is a long-standing challenge for ab initio electronic structure methods. In typical semiconductors, insulators, metals, and semiconductors, the key mechanism leading to temperature dependent properties is the thermal motion of the atoms in the crystal lattice, and the effect of this motion on the electronic structure of the system. Recent advances have made it possible to study these effects from first principles with predictive accuracy [1].
In the case of semiconductors and insulators, one problem that attracted considerable attention during the past decade is the electron-phonon renormalization of band structures, including both quantum zero-point effects and temperature dependence . This problem is becoming increasingly important as we strive to achieve close quantitative agreement between ab initio calculations and experimental data. Furthermore this problem underpins calculations of many important properties, from temperature-dependent optical absorption [24][25][26] and emission spectra [27] to temperature-dependent * fgiustino@oden.utexas.edu transport coefficients [28,29].
The calculation of temperature-dependent electronic and optical properties has been demonstrated using both perturbative approaches based on density functional perturbation theory [30] (DFPT) calculations in the crystal unit cell, and non-perturbative approaches based on density functional theory (DFT) calculations in large supercells. In perturbative methods the key ingredient of the calculations is the electron-phonon matrix element, that is the scattering matrix element between two Kohn-Sham states which results from the variation of the self-consistent potential associated with a vibrational eigenmode. The matrix elements are employed to obtain temperature-dependent band structures in the Allen-Heine (AH) method [31,32], and to compute indirect optical absorption in the Hall, Bardeen, and Blatt theory [33,34]. These approaches have enjoyed considerable success during the past decade across a broad range of materials [3-5, 7-14, 24, 35-38].
Non-perturbative supercell-based methods offer an alternative approach to calculations of electronic and optical properties at finite temperature. The central idea of these methods is that the interactions between electrons and phonons can be described using large supercells where the atoms are displaced to capture thermal disorder [39]. Supercell methods derive from earlier frozenphonon approaches [2,40], and do not require the ex-plicit evaluation of electron-phonon matrix elements. In order to sample thermal disorder, the first attempts in this area relied on the Monte Carlo sampling of quantum nuclear wavefunctions, either in the atomic configuration space [25,41,42] or in the space of normal vibrational modes [19,20,43]. The formal basis for this approach is provided by a theory developed in the 1950s by Williams [44] and Lax [45] to study defects in solids, and recently extended by Zacharias et al. to finitetemperature optical spectra, phonon-assisted optical processes, and band gap renormalization [25].
The central idea behind Monte Carlo supercell methods is that the temperature dependence of the dielectric function (ω; T ) results from a multivariate Gaussian integral given by [26]: where the product runs over the normal coordinates x ν , and each width σ ν,T of the independent Gaussians represents the root mean square displacement of the atoms at temperature T along a phonon mode ν. The interpretation of Eq. (1) is that the temperature dependence of the spectrum can directly be obtained as an ensemble average over the spectra calculated at fixed nuclear coordinates x, whose probability density function is a multivariate Gaussian.
The disadvantages of stochastic supercell approaches are that (i) it is difficult to control the rate of convergence with the number of random samples, and (ii) the final outcome of the calculations is a distribution of values at each temperature, e.g. band gaps, rather than a unique value as in perturbative methods. Furthermore, performing stochastic sampling on large supercells can be prohibitively costly from the computational standpoint. In order to overcome these limitations, in Ref. [26] we demonstrated that it is possible to replace the stochastic sampling of the nuclear wavefunctions by a single supercell calculation with a precise choice of atomic displacements. The possibility of studying electronic and optical spectra at finite temperature using a single supercell calculation enabled several interesting applications. For example, in Ref. [26] we reported temperature-dependent phonon-assisted optical spectra of Si, diamond, and GaAs based on this method; Refs. [46], [47] and [48] applied this method to the dielectric function of Zn 2 Mo 3 O 8 , metallic hydrogen, and BaSnO 3 , respectively; Ref. [23] used this technique to compute GW quasiparticle band gaps at finite temperature for diamond, BN, SiC, Si, AlP, ZnO, GaN, ZnS, GaP, AlAs, ZnSe, CdS, GaAs, Ge, AlSb, CdSe, ZnTe, and CdTe; Ref. [49] calculated exciton-phonon couplings in hexagonal boron nitride; and Refs. [28,29] demonstrated calculations of finite-temperature carrier mobilities in silicon n-i-n and p-n junctions using this approach. In retrospect, these successes across a broad range of applications are not too surprising, since the method of Ref. [26] is designed to provide, in the limit of large supercell, the exact thermodynamic average of any property that can be expressed in the form of Fermi's golden rule: this includes in principle photoelectron spectra (hence band gaps and band structures), optical spectra, tunneling spectra, and transport coefficients.
One potential limitation of Ref. [26] is that the theory was developed using a Γ-point formalism, therefore phonon calculations to determine the vibrational eigenmodes in the supercell are demanding. Furthermore, by addressing only Γ-point properties, the calculations are limited to angle-integrated spectra, such as density-ofstates (DOS) and optical absorption spectra. Finally, it has been pointed out that the choice of the special displacement employed in Ref. [26] might be improved to reduce the size of the supercell required to achieve convergence [23].
In this manuscript, we significantly expand the scope of the methodology introduced in Ref. [26] by formulating the theory within a compact reciprocal space formulation, and exploiting translational invariance and timereversal symmetry. This upgrade allows us to determine the special supercell displacement using quantities that are computed in a crystal unit cell via DFPT. In a nutshell, we demonstrate that a supercell calculation where the atoms are displaced according to: (2n qν,T + 1) 1 2 × 2 Re e iq·Rp e κ,ν (q) .
yields the exact thermodynamic average of electronic and optical properties at the temperature T in the adiabatic and harmonic approximations. In the above expressions ∆τ pκ indicates the displacement of the atom κ with mass M κ in the unit cell with lattice vector R p , and N p is the number of unit cells in the supercell. e κ,ν (q) is the phonon polarization vector of the normal mode (normalized within the unit cell) with wavevector q, branch index ν, frequency ω qν , and Bose-Einstein occupation n qν,T . The quantities S qν are signs, + or −, which depend on the normal mode, as specified in Sec. V. The summation is restricted to phonon wavevectors that are not timereversal partners. In particular, we define this group of phonons as set B. The real part in the equation arises from grouping together a phonon with wavevector q ∈ B and its partner −q. Phonon wavevectors that coincide with their time-reversal partners are grouped in a finite set A, and their contribution to the atomic displacements vanishes in the limit of dense Brillouin-zone sampling. This partitioning is discussed in greater detail in Sec. II. Equation (2) reduces to our previous prescription provided in Ref. [26] if we perform Γ-point sampling, but is much more powerful since the construction of the special displacement relies on DFPT calculations in the unit cell. Therefore the only expensive step is one calculation of the desired property in a large supercell. To demonstrate the power of this approach, we report for the first time complete temperature-dependent band structures based on Eq. (2) for both three-dimensional semiconductors, Si and GaAs, and two-dimensional semiconductors (MoS 2 ), and we show that this method delivers an accuracy comparable to perturbative techniques, without requiring the so-called rigid-ion approximation of the Debye-Waller self-energy [12,14]. We also demonstrate that Eq. (2) is an exact single-point approximant to an imaginary-time path integral in the space of nuclear configurations, and that it reproduces very accurately the thermal displacement ellipsoids measured by X-ray diffraction (XRD).
The method of Ref. [26] was originally named 'oneshot configuration', but subsequently it has variably been referred to as the 'WL-HA' method [47] or the 'STD' method [28]. In order to avoid a proliferation of acronyms, in this manuscript we will refer to the general theory as the special displacement method (SDM), and to Eq. (2) as the ZG displacement.
The manuscript is organized as follows. In Sec. II we describe the Williams-Lax theory that underpins our method, and we write the key equations using a reciprocal-space formulation. In Sec. III we establish the connection between the ZG displacement and thermal ellipsoids. In Sec. IV we show that the Williams-Lax theory can be recast in the form of an imaginary-time Feynman's path integral, and that the ZG displacement enables an efficient single-point evaluation of this path integral. Section V is devoted to the construction of the ZG displacement and the determination of the mode signs S qν in Eq. (2). In Sec. VI we provide all computational details of the present calculations, including the unfolding of bands and spectral functions. In Sec. VII we briefly outline the procedure that we follow to determine the ZG displacement, and we demonstrate our main computational results. We show the temperature-dependent thermal ellipsoids, band structures and band gaps of Si, GaAs, and monolayer MoS 2 using the SDM. Here we also provide detailed convergence tests and show that accurate results can be obtained even with relatively small supercells. In Sec. VIII we summarize our key findings and discuss future directions. More technical aspects are left to the Appendices.

A. General remarks
The theoretical framework underlying our present methodology was laid out in two works by Williams [44] and Lax [45] in the 1950s, which we refer to collectively as the Williams-Lax theory. This theory starts from the Herzberg-Teller rate [42,50] that describes transitions between coupled electron-phonon states driven by an external field, and replaces the final quantum nuclear states by a semiclassical continuum. Formally this step corresponds to neglecting commutators involving the nuclear kinetic energy operator, and is related (albeit not identical) to the adiabatic Born-Oppenheimer approximation [42]. This approach is also closely related to Feynman path integrals [51], as we discuss in detail in Sec. IV.
If we denote a Born-Oppenheimer quantum state using the ket |αn , with the Greek letter referring to the electronic part and the integer to the nuclear part [42], then the Williams-Lax theory provides a semiclassical approximation for the transition rate from an initial state |αn to all final states |βm with the same β and every possible vibrational state m, evaluated using Fermi's golden rule: In this expression ω is the frequency of the driving external field, the χ αn ({τ }) are the quantum nuclear wavefunctions for the potential energy surface associated with the electronic state α, and {τ } denotes the set of all atomic coordinates. Γ {τ } α→β (ω) is the transition rate evaluated with the atoms clamped in the positions {τ }: α→β denote the energies of the initial and final electronic states and the associated transition matrix elements, respectively, all evaluated at clamped atoms. For example Γ {τ } α→β (ω) can represent the optical transition rates for light absorption, or the charge current in electron tunneling. Equation (3) formalizes the intuitive concept that the transition rate including quantum nuclear effects can be obtained by averaging 'static' clamped-ion rates over the nuclear probability distributions |χ αn ({τ })| 2 of the initial state.
The transition rate at a finite temperature T is obtained from Eq. (3) by carrying out a canonical average over the vibrational quantum numbers of the initial state [26,42,45]: where k B stands for the Boltzmann constant and Z = n exp(−E αn /k B T ) is the canonical partition function. Here E αn denotes the vibrational energy of the state |αn . This approach has been used in a number of investigations and is well established by now [23,25,26,42,[44][45][46][52][53][54].
Since Eqs. (3)-(5) involve the calculations of transition rates for displaced atomic configurations {τ }, the results automatically incorporate the effect of electronphonon couplings, because the method probes the change in the electronic energies and wavefunctions upon displacing the atoms. We emphasize that Eqs. (3)- (5) apply to any property that can be obtained from the Fermi golden rule, in particular optical spectra, photoemission spectra, transport coefficients, and tunneling spectra. Furthermore, these equations provide the formal basis for computing temperature-dependent electronic eigenvalues. Indeed, if we consider transitions from an electronic state |αn of a solid into free electron states, as in photoemission experiments, the final free electron states are independent of atomic positions, therefore the first frequency moment of the transition rates, dω ω Γ α→β (ω, T ), yields the temperaturedependent electronic energy of the initial state: (6) This relation states that, in the Williams-Lax theory, the temperature-dependent electronic energy E α (T ) can be obtained by averaging the position-dependent energy E {τ } α over thermal fluctuations. This observation forms the basis for all non-perturbative supercell calculations of temperature-dependent electronic properties, and is also at the basis of the adiabatic formulation of the Allen-Heine theory of temperature-dependent band structures [31].
An important implication of Eq. (6) is that, if we describe ions classically, then E α (T ) corresponds to averaging electronic energies over the snapshots of a molecular dynamics trajectory. This concept has been used extensively in the literature, but to the best of our knowledge the connection to the Williams-Lax theory is still not fully appreciated.

B. Reciprocal space formulation
In this section we recast Eqs. (3)-(5) in a language appropriate for ab initio calculations in periodic crystals. From now on we assume the harmonic approximation, so that we can define phonons in the usual way. We employ the same standard notation as in Ref. [1], and we choose the convention for the phase of vibrational eigenmodes as in Ref. [55].
We consider a Born-von Kármán supercell containing N p unit cells with M atoms each. The position vector of the atom κ in the unit cell p at zero temperature is τ κp = R p + τ κ , where R p indicates a direct lattice vector, and τ κ is the position vector in the primitive unit cell, with Cartesian components τ pκα . The atomic displacements from equilibrium are written as a linear combination of normal vibrational modes [1]: where the z qν are the complex-valued normal coordi-nates [56]. The inverse relation is: In both equations M 0 is an arbitrary reference mass, usually chosen to be the proton mass, and the summations run over all the 3M N p atomic degrees of freedom. For completeness, in Appendix A we give the standard textbook relations between normal modes. If we write z qν = x qν +iy qν , with x qν and y qν being the real normal coordinates, Eq. (8) and time-reversal symmetry imply x −qν = x qν and y −qν = −y qν . Therefore only half of the real normal coordinates are independent. The wavevectors q of these independent coordinates occupy half of the first Brillouin zone. This is illustrated in Fig. 1, where we consider a regular q-grids of size 8×8. This grid can be partitioned in three sets, A, B, and C, following Appendix B of Ref. [1]. Set A includes the q-points that are invariant under inversion, modulo a reciprocal lattice vector. Therefore the q-points in this set correspond to the center of the Brillouin zone, the centers of its faces, and the corners, as shown by the filled disks in Fig. 1. Set B includes all the q-points that are not inversion partners (modulo a reciprocal lattice vector), as shown by the empty circles in Fig. 1. Set C is obtained by changing the sign of all the points in B, and is denoted by crossed filled circles in Fig. 1. Using this partitioning we can write the atomic displacements in terms of independent real normal coordinates as: q∈A,ν e κα,ν (q)x qν cos(q · R p ) + 2Re q∈B,ν e iq·Rp e κα,ν (q)(x qν + iy qν ) .
Equation (9) allows us to write the total quantum nuclear wavefunction of the state |αn in the harmonic approximation [57]: Here we omitted the subscript α for notational simplicity. χ nqν (x) represents the wavefunction of a quantum harmonic oscillator with frequency ω qν and quantum number n qν ; in particular, for q ∈ A we have: while for q ∈ B we have: qν /2l 2 qν H nqν (x/l qν ). (12) Here l qν = ( /2M 0 ω qν ) 1/2 is the zero-point vibrational amplitude, and H m (x) denotes the Hermite polynomial of order m. The total energy of the state χ n is given by the standard expression: By substituting Eqs. (10)-(13) inside Eqs. (3) and (5), and using Mehler's sum rule for Hermite functions [58], we obtain the compact expression [42]: where σ qν is the mean square displacement of the normal mode at the temperature T and is given by: Here n qν,T = [exp( ω qν /k B T )−1] −1 is the Bose-Einstein occupation of the mode with frequency ω qν . Equation (14) states that the temperature-dependent transition rate Γ α→β (ω, T ) is obtained by averaging the transition rates calculated at clamped atoms for a variety of atomic configurations specified by the normal coordinates {x qν , y qν }, and that the average is to be taken using a multidimensional Gaussian importance function. In this context the temperature T sets the width of each Gaussian via Eq. (15). As a sanity check, we note that if Γ {τ } α→β (ω) does not depend on the atomic coordinates, Eq. (14) yields the correct temperature-independent rate.
The core of the special displacement method described in this manuscript is to identify one set of atomic displacements so that a single evaluation of Γ {xqν ,yqν } α→β (ω) yields the same result as the multidimensional integral in Eq. (14). From this perspective, the task of finding the ZG displacement is similar to finding the mean value point of a definite integral; the only difference is that we are dealing with a multi-dimensional integral with hundreds to thousands of variables.

III. ZG DISPLACEMENT AND THERMAL ELLIPSOIDS
In this section we discuss the physical meaning of the ZG displacement in Eq. (2). In particular we prove rigorously that the ZG displacement reproduces the displacement autocorrelation function obtained from the quantum canonical average [59], and therefore encodes information about the thermal ellipsoids measured via XRD. Furthermore, using explicit calculations, we demonstrate that the ZG displacement yields the correct thermal distribution of the atomic coordinates to all orders.
When a harmonic crystal is in thermodynamic equilibrium at the temperature T , the atomic displacements ∆τ pκα follow a normal distribution. This a direct consequence of the fact that the marginal distribution of a multivariate normal is also normal. In particular we have (cf. Eq. 7.2.21 of Ref. [60]): where the width σ κα (T ) of the Gaussian is given by (cf. Eq. 7.2.5 of Ref. [60]): Here the limit of dense Brillouin-zone sampling is implied; in this limit the contribution to the sum of phonons with q ∈ A vanishes. The width σ κα (T ) is a particular case of the tensor of anisotropic displacement parameters (ADPs), defined as [61]: where · T denotes the canonical average over vibrational quantum states. In fact, by combining Eqs. (18) and (16) it follows directly that σ 2 κα (T ) = U κ,αα (T ). The ADPs of Eq. (18) are the values employed to generate thermal ellipsoids when visualizing crystal structures at finite temperature.
If we start from the ZG displacement, and we take the mean square displacement of the atom κ over all the unit cells of the supercell, then we obtain precisely the thermodynamic average given by Eq. (17). In fact, by replacing the ZG displacement of Eq. (2) inside the sum p ∆τ 2 pκα /N p , and using the sum rule in Eq. (A5), we find immediately: The third line contains terms associated with q ∈ A phonons. As we prove below in Sec. V D, all contributions arising from q ∈ A phonons vanish in the limit of large supercell. Furthermore, the choice of signs S qν in the ZG displacement guarantees that also the second line of Eq. (19) vanishes in the thermodynamic limit. Therefore in this limit we recover precisely the mean-square displacements σ 2 κα (T ) of Eq. (17). In Sec. VII A we show that these results are confirmed by explicit calculations on Si, GaAs, and MoS 2 .

IV. THE ZG DISPLACEMENT AND FEYNMAN'S PATH INTEGRALS
An important property of the ZG displacement is that it can be understood as an exact single-point approximant to a Feynman's path integral [51].
In path integral approaches, such as path integral Monte Carlo simulations [52,62,63] and path integral molecular dynamics [47,64,65], the electronic and ionic degrees of freedom are decoupled using the adiabatic approximation, and the quantum nature of the atomic nuclei is taken into account by averaging the electronic properties over all possible trajectories of the atomic nuclei. In order to obtain the correct quantum statistics, each trajectory is weighted by exp(−S E / ), where S E is the Euclidean action associated with a path [51], and is defined in Eq. (22) below.
We now show that the Williams-Lax rate can be recast exactly in the language of path integrals. The equivalence between the two approaches was originally identified in Ref. [52]. By combining Eqs. (5) and (3) and using the relations χ αn ({τ }) = τ |αn , n |αn αn| = 1, and H n |αn = E αn |αn , we find: In these expressions |τ represents position eigenstates for all atomic coordinates, andĤ n is the nuclear Hamiltonian for the adiabatic potential energy surface E {τ } α . Now by combining together Eqs. (B1)-(B4), the Williams-Lax transition rate assumes the form of a thermodynamic Feynman path integral [47,51,52,63]: where the notation Dτ indicates the sum over all paths that begin from the set of atomic coordinates {τ } at the time t = 0 and go back to the same initial coordinates at the time t = /k B T . S E is the Euclidean action evaluated along each one of these paths: whereτ pκα represents the classical velocities of the nuclei and U α is the potential energy of the nuclear Hamiltonian. The partition function Z F is equal to the second line of Eq. (21). Equation (21) states that the Williams-Lax transition rates can be obtained by averaging the electronic transitions at clamped ions using thermodynamic path integrals as weighting coefficients. For each atomic configuration 'snapshot' {τ }, we must evaluate all closed-loop path integrals that begin and terminate at {τ }, see Fig. 2. At high temperature, k B T ω (with ω being the highest vibrational frequency), Eq. (21) assigns the largest weights to paths that minimize the classical action, therefore the William-Lax theory reduces to a standard configurational average using classical Boltzmann weights. Based on these observations, we can regard Eq. (21) as the link between the Williams-Lax theory, path-integral molecular dynamics calculations of electronic and optical properties at finite temperature, and finite-temperature calculations using classical molecular dynamics simulations [52]. We note that Eq. (21) does not require the harmonic approximation, and maintains its validity even for strongly anharmonic systems.
Since in the case of harmonic systems the special displacement method presented here yields the same result as the Williams-Lax theory (as we prove in Sec. V), and since we made no assumption about the form of the transition rates Γ {τ } α→β (ω) in Eq. (21), it follows that any calculation of finite-temperature properties using path integrals can be replaced altogether by a single evaluation with the ZG displacement. In symbols, for a generic observable O({τ }) which is a parametric function of the atomic coordinates {τ }, we have: (23) In other words, the ZG displacement provides an exact single point approximant to an thermodynamic Feynman path integral, as is shown schematically in Fig. 2. This result constitutes a significant computational simplification. For example, calculations of temperaturedependent band gaps using ring-polymer molecular dynamics require simulating hundreds of bids and averaging over thousands of snapshots [64], while the SDM can perform the same operation by performing a single supercell calculation.

V. THE ZG DISPLACEMENT
In this section, using a rigorous reciprocal space formulation, we prove that the ZG displacement given by Eq. (2) yields the correct Williams-Lax transition rates of Eq. (5) in the limit of large supercell. The key ingredients of our derivation are (i) the translational invariance of the lattice, (ii) time-reversal symmetry, and (iii) a smooth connection between vibrational eigenmodes at nearby wavevectors.
The strategy that we follow is similar to our previous work [26], that is we choose a displacement defined by normal mode coordinates of magnitude |x qν | = √ 2 σ qν for q ∈ A, and |x qν | = |y qν | = σ qν / √ 2 for q ∈ B. These choices leave us with the freedom to select appropriate signs that define in which direction the normal coordinates are being displaced (±). In Ref. [26] we considered a Γ-point formalism, the vibrational modes were ordered by increasing energy, and the choice of signs was simply the sequence +, −, +, · · · , see Eq. (5) of Ref. [26]. In the present work the situation is more complicated, because we now have to deal with phonon wavevectors q and phonon branch indices ν. As we discuss below, in this case a simple sequence of alternating signs is not sufficient to recover the correct Williams-Lax rate, and a more structured choice is necessary. The following sections describe how we make this choice and why.
A. The simplest case: one vibrational mode Before proceeding with the derivations, it is useful to examine the key idea of this method using the simplest example: a hypothetical system with only one vibrational mode. We rewrite Eq. (14) by replacing the transition rate with a generic 'observable' O {x} , and we call the WL thermal average of this quantity O(T ). For definiteness we consider the form of the integral corresponding to q ∈ A in Eq. (14): By expanding O {x} in powers of x near the equilibrium coordinates x = 0 and evaluating the integrals, we find: which is correct to fourth order in σ. Alternatively, we can evaluate O {x} for the normal coordinate x = √ 2 σ. This procedure yields: (26) If we now average the expansions for x = ± √ 2 σ we obtain: up to fourth order in σ. By comparing Eq. (25) and (27) we see that two evaluations of the observable O {x} at clamped nuclei are sufficient to reproduce the thermal average given by Eq. (24).
If we now consider the form of the integral corresponding to q ∈ B in Eq. (14), we arrive at a similar result. In this case the integral yielding the thermal average is to be replaced by two evaluations of the property O {x} at x = ±σ/ √ 2. In the case of a real system with many vibrational modes, our strategy is to exploit the same mechanism leading to Eq. (27), by leveraging the cancellation between 'similar' modes, i.e. modes of the same branch at adjacent q-points. In the case of many vibrational modes, the displacements can be chosen so as to reproduce the thermal average with only one atomic configuration.

B. Thermal average of an observable in the Williams-Lax formalism
In this section we derive the expression for the WL thermal average of an observable O {τ } which depends parametrically on the atomic positions {τ }, starting from Eq. (14) and following the same steps that lead from Eq. (24) to Eq. (25).
To second order in the displacements ∆τ pκα from equilibrium, the observable O {τ } reads: where O {0} indicates the observable evaluated at the equilibrium configuration. We can express the displacements in term of the real normal coordinates using Eq. (9) and the chain rule. The result for the linear variation is: while the result for the quadratic variation is: The derivation of Eqs. (29) and (30) (30), and we evaluate the integrals in the real normal coordinates x qν and y qν . There are only two types of integrals, those that are odd in x qν or y qν , which vanish identically, and those that involve the powers x 2 qν or y 2 qν . The latter are standard Gaussian integrals of the form dx x 2 e −x 2 = √ π/2. The resulting expression for the WL thermal average of O {τ } , correct to third order in σ qν , is: This result constitutes the generalization to many qpoints and many phonon branches of the result in Eq. (25) for the single-mode model.
Here we see that all odd powers in the normal coordinates and all second order cross-terms with q = q or ν = ν vanish upon integration.
C. Taylor expansion of an observable in normal coordinates

Linear variations
Now we perform a step similar to Eq. (26), but for the general case of many q-points and many phonon branches. Although we already derived the required Taylor expansion in Eqs. (28)- (30), it is convenient to find equivalent expressions which are easier to handle and where the translational invariance and time-reversal symmetry are built in from the outset.
Translational invariance implies that the derivative of O {τ } with respect to an atomic displacement be the same for every unit cell: This property can be used to prove that: To see this we write the derivatives using the chain rule and we employ Eq. (9). For q ∈ B we find: The sum rule in Eq. (A5) requires q = 0 for the sum over p to be nonzero, but this condition cannot be fulfilled when q ∈ B. The same argument applies to derivatives with respect to y qν . This proves the result in Eq. (33). The main consequence of this result is that the linear variation of the observable O {τ } in Eq. (29) takes the simple form: where the right-hand side contains only phonons with q ∈ A. This result indicates that, if we perform a supercell calculation with atoms displaced away from equilibrium, then the contribution to the observable O {τ } which is linear in the displacements comes entirely from phonons with wavevectors q ∈ A.

Quadratic variations
Here we employ again translational invariance to simplify the expression for the quadratic variations of the observable O {τ } appearing in Eq. (30). In this case translational invariance dictates: where the unit cell p − p corresponds to the lattice vector R p −R p . We now rewrite the second derivatives in Eq. (30) in terms of the real normal coordinates x qν and y qν using the chain rule and Eq. (9). After some lengthy but straightforward algebra we obtain: where δ q∈A = 1 for q ∈ A and 0 otherwise, and similarly for δ q∈B . From these expressions we see that, as a consequence of translational invariance, all derivatives with q = q vanish identically. By comparing Eqs. (38) and (40) we have: Furthermore, by using translational invariance in Eq. (39) it can be verified that: These properties allow us to simplify Eq. (30) as follows: As in the case of the linear variations in Sec. V C 1, also for the quadratic variations the contribution from phonons with q ∈ A vanishes in the limit of dense Brillouin zone sampling.

D. Choice of normal coordinates for the ZG displacement
The last step of the procedure outlined in Sec. V A is to choose the values of the normal coordinates x qν and y qν so that Eqs. (36) and (43) reproduce the WL average given by Eq. (31) in the limit of dense Brillouin-zone sampling.
By comparing Eq. (31) with Eq. (43) it is evident that the magnitude of the normal coordinates must be: |x qν | = √ 2 σ qν for q ∈ A, and |x qν | = |y qν | = σ qν / √ 2 for q ∈ B. The remaining degrees of freedom are the signs of the normal coordinates, say S qν,x = ±1 and S qν,y = ±1, which are to be determined. If we evaluate the observable O {τ } at the atomic positions specified by the normal coordinates x qν = √ 2 S qν,x σ qν for q ∈ A, and x qν = S qν,x σ qν / √ 2, y qν = S qν,y σ qν / √ 2 for q ∈ B, then the combination of Eq. (28) with Eqs. (36) and (43) yields, to second order in the displacements: Using the WL average in Eq. (31), the last expression takes the form: where ∆ A and ∆ B represent the deviation of O {τ } with respect to the WL average O(T ), and are given by: The term ∆ A only contains contributions from phonons with q ∈ A. In the limit of dense sampling of the Brillouin zone, the number of elements of A remains finite, while the number of elements in B goes to infinity. Therefore the Lesbegue measure of set A vanishes in this limit, and we have the result: Given that phonons with q ∈ A do not contribute for large supercells, to reproduce the WL average using a single atomic configuration, we only need to determine the signs S qν,x and S qν,y so that the term ∆ B in Eq. (47) vanishes. In order to simplify the task, we choose to assign the same sign to x qν and y qν : S qν,x = S qν,y = S qν . If we replace this choice of normal coordinates inside Eq. (7) and ignore the contributions from q ∈ A [which vanish in the limit of large supercell according to Eq. (48)], then we obtain precisely the ZG displacement in Eq. (2). Strictly speaking the result carries an additional phase e iπ/4 , but this can be absorbed in the eigenmodes (adding this extra phase does not pose any problem because the set q ∈ B does not contain time-reversal partners). Using the choice S qν,x = S qν,y = S qν inside Eqs. (38) and (39), we can rewrite ∆ B as follows: having defined: Clearly the quantity A νν (q) in the last expression is a relative of the lattice dynamical matrix, and the second derivatives ∂ 2 O {τ } /∂τ 0κα ∂τ p,κ α are relatives of the interatomic force constants, but for the observable O {τ } instead of the total energy. We note that the construction of the ZG displacement does not require the explicit evaluation of the second derivatives ∂ 2 O {τ } /∂τ 0κα ∂τ p,κ α to achieve the minimization of ∆ B . In order to determine the ZG displacement, we need to choose the signs S qν in such a way that ∆ B in Eq. (49) vanishes in the limit of large supercells. To this aim we note that in the limit of dense Brillouin-zone sampling, the vibrational eigenmodes e κα,ν (q) and eigenfrequencies ω qν can be chosen to be nearly the same between adjacent q-points: e κα,ν (q + ∆q) = e κα,ν (q).
(51) These conditions are always true for the frequencies, but are usually not fulfilled by the eigenmodes obtained from the diagonalization of the dynamical matrix. In fact nondegenerate eigenmodes may carry an arbitrary complex phase, and in case of degeneracy any rotation in the degenerate subspace is admissible. In order to enforce Eq. (51), we set up a smooth Bloch gauge in the Brillouin zone by performing a unitary rotation of all eigenmodes. This is described in detail in Sec. VI C. Then, by combining Eqs. (51) and Eq. (50) we have: Owing to this relation, if we select a set D of adjacent qpoints in the Brillouin zone, in the limit of dense sampling we can rewrite Eq. (49) as: whereq is the centroid of D. For the sum on the righthand side to vanish, we need to determine a set of signs S qν so that half of the products S qν S qν are positive and the other half is negative.
If we denote by n the number of phonon branches, we can assign 2 n unique combinations of signs to these modes. Let us call such combinations S The proof proceeds by induction: the equality is trivially verified for n = 2; for n > 2 assume that the result holds for a given n; when considering n+1 all the possible 2 n+1 combinations of signs are obtained by duplicating the previous 2 n combinations of n signs, and appending an extra sign at the end of each duplicate sequence. By construction, when ν, ν ≤ n we have i) qν = 0; when ν = n + 1 or ν = n + 1 every term of the sequence for n modes appears twice and with opposite signs, therefore also in this case the sum vanishes.
At this point we can perform a limiting procedure and consider a dense sampling of the Brillouin zone. If we partition the set B of q-points into disjoint subsets containing 2 n elements each, and we attach to each of these elements one of the 2 n unique sequences of signs S (i) qν , then the right-hand side of Eq. (53) vanishes identically, and so does the error ∆ B in Eq. (49): Taken together, Eqs. (48) and (54) constitute the formal proof that the ZG displacement provided by Eq. (2) yields exactly the Williams-Lax thermal average, to second order in the mean-square displacements. While we have not proven formally that the equivalence between the ZG displacement and the Williams-Lax theory extends beyond second order in the displacements, in Sec. VII A we show numerically that Eqs. (48) and (54) also hold for higher orders. A proof of equivalence of the two formalisms to all orders in the displacements was provided in Ref. [26] for the simpler case of Γ-point calculations. Therefore the present method does not rely only on the quadratic expansion of the relevant operator, and the equivalence between the ZG displacement and the multivariate Gaussian integral in Eq. (1) holds to all orders in the atomic displacements.
We also note that the same proof leading to Eq. (54) can be employed to demonstrate the equivalence between the ZG displacement and the mean-square displacements in Eq. (19). To this aim it is sufficient to replace the first line of Eq. (50) by a constant, absorb the imaginary factor e iπ/4 inside the eigenmodes, and observe that with these changes ∆ B of Eq. (49) reduces precisely to the second line of Eq. (19).

E. Additional considerations for calculations using small supercells
The proof outlined in the previous section considers the limit of very large supercells. For practical calculations it is important to ensure that the ZG displacement delivers good accuracy also in the case of computationally tractable, smaller supercells. This improvement can be achieved as follows.
The combinations of signs necessary to achieve the cancellation in Eq. (54) carries some redundancy: for every set S Obviously the latter combination yields the same products as the former, −S qν , therefore it is not useful to achieve the cancellation in Eq. (53). We can call the latter combination 'antithetic' to the former. In order to reduce the size of the supercell required to achieve convergence, we can eliminate the antithetic combinations.
This reasoning implies that the minimum number of q-points required to achieve exact cancellation of ∆ B [in the limit where Eq. (53) is satisfied] is precisely 2 n−1 . For example, in a tetrahedral semiconductor like Si or GaAs we would need at least 2 6−1 = 32 q-points and therefore supercells of size at least 4 × 4 × 4 are necessary to enable such cancellation.
To see how the choice of signs without antithetic sets works in a simple example, let us consider a system with n = 3 phonon branches. In this case we have 2 3 = 8 distinct combinations of signs as follows: Here the combinations i = 5, · · · , 8 are antithetic to i = 1, · · · , 4 and we can discard them. If we now consider a set D of 2 3−1 = 4 adjacent q-points, say q 1 , · · · , q 4 , we can choose the 'sign matrix' as follows: With this choice the summation in Eq. (53) becomes: By inspecting each column we see that terms with the same superscript ν, ν appear the same number of times with positive and negative signs, therefore the sum vanishes.
In the above example we selected 4 rows from the complete set of combinations in Eq. (55). How to perform this selection and in which order is not specified by the theory leading to Eq. (54). For example we could have selected rows 7, 8, 1, 2, and this choice would have led again to the same cancellation of Eq. (57). In practice, however, since we sample the Brillouin zone on a discrete grid of q-points, the property A νν (q) A νν (q) used in Eq. (53) does not hold for small supercells, and the cancellation is incomplete. This observation can be exploited to identify an optimal choice of signs from the complete set in Eq. (55). Indeed, we can select 2 n−1 inequivalent rows in such a way as to numerically minimize the error ∆ B .
A global optimization would require us to evaluate the second derivatives of the observable O {τ } for every atomic displacement ∆τ pκα , as it can be seen from Eq. (50). However, this step would be computationally costly and would make the special displacement method equivalent to standard frozen-phonon techniques [2,40].
Instead we make the observation that, as for the matrix of interatomic force constants, the second derivatives ∂ 2 O {τ } /∂τ 0κα ∂τ p,κ α must vanish for |R p | → ∞. Therefore the leading components of A νν (q) in Eq. (50) must be those for small |R p |. This suggests to minimize Eq. (49) by retaining only the second line of Eq. (50), multiplied by the respective signs, for every κα and κ α when R p = 0. In practice we chose to minimize the following normalized function of the signs S qν : The minimization of E({S qν }) with respect to {S qν } only involves the phonon eigenmodes and eigenfrequencies, and does not require any evaluation of the observable O {τ } . Hence the resulting ZG displacement is agnostic to the specific temperature-dependent property that is being computed, and can be constructed from quantities that are easily evaluated via standard DFPT in the crystal unit cell. When the signs are obtained by minimizing Eq. (58), the ZG displacements of a given crystal becomes uniquely defined for each temperature and each supercell.
The physical meaning of selecting the signs in such a way as to minimize Eq. (58) is that this choice leads precisely to the ZG displacements that best approximates the exact thermal mean-square displacements given by Eq. (17). Therefore this choice constitutes a way to construct displacements that reproduce XRD spectra as accurately as possible, as we show in Sec. VII A.

A. Computational setup
All calculations were performed using the local density approximation [66,67] to DFT for Si and GaAs, and the PBE generalized gradient approximation [68] for monolayer MoS 2 . We employed planewaves basis sets, norm-conserving pseudopotentials [69] for Si and GaAs, and ultrasoft pseudopotentials [70] for monolayer MoS 2 , as implemented in the Quantum ESPRESSO package [71]. The planewaves kinetic energy cutoff was set to 40 Ry for Si and GaAs, and 50 Ry for MoS 2 . To minimize interactions between periodic images of the monolayer, we used an interlayer separation of 15Å. The ZG displacement was constructed using vibrational eigenmodes and eigenfrequencies obtained from DFPT [30], starting from Brillouin zone grids of 8×8×8 points (Si), 8×8×8 points (GaAs), 8×8×1 points (monolayer MoS 2 ), and then using standard Fourier interpolation to generate dynamical matrices for coarser or denser grids.
In the case of polar materials (GaAs and monolayer Mo 2 ) our calculations correctly include the long-range component of the interatomic force constants via the nonanalytic correction to the dynamical matrix [72].
The algorithm used to construct the ZG displacement, the generation of a smooth connection between vibrational eigenmodes in the Brillouin zone, and the unfolding of spectral functions and band structures from the supercell to the unit cell are described in the following three sections, respectively.

B. Generation of the ZG displacement
To compute the ZG displacement in Eq. (2) we proceed as follows. First we perform a DFPT calculation of phonons in the crystal unit cell, using a coarse uniform grid of q-points in the Brillouin zone. Then we decide the size of the supercell for the ZG displacement, say N 1 × N 2 × N 3 unit cells, so that N p = N 1 N 2 N 3 . This choice sets the grid of q-points that we need to use in Eq. (2), namely q = From this grid we extract the sets A and B as illustrated in Fig. 1, and discard all remaining points. Using the DPFT results from the coarse Brillouin-zone grid, we perform standard Fourier interpolation to obtain the eigenmodes e κα,ν (q) and eigenfrequencies ω qν for q-points in sets A and B. This operation is identical to the procedure for calculating phonon dispersion relations using the matrix of interatomic force constants.
In order to enforce a smooth Berry connection between phonon eigenmodes as dictated by Eq. (51), we determine unitary rotations for adjacent q-points using a singularvalue decomposition (SVD), as described in Sec. VI C. This calculation only requires evaluating scalar products between eigenmodes at adjacent q-points.
At this point we are left with the determination of the signs S qν . Here we could proceed with a global optimization of all the signs, using Eq. (58). However, in order to keep the algorithm as simple as possible, we proceed with a partial optimization as follows. We order the qpoints along a path that is designed to span the entire set B in the Brillouin zone. A simple representative path in two dimensions is shown in Fig. 3. More complex choices such as the Peano-Hilbert space-filling curve are possible [73], but we have not explored these alternatives. The only requirement for the construction of this path is that it should exhibit 'locality', in the sense that pairs of q-points which are close in the Brillouin zone should also be close along the path, so that Eq. (53) is fulfilled. We then group the q-points along the path in sets of neighbors, with 2 n−1 q-points per set. The signs in each set are determined by extracting 2 n−1 sequences from the 2 n possible combinations, as explained in Sec. V D, excluding antithetic sets. Then we consider 2 n−1 consecutive sets like D, and we choose the signs as cyclic permutations of those from the first set. This procedure allows us to assign S qν for a total of 2 2(n−1) q-points. At this stage we evaluate the error E({S qν }) of Eq. (58) for these q-points. If the error is above a pre-defined threshold (defined as an external parameter, say δ = 5%) then we restart the procedure by performing a new selection of the 2 n−1 sign sequences and their order in the first D set.
We emphasize that the above optimization is unnecessary for large supercells, and many other choices of signs will lead to essentially the same results. This procedure is only advantageous when trying to obtain temperaturedependent properties using small supercells (e.g. 4×4×4 supercells). If the set B contains more than 2 2(n−1) q-points, then we continue the sequence by restarting from the beginning. We note that, while this procedure does not necessarily respect the periodic gauge across the Brillouin-zone boundaries, this does not constitute a limitation since we are only interested in minimizing E({S qν }).
Having established the signs S qν , we finally compute the ZG displacement using Eq. (2). In this expression the temperature T is an external parameter and enters via the Bose-Einstein factors n qν .

C. Smooth gauge of phonon eigenmodes along a path in reciprocal space
In order to satisfy Eq. (51), we perform unitary transformations of phonon eigenmodes at adjacent q-points on the space-filling curve described in Sec. VI B. The transformation is defined in such a way that similar eigenmodes at different q-points have a similar complex phase, and the ordering of eigenmodes is preserved in the case of branch crossing. This is equivalent to enforcing a smooth Berry connection across the Brillouin zone [74]. These ideas are related to standard concepts used in the theory of maximally-localized Wannier functions [75].
Given two adjacent reciprocal-space vectors q and q + ∆q, we seek for a transformation such that e κα,ν (q) and e κα,ν (q + ∆q) be as similar as possible. We can define similarity between eigenmodes using the overlap matrix: Using this definition and using the orthonormality relations in Eq. (A2) we have: e κα,ν (q + ∆q) = ν M νν e κα,ν (q).
For the modes e κα,ν (q) and e κα,ν (q + ∆q) to be as similar as possible, the overlap matrix M νν needs to be as close as possible to the identity matrix. Generally this is not the case since the diagonalization of the dynamical matrix at different q-points introduces an arbitrary gauge in the normal modes. To address this problem we perform a transformation of e κα,ν (q + ∆q) as follows: e κα,ν (q + ∆q) = ν U νν e κα,ν (q + ∆q), where U νν is a unitary rotation among the vibrational eigenmodes. After this transformation, the new overlap matrix reads: Minimizing the left-hand side with respect to U is equivalent to mazimizing Tr(M † U † + U M ). Using the properties of the matrix trace this is the same as maximizing Re Tr (U M ).
The matrix M is not Hermitian in general, but it can be decomposed via SVD as: where L and R are unitary matrices and S is a diagonal matrix with non-negative real values on the diagonal. With these definitions, it can be shown that the matrix U that maximizes Re Tr (U M ) is precisely U = RL † [76]. In order to use this strategy in practical calculations, we determine unitary transformations for each q-point along a space-filling curve, by evaluating overlap matrices M between each pair of successive q-points, say q 1 and q 2 . Then we apply the transformation U to the modes in q 2 , and repeat the procedure for q 2 and q 3 , and so on for all wavevectors. Figure 4 shows the eigenmodes e κα,ν (q) of silicon before and after the procedure just described. We can see that the resulting eigenmodes vary continuously between adjacent q-points, as desired.

D. Generation of temperature-dependent band structure by Brillouin-zone unfolding
In this section we present the recipe for calculating spectral functions that represent momentum-resolved density of states of the electrons at a given temperature. In order to obtain the spectral functions in the Brillouin zone of the primitive unit cell we employ the following unfolding procedure.
The spectral function is defined as [77,78]: where P mK,k (T ) are temperature-dependent spectral weights given by: In these expressions ε mK (T ) and ψ SC mK (T ) represent eigenvalue and wavefunction of a Kohn-Sham state in the supercell, respectively, obtained from a calculation with the ZG displacement at the temperature T ; ψ PC nk denotes a state in the primitive unit cell. We employ lower (upper) case bold fonts to indicate wavevectors in the primitive cell (supercell). The integral is over a volume that encompasses the supercell and is commensurate both with K and k. The equivalence between Eq. (65) and the standard definition of the electron spectral function using the Lehmann representation is demonstrated in Ref. [79].
We now expand Kohn-Sham states in a planewaves basis set as follows: where V is the volume where the wavefunctions are normalized. By replacing these expansions inside Eq. (66) and using the resolution of identity n c PC nk (g)c PC, * nk (g ) = δ gg one obtains: We note that in the last expression only the planewaves coefficients of the supercell state appear; therefore the calculation of the spectral function at finite temperature does not require explicit projections onto the states of the primitive cell. In actual calculations we select a k-path in the Brillouin zone of the primitive cell, and for each kpoint we proceed as follows: we identify all the supercell K-points that unfold into k via a reciprocal lattice vector G of the supercell; for each of these points we evaluate the weights P mK,k (T ) using Eq. (67); then we use the calculated weights inside Eq. (65). In the case of ultrasoft pseudopotentials, which we employed for calculations on monolayer MoS 2 in Sec. VII D, we use a slightly modified version of Eq. (67) to account for the augmentation charge [70]. Starting from the spectral function A k (ε; T ), we extract the renormalized band structure by numerically identifying the quasiparticle peaks along the energy axis.
In principle one could also evaluate a momentum-and band-resolved spectral function, by considering scalar products like in Eq. (66) but without the summation over the states n of the primitive unit cell. In this case Eq. (67) must be replaced by a more complicated expression which requires an explicit evaluation of wavefunctions in the primitive unit cell. We explored this possibility, but we found that in order to achieve numerical convergence one would need an impractically large number of bands in the supercell.

VII. SUMMARY OF PROCEDURE AND NUMERICAL RESULTS
In this section we summarize the procedure for the calculation of temperature-dependent band structures using the ZG displacement, and we discuss our results for silicon, gallium arsenide, and monolayer molybdenum disulfide. The special displacement method consists of the following steps: (i) We compute the vibrational eigenmodes e κα,ν (q) and eigenfrequencies ω qν in the unit cell using DFPT on a coarse uniform grid of q-points.
(ii) In preparation for calculations on a supercell of size N 1 ×N 2 × N 3 , we interpolate the vibrational eigenmodes and frequencies on a finer q-points grid with the same size as the supercell, N 1 ×N 2 × N 3 . We partition this grid into the sets A, B, and C.
(iii) We enforce a smooth Berry connection for the vibrational eigenmodes. To this aim we perform unitary transformations that makes eigenmodes at nearby q-points as similar as possible, as described in Sec. VI C.
(iv) We build an N 1 ×N 2 × N 3 supercell with the atoms displaced by ∆τ pκ given in Eq. (2). The signs S qν in Eq. (2) are determined using the procedure described in Sec. V D.
(v) To compute band structures, we first set up the k-point path in the Brillouin zone of the primitive cell, then we obtain the folded k-points in the Brillouin zone of the supercell. We perform a DFT band structure calculation in the supercell, and we unfold the result using the method of Ref. [77]. This procedure is discussed in Sec. VI.
In this manuscript we focus on temperature-dependent band structures to limit the length of the presentation. For calculations of temperature-dependent optical spectra, photoelectron spectra, tunneling spectra, or transport coefficients, steps (i)-(iv) above remain the same, while step (v) shall be replaced by the calculation of the required property. We emphasize that the ZG displacement in Eq. (2) does not contain eigenmodes with q ∈ A. These eigenmodes correspond to stationary lattice waves, and break the crystal symmetry. In the following sections we show how this observation can be used to analyze the convergence and the accuracy of the calculations as a function of supercell size.
The main difference between Eq. (2) and the displacement provided in our previous work [Eq. (5) of Ref. [26]] is that here a more structured choice of signs allow us to better control the convergence rate and to perform accurate calculations using relatively small supercells.
A. Thermal mean-square displacements of Si, GaAs, and MoS2 In Fig. 5(a) we compare the probability distribution P (∆τ pκα ; T ) in Eq. (16) evaluated for silicon at T = 0 K with an histogram of the displacements obtained from the ZG formula in Eq. (2). The histogram is obtained numerically by binning the values of ∆τ pκα for all atoms along the [100] direction. It is remarkable that the distribution provided by the ZG displacement follows the exact thermodynamic average with very high precision. The choice of the Cartesian direction is not important in this case, since silicon is isotropic. In fact the inset of Fig. 5(a) shows that the distribution in the (100) plane is also isotropic. For completeness we also show in Fig. 5(b) how the ZG displacement appears in a three-dimensional rendering.
In Fig. 6(a)-(c) we show the thermal mean-square displacements of Si, GaAs, and MoS 2 calculated using the ZG displacements (colored disks), the mean-square displacements evaluated using the exact expression in Eq. (17) (grey disks), and experimental data from XRD where available (triangles) [80][81][82][83]. We can see that the ZG displacement yields mean-square displacements in excellent agreement with Eq. (17), and that the agreement between our calculations and experiments is also very good. These successful comparisons reinforce the notion that the ZG displacement provides a very accurate classical representation of a thermodynamic average over the quantum fluctuations of the atomic positions.
The ZG displacement can also be employed to obtain the thermal displacement ellipsoids and compare the complete ADP tensor U κ,αβ (T ) of Eq. (18) with experiments. In Fig. 6(d)-(f) we show the computed thermal ellipsoids of Si, GaAs, and MoS 2 at T = 300 K. The ellipsoids of Si, Ga, and As are found to be spheres with radius 0.49Å 2 , 0.59Å 2 , and 0.50Å 2 , respectively. These findings are consistent with the cubic symmetry of the Si and GaAs lattices. In the case of MoS 2 , the ellipsoids reflect the two-dimensional nature of the material: the in-plane parameters U are 0.22Å 2 and 0.42Å 2 for Mo and S atoms, respectively; the out-of-plane parameters U ⊥ are 0.82Å 2 and 0.86Å 2 for Mo ans S, respectively. In all cases the ADPs obtained from the ZG displacements are within 25% of the corresponding experimental values [80][81][82][83].
B. Temperature dependent band structure of Si Figure 7 shows our results and analysis for the band structure of Si at finite temperature. Full computational details, including the evaluation of the spectral function and band energies in the primitive unit cell, are provided in Sec. VI; here we only mention that the calculations are based on DFT in the local density approximation (LDA). In Fig. 7(a), we show the spectral function A k (ε; T ) obtained from the ZG displacement in a 8×8×8 supercell. We consider T = 0 K to focus on the effect of zeropoint renormalization. The spectral function provides the momentum-resolved electronic density of states [84], and it is shown as a colormap. In Fig. 7(b) we compare the bands ε nk (T ) extracted from the spectral function A k (ε; T ) with the usual DFT band structure evaluated at clamped ions; in particular we consider the dispersions along the LΓX path of the Brillouin zone for T = 0 K (blue) and T = 300 K (green); the DFT bands at clamped ions are in black. This calculation indicates that, as a result of zero-point motion, the energy of the valence band maximum (VBM) at Γ increases by 35 meV, while the energy of the indirect conduction band minimum (CBM) at 0.87 ΓX decreases by 22 meV. These values are in excellent agreement with Ref. [12], which obtained 35 meV and 21 meV, respectively, when using the perturbative Allen-Heine approach and the adiabatic approximation.
From the band structure in Fig. 7(b) we can also extract the phonon-induced mass-enhancement. For simplicity we focus on the longitudinal electron mass of silicon, and we define the coupling strength λ T such that m * l (T ) = (1 + λ T )m * l where m * l = 0.95 m e is the DFT mass at clamped ions. From the calculated bands we obtain λ T = 0.03 and 0.04 for T = 0 K and T = 300 K, respectively. This finding is in agreement with the mass renormalization reported in Ref. [85].
Our calculated band gap narrowing due to zero-point effects is 57 meV. This value is in very good agreement with previous calculations based on non-perturbative adiabatic approaches, yielding 56-65 meV [15,19,23,25,26] as well as with experimental values, in the range of 62-64 meV [86,87]. Ref. [12] showed that non-adiabatic corrections within the Allen-Heine theory increase the zero-point renormalization by 8 meV [12]. This effect is not captured by the present special displacement method, which is in essence an adiabatic approach. We also point out that the most recent GW calculations using an earlier version of the present approach [23] yield a slightly larger gap renormalization of 66 meV. This is expected given that the electron-phonon interaction is overscreened in DFT/LDA due to the band gap problem.
In order to analyze the convergence behavior of the SDM, in Fig. 7(c) we plot the dependence of the zeropoint band gap renormalization on the supercell size. Two sets of data are shown. The green lines and datapoints correspond to calculations performed using Eq. (2). The black lines and grey datapoints were obtained after modifying Eq. (2) to include the contributions of phonons with q ∈ A. Phonons with q ∈ A correspond to stationary waves in the primitive unit cell (e.g. q = 0 phonons). As we demonstrate in Sec. V D, in the thermodynamic limit of infinite supercell the con-tribution of these modes vanishes identically. Therefore, the calculations performed by including or excluding phonons with q ∈ A in the ZG displacement given by Eq. (2) converge to the same limit. However, by considering only q ∈ B phonons as in Eq. (2) we reach convergence from the bottom (in terms of magnitude); while by including phonons with q ∈ B and q ∈ A we reach convergence from the top. By construction the converged answer must lie in between these two limits, therefore the analysis presented in Fig. 7(c) can be used as a way to quantify the convergence error of the calculations. For example the data obtained for a 10×10×10 supercell in Fig. 7(c) show that the fully converged result for an infinitely-large supercell will be in the interval 55-65 meV. This observation carries general validity for all the systems considered in this work.
The inset of Fig. 7(c) shows the spectral function near the threefold degenerate VBM of silicon, as computed using the ZG displacement for a 4×4×4 supercell at T = 0 K. If we include q ∈ A points in the ZG displacement, then the crystal symmetry is broken, and by consequence we observe a band splitting (black line). In contrast, when we use the pure ZG displacement from Eq. (2), i.e. without phonons with q ∈ A, the band degeneracy is correctly preserved. This analysis indicates that, when performing electron-phonon calculations using non-perturbative supercell approaches, q = 0 phonons as well as all phonons with q ∈ A are the least representative since they break crystal symmetry, which may lead to calculation artifacts. This issue is resolved by the present formulation of the ZG displacement as provided by Eq. (2).
In Fig. 7(d) we compare our calculations of the indirect band gap of silicon using the SDM with experiments [88], up to T = 500 K. To facilitate comparison we scissorshifted our DFT/LDA results by 0.73 eV, which is close to the GW correction reported in Ref. [89]. The agreement between our calculation and experiment is very good, except that we underestimate slightly the temperature slope. This effect is a well-known consequence of the fact that the strength of the electron-phonon interaction is underestimated by DFT/LDA; the slope can be improved by using GW calculations in combination with the SDM, as demonstrated in Ref. [23].
C. Temperature dependent band structure of GaAs Figure 8 shows our calculated band structure of GaAs at finite temperature. Also in this case we employ DFT and the LDA, as described in Sec. VI. In Fig. 8(a) we have the spectral function A k (ε; T ) at T = 0 K as a color map, and in Fig. 8(b) we have the bands extracted from the spectral function at T = 0 K (blue) and T = 300 K (green). The bands at clamped ions are shown in black for comparison. For simplicity we are not including spinorbit coupling in the calculations, therefore we do not see the spin-orbit splitting of the split-off holes in the valence bands [90][91][92]. Since the holes at the top of the valence band have the same orbital character, we expect a similar zero-point renormalization for the split-off holes.
From the data in Fig. 8(b) we obtain zero-point corrections to the valence and conduction band edges at Γ of +21 meV and −11 meV, respectively. The resulting band gap narrowing, ∆E g = 32 meV, lies within the experimental range of 57 ± 29 meV [93]. The calculations in Fig. 8(b) are performed using an 8×8×8 supercell. To better compare with the finite-differences results of Ref. [10], we repeated the calculations with a smaller, 4×4×4 supercell. In this case we obtain ∆E g = 25 meV, which matches the value of 25 meV reported in Ref. [10]. It is well established that GW quasiparticle corrections change these results by strengthening the electron-phonon coupling [10,23]: in order to incorporate this effect it will be sufficient to perform a GW calculation using the ZG displacement.
From the band structure in Fig. 8(b) we can determine the phonon-induced mass-enhancement. Since we are not including spin-orbit effects, we only report the value for the CBM at Γ. By denoting m * e = 0.056 m e the conduction band mass at clamped ions, we find m * e (T ) = (1 + λ T )m * e with λ T = 0.005 and 0.007 for T = 0 K and T = 300 K, respectively. We are unaware of previous calculations of the mass renormalization in GaAs as a function of temperature.
In Fig. 8(c) we show the convergence of the zero-point renormalization of the direct gap of GaAs with respect to the supercell size. The converged value obtained with Eq. (2) is ∆E g = 32 meV and is obtained using an 8×8×8 supercell. As for the case of Si in Fig. 7(c), also for GaAs the threefold degeneracy of the VBT is preserved when using the ZG displacement. And also in this case, if we include phonons with q ∈ A in Eq. (2), the degeneracy is lifted due to symmetry breaking, as it can be seen in the inset of Fig. 8(c).
In previous work it has been suggested that in the case of polar materials it might be impossible to achieve convergence within the adiabatic approximation [12]. Here we wish to emphasize that the divergence of perturbative approaches for polar materials is not a consequence of the adiabatic approximation, but it results from taking the limit ω qν → 0 before performing the integration over phonon wavevectors to obtain the Fan-Migdal self-energy [1]. When the limit ω qν → 0 is correctly performed after evaluating the integral over q, there is no divergence. The correct procedure can be found in Sec. IV of the early work by Fan [94] and is summarized in Appendix C. Figure 8(c) shows indeed that adiabatic calculations using the SDM converge smoothly as a function of supercell size, and that the fully converged band gap renormalization lies in a very narrow bracket between 32 meV (ZG displacement without q ∈ A phonons) and 36 meV (ZG displacement with q ∈ A phonons) already for a 8×8×8 supercell.
Generally speaking, the smooth and fast convergence of the SDM can be ascribed to the fact that the formalism relies on a standard DFT calculation for a supercell with displaced atoms. This calculation is intrinsically easier to converge than perturbative approaches. In fact perturbative methods require stringent q-point sampling to evaluate principal-value integrals of first-order poles that appear in the Fan-Migdal and Debye-Waller self-energies; furthermore these integrals yield large and opposite contributions, so the final result is obtained by subtracting two large numbers. The SDM method circumvents these numerical challenges.
In Fig. 8(d) we show our results for the direct gap of GaAs using the ZG displacement, in the temperature range 0-500 K, and we compare with experimental data from Ref. [93]. In order to take into account the non-negligible expansion of the GaAs lattice with temperature, we employ the quasi-harmonic approximation. We use a scissor-shift of 0.53 eV to mimic GW corrections as reported in Ref. [95]. The slight underestimation of the experimental slope with temperature can be corrected by performing GW calculations instead of standard DFT/LDA [10,23], but overall the agreement between the present calculations and experiments is very good.
For completeness, we also mention that the results shown in Fig. 8(d) and obtained from the analysis of temperature-dependent band structures are in excellent agreement with the values that we obtained in an earlier work by analyzing the joint density of states [26].
D. Temperature dependent band structure of monolayer MoS2 Figure 9 shows our calculated temperature-dependent band structure of monolayer MoS 2 , as an example of application of the SDM to two-dimensional materials. In this case we employed the PBE exchange and correlation functional [68], and fully relativistic ultrasoft pseudopotentials [70] to take spin-orbit coupling into account (Sec. VI).
In Fig. 9(a) we have the spectral function of MoS 2 at T = 0 K obtained from the ZG displacement in a 10×10×1 supercell. By extracting the corresponding temperature-dependent bands we obtain a spin-orbit splitting of 136 meV for the valence band states at K, to be compared to the clamped-ion value of 142 meV. Our calculation is in agreement with the spliting of 135 meV obtained in Ref. [14] within the Allen-Heine theory. We also determined the electron effective mass renormalization at the K point, and found λ T 0.04 and 0.06 for T = 0 K and T = 300 K, respectively. These latter two values are not fully converged since we used finite differences with a coarse k-point spacing of 6 · 10 −3 2π/a. In Fig. 9(b) and (c) we show convergence tests for the gap renormalization. In Fig. 9(b) we have the renormalization as a function of interlayer separation by keeping the in-plane supercell size fixed; in (c) we vary the supercell size, keeping the interlayer separation constant. As in the case of the three-dimensional materials Si and GaAs, the band gap narrowing converges relatively quickly with the size of the cell; furthermore the calculations are not very sensitive to the size of the vacuum region provided periodic replica are separated by more than 10Å.
For a 10×10×1 supercell we obtain a zero-point gap renormalization of 64 meV using the ZG displacement. By including also q ∈ A phonons in Eq. (2) the value changes only slightly to 65 meV, therefore we estimate an error with respect to the fully converged result of less than 1 meV. Our result is in good agreement with the perturbative calculations of Ref. [14], which reported 75 meV. The residual difference may be due to the 'rigidion' approximation used in the Allen-Heine approach of Ref. [14], and to differences in the DFT exchange and correlation functionals employed.
In Fig. 9(d) we compare our calculated temperaturedependent band gap of MoS 2 (green) with experimental data from Ref. [96] (black). To facilitate comparison we introduced a scissor correction of 0.34 eV to match the measured band gap at 4 K. This correction is similar to that obtained using GW and the Bethe-Salpeter equation (BSE) [97]. Unlike in Si and GaAs, here the calculated data follow the experimental trend very closely. This is unexpected since we are computing electron-phonon effects at the DFT/PBE level, and suggests that quasiparticle corrections of the electron-phonon coupling [98] are not significant in MoS 2 .

E. General remarks on the SDM results
The applications to Si, GaAs, and MoS 2 described in the previous sections show that the SDM yields temperature-dependent bands of the same quality as perturbative approaches based on the Allen-Heine theory. In particular, the capability to access momentum-resolved quantities such as temperature-renormalization at specific k-points and the phonon-induced enhancement of the band mass are entirely novel, and considerably expand the range of applicability of the ZG displacement.
One important point of note is that the present reciprocal-space formulation of the ZG displacement allowed us to identify and remove the contributions of vibrational modes that break crystal symmetry. This improvement is important in order to avoid spurious and uncontrolled lifting of electronic degeneracy. Since the symmetry-breaking modes are those with wavevector q ∈ A, that is phonons corresponding to standing Bloch waves, the present analysis and results demonstrate that these modes (especially q = 0 phonons) are the least representative in a thermodynamic average, and should be avoided for accurate calculations.
Another interesting point is that we can bracket the fully converged results by performing two calculations: one with the pure ZG displacement of Eq. (2), and one with the modified version including q ∈ A phonons. This provides a simple and effective strategy to quantify the convergence error as a function of supercell size.
We also found that, when used in conjunction with the ZG displacement, the adiabatic approximation does not suffer from the convergence problems that are encountered in perturbative approaches for polar materials. This advantage arises from the fact that the special displacement method does not require integrating over poles as in perturbative approaches, and the calculation is as easy as a standard calculation at clamped ions.
For the systems considered in this work, the supercells required to achieve relatively accurate results are surprisingly small. For example, the ZG displacement in a 4×4×4 Si supercell (128 atoms) yields a zero-point renormalization which differs by less than 10% from the fully converged value; for GaAs we need a 6×6×6 supercell (432 atoms) to achieve similar accuracy; for MoS 2 a 6×6×1 supercell (72 atoms) is enough to converge the results within 10%. This indicates that the ZG displacement can be used to perform calculations with relatively small supercells, and this may open the way to post-DFT calculations at finite temperature, including GW quasiparticle calculations and exciton calculations via the BSE method.

VIII. CONCLUSIONS AND OUTLOOK
In this manuscript we develop a new methodology for performing electronic structure calculations at finite temperature. In a nutshell, this method consists of performing a single calculation in a supercell where the atoms have been displaced according to Eq. (2). We refer to this displacement as the ZG displacement, and to the methodology as the special displacement method.
This work follows our earlier study in Ref. [26], where the original idea was first proposed. The key novelty of the present study is that we reformulate the entire theory using a compact and rigorous reciprocal space formulation, building on density functional perturbation theory. This new formulation allows us to go beyond angle-integrated electronic and optical spectra, and to compute complete band structures at finite temperature. We demonstrate this concept for three-dimensional bulk semiconductors, silicon and gallium arsenide, and for a two-dimensional semiconductor, monolayer molybdenum disulfide. In all cases our results match the accuracy of established perturbative techniques based on the Allen-Heine theory. The added advantage of the present approach over perturbative methods is that it does not require the evaluation of self-energy energy poles, and it does not require the rigid-ion approximation for the Debye-Waller self-energy. As a consequence, the method is robust, reliable, and simple to use.
Beyond demonstrating calculations of band structures at finite temperature, we show that the ZG displacement accurately reproduces the characteristic anisotropic displacement parameters measured by X-ray crystallography, and can be used to determine thermal ellipsoids as a function of temperature. Therefore the ZG displacement represents an accurate classical snapshot of thermal disorder in solids, and eliminates the need for the stochastic sampling of the nuclear wavefunctions.
More fundamentally, we show that the ZG displacement can be understood as the best single-point approximant to a thermodynamic Feynman path-integral. This link may open new avenues in the study of path integrals using the quantum-classical analogy and special displacements.
In the present approach, the choice of the atomic displacements is carefully performed by first establishing a smooth Berry connection among vibrational eigenmodes in the Brillouin zone, and then by choosing the phase of each eigenmode using rigorous sum rules. This approach allows us to prove that the ZG displacement yields the exact Williams-Lax average of an electronic observable in the thermodynamic limit of infinite supercell. We also provide procedures for accelerating the convergence of the calculations for those cases where only small supercells are within reach.
For reasons of space we only discussed applications to finite-temperature band structures, but we emphasize that this methodology is much more general. For example the earlier version of this method [26] has already been applied to compute phonon-assisted optical properties, dielectric functions, GW quasiparticle corrections, zero-point renormalization of band gaps, exciton-phonon couplings, and charge transport. This broad range of applications is possible because the Williams-Lax theory can be employed to compute any property which can be expressed as a Fermi's golden rule.
The present extension of the special displacement method to compute entire band structures is particularly appealing for testing quasiparticle corrections at finite temperature. Since only one supercell calculation is required, and since we developed a new algorithm to accelerate convergence with respect to the supercell size, this approach opens the way to systematic many-body calculations of band structures at finite temperature.
Up to this point the SDM relied on the harmonic approximation. It would be interesting to consider extensions to anharmonic systems. We expect that the method will work seamlessly in conjunction with the quasi-harmonic approximation [99,100] and with the selfconsistent harmonic approximation [101,102]. In fact, in both cases the original anharmonic potential is replaced by its 'best' harmonic approximation; after this replacement the SDM method can be used without changes. In the case of strongly anharmonic systems, for example in the presence of double-well potential energy surfaces, it should be possible to modify the present method by treating all the harmonic modes as described in this manuscript, and by adapting the ZG configuration to describe averages over double-well quantum nuclear wavefunctions. The feasibility and accuracy of this approach will have to be demonstrated in future work. Here we summarize the relations between normal modes e κα,ν (q) obtained from the diagonalization of the dynamical matrix at each q wavevector, and the standard sum rules resulting from the Born-von Kármán boundary conditions.
In the harmonic approximation the dynamical matrix D κα,κ α (q) at each q-point and the vibrational eigenmodes and frequencies satisfy the equation: [1] κ α D κα,κ α (q)e κ α ,ν (q) = ω qν e κα,ν (q). (A1) Since the dynamical matrix is Hermitian, the eigenmodes form an orthonormal basis: Equation (A1) implies the relations: [1] e κα,ν (−q) = e * κα,ν (q), where we followed the convention of Ref. [55]. We denote the lattice vector pointing to the p-th unit cell as R p = n p1 a 1 +n p2 a 2 +n p3 a 3 , where the a i represent the primitive lattice vectors and n pi are integers between 0 and N i − 1. The Born-von Kármán supercell contains N p = N 1 × N 2 × N 3 unit cells. For the uniform grid of phonon wavevectors q in the Brillouin zone we choose where the b i represent the primitive reciprocal lattice vectors, and m i are integers between 0 and N i − 1. The vectors b i are such that a i · b j = 2πδ ij . With these conventions we have the sum rules: Appendix B: Link between Williams-Lax theory and Feynman's path integrals In this appendix we provide some of the mathematical background required to link the Williams-Lax theory with thermodynamic path integrals. The following equations are used in Sec. IV.
Starting from Eq. (20), using the Trotter decom- [103] and inserting the resolutions of identity dτ i |τ i τ i | = 1 with i = 1, · · · , N −1 in between each product, Eq. (20) becomes: where we introduced |τ 0 = |τ N = |τ to make the notation more compact. The nuclear Hamiltonian can be expressed as a sum of kinetic (T ) and potential (Û α ) energies, and the exponential in Eq. (B1) is rewritten for large N using the Baker-Campbell-Hausdorff formula [104]: In this limit the matrix elements appearing in Eq. (B1) take the form: having defined ∆t = /N k B T . The matrix element with the kinetic energy is simplified by inserting the resolution of identity for the eigenstates of the nuclear momentum operators, dk|k k| = 1, where {k} represent a complete set of planewaves associated with the atomic positions {τ }. After evaluating terms like k|τ i and solving the resulting Gaussian integral, this procedure leads to: where c is a normalization prefactor. By taking the limit of large N , which is equivalent to taking ∆t → 0, the terms (τ pκα,i+1 − τ pκα,i )/∆t become the classical velocities of the nuclei, and the exponent of Eq. (B4) yields the classical kinetic energy of the nuclei.
The combination of Eqs. (B1)-(B4) with Eq. (20) yields the link of the Wiliams-Lax transition rate to a thermodynamic Feynman path integral, as shown in Eq. (21). divergence of the Fan-Migdal self-energy is an artifact of the integration procedure. Mathematically the divergence arises from collapsing the two imaginary poles at ±iq LO into a double pole at the origin. A similar problem arises in the textbook Fourier transform of the Coulomb potential.
In practical ab initio calculations the adiabatic approx-imation can be retained without incurring into a divergence as follows. First we perform calculations where all the phonon frequencies are set to a small constant, say ω 0 . After converging the summation over the Brillouin zone, we repeat the calculations for smaller values of ω 0 , so as to take the adiabatic limit ω 0 → 0. This procedure will yield a finite self-energy.   The same eigenmodes as in (a), this time after using the algorithm described in Sec. VI C to enforce a smooth Berry connection between eigenmodes at adjacent q-points. In this case the eigenmodes vary smoothly along the Brillouin-zone path.  6. (a) Mean-square thermal displacements of silicon as a function of temperature, evaluated using Eq. (17) (gray disks) and the ZG displacement (green disks). Experimental data from Ref. [80] are shown as black triangles. (b) Thermal ellipsoids of silicon at T = 300 K, as obtained from the ZG displacement. (c) Mean-square thermal displacements of GaAs as a function of temperature, evaluated using Eq. (17) (gray disks) and the ZG displacement (green and blue disks). We also report experimental data from Ref. [81] (filled triangles) and Ref. [82]. These data correspond to weighted averages of the displacements of Ga and As. (d) Thermal ellipsoids of GaAs at T = 300 K, as obtained from the ZG displacement. (e) In-plane mean-square thermal displacements of monolayer MoS2 as a function of temperature, evaluated using Eq. (17) (gray disks) and the ZG displacement (green and blue disks). We also report experimental data from Ref. [83]. (f) Thermal ellipsoids of monolayer MoS2 at T = 300 K, as obtained from the ZG displacement. (a) Spectral function of silicon calculated using the ZG displacement at T = 0 K. The calculation was performed using an 8×8×8 supercell and the unfolding procedure described in Sec. VI D. We sampled the L-Γ-X path on 280 equally-spaced k-points, and the zero of the energy axis is referred to the valence band top. (b) Band structures of silicon at clamped ions (black), T = 0 K (blue), and T = 300 K (green). The bands were extracted from spectral functions like the one in (a). We also report the calculated zero-point renormalization (ZPR) and a ball-stick model. (c) Sensitivity of the calculated ZPR of silicon to the size of the supercell. The horizontal axis indicates the linear size N of the N×N×N supercell. We show both the calculations performed using the ZG displacement (green), and the results obtained by also including q-points in set A (grey). In the latter case the threefold degeneracy of the valence band top is lifted (inset), and the band gap is evaluated by considering the topmost valence state. (d) Temperature dependence of the indirect band gap of silicon up to 500 K. We show the results of the special displacement method (green circles) and experimental data from Ref. [88] (black triangles). The calculated band gaps were scissor-shifted by 0.73 eV to match the experimental value at 4 K. The straight line is the high-temperature limit and intercepts the T = 0 K axis at the clamped-ion band gap (1.23 eV, empty circle). (a) Spectral function of GaAs calculated using the ZG displacement at T = 0 K. The calculation was performed using an 8×8×8 supercell and the unfolding procedure described in Sec. VI D. We sampled the L-Γ-X path on 280 equally-spaced k-points, and the zero of the energy axis is referred to the valence band top. (b) Band structures of GaAs at clamped ions (black), T = 0 K (blue), and T = 300 K (green). The bands were extracted from spectral functions like the one in (a). We also report the calculated zero-point renormalization (ZPR) and a ball-stick model. (c) Sensitivity of the calculated ZPR of GaAs to the size of the supercell. The horizontal axis indicates the linear size N of the N ×N ×N supercell. We show both the calculations performed using the ZG displacement (green), and the results obtained by also including q-points in set A (grey). In the latter case the threefold degeneracy of the valence band top is lifted (inset), and the band gap is evaluated by considering the topmost valence state. (d) Temperature dependence of the indirect band gap of GaAs up to 500 K. We show the results of the special displacement method (green circles) and experimental data from Ref. [93] (black triangles). The calculated band gaps were scissor-shifted by 0.53 eV to match the experimental value at 25 K. The straight line is the high-temperature limit and intercepts the T = 0 K axis at the clamped-ion band gap (1.56 eV, empty circle). The horizontal axis indicates the linear size N of the N ×N ×1 supercell. We show both the calculations performed using the ZG displacement (green), and the results obtained by also including q-points in set A (grey). (d) Temperature dependence of the indirect band gap of monolayer MoS2 up to 400 K, evaluated using a 10×10×1 supercell. We show the results of the special displacement method (green circles) and experimental data from Ref. [96] (black triangles). The calculated band gaps were scissor-shifted by 0.34 eV to match the experimental value at 4 K. The straight line is the high-temperature limit and intercepts the T = 0 K axis at the clamped-ion band gap (1.93 eV, empty circle).