All-electron, real-space perturbation theory for homogeneous electric fields: theory, implementation, and application within DFT

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 sum-frequency generation spectra. In first-principles frameworks, these quantities are typically computed via time-dependent density-functional theory (DFT) [1, 2] or via analytical perturbation theory (PT) in either its density-functional perturbation theory (DFPT) [3–6] or coupled perturbed self-consistent field (CPSCF) formulation [7–12]. Within these linear-response approaches, an additional complexity arises for the treatment of solids: as discussed in more detail in section 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 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 section 6. Our implementation uses the all-electron, numeric atom-centered orbital based framework of the FHI-aims code [13–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 sections 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 section 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 PT framework is discussed for finite and extended periodic systems in sections 2 and 3, respectively. In section 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 section 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 section 6 by applying the developed formalism to compute harmonic and anharmonic Raman spectra for different polymorphs of the paracetamol crystal.

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 section 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 section 3.

In Kohn–Sham DFT, the total energy functional is given by

$$\begin{eqnarray}&&{E}_{{\rm{KS}}}[n]={T}_{{\rm{s}}}[n]+{E}_{{\rm{ext}}}[n]+{E}_{{\rm{H}}}[n]+{E}_{{\rm{xc}}}[n]+{E}_{{\rm{nuc-nuc}}}.\end{eqnarray} \tag{ 1 }$$

Here, n(r) is the electron density, ${T}_{{\rm{s}}}$ the kinetic energy of non-interacting electrons, ${E}_{{\rm{ext}}}$ the external energy due to the electron-nuclear attraction, ${E}_{{\rm{H}}}$ the Hartree energy, ${E}_{{\rm{xc}}}$ the exchange-correlation energy, and ${E}_{{\rm{nuc-nuc}}}$ the repulsion energy of the nuclei. The ground-state electron density ${n}_{0}({\bf{r}})$ (and the associated ground-state total energy) is obtained by variationally minimizing equation (1) under the constraint that the number of electrons N_e is conserved. This yields the chemical potential $\mu =\delta {E}_{{\rm{KS}}}/\delta n$ of the electrons and the Kohn–Sham single particle equations

$$\begin{eqnarray}&&{\hat{h}}_{{\rm{KS}}}{\psi }_{p}=[{\hat{t}}_{{\rm{s}}}+{\hat{v}}_{{\rm{ext}}}+{\hat{v}}_{{\rm{H}}}+{\hat{v}}_{{\rm{xc}}}]{\psi }_{p}={\epsilon }_{p}{\psi }_{p}\end{eqnarray} \tag{ 2 }$$

for the Kohn–Sham Hamiltonian ${\hat{h}}_{{\rm{KS}}}$. In equation (2), ${\hat{t}}_{{\rm{s}}}$ denotes the kinetic energy operator, ${\hat{v}}_{{\rm{ext}}}$ the external potential, ${\hat{v}}_{{\rm{H}}}$ the Hartree potential, and ${\hat{v}}_{{\rm{xc}}}$ the exchange-correlation potential. Solving equation (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

$$\begin{eqnarray}&&n({\bf{r}})=\displaystyle \sum _{p}f({\epsilon }_{p})| {\psi }_{p}({\bf{r}}){| }^{2},\end{eqnarray} \tag{ 3 }$$

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 equation (2) in numerical implementations, the Kohn–Sham states are expanded in a finite basis set ${\chi }_{\mu }({\bf{r}}-{{\bf{R}}}_{I(\mu )})$ as

$$\begin{eqnarray}&&{\psi }_{p}({\bf{r}})=\displaystyle \sum _{\mu }{C}_{\mu p}\ {\chi }_{\mu }({\bf{r}}-{{\bf{R}}}_{I(\mu )}),\end{eqnarray} \tag{ 4 }$$

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 ${{\bf{R}}}_{I(\mu )}$. In such a basis set, equation (2) becomes a generalized eigenvalue problem

$$\begin{eqnarray}&&\displaystyle \sum _{\nu }{H}_{\mu \nu }{C}_{\nu p}={\epsilon }_{p}\displaystyle \sum _{\nu }{S}_{\mu \nu }{C}_{\nu p}.\end{eqnarray} \tag{ 5 }$$

Using the bra-ket notation $\langle .| .\rangle$ for the inner product in Hilbert space, H_μν denotes the elements $\langle {\chi }_{\mu }| {\hat{h}}_{{\rm{KS}}}| {\chi }_{\nu }\rangle$ of the Hamiltonian matrix and S_μν the elements $\langle {\chi }_{\mu }| {\chi }_{\nu }\rangle$ of the overlap matrix. Accordingly, the variation with respect to the density becomes a minimization with respect to the expansion coefficients C_νp

$$\begin{eqnarray}&&{E}_{{\rm{tot}}}^{0}={E}_{{\rm{KS}}}[{n}_{0}({\bf{r}})]=\mathop{\min }\limits_{{C}_{\nu p}}\left[{E}_{{\rm{KS}}}-\displaystyle \sum _{p}{\epsilon }_{p}(\langle {\psi }_{p}| {\psi }_{p}\rangle -1)\right],\end{eqnarray} \tag{ 6 }$$

in which the eigenstates ψ_p are constrained to be orthonormal. Typically, the ground-state density ${n}_{0}({\bf{r}})$ and the associated total energy ${E}_{{\rm{tot}}}$ are determined by solving equation (6) iteratively, until self-consistency is achieved.

If an external electric field ${\bf{E}}$ is applied to an isolated system, the KS Hamiltonian gains an additional term ${\hat{h}}_{{\rm{E}}}=-{\bf{r}}\cdot {\bf{E}}$. If this electric field ${\bf{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 ${\hat{h}}_{{\rm{E}}}$ contributes

$$\begin{eqnarray}&&{E}_{{\rm{E}}}[n]=-\displaystyle \sum _{\gamma }\int {e}_{\gamma }{r}_{\gamma }\,n({\bf{r}})\,{\rm{d}}{\bf{r}}\end{eqnarray} \tag{ 7 }$$

to the total energy functional in equation (1). A perturbative Taylor expansion of the total energy in the zero-field limit then gives

$$\begin{eqnarray}&&{E}_{{\rm{tot}}}({\bf{E}})\approx {E}_{{\rm{tot}}}^{0}+\displaystyle \sum _{\gamma }{\mu }_{\gamma }{e}_{\gamma }+\displaystyle \frac{1}{2}\displaystyle \sum _{\gamma ,\delta }{\alpha }_{\gamma \delta }{e}_{\gamma }{e}_{\delta }\,+\cdots ,\end{eqnarray} \tag{ 8 }$$

where η, γ are Cartesian directions. For isolated systems, the coefficient in the linear term

$$\begin{eqnarray}&&{\mu }_{\gamma }={\left.\displaystyle \frac{\partial {E}_{{\rm{E}}}[{n}_{0}]}{\partial {e}_{\gamma }}\right|}_{{\bf{E}}=0}=-\int {n}_{0}({\bf{r}}){r}_{\gamma }{\rm{d}}{\bf{r}},\end{eqnarray} \tag{ 9 }$$

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

$$\begin{eqnarray}&&{\alpha }_{\gamma \delta }={\left.\displaystyle \frac{{\partial }^{2}{E}_{{\rm{E}}}[{n}_{0}]}{\partial {e}_{\gamma }\partial {e}_{\delta }}\right|}_{{\bf{E}}=0}={\left.\displaystyle \frac{\partial {\mu }_{\gamma }}{\partial {e}_{\delta }}\right|}_{{\bf{E}}=0}=-\int {r}_{\gamma }\left({\left.\displaystyle \frac{\partial {n}_{0}({\bf{r}})}{\partial {e}_{\delta }}\right|}_{{\bf{E}}=0}\right){\rm{d}}{\bf{r}},\end{eqnarray} \tag{ 10 }$$

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–23] cannot be directly 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 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

$$\begin{eqnarray}&&{M}^{(1)}=\displaystyle \frac{{{\rm{d}}{M}}^{(0)}}{{{\rm{d}}{e}}_{\gamma }}\end{eqnarray} \tag{ 11 }$$

used in the following for ground-state M⁽⁰⁾ and response properties M⁽¹⁾ is always well-defined, e.g.,

$$\begin{eqnarray}&&{n}^{(1)}=\displaystyle \frac{{{\rm{d}}{n}}^{(0)}}{{{\rm{d}}{e}}_{\gamma }}=\displaystyle \frac{{{\rm{d}}{n}}_{0}}{{{\rm{d}}{e}}_{\gamma }}.\end{eqnarray} \tag{ 12 }$$

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:

$$\begin{eqnarray}&&{\hat{h}}_{{\rm{KS}}}({e}_{\gamma })\approx {\hat{h}}_{{\rm{KS}}}^{(0)}+{\hat{h}}_{{\rm{KS}}}^{(1)}{e}_{\gamma }\,+\cdots ,\end{eqnarray} \tag{ 13 }$$

where the response of the Hamiltonian operator is

$$\begin{eqnarray}&&{\hat{h}}_{{\rm{KS}}}^{(1)}={\hat{v}}_{{\rm{ext}}}^{(1)}(r)+{\hat{v}}_{{\rm{H}}}^{(1)}+{\hat{v}}_{{\rm{xc}}}^{(1)}-{r}_{\gamma }.\end{eqnarray} \tag{ 14 }$$

Introducing the analogous expansions

$$\begin{eqnarray}&&{\psi }_{p}({e}_{\gamma })\approx {\psi }_{p}^{(0)}+{\psi }_{p}^{(1)}{e}_{\gamma }+\cdots \quad \quad \quad {\epsilon }_{p}({e}_{\gamma })\approx {\epsilon }_{p}^{(0)}+{\epsilon }_{p}^{(1)}{e}_{\gamma }+\cdots \end{eqnarray} \tag{ 15 }$$

for the single particle states ψ_p(e_γ) and their eigenvalues _p(e_γ), rearranging the linear-order terms in the KS equation ${\hat{h}}_{{\rm{KS}}}({e}_{\gamma }){\psi }_{p}({e}_{\gamma })={\epsilon }_{p}({e}_{\gamma }){\psi }_{p}({e}_{\gamma })$, and applying the normalization condition $\langle {\psi }_{p}({e}_{\gamma })| {\psi }_{p}({e}_{\gamma })\rangle =1$, yields the Sternheimer equation

$$\begin{eqnarray}&&({\hat{h}}_{{\rm{KS}}}^{(0)}-{\epsilon }_{p}^{(0)})| {\psi }_{p}^{(1)}\rangle =-({\hat{h}}_{{\rm{KS}}}^{(1)}-{\epsilon }_{p}^{(1)})| {\psi }_{p}^{(0)}\rangle ,\end{eqnarray} \tag{ 16 }$$

as well as the condition

$$\begin{eqnarray}&&\langle {\psi }_{p}^{(1)}| {\psi }_{p}^{(0)}\rangle +\langle {\psi }_{p}^{(0)}| {\psi }_{p}^{(1)}\rangle =0.\end{eqnarray} \tag{ 17 }$$

By multiplying equation (16) with $\langle {\psi }_{q}^{(0)}|$ from the left one obtains

$$\begin{eqnarray}&&({\epsilon }_{q}^{(0)}-{\epsilon }_{p}^{(0)})\langle {\psi }_{q}^{(0)}| {\psi }_{p}^{(1)}\rangle =-(\langle {\psi }_{q}^{(0)}| {\hat{h}}_{{\rm{KS}}}^{(1)}| {\psi }_{p}^{(0)}\rangle -{\epsilon }_{p}^{(1)}{\delta }_{{qp}}).\end{eqnarray} \tag{ 18 }$$

To solve this equation numerically, we expand the response of the wave functions

$$\begin{eqnarray}&&{\psi }_{p}^{(1)}({\bf{r}})=\displaystyle \sum _{q}{U}_{{qp}}^{(1)}{\psi }_{q}^{(0)}({\bf{r}})=\displaystyle \sum _{\mu }\mathop{\underbrace{\displaystyle \sum _{q}{U}_{{qp}}^{(1)}{C}_{\mu q}^{(0)}}}\limits_{{C}_{\mu p}^{(1)}}\ {\chi }_{\mu }({\bf{r}})\end{eqnarray} \tag{ 19 }$$

in terms of the unperturbed states ${\psi }_{q}^{(0)}({\bf{r}})$. Here, we chose ${U}_{{pp}}^{(1)}=0$ for all p to fulfill equation (17) and hence obtain an algebraic expression for equation (18)

$$\begin{eqnarray}&&({\epsilon }_{q}^{(0)}-{\epsilon }_{p}^{(0)}){U}_{{qp}}^{(1)}=-\displaystyle \sum _{\mu \nu }{({C}_{\mu q}^{(0)})}^{* }{C}_{\nu p}^{(0)}\langle {\chi }_{\mu }| {\hat{h}}_{{\rm{KS}}}^{(1)}| {\chi }_{\nu }\rangle +{\epsilon }_{p}^{(1)}{\delta }_{{qp}}.\end{eqnarray} \tag{ 20 }$$

The expansion using the matrix ${U}_{{qp}}^{(1)}$ employed in this work is typical for the CPSCF formulation [7–12] of PT. For our implementation described in section 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}_{\mu \nu }^{(1)}=\langle {\chi }_{\mu }| {\hat{h}}_{{\rm{KS}}}^{(1)}| {\chi }_{\nu }\rangle$ are defined in the same way as for unperturbed calculations, i.e., using the numeric atomic orbitals introduced in equation (4). This allows us to directly compute the non-diagonal elements ($q\ne p$) of

$$\begin{eqnarray}&&{U}_{{qp}}^{(1)}=\displaystyle \frac{{\displaystyle \sum }_{\mu \nu }{({C}_{\mu q}^{(0)})}^{* }\ {H}_{\mu \nu }^{(1)}\ {C}_{\nu p}^{(0)}}{{\epsilon }_{p}^{(0)}-{\epsilon }_{q}^{(0)}}.\end{eqnarray} \tag{ 21 }$$

The matrix ${U}_{{qp}}^{(1)}$, which fulfills ${U}_{{qp}}^{(1)}=-{({U}_{{pq}}^{(1)})}^{* }$, plays a central role in our implementation: As discussed in detail in section 4, it allows us to directly determine the response of the density

$$\begin{eqnarray}&&{n}^{(1)}=\displaystyle \sum _{p}f({\epsilon }_{p})[{\psi }_{p}^{(1)}{\psi }_{p}^{(0)}+{\psi }_{p}^{(0)}{\psi }_{p}^{(1)}]\end{eqnarray} \tag{ 22 }$$

in a density matrix formalism, i.e., without explicitly computing the response of the eigenvalues ${\epsilon }_{p}^{(1)}$, of the wavefunction ${\psi }_{p}^{(1)}({\bf{r}})$, or its coefficients ${C}_{\mu p}^{(1)}$, which is computationally advantageous. Using ${U}_{{qp}}^{(1)}$, one can then directly evaluate the polarizability tensor α_γδ defined in equation (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 equation (19) via the ${U}_{{qp}}^{(1)}$ matrix, but rather compute the coefficients ${C}_{\mu p}^{(1)}$ by directly solving equation (18) in the space spanned by the KS states using the DFPT formalism [3–6]. In such codes, in which thousands of orbitals, i.e., plane waves, need to be considered, the DFPT approach is advantageous.

For periodic boundary conditions (PBCs), the main physical reasoning behind the derivation of equations (1)–(22) still remains valid. However, three specific adaptations have to be made:

First, the basis set expansion introduced in equation (4) is slightly different, as described in detail in [13, 16, 24]: The periodic images of the nuclei ${{\bf{R}}}_{{Im}}={{\bf{R}}}_{I}+{{\bf{R}}}_{m}$ are accounted for by summing over the lattice vectors ${{\bf{R}}}_{m}$, i.e., over linear combinations of the unit cell (u.c.) lattice vectors ${{\bf{a}}}_{1},{{\bf{a}}}_{2},{{\bf{a}}}_{3}$. Analogously, also the numeric atomic orbitals associated with such periodic images, e.g., ${\chi }_{\mu m}({\bf{r}})={\chi }_{\mu }({\bf{r}}-{{\bf{R}}}_{I(\mu )}-{{\bf{R}}}_{m})$ associated with the periodic image m of nucleus ${{\bf{R}}}_{I}$, gain an additional index m that describes their relative position to the u.c. equivalent. To account for translational symmetry and exploit Bloch's theorem, Bloch-like generalized basis functions

$$\begin{eqnarray}&&{\varphi }_{\mu }({\bf{k}},{\bf{r}})=\displaystyle \sum _{m}{\chi }_{\mu m}({\bf{r}})\exp (-{\rm{i}}{\bf{k}}\cdot {{\bf{R}}}_{m})\end{eqnarray} \tag{ 23 }$$

are constructed from the local atomic orbitals and then used in the basis set expansion

$$\begin{eqnarray}&&{\psi }_{p}^{(0)}({\bf{k}},{\bf{r}})=\displaystyle \sum _{\mu }{C}_{\mu p}^{(0)}({\bf{k}})\ {\varphi }_{\mu }({\bf{k}},{\bf{r}}).\end{eqnarray} \tag{ 24 }$$

Accordingly, all relevant physical quantities such as the KS Hamiltonian

$$\begin{eqnarray}&&{H}_{\mu \nu }^{(0)}({\bf{k}})=\displaystyle \sum _{m,n}{{\rm{e}}}^{-{\rm{i}}{\bf{k}}({{\bf{R}}}_{n}-{{\bf{R}}}_{m})}{\int }_{{\rm{u.c.}}}{\chi }_{\mu m}({\bf{r}}){\hat{h}}_{{\rm{KS}}}\,{\chi }_{\nu n}({\bf{r}}){\rm{d}}{\bf{r}}\end{eqnarray} \tag{ 25 }$$

gain an additional dependence on the wavevector ${\bf{k}}$, so that equation (5) becomes

$$\begin{eqnarray}&&\displaystyle \sum _{\nu }{H}_{\mu \nu }^{(0)}({\bf{k}}){C}_{\nu p}^{(0)}({\bf{k}})={\epsilon }_{p}^{(0)}({\bf{k}})\displaystyle \sum _{\nu }{S}_{\mu \nu }({\bf{k}}){C}_{\nu p}^{(0)}({\bf{k}}).\end{eqnarray} \tag{ 26 }$$

Therefore, the summations over electronic states appearing in equations (1)–(22) now feature an additional analytical integration over the Brillouin zone that is approximated numerically by a sum over a finite ${\bf{k}}$-grid with N_k points. Similarly, the real-space integrals in equations (1)–(22) are no longer indefinite, but definite and limited to the u.c., as it is the case in equation (25).

Second, it is necessary to consider the screened electric field ${\bf{E}}={\bf{D}}-4\pi {\bf{P}}$ in the solid, where ${\bf{D}}$ is the electric displacement [25] and the polarization in the u.c. volume V is given by [3, 6]:

$$\begin{eqnarray}&&{P}_{\gamma }=-\displaystyle \frac{1}{V}{\int }_{{\rm{u.c.}}}{r}_{\gamma }{n}_{0}({\bf{r}}){\rm{d}}{\bf{r}}.\end{eqnarray} \tag{ 27 }$$

The relationship between the components of the electric displacement and the screened field defines the high-frequency dielectric constant [25]

$$\begin{eqnarray}&&{\varepsilon }_{\gamma \delta }^{\infty }=\displaystyle \frac{\partial {D}_{\gamma }}{\partial {E}_{\delta }}={\delta }_{\gamma \delta }+4\pi \displaystyle \frac{\partial {P}_{\gamma }}{\partial {E}_{\delta }},\end{eqnarray} \tag{ 28 }$$

where η, γ are Cartesian directions and δ is the Kronecker delta symbol. For a screened field ${\bf{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

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial {E}_{{\rm{E}}}[{n}_{0}]}{\partial {e}_{\gamma }}\right|}_{{\bf{E}}=0}=-{\int }_{{\rm{u.c.}}}{r}_{\gamma }{n}_{0}({\bf{r}}){\rm{d}}{\bf{r}}={{VP}}_{\gamma }\to {\mu }_{\gamma },\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{\partial }^{2}{E}_{{\rm{E}}}[{n}_{0}]}{\partial {e}_{\gamma }\partial {e}_{\delta }}\right|}_{{\bf{E}}=0}=-{\int }_{{\rm{u.c.}}}{r}_{\gamma }\left({\left.\displaystyle \frac{\partial {n}_{0}({\bf{r}})}{\partial {e}_{\delta }}\right|}_{{\bf{E}}=0}\right){\rm{d}}{\bf{r}}=V{\left.\displaystyle \frac{\partial {P}_{\gamma }}{\partial {e}_{\delta }}\right|}_{{\bf{E}}=0}\to {\alpha }_{\gamma \delta }.\end{eqnarray} \tag{ 30 }$$

A comparison with equations (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 ${\bf{E}}$ is not periodic, as alluded to in the introduction. As a consequence, the definite integral over the u.c. required to determine the Hamiltonian response ${H}_{\mu \nu }^{(1)}({\bf{k}})$ is ill-defined in PBCs, since ${\hat{h}}_{{\rm{KS}}}^{(1)}$ given in equation (14) contains the position operator r_γ, which is itself ill-defined in this case. The same problem affects equation (30). In reciprocal-space implementations, the Berry-phase formalism [4, 5, 21, 23, 27, 28] is typically the method of choice. A tutorial introduction to this approach can be found in [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 $[{\hat{h}}_{{\rm{KS}}}^{(0)}({\bf{k}}),{\bf{r}}]=-{\rm{\nabla }}$. With that, one gets the well-known expression

$$\begin{eqnarray}&&\langle {\psi }_{q}^{(0)}({\bf{k}})| {{\rm{\nabla }}}_{\gamma }| {\psi }_{p}^{(0)}({\bf{k}})\rangle =-\langle {\psi }_{q}^{(0)}({\bf{k}})| [{\hat{h}}_{{\rm{KS}}}^{(0)}({\bf{k}}),{r}_{\gamma }]| {\psi }_{p}^{(0)}({\bf{k}})\rangle \end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&=({\epsilon }_{p}^{(0)}({\bf{k}})-{\epsilon }_{q}^{(0)}({\bf{k}}))\langle {\psi }_{q}^{(0)}({\bf{k}})| {r}_{\gamma }| {\psi }_{p}^{(0)}({\bf{k}})\rangle ,\end{eqnarray} \tag{ 32 }$$

that can be used to evaluate the non-diagonal matrix elements ($q\ne p$)

$$\begin{eqnarray}&&\langle {\psi }_{q}^{(0)}({\bf{k}})| {r}_{\gamma }| {\psi }_{p}^{(0)}({\bf{k}})\rangle =\displaystyle \frac{\langle {\psi }_{q}^{(0)}({\bf{k}})| {{\rm{\nabla }}}_{\gamma }| {\psi }_{p}^{(0)}({\bf{k}})\rangle }{{\epsilon }_{p}^{(0)}({\bf{k}})-{\epsilon }_{q}^{(0)}({\bf{k}})}.\end{eqnarray} \tag{ 33 }$$

Using equations (23) and (24) we obtain the representation

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{qp}}({\bf{k}})=-\langle {\psi }_{q}^{(0)}({\bf{k}})| {r}_{\gamma }| {\psi }_{p}^{(0)}({\bf{k}})\rangle =-\displaystyle \sum _{\mu \nu }\displaystyle \frac{{({C}_{\mu q}^{(0)}({\bf{k}}))}^{* }{C}_{\nu p}^{(0)}({\bf{k}})}{{\epsilon }_{p}^{(0)}({\bf{k}})-{\epsilon }_{q}^{(0)}({\bf{k}})}{R}_{\mu \nu }^{(0)}({\bf{k}})\end{eqnarray} \tag{ 34 }$$

with

$$\begin{eqnarray}&&{R}_{\mu \nu }^{(0)}({\bf{k}})=\displaystyle \sum _{{mn}}{{\rm{e}}}^{-{\rm{i}}{\bf{k}}({{\bf{R}}}_{n}-{{\bf{R}}}_{m})}{\int }_{{\rm{u.c.}}}{\chi }_{\mu m}({\bf{r}}){{\rm{\nabla }}}_{\gamma }\,{\chi }_{\nu n}({\bf{r}}){\rm{d}}{\bf{r}}.\end{eqnarray} \tag{ 35 }$$

We can then recast the expectation value ${H}_{{qp}}^{(1)}({\bf{k}})=\langle {\psi }_{q}^{(0)}({\bf{k}})| {\hat{h}}_{{\rm{KS}}}^{(1)}| {\psi }_{p}^{(0)}({\bf{k}})\rangle$ appearing in equation (18):

$$\begin{eqnarray}&&\langle {\psi }_{q}^{(0)}({\bf{k}})| {\hat{h}}_{{\rm{KS}}}^{(1)}| {\psi }_{p}^{(0)}({\bf{k}})\rangle =\langle {\psi }_{q}^{(0)}({\bf{k}})| {\hat{v}}_{{\rm{ext}}}^{(1)}(r)+{\hat{v}}_{{\rm{H}}}^{(1)}+{\hat{v}}_{{\rm{xc}}}^{(1)}| {\psi }_{p}^{(0)}({\bf{k}})\rangle -\langle {\psi }_{q}^{(0)}({\bf{k}})| {r}_{\gamma }| {\psi }_{p}^{(0)}({\bf{k}})\rangle \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&=\displaystyle \sum _{\mu \nu }{({C}_{\mu q}^{(0)}({\bf{k}}))}^{* }{C}_{\nu p}^{(0)}({\bf{k}}){V}_{\mu \nu }^{(1)}({\bf{k}})+{{\rm{\Omega }}}_{{qp}}({\bf{k}}).\end{eqnarray} \tag{ 37 }$$

Here, the matrix elements ${V}_{\mu \nu }^{(1)}({\bf{k}})=\langle {\varphi }_{\mu }({\bf{k}})| {\hat{v}}_{{\rm{ext}}}^{(1)}(r)+{\hat{v}}_{{\rm{H}}}^{(1)}+{\hat{v}}_{{\rm{xc}}}^{(1)}| {\varphi }_{\nu }({\bf{k}})\rangle$ can be directly evaluated as done in equation (25)

$$\begin{eqnarray}&&{V}_{\mu \nu }^{(1)}({\bf{k}})=\displaystyle \sum _{m,n}{{\rm{e}}}^{-{\rm{i}}{\bf{k}}({{\bf{R}}}_{n}-{{\bf{R}}}_{m})}{\int }_{{\rm{u.c.}}}{\chi }_{\mu m}({\bf{r}})({\hat{v}}_{{\rm{ext}}}^{(1)}(r)+{\hat{v}}_{{\rm{H}}}^{(1)}+{\hat{v}}_{{\rm{xc}}}^{(1)})\,{\chi }_{\nu n}({\bf{r}}){\rm{d}}{\bf{r}},\end{eqnarray} \tag{ 38 }$$

since they only depend on lattice periodic operators. Now, the matrix ${U}_{{qp}}^{(1)}({\bf{k}})$ introduced in equation (21) is computed as

$$\begin{eqnarray}&&{U}_{{qp}}^{(1)}({\bf{k}})=\displaystyle \frac{{\displaystyle \sum }_{\mu \nu }{({C}_{\mu q}^{(0)}({\bf{k}}))}^{* }\ {V}_{\mu \nu }^{(1)}({\bf{k}})\ {C}_{\nu p}^{(0)}({\bf{k}})}{{\epsilon }_{p}^{(0)}({\bf{k}})-{\epsilon }_{q}^{(0)}({\bf{k}})}+\displaystyle \frac{{{\rm{\Omega }}}_{{qp}}({\bf{k}})}{{\epsilon }_{p}^{(0)}({\bf{k}})-{\epsilon }_{q}^{(0)}({\bf{k}})}.\end{eqnarray} \tag{ 39 }$$

Similarly, the polarizability tensor components appearing in equation (10) can be rewritten as

$$\begin{eqnarray}{\alpha }_{\gamma \delta } & = & -{\displaystyle \int }_{{\rm{u.c.}}}{r}_{\gamma }\left({\left.\displaystyle \frac{\partial {n}_{0}({\bf{r}})}{\partial {e}_{\delta }}\right|}_{{\bf{E}}=0}\right){\rm{d}}{\bf{r}}\\ & = & -\displaystyle \frac{1}{{N}_{k}}\displaystyle \sum _{p,{\bf{k}}}f({\epsilon }_{p}({\bf{k}}))[\langle {\psi }_{p}^{(1,\delta )}({\bf{k}})| {r}_{\gamma }| {\psi }_{p}^{(0)}({\bf{k}})\rangle +\langle {\psi }_{p}^{(0)}({\bf{k}})| {r}_{\gamma }| {\psi }_{p}^{(1,\delta )}({\bf{k}})\rangle ]\\ & = & \displaystyle \frac{1}{{N}_{k}}\displaystyle \sum _{q,p,{\bf{k}}}f({\epsilon }_{p}({\bf{k}}))[{({U}_{{qp}}^{(1,\delta )}({\bf{k}}))}^{* }{{\rm{\Omega }}}_{{qp}}^{(\gamma )}({\bf{k}})\ +\ {U}_{{qp}}^{(1,\delta )}({\bf{k}}){{\rm{\Omega }}}_{{pq}}^{(\gamma )}({\bf{k}})],\end{eqnarray} \tag{ 40 }$$

using the matrix elements defined in equation (34). As explicitly highlighted in the notation, the matrix ${U}_{{qp}}^{(1,\delta )}({\bf{k}})$ associated with a perturbation along the Cartesian axis δ has to be used in this case, whereas the matrix ${{\rm{\Omega }}}_{{qp}}^{(\gamma )}({\bf{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 ${\varphi }_{\mu }({\bf{k}})$ and wave vectors ${\bf{k}}$ is used, since a simplification to finite systems is straightforward.

Our implementation closely follows the flowchart shown in figure 1: After a ground-state DFT calculation (see [13]) is completed, the matrix ${{\rm{\Omega }}}_{{qp}}({\bf{k}})$ is computed. If ${U}_{{pq}}^{(1)}({\bf{k}})=0$ is used as initial guess, one obtains

$$\begin{eqnarray}&&{U}_{{qp}}^{(1)}({\bf{k}})=\displaystyle \frac{{{\rm{\Omega }}}_{{qp}}({\bf{k}})}{{\epsilon }_{p}^{(0)}({\bf{k}})-{\epsilon }_{q}^{(0)}({\bf{k}})}\end{eqnarray} \tag{ 41 }$$

in the first iteration, which can then be fed back to the self-consistency loop to determine the first-order density response ${n}^{(1)}({\bf{r}})$ in a density matrix formalism (see section 4.1). As detailed in section 4.2, we then use ${n}^{(1)}({\bf{r}})$ to compute the remaining, individual ingredients that enter $\langle {\psi }_{q}^{(0)}({\bf{k}})| {\hat{h}}_{{\rm{KS}}}^{(1)}| {\psi }_{p}^{(0)}({\bf{k}})\rangle$, i.e., the matrix elements ${V}_{\mu \nu }^{(1)}({\bf{k}})$ defined in equation (38). The Sternheimer equation then provides a new matrix ${U}_{{qp}}^{(1)}({\bf{k}})$, as discussed in section 4.3. We iteratively restart the PT loop using a Pulay mixer [30], until self-consistency is reached, i.e., until the changes in the response of the density matrix ${{\bf{P}}}^{(1)}$ become smaller than a user-given threshold. In the last step, the polarizability and the dielectric constant are computed, as discussed in section 4.4. Atomic units are used in the complete workflow.

Zoom In Zoom Out Reset image size

Both the ground-state density n⁽⁰⁾(r) and the response of the density n⁽¹⁾(r) are periodic, i.e., invariant against translations

$$\begin{eqnarray}&&{n}^{(0)}({\bf{r}}+{{\bf{R}}}_{m})={n}^{(0)}({\bf{r}})\quad \quad {n}^{(1)}({\bf{r}}+{{\bf{R}}}_{m})={n}^{(1)}({\bf{r}})\end{eqnarray} \tag{ 42 }$$

by a lattice vector ${{\bf{R}}}_{m}$, as shown in figure 2. Accordingly, we can use the algorithms used in ground-state calculations and discussed in detail in [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.

Zoom In Zoom Out Reset image size

4.1. Response of the electronic density

To numerically compute the electronic density $n({\bf{r}})$ in ground-state calculations [13], we use a density matrix formalism

$$\begin{eqnarray}&&{n}^{(0)}({\bf{r}})=\displaystyle \sum _{\mu m,\nu n}{P}_{\mu m,\nu n}^{(0)}{\chi }_{\mu m}^{(0)}({\bf{r}}){\chi }_{\nu n}^{(0)}({\bf{r}}),\end{eqnarray} \tag{ 43 }$$

which is obtained by inserting equations (24) and (23) into equation (3). Hence, the density matrix is given by

$$\begin{eqnarray}&&{P}_{\mu m,\nu n}^{(0)}=\displaystyle \frac{1}{{N}_{k}}\displaystyle \sum _{{\bf{k}}}\left[{{\rm{e}}}^{-{\rm{i}}{\bf{k}}({{\bf{R}}}_{n}-{{\bf{R}}}_{m})}\displaystyle \sum _{o}f({\epsilon }_{o}({\bf{k}})){({C}_{\mu o}^{(0)}({\bf{k}}))}^{* }{C}_{\nu o}^{(0)}({\bf{k}})\right].\end{eqnarray} \tag{ 44 }$$

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({\epsilon }_{o}({\bf{k}}))\ne 0$ in practice. Similarly, the response of the electronic density can thus be expressed as

$$\begin{eqnarray}&&{n}^{(1)}({\bf{r}})=\displaystyle \sum _{\mu m,\nu n}{P}_{\mu m,\nu n}^{(1)}{\chi }_{\mu m}^{(0)}({\bf{r}}){\chi }_{\nu n}^{(0)}({\bf{r}})\end{eqnarray} \tag{ 45 }$$

using the response of the density matrix given by

$$\begin{eqnarray}{P}_{\mu m,\nu n}^{(1)} & = & \displaystyle \frac{1}{{N}_{k}}\displaystyle \sum _{{\bf{k}}}\left\{{{\rm{e}}}^{-{\rm{i}}{\bf{k}}({{\bf{R}}}_{n}-{{\bf{R}}}_{m})}\displaystyle \sum _{o}f({\epsilon }_{o}({\bf{k}}))[{({C}_{\mu o}^{(1)}({\bf{k}}))}^{\ast }{C}_{\nu o}^{(0)}({\bf{k}})\right.+{({C}_{\mu o}^{(0)}({\bf{k}}))}^{\ast }{C}_{\nu o}^{(1)}({\bf{k}})]\}.\end{eqnarray} \tag{ 46 }$$

In the sum over states o, we express ${{\bf{C}}}^{(1)}({\bf{k}})$ in terms of ${{\bf{U}}}^{(1)}({\bf{k}})$ via

$$\begin{eqnarray}{C}_{\mu p}^{(1)}({\bf{k}}) & = & \displaystyle \sum _{q}{U}_{{qp}}^{(1)}({\bf{k}}){C}_{\mu q}^{(0)}({\bf{k}})\\ & = & \displaystyle \frac{1}{2}\displaystyle \sum _{o^{\prime} }f({\epsilon }_{o^{\prime} }({\bf{k}})){U}_{o^{\prime} p}^{(1)}({\bf{k}}){C}_{\mu o^{\prime} }^{(0)}({\bf{k}})+\displaystyle \frac{1}{2}\displaystyle \sum _{u}(2-f({\epsilon }_{u}({\bf{k}}))){U}_{{up}}^{(1)}({\bf{k}}){C}_{\mu u}^{(0)}({\bf{k}}),\end{eqnarray} \tag{ 47 }$$

whereby we split the sum over q into two separate sums over $o^{\prime}$ and u. In practice, these two sums can then be later evaluated by restricting the sum over $o^{\prime}$ to occupied and the sum over u to unoccupied states, respectively. Accordingly, also the sum over o appearing in equation (46) can be split into two double sums, one over $o,o^{\prime}$ and one over $o,u$. For the first one, we obtain

$$\begin{eqnarray}&&\displaystyle \sum _{o,o^{\prime} }\displaystyle \frac{f({\epsilon }_{o}({\bf{k}}))f({\epsilon }_{o^{\prime} }({\bf{k}}))}{2}[{({C}_{\mu o^{\prime} }^{(0)}({\bf{k}}){U}_{o^{\prime} o}^{(1)}({\bf{k}}))}^{* }{C}_{\nu o}^{(0)}({\bf{k}})+{({C}_{\mu o}^{(0)}({\bf{k}}))}^{* }{U}_{o^{\prime} o}^{(1)}({\bf{k}}){C}_{\nu o^{\prime} }^{(0)}({\bf{k}})]\\ &&=\displaystyle \sum _{o,o^{\prime} }\displaystyle \frac{f({\epsilon }_{o}({\bf{k}}))f({\epsilon }_{o^{\prime} }({\bf{k}}))}{2}[{({C}_{\mu o^{\prime} }^{(0)}({\bf{k}}){U}_{o^{\prime} o}^{(1)}({\bf{k}}))}^{* }{C}_{\nu o}^{(0)}({\bf{k}})+{({C}_{\mu o^{\prime} }^{(0)}({\bf{k}}))}^{* }{U}_{{oo}^{\prime} }^{(1)}({\bf{k}}){C}_{\nu o}^{(0)}({\bf{k}})]=0,\end{eqnarray} \tag{ 48 }$$

which vanishes due to ${U}_{{oo}^{\prime} }^{(1)}({\bf{k}})=-{({U}_{o^{\prime} o}^{(1)}({\bf{k}}))}^{* }$, see equation (39). For the second double sum, we obtain

$$\begin{eqnarray}&&\displaystyle \sum _{o,u}\displaystyle \frac{f({\epsilon }_{o}({\bf{k}}))(2-f({\epsilon }_{u}({\bf{k}})))}{2}[{({C}_{\mu u}^{(0)}({\bf{k}}){U}_{{uo}}^{(1)}({\bf{k}}))}^{* }{C}_{\nu o}^{(0)}({\bf{k}})+{({C}_{\mu o}^{(0)}({\bf{k}}))}^{* }{U}_{{uo}}^{(1)}({\bf{k}}){C}_{\nu u}^{(0)}({\bf{k}})]\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&=\displaystyle \sum _{o,u}(f({\epsilon }_{o}({\bf{k}}))-f({\epsilon }_{u}({\bf{k}}))){({C}_{\mu u}^{(0)}({\bf{k}}){U}_{{uo}}^{(1)}({\bf{k}}))}^{* }{C}_{\nu o}^{(0)}({\bf{k}})\ \end{eqnarray} \tag{ 50 }$$

by switching the summation indices u, o in the second term and using ${U}_{{ou}}^{(1)}({\bf{k}})=-{({U}_{{uo}}^{(1)}({\bf{k}}))}^{* }$, as done already for equation (48). By this means, the response of the density matrix can be written as

$$\begin{eqnarray}&&{P}_{\mu m,\nu n}^{(1)}=\displaystyle \frac{1}{{N}_{k}}\displaystyle \sum _{{\bf{k}}}{e}^{-i{\bf{k}}({{\bf{R}}}_{n}-{{\bf{R}}}_{m})}\displaystyle \sum _{o,u}(f({\epsilon }_{o}({\bf{k}}))-f({\epsilon }_{u}({\bf{k}})))[{({C}_{\mu u}^{(0)}({\bf{k}}){U}_{{uo}}^{(1)}({\bf{k}}))}^{* }{C}_{\nu o}^{(0)}({\bf{k}})].\end{eqnarray} \tag{ 51 }$$

In practice, the evaluation of equation (51) can thus be restricted to the double sum over occupied o and unoccupied states u.

4.2. Response of the Kohn–Sham Hamiltonian

As discussed in section 3 for equation (37), the computation of $\langle {\psi }_{q}^{(0)}({\bf{k}})| {\hat{h}}_{{\rm{KS}}}^{(1)}| {\psi }_{p}^{(0)}({\bf{k}})\rangle$ is split into different steps: The matrices ${R}_{\mu \nu }^{(0)}({\bf{k}})$ and ${{\rm{\Omega }}}_{\mu \nu }^{(0)}({\bf{k}})$, which are defined in equations (34), (35) and which are required to calculate $\langle {\psi }_{q}^{(0)}({\bf{k}})| -{r}_{\gamma }| {\psi }_{p}^{(0)}({\bf{k}})\rangle$, are computed before the self-consistency loop, since they only depend on unperturbed properties. The definite u.c. integral appearing in equation (35) is integrated on a real-space grid using the formalisms described in [13, 14]. Conversely, the matrix ${V}_{\mu \nu }^{(1)}({\bf{k}})$, which is defined in equation (38) and which is required to compute $\langle {\psi }_{q}^{(0)}({\bf{k}})| {\hat{v}}_{{\rm{ext}}}^{(1)}(r)+{\hat{v}}_{{\rm{H}}}^{(1)}+{\hat{v}}_{{\rm{xc}}}^{(1)}| {\psi }_{p}^{(0)}({\bf{k}})\rangle$, explicitly depends on the response of the density ${n}^{(1)}({\bf{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 ${\hat{v}}_{{\rm{ext}}}^{(1)}({\bf{r}})$ and ${\hat{v}}_{{\rm{H}}}^{(1)}({\bf{r}})$ as well as the response of the exchange-correlation potential ${\hat{v}}_{{\rm{xc}}}^{(1)}({\bf{r}})$, as discussed below. The matrix elements ${V}_{\mu \nu }^{(1)}({\bf{k}})$ are then again obtained by performing the real-space u.c. integral appearing in equation (38) with the aforementioned techniques.

4.2.1. Response of the electrostatic potentials

As discussed in detail in [13, 16, 24], 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)}({\bf{r}})$ is decomposed into two terms

$$\begin{eqnarray}&&{n}^{(0)}({\bf{r}})=\displaystyle \sum _{{Im}}{n}_{{Im}}^{{\rm{free}}}({\bf{r}}-{{\bf{R}}}_{{Im}})+\delta n({\bf{r}}).\end{eqnarray} \tag{ 52 }$$

The first term describes the density associated with a superposition of 'free', i.e., completely isolated, spherically symmetric atoms ${n}_{{Im}}^{{\rm{free}}}({\bf{r}})$ located at the positions of the nuclei and of their periodic images ${{\bf{R}}}_{{Im}}$. The potentials of ${n}_{{Im}}^{{\rm{free}}}({\bf{r}})$ and $\delta n({\bf{r}})$ are computed independently and then reassembled to get the full electrostatic potential that enters the Kohn–Sham equations. For this purpose, $\delta n({\bf{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, see [13, 24, 31]. Given that the density response ${n}^{(1)}({\bf{r}})$ is also periodic in the perturbed case, see figure 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.

4.2.2. Response of the exchange-correlation potential

In semi-local approximations, the exchange-correlation potential ${\hat{v}}_{{\rm{xc}}}({\bf{r}})$ entering the Kohn–Sham Hamiltonian in equation (2) is given by

$$\begin{eqnarray}&&{\hat{v}}_{{\rm{xc}}}({\bf{r}})=\displaystyle \frac{\partial {E}_{{\rm{xc}}}[n({\bf{r}})]}{\partial n({\bf{r}})}.\end{eqnarray} \tag{ 53 }$$

Accordingly, its response ${\hat{v}}_{{\rm{xc}}}^{(1)}({\bf{r}})$ can be obtained via

$$\begin{eqnarray}&&{\hat{v}}_{{\rm{xc}}}^{(1)}({\bf{r}})=\int {\rm{d}}{\bf{r}}^{\prime} \displaystyle \frac{{\partial }^{2}{E}_{{\rm{xc}}}}{\partial n({\bf{r}})\partial n({\bf{r}}^{\prime} )}\displaystyle \frac{\partial n({\bf{r}}^{\prime} )}{\partial {e}_{\gamma }}=\int {\rm{d}}{\bf{r}}^{\prime} {f}_{{\rm{xc}}}({\bf{r}},{\bf{r}}^{\prime} )\ {n}^{(1)}({\bf{r}}^{\prime} )\end{eqnarray} \tag{ 54 }$$

by integrating over the exchange-correlation kernel ${f}_{{\rm{xc}}}({\bf{r}},{\bf{r}}^{\prime} )$, i.e., the second functional derivative of the exchange-correlation energy ${E}_{{\rm{xc}}}[n({\bf{r}})]$, and the density response ${n}^{(1)}({\bf{r}}^{\prime} )$. 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}_{{\rm{xc}}}({\bf{r}},{\bf{r}}^{\prime} )$. Additionally, many more exchange-correlation 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

$$\begin{eqnarray}&&{[{V}_{{\rm{HFX}}}^{(1)}]}_{\mu ,\nu }=-\displaystyle \frac{1}{2}\displaystyle \sum _{\lambda ,\sigma }{P}_{\lambda ,\sigma }^{(1)}({\chi }_{\mu }{\chi }_{\lambda }| {\chi }_{\nu }{\chi }_{\sigma })\ \end{eqnarray} \tag{ 55 }$$

needs to be added to the entries ${H}_{\mu ,\nu }^{(1)}$ of the Hamiltonian response matrix. Here, $({\chi }_{\mu }{\chi }_{\lambda }| {\chi }_{\nu }{\chi }_{\sigma })$ is the two-electron, four-index Coulomb integral defined and discussed in [15, 37, 38] and ${{\bf{P}}}^{(1)}$ is the first-order density matrix defined in equation (51).

4.3. Stable evaluation of the expansion matrix ${{\bf{U}}}^{(1)}({\bf{k}})$

To compute ${{\bf{U}}}^{(1)}({\bf{k}})$, one can in principle just evaluate equation (39) as discussed in the beginning of section 4. Thereby, only the entries

$$\begin{eqnarray}&&{U}_{{uo}}^{(1)}({\bf{k}})=\displaystyle \frac{1}{{\epsilon }_{o}^{(0)}({\bf{k}})-{\epsilon }_{u}^{(0)}({\bf{k}})}\left[\displaystyle \sum _{\mu \nu }{({C}_{\mu u}^{(0)}({\bf{k}}))}^{* }\ {V}_{\mu \nu }^{(1)}({\bf{k}})\ {C}_{\nu o}^{(0)}({\bf{k}})+{{\rm{\Omega }}}_{{uo}}({\bf{k}})\right]\end{eqnarray} \tag{ 56 }$$

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 ${{\bf{P}}}^{(1)}$, as shown and discussed for equation (51). Obviously, equation (56) becomes numerically unstable when quasi-degenerate eigenvalues are present close to the Fermi energy ${\epsilon }_{{\rm{F}}}$, since the denominator ${\epsilon }_{o}^{(0)}({\bf{k}})-{\epsilon }_{u}^{(0)}({\bf{k}})$ approaches zero in that case. In order to overcome this difficulty, we employ a technique originally proposed by de Gironcoli [6, 39] for DFPT-based lattice dynamics calculations in metals. We use a Fermi function with a small smearing width σ

$$\begin{eqnarray}&&\tilde{\theta }(\epsilon )=\displaystyle \frac{1}{1+{e}^{\epsilon /\sigma }}=\displaystyle \frac{1}{2}\left[1-\tanh \left(\displaystyle \frac{\epsilon }{2\sigma }\right)\right]\end{eqnarray} \tag{ 57 }$$

for the occupation numbers f(_o) and f(_u) appearing as difference in equation (51). We then pull this difference f(_o) − f(_u) inside the evaluation of ${U}_{{uo}}^{(1)}$ and re-write the problematic prefactor in equation (56) as

$$\begin{eqnarray}&&\displaystyle \frac{f({\epsilon }_{o})-f({\epsilon }_{u})}{{\epsilon }_{o}^{(0)}({\bf{k}})-{\epsilon }_{u}^{(0)}({\bf{k}})}\quad \longrightarrow \quad 2\displaystyle \frac{\tilde{\theta }({\epsilon }_{{\rm{F}}}-{\epsilon }_{u}({\bf{k}}))-\tilde{\theta }({\epsilon }_{{\rm{F}}}-{\epsilon }_{o}({\bf{k}}))}{{\epsilon }_{o}({\bf{k}})-{\epsilon }_{u}({\bf{k}})}\tilde{\theta }({\epsilon }_{o}({\bf{k}})-{\epsilon }_{u}({\bf{k}})),\end{eqnarray} \tag{ 58 }$$

as detailed in [39]. This has virtually no influence in the regime ${\epsilon }_{o}({\bf{k}})-{\epsilon }_{u}({\bf{k}})\gt \sigma$. For ${\epsilon }_{o}({\bf{k}})-{\epsilon }_{u}({\bf{k}})\ll \sigma$, we replace and evaluate the rewritten problematic factor by its analytic limit for ${\epsilon }_{u}\to {\epsilon }_{o}$:

$$\begin{eqnarray}&&\displaystyle \frac{\tilde{\theta }({\epsilon }_{{\rm{F}}}-{\epsilon }_{u}({\bf{k}}))-\tilde{\theta }({\epsilon }_{{\rm{F}}}-{\epsilon }_{o}({\bf{k}}))}{{\epsilon }_{o}({\bf{k}})-{\epsilon }_{u}({\bf{k}})}\tilde{\theta }({\epsilon }_{o}({\bf{k}})-{\epsilon }_{u}({\bf{k}}))\to -\displaystyle \frac{1}{2\sigma \left[1+\cosh \left(\tfrac{{\epsilon }_{{\rm{F}}}-{\epsilon }_{o}({\bf{k}})}{\sigma }\right)\right]}.\end{eqnarray} \tag{ 59 }$$

This expression is always finite and therefore numerically stable, even in the case of vanishingly small energy differences.

4.4. Evaluation of polarizabilities

In the last step, we evaluate the polarizability by rewriting equation (40):

$$\begin{eqnarray}&&{\alpha }_{\gamma \delta }=\displaystyle \frac{2}{{N}_{k}}\displaystyle \sum _{q,o,{\bf{k}}}f({\epsilon }_{o}({\bf{k}}))\ {\rm{Re}}\{{U}_{{qo}}^{(1,\delta )}({\bf{k}}){{\rm{\Omega }}}_{{oq}}^{(\gamma )}({\bf{k}})\}\end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{2}{{N}_{k}}\displaystyle \sum _{u,o,{\bf{k}}}\displaystyle \frac{f({\epsilon }_{o}({\bf{k}}))(2-f({\epsilon }_{u}({\bf{k}})))}{2}\ {\rm{Re}}\{{U}_{{uo}}^{(1,\delta )}({\bf{k}}){{\rm{\Omega }}}_{{ou}}^{(\gamma )}({\bf{k}})\}\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{2}{{N}_{k}}\displaystyle \sum _{u,o,{\bf{k}}}(f({\epsilon }_{o}({\bf{k}}))-f({\epsilon }_{u}({\bf{k}})))\ {\rm{Re}}\{{U}_{{uo}}^{(1,\delta )}({\bf{k}}){{\rm{\Omega }}}_{{ou}}^{(\gamma )}({\bf{k}})\}.\end{eqnarray} \tag{ 62 }$$

In the first step, the use of ${{\rm{\Omega }}}_{{qp}}({\bf{k}})={{\rm{\Omega }}}_{{pq}}^{* }({\bf{k}})$ reduces the summands to a real part ${\rm{Re()}}$, while in the second step the same procedure as used to obtain equation (48) is applied. In this way, the double sum can be limited in practice to only run over unoccupied u and occupied o states. The same strategy introduced for equation (58) and discussed in the previous section can be applied to deal with quasi-degenerate eigenvalues here. Again, the matrix ${{\bf{U}}}^{(1,\delta )}({\bf{k}})$ appearing in equation (62) is associated with a perturbation along the Cartesian axis δ, while the ${{\boldsymbol{\Omega }}}^{(\gamma )}({\bf{k}})$ matrix is associated to a perturbation along the axis γ.

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

5.1. 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 equation (4). As an example, we use ethylene (C₂H₄), 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⁻⁴ eV ${{\rm{\mathring{\rm A} }}}^{-1}$) 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 [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}_{0}^{3}$ (11%) at a tier 2 level. Only at the tier 3 level we get a maximum absolute (relative) error of 0.23 ${a}_{0}^{3}$ (1%). For semi-infinite systems, the dielectric constant, which is directly proportional to the polarizability as noted in equation (28), converges much faster with increasing basis set size, as also shown in figure 3 for bulk silicon. Even at a tier 1 level we essentially achieve convergence with an absolute (relative) error of 0.007 (0.05%).

Zoom In Zoom Out Reset image size

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 [40] for an example of basis sets adapted for polarizabilities [41], for hyperpolarizabilities, or [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₂H₄ 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 figure 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.

Zoom In Zoom Out Reset image size

Table 1.  Influence of H- and C-like 'ghost' atoms on the diagonal elements of the polarizability of C₂H₄, using light settings and LDA. Numbers in brackets indicate the mean polarizability.

Basis set	Polarizability				Time (s)
tier1	21.8	32.8	13.2	(22.6)	7.4
tier1 + 2 H-ghosts	22.3	33.2	18.6	(24.7)	12.9
tier1 + 2 C-ghosts	24.4	33.5	19.6	(25.8)	18.4
tier2	23.9	35.0	19.7	(26.2)	36.3

Finally, to study the sensitivity of the polarizability tensor on the k-point grid density in periodic systems, we also use silicon as example. Figure 5 displays the convergence behavior with respect to the size of the reciprocal-space ${\bf{k}}$-mesh in the primitive Brillouin zone. We observe a maximum absolute (relative) error of 0.12 (0.15%) when using 16 × 16 × 16 ${\bf{k}}$-points with respect to the converged result. This convergence behavior is comparable or slightly slower than what is observed for the total energy.

Zoom In Zoom Out Reset image size

5.2. 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 table 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.

Table 2.  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.

$| {\alpha }^{{\rm{FD}}}-{\alpha }^{{\rm{PT}}}|$	MAE (${a}_{0}^{3}$)	MAPE
Dimers	0.0004	0.0007%
Trimers	0.0002	0.001%
Molecules	0.0002	0.0008%

5.3. Extended systems: high-frequency dielectric constant

5.4. Performance and scaling of the implementation

To demonstrate the performance and scaling of our implementation, we show timings for the H(C₂H₄)_nH 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 u.c. from (C₂)₃₂ to (C₂)₅₁₂. 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 figure 6 (molecules), we find that the integration of the Hamiltonian response matrix ${{\bf{H}}}^{(1)}({\bf{k}})$ determines the computational time for small system sizes, i.e., for less than 200 atoms. Like for the update of the response density n⁽¹⁾, which involves similar numerical operations, we find a scaling of nearly O(N) for this step (see table 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 ${{\bf{P}}}^{(1)}$ becomes dominant, since it scales with O(N^2.5) in this regime. As discussed in section 4, the computation of ${{\bf{P}}}^{(1)}$ requires matrix multiplication operations, which traditionally scale O(N³). For bulk diamond we find a similar behavior and fit similar exponents, as shown in figure 6 and table 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.

Zoom In Zoom Out Reset image size

Table 4.  Fitted CPU time exponents α for the H(C₂H₄)_nH molecules (n = 8–128) and the periodic diamond discussed in the text. The fits were performed using the expression t = c N^α for the CPU time as function of the number of atoms N.

	H(C2H4)nH	Diamond
${n}^{(1)}({\bf{r}})$	1.1	1.4
${{\bf{H}}}^{(1)}({\bf{k}})$	1.4	1.5
${{\bf{P}}}^{(1)}$	2.5	2.6
Total	1.3	1.4

In summary, we find an overall scaling behavior that is always smaller than O(N²) 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 figure 7 for a u.c. 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 u.c. are in principle treatable within the relatively moderate CPU and memory resources offered by a single state-of-the-art workstation.

Zoom In Zoom Out Reset image size

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 [18, 53, 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 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 [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 u.c., 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 MD and compute the polarizabilities along these trajectories via PT. As explained in, e.g., [57], the polarizability tensor ${\boldsymbol{\alpha }}$ can be divided into an isotropic $\bar{\alpha }$, and an anisotropic component ${\boldsymbol{\beta }}$,

$$\begin{eqnarray}{\boldsymbol{\alpha }} & = & \bar{\alpha }{\boldsymbol{I}}+{\boldsymbol{\beta }}\\ \bar{\alpha } & = & \displaystyle \frac{1}{3}({\alpha }_{{xx}}+{\alpha }_{{yy}}+{\alpha }_{{zz}})\,,\,{\rm{Tr}}[{\boldsymbol{\beta }}]=0.\end{eqnarray} \tag{ 63 }$$

As shown in [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

$$\begin{eqnarray}I(\omega ) & = & {I}_{{\rm{iso}}}+\displaystyle \frac{4}{3}{I}_{{\rm{aniso}}}\\ {I}_{{\rm{iso}}} & = & \displaystyle \frac{N}{2\pi }{\displaystyle \int }_{-\infty }^{+\infty }{\rm{d}}{t}{{\rm{e}}}^{-{\rm{i}}\omega t}\langle \bar{\alpha }(0)\bar{\alpha }(t)\rangle \\ {I}_{{\rm{aniso}}} & = & \displaystyle \frac{N}{2\pi }{\displaystyle \int }_{-\infty }^{+\infty }{\rm{d}}{t}{{\rm{e}}}^{-{\rm{i}}\omega t}\displaystyle \frac{1}{10}\langle {\rm{Tr}}[{\boldsymbol{\beta }}(0)\cdot {\boldsymbol{\beta }}(t)]\rangle .\end{eqnarray} \tag{ 64 }$$

Here, N is again the number of atoms in the system. Furthermore, since the autocorrelation functions $\langle \cdot \rangle$ 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 equation (64) by $\beta {\hslash }\omega /(1-{{\rm{e}}}^{-\beta {\hslash }\omega })$, where $\beta =1/{k}_{{\rm{B}}}T$ [58]. Further frequency-dependent factors that multiply the vibrational Raman lineshapes are experiment-dependent [59–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].

We first analyze the sensitivity of the harmonic Raman intensities of the paracetamol molecule to the employed exchange-correlation functional. In figure 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.

Zoom In Zoom Out Reset image size

Conversely, the xc-functional chosen for the assessment of the PES has a large impact on the position of the peaks. In figure 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 blueshifts of the peak positions, especially above 1000 cm⁻¹, which means that these vibrational modes become more stiff [65, 66].

Zoom In Zoom Out Reset image size

In figures 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 quantum-nuclear effects in the dynamics [67, 68]. We have used 16 replicas of the system and the generalized-Langevin equation thermostat proposed in [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 redshift (with respect to the anharmonic classical spectrum) of about 70–100 cm⁻¹ 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 [68–70] 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⁻¹ solely due to nuclear quantum effects [68, 70].

Zoom In Zoom Out Reset image size

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 theoretical 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 blueshifts of up to 180 cm⁻¹ for frequencies above 3000 cm⁻¹ (see figure 9), while considering the quantum nature of the nuclei would redshift these frequencies by up to 150 cm⁻¹ 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⁻¹ region between calculated and experimental spectra. However, the calculated spectra are still blue-shifted with respect to experiment above 2500 cm⁻¹, even though lineshapes are well reproduced. The inclusion of nuclear quantum effects in the simulation would most likely solve this discrepancy, but the cost of such a simulation would be prohibitive at this point.

Zoom In Zoom Out Reset image size

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–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 [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.

In this paper, we derived and implemented a real-space formulation of PT 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²) 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 ${{\bf{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 blueshifts 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 PES, 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 the LDA functional 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] on the NOMAD Repository.

HS acknowledges Vivekanand Gobre and Wanzhen Liang for inspiring discussions. MR and NR acknowledge computer time from the Swiss National Supercomputing Centre (CSCS), under Project No. s786. We further acknowledge Volker Blum for his continued support during this project, Raul Laasner for providing the first version of the Pulay mixer we use, and Florian Knoop for critical proofreading. The project received funding from the Einstein foundation (project ETERNAL), the Deutsche Forschungsgemeinschaft (DFG) through SFB 951, and the European Unions Horizon 2020 research and innovation program under grant agreement no. 676580 with The Novel Materials Discovery (NOMAD) Laboratory, a European Center of Excellence. PR acknowledges financial support from the Academy of Finland through its Centres of Excellence Program (Project No. 251748 and 284621).

To validate our implementation for isolated systems, we compared the PT polarizabilities of 16 dimers (see table A1), 5 trimers (see table A2), and 11 molecules (see table A3) with those obtained by finite differences. In the latter case, the polarizability tensors were calculated using a finite, external electric field perturbation of  ±0.01 V Å⁻¹. All calculations were performed for fully relaxed geometries (remaining maximum force components smaller than 10⁻⁴ eV Å⁻¹) at the LDA level of theory using tier 2 basis sets and really tight defaults for all other numerical parameters such as integration grids. In all cases, we find that the observed deviations between the polarizabilities obtained via PT and via finite differences are orders of magnitude smaller than the polarizabilities themselves, as also substantiated by the respective MAEs and MAPEs given in tables A1–A3. Please note that even the largest observed absolute error (0.001 8 ${a}_{0}^{3}$ for HCN) corresponds to a very small relative error of only ∼0.008%.

Table A1.  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 errors (APE) with respect to finite difference calculations are given. Please note that these errors are several orders of magnitude smaller than the relevant digits in α_ii.

		PT (a03)	AE · 103 (a03)	APE · 103 (%)
Cl2	${\alpha }_{{xx}}$	24.123	0.07	0.29
	αyy	24.123	0.07	0.29
	αzz	41.309	0.10	0.24
ClF	αxx	15.927	0.05	0.31
	αyy	15.927	0.05	0.31
	αzz	22.292	0.06	0.27
CO	αxx	11.656	0.00	0.00
	αyy	11.656	0.00	0.00
	αzz	15.493	0.01	0.06
CS	αxx	22.239	0.03	0.13
	αyy	22.239	0.03	0.13
	αzz	37.652	0.07	0.19
F2	αxx	6.170	0.11	1.78
	αyy	6.170	0.11	1.78
	αzz	11.684	0.16	1.37
H2	αxx	3.902	0.01	0.26
	αyy	3.902	0.01	0.26
	αzz	7.532	0.03	0.40
HCl	αxx	16.815	0.00	0.00
	αyy	16.815	0.00	0.00
	αzz	18.868	0.05	0.26
HF	αxx	4.964	0.08	1.61
	αyy	4.964	0.08	1.61
	αzz	6.410	0.08	1.25
Li2	αxx	120.631	1.16	0.96
	αyy	120.631	1.16	0.96
	αzz	231.987	3.83	1.65
LiF	αxx	11.162	0.19	1.70
	αyy	11.162	0.19	1.70
	αzz	11.064	0.06	0.54
LiH	αxx	29.868	0.43	1.44
	αyy	29.868	0.43	1.44
	αzz	30.634	0.97	3.17
N2	αxx	9.923	0.04	0.40
	αyy	9.923	0.04	0.40
	αzz	15.033	0.14	0.93
Na2	αxx	121.132	0.58	0.48
	αyy	121.132	0.58	0.48
	αzz	283.915	7.22	2.54
NaCl	αxx	28.156	0.22	0.78
	αyy	28.156	0.22	0.78
	αzz	40.558	0.00	0.00
P2	αxx	34.724	0.02	0.06
	αyy	34.724	0.02	0.06
	αzz	67.280	0.09	0.13
SiO	αxx	24.570	0.40	1.63
	αyy	24.570	0.40	1.63
	αzz	34.021	0.13	0.38
Mean			0.41	0.77

Table A2.  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 errors (APE) with respect to finite difference calculations are given. Please note that these errors are several orders of magnitude smaller than the relevant digits in α_ii.

		PT (a03)	AE · 103 (${a}_{0}^{3}$)	APE · 103 (%)
CO2	αxx	12.041	0.14	1.16
	αyy	12.041	0.14	1.16
	αzz	26.559	0.23	0.87
H2O	αxx	8.576	0.06	0.70
	αyy	9.795	0.11	1.12
	αzz	9.191	0.03	0.33
HCN	αxx	13.101	0.08	0.61
	αyy	13.101	0.08	0.61
	αzz	23.102	1.88	8.14
SH2	αxx	23.169	0.03	0.13
	αyy	24.109	0.24	1.00
	αzz	24.052	0.04	0.17
SO2	αxx	18.869	0.03	0.16
	αyy	33.634	0.05	0.15
	αzz	22.710	0.01	0.04
Mean			0.21	1.09%

Table A3.  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 errors (APE) with respect to finite difference calculations are given. Please note that these errors are several orders of magnitude smaller than the relevant digits in α_ii.

		PT (${a}_{0}^{3}$)	AE · 103 (${a}_{0}^{3}$)	APE · 103 (%)
C2H2	αxx	16.323	0.09	0.55
	αyy	16.323	0.09	0.55
	αzz	31.802	0.22	0.69
C2H4	αxx	20.208	0.10	0.49
	αyy	24.666	0.34	1.38
	αzz	35.705	0.11	0.31
CH3Cl	αxx	26.327	0.02	0.08
	αyy	26.327	0.03	0.11
	αzz	35.999	0.10	0.28
CH4	αxx	16.974	0.62	3.65
	αyy	16.974	0.62	3.65
	αzz	16.974	0.62	3.65
H2CO	αxx	11.993	0.04	0.33
	αyy	18.333	0.09	0.49
	αzz	23.032	0.11	0.48
H2O2	αxx	13.596	0.12	0.88
	αyy	17.595	0.17	0.97
	αzz	12.362	0.13	1.05
N2H4	αxx	20.997	0.08	0.38
	αyy	25.882	0.07	0.27
	αzz	21.209	0.07	0.33
NH3	αxx	13.339	0.02	0.15
	αyy	13.339	0.00	0.00
	αzz	14.608	0.07	0.48
PH3	αxx	30.003	0.70	2.33
	αyy	30.003	0.70	2.33
	αzz	31.116	0.01	0.03
Si2H6	αxx	57.444	0.20	0.35
	αyy	57.444	0.23	0.40
	αzz	77.035	0.36	0.47
SiH4	αxx	31.967	0.14	0.44
	αyy	31.967	0.14	0.44
	αzz	31.967	0.14	0.44
Mean			0.20	0.86