Lattice Dynamics Calculations based on Density-functional Perturbation Theory in Real Space

A real-space formalism for density-functional perturbation theory (DFPT) is derived and applied for the computation of harmonic vibrational properties in molecules and solids. The practical implementation using numeric atom-centered orbitals as basis functions is demonstrated exemplarily for the all-electron Fritz Haber Institute ab initio molecular simulations (FHI-aims) package. The convergence of the calculations with respect to numerical parameters is carefully investigated and a systematic comparison with finite-difference approaches is performed both for finite (molecules) and extended (periodic) systems. Finally, the scaling tests and scalability tests on massively parallel computer systems demonstrate the computational efficiency.


Introduction
Density-functional theory (DFT) [1,2] is to date the most widely applied method to compute the ground-state electronic structure and total energy for polyatomic systems in chemistry, physics, and material science. Via the Hellmann-Feynman [3,4] theorem the DFT ground state density also provides access to the first derivatives of the total energy, i.e., the forces acting on the nuclei and the stresses acting on the lattice degrees of freedom. The forces and stress in turn can be used to determine equilibrium geometries with optimization algorithms [5], to traverse thermodynamic phase space with ab initio molecular dynamics [6], and even to search for transition states of chemical reactions or structural transitions [7]. Second and higher order derivatives, however, cannot be calculated on the basis of the ground state density alone, but also require knowledge of its response to the corresponding perturbation: The 2n + 1 theorem [8] proves that the nth order derivative of the density/wavefunction is required to determine the 2n + 1th derivative of the total energy. For example, for the calculation of vibrational frequencies and phonon band-structures (second order derivative) the response of the electronic structure to a nuclear displacement (first order derivative) is needed. These derivatives can be calculated in the framework of density-functional perturbation Sun and Bartlett [51] have analytically generalized the formalism to account for non-commensurate perturbations (corresponding to non-Γ periodicity in reciprocal-space), but no practical implementation has been reported.
In the aforementioned reciprocal-space implementations, each perturbation characterized by its reciprocal-space vector q requires an individual DFPT calculation. Accordingly, this formalism can become computationally expensive quite rapidly, whenever the response to the perturbations is required to be known on a very tight q grid. To overcome this computational bottleneck, various interpolation techniques have been proposed in literature: For instance, Giustino et al. [52] suggested to Fourier-transform the reciprocal-space electron-phonon coupling elements to realspace. The spatial localization of the perturbation in real-space (see Fig. 1) allows an accurate interpolation by using Wannier functions as a compact, intermediate representation. In turn, this then enables a back-transformation onto a dense q grid in reciprocalspace.
To our knowledge, however, no real-space DFPT formalism that directly exploits the spatial localization of the perturbations under periodic boundary conditions has been reported in the literature, yet. This is particularly surprising, since real-space formalisms have attracted considerable interest for standard ground-state DFT calculations [53][54][55][56][57][58][59] in the last decades due to their favorable scaling with respect to the number of atoms and their potential for massively parallel implementations. Formally, one would expect a real-space DFPT formalism to exhibit similar beneficial features and thus to facilitate calculations of larger systems with less computational expense on modern multi-core architectures.
We here derive, implement, and validate a real-space formalism for DFPT. The inspiration for this approach comes from the work of Giustino et al. [52], who demonstrated that Wannierization [60] can be used to map reciprocal-space DFPT results to realspace, which in turn enables numerically efficient interpolation strategies [61]. In contrast to these previous approaches, however, our DFPT implementation is formulated directly in real space and utilizes the exact same localized, atom-centered basis set as the underlying ground-state DFT calculations. This allows us to exploit the inherent locality of the basis set to describe the spatially localized perturbations and thus to take advantage of the numerically favorable scaling of such a localized basis set. In addition, all parts of the calculation consistently rely on the same real-space basis set. Accordingly, all computed response properties are known in an accurate real-space representation from the start and no potentially error-prone interpolation (reexpansion) is required. However, this reformulation of DFPT also gives rise to many non-trivial terms that are discussed in this paper. For instance, the fact that we utilize atom-centered orbitals require accounting for various Pulay-type terms [62]. Furthermore, the treatment of spatially localized perturbations that are not translationally invariant with respect to the lattice vectors requires specific adaptions of the algorithms used in ground-state DFT to compute electrostatic interactions, electronic densities, etc. We also note that the proposed approach facilitates the treatment of isolated molecules, clusters, and periodic systems on the same footing. Accordingly, we demonstrate the validity and reliability of our approach by using the proposed real-space DFPT formalism to compute the electronic response to a displacement of nuclei and harmonic vibrations in molecules and phonons in solids.
The remainder of the paper is organized as follows. In Section 2 we succinctly summarize the fundamental theoretical framework used in DFT, in DFPT, and in the evaluation of harmonic force constants. Starting from the established real-space formalism for ground-state DFT calculations, we derive the fundamental relations required to perform DFPT and lattice dynamics calculations in Section 3. The practical and computational implications of these equations are then discussed in Section 4 using our own implementation in the all-electron, full-potential, numerical atomic orbitals based code FHI-aims [55,63,64] as an example. In Section 5 we validate our method and implementation for both molecules and extended systems by comparing vibrational and phonon frequencies computed with DFPT to the ones computed via finitedifferences. Furthermore, we exhaustively investigate the convergence behavior with respect to the numerical parameters of the implementation (basis set, system sizes, integration grids, etc.) and we discuss the performance and scaling with system size. Eventually, Section 6 summarizes the main ideas and findings of this work and highlights possible future research directions, for which the developed formalism seems particularly promising.

Density-functional theory
In DFT, the total energy is uniquely determined by the electron density n(r) (1) in which T s is the kinetic energy of non-interacting electrons, E ext the electron-nuclear, E H the Hartree, E xc the exchange-correlation, and E ion−−ion the ion-ion repulsion energy. All energies are functionals of the electron density. Here we avoid an explicitly spin-polarized notation, a formal generalization to collinear (scalar) spin-DFT is straightforward.

E KS = T s [n] + E ext [n] + E H [n] + E xc [n] + E ion−−ion ,
The ground state electron density n 0 (r) (and the associated ground state total energy) is obtained by variationally minimizing Eq. (1) δ δn whereby the chemical potential µ = δE KS /δn ensures that the number of electrons N e is conserved. This yields the Kohn-Sham single particle equationŝ for the Kohn-Sham Hamiltonianĥ KS . In Eq. (3)t s is the single particle kinetic operator,v ext the (external) electron-nuclear potential, v H the Hartree potential, andv xc the exchange-correlation potential. Solving Eq. (3) yields the Kohn-Sham single particle states ψ i and their eigenenergies ϵ i . The single particle states determine the electron density via in which f (ϵ i ) denotes the Fermi-Dirac distribution function.
To solve Eq. (3) in numerical implementations, the Kohn-Sham states are expanded in a finite basis set χ µ (r) using the expansion coefficients C µi . In this expansion, Eq. (3) becomes a generalized algebraic eigenvalue problem Using the bra-ket notation ⟨.|.⟩ for the inner product in Hilbert space, H µν denotes the elements ⟨χ µ |ĥ KS |χ ν ⟩ of the Hamiltonian matrix and S µν the elements ⟨χ µ |χ ν ⟩ of the overlap matrix.
Accordingly, the variation with respect to the density in Eq. (2) becomes a minimization with respect to the expansion coefficients C νi in which the eigenstates ψ i are constrained to be orthonormal. Typically, the ground state density n 0 (r) and the associated total energy E tot are determined numerically by solving Eq. (7) iteratively, until self-consistency is achieved.
To determine the force F I acting on nucleus I at position R I in the electronic ground state, it is necessary to compute the respective gradient of the total energy, i.e., its total derivative [65][66][67] In Eq. (8) we have used the notation ∂/∂R I to highlight partial derivatives. The first term in Eq. (8) describes the direct dependence of the total energy on the nuclear degrees of freedom. The second term, the so-called Pulay term [62], captures the dependence of the total energy on the basis set chosen for the expansion in Eq. (5). It vanishes for a complete basis set or if the chosen basis set does not depend on the nuclear coordinates, e.g., in the case of plane-waves. The last term vanishes, if Eq. (7) has been variationally minimized with respect to the expansion coefficients C µi to obtain the ground state total energy and density. That this holds true also in practical numerical implementations is demonstrated in Appendix A. However, for higher order derivatives of the total energy, e.g., the Hessian, the last term no longer vanishes since the forces are not variational with respect to the expansion coefficients C µi . Accordingly, a calculation of the Hessian does not only require the analytical derivatives appearing in the first two terms, but also the response of the expansion coefficients and the basis functions to a nuclear displacement (∂C µi /∂R J and ∂χ µ /∂R J , respectively). More generally, according to the (2n + 1) theorem, knowledge of the nth order response (i.e. the nth order total derivative) of the electronic structure with respect to a perturbation is required to determine the respective (2n + 1)th total derivatives of the total energy [8].
These response quantities are, however, not directly accessible within DFT, but require the application of first order perturbation theory.

Density-functional perturbation theory
To determine the ∂C µi /∂R J and ∂χ µ /∂R J needed for the computation of the Hessian (Eq. (9)), we assume that the displacement from equilibrium ∆R J only results in a minor perturbation (linear response) of the original Hamiltonianĥ (0) KS . We then expand the wave func- (1) i (∆R J ) linearly and apply the normalization condition we then immediately obtain the Sternheimer equation [68] ( The corresponding first order density is given by To solve the Sternheimer equation (Eq. (12)), we use the DFPT formalism [9,11] and thus the same expansion for ψ (1) i as used in To determine the unknown coefficients C (1) µi , it is necessary to iteratively solve Eq. (12) until self-consistency is achieved. This is best done in matrix form: Formally, DFPT and CPSCF are equivalent and only differ in the way the first order wave function coefficients C (1) are obtained. In the DFPT formalism, C (1) is calculated directly by solving Eq. (15) self-consistently. In the CPSCF formalism, the coefficients C (1) are further expanded in terms of the coefficients of the unperturbed system [12,13] whereby the respective expansion U pi coefficients are given by Here, the Ď is used to denote the respective Hermitian conjugate of the matrices, and E (0) denotes the diagonal matrices containing the eigenvalues ϵ i .

The harmonic approximation: Molecular vibrations and phonons in solids
DFPT is probably most commonly applied to calculate molecular vibrations or phonon dispersions in solids in the harmonic approximation, although its capabilities extend much beyond this [42]. Since we will later use vibrational and phonon frequencies to validate our implementation, we will now briefly present the harmonic approximation to nuclear dynamics.
To approximately describe the dynamics for a set of nuclei {R I }, the total energy Eq. (7) is Taylor-expanded up to second order around the nuclei's equilibrium positions {R 0 The linear term in this expansion is not noted because it vanishes at the equilibrium positions. The Hessian in the second term (often referred to as force constants) can be determined with DFPT as described in the previous section. The equations of motions for the nuclei in this potential E harm tot ({R I }) are analytically solvable and yield a superposition of independent harmonic oscillators for the displacements from equilibrium ∆R I (t) = R I (t) − R 0 I . In the complex plane, these displacements correspond to the real part of in which the complex amplitudes (and phases) A λ are dictated by the initial conditions; the eigenfrequencies ω λ and the individual components [e λ ] I of the eigenvectors e λ are given by the solution of the eigenvalue problem: for the dynamical matrix A technical complication arises for periodic solids, which are characterized by a translationally invariant unit cell defined by the lattice vectors a 1 , a 2 , and a 3 . Each of the N atoms R I in the primitive unit cell thus has an infinite number of periodic replicas whereby R m denotes an arbitrary linear combination of a 1 , a 2 , and a 3 (see Fig. 2). Accordingly, also the size of the Hessian becomes in principle infinite, since also vibrations that break the perfect translational symmetry need to be accounted for. This problem can be circumvented by transforming the harmonic force constants Φ harm Im,J into reciprocal space. Formally, this transforms this problem of infinite size into an infinite number of problems of finite size [69] since the finite (3N × 3N) dynamical matrix D(q) would in principle have to be determined for an infinite number of qpoints in the Brillouin zone. Its diagonalization would produce a set of 3N q-dependent eigenfrequencies ω λ (q) and -vectors e λ (q).
Furthermore, the displacements defined in Eq. (19) acquire an additional phase factor: In reciprocal-space DFPT implementations [9,10,36,47,70], perturbations that are incommensurate with the unit cell (q ̸ = 0) are typically directly incorporated into the DFPT formalism itself. For instance, a perturbation vector leads to a density response that is not commensurate with the primitive unit cell. By adding an additional phase factor to the perturbation the translational periodicity of the unperturbed system can be restored so that also q ̸ = 0 perturbations become tractable within the original, primitive unit cell, which is computationally advantageous. However, one DFPT calculation for each q point is required in such cases. In our implementation, we take a different route by choosing a real-space representation, as discussed in detail in the next section.

Total energies and forces in a real-space formalism
In practice, FHI-aims uses the Harris-Foulkes total energy functional [71,72] In both Eq. (29) and here, Z I is the nuclear charge, and n MP (r) the multipole density obtained from partitioning the density n(r) into individual atomic multipoles to treat the electrostatic interactions in a computationally efficient manner. Accordingly, is the full electrostatic potential stemming from atom I, which includes the electronic and nuclear contributions. The respective forces can be split into three individual terms. The Hellmann-Feynman force is The Pulay force can be written as and the force arising from the multipole correction is

Periodic boundary condition
To treat extended systems with periodic boundary conditions in a real-space formalism, the equations for the total-energy and the forces given in the previous section need to be slightly adapted. The general idea follows this line of thought: A periodic solid is characterized by a (not-necessarily primitive) unit cell that contains atoms at the positions R I , whereby the lattice vectors a 1 , a 2 , a 3 characterize the extent of this unit cell and impose translational invariance. To compute the properties of such a unit cell, it is not sufficient to only consider the mutual interactions between the electronic density n(r) and atoms R I in the unit cell, but it is also necessary to account for the interactions of the N uc atoms in the unit cell with the respective periodic images of the atoms R Im and of the density n(r +R m ) = n(r), as introduced and discussed in Eq. (22). Accordingly, the double sum in Eq. (29) and the single sum in Eq. . (37) Given that the extent of our atom-centered basis set is confined [55], only a finite number of periodic images needs to be . The blue dashed line shows the maximum extent of its orbitals. To treat periodic boundary conditions in DFT in real space, it is necessary to construct a supercluster (red solid line) which includes all periodic images that have nonvanishing overlap with the orbitals of the atoms in the original unit cell, as exemplarily shown here for atom A and B. In practice, it is sufficient to carry out the integration in the unit cell alone, since translational symmetry then allows to reconstruct the full information, as discussed in more detail in Sections 3.2 and 4. In turn, only the dark gray atoms that have non-vanishing overlap with the unit cell need to be accounted for in the integration, as exemplarily shown here for atom C. The DFPT supercell highlighted in black is the smallest possible supercell that encompasses the DFT supercluster and exhibits the same translational Born-von Kármán periodicity as the original unit cell. Accordingly, it contains slightly more atoms than the DFT supercluster, e.g., atom D. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) accounted for in these sums, since only a finite number of periodic images feature atomic orbitals that have non-zero overlap with the orbitals of the atoms in the unit cell, as sketched in Fig. 3. In practical calculations, these periodic images are accounted for explicitly by the construction of superclusters that encompass all N sc atoms with non-vanishing overlap with the orbitals of the N uc atoms in the original unit cell (see Fig. 3). As discussed in detail in Ref. [55,73], also the basis set needs to be adapted to reflect the translational symmetry. Since each local atomic orbital χ µ (r) in Eq. (5) is associated with an atom I(µ), we first introduce periodic images χ µm (r) = χ µ (r − R I(µ) + R m ) for them as well. Following the exact same reasoning as in Section 2.3, the atomic orbitals used for the expansion of the eigenstates (5) are then replaced by Bloch-like generalized basis functions so that all matrix elements ⟨.|.⟩ become k-dependent, e.g.,  uc χ µm (r)ĥ ks χ νn (r)dr.
Please note that for practical reason the integration has been restricted to the unit cell (uc) in this case. To reconstruct the full information, e.g., of the N uc × N sc overlap matrix, the double sum and the associated phase factors run over all periodic images N sc × N sc , whereby only atoms with non-vanishing overlap in the unit cell contribute (see Fig. 3 and Ref. [73]). These sums are finite, since Table 1 Number of atoms in the unit cell and the corresponding number of atoms in the supercluster used in the ground-state DFT calculations (atoms in the red box in Fig. 3) and in the DFPT supercell (black box in Fig. 3). Please note that in the case of Si the increased number of atoms in the DFPT supercell originates from the fact that in this case the circle-like DFT supercluster is encompassed by an oblique DFPT supercell with the same shape as the primitive unit cell of the diamond structure. all basis functions are bounded by a confinement potential [55]. In the expression for the Kohn-Sham energy (Eq. (29)) and the Pulay force (Eq. (35)), the sum over electronic states now also runs over the N k k-points Formally, the infinite periodic solid is thus treated in real-space within a finite DFT ''supercluster'' (see Fig. 3), which explicitly includes all periodic images R Im that have non-vanishing orbital overlap with the unit cell. Thereby, Eq. (38) enforces the translational symmetries to be retained. Accordingly, this real-space formalism for periodic solids leads to a notable, but reasonably tractable computational overhead for DFT calculations, e.g., when comparing calculations with N primitive atoms in a unit cell to calculations with N atoms in an isolated molecule. This becomes immediately evident from Table 1, which lists some typical supercell sizes that are used in the ground state total energy calculations at the DFT level for representative 1D, 2D, and 3D systems. However, the fact that the underlying DFT formalism explicitly accounts for all periodic images R Im turns out to even be advantageous in DFPT calculations. For instance, the computation of the dynamical matrix in Eq. (23) explicitly requires the derivatives with respect to all periodic replicas R Im . As discussed in detail in the Section 3.3, the real-space formalism allows to reconstruct all the necessary, non-vanishing elements of the Hessian that enter Eq. (23) within one DFPT run. In turn, this allows us to exactly compute the dynamical matrix (Eq. (23)) -and thus all eigenvalues ω 2 λ (q) andvectors e λ (q) -at arbitrary q-points by simple Fourier transforms.
In practice, we achieve this goal by computing the Hessian in a slightly larger Born-von Kármán [69] DFPT supercell that encompasses the supercluster used for DFT ground state calculations (cf. Fig. 3). By these means, the minimum image convention associated with translational symmetry can be straightforwardly exploited also in the case of perturbations that break the original symmetry of the crystal.
It should be noted that, for semiconductors and insulators, the size of the DFPT supercell is typically determined by the extent of the orbitals. However, for metals, this may not be enough since a large number of k-points is required for convergence. To be consistent with this finer k-mesh, the DFPT supercell would have to be extended to a much larger size for metals. The traditional reciprocal space approach [9][10][11] might therefore be computationally advantageous for metal. For this reason, we only apply our real-space formalism to semiconductors and insulators in the following sections.

Real-space force constants calculations
To derive the expressions for the force constants in realspace, we will directly use the general case of periodic boundary conditions, as introduced in the previous section. Analogously to Eq. (33) we can split the contributions to the Hessian (or to the force constants) defined in Eq. (9) into the respective derivatives of the contributions to the force Φ harm Please note that we have omitted the multipole term here, since its contribution is already three orders of magnitude smaller at the level of the forces.
Due to the permutation symmetry (Φ Is,J = Φ J,Is ) of the force constants, the order in which the derivatives are taken is irrelevant. The formulas given above for the forces F I acting on the atoms in the unit cell are equally valid for the forces F Is acting on its periodic images R Is , as long as the sums and integrals in the supercell (see Fig. 3) are performed using the minimum image convention. In the following, we will exploit this fact so that only total derivatives with respect to the atoms in the primitive unit cell need to be taken. Consequently, the total derivative of the Hellmann-Feynman force yields in which δ Is,J0 = δ IJ δ s0 denotes a multi-index Kronecker delta.
To determine the total derivative of the Pulay force, we first split Eq. (35) into two terms using the density matrix and the energy weighted density matrix (46) which also incorporate the phase factors arising due to periodic boundary conditions. Using this notation, the total derivative of the Pulay term can be split into four terms for the sake of readability: The first term account for the response of the Hamiltonianĥ ks (k), while the third and fourth term account for the response of the energy weighted density matrix W µm,νn and the overlap matrix S µm,νn , respectively (cf. Section 4.1). Please note that in all four contributions many terms vanish due to the fact that the localized atomic orbitals χ µm (r) are associated with one specific atom/periodic image R J(µ)m , which implies, e.g., This allows us to re-index the sums over (µm, νn) in a computationally efficient, sparse matrix formalism (cf. Ref. [74]). Similarly, it is important to realize that all partial derivatives that appear in the force constants can be readily computed numerically, since the χ µm are numeric atomic orbitals, which are defined using a splined radial function and spherical harmonics for the angular dependence [55].

Details of the implementation
The practical implementation of the described formalism closely follows the flowchart shown in Fig. 4. For the sake of readability we use the notation to highlight that in each step of the flowchart a loop over all atoms in the unit cell R I viz. all periodic replicas R Is is performed to compute all associated derivatives. In the following chapters, we will use subscripts i, j for occupied KS orbitals in the DFPT supercell, and a for the corresponding unoccupied (virtual) KS orbitals, and p, q for the entire set of KS orbitals in the DFPT supercell.
After the ground state calculation (see Section 2.1 and Ref. [55]) is completed, the first step is to compute the response of the overlap matrix S (1) . We then use U (1) ai = 0 Appendix B as the initial guess for the response of the expansion coefficients and determine the response of the density matrix P (1) , which then allows to construct the respective density n (1) (r). Using that, we compute the associated response of the electrostatic potential and of the Hamiltonian H (1) . In turn, all these ingredients then allow to set up the Sternheimer equation, the solution of which allows to update the response of the expansion coefficients C (1) . Using a linear mixing scheme, we iteratively restart the DFPT loop until self-consistency is reached, i.e., until the changes in C (1) become smaller than a user-given threshold. In the last steps, the response of the energy weighted density matrix W (1) , the force-constants Φ Im,J , and the dynamical matrix D(q) are computed and diagonalized on userspecified paths and grids in reciprocal space.

Response and Hessian of the overlap matrix
The first step after completing the ground state DFT calculation is to compute the first order response of the overlap matrix, a quantity that is not required in plane-wave implementations, but that needs to be accounted for when using localized atomic orbitals [62]. Using the definition of the overlap matrix S given in Therefore, it is possible to restrict the integration to the unit cell (uc) and to reconstruct the whole integral by summing over all periodic replicas n, as illustrated in Fig. 5.
For the response of the overlap matrix, translational symmetry enables us again to restrict the integration to the unit cell as illustrated in Fig. 6. Please note that only very few nonvanishing contributions exist, since every orbital only depends on the position of one specific atom or replica Following the same strategy, also the second order derivatives of the overlap matrix required in Eq. (54) can be computed using: Again, only a few contributions exist for the first term and for the second term

Response of the density matrix
The first step in the DFPT self-consistency cycle is to calculate of the response of the density matrix using the given expansion coefficients C (0) and C (1) . Using the discrete Fourier transform to get real-valued coefficients C (0) µm,i , the density matrix defined in Eq. (45) becomes: Accordingly, its response is In the practical solution of the Sternheimer equation (cf. Section 4.6), we use the CPSCF approach (Eq. (16)) and use matrix U (1) to expand the response of the expansion coefficients C (1) We have also solved the Sternheimer equation use DFPT approach (Eq. (15)) directly, and obtained exactly the same results Integration strategy for the computation of the response matrix elements, here shown for the first order overlap matrix S (1) in Eq. (60). Please note that to be able to restrict the integration to the unit cell, the derivative has to be translated together with the orbital as shown in Eq. (59).
as with Eq. (16) for the systems (e.g. molecules) discussed in this paper. In praxis, the density matrix can then be directly evaluated in terms of U (1) , as shown in Appendix B.

Response of the electronic density
To determine the electronic density n(r), we use a density matrix based formalism Similarly, the response of the electronic density can thus be expressed as µm,νn χ (0) µm (r)χ (0) νn (r) Please note that the ground state density is periodic (translationally invariant) but its response is not As already discussed for the response of the overlap matrix in Section 4.1, the individual contributions to the response are however related to each other via their translation property

Response of the total electrostatic potential
In a real-space formalism [53,55] such as FHI-aims it is necessary to treat the electrostatic interactions (electronic Hartree potential v es and nuclear external potential v ext in a unified formalism [55,73]. Using Eq. (31), the electrostatic potential entering the zero-order Kohn-Sham Hamiltonianĥ The contribution of each atom R Jn consists of two contributions V es,tot Jn In this expression denotes the electrostatic potential associated with an isolated (''free'') atom of the same species with the electron density n free (r− R Jn ). Both the free-atom electron densities n free (r − R Jn ) and the electrostatic potential V free Jn (r − R Jn ) are accurately known as cubic spline functions on dense grids. The second term in the total electrostatic potential V es,tot Jn is computed by partitioning [73] the difference density δn(r) = n(r) −  J,n n free (r − R Jn ) into individual contributions δ I n(r). Their contribution δV Jn (r − R Jn ) to the translationally invariant and periodic electrostatic potential is computed using a combined multipole expansion and Ewald summation formalism proposed by Delley [53].
As the perturbations break the local periodicity of the crystal, also, their response is localized in non-polar materials [52]. Accordingly, no Ewald summation is needed for the response potential. Instead, we use a real-space multipole expansion for the computation of the first order potential V (1) es,tot (r). From the given first-order density n (1) (r), we first construct whereby n free (r − R Is ) and its first derivative is available by splines [55]. The respective first order potential thus becomes The first term is readily accessible, given that V free (r − R Is ) is accurately known as a cubic spline. For the second term, we first partition δn (1) into individual contributions stemming from the different atoms and periodic replicas R Is , so we have the radial part of density: Here the upper index (lm) refers to the quantum numbers of the spherical harmonics. The p J (r) are the atom-centered partition functions [55]. From that, we get the radial part of the electrostatic potential: Here, g l (r < , r > ) = r l < /r l+1 > is the Green function for the unscreened Hartree potential [55]. Then the full electrostatic potential is reassembled using and δV (1) Please note that the chosen approach is valid to describe the electrostatics in non-polar materials, in which the perturbation of the electrostatic potential is indeed spatially localized [52]. Accordingly, it can be treated accurately within the finite supercells used in our real-space DFPT approach (see Section 3). Exemplarily, this is demonstrated in Fig. 7 for the response of the electrostatic potential computed in a one-dimensional, infinite chain of polyethylene (C 2 H 4 ). In polar materials, long-ranged dipole interactions can arise, which would extend beyond the boundaries of the DFPT supercells used in the real-space formalisms. In that case, additional correction terms to the electrostatic perturbation potential [75] need to be accounted for.

Response of the Kohn-Sham Hamiltonian
To determine the Hamiltonian matrix and its response, we again exploit their properties under translations already discussed for the overlap matrix in Section 4.1: The response of the Hamiltonian operator includes the response of the total electrostatic potential V (1) es,tot discussed in the previous section and the response of the exchangecorrelation potential V (1) xc . In the case of the LDA [76,77] functional considered in this work, evaluating the functional derivative in the latter term yields: A sketch of the employed integration strategy to compute the response of the Hamiltonian is shown in Fig. 8.

Solution of the Sternheimer equation
Using the notations introduced in this section, the Sternheimer equation given defined in Eq. (15) becomes The calculation was performed at the LDA level of theory using fully converged numerical parameters (cf. Section 5.1). In this non-polar system, the response of the electrostatic potential is strongly localized at the perturbation and thus contained in the DFPT supercell used in the calculation (cf. Fig. 3 and Table 1).
µm,ν0 . The first row (a) shows the ground-state Kohn-Sham Hamiltonian, which -due to its periodicity -can be integrated using the exact same strategy used for the overlap matrix S (0) (see Fig. 5). The remaining rows (b) highlight that the response H (1) µm,ν0 requires to account for derivatives of the Kohn-Sham Hamiltonian dĥ KS /dR Is , which is not periodic. To restrict the integration to the unit cell, it is thus necessary to translate also this perturbation accordingly. For this exact reason, a Born-von Kármán supercell [69] supercell is needed in DFPT, but not in the case of a periodic Hamiltonian as in DFT.
More conveniently, it can be written in matrix form as whereby E (0) and E (1) denote the diagonal matrices containing the eigenvalues ϵ i and their responses respectively. By multiplying with the Hermitian conjugate C (0)Ď and by expanding the response C (1) in terms of the zero-order expansion coefficients C (0) using we get Thereby, we have used the orthonormality relation: Due to the diagonal character of E (0) and E (1) , this matrix equation contains the response of the eiqenvalues on its diagonal Conversely, the off-diagonal elements determine the response of the expansion coefficients for p ̸ = q The orthogonality relation then also yields the missing diagonal elements

Response of the energy weighted density matrix
After achieving self-consistency in the DFPT loop, the last task is to determine the response of the energy weighted density matrix that is required for the evaluation of Eq. (53). In close analogy to the density matrix formalism discussed in Section 4.2, the response of the energy weighted density matrix can be expressed as: In close analogy to our discussion of the density matrix, the energy weighted density matrix is also evaluated in practice directly in terms of U (1) , as detailed in Appendix C.

Symmetry of the force constants
As mentioned above, the individual force constant elements are related to each other by translational symmetry Φ Is,J0 = Φ I(s+m),Jm , and permutation symmetry Φ Is,Jm = Φ Jm,Is .
Due to these symmetries, only a subset N uc × N sc of the complete N sc ×N sc force constant matrix needs to be computed for a supercell containing N sc atoms (see Fig. 3 and Table 1). Similarly, invariance which enables us to determine the entries on the diagonal Φ J0,J0 from the off-diagonal elements. For our implementation, this is computationally favorable, since no special treatment of ''on-site'' terms, i.e., contributions stemming from one individual atom, is required, e.g., in Eq. (42) or for the integration of ''on-site'' matrix elements [73]. Please note that space and point group symmetries [69], which would allow to further reduce the amount of force constants that need to be computed, are not exploited in the implementation, yet.

Validation and results
To validate our implementation we have specifically investigated the convergence of vibrational frequencies with respect to the numerical parameters used in the calculation in Section 5.1. Furthermore, a systematic validation of the implementation by comparing to vibrational frequencies obtained from finitedifferences is presented in Section 5.2; these tests are extended to periodic systems in Section 5.3. All benchmark data is available in the NoMaD Repository (https://repository.nomad-coe.eu) via http://dx.doi.org/10.17172/NOMAD/2017.02.19-1. Eventually, the computational performance of the implementation is discussed in Section 5.4.

Convergence with respect to numerical parameters
First, we analyze the convergence behavior of our DFPT implementation with respect to the numerical parameters used in the calculation, i.e., the basis set size used in the expansion (Eq. (5)) of the Kohn-Sham states in numerical, atom centered orbitals and the (radial and angular) grids used for the numerical integration. Exemplarily, we discuss these effects using the six infrared active frequencies of ethane (C 2 H 6 ), which in all cases are computed using a local approximation for exchange and correlation (LDA parametrization of Perdew and Zunger [76] for the correlation energy density of the homogeneous electron gas based on the data of Ceperley and Alder [77]). In all cases, the DFPT calculations were performed for the respective equilibrium geometry, i.e., the structure obtained by relaxation (maximum force < 10 −4 eV/Å) using the exact same computational settings. Due to the fact that the exact same formalism is used for both for finite systems and We use a ''tier 2'' basis set and N ang,max = 590 here. The benchmark values are calculated using N r,mult = 3. Fig. 11. Convergence of the infrared-active vibrational frequencies of ethane with respect to the angular integration grid, as controlled by the parameter N ang,max (see text). We use a ''tier 2'' basis set and N r,mult = 2 here. The benchmark values are calculated using N ang,max = 590. periodic materials, the presented convergence studies are also valid for both cases. Fig. 9 shows the absolute change in these vibrational frequencies if the basis set size is increased. Here, a minimal basis (half a basis function per electron in the spin-unpolarized case) includes the orbitals that would be occupied orbitals in a free atom following the Aufbau principle. Additional sets of basis functions are added in ''tier 1'', ''tier 2'', . . . , calculations, see Ref. [55] for more details. The vibrational frequencies converge quickly with the basis set size. Already at a ''tier 1'' level we get qualitatively correct results with a maximal absolute/relative error of 18 cm −1 /0.6%.
Fully quantitatively converged calculations are achieved with the ''tier 2'' basis set.
Atom-centered grids are used for the numerical integrations in FHI-aims [55]: Radially, each atom-centered grid consists of N r spherical integration shells, the outermost of which lies at a distance r outer from the nucleus. The shell density can be controlled by means of the radial multiplier N r,mult . For example, N r,mult = 2 results in a total of 2N r + 1 radial integration shells. On these shells, angular integration points are distributed in such a way that spherical harmonics up to a certain order are integrated exactly by the use of Lebedev grids as proposed by Delley [78]. Here, we characterize the angular integration grids by the maximum number of angular integration points N ang,max used in the calculation. Figs. 10 and 11 show our convergence tests with respect to N r,mult and N ang,max , respectively. In both cases, we find that the computed vibrational frequencies depend only weakly on the chosen integration grids: For N r,mult , even the most sparse radial integrations grids yields qualitative and almost quantitatively correct frequencies, since the maximum absolute and relative errors are 5.5 cm −1 and 1.8%, respectively.
Quantitatively converged results are achieved at the N r,mult = 2 level with absolute and relative errors of 0.2 cm −1 and 0.08%. As Fig. 11 shows, the vibrational frequencies are virtually unaffected by the angular integration grids; the maximum absolute error is always smaller than 0.01 cm −1 .

Validation against finite-differences
To validate our DFPT implementation, we have compared the obtained vibrational frequencies to finite-difference calculations, in which the Hessian was obtained via a first order finite-difference expression for the forces and dipole moments (see below) using an atomic displacement of 0.0025 Å. Exemplarily, we discuss the performance of our implementation using the infrared (IR) spectrum of the C 60 molecule. The IR intensity for a given vibrational eigenmode e λ can be computed both with finite-differences and DFPT by inspecting the changes induced in the dipole moment µ =  n(r) r dr by the displacements associated with the vibrational mode λ. As Fig. 12 illustrates, both the IR frequencies and intensities agree very well between the finite-difference approach and our DFPT implementation. To validate our DFPT implementation in a more systematic way, we have also compared the vibrational frequencies of 32 selected molecules with finite-difference calculations, utilizing the exact same first order finite-difference formalism used for the C 60 molecule. All calculations were performed at the LDA level of theory using fully converged numerical parameters 2 for the equilibrium geometry determined by relaxation (maximum force < 10 −4 eV/Å). A detailed list of results for these calculations is given in the Appendix D. For the sake of readability, we here only discuss the difference between the vibrational frequencies obtained via DFPT and via finite-differences, which we quantify by the mean absolute error (MAE), the maximum absolute error (MaxAE), the mean absolute percentage error (MAPE) and the maximum absolute percentage error (MaxAPE) for each molecule. These statistical data is succinctly summarized in Table 2: Overall, we find an excellent agreement between our DFPT implementation and the finite-difference results (average MaxAE of 1.40 cm −1 and average MaxAPE of 0.16%). Please note that the largest occurring absolute error (10.13 cm −1 in P 2 ) and the largest occurring relative error (1.46% in H 2 O 2 ) still correspond to relatively moderate relative and absolute errors (1.26% and 5.73 cm −1 , respectively). The occurrence of these deviations are in part caused by numerical errors, e.g., the ones arising due to the moving integration grid [55] and due to the finite multipole expansion [55] (the multipole term in force constants calculation Eq. (41) has been omitted). Such errors affect these two approaches (finite difference and DFPT) differently. To a large extent, this is mitigated in these 2 So called ''tier 2'' basis sets and ''really tight'' defaults were used for the numerical settings. Additionally, we increased the order of the multipole expansion to l = 12 and the radial integration grid to N r,mult = 4 for all systems except LiF, NaCl, and P 2 . An atomic displacement of 0.013 Å was used in the finite-difference calculations. benchmark calculations by choosing highly-accurate settings. Still, the finite-difference reference calculations themselves exhibit a certain uncertainty, since they can be sensitive to the atomic displacement chosen for evaluating the numeric derivatives. For instance, this is the case for the P 2 molecule, which exhibits the largest absolute error in Table 2. For this reason, we have also compared our DFPT calculations with benchmark results (Gaussian code, aug-cc-pVTZ basis set) reported in the ''NIST Computational Chemistry Comparison and Benchmark Database'' [79]. For the 15 dimers contained both in Table 2 and in this database, the mean absolute percentage errors is only 0.5%.

Extended systems: Phonons
To showcase the ability of our implementation to treat finite systems and periodic solids on the same footing, we compare the vibrational frequencies of various polyethylene chains H(C 2 H 4 ) n H with different lengths (n from 1 to 8) to the respective periodic, infinite chain of C 2 H 4 . In the latter case, we compute the vibrational/phonon density of states (DOS) whereby a normalized Gaussian function with a width σ of 5 cm −1 is used to approximate the Delta-distribution δ[ω − ω λ (q)]. It should be noted that the phonon DOS of an infinite C 2 H 4 chain is not zero at the Γ -point, because it is a one-dimensional system [80]. All calculations have been performed for relaxed equilibrium geometries (maximum force < 10 −4 eV/Å) with fully converged numerical parameters, i.e., using the aforementioned really-tight integration grids and ''tier 2'' basis sets. For the periodic chain, a reciprocal-space grid of 11 × 1 × 1 electronic k-points and a grid of 200 × 1 × 1 vibrational q-points (in the primitive Brillouin zone) has been utilized to converge the density of states g(ω), as substantiated in Figs. 13 and 14. Whereas the convergence with respect to electronic k-points is reasonably fast, a large amount of vibrational q-points is required to sample the Brillouin zone, especially for the relatively moderate broadening σ of 5 cm −1 . In this context, it is important to realize that the Table 2 Mean absolute error (MAE), maximum absolute error (MaxAE), mean absolute percentage error (MAPE) and max absolute percentage error (MaxAPE) for the difference between the vibrational frequencies obtained via DFPT and via finitedifferences using an atomic displacement of 0.013 Å for a set of 32 molecules. All calculations are performed at the LDA level of theory with fully converged numerical settings and relaxed geometries (see text and respective footnote). actual number of q-points used is not at all computationally critical in our implementation: As discussed in Section 2.3, our implementation involves determining all non-vanishing forceconstants in real-space; the respective q-dependent properties can then be determined exactly by a simple Fourier transform with minimal numerical effort. For instance, using q = 2000 only requires ∼1 s more computational time than the q = 20 case.

MAE (cm
The outcome of these investigations is summarized in Fig. 15, in which the vibrational density of states (σ = 1 cm −1 ) for the isolated H(C 2 H 4 ) n H chains with variable length (n from 1 to 8) is compared to the vibrational density of states (σ = 5 cm −1 ) of the extended, infinitely long polyethylene (C 2 H 4 ) chain. With increasing length n, the vibrational frequencies of the isolated chain start to resemble the density of states g(ω) of the infinitely long polyethylene chain. Still, some features, e.g., the low frequency modes that stem from long-wavelength phonons can only be correctly captured in the periodic DFPT calculation. Please note that the differences between the vibrational density of state of the H(C 2 H 4 ) 8 H molecule (50 atoms) and the C 2 H 4 chain (66 atoms in the DFPT supercell) are to a large extend not caused by the additional force-constants accounted for in the periodic case. Rather, the differences stem from the fact that the molecular vibrational density of states effectively corresponds to a reciprocalspace sampling of q ≈ 8, which -as Fig. 14 shows -is not sufficient to capture the contributions of long-range wavelengths to the density of states. Eventually, we have also validated our real-space implementation against finite-difference calculations performed using phonopy [81,82] for two realistic periodic systems. As a twodimensional example, we use graphene, the vibrational properties of which have been controversially debated in the literature e Fig. 13. Convergence of the phonon density of states of polyethylene with respect to the number of k-points utilized in the primitive Brillouin zone for DFPT calculations of the C 2 H 4 chain. The top panel shows the density of states for 18 kpoints and the bottom panel shows the difference with respect to this converged reference. A Gaussian broadening of 5 cm −1 and 200 q points was used in the computation of g(ω). [83,84], especially regarding the role of long-ranged interactions that are not treatable in real-space. As discussed in Section 4.4 already, correction terms that can account for such interactions are not yet part of the implementation discussed in this work. To avoid possible artifacts due these effects, we have thus performed finitedifference calculations (displacement 0.008 Å) in the exact same 11 × 11 × 1 supercell (242 atoms) that is also inherently used in the DFPT calculations itself (see Fig. 3). In both the case of DFPT and finite-differences, all calculations have been performed for relaxed equilibrium geometries (maximum force <10 −4 eV/Å) with 11 × 11 × 1 k-points in the primitive unit cell, tight settings, the ''tier 1'' basis set, and the LDA functional. As Fig. 16 shows, we find an excellent agreement between our DPFT implementation and the finite-difference calculations: By using such extended supercells, even the parabolic dispersion [85] in the lowest acoustic branch and the Kohn anomalies at Γ and K are captured in a qualitatively correct fashion by both our real-space DFPT and the finite-difference approach, as shown by Maultzsch et al. [83] before. Our implementation is thus ideally suited to further investigate to which extent long-range corrections to the perturbation potential will alter these effects.
For a three-dimensional system, we have used silicon in the diamond structure as an example. All calculations have been performed for relaxed equilibrium geometries (maximum force  Fig. 13, but for the convergence with respect to the number of q-points in the primitive Brillouin zone. A Gaussian broadening of 5 cm −1 is used. < 10 −4 eV/Å) with 7 × 7 × 7 k points in the primitive Brillouin zone, tight settings for the integration, a ''tier 1'' basis set, and the LDA functional. Finite-difference calculations have been performed again using phonopy [81,82] with a 5 × 5 × 5 supercell of the conventional cubic fcc cell (1000 atoms) and a finite displacement of 0.01 Å, which yields fully converged vibrational band structures (error < 1 cm −1 ). This was systematically checked by running finite-difference calculations for up to 9 × 9 × 9 supercells of the primitive unit cell (1458 atoms). As shown in Fig. 17, our DPFT implementation again yields an excellent agreement with the respective finite-difference calculations.

Performance and scaling of the implementation
To systematically investigate the performance and scaling of our implementation, we here show timings for the H(C 2 H 4 ) n H molecules with variable length n = 1-90 and the polyethylene chain C 2 H 4 . In the latter case, we have systematically increased the number of building units in the unit cell from (C 2 H 4 ) 1 to (C 2 H 4 ) 12 . All calculations use a ''tier 1'' basis set, light settings for the integrations, and the LDA functional. 11 × 1 × 1 k-points were used to sample the primitive Brillouin zone in the periodic case. We performed all these calculations on a single node featuring two Intel Xeon E5-2698v3 CPUs (32 cores) and 4 Gb of RAM per core.
For the timings of the finite molecules shown in Fig. 18, we find that the integration of the Hamiltonian response matrix H (1) determines the computational time for small system sizes i.e., for Density of States (a.u.) less than 200 atoms. As it is the case for the update of the response density n (1) , which involves similar numerical operations, we find a scaling of O(N 2 ) for this step (see Table 3). This is not too surprising, since these operations, which scale with O(N) at the ground-state DFT level [55], need to be performed 3N times when assessing the Hessian at the DFPT level, i.e., once for each cartesian perturbation of each atom. For the exact same reasons, the treatment of electrostatic effects, which scales as ∼O(N 1.6 ) at the ground-state DFT level [55], scales as O(N 2.4 ) for the computation of the electrostatic response potential V     Fig. 4, also the timings required for the computation of the individual response properties (density n (1) , electrostatic potential V (1) es,tot , Hamiltonian matrix H (1) , density matrix P (1) ) are given. Here we use light settings for the integration, a ''tier 1'' basis set, and the LDA functional. systems (N ≪ 1000), it would thus be beneficial to switch to a more advanced formalism for this computational step [16,17].
To understand the timings shown in Fig. 19 for the periodic linear chain, it is important to realize that such periodic calculations do not directly scale with the number of atoms N, as it was the case in the finite system, in which an N × N Hessian was computed. Rather, the calculations are inherently performed in a supercell (see Fig. 3) that features N sc atoms in total. As discussed in Section 2.3, only an N × N sc subsection of the Hessian needs to be determined. Accordingly, the scaling is thus best rationalized as function of the effective number of atoms N eff = √ N · N sc , as shown in Fig. 19 and Table 3. In this representation, the scaling and the respective exponents closely follow the behavior discussed for the finite systems already with one exception: Due to the fact that a sparse matrix formalism is used in the periodic implementation (see Section 3.3 and Ref. [74]), a more favorable scaling for the construction of the density matrix response P (1) is found.
As also shown in the lower panel of Fig. 19 and Table 3, the scaling does however not follow these intuitive expectations if plotted with respect to the number of atoms N present in the primitive unit cell, since N eff , N sc , and N are not necessarily linearly related. For the case of the linear chain, the number of periodic images N sc −N with atomic orbitals that reach into the unit cell should be a constant that is independent of the chain length viz. number of atoms N present in the unit cell. Accordingly the ratio N sc /N decreases from a value of 9 in the primitive C 2 H 4 unit cell (6 atoms) to a value of N sc /N = 3, if a (C 2 H 4 ) 4 unit cell with 24 atoms is used.
In this regime, in which N eff is approximately proportional to √ N, we find a very favorable overall scaling of O(N 1.3 ), whereby neither of the involved steps scales worse than O(N 1.7 ).
For larger system sizes (N > 24), however, the scaling deteriorates. The reason for this behavior is the rather primitive and simple strategy that we have employed in the generation of the DFPT supercells to facilitate the treatment of integrals using the minimum image convention, as discussed in Section 3.2. Effectively, these supercells are constructed using fully intact, translated unit cells -even if a considerable part of the periodic atomic images contained in this translated unit cell do not overlap with the original unit cell. For the case of the linear chain, the minimal possible ratio N sc /N = 3 is thus reached in the N = 24 case and retained for all larger systems N > 24. In this limit, N eff becomes proportional to N, so that we effectively recover the scaling exponents found for N eff and for finite molecular systems (cf. Table 3). In summary, we find an overall scaling behavior that is always clearly smaller than O(N 3 ) for the investigated system sizes both in the molecular and the periodic case. For the periodic case, we find a particularly favorable scaling regime of O(N 1.3 ) for small to medium sized unit cells N 24. As discussed in more detail in the outlook, this regime can be potentially improved and extended to larger unit cell sizes. Please note that the scaling relations discussed above for the linear chain are qualitatively also found in the case of 2D and 3D materials. Given that the utilized atomic orbitals are spatially confined within a cut-off radius [55], similar relations between N sc and N are effectively found in the case of graphene and silicon. Although the prefactors depend on the shape and dimensionality of the unit cell, the relation N eff ∝ √ N also approximately holds in these cases. In this context it is very gratifying to see that even quite extended systems (molecules with more than 100 atoms and periodic solids with more than 50 atoms in the unit cell) are in principle treatable within the relatively moderate CPU and memory resources offered by a single state-ofthe-art workstation.
Eventually, let us note that a parallelization over cores viz. nodes is already part of the presented implementation, given that the discussed real-space DFPT formalism closely follows the  Fig. 4, also the timings required for the computation of the individual response properties (density n (1) , electrostatic potential V (1) es,tot , Hamiltonian matrix H (1) , density matrix P (1) ) are given. Here we use light settings for the integration, a ''tier 1'' basis set, and the LDA functional. strategies used for the parallelization of ground-state DFT calculations in FHI-aims [55,63]: The parallelization of the operations performed on the real-space grid closely follows the strategy described in [63]; For the matrix operation, MPI based ScaLapack routines have been used to achieve a reasonable performance both regarding computational and memory parallelization.
The parallel scalability for a unit cell containing 1024 Si atoms is shown in Fig. 20. All calculations use a ''tier 1'' basis set, light settings for the integrations, and the LDA functional. One k-point is sufficient to sample the reciprocal space due to the large unit cell. Here we give the CPU time required for one single perturbation (one atom and one cartesian coordinate). Clearly, almost ideal scaling is achieved.

Conclusion and outlook
In this paper, we have derived and implemented a reformulation of density-functional perturbation theory in real-space and validated the proposed approach by computing vibrational properties of molecules and solids. In particular, we have shown that these calculations can be systematically converged with respect to the numerical parameters used in the computation. Also, we have demonstrated that the computed vibrational frequencies are essentially equal to those obtained from finite-differences -both for finite molecules and extended, periodic systems. Comparison of our results with vibrational frequencies stemming from different  (1) , density matrix P (1) ) are also given. Then red line corresponds to ideal scaling. The parallel efficiency is shown in the lower panel. Here we use light settings for the integration, a ''tier 1'' basis set, and the LDA functional. codes and implementations is urgently needed, but would go beyond the scope of this work.
The key idea of the proposed approach relies on the localized nature of the response density in non-polar materials, which enables the treatment of perturbations directly in real-space. On the one hand, this allows utilizing the computationally favorable real-space techniques developed over the last decades, e.g., massively parallel grid operations that scale O(N) [55,63]. On the other hand, the proposed approach allows us to determine the full, nonvanishing response in real-space in one DFPT run. In turn, simple and numerically cheap Fourier transforms -without the need of invoking any Fourier interpolation -give access to the exact associated response properties in reciprocal-space. We have explicitly demonstrated the viability of this approach for lattice dynamics calculations in periodic systems: In that case, we get fully q-point converged densities of states and vibrational band structures along arbitrary paths from one DFPT run in real-space. Conversely, traditional reciprocal-space implementations would in principle have required a single DFPT run for each individual value of q. In practice, this is often circumvented in reciprocalspace implementations, since efficient and accurate interpolation schemes for vibrational frequencies exist [86]. For the exact same reasons, finite-difference strategies can yield accurate results even in very limited supercells [81,82]. However, this is no longer the case if more complex response properties such as the electron-phonon coupling [33,52] need to be assessed. In that case, reciprocal-space formalisms either need to sample the Brillouin zone by brute-force [33] or to rely on approximate interpolation strategies, e.g., using a Wannierization of the interactions in realspace [52]. The approach discussed in this work allows to overcome these limitations and to consistently assess all these properties using the well-controlled wavefunction expansion already used in the ground-state DFT and thus potentially lays the foundation for future research directions in this field. This is further substantiated by the scaling behavior discussed in the previous section. Despite being a proof-of-concept implementation that has not undergone extensive numerical optimization, we find the code to exhibit quite favorable scaling properties and a promising performance that can be even improved further. For instance, the exploitation of space and point group symmetry would straightforwardly lead to significant savings in computational time, especially for high-symmetry periodic systems. Along these lines, symmetry can also be used to optimize the construction of the supercell used in the DFPT calculations. For the sake of simplicity, this procedure so far relies on translated images of the complete and intact unit cell. For particularly large and/or oblique unit cells this can result in a significant computational overhead, since the supercell can contain periodic images of atoms that do not interact with the unit cell at all. Accordingly, optimizing the supercell construction procedure can immediately lead to computational savings without loss of accuracy. Following these strategies, linear scaling should be achievable [87] for large system sizes (hundred and more atoms per unit cell). This would facilitate DFPT calculations of vibrational properties and of the electron-phonon coupling for fully converged q-grids in complex systems, such as organic molecules adsorbed on surfaces. For such kind of applications, additional computational savings can be gained in our proposed real-space approach by artificially restricting the calculation to the actual degrees of freedom of interest, e.g., the ones of the absorbed molecule.
The formalism described in this paper could also be extended to all type of perturbations, e.g. homogeneous electric field perturbations, in this case only one perturbation per cartesian direction needs to be considered regardless of the system size.

Acknowledgments
H.S. acknowledges Wanzhen Liang and Xinguo Ren for inspiring discussions. We further acknowledge Volker Blum for his continued support during this project. The project received funding from the Einstein foundation (project ETERNAL) and the European Union's Horizon 2020 research and innovation program under grant agreement no. 676580 with The Novel Materials Discovery (NOMAD) Laboratory, a European Center of Excellence. P.R. acknowledges financial support from the Academy of Finland through its Centres of Excellence Program (Project No. 251748 and 284621).

Appendix A. Convergence behavior of forces with respect to the degree of self-consistency
To investigate to which extent the last term of Eq. (8) really vanishes in practice, we have chosen Si (diamond structure) and Al (fcc) as examples. In both cases, one atom was displaced by 0.1 Å, which results in forces on this atom in the order of 10 0 eV/Å and 10 −1 eV/Å, respectively. To investigate what happens in calculations, in which full self-consistency has not yet been reached, we have then run a series of calculations with different break conditions for the self-consistency cycle. We only used the maximally allowed change in charge density as break condition and varied its value between 10 −2 and 10 −8 electrons. For the last setting, full self-consistency is achieved: Indeed, the observed change in energy/eigenvalues in the last iteration of such fully converged calculations is 10 −11 eV/10 −7 eV for Si, and 10 −12 eV/10 −7 eV for Al. In Fig. A  charge density < 10 6 ), the error in the force due to the non-fullyachieved self-consistency is thus typically smaller than 1 meV/Å. In this context, it is however important to note that in relaxation or MD calculations FHI-aims requires to specify a self-consistency break condition also for the maximum change in the forces, so that in practice these errors are well-controlled.

Appendix B. First order density matrix
The sum over states in the first order density matrix can be divided into sums over occupied-occupied states, occupied-unoccupied states, and unoccupied-occupied states: µm,i C This means that we only need to calculate the occupiedunoccupied sum for U (1) ai in the CPSCF equation, Eq. (95), while the occupied-occupied part is computed from the first order overlap matrix.

Appendix C. First order energy weighted density matrix
Similar like the first order density matrix, by using Eqs. (93) and (94), we can rewrite Eq. (99) into an occupied-occupied part, occupied-unoccupied part, and an unoccupied-occupied part:

Appendix D. 32 molecules' frequencies
Tables D.4-D.9 list the vibrational frequencies obtained via DFPT and via finite-differences for systems containing two, three, four, five, six, and eight atoms, respectively. In these comparisons, the atomic displacement is set to 0.013 Å in finite-difference