Phonons from density-functional perturbation theory using the all-electron full-potential linearized augmented plane-wave method FLEUR

Phonons are quantized vibrations of a crystal lattice that play a crucial role in understanding many properties of solids. Density functional theory provides a state-of-the-art computational approach to lattice vibrations from first-principles. We present a successful software implementation for calculating phonons in the harmonic approximation, employing density-functional perturbation theory within the framework of the full-potential linearized augmented plane-wave method as implemented in the electronic structure package FLEUR. The implementation, which involves the Sternheimer equation for the linear response of the wave function, charge density, and potential with respect to infinitesimal atomic displacements, as well as the setup of the dynamical matrix, is presented and the specifics due to the muffin-tin sphere centered linearized augmented plane-wave basis-set and the all-electron nature are discussed. As a test, we calculate the phonon dispersion of several solids including an insulator, a semiconductor as well as several metals. The latter are comprised of magnetic, simple, and transition metals. The results are validated on the basis of phonon dispersions calculated using the finite displacement approach in conjunction with the FLEUR code and the phonopy package, as well as by some experimental results. An excellent agreement is obtained.


Introduction
Phonons are quantized collective lattice vibrations featuring a discrete spectrum of frequencies.They are also described and understood as quasiparticles using the framework of quantum field theory.The basic theory of phonons is well understood and has been described in detail in text books [1,2].In the harmonic approximation, the phonon frequencies ω are determined by the eigenvalues det D − ω 2 1 = 0 (1) of the dynamical matrix (DM) D. 1 denotes the identity matrix.The dynamical matrix, here presented in terms of matrix elements with the atomic masses M γ( ′ ) is the reduced Fourier transform of the harmonic force constant matrix ϕ, also known as Hesse matrix of second-order derivatives E tot with matrix elements that describes consistent with the harmonic approximation the second-order expansion of the Born-Oppenheimer energy E tot with respect to the positions τ of the atoms γ, and γ ′ in some unit cells.τ i with i ∈ {1, 2, 3} are the cartesian components of τ .γ ′ R describes the atom γ ′ in the unit cell with lattice vector R at τ γ ′ R = τ γ ′ + R. The symbol q denotes the phonon wave vector defined within the Brillouin zone (BZ) of the crystal lattice.It represents the momentum associated with a phonon and describes the propagation direction and the wavelength of the lattice wave.The dimension of the hermitian dynamical matrix scales with the dynamical degrees of freedom of the lattice, i.e. the number of atoms in the unit cell, N A , along the three cartesian coordinates as dim(D) = 3N A × 3N A .The solutions of (1) are displacement modes q)e i(q•R+ωµ(q)t) + P * γ µ (q)e −i(q•R+ωµ(q)t) , where the normalized polarization vector P γ µ ∈ C 3 and w γ R µ (0) denotes the direction and the arbitrary amplitude of the displacement for the 3N A phonon modes µ of atom γ in unit cell R with wave vector q in Brillouin zone of volume Ω BZ .
Phonons play a crucial role in understanding a vast number of material phenomena.They lie at the heart of thermodynamical properties of solids, heat and sound propagation, reveal elastic properties of materials, and contribute to electrical resistivity [1,[3][4][5].In conventional superconductors, the interaction of electrons with phonons is the primary mechanism responsible for the attractive pairing of electrons leading to the superconducting state [6].Phonons are of continuous interest due to their role in engineering acoustic metamaterials [7,8], as driving force for charge-density waves [9], for the optimization of the phonon transport in thermoelectrics [10], and in the context of ferroelectric [11], 2D [12][13][14], and magnetic [15] materials.In the latter they contribute to spin-relaxation [16], Gilbert damping [17] and -equilibration [18], assist magnetization switching by linearly- [19] and circularly-polarized [20], or chiral phonons [21], influence the temperature dependence of the magnetocrystalline anisotropy [22], and can be of interest in the fields of orbitronics [23] and thermal Hall physics [24].
Several theoretical approaches are employed to evaluate phonon properties in condensed matter systems [25].Among those, the Kohn-Sham (KS) density functional theory (DFT) [26][27][28][29][30] has established itself as the method of choice for providing materials specific information directly from the electronic structure without adjustable parameters.Concerning the computation of phonons and related quantitites, there are basically two established DFT approaches in use: (i) the finite displacement (FD) method [31][32][33], and (ii) the density-functional perturbation theory (DFPT) [34][35][36][37][38].The FD method emerged first and has been preferred in a wide spectrum of the literature to this day (see e.g.[39]), while the number of publications using DFPT is constantly increasing.FD and DFPT are complementary to each other and, given the same input, are able to deliver quantitatively comparable results.They are often applied in parallel, e.g. to test the degree of unharmonicity in teh FD results.Beneficial for both methods is the 2n + 1-theorem [36], which gives access to quantities of order 2n + 1, while only having input quantities of order n at hand.Both methods do not deliver continuous dispersion relations, which is why usually interpolation methods are performed as a postprocessing step.However, the DFPT method excels in improving the interpolation at specific points of interest in the Brillouin zone, since it offers access to them at affordable numerical costs.A comprehensive overview is given in the review of Baroni [37] or the book of Martin [40].
In the FD method, one takes advantage that the force acting on an atom due to a displacement is the first derivative of the total energy, F = −∂E tot /∂τ , and the force-constant matrix elements are calculated by a difference quotient of the i-th Cartesian component of the force F γi acting on atom γ when another atom γ ′ of the solid is displaced by a small displacement ∆τ γ ′ i ′ from the equilibrium position τ γ ′ into direction i ′ .At equilibrium F γi (τ γ ′ i ′ ) is usually zero.For crystalline solids, on which we focus throughout the paper, symmetry can usually be exploited, reducing the number of necessary force-vector components and displacements.Additionally, the combinations of (γ ′ , γ) reduce to the N 2 A pairs of atoms (β, α) in the representative unit cell.This is normally automated by software packages such as phonopy [41,42], which provide phonon calculations at harmonic and quasiharmonic levels.Overall, the implementation of the FD method is quite simple, provided the DFT code delivers reliable forces.Nevertheless, the supercells must be chosen to include different periods determined by the phonon vector q.As a consequence, the qvector must be commensurate to the supercell, restricting this method to rational qvector components, and making the calculation of phonons with a q-vector exhibiting a small absolute value very expensive.
In DFPT, we take an analytical approach to the second derivates of the total energy.Then, the DM contains, among other terms, linear responses of the charge density and the effective potential to the change of the external potential caused by the phonon.In DFPT, the first-order response functions of the electronic structure to small perturbations of the atom positions without the need to perform completely new calculations for each perturbation are calculated in a self-consistent way using the Sternheimer equation [43], which is a first-order version of the KS eigenvalue equation.As will be outlined in subsection 2.2, this gives access to E (2) tot avoiding supercell calculations.The costs of a DFPT calculation are equally distributed among arbitrary q-vectors and comparable to a DFT self-consistency procedure.
Most phonon studies using DFPT have been performed with norm-conserving pseudopotentials [44][45][46][47], but there are now also a plethora of studies using ultrasoft (US) pseudopotentials [48][49][50][51] and the projector-augmented wave (PAW) method [52,53].Publications and codes combining DFPT and all-electron muffin-tin based electronic structure methods such as the augmented spherical wave (ASW) [54], linear muffin-tin orbital techniques (LMTO) [55], Korringa-Kohn-Rostoker (KKR) Green function [56], and full-potential linearized augmented plane-wave (FLAPW) method [57][58][59][60] are scarce [61][62][63] and the technicalities of the implementation are not well explored.All-electron methods treat core and valence electrons on the same footing.In order to deal with the Coulomb singularity produced by the nuclear charge and the associated rapid variation of the core and valence electron wave functions and charge densities in the vicinity of the nucleus, all-electron methods partition the space of the unit cell into muffin-tin spheres in which wave functions, charge densities and potentials are represented in real space.A Fourier representation of these quantities would hardly converge.
In this paper, we present a successfully working implementation of DFPT in the context of the all-electron FLAPW method.The FLAPW method is frequently considered a reference for electronic structure (DFT) calculations [64,65], especially when dealing with magnetism, systems with localized electrons such 2p, 3d, and 4f electrons, or open systems, and systems in lower dimensions.The FLAPW methodology is well-developed [66] and first-order changes of the total energy such as forces [67][68][69] or the stress-tensor [70] are well-established.The second-order changes, however, are at a different scale, since they also require the density response in the form of first-order changes of the density and second-order derivatives, which require greater numerical attention as differentiation acts numerically as a roughening operator.Here, we present solutions to known numerical challenges of muffin-tin based electronic structure methods in general, and the FLAPW method in particular, in the context of the DFPT approach such as: The Madelung summation, the Coulomb singularity of the potential, the rapidly varying wave functions and charge densities in the vicinity of the nucleus, the calculations of gradients of the all-electron potential, the presence of the core electrons, the incompleteness and the position dependence of the basis-set, the different representations of the basis-set in muffin-tin-spheres and the interstitial region and their match at the muffin-sphere boundary.On a more general level, this implementation allows to gain insight in response properties of highly complex materials.
We implemented our approach in the open source electronic structure package FLEUR [71,72], more precisely in the bulk version of general symmetry.In the context of this work, it is worth mentioning that an emphasis was placed on the implementation of a numerically accurate force formalism [69,73] (to which the DFPT implementation is very alike to), and on the choice of local orbitals [74] to reduce the linearization error [75,76] and to improve the LAPW basis set [77] towards unoccupied states of higher energies [78].We show that the implementation of the DFPT presented here, which is based on the dissertations of Klüppelberg [79] and Gerhorst [80], in which further nitty-gritty details can be found, provides a solid foundation for calculating the phononic properties, and charge density response properties in general of complex materials with the FLAPW method according to first principles.
This paper is organized as follows: We briefly recapitulate the central theoretical background of DFT, DFPT, and the FLAPW method, also to establish a consistent notation.We then explain the technical details of our implementation.We present the general concept and the workflow for DFPT calculations, discuss the challenges related to the choice of the LAPW basis, the implementation of the Sternheimer equation and the dynamical matrix, and solutions to the challenges.In order to guarantee a good reading flow of the paper and not to be overloaded with details, we have separated additional technical details into Appendix A to Appendix H.Although for clarity and simplicity the implementation is presented based on electronic charge density only as in the context of non-spin-polarized DFT, we also apply this method to collinear magnets by doubling the formalism and incorporating the magnetization density, thus replacing the charge density by spin-densities of spin-up and -down electrons.Finally, we validate our DFPT framework with respect to the quality of Goldstone modes and phonon dispersion relations against the FD approach on a selection of materials and conclude with an outlook to future developments.

Density Functional Theory
According to the Kohn-Sham DFT, the total energy of a system of interacting electrons is uniquely determined by its ground-state charge-density distribution and the problem of a system of interacting electrons is mapped onto an equivalent non-interacting problem with same ground-state density.This is made possible by expressing the unknown density functional in the form where the first term is the kinetic energy of the non-interacting system, the second term is the classical electrostatic self-interaction of the electron charge-density distribution, known as Hartree energy E H , the third is the unknown and well-approximated exchangecorrelation (xc) energy E xc , and the fourth term describes the interaction of electrons with the potential external to the electrons, e.g. of the nuclei positioned at τ .The final term describes the electrostatic interaction among the nuclei.This approach is in principle exact, but the aforementioned exchange-correlation energy is not known explicitly and there exists a large variety of approximations [81,82].In this paper we work with the local-density approximation [83], a simple representative of the xcfunctionals, which leads to good results for a large class of materials.For a set of atoms located at {τ }, the ground-state density is obtained through the Kohn-Sham equations that are solved self-consistently and comprise the effective potential V eff , being a functional of the ground-state charge density n (0) (r), subdivided into a sum of the external (ext), Hartree (H), and exchange-correlation potential (xc).The external and Hartree apart are often grouped as the Coulomb potential V C .Ψ ν are the eigenstates and ϵ ν the eigenenergies.The index ν (o) denotes (occupied) spin-degenerate states.
For simplicity, throughout this work the spin index is omitted.For magnetic systems, we switch to the well-established spin-density functional theory [84], where the treatment of collinear magnets is straightforward: The spin-degeneracy is lifted and (7a) is solved separately for the spin-up (↑) and -down (↓) states, Ψ ν↑(↓) , solutions of a spin-dependent potential ↓ ].The latter term is obtained by generalizing the ground-state density to the ground-state spin densities n (0) ↑(↓) , calculated via the summations (7c) of spin-up and -down states, separately, with is related to the spin-independent exchange correlation potential V xc and a magnetic exchange-correlation field ↓ ].These generalizations hold also true for the phonon calculations below.

Density-Functional Perturbation Theory
Depending on the energy scales or phenomena of interest, quantities of certain orders in a perturbation are to be determined.Given a phonon, the dynamical matrix and consequently the second-order changes in the total energy, E tot , with respect to atomic displacements ∂τ turn out to be pivotal.Applying the Hellmann-Feynman Theorem to the second order derivative of the energy E in (6) and restricting ourselves here for simplicity to one displacement component of one atom λ = τ γi the second-order change of the energy reads where the basis-set independent variation of the ion-ion interaction E (2) ion−ion is included.The integral spans over the volume of the unit cell Ω.In our nomenclature, quantities with the superscript (1) (or (2)) are defined as perturbed quantities to first (second) order, while ones without a superscript are the unperturbed quantities of the groundstate system.Unlike first-order changes of the energy, such as forces or stress tensors, for second-order changes, the terms involving the derivative of the density do not vanish.This means that it is necessary to compute the electronic response of the system to the displacement of atoms to perform ab initio lattice dynamics calculations.The requirement of a first-order density change makes the calculation of quantities requiring second-order energy derivatives qualitatively very different from the evaluation of quantities requiring only first-order energy changes.The first-order change in the density, n (1) , constitutes a key quantity of the DFPT and reads o (r) with N (1) = Ω n (1) (r) dr = 0 , (9) relating to the first-order change in the eigenfunctions Ψ o .
For our current implementation, exploiting symmetry this can be simplified further (see Appendix A).The relationship between the second order change in the energy and the first-order change in the eigenvalues and eigenstates of the underlying Hamiltonian satisfies the well-known 2n + 1-theorem [36], which states that a (2n + 1)-order derivative of the energy of some Hamiltonian can be calculated from the knowledge of the eigenfunction and its derivatives up to order n.
Access to the aforementioned first-order change of the electronic quantities is provided by the solution of the Sternheimer equation in a self-consistent fashion, as the change in the charge density creates a change in the effective potential and vice versa.Assuming non-degenerate states, the following basic form of the Sternheimer equation holds where H is the Hamiltonian and V eff is the first-order change of the effective potential, which contains not only V (1) ext but also the Hartree and exchange correlation kernel, δ(V H + V xc )/δn, and the density response n (1) .The projector onto the unoccupied subspace of states denoted by subscript u is explicitly included.This gives rise to a self-consistency calculation: The potential response determines the response of the eigenstates, which are used to calculate the density response, that in turn is used to construct the potential response.This Sternheimer self-consistency cycle is very similar to that of a DFT ground-state calculation.The generation of the density and potential is replaced by the generation of their respective responses, the starting perturbation is only that of the external potential, instead of the original Hamiltonian and overlap matrices the corresponding response matrices are set up and the density response is mixed to achieve self-consistency instead of the density itself.A key difference is that there is no diagonalization step as for the Schrödinger equation (solving the Sternheimer equation is purely matrix-vector multiplication), and that we need access to the full eigenspectrum of each k-point, not only the occupied states.After self-consistency is reached, the variational solution can be used to calculate the density response (9) and subsequently the force constant matrix (8).The ramifications of applying the formalism in the LAPW basis [77] will be explored in section 3.In summary, the DFPT for phonons requires the first-order changes in the density, of the wave function, the external and effective potential, as well as the second order changes in the external potential (whose evaluation will be avoided in practical implementation) and the ion-ion energy.In reality, all first-order changes are vector quantities and all second-order changes are matrices, and the product n (1) V (1) ext in (8) turns into a direct product n (1) ⊗ V (1) ext .Analogously to the discussions in 2.1, for collinear magnetic systems the Sternheimer equation ( 10) is solved for the changes of the spin-up and -down states used to synthesize the spin-density response (9), which sum to the required density response.In principle, (8) and (10) are sufficient in a plane-wave ansatz and together with (9) they make up the concept of this DFPT implementation.However, the position dependent, incomplete and multi-domain represented basis-set of the LAPW basis gives rise to a multitude of additional terms, each of them to be carefully taken into account.

Full-Potential Linearized Augmented Plane-Wave Method
When describing wave functions, Ψ kν (r), in a periodic lattice, their natural form are Bloch waves characterized by a crystal momentum vector k restricted to the first Brillouin zone (BZ) of the reciprocal lattice, and a band index ν.They are typically expanded into basis functions, e.g.plane waves, with reciprocal lattice vectors G, where z k+G,ν are the corresponding expansion coefficients.The maximum length of the reciprocal lattice vectors, K max , determines the number of basis functions N B and controls the numerical effort and precision of the results.Care has to be taken when selecting the k-point set to maintain the symmetry of the lattice: here, we always choose an equidistant mesh containing the Γ-point with N kx × N ky × N kz k-points in the reciprocal space.For odd N k i this corresponds to a Monkhorst-Pack mesh [85].For further details on the choice of the k-point mesh see also section 4.2.
To deal with the Coulomb singularity ∼ 1/r at the center of the atoms due to the positively charged nuclei and the rapidly oscillating core and valence electron wave functions in the vicinity of the nuclei, as typical for all-electron methods, in the FLAPW method [86], the computational domain is divided into spheres MT γ around the centers of each atom γ -the union of all these spheres is called the muffin-tin (MT) region -and into an interstitial (IR) region.The basic plane-wave approach is kept in the interstitial region of the unit-cell with volume Ω, but it is augmented by radial functions u ℓ (r) and spherical harmonics Y L (r) with the angular momentum and magnetic quantum numbers L = (ℓ, m), and the unit vector r = r/r in the MT region.To guarantee sufficient variational flexibility of these LAPW basis functions, different "orders" (denoted by the index p) of radial functions are used.The zeroth order u ℓ0 (r) corresponds to the solution of the radial Schrödinger equation for a spherical potential (containing the full atomic ∼ 1/r singularity) in the MT spheres to a given energy parameter characteristic of the valence electrons, and the first order functions u ℓ1 (r) correspond to their first order energy derivatives [77].Extending this logic, in certain cases, we supplement the LAPW basis [77] with local orbitals [78,87].These are used to give more variational freedom, to accurately describe semicore states, high-lying unoccupied states, as well as to reduce the linearization error [76].These are exclusively present in the MT.The LAPW basis functions are thus with the unit cell volume Ω and coefficients a k+G,γ

Lp
, with ℓ ≤ ℓ max ≃ R MT γ K max , determined such as to guarantee continuity and smoothness of the basis function at the muffin boundary.ℓ max is a numerical cut-off parameter often set by the muffin-tin radius R MT γ and the largest reciprocal lattice vector controlling the quality of the basis, K max .The radial functions are represented at a set of N MT radial mesh points.
Consequently, it is natural to choose the computational domain consisting of MT and IR also for the charge density, n(r), and potential, V eff (r), and expand both into plane waves (up to a maximal wave vector length G max ) and radial functions times spherical harmonics (up to an angular quantum number L max ≤ 2ℓ max ), as exemplified here for the densities: In practice, the support of the different regions is mediated by step functions Θ γ and Θ IR .Θ γ are 1 in the respective MT γ sphere of atom γ and 0 everywhere else.The step function Θ IR = 1 − γ Θ γ removes the MT region from its integration.The Fourier representation Θ γ (G) can be found in equation 5.41 of reference [66].As the representation of the electronic structure in the MT spheres now explicitly depends on the atomic positions, several amendments to the previously outlined theory of phononic perturbations become necessary.These will be discussed in the next chapter.
For practical purposes, the symmetry properties of the crystal lattice are used and in the FLEUR code [71,72], the charge density and potential is represented in terms of symmetrized plane waves, so-called star-functions [66], and symmetrized spherical harmonics, so-called lattice harmonics [88].For the sake of readability, we largely omit this additional layer of representation in this work.
In the LAPW basis, the Kohn-Sham equation (7a) turns into a generalized eigenvalue problem with a Hermitian Hamiltonian and overlap matrix of dimension respectively.The setup of the Hamiltonian matrix due to the non-spherical potential V L ′ in each MT sphere expressed as in (12) comprises in the order of O(N A (ℓ max + 1) 4 ) matrix elements of the type ⟨L|V L ′ |L ′′ ⟩, i.e. between the basis functions in spherical representation with the angular and magnetic quantum number L and L ′′ and the non-spherical potential.
In order to exclude matrix elements of irrelevant magnitude we introduce an angular momentum cut-off ℓ max,nsph ≤ ℓ max for the basis functions ℓ and ℓ ′′ contributing to the Hamiltonian setup.The diagonalization of the eigenvalue problem (14) is the runtime determining step of a self-consistent determination of the ground-state charge density, making the runtime of LAPW methods scale ∝ N k × N 3 B , where scaling with respect to the number of basis functions, N B , stands also for the precision scaling of the physical properties as well as the volume scaling O(N 3 A ), as the number of basis functions scale linearly with the number of atoms, N B ∝ N A .

Implementation
The central motivation is to determine the phonon dispersion (1) in the harmonic approximation by means of the force constant matrix (8), which requires the determination of the charge-density response, the wave function response, the response of the external and effective potential, as well as the second order changes in the external potential and the ion-ion energy.In the following outline of the implementation we deal with vectors and matrices in the space of the 3N A dynamical degrees of freedom and of the electronic degrees of freedom determined by the number of basis functions, N B , in which the Kohn-Sham orbitals are expanded.Some quantities are vectors in one space and matrices in the other.In general, we do not distinguish both by different types of vector or matrix symbols for different types of spaces, but depending on the context one space is emphasized over the other by the relevant vector or matrix symbol.

General Concept
The DFPT formalism to phonon properties consists of three parts: (i) a single groundstate DFT calculation, (ii) the setup and convergence of the Sternheimer equation to obtain the first-order density response n (1) upon the displacement of atoms, as well as (iii) the setup and diagonalization of the dynamical matrix.Since the phonon is a wavevector q dependent displacive perturbation, many response or perturbed quantities carry naturally the superscript q (for more details see section 3.2).The dynamical matrix is the central input to any phonon property calculator.It is set-up and diagonalized q-vector by q-vector in a totally sequential but parallelizable fashion.The basic algorithm is sketched in Figure 1.After the ground-state DFT calculation (red box) was carried out, we start the DFPT part.For this we need to set-up some additional quantities (blue box) beyond those we determined during the ground state run.These are calculated once before the start of the outer q-loop (blue frame), as they are not q-dependent.They consist of (i) the gradients of the ground-state density and all constituents of the potential (external, Coulomb, xc), as they are needed to determine the basis set corrections and the corrections of the discontinuities at the surface of the MT sphere, and (ii) of the complete set of eigenvalues and eigenvectors of the unperturbed Hamiltonian, which are needed for the quasi-analytical inversion of the Hamiltonian and overlap matrices in the Sternheimer equation.This yields a set of self-consistency equations (red frame) for each q-point, as well as atom in the unit cell and displacement direction (cyan frame).By treating each of the 3N A perturbations sequentially for each atom β and cartesian displacement coordinate j, the response quantities are just scalar components of vectorial quantities of dimension 3N A .To keep the overview, we have omitted the indices β, j in workflow Figure 1.We construct the (βj)-components of the first-order external potential (V (1)q ext ) as our initial perturbation and establish then the Hamiltonian response (H (1)q ) together with the overlap matrix (S (1)q ).From the Sternheimer equation we thus determine for a given k-vector the response of the wave function expansion coefficients for all electronic eigenstates ν.After the k-point loop is completed, these enter the charge density generator to construct the (βj)-component of the density response (n (1)q ), which is then used to construct the (βj)-component of the effective potential response (V (1)q eff ).This accounts for the self-consistent nature of the problem, since the effective potential response in turn requires a new Hamiltonian and overlap matrix response.We repeat the calculation until the response-density changes less than a given threshold between iterations (red frame).After convergence is reached, a final iteration of the loop is started to construct additional quantities needed to compute the (βj)-row of the dynamic matrix (D βq⊤ j ).Finally, for this we also need the second-order variations of the external potential (V (2)q ext ) and the ion-ion interaction (E (2)q ion−ion ).The individual parts of the workflow are laid out in section 3.4.We deal first with the density response (section 3.4.1),then with the exact form of the Sternheimer equation and the related matrices (section 3.4.2), the generation of the potential responses (section 3.4.3),and finally with the mixing procedure (section 3.4.4).The dynamical matrix setup is found in the succeeding section 3.5.In all calculation steps apart from the mixing, the first and second order quantities are determined in part by quantities obtained from the ground-state calculation, which will be highlighted accordingly.

Definitions
First we define the perturbed quantities: provided the periodic displacement ∆τ of atom α in unit cell R with amplitude Q α and phonon wave vector q offsetting the atom from the equilibrium position R + τ α is expressed as any perturbed or response quantity X (n) of order n, e.g.think of the charge density response n (1) , will be a sum of terms with n different ±q-dependent phases.So, properties originally periodic according to the translation symmetry of the crystalline lattice now carry an additional plane-wave factor, potentially altering the Bloch character.For first and second order quantities, X (1) and X (2) , respectively, of the general quantity X, we subsume the N A atoms (α) and the three Cartesian coordinates (i) into 3N A dimensional vectors, X (1) , and 3N A ×3N A dimensional matrices (underlined quantities), X (2) , respectively: where superscript ⊤ stands for the transpose operation, and use of the relation Q(−q) = Q * (q) was made.The following notation was introduced: Response quantities, either represented as vector or matrices or with indices, represent the direct derivatives with respect to the phononic perturbation, as e.g. the Hesse matrix to second order, while a absence of indices corresponds to the full scalar perturbation of that order, i.e. for the same example the contraction of the Hesse matrix with the displacement vectors.X (•) (q) refers the Fourier transform of X (•) over the unit cells.Working with these projections enables us to suppress the dimensional character of the involved quantities.It can easily be shown, by requiring integrals over the unit cell to be non-vanishing, that (i) to second order only combinations of +− or −+ contribute [61] and (ii) those terms are the Hermitian conjugates of each other.To that end, it is sufficient to solely calculate the + part to first order in the eigenstates, charge density, and potential perturbations as well as to only calculate the DM for the combination +−.By denoting the atomic and directional indices for the columns by αi and the rows by βj, we need to calculate quantities such as:

Challenges
The LAPW basis function adds a position-dependent basis set in the MT region to the basis consisting of plane waves in the IR, which leads to non-trivial additional complexities.When the wave function is varied due to a phonon perturbation, owing to the displacement sensitivity of the basis, the first-order change of the wave function is not expressed only by the first-order change of the expansion coefficients to a momentum increased by the phonon wave vector q, z k+G+q,ν , Ψ but also by a term that is the derivative ϕ (1)R k+G (r) of the basis function with respect to the atomic displacement in a unit cell R.This latter term (and all related quantities) lies typically outside the Hilbert space (HS) spanned by the original LAPW basis functions and cannot be included efficiently by increasing the number of functions.We distinguish between two out-of-HS contributions: (i) the Pulay terms [89], akin to a set of corrections in atomic-force calculations [69], and (ii) the variation of the effective potential, affecting the radial solutions u ℓ in the MT part of the LAPW basis (basis response) [90].Both terms together are known as the incomplete basis set correction (IBC).While the Pulay terms are indispensable for a successful DFPT calculation in the FLAPW method, the basis response is (according to literature for force calculations) assumed to be small [69].This results in the so-called frozen-augmentation approximation, where this contribution is neglected and that we likewise adopt.The remaining terms stem from the differentiation of the matching coefficients, yielding an imaginary prefactor, and the direct differentiation of the position dependence, which is expressed as a gradient, ∇, with respect to the space coordinate r.Written with the full set of response indices, we find The latter gradient operator poses a numerical challenge, as we have to calculate the corresponding radial derivative of u ℓ on a finite radial grid that extends to the atomic nucleus, where this can become inaccurate.Another challenge lies in the numerical discontinuity of the LAPW basis functions at the MT boundary (see also discussion in Appendix C).Although the LAPW basis is in principle differentiable throughout the unit cell when the real-space basis function in the MT region is expanded into an unlimited number of angular momentum coefficients ℓ, reasonable cutoffs lead to slight discontinuities in zeroth and first order or even severe discontinuities in higher orders.Especially in the case of phonons, these higher orders become relevant (e.g. when applying a dyadic product of gradients to the basis function).To mitigate this and the dependence of integrals over quantities X involving wave functions on the positions of the moving nuclei, additional surface integrals need to be considered [79]: Terms of the form [X] SF are to be understood as the difference of the function taken in the MT domain (∂β + ) and in the interstitial domain (∂β − ).The former can be used to cancel one problem against another.In practice, each response of a r-dependent quantity has a certain resemblance to the gradient of the same quantity in the displaced MT-sphere.For the response of the basis function, this similarity is analytically explicit, while for quantities like the potential responses, it is implicit.Overall, it is beneficial to sum up response quantities with the corresponding gradient.The gradients can be readily obtained at various points in the calculation, due to the fact that any surface integral of a function over a closed surface can be rewritten as an integral over the enclosed volume of the function's gradient.We make heavy use of this and hence eliminate the necessity to deal with gradient terms as often as possible, which leads to the regrouping of terms and the existence of surface integrals of IR quantities, that are not paired with their corresponding MT representation anymore.As a last aside, we need a numerically stable representation of the various Coulomb potentials and their derivatives, as well as the Coulomb energy between the nuclei, despite the 1/r -singularities at the nuclei, in particular when a gradient or a dyadic product of gradients is involved.Here the method of Weinert [57] for solving Poisson's equation without shape approximation for an arbitrary periodic charge distribution can be used and generalized.The Coulomb terms are then perfectly continuous by construction (i.e.chapter 3.4.3).It turns out that the standard FLEUR integration scheme (6-point Simpson integrator) for radial integrals in the MT sphere has trouble handling the density gradient and related quantities, due to wild oscillations at the core that propagate outwards, and it was thus replaced by a 4-point spline integrator in the DFPT part of the code.

Sternheimer Equation
We start with elaborating on the first-order density variation, then discuss the setup of the Sternheimer equation, shortly introduce how we calculate the linear potential variations and close with a brief overview of our strategy to achieve self-consistency in the Sternheimer equation.

Linear Density Response
The representation of the first-order density response n (1) around the unperturbed density n in the FLAPW method comprises various terms.Ultimately, it comes back to the wave function response.The wave function part related to the temperature dependent occupation function ( fkν ) and expansion coefficient (z k+G,ν ) responses are very similar to the plane-wave part in the pseudopotential method [37], but with the basis-function response we instead find (using the simplification of Appendix A) where a possible spin index is omitted again.Explicit terms are induced as (i) the "direct" response of the expansion coefficients, (ii) the consequence of the basis set variation that creates a term with an imaginary prefactor in the displaced MT sphere and a gradient term (that is in practice grouped with its complex conjugate into a ∇ j nterm), as well as (iii) a term dependent on the perturbed occupation function of the electronic state (for materials like metals where the occupation is fractional).Thus, we need to precalculate the gradient of the unperturbed density (as well as gradients of the potential for steps later in the calculation).The corresponding formulae and details on the numerical accuracy can be found in reference [80].It is noteworthy to realize that the density response field is a functional of the unperturbed density, its gradient, the first order change of the expansion coefficient, which is a functional of the first order change of the potential, which depends on the first order change of density, n (1) [n, ∇n, z (1) [V (1) [n, n (1) ]]].Evaluating (24) for the part independent of f (1)βj kν in the IR results in whereas in the MT sphere of atom γ we implement where the brackets [ • ] in [∇n] MT γ ,L (r γ ) denote the expansion of the density gradient into lattice harmonics [88].The coefficients d, containing the linear response of the matching coefficients, expansion coefficients, and occupation numbers, are defined in Appendix B. Within these interstitial and muffin-tin representations, the central quantity not given by a preceding DFT calculation is the first-order variation of the wave function coefficients z (1)βj k+G+q,ν .They are determined by a self-consistent solution of the Sternheimer equation.
Core electrons In the charge density response quantity, n (1) , we explicitly consider the valence states only, i.e. we apply the frozen-core approximation, stating the core electrons not to be perturbed by a shift of the atomic positions, while in the ground state density and the gradient of the ground state density, ∇n, the full density enters, i.e. including the core electrons.In FLEUR, there is furthermore an option to explicitly consider electrons leaking out of the MT spheres -they originate from and permeate other muffin-tins -by applying so-called core-tail corrections [66,91], which would warrant additional terms to first order.We postpone the implementation of their perturbation [73], because core-tails can be suppressed by using local orbitals [87].The core electrons do, however, contribute to the gradient of the all-electron charge density in the same equation.

Setup of the Sternheimer Equation
As mentioned before, the Sternheimer equation takes a more lengthy form in LAPW methods as compared to pure plane-wave formulations.Inserting ( 21) into (10) and explicitly accounting for the basis variation and the surface terms yields in the space of LAPW basis functions The first line constitutes the Hellmann-Feynman contribution, the second line contains Pulay terms, and the third line contains surface terms.In this representation, the Pulay terms are of significant value.They consist of a prefactor part and a part containing the gradient of the basis function.The latter are not so numerically well-behaved, especially in the core region, that their numerical integration over the muffin-tin sphere guarantees sufficient accuracy for reliable phonon properties.
To solve this equation, we exploit that the left side is a matrix-vector product, and we need to invert the matrix.Instead of working in the space spanned by the LAPW basis functions, we switch to a representation of the Sternheimer equation where the space is spanned by the Kohn-Sham wave functions by multiplying (27) from the left with z * k+G ′ +q,ν and contracting over G ′ .This procedure avoids a costly inversion of a matrix that is nearly singular at certain eigenvalues.Given the definitions of both the band representation of the perturbed expansion coefficients and the prefactor part of the perturbed wave function we rewrite (27) into with the kinetic energy operator T (on which some notes are found in Appendix C).Since V eff depends on n (1) and, therefore, on z (1) , this Sternheimer equation must be solved self-consistently according to the scheme in Figure 1.Comparing this form with the initial Sternheimer equation ( 27), which subdivides into the Hellmann-Feynman (first line), the Pulay (2nd line) and the surface terms (last line), we now group the contributions differently (as discussed in section 3.3).We highlight (i) the complete representation of the Sternheimer equation in the Kohn-Sham wave function spanned Hilbert space (and consequently the contraction of the G-vectors), (ii) the summation of the first-order effective potential and the gradient of the unperturbed potential in the muffin-tin matrix element of atom γ, that avoids the integration over large terms around the center of the MT spheres, (iii) the overall avoidance of gradients of wave functions and thus contributions outside the established LAPW Hilbert space by cancelling them with the MT surface terms, and (iv) the grouping of IR terms into a combined perturbation of the interstitial potential and step function in Fourier representation.Considering (ii), we introduce a shorthand notation for combinations of perturbations and gradients in the displaced MT.We write to streamline further equations, as such combinations reappear frequently in the dynamical matrix setup.Details on the general evaluation of IR or MT matrix elements are pointed out in reference [80].It becomes obvious that knowledge about the unperturbed system at shifted Bloch vectors k + q is required.We choose to calculate all that information once in the beginning, before the Sternheimer loop of a particular q-point.As a test of the implementation, one can show that the analytical solution of the Sternheimer equation for q = 0 and one atom is [80] z There are two more things to consider.Firstly, (30) only holds for non-vanishing energy differences δ qν ′ ,kν := |ϵ k+q,ν ′ − ϵ kν |.We group terms with and without a prefactor ϵ kν together and identify them as a perturbed overlap matrix S (1) and Hamiltonian H (1) , respectively, also referred to as overlap matrix response and Hamiltonian response in the space of the N B × N B Kohn-Sham states and at the same time 3N A dimensional vectors in the space spanning the dynamical matrix, to rewrite (30) as If the energy difference δ qν ′ ,kν is close to zero (a threshold of 10 −12 htr was used for the calculations in this work), a special treatment is required in order to avoid an explicit division by very small numbers, and the following reformulated expression (see derivation in Appendix D) is applied.Secondly, arithmetic and sum reformulations yield an individual procedure for the case in which the energy difference is finite and both ν ′ , ν represent occupied states, meaning their occupation-number prefactor f is larger than a certain threshold (set by default to 10 −8 in FLEUR).Provided this condition, we can derive where F (ϵ k+q,ν ′ ) is the Fermi smearing for the respective eigenenergy (see Appendix H).Its smoothness can be controlled by the Fermi smearing parameter k B T .Using this modification can improve the stability of the self-consistency calculation.

Potential Responses V
eff and ∇V eff , the first order response and the gradient of the effective potential V eff both enter the Sternheimer equation (30) as well as the set-up of the dynamical matrix.For the latter, we also need the response and gradient of the external V ext and of the Coulomb potential V C , due to various correction terms that occur as a consequence of the LAPW basis (see section 3.5).Since we need each term up to first order and also the corresponding real-space gradients.We briefly discuss the calculations of these terms in the following.
Hartree and external potential -The Hartree potential response, V H , is basically the Hartree potential of the response charge density n (1) , and the Hartree potential gradient, ∇V C , is basically the Hartree potential of the gradient of the charge density ∇n, plus additional surface integrals introduced in ( 23).These surface corrections apply for the displaced atom in the potential response calculation and for all atoms in the gradient case and correct possible discontinuities at the muffin-tin boundary.Therefore, we use the Weinert algorithm [57] for solving the Poisson equation to obtain both V (1) C and ∇ j V C from the charge density response to first order and the gradient of the unperturbed charge density, respectively.The procedure is similar to that of the ground-state calculation, where the radially dependent MT densities are replaced by a smooth Fourier transformable pseudo-density valid in the whole unit cell, so that the interstitial potential can be directly expressed by the pseudo-density components, while the MT potential is obtained by solving a boundary value problem inside the MT sphere.However, instead of using the density, we use the first-order response density or the gradient of the density.We employ the Weinert algorithm not only for the response, but also for the gradient of the Coulomb potential, since continuity at the muffin-tin boundary (being important for well-behaving numerics) is then ensured by construction, which is not the case if we straight-forwardly differentiate the ground-state Coulomb potential across the sphere boundary.
Basically, we employ (with small modifications) the same Coulomb solver routines for the potential response as for the unperturbed potential, but there are some points to consider.Firstly, additional care has to be taken with respect to the radial integration of the density response or density gradient in the Coulomb solver, as they can be less smooth than the typically used ground-state density.A second point is the emergence of surface terms in (37) and (38).We express them here as a correction to the basic multipole moments q lm of degree (l, m) as described in [57].Keeping to Weinert's original notation (Equation (11) of his paper), the effective multipole moment of the response charge density inside the ith MT sphere can be written as qi lm = q i lm − q iI lm + q i,SF lm − q iI,SF lm .
The first two terms describe the multipole moments of the true response charge density in the sphere of atom i subtracted by the multipole moments of the plane-wave response charge in the ith atom.The second two terms, denoted by the superscript SF, correspond to corrections from the surface integral.Lastly, the infinitesimal displacement of the Coulomb singularity of the atoms contributes in first order response term by an (ℓ = 1)character instead of being spherical with (ℓ = 0, m = 0).Aside from this, the Weinert procedure is used as in his seminal paper.The specifics of the modified terms can be found in Appendix E, while the full derivation of the adapted method was elaborated on in reference [79].Like in the original method, where the combination of the Hartree and external potential mitigates the ∼ 1/r singularity, here the ∼ 1/r 2 contributions with large absolute values compensate each other and lead to a better controllable numerical behaviour.
Exchange correlation potential -In order to calculate the first-order variation and the gradient of the xc potential both the first-order variation of the charge density and the gradient of the unperturbed density are multiplied with the exchange-correlation kernel that is the functional derivative of xc-potential with respect to the charge density evaluated at the DFT ground state density n (0) (r) of the unperturbed system.Algorithmically, all operations are carried out in real space after the respective IR and MT coefficients of the density response and density gradients in the coefficient space (see definition (13) for clarity) have been transformed to real space by Fourier transformation and the evaluation of the lattice harmonics [88] on a spherical grid, respectively.The results of the multiplication are then transformed back to coefficient space.
For the sake of algorithmic locality we recalculate the xc-kernal at each Sternheimer iteration again, although it depends only on the ground-state density and does not change with iteration.For the sake of convenience, we employ the libxc library of functionals [92] making all necessary quantities readily available when provided with the real-space density.Currently, we are limited to LDA functionals.In the future we plan also an extension to GGA functionals, for which the evaluation of the xc kernel is significantly more involved.

Achieving Self-Consistency
The self-consistent solution for the charge-density response field by means of the Sternheimer equation bares a lot similarities to the selfconsistent solution of the charge density by means of the Kohn-Sham equation in a conventional DFT calculation.In both cases we deal with a nonlinear problem that is solved iteratively.The output response density n (1) (m+1) , here and below written as scalar quantity as each atom and displacement coordinate is converged independently, obtained after completing iteration step m is a functional of the first-order change of the potential V (1) (m+1) is calculated.Therefore, we adopted the existing charge-density mixing technology [93] and mixed the charge-density response (m) , α mix (42) according to the Broyden-like scheme of Anderson [94] with the same mixing parameter α mix as in the ground-state calculation.δn (m) is the Anderson-preconditioned residual response density F n in with the preconditioner synthesized from the history of all charge density responses letting n the input and output response charge density at some iteration m, respectively.To cope with the fact that in the MT region the density response is complex valued while the original charge density is real valued, we mix the real and imaginary part of the response density independently by mapping both onto the mixing scheme for a magnetic system relating the real and imaginary part of the density response to spin-up and -down densities.
In reality, the response density depends also on the ground-state charge densities.Throughout the self-consistency cycles, however, all ground-state properties remain unchanged and thus this charge-density field is taken off from the mixing procedure, as it provides a constant static offset that might contribute to an instability of the procedure.After the mixing is completed these terms are added again.As a measure of convergence we use the L 2 -norm induced metric distance n (1)βj+ out , n and require it to be smaller than a preset threshold of ϵ scf .All tests done to this point indicate very stable convergence behaviour for any material which converged properly in its ground-state calculation.But it should be noted that for the first few iterations the distance can start from very large values, especially when dealing with small q-vectors.Some additional details on the mixing are found in Appendix F.

Dynamical Matrix
We recall, that according to (2) and ( 3), the DM is related to the second derivative of the Born-Oppenheimer energy surface.From (8) we have seen, that this second derivative is related to the Coulomb interaction between the charge density response n (1)  and the perturbed external potential V ext generated by the nuclear charge, the Coulomb interaction interaction between the ground-state density n and the second order external potential V (2) ext , and the second derivative of the repulsive Coulomb energy generated by the nuclear charges.The Hellmann-Feynman force constant ( 8) is an important contribution to the DM, but it is incomplete for many electronic structure methods, in particular for the LAPW basis set.In the following we derive the DM step by step starting with the first derivative of the Born-Oppenheimer energy surface.Of course, by this we find the Hellmann-Feynman terms again, but also the Pulay terms, and the terms due to the discontinuity at the MT-sphere boundary.
The ground-state energy (6) per unit cell of volume Ω of the unperturbed system can equivalently be expressed in terms of the Kohn-Sham eigenvalues as where the first term together with the third one corresponds to the kinetic energy T 0 , the fifth term is the Hartree energy E H , the sixth term the exchange-correlation energy E xc and the last term the Coulomb energy E ion−ion between nuclei of atoms α and β with atomic numbers Z α , Z β of the energy functional (6).We introduced a term dependent on the temperature T and electronic entropy S as proposed by Weinert and Davenport [95] to deal consistently with the temperature dependent Fermi-Dirac distribution of the occupation of electron states in case of metals.
From this we derive an optimized representation of the first-order total energy variation.All contributions related to the first-order occupation numbers cancel between the sum of the single particle energies ϵ kν and entropy terms, and we find for a displacement of atom α along coordinate i E tot is, aside from the explicit q-dependence in the first-order quantities and the fact that it is not evaluated for a finite displacement with amplitude Q (see (15)), reminiscent of the LAPW force expression introduced by Yu and Krakauer [69] with the discontinuity extension of Klüppelberg et al. [73], and thus the i-th force component acting on atom α, F α i , is related to E (1)αi tot (q = 0) as . The first line of (44a) corresponds to the well-known Hellmann-Feynman force and, like in the implemented form of the Sternheimer equation ( 30), the Pulay and the MT surface-term contributions are smartly rearranged to discard gradients applied to e.g.wave functions (by reformulation of the MT surface integrals into volume integrals of gradients).We arrive at (i) state-dependent correction terms C (1)αi− kν , which are a sum of typical Pulay-type and MT surface terms evaluated in the MT-sphere and IR, respectively, (ii) the potential energy of the ground-state charge density in the field of the gradient of the Coulomb potential ∇V C in the displaced muffin-tin sphere, as well as (iii) the electrostatic energy of the charge density in the field of the Coulomb potential and the exchange correlation energy E xc both evaluated in the IR with the perturbed step function.
Based on the same reformulation ideas, we obtain the following collection of terms for the second order change of the total energy per unit-cell volume with respect of the displacement of atom α into direction i and atom β into direction j: (r) dr This lengthy expression is the complete FLAPW analogue of the Hellmann-Feynman expression of the Hessian matrix (8).The first four lines constitute the Hellmann-Feynman part, where parts of the integral terms stem from rearrangements by partial integration to avoid second order dyadic gradient terms (∇ j ∇ i ) and gradient terms of the perturbed quantities at the expense of additional surface integrals.This was done on account of the observation that such terms (resulting here from the double direct differentiation of the external potential) are very demanding for the radial integration and are a major source of numerical inaccuracies.The same rationale holds for the various integral terms in the bottom four lines.The (kν)-dependent terms contain only the part of the IBC that is directly basis dependent and mixes Pulay and surface contributions.The composition of the coefficients C (2)βj+αi− kν can be found in Appendix G.
We derive the second-order variation of the ion-ion interaction, E ion−ion , following a scheme for the ground state energy already published by Weinert [59], bearing similarities to the perturbed electronic potentials.Ultimately, we use The parameter N appearing in the pseudodensity, n ps , is chosen for its optimal convergence and that of the IR potential according to the given choice of ℓ max and G max R MT,max [57].We choose this parameter in the same way as for the calculation of the ground state.The expression for n ps results from the evaluation of equation ( 28) in [57] for ℓ = 2 with multipole coefficients representing the second order atomic displacements already calculated and expressed as factors containing the reciprocal wave vector components (G + q) i/j .The power N + 3 in the denominator is two orders higher than in the reference, which is compensated again by the factors G + q occurring in the nominator resulting from the second-order differentiation of the energy.The numerical quality of this formalism is in good agreement with the results obtained from the ABINIT [44,96,97] code, where the algorithm is based on an Ewald approach [80].
We then set up the DM by symmetrizing the energy perturbation to ensure its hermiticity and dividing each element by a factor of M α M β .Then, we calculate all eigenvalues {λ µ } µ=1,...,3N A and eigenvectors of the Hermitian matrix using a standard eigenvalue solver [98].The actual phonon frequencies {ω µ } µ=1,...,3N are the square roots of these eigenvalues.As described in the literature, if an eigenvalue is positive, we take the resulting square root with a positive sign, and if it is negative and thus the square root would yield an imaginary frequency, we represent it as real with a negative value in our calculated phonon dispersions.At negative frequencies, the phonon spectrum thus indicates instabilities in the crystal lattice.A deeper insight into the technical nuances of the implementation (such as integral evaluations, pseudodensity coefficients, and gradient calulations) is provided as an integral part of references [79] and [80].

Scaling behaviour
The DFPT algorithm comes on top of a ground-state calculation whose computational effort was briefly discussed in section 2.3 and for which a detailed discussion can be found in [99].The runtime determining step of the DFPT algorithm is the iterative solution of the Sternheimer equation (30) for each wave vector q, for all three Cartesian coordinates of the displacement perturbation, all N A atoms in the unit cell, and all N k k-points in the BZ.In practice, this is done by a series of matrix multiplications and thus the computational effort is bounded by the largest among them.This is already the first one, where we multiply the perturbed Hamiltonian and overlap matrices (N B × N B , where N B is the number of basis functions as determined by K max ) with the matrix of unperturbed expansion coefficients in the occupied subspace (N B ×N o ), with the number of occcupied states N o .The order of operations for this multiplication is O(N o N 2 B ).The other matrix multiplications are of the same order, as the dimension of the occupied subspace gets passed on with each product, and there is no proper matrix inversion necessary for the initial Hamiltonian and overlap, as we use the spectral representation for a quasi-analytic inversion.This is of the order O(N o N B ). Summarizing, for each wave vector q, the runtime of the DFPT algorithm scales as . Since the number of occupied states as well as the number of basis functions scale with the number of atoms, the DFPT has a volume scaling of O(N 4 A ) and the precision scaling is of O(N 2 B ) in the number of basis functions.Although in the DFPT approach, the volume scaling is worse than for the conventional DFT self-consistency cycle (∝ N k N 3 B ), in the FLAPW method the number of occupied states are only a fraction of all N B , e.g. in fcc Ne we find 4 occupied states for 162 to 177 states overall (depending on the k-point).This is at most 2.5%.In general we expect a maximum occupancy in the order of 5-10%.Currently we use all available unoccupied states (N B − N o ≃ N B ) in calculating the response matrix.Thus, N o produces a prefactor that is a fraction of N B and an iteration of the Sternheimer loop is faster than that of a conventional DFT calculation with no symmetry.The memory requirement, as opposed to the groundstate calculation, is more than tripled.This is due to the necessity of not only keeping the occupied unperturbed eigenvalues ϵ kν and eigenvectors z kν in storage, but also the full set of unperturbed ϵ k+q,ν ′ and z k+q,ν ′ , as well as the occupied perturbed quantities ϵ (1)βj k+q,ν and z (1)βj k+q,ν .The q-dependent quantities, however, can be deleted once a specific q-point calculation is finished.

Results and Discussion
In this section we validate our DFPT framework with respect to the quality of Goldstone modes and phonon dispersion relations against the FD approach.We choose a set of six distinct elemental materials, none of which share the same attributes.We cover a simple alkali metal (Na), several magnetic (Fe, Ni) transition and noble metals (Cu) exhibiting different crystal structures, a semiconductor (Si) and an insulating noble gas crystal (Ne).The alkali metal is distinguished by its rather simple Fermi surface with a low number of electrons, and the noble gas crystal by its low-energy phonon modes.Our strategy is the following: We set up the unperturbed unit cell of the material under study and optimize its volume (fit the total energy curve as a function of different lattice constants to the Birch-Murnaghan equation of states [100]) and internal degrees of freedom of the atom positions if necessary.We use the resulting structure as input for comparative FD calculations with phonopy and DFPT in FLEUR.We will first give a short summary on how the FD calculations are conducted.To begin a FD calculation, we provide the unit cell optimized by FLEUR as input to phonopy along with a 3×3 matrix of integers, M S , that extends the original Bravais lattice, A u , to a supercell with lattice vectors A S , by the matrix multiplication: The supercell is subsequently filled with copies of the original set of atoms at appropriate positions.To ensure the best possible comparability between our benchmarks, we set a list of computational parameters identically for all materials considered (Table 1) and we work with the same k-point densities across all different Brillouin zones in use.We have chosen a k-point set of 16 × 16 × 16 for the ground-state calculations performed in the primitive unit cell and adjusted the k-point set for the supercells accordingly.As default size of the supercell we chose 2 × 2 × 2 times the primitive unit cell with a k-point set reduced to 8 × 8 × 8.For the 4 × 4 × 4 supercell we work with a 4 × 4 × 4-k-point set.Aside from parameters previously mentioned, there is the force convergence criterion ϵ force (similar to ϵ scf , but for the difference between the forces in two iterations), and the force level (0 means there will be no corrections as described in reference [73], as we do not expect significant drift forces emerging for the selected materials).The same parameters are used for all calculations, i.e. the groundstate calculation, the supercell ground-state and force calculation for the FD supercell, and the DFPT run for all systems discussed here.
Along with a perfect supercell, phonopy analyzes the symmetry of the system and gives a list of supercell inputs with displacements that include all information needed to construct the force constant matrix and thereby the dynamical matrix.This list of inputs goes back to the FLEUR code, which calculates the corresponding forces upon each suggested displacement.The force calculations for the different displacements are fully independent of each other, so the process can be run in parallel.Based on the set of force and displacement information, phonopy delivers the force constant matrix and the final output is a continuous phonon dispersion relation based on a Fourier transform of this matrix onto reciprocal space.We acknowledge that the FD method contains harmonic and anharmonic contributions to the phonon-dispersion.The anharmonic contribution depends on the magnitude of the displacement amplitude.To compare our results with the DFPT, which contains strictly only the harmonic terms, we have carefully monitored the role of the displacement amplitude.Finally, we use a displacement amplitude of 0.02 a 0 for each of the structures.

Computational Details: Density Functional Perturbation Theory Phonons
From the optimized FLEUR input cell, the DFPT calculation is started with the same computational parameter set as for the FD benchmark.It is important to note that the cutoff K max , which limits the number of reciprocal lattice vectors G for every k-point according to |k + G| ≤ K max is also applied to the q-shifted k-points k + G + q, just as the cutoff G max is applied to the density and potential responses.Practical experience has shown that we numerically obtain the best results when the q-points for which the phonon properties are calculated are part of the k-point mesh.Thus the choice of the selected q-points impacts also the choice of the equidistant k-point mesh.We would also like to point out that the differentiation of a function expanded into an angular momentum representation with angular momentum index ℓ, also generates contributions in the angular momentum components of index ℓ±1.To fully capture these components, we increase the maximum angular momentum of the LAPW basis set from ℓ max in a DFT calculation, which is typically an even number, to ℓ max + 1 in the DFPT calculation, which explains the odd values of ℓ max in Table 1.Analogously we proceed for the response charge density and potential.Also here we increase the angular expansion to L max + 1.By cubic point group symmetry, these angular momentum components are not occupied for ground state calculations and thus these quantities are not altered by increasing the cutoff by 1.
We first run a standard ground-state DFT calculation (red box in Figure 1) and, after the density is converged, modify the inp.xml file so that all states are taken into account in the eigenvalue determination (numbands="all"), and add a path with all q-points we want to evaluate in the juPhon tag.The calculations were performed for the FLEUR version that can be found on the repository under the Git tag juBranch before DFPT merge.A comprehensive description of the full workflow can additionally be found under the tag Phonon README for paper.Starting the FLEUR calculation with this modified input will calculate dynamical matrices for each q-point provided.

Quality of the Goldstone Modes
For any crystal with N A atoms in the unit cell, the phonon spectrum will have 3N A distinct branches, three acoustic and 3(N A − 1) optical ones, some of which might be degenerate depending on the crystal symmetry.Especially near the Γ-point, q = 0, the acoustic branches are related to the speed of sound in a material by their slope.At the Γ-point, i.e. at the infinite-wavelength limit, the phonon reduces to a rigid translation of the solid, which does not cost any energy and the lowest three frequencies are required to be exactly zero summarized by the acoustic sum rule.In a FD calculation, this corresponds to a vanishing net force summed over all atoms (drift force) [79].In DFPT, Table 2. Overview of the lattice constants, a, in a 0 , MT radii, R MT , in a 0 , and acoustic Γ-point modes, ω, in 1/cm for the selected materials.No frequency exceeds an absolute value of 1.0 in units of 1/cm (or 0.12398 meV -0.029979 THz, respectively) providing an upper bound for the error of the Goldstone mode ω(0) = 0.  with the analytical solution of the Sternheimer equation for monoatomic materials (33), one can show [80] that the dynamical matrix itself must vanish for q = 0, hence making the acoustic phonons gapless Goldstone modes [101].This is not the case for polyatomic solids, where the acoustic branches have finite value and the matrix is not 0 in every element.With respect to the numerical approach taken here, which results to finite accuracy in the evaluation of all equations, this zero condition required by physics is usually not perfectly realized in practice, and in the development of many phonon codes one has chosen to explicitly enforce it by subtracting either the drift force for an FD approach or a diagonal matrix with the three lowest eigenvalues for the DFPT implementations.Thus, evaluating the eigenvalue spectrum for the acoustic modes at Γ-point and in particular their deviation from zero, is a numerical check of the quality of the Goldstone modes and constitutes a very good test for the accuracy of our calculations and whether such corrections are warranted.Table 2 summarizes the Goldstone mode quantities for each of our test systems as well as the material specific parameters.Each of the supercell calculations was carried out with a 2 × 2 × 2 supercell.For each of the systems phonopy suggests excactly one displacement pattern.

Na
The frequencies are overall small, though in general a bit larger for the DFPT case.It can also be seen, that in the FD case the modes are closest to zero for the simplest materials, fcc Cu and fcc Ne.The other materials contain either local orbitals, a spinpolarization or more than one atom in the unit cell.Since the above Goldstone-mode requirement is well met for all systems, we have come to the decision not to correct our spectrum by applying the acoustic sum rule.Furthermore, the deviation from 0 can be seen as a measure of accuracy for the overall frequencies.
We note that the convergence behaviour of the density response is directly linked to that of the DFT ground-state calculation.For fast converging materials, the Γ-point calculation will converge with similar speed.The calculations for other high-symmetry points in the phonon BZ require some more self-consistency iterations and start with higher initial distances (43).The calculation of intermediary q-vectors takes even longer, with the iteration count growing noticeably with decreasing magnitude of q.Overall, the calculations tend to finish in at most 15 iterations.

Comparison of Phonon Dispersion Relations
Here we validate our implementation of the DFPT by comparing the phonon dispersion relations of the materials introduced above with results from FD calculations along high-symmetry lines of the BZ.The results of the FD are shown as red dashed lines and the DFPT data points as blue squares.Since we deal with monoatomic systems (with the exception of Si), we find 3 acoustic modes that are partly degenerate.Overall we find an excellent agreement between the DFPT and FD approach.For some systems we find unsatisfactory convergence at certain q-points.For these cases we investigate the convergence of the dispersion relation with respect to the increase of the supercell size for FD calculations.In the following we first present the alkali metal Na and the noble metal Cu, both having one valence s electron, then we turn to the magnetic transition metals Fe and Ni, and finally we present the covalently bonded semiconductor Si and the van-der-Waals bonded insulating noble gas crystal Ne.Since our main emphasis is on the numerical validation of our results, we do not discuss the physics of the lattice dynamics of the individual systems, but rather try to cover different classes of materials with our examples.
We begin with bcc Na and fcc Cu to test the implementation for simple, nonmagnetic metals.We restrict our DFPT calculations to q-vectors that mediate between the k(')-points of the set sampling the first Brillouin zone, i.e. k + q = k ′ + G, where G is an arbitrary reciprocal lattice vector.q-points unrelated to the k-point grid show a more erratic convergence behaviour and generally lead to unfavourable results.This gives us 24 distinct data points to compare our phonopy curves to: From Figure 2 it can be seen that the overall agreement between FD and DFPT is good, but Cu matches more closely.In this context, it is useful to point out that the frequency scale of Na has twice the resolution of Cu.The Na DFPT data points deviate slightly from the FD curve for the upper (longitudinal) branch along the N-Γ-path, along the Γ-H path the degeneracy between the longitudinal and transversal branch is lifted, which can be recognized by two little blue squares at different frequency for each k-point, i.e. a gap between both branches opens, which is a bit too big, and the high symmetry point P is not reproduced perfectly.A similar picture emerges for two ferromagnetic metals, fcc Ni and bcc Fe (Figure 3).Once again, the agreement for the face-centered cubic material is better than for the body-centered one.Especially the peak left of the H-point and the area right of it are not described well.We speculate that the discrepancy is caused by the FD curve.
To check for both bcc metals Na and Fe, whether the FD curves are not sufficiently converged in some regard and whether these discrepancies originate from the FD or  DFPT part of the data, convergence tests were made.Differences, e.g. between a 8×8×8 and a 16 × 16 × 16 k-point set were marginal.However, rerunning the FD calculations with a bigger 4 × 4 × 4 supercell (this already constitutes a cell with 64 instead of 8 atoms) leads to visibly improved results (Figure 4).It is evident that the match between the curves and data points is neatly improved.We take away that certain materials may require larger supercells, but assume that they will converge slowly with respect to the supercell size and therefore omit further calculations here, as their compute time grows disproportionately.
Finally, we show Si alongside fcc Ne, both FD calculations are carried out in the previous 2 × 2 × 2-supercell, to have an example for a covalently bonded semiconductor and a van-der-Waals bonded insulator with low phonon frequencies.For Si we came  again across the effect of an insufficient supercell size.The mismatch for certain Si branches is clearly visible, so to improve the fit we enlarge the supercell again.This time, we opt to use a more complex one, that reads (M S ) ij = 2 for i ̸ = j and (M S ) ii = −2 for i ∈ {1, 2, 3}.Again, we reduce the corresponding k-point set to 8 × 8 × 8 points.It is equivalent to unfolding the diamond structure, fcc with a 2-atom basis, into a simple cubic supercell with 8 atoms and then duplicating it in each direction.This is computationally much cheaper than calculating a 4 × 4 × 4 supercell, as the number of atoms in the unit cell is just half and the symmetry is reduced less by the single necessary perturbation in one atom.The result is shown in Figure 6.It can easily be seen, that the larger supercell improves the overall match nicely.The results of both methods give a good fit to the experimental data taken from various sources [102][103][104] that we show together with our computational results, making both methods equivalently viable.A good agreement is obtained with reference calculations carried out with the normconserving pseudopotential method [105] in combination with the LDA functional.Comparing the FD and DFPT results for fcc Ne, we find them in very good agreement, especially considering the small overall magnitude of the phonon dispersion.Although we focus in this paper on the internal consistency of the implementation of the DFPT, it is worth mentioning that for Ne the computational results do not agree well with the experimental data.Experiments at low temperatures [106] show a phonon dispersion that, when scaled to 1/cm, is roughly half as high in frequency at its maximum as the dispersion in Figure 5.This is understandable.Since Ne is a van-der-Waals bonded solid, we should have applied a van-der-Waals functional [107].Using the conventional LDA, the Ne bonding becomes too strong and the phonon energy too high.This is consistent with the computationally optimized lattice constant, which is around 3.924 Å, while the experimental data taken at 6 K give a lattice constant of 4.466 ± 0.002 Å.Just for comparison, the theoretical lattice constant of Si (5.401 Å) matches the experimental one (5.431Å [108]) quite well.To include the van-der-Waals functionals into the DFPT algorithm is part of our future plans.

Conclusion and Outlook
We presented an implementation of density-functional perturbation theory (DFPT) in the all-electron full-potential linearized augmented plane-wave (FLAPW) method FLEUR for the calculation of phonons, that is computationally stable and efficient.This complements the DFPT calculations of phonons, which are typically performed using pseudopotential methods with an all-electron approach, and extends an effective application of the DFPT to magnetic systems and systems of localized electrons.The research software is built up modularly and can be extended in the future.We developed and implemented algorithmic concepts to overcome or bypass numerical challenges inherent to the FLAPW concept, which are provided by the Madelung summation, the Coulomb singularity of the potential, the rapidly varying wave functions and charge densities in the vicinity of the nucleus, the calculations of gradients of the all-electron potential, the presence of the core electrons, the incompleteness and the position dependence of the basis-set, the different representations of the basis-set in muffintin-spheres and the interstitial region and their match at the muffin-sphere boundary to a point that the criterion for the Goldstone mode is satisfied to better than 0.125 meV.We highlighted some particularly challenging points and provided nitty-gritty details in how we dealt with them, leading to a collection of stable and accurate results validated by the finite difference (FD) method relying on accurate force calculations with respect to atomic displacements orchestrated by the phonopy software package [41,42].To achieve agreement between the FD and DFPT approaches, we noted the necessity of converging the FD calculations with respect to the supercell size, again confirming the quality of our DFPT results.Considering the calculation of the phonon energy for the same three-dimensional grid of phonon wave vectors, at present, the FD approach shows a lower computational effort and takes less computer time than the DFPT method.This is also due to missing optimizations in the latter case, while for FD, the full symmetry of the atoms and forces can be exploited by the FLEUR code.It should be noted though, that the convergence of DFPT w.r.t. the k-point grid is much better than that of FD w.r.t. the supercell size.An in-depth optimization of the computational parameters with respect of the convergence of both methods in FLEUR will be conducted in the future.
This paper serves as evidence that reliable and efficient phonon calculations with DFPT are possible in the FLAPW method.The computational efficiency can be further advanced by the full implementation of phonon symmetries [109] as well as the implementation of effective parallelization strategies.The extension to polar materials [110], and the implementation of the spin-orbit coupling [111], non-collinear magnetism [112], different exchange correlation functionals such as the generalized gradient approximation (GGA) [113], a van-der-Waals functional [107] or the extension to strongly correlated electrons systems using Hubbard U (DFPT+U) [114] are straightforward and are subject to future work.
Finally, the successful implementation of the DFPT formalism for phonons provides a motivation to translate the algorithm to other types of perturbations, e.g. to gain insight about responses to external electric [115] or magnetic fields [112].
Platform for Research Software Engineering -Preparatory Study (HIRSE PS), the Joint Lab Virtual Materials Design of the Forschungszentrum Jülich funded through the Innovation Fonds of the Federal Ministry of Education and Research (BMBF), the Joint Virtual Laboratory of the Forschungszentrum Jülich and the French Alternative Energies and Atomic Energy Commission -AI, Data Analytics and Scalable Simulation (AIDAS), and the Bavarian Ministry of Economic Affairs, Regional Development and Energy for financial support within the High-Tech Agenda Project "Bausteine für das Quantencomputing auf Basis topologischer Materialien mit experimentellen und theoretischen Ansätzen".We gratefully acknowledge computing time on the supercomputer JURECA [116] at Forschungszentrum Jülich under grant no.jiff13.
Finally, the authors dedicate this work to the memory of Henry Krakauer one of the original developer of the FLAPW method, teacher, advisor and mentor.one atom is present in the unit cell (see Appendix H for further details).Provided the aforementioned conditions, the perturbed occupation-number terms are trivial in the sense that they couple to the unperturbed basis and produce charge density contribution to the response density.In the IR, we consequently add In the MT spheres, the contributions can directly be absorbed into the d-coefficients that were referenced before.So, their full form is where the band-dependent matching coefficients enter as βj k+G+q,ν a k+G+q,γ

Appendix C. Evaluating the Kinetic Energy Operator
There are several different ways of applying the kinetic energy operator T in an APW context.In deriving the Kohn-Sham equations, the variational expression of the kinetic energy of state ν reads with first derivatives acting on the Kohn-Sham orbital Ψ ν of state ν.Conceptually and numerically, it is very convenient to determine the radial basis functions in the MT region as solutions of the Schrödinger equation.Therefore, by applying Green's theorem, one converts the representation of the kinetic energy in terms of the scalar product of two gradient terms into the Schrödinger form with the well-known Laplace operator, as in (7a), The latter term is the integral over the boundary ∂MT of each MT sphere, with the surface element dS pointing outwards of the enclosed domain.Obviously, the surface term is zero if the wave function or its derivative is zero at domain boundary.This is in general not the case if the domain boundary is the surface between the MT and IR region.Applying the expression of the kinetic energy for the MT and IR region, we get the representation of the kinetic energy by the Laplace operator over the entire unit cell plus the difference of the surface terms at the muffin-tin spheres taken once from the domain of the MT and once from the domain of the IR region (for the definition of ∂MT [X]dS see ( 23)) In the limit of increasingly higher angular momentum ℓ max of the radial basis set in the muffin-tin sphere, the difference of the surface intergals converges to zero.In practice, we use finite ℓ max cutoffs and the surface integrals at the boundary discontinuity are finite and not negligible.In the FLEUR code, we go one step further and symmetrize the form (C.4) by applying the Laplace operator to both Ψ ν and Ψ * ν . T The remaining symmetrized average surface contribution T SF,sym is then negligible.It was tested for the DFPT implementation that there is no significant difference for calculations with the mixed form (C.3) as opposed to the symmetrized form (C.5).We opt to use the latter for conformity with the base calculation.

Appendix D. Modifying the Perturbed Expansion Coefficients
One may naively think to ignore expression (34) in case of tiny energy differences δ qν ′ ,kν .However, this is theoretically not correct and can cause numerical trouble at particular q vectors.Instead, in order to derive numerically stable forms of z (1)βj qν ′ ,kν we exploit the T and P symmetry between pairs of occupied states that enter the sum.We inspect the respective part of the first order density response: We take a closer look at the occupied-occupied subspace and introduce a factor 1 = 1 − F (ϵ k+q,ν ′ ) + F (ϵ k+q,ν ′ ) to find This directly corresponds to evaluating (36) for the expansion coefficients.A similar train of thought (without the inserted factor) can be followed for δ qν ′ ,kν ≈ 0 by recognizing that the occupation prefactor will then be the same for both the original and the shifted Bloch vector k.This leads to C , and the Coulomb potential gradient, ∇V C , is largely analogous to the description in Ref. [57], when the density is replaced by the density response or the charge-density gradient, respectively.This appendix serves to outline the differences to the conventional generation of the Coulomb potential in a ground state calculation.
Firstly, there are the surface corrections to the multipole moments q γR,SF ℓm in the MT sphere and in the IR, q γRI,SF ℓm .Concerning the Coulomb potential response, V (1)βj+ C , they result from the displacement of atom β into the direction j by a phonon with wave vector q.For an atom at τ γ in unit cell R, the MT contribution reads q βj+|γR,SF ℓm , (E.1a) where we omit here explicitly the transformation of the density representation in the sphere from lattice harmonics denoted as [ • ], to spherical harmonics, and the interstitial contribution reads q βj+|γRI,SF ℓm .
In both cases we need Gaunt coefficients, G m,m ′ ,m ′′ ℓ,ℓ ′ ,1 = Y * ℓm (Ω)Y ℓ ′ m ′ (Ω)Y 1m ′′ (Ω)dS, and a matrix ζ, that links the natural spherical tensorial coordinates of the magnetic quantum number m ′′ with indices {−1, 0, 1} to the Cartesian ones, The structure factor e iG•τ γ in (E.1b) results in the expression for the pseudo-density and Coulomb potential being evaluated with a reciprocal vector G + q instead of G.The same holds true for the multipole moments of the density response, which is the second main difference to the ground-state procedure.
The nuclear charge contribution Z γ of atom γ to the multipole moments reads q βj+|γR,ext 1m It replaces the standard contribution to q 00 from the spherical Coulomb potential Z γ /r of the positively charged nuclei.Aside from these deviations, the procedure from the seminal paper [57] can be followed.
In the case of the Coulomb potential gradient, ∇ j V C , no additional vector q appears and in comparison to (E.1a), the structure factor vanishes, there is no restriction δ γβ to the displaced MT sphere β, and the expression changes sign.We find then for the MT part of the surface correction: q j|γ,SF ℓm . (E.4) The changes to the IR part and to the nuclear term are analogous.

Appendix F. Peculiarities of the Sternheimer Mixing
Here are two technical notes about the mixing of the density perturbation during the Sternheimer self-consistency loop: Firstly, we decided to mix only the density response without the gradient part of the density that appears in the displaced MT spheres, as we then deal with a more well-behaved quantity, and the gradient does not change between iterations anyway.Secondly, before the mixing starts, two initial cycles of the Sternheimer loop are performed already.The first one with only the external part of the potential perturbation in the Hamiltonian, which can be understood as constructing a "starting perturbation", and the second with the first full effective potential perturbation.This is the first density designated to be mixed.We thereby ensure that all density perturbations coming into the mixing procedure are constructed in the same way with the same kind of potential.

Appendix G. State-Dependent Terms of the Dynamical Matrix
Due to the complexity of the second derivative, a bunch of state-dependent terms appear in the calculation of the DM.With the introduction of matrix-vector products of the (perturbed) expansion coefficients with matrices akin to the Hamiltonian and overlap, the C-coefficients from (45) can be rearranged into a somewhat compact form that looks as follows: C The auxiliary matrices we introduce (with omitted superscripts referring to the perturbations) are modified forms of the unperturbed Hamiltonian and overlap.To first order they are where the main modification is given by a prefactor stemming from the basis variations.The same prefactor, albeit in the other perturbation direction, again modifies the matrices to second order: where, again, for δ qν ′ ,kν ≈ 0 we find F (ϵ k+q,ν ′ ) → 1 and ϵ k+q,ν ′ → ϵ kν .

Appendix H. Calculating the Perturbed Occupation Numbers
The perturbed occupation numbers f , there are two options.A straight-forward way is to iteratively determine the Fermi energy derivative in the same vein as the Fermi energy itself resulting from the ground-state calculation.We instead aim for another analytical scheme that stems from the requirement of a conserved electron count N = kν fkν . (H.4) Once again differentiating both sides using (H.3), and rearranging terms leads to It was taken into account that according to (9), the variation of the left side of (H.4) in terms of an atomic displacement is zero, N (1) = 0.In the case of low smearing, as for insulators, the Fermi energy derivative is taken to be 0.

Figure 1 .
Figure 1.Sketch of workflow for the DFPT calculation leading to phonon properties.The colored frames highlight the loop structure of the calculation.V

( 1 )
eff , which depends on the input response density n (1) (m) , which enters the Sternheimer equation to generate the wave function response from which the chargedensity response n

Figure 2 .
Figure 2. Phonon dispersions for (a) bcc Na and (b) fcc Cu.The red dashed curve shows the FD reference and the blue squares show the DFPT data.The MT-radii and lattice constants are (a) 2.6/7.651 a 0 and (b) 2.2/6.651 a 0 .Both FD calculations are performed with a 2×2×2 supercell.In the case of Na, the LAPW basis is supplemented with local orbitals (LO) to treat the 2s and 2p semicore states as valence states.

Figure 3 .
Figure 3. Phonon dispersions for ferromagnetic (a) bcc Fe and (b) fcc Ni.Colors and symbols like in Figure 2. The MT-radii and lattice constants are (a) 2.2/5.209 a 0 and (b) 2.2/6.466 a 0 .Both FD calculations are performed with a 2 × 2 × 2 supercell.In the case of Fe, the LAPW basis is supplemented with local orbitals (LO) to treat the 3s and 3p semicore states as valence states.

Figure 4 .
Figure 4. Improved phonon dispersions for (a) bcc Na and (b) bcc Fe.Colors and symbols like in Figure 2.Both FD calculations are performed with a 4×4×4 supercell.The other parameters are left unchanged.

Figure 6 .
Figure 6.Improved phonon dispersion for Si with the modified supercell.The red curve shows the FD reference and the blue squares show the DFPT data.Additionally, we show three sets of experimental data.The black diamonds, brown upward triangles, and yellow downward triangles belong to references [102], [103], and [104] respectively.

. 2 F
An additional derivation akin to Appendix D holds true for the usage of the coefficients in the second line of (G.1b), where the exact same coefficients from the Sternheimer loop can be used, albeit resulting in one last additional termC (ϵ k+q,ν ′ ) H (1)βj+ qν ′ ,kν − ϵ k+q,ν ′ S

F
(1)βj kν are analytically derived from their original definition fkν= f k F (x), x = ϵ kν − E F k B T , (H.1)where the smearing function F (x) is taken as the Fermi-Dirac-distribution temperature k B T in units of the Boltzmann constant k B .By taking the derivative of (H.1) and doing some arithmetic, one can find f(1)βj kν = − fkν F (−x) ϵ

Table 1 .
Overview of the shared calculational parameters for the selected materials.Parameters not contained in the table are kept at the FLEUR default or are explicitly mentioned when presenting the respective calculation results.