All-Electron, Real-Space Perturbation Theory for Homogeneous Electric Fields: Theory, Implementation, and Application within DFT

Within density-functional theory, perturbation theory~(PT) is the state-of-the-art formalism for assessing the response to homogeneous electric fields and the associated material properties, e.g., polarizabilities, dielectric constants, and Raman intensities. Here we derive a real-space formulation of PT and present an implementation within the all-electron, numeric atom-centered orbitals electronic structure code FHI-aims that allows for massively-parallel calculations. As demonstrated by extensive validation, this allows the rapid computation of accurate response properties of molecules and solids. As an application showcase, we present harmonic and anharmonic Raman spectra, the latter obtained by combining hundreds of thousands of PT calculations with \textit{ab initio} molecular dynamics. By using the PBE exchange-correlation functional with many-body van der Waals corrections, we obtain spectra in good agreement with experiment especially with respect to lineshapes for the isolated paracetamol molecule and two polymorphs of the paracetamol crystal.


Introduction
The response of molecules and solids to an applied electric field is a fundamental physical mechanism of prime importance, since it determines significant properties and spectroscopic signals, such as dielectric constants, Raman spectra, and sumfrequency generation spectra.In first-principles frameworks, these quantities are typically computed via time-dependent density-functional theory [1,2] or via analytical perturbation theory (PT) in either its density-functional perturbation theory (DFPT) [3,4,5,6] or coupled perturbed self-consistent field (CPSCF) formulation [7,8,9,10,11,12].Within these linear-response approaches, an additional complexity arises for the treatment of solids: As discussed in more detail in Sec. 3, the position operator appearing in the respective equations is not well defined and the expressions need to be recast into a more suitable form.Practical implementations of these methods within Kohn-Sham density-functional theory (DFT) differ substantially by their choice of basis sets (e.g., plane-waves or localized basis sets) and by their treatment of the core electrons (e.g., all-electron or pseudopotentials).In this paper, we address an implementation of PT for the response to a homogeneous electric field targeted towards handling large periodic systems, which will be the subject of our showcase in Sec. 6.Our implementation uses the all-electron, numeric atom-centered orbital based framework of the FHI-aims code [13,14,15], which also features implementations of PT for vibrational [16] and magnetic properties [17].Notably, this infrastructure allows us to treat isolated systems (such as molecules) and extended systems (such as periodic, crystalline solids) on the same footing, as discussed in Sec. 2 and 3, respectively.
As an application, we focus on the important task of calculating anharmonic vibrational Raman spectra of molecular crystals.These spectra are able to provide information about differences in the polymorphic structure of these crystals, the presence of impurities, and the onset of phase transitions.Importantly, they are quantities that can be readily accessed experimentally under different thermodynamic conditions, which can also be unambiguously simulated.In that respect, the calculation of these spectra in an anharmonic fashion using time-correlation functions [18,19], as further detailed in Sec.6, represents an important link between computer simulations and experiments.It can help to gauge the impact of anharmonicities in different phonon modes, which opens the path for a better understanding and control of the polymorphic forms of molecular crystals.The particular character of our implementation results in an efficient scaling with respect to system size (due to a sparse representation of the density matrices) and efficient numerical scaling with respect to the number of cores used on modern, massively parallel architectures (due to the use of local real-space operations).This facilitates the calculations of tens of thousands of polarizability tensors along ab initio molecular dynamics (MD) trajectories and thus enables the evaluation of anharmonic Raman spectra.We discuss how these spectra depend on different functionals and propose ways to obtain them at minimal cost.Our ab initio spectra computed at room temperature compare very well with experimental data obtained at the same conditions.The remainder of this paper is organized as follows.The fundamental perturbation theory framework is discussed for finite and extended periodic systems in Sec. 2 and 3, respectively.In Sec. 4, a detailed derivation of the respective equations and their implementation in a real-space, all-electron, numeric atom-centered orbitals based framework is presented.In Sec. 5, our approach and implementation is validated by comparing the calculated analytical polarizabilities and dielectric constants to literature values or to ones computed via finite-differences.Furthermore, we discuss the convergence behavior of our implementation, the scaling with system size, and the parallel performance when a large number of cores is used.We finish in Sec.6 by applying the developed formalism to compute harmonic and anharmonic Raman spectra for different polymorphs of the paracetamol crystal.

Fundamental Theoretical Framework
Before addressing the implementation in the FHI-aims code, we recall the basic equations used in this work.Throughout the text, we use a spin-unpolarized notation for the sake of simplicity, but a formal generalization to a collinear (scalar) spin treatment is straightforward.Moreover, we focus on systems with a non-vanishing energy gap for electronic excitations, because electric fields in metals are fully screened (Our numerical strategy to deal with quasi-degenerate electronic states in non-metallic systems is discussed in Sec.4.3.).In this section, a detailed derivation of the equations for finite, molecular systems is given; a generalization to extended periodic solids follows in Sec. 3.
In Kohn-Sham DFT, the total-energy functional is given by Here, n(r) is the electron density, T s the kinetic energy of non-interacting electrons, E ext the external energy due to the electron-nuclear attraction, E H the Hartree energy, E xc the exchange-correlation energy, and E nuc-nuc the repulsion energy of the nuclei.The ground state electron density n 0 (r) (and the associated ground-state total energy) is obtained by variationally minimizing Eq. ( 1) under the constraint that the number of electrons N e is conserved.This yields the chemical potential µ = δE KS /δn of the electrons and the Kohn-Sham single particle equations ĥKS for the Kohn-Sham Hamiltonian ĥKS .In Eq. ( 2), ts denotes the kinetic energy operator, vext the external potential, vH the Hartree potential, and vxc the exchange-correlation potential.Solving Eq. ( 2) yields the Kohn-Sham single particle states ψ p and their eigenenergies p .For a spin-unpolarized system, these states determine the electron density via whereby the occupation numbers f ( p ) are chosen in such a way that the N e /2 states with the lowest eigenvalues p are doubly occupied.
To solve Eq. ( 2) in numerical implementations, the Kohn-Sham states are expanded in a finite basis set χ µ (r − R I(µ) ) as with the expansion coefficients C µp .The chosen notation highlights that in a numerical atom-centered basis set, each basis function µ is associated to an atom I(µ) situated at R I(µ) .In such a basis set, Eq. ( 2) becomes a generalized 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 becomes a minimization with respect to the expansion coefficients C νp in which the eigenstates ψ p are constrained to be orthonormal.Typically, the ground state density n 0 (r) and the associated total energy E tot are determined by solving Eq. ( 6) iteratively, until self-consistency is achieved.
If an external electric field E is applied to an isolated system, the KS Hamiltonian gains an additional term ĥE = −r • E. If this electric field E = (e x , e y , e z ) is a superposition of homogeneous electrical fields with strengths e γ aligned along the different cartesian axes γ, the additional term ĥE contributes to the total energy functional in Eq. (1).A perturbative Taylor-expansion of the total energy in the zero-field limit then gives where η, γ are Cartesian directions.For isolated systems, the coefficient in the linear term which corresponds to the γ-component of the dipole moment of the system in its ground state, can be directly evaluated at the DFT level of theory due to the Hellmann-Feynman theorem.However, this is not possible for the coefficient in the second-order term, i.e., the polarizability since the derivative (or response) of the ground-state density with respect to the field strength is explicitly required.More generally, this is formalized within the 2n + 1 rule [20], which states that second-order derivatives of the total energy [21,22,23] cannot be direcly calculated from the ground-state electron density or wavefunction alone, but also require the respective first-order derivatives of the electron density or wavefunction, i.e., their linear response to the perturbation.We use perturbation theory (PT) to obtain the required derivatives.In this formalism, the response to perturbations along different Cartesian axes γ can be treated independently viz.subsequently, so that the short-hand notation used in the following for ground-state M (0) and response properties M (1) is always welldefined, e.g., In this way, we can express the linear Taylor-expansion of the Kohn-Sham Hamiltonian in the limit of vanishing field along the γ-axis as: ĥKS (e γ ) ≈ ĥ(0) KS + ĥ( 1) where the response of the Hamiltonian operator is ĥ(1) Introducing the analogous expansions for the single-particle states ψ p (e γ ) and their eigenvalues p (e γ ), rearranging the linearorder terms in the KS equation ĥKS (e γ )ψ p (e γ ) = p (e γ )ψ p (e γ ), and applying the normalization condition ψ p (e γ )|ψ p (e γ ) = 1, yields the Sternheimer equation ĥ as well as the condition By multiplying Eq. ( 16) with ψ (0) q | from the left one obtains To solve this equation numerically, we expand the response of the wave functions in terms of the unperturbed states ψ (0) q (r).Here, we chose U (1) pp = 0 for all p to fulfill Eq. ( 17) and hence obtain an algebraic expression for Eq. ( 18) The expansion using the matrix U qp employed in this work is typical for the coupled perturbed self-consistent field (CPSCF) formulation [7,8,9,10,11,12] of PT.For our implementation described in Sec. 4, such an expansion in terms of orbitals is advantageous, since it allows leveraging the already existing algorithms for the massively-parallel evaluation of matrix elements in this representation [13,14,16].Accordingly, the matrix elements H (1) KS |χ ν are defined as for unperturbed calculations, i.e., using the numeric atomic orbitals introduced in Eq. ( 4).This allows us to directly compute the non-diagonal elements (q = p) of The matrix U qp , which fulfills U , plays a central role in our implementation: As discussed in detail in Sec. 4, it allows us to directly determine the response of the density n (1) = p f ( p ) ψ (1)  p ψ (0) p + ψ (0) p ψ (1)   p (22) in a density-matrix formalism, i.e., without explicitly computing the response of the eigenvalues p , of the wave function ψ p (r), or its coefficients C µp , which is computationally advantageous.Using U (1) qp , one can then directly evaluate the polarizability tensor α γδ defined in Eq. (10) in finite, isolated systems.Let us note that implementations of PT in plane wave codes typically do not use the expansion in terms of orbitals defined in Eq. ( 19) via the U (1) qp matrix, but rather compute the coefficients C (1) µp by directly solving Eq. ( 18) in the space spanned by the KS states using the density-functional perturbation theory formalism (DFPT) [3,4,5,6].In such codes, in which thousands of orbitals, i.e., plane waves, need to be considered, the DFPT approach is advantageous.

Generalization to Periodic Solids
For periodic boundary conditions (PBCs), the main physical reasoning behind the derivation of Eqs.(1)-( 22) still remains valid.However, three specific adaptations have to be made: First, the basis set expansion introduced in Eq. ( 4) is slightly different, as described in detail in Refs.[13,24,16]: The periodic images of the nuclei R Im = R I + R m are accounted for by summing over the lattice vectors R m , i.e., over linear combinations of the unit cell lattice vectors a 1 , a 2 , a 3 .Analogously, also the numeric atomic orbitals associated with such periodic images, e.g., χ µm (r) = χ µ (r − R I(µ) − R m ) associated with the periodic image m of nucleus R I , gain an additional index m that describes their relative position to the unit-cell equivalent.To account for translational symmetry and exploit Bloch's theorem, Bloch-like generalized basis functions are constructed from the local atomic orbitals and then used in the basis set expansion Accordingly, all relevant physical quantities such as the KS Hamiltonian gain an additional dependence on the wavevector k, so that Eq. ( 5) becomes ν Therefore, the summations over electronic states appearing in Eqs. ( 1)-( 22) now feature an additional analytical integration over the Brillouin zone that is approximated numerically by a sum over a finite k-grid with N k points.Similarly, the real-space integrals in Eqs. ( 1)-( 22) are no longer indefinite, but definite and limited to the unit cell (u.c.), as it is the case in Eq. (25).Second, it is necessary to consider the screened electric field E = D − 4πP in the solid, where D is the electric displacement [25] and the polarization in the unit cell volume V is given by [3,6]: The relationship between the components of the electric displacement and the screened field defines the high-frequency dielectric constant [25] ε where η, γ are Cartesian directions and δ is the Kronecker delta symbol.For a screened field E = (e x , e y , e z ) that consists of a superposition of homogeneous electrical fields with strengths e γ aligned along the different cartesian axes γ, one can follow the derivation given in the previous section to obtain A comparison with Eqs. ( 9) and ( 10) reveals the formal relationship between the dipoles µ γ and the polarizabilities α γδ discussed in the previous section for molecules and the polarization P γ , i.e., a dipole density [26], and its derivative with respect to the screened field in solids.For the sake of notational clarity, the "molecular" notation with µ γ and α γδ is used for the remainder of this paper.Third, complications arise due to the fact that the superposition of homogeneous electric fields E is not periodic, as alluded to in the introduction.As a consequence, the definite integral over the unit cell required to determine the Hamiltonian response H (1) µν (k) is ill-defined in PBCs, since ĥ(1) KS given in Eq. ( 14) contains the position operator r γ , which is itself ill-defined in this case.The same problem affects Eq. ( 30).In reciprocal space implementations, the Berry-phase formalism [27,28,5,21,4,23] is typically the method of choice; a tutorial introduction to this approach can be found in Ref. [29].In real space implementations, the position operator can be rewritten in a boundary-insensitive form [22] by exploiting the properties of the commutator between the KS-Hamiltonian and the position operator ĥ(0) KS (k), r = −∇.With that, one gets the well-known expression that can be used to evaluate the non-diagonal matrix elements (q = p) Using Eqs. ( 23) and ( 24) we obtain the representation with With that we can recast the expectation value H (1) p (k) appearing in Eq. ( 18): Here, the matrix elements V (1) xc |ϕ ν (k) can be directly evaluated as done in Eq. ( 25) since they only depend on lattice periodic operators.Now, the matrix U qp (k) introduced in Eq. ( 21) is computed as Similarly, the polarizability tensor components appearing in Eq. ( 10) can be rewritten as using the matrix elements defined in Eq. ( 34).As explicitly highlighted in the notation, the matrix U (1,δ) qp (k) associated with a perturbation along the Cartesian axis δ has to be used in this case, whereas the matrix Ω (γ) qp (k) is associated with a perturbation along the Cartesian axis γ.Throughout the remainder of this work, the more general formulation in terms of Bloch-functions ϕ µ (k) and wave vectors k is used, since a simplification to finite systems is straightforward.

Details of the Implementation
Our implementation closely follows the flowchart shown in Fig. 1: After a ground state DFT calculation (see Ref. [13]) is completed, the matrix pq (k) = 0 is used as initial guess, one obtains in the first iteration, which can then be fed back to the self-consistency loop to determine the first-order density response n (1) (r) in a density matrix formalism (see Sec. 4.1).As detailed in Sec.4.2, we then use n (1) (r) to compute the remaining, individual ingredients that enter ψ p (k) , i.e., the matrix elements V µν (k) defined in Eq. (38).The Sternheimer equation then provides a new matrix U (1) qp (k), as discussed in Sec.4.3.We iteratively restart the PT loop using a Pulay-mixer [30], until selfconsistency is reached, i.e., until the changes in the response of the density-matrix P (1)  become smaller than a user-given threshold.In the last step, the polarizability and the dielectric constant are computed, as discussed in Sec.4.4.Atomic units are used in the complete workflow.
Both the ground-state density n (0) (r) and the response of the density n (1) (r) are periodic, i.e., invariant against translations n (1) (r) by a lattice vector R m , as shown in Fig. 2. Accordingly, we can use the algorithms used in ground-state calculations and discussed in detail in Refs.[13,24] for many aspects of our implementation.In the following, we thus mainly focus on the practical details that are specifically needed for the computation of the response to a homogeneous electric field.

Response of the Electronic Density
To numerically compute the electronic density n(r) in ground-state calculations [13], we use a density matrix formalism which is obtained by inserting Eqs. ( 24) and ( 23) into Eq.(3).Hence, the density matrix is given by Here, the chosen notation using the index o highlights that the sum over all states only needs to be performed over occupied states with f ( o (k)) = 0 in practice.Similarly, the  response of the electronic density can thus be expressed as n (1) (r) = µm,νn using the response of the density matrix given by In the sum over states o, we express C (1) (k) in terms of U (1) (k) via whereby we split the sum over q into two separate sums over o and u.In practice, these two sums can then be later evaluated by restricting the sum over o to occupied and the sum over u to unoccupied states, respectively.Accordingly, also the sum over o appearing in Eq. ( 46) can be split into two double sums, one over o, o and one over o, u.
For the first one, we obtain , cf.Eq. (39).For the second double sum, we obtain by switching the summation indices u, o in the second term and using U (1) , as done already for Eq.(48).By this means, the response of the density matrix can be written as P (1)  µm,νn = 1 In practice, the evaluation of Eq. ( 51) can thus be restricted to the double sum over occupied o and unoccupied states u.

Response of the Kohn-Sham Hamiltonian
As discussed in Sec. 3 for Eq. ( 37), the computation of ψ µν (k) and Ω (0) µν (k), which are defined in Eqs. ( 34)- (35) and which are required to calculate ψ p (k) , are computed before the self-consistency loop, since they only depend on unperturbed properties.The definite unit-cell integral appearing in Eq. ( 35) is integrated on a real-space grid using the formalisms described in Refs.[13,14].Conversely, the matrix V (1) µν (k), which is defined in Eq. (38) and which is required to compute ψ (0) p (k) , explicitly depends on the response of the density n (1) (r) and thus needs to be updated each cycle.For that purpose, we first compute its ingredients on a real-space grid, i.e., the response of the electrostatic potentials v(1) ext (r) and v(1) H (r) as well as the response of the exchange-correlation potential v(1) xc (r), as discussed below.The matrix elements V (1) µν (k) are then again obtained by performing the real-space unit-cell integral appearing in Eq. ( 38) with the aforementioned techniques.

Response of the Electrostatic Potentials
As discussed in detail in Refs.[13,24,16], the electrostatic potential generated by the nuclei and the electrons is computed in FHI-aims ground-state calculations using a scheme proposed by Delley [31]: The ground-state density n (0) (r) is decomposed into two terms The first term describes the density associated with a superposition of "free", i.e., completely isolated, spherically symmetric atoms n free Im (r) located at the positions of the nuclei and of their periodic images R Im .The potentials of n free Im (r) and δn(r) are computed independently and then reassembled to get the full electrostatic potential that enters the Kohn-Sham equations.For this purpose, δn(r) is further decomposed into atom-specific multipoles, the contributions of which are added up in an Ewald-like summation to account for long-range interactions, cf.Refs.[13,24,31].Given that the density response n (1) (r) is also periodic in the perturbed case, see Fig. 2, we can use the exact same formalism to obtain the electrostatic potential associated with it.There is only one small difference: In this case, the "free", spherically symmetric atoms do not contribute to the associated electrostatic potential at all.

Response of the Exchange-Correlation Potential
In semi-local approximations, the exchange-correlation potential vxc (r) entering the Kohn-Sham Hamiltonian in Eq. ( 2) is given by Accordingly, its response v( 1) xc (r) can be obtained via v( 1) by integrating over the the exchange-correlation kernel f xc (r, r ), i.e., the second functional derivative of the exchange-correlation energy E xc [n(r)], and the density response n (1) (r ).For the local-density approximation (LDA) [32,33] and the PBE functional [34,35] in the generalized-gradient approximation (GGA), we have implemented the standard expressions for f xc (r, r ).Additionally, many more exchangecorrelation kernels are accessible in our implementation via the Libxc library [36].
For isolated systems, we have also implemented the response of the exact-exchange potential.For Hartree-Fock and hybrid functionals, an additional exchange term needs to be added to the entries H µ,ν of the Hamiltonian response matrix.Here, (χ µ χ λ |χ ν χ σ ) is the two-electron, four-index Coulomb integral defined and discussed in Refs.[15,37,38] and P (1) is the first order density matrix defined in Eq. (51).

Stable Evaluation of the Expansion Matrix
To compute U (1) (k), one can in principle just evaluate Eq. (39) as discussed in the beginning of Sec. 4. Thereby, only the entries associated to unoccupied-occupied (uo) orbital pairs need to be computed, since these are the only entries that enter the response of the density matrix P (1) , as shown and discussed for Eq. ( 51).Obviously, Eq. ( 56) becomes numerically unstable when quasi-degenerate eigenvalues are present close to the Fermi energy F , since the denominator u (k) approaches zero in that case.In order to overcome this difficulty, we employ the technique originally proposed by de Gironcoli [39,6] for DFPTbased lattice dynamics calculations in metals.For this purpose, we use a Fermi function with a small smearing width σ for the occupation numbers f ( o ) and f ( u ) appearing as difference in Eq. ( 51).We then pull this difference f ( o ) − f ( u ) inside the evaluation of U uo and re-write the problematic prefactor in Eq. ( 56) as as detailed in Ref. [39].This has virtually no influence in the regime o (k) − u (k) > σ.
For o (k) − u (k) σ, we replace and evaluate the rewritten problematic factor by its analytic limit for u → o : This expression is always finite and therefore numerically stable, even in the case of vanishingly small energy differences.

Evaluation of Polarizabilities
In the last step, we evaluate the polarizability by rewriting Eq. ( 40): In the first step, the use of Ω qp (k) = Ω * pq (k) reduces the summands to a real part extraction Re(), while in the second step the same procedure as used to obtain Eq. ( 48) is applied.By this means, the double sum can be limited in practice to only run over unoccupied u and occupied o states.The same strategy introduced with Eq. ( 58) and discussed in the previous section can be applied to deal with quasi-degenerate eigenvalues here.Again, the matrix U (1,δ) (k) appearing in Eq. ( 62) is associated with a perturbation along the Cartesian axis δ, while the Ω (γ) (k) matrix is associated to a perturbation along the axis γ.

Validation and Performance
To validate our implementation we show how our simulations converge with respect to the numerical parameters used in the calculation in Sec.5.1.Furthermore, we compare our PT polarizabilities to those obtained from finite differences in Sec.5.2.These tests are then extended to periodic systems in Sec.5.3.The computational performance of the implementation is discussed in Sec.5.4.

Convergence with respect to Basis Set Size and k-Point Grid Density
We observe that our calculated polarizability tensors are most sensitive to the basis set size and the amount of k-points used in the simulation, as shown below.All other numerical parameters either influence the results very little or show a similar convergence behavior as in ground-state DFT calculations.
First, we discuss the convergence of polarizabilities in our implementation with respect to the basis set size used for the expansion of the Kohn-Sham states in Eq. (4).As an example, we use ethylene (C 2 H 4 ), for which we compute the three diagonal components α γγ of the polarizability tensor using LDA [32,33].In all cases, the PT calculations were performed for the same geometry, i.e., the structure obtained by geometry optimization (maximum residual force < 10 −4 eV/ Å) with tight basis sets and numerical settings.The C-C bond of the molecule is oriented along the Y-axis.
Figure 3 shows the absolute error in the diagonal components of the polarizability tensor with increasing basis-set size.Here, a minimal basis includes exactly one basis function per electron; additional functions are then added in groups, so-called tier 1, tier 2, etc., basis sets (see Ref. [13] for more details).The polarizabilities converge slowly with the basis set size in finite molecular systems as ethylene: Although getting qualitatively correct results, the maximum absolute (relative) error is for instance still 2.44 a 3 0 (11 %) at a tier 2 level.Only at the tier 3 level we get a maximum absolute (relative) error of 0.23 a 3 0 (1 %).For semi-infinite systems, the dielectric constant, which is directly proportional to the polarizability as noted in Eq. 28, converges much faster with increasing basis set size, as also shown in Fig. 3 for bulk silicon.Even at a tier 1 level we essentially achieve convergence with an absolute (relative) error of 0.007 (0.05 %).
The slower convergence observed for molecular systems arises from the inhomogeneous distribution of the localized basis sets in isolated systems.The standard basis sets in FHI-aims have been optimized to obtain converged ground state energies, but are not necessarily even-tempered for the calculation of polarizabilities, which can create an imbalance in the extent of the polarization that is possible in different directions.One possibility to improve convergence would be the construction of basis sets that are specifically tailored for the calculation of polarizabilities, see for instance Ref. [40] for an example of basis sets adapted for polarizabilities, Ref. [41] for hyperpolarizabilities, or Ref. [17] for magnetic response properties.Alternatively, it is possible to include extra basis functions in otherwise empty regions to span the space much more efficiently.As shown below, this allows to reduce the computational cost by using much smaller overall basis sets without sacrificing accuracy.The difficult task in this procedure is to determine in which region of space the original basis sets are not sufficient, in order to determine where to best place the extra basis functions.In general, the symmetries of the molecule are helpful in this task and thus need to be considered as well.We illustrate this procedure for the polarizabilities of the C 2 H 4 molecule with LDA (see Table 1).It is clear that the addition of 2 carbon-like ghost atoms (i.e.only the tier 1 basis set of a carbon atom), which we positioned below and above the molecular plane on the bisection of the C-C segment (see Fig 4), significantly improves the convergence, almost to the level of tier 2, but at only half the computing time.Please note that simply increasing the onset of the cutoff potential for the usual basis sets in FHI-aims does not improve the performance of our results.Finally, to study the sensitivity of the polarizability tensor on the k-point grid density in periodic systems, we also use silicon as example.Fig. 5 displays the convergence behavior with respect to the size of the reciprocal-space k-mesh in the primitive Brillouin zone.We observe a maximum absolute (relative) error of Si, tier 2 Si, tier 3

Validation against Finite Differences
To validate our PT implementation, we also compared the obtained polarizabilities of 32 selected molecules to the ones obtained via finite difference calculations, as detailed in the Appendix.There, the details for each individual molecule can be found; here, this data is succinctly summarized in Tab. 2, where we list the mean absolute percentage error (MAPE) and the mean absolute error (MAE) for all tested molecules.Overall, we find an excellent agreement between our implementation and the finite-difference results.Mean absolute error (MAE) and mean absolute percentage error (MAPE) for the difference between the polarizabilities obtained via PT and finite differences (FD) for a set of 16 dimers, 5 trimers, and 11 molecules.All calculations are performed at the LDA level of theory with fully converged numerical settings and relaxed geometries.Detailed informations including the values for each individual molecule can be found in the Appendix.

Extended Systems: High-Frequency Dielectric Constant
In order to validate our implementation for extended systems, we have calculated the dielectric constant of several semiconductors using the local-density approximation (LDA [42]) and the generalized gradient approximation (GGA-PBE [34,35]) and compared it with experimental and theoretical data compiled from literature [22,1], see Tab. 3.All calculations have been performed at the theoretical equilibrium lattice constant using 16×16×16 k-points in the primitive unit cell and "tight" basis set and integration settings.Also, we list LDA/GGA literature results obtained using a plane wave basis set and norm-conserving pseudopotentials (NCPP) [22,1] or the projector augmented wave method (PAW [43,44]).With respect to experiment, we note that all LDA and GGA calculations overestimate the electronic dielectric constant by roughly 10% due to the well-known fact that these functionals yield too small band gaps [45,22].
With respect to theoretical results, the most recent literature data computed with the PAW method (LDA [43]; PBE [44]) is in excellent agreement with our implementation.Slightly larger deviations are observed with respect to earlier calculations that rely on norm-conserving pseudopotentials (NCPP): The agreement is generally better with literature results obtained using non-linear core corrections [46,47].For instance, this can be observed for GaSb: The work of Dal Corso et al. [1] made use of non-linear core corrections, but not the earlier one of Giannozzi et al. [22].For the latter work, the use of a smaller k-point grid may also be partially responsible for the observed deviations.
Exp. [48,49,50,51,52] this ) n H molecules (left) and periodic diamond (right) as a function of the number of atoms (diamond: in the unit cell) on 32 cores (see text).Following the flowchart in Fig. 1, also the timings required for the computation of the individual responses, i.e., the ones of the density n (1) (r), of the Hamiltonian matrix H (1) (k), and of the density matrix P (1) , are given.
To demonstrate the performance and scaling of our implementation, we show timings for the H(C 2 H 4 ) n H molecules with variable n = 8 − 128 and diamond.In the latter case, different supercell sizes were considered by increasing the number of building units in the unit cell from (C 2 ) 32 to (C 2 ) 512 .All calculations use light settings and the LDA functional.Only the Γ point is considered in the periodic case.Calculations were performed on a single node featuring two Intel Xeon E5-2698v3 CPUs (32 cores) and 4 Gb of RAM per core.
For the timings shown in Fig. 6 (molecules), we find that the integration of the Hamiltonian response matrix H (1) (k) determines the computational time for small 1.4 1.5 P (1)  2.5 2.6 Total 1.3 1.4 and the periodic diamond discussed in the text.The fits were performed using the expression t = cN α for the CPU time as function of the number of atoms N .
system sizes, i.e., for less than 200 atoms.Like for the update of the response density n (1) , which involves similar numerical operations, we find a scaling of nearly O(N ) for this step (see Tab. 4), as it is the case in ground-state DFT calculations [13].For very large system sizes (N 1000), the update of the response density matrix P (1)  becomes dominant, since it scales with O(N 2.5 ) in this regime.As discussed in Sec. 4, the computation of P (1) requires matrix multiplication operations, which traditionally scale O(N 3 ).For bulk diamond we find a similar behavior and fit similar exponents, as shown in Fig. 6 and Tab. 4 as well.We note that the prefactors to these timings are higher for dense 3D systems than for 1D systems and that they also are system dependent.Our real-space PT implementation thus exhibits a similar scaling as the underlying DFT calculations, as it is generally the case for DFPT/CPSCF codes.
In summary, we find an overall scaling behavior that is always smaller than O(N 2 ) for the investigated system sizes both in the molecular and the periodic case.Note that a parallelization over cores is already part of the presented implementation, given that the discussed real-space formalism closely follows the strategies used for the parallelization of ground-state DFT calculations in FHI-aims [13,14].As shown in Fig. 7 for a unit cell of the paracetamol crystal containing 160 atoms, almost ideal scaling is achieved for the time per PT iteration when different number of CPU-cores (same specifications as in the previous tests) are used.Still, it is very gratifying to see that even quite extended systems with more than 100 atoms in the unit cell are in principle treatable within the relatively moderate CPU and memory resources offered by a single state-of-the-art workstation.

Application: Harmonic and Anharmonic Raman Spectra
In order to showcase the usefulness and efficiency of our implementation, we calculate the non-resonant Raman spectra of paracetamol in its molecular form, as well as in its first (monoclinic) and second (orthorhombic) crystalline polymorph.More specifically, we wish to investigate the impact of anharmonicities on these spectra, due to their acknowledged importance in H-bonded, flexible systems [53,18,54].Focusing on molecular crystals, which often exist in multiple competing polymorphs with very different physicochemical properties, makes it necessary to have an accurate and efficient (r) Total (DFPT) Ideal scaling n (1) (r) P (1)  Total time Ideal scaling model to characterize such structures.Taking into account anharmonicities in the computation of Raman spectra is of crucial importance, as has already been proven in the past, for example for the characterization of phase transitions in high-pressure ice [18].
Vibrational Raman spectra are typically computed in the harmonic approximation, where the Raman intensities are proportional to the derivatives of the polarizability with respect to atomic displacements, as detailed, for example, in Refs.[55,56].In this work, we calculate these harmonic Raman intensities through finite-differences by numerically computing the derivatives of the polarizability tensor via finite displacements of the nuclear coordinates.These displacements are performed in the unit cell, since only phonons at the Γ point of the lattice contribute to the Raman intensity.Additionally, we also compute anharmonic Raman spectra through the calculation of polarizability autocorrelation functions in thermodynamic equilibrium.We simulate the nuclear dynamics using ab initio molecular dynamics and compute the polarizabilities along these trajectories via PT.As explained in, e.g., Ref [57], the polarizability tensor α can be divided into an isotropic ᾱ, and an anisotropic component β, As shown in Refs.[18,19] where pioneering simulations of this type were presented, the (non-resonant) Raman intensity I(ω) is then expressed as a sum of isotropic and anisotropic parts as Here, N is again the number of atoms in the system.Furthermore, since the autocorrelation functions • are computed classically, a quantum correction factor is usually applied.Due to the fact that the classical correlation function better approximates the Kubo transform of the quantum autocorrelation function, we multiply I(ω) in Eq. 64 by β ω/(1 − e −β ω ), where β = 1/k B T [58].Further frequencydependent factors that multiply the vibrational Raman lineshapes are experimentdependent [59,60,61].Here, we normalized experimental and theoretical spectra by their areas for comparison.All MD trajectories used in this paper have been obtained using the PBE functional in combination with many-body van der Waals interactions [62] (PBE+MBD), which have been previously shown to play an important role for the accurate assessment of potential-energy surfaces (PES) and free energies of molecular crystals [19,53,63,64].
600 900 1200 1500 Wavenumber (cm -1 ) I(ω) (normalized) 3000 3200 3400 3600 3800 PBE+MBD/LDA (~1h) PBE+MBD/PBE (~3h) PBE+MBD/PBE0 (~4h) Figure 8. Harmonic Raman spectrum of the paracetamol molecule: The notation XX/YY denotes that the PES (energy and forces) were calculated with the XX functional, while the YY functional was used for the polarizabilities in the PT part.In this figure, the PES is always obtained at the PBE+MBD level, while different functionals are used for the polarizabilities.Tight settings and basis sets were used and the calculated (finite-difference) Raman intensities were convoluted with a Gaussian function of fixed width for better visualization.Computational time required for each simulation is also denoted.
We first analyze the sensitivity of the harmonic Raman intensities of the paracetamol molecule to the employed exchange-correlation functional.In Fig. 8, the PES is always treated at the PBE+MBD level, but different functionals are used to calculate the polarizabilities.We observe that the Raman spectra are essentially insensitive to the choice of xc-functional (LDA, PBE, and hybrid functional PBE0) used in the PT part.This is due to the fact that the magnitude of the Raman peaks are proportional to the polarizability derivatives with respect to atomic displacements and these derivatives are very similar in all functionals.Since evaluating the polarizabilities at the PBE0 level is four times more expensive than with LDA, we can decrease the cost of these simulations without sacrificing accuracy by evaluating the PT portion at the LDA level.Conversely, the xc-functional chosen for the assessment of the PES has a large impact on the position of the peaks.In Fig. 9, we highlight this fact by showing the harmonic Raman spectra of the paracetamol molecule obtained using the LDA xc-functional for the polarizabilities, but different xc-functionals for the PES (energy and forces including full geometry relaxation).Switching from LDA to PBE (or to PBE0) for probing the PES results in noticeable changes in the harmonic Raman spectrum, as can be seen from the shifts in the peak positions.We also note that the main differences introduced by Hartree-Fock exchange in the spectrum are rigid blue shifts of the peak positions, especially above 1000 cm −1 , which means that these vibrational modes become more stiff [65,66].In (a) we also show a spectrum obtained from thermostatted ring polymer molecular dynamics (TRPMD), which accounts for the quantum nature of the nuclei.In all cases, the PES was probed with the PBE+MBD functional, while the polarizabilities were calculated with the LDA functional.Harmonic Raman intensities were convoluted with Gaussian functions for better visualization.

Wavenumber (cm
In Fig. 10 (a) and (b), we show our calculated harmonic and anharmonic Raman spectra for the isolated paracetamol molecule and the paracetamol crystal in its monoclinic form I. For the molecule, we show anharmonic spectra obtained from an ab initio MD trajectory at 300 K with classical nuclei and also spectra obtained from a thermostatted ring polymer MD trajectory at 300 K, which accounts for quantumnuclear effects in the dynamics [67,68].We have used 16 replicas of the system and the Generalized-Langevin Equation thermostat proposed in Ref. [68] for the internal modes of the ring polymer.While the general shape of the harmonic and anharmonic spectra are similar (both for the molecule and the crystal), several peaks are shifted with respect to one another and feature different relative magnitudes, which substantiates the importance of anharmonic effects.In the molecular case, nuclear quantum effects induce a red-shift (with respect to the anharmonic classical spectrum) of about 70-100 cm −1 in the high-frequency range.The effect is much less pronounced in the lower frequency regions, as expected [66].In the crystal, the lineshapes of the harmonic and anharmonic approximations are quite different, which highlights the fact that in our anharmonic spectrum we are able to capture the Raman peak lifetimes, while in the harmonic approximation we are simply convoluting the Raman intensities with Gaussian functions of a fixed width (and not explicitly calculating lifetimes).For periodic and condensed phase systems, we have previously shown [69,70,68] that nuclear quantum effects would have a similar impact on the spectrum as for the molecular case.For water at room temperature for instance, the OH-stretch peaks are red-shifted by 150 cm −1 solely due to nuclear quantum effects [70,68].
In order to further evaluate the quality of our simulations we turn to a comparison to experimental data: Figure 11 shows our computed anharmonic Raman spectra for polymorphs I and II of the paracetamol crystal, respectively, compared to experimental spectra from literature.Both spectra were calculated from 2 independent MD runs of 15 picoseconds each.A time step of 0.5 femtosecond was used and the polarizability was computed every femtosecond.Our results show a very good agreement with experiment, especially in terms of lineshapes for both crystalline forms.As previously discussed in literature and above, the observed rigid shifts between experimental and theroretical spectra originate from the choice of functional and the lack of nuclear quantum effects in the simulations.Employing a higher-level hybrid functional can be estimated to lead to blue-shifts of up to 180 cm −1 for frequencies above 3000 cm −1 (see Fig. 9), while considering the quantum nature of the nuclei would red-shift these frequencies by up to 150 cm −1 at high frequencies, as discussed above.To some extent, these effects hence cancel each other and are less pronounced at lower frequencies, which explains the good agreement observed in the 600-1800 cm −1 region between calculated and experimental spectra.However, the calculated spectra are still blue-shifted with respect to experiment above 2500 cm −1 , even though line-shapes are well reproduced.The inclusion of nuclearquantum effects in the simulation would most likely solve this discrepancy, but the cost of such a simulation would be prohibitive at this point.
It is interesting to note that experimental spectra may sometimes differ slightly from one another as well.These differences are noticeable in the relative intensities of peaks or appearance/disappearance of low-intensity peaks [71,72,73,74].These differences reflect the difficulty to control the experimental setup for a wide range of frequencies and to synthesize a pure sample especially in the case of those polymorphic crystals, which undergo phase transitions under specific thermodynamic conditions.In particular, as explained in Ref. [72], the crystallization of paracetamol in form II is often not perfect, as some traces of metacetamol may remain present, leading to partially mixed Raman spectra.

Conclusions
In this paper, we derived and implemented a real-space formulation of perturbation theory for homogeneous electric fields within an all-electron, numeric atom-centered orbitals DFT framework.We validated the approach by computing polarizabilities (and dielectric constants) 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.Due to the slow convergence of polarizabilities with respect to the basis set size for isolated systems, we propose a simple solution based on the addition of so-called "ghost" atoms (i.e.only basis functions) in parts of space that are not densely populated.Also, we show how to stabilize our implementation for situations where small differences between occupied and unoccupied eigenvalues are present, arriving at a formulation which proved always stable.The scaling behavior of our implementation for the calculation of polarizabilities is between O(N ) and O(N 2 ) for both non-periodic (O(N 1.3 )) and periodic systems (O(N 1.4 )).In order to reduce the total time to O(N ), more advanced algorithms [11,12] for the evaluation of the density-matrix response P (1) could be pursued in the future.
We have tested our approach for the computation of dielectric constants by comparing theoretical and experimental literature data for a variety of semiconductors, obtaining very good agreement.To highlight the power of our PT implementation, we applied it to the calculation of anharmonic Raman spectra of the isolated molecule of paracetamol, as well as two of its polymorphic crystal forms, which involved the computation of hundreds of thousands of polarizability tensors in order to build the time series.We obtained good agreement with experiment in all cases especially for the lineshapes, which highlights the power of ab initio MD to capture anharmonic phonon frequencies and lifetimes, as well as the respective material properties [75].Regarding the calculated peak positions, we observe blue-shifts in the NH and CH stretching regions that stem from the lack of nuclear quantum effects in the MD simulations, as we explicitly show for the isolated molecule.We also found that these spectra are very sensitive to the xc-functional employed for the assessment of the potentialenergy surface, but that they are rather insensitive to the xc-functional employed for the calculation of the polarizabilities.In fact, we obtain correct spectra in a computationally efficient manner by using LDA for the polarizability tensors, but the PBE functional with many-body van der Waals corrections for the PES.We have shown that having such an efficient implementation that gives access to anharmonic Raman signals will be extremely useful for the analysis of experimental Raman spectra, which are often used to characterize new polymorphic forms of (molecular) crystals.
The data presented in this work as well as the input and output files used to produce it are publicly available as a dataset [76]

Figure 1 .
Figure 1.Flowchart for the calculation of the polarizability.Loops are performed over the different Cartesian coordinates and, in the case of periodic boundary conditions, over the finite k-grid.

Figure 2 .
Figure2.Ground-state electronic density n (0) (r) and its response n(1) (r) to an electric field, as exemplarily computed for an infinite, periodic H 2 chain.

Figure 3 .
Figure 3. Convergence behaviour of the polarizabilities α xx , α yy , α zz of ethylene and of the high-frequency dielectric constant ε ∞ xx of bulk silicon (16 × 16 × 16 k-points) with respect to the basis set size (see text).

Figure 4 .
Figure 4. Sketch of the C 2 H 4 molecule and its two ghost atoms used to improve the convergence.Ghost atoms are pictured at the top and bottom of the molecule from this perspective.

Figure 5 .
Figure 5. Convergence of the diagonal components of the high-frequency dielectric constant ε ∞ xx of bulk silicon with respect to the k-point density.The size of the k-grid is N k × N k × N k .Tight grid settings and tier 2 as well as tier 3 basis sets are used.The benchmark value is calculated using N k =24.
s.c.f.iteration (s) n

Figure 7 .
Figure 7. Scaling of the CPU time per PT cycle with the number of cores (parallel scalability) for the paracetamol crystal (form II) containing 160 atoms in the unit cell.Tier 1 basis sets and a 2×2×2 k-grid are used.The time required for the computation of the individual response properties is also shown.

Figure 9 .
Figure 9. Harmonic Raman spectrum of the paracetamol molecule: The notation XX/YY denotes that the PES (energy and forces) were calculated with the XX functional, while the YY functional was used for the polarizabilities in the PT part.In this figure, the polarizabilities are always obtained at the LDA level, while different functionals are used for the PES.Tight settings and basis sets were used and the calculated (finite-difference) Raman intensities were convoluted with a Gaussian function of fixed width for better visualization.

Figure 10 .
Figure 10.Comparison of harmonic and anharmonic (300 K) Raman spectra of (a) the paracetamol molecule, and (b) the paracetamol crystal (form I).In (a) we also show a spectrum obtained from thermostatted ring polymer molecular dynamics (TRPMD), which accounts for the quantum nature of the nuclei.In all cases, the PES was probed with the PBE+MBD functional, while the polarizabilities were calculated with the LDA functional.Harmonic Raman intensities were convoluted with Gaussian functions for better visualization.

Figure 11 .
Figure 11.Raman spectra of paracetamol-form I (top) and II (bottom) calculated at 300 K. Experiment from Ref.[71] at room temperature.The spectra have been normalized to one in each panel.

Table 1 .
Influence of H-and C-like "ghost" atoms on the diagonal elements of the polarizability of C 2 H 4 , using light settings and LDA.Numbers in brackets indicate the mean polarizability.0.12 (0.15 %) when using 16 × 16 × 16 k-points with respect to the converged result.This convergence behavior is comparable or slightly slower than what is observed for the total energy.

Table 3 .
Comparison of the high-frequency dielectric constants of various semiconductors computed at the LDA/PBE level with literature values: Tight-default settings and basis sets as well as a 16× 16× 16 k-point mesh are used.
Figure 6.H(C 2 H 4 ) n H molecules: CPU time per PT cycle required for finite H(C 2 H 4

Table A1 .
on the NOMAD Repository.P.R. acknowledges financial support from the Academy of Finland through its Centres of Excellence Program (Project No. 251748 and 284621).Polarizability tensor elements α ii for 16 dimers, as computed with the presented PT implementation at the LDA level of theory.Additionally, absolute errors (AE) and absolute percentage erros (APE) with respect to finite difference calculations are given.Please note that this errors are several orders of magnitude smaller than the relevant digits in α ii .

Table A2 .
(continued) Polarizability tensor elements α ii for 16 dimers, as computed with the presented PT implementation at the LDA level of theory.Additionally, absolute errors (AE) and absolute percentage erros (APE) with respect to finite difference calculations are given.Please note that this errors are several orders of magnitude smaller than the relevant digits in α ii .

Table A3 .
Polarizability tensor elements α ii for five trimers, as computed with the presented PT implementation at the LDA level of theory.Additionally, absolute errors (AE) and absolute percentage erros (APE) with respect to finite difference calculations are given.Please note that this errors are several orders of magnitude smaller than the relevant digits in α ii .

Table A4 .
Polarizability tensor elements α ii for eleven molecules, as computed with the presented PT implementation at the LDA level of theory.Additionally, absolute errors (AE) and absolute percentage erros (APE) with respect to finite difference calculations are given.Please note that this errors are several orders of magnitude smaller than the relevant digits in α ii .