Abstract
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, we achieve a 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 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.
Export citation and abstract BibTeX RIS
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1. 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 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.
2. 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 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
Here, n(r) is the electron density, the kinetic energy of non-interacting electrons, the external energy due to the electron-nuclear attraction, the Hartree energy, the exchange-correlation energy, and the repulsion energy of the nuclei. The ground-state electron density (and the associated ground-state total energy) is obtained by variationally minimizing equation (1) under the constraint that the number of electrons Ne is conserved. This yields the chemical potential of the electrons and the Kohn–Sham single particle equations
for the Kohn–Sham Hamiltonian . In equation (2), denotes the kinetic energy operator, the external potential, the Hartree potential, and 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
whereby the occupation numbers f(p) are chosen in such a way that the Ne/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 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 . In such a basis set, equation (2) becomes a generalized eigenvalue problem
Using the bra-ket notation for the inner product in Hilbert space, Hμν denotes the elements 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 and the associated total energy are determined by solving equation (6) iteratively, until self-consistency is achieved.
If an external electric field is applied to an isolated system, the KS Hamiltonian gains an additional term . If this electric field is a superposition of homogeneous electrical fields with strengths eγ aligned along the different cartesian axes γ, the additional term contributes
to the total energy functional in equation (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 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
used in the following for ground-state M(0) and response properties M(1) is always well-defined, 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:
where the response of the Hamiltonian operator is
Introducing the analogous expansions
for the single particle states ψp(eγ) and their eigenvalues p(eγ), rearranging the linear-order terms in the KS equation , and applying the normalization condition , yields the Sternheimer equation
as well as the condition
By multiplying equation (16) with from the left one obtains
To solve this equation numerically, we expand the response of the wave functions
in terms of the unperturbed states . Here, we chose for all p to fulfill equation (17) and hence obtain an algebraic expression for equation (18)
The expansion using the matrix 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 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 () of
The matrix , which fulfills , 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
in a density matrix formalism, i.e., without explicitly computing the response of the eigenvalues , of the wavefunction , or its coefficients , which is computationally advantageous. Using , 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 matrix, but rather compute the coefficients 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.
3. Generalization to periodic solids
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 are accounted for by summing over the lattice vectors , i.e., over linear combinations of the unit cell (u.c.) lattice vectors . Analogously, also the numeric atomic orbitals associated with such periodic images, e.g., associated with the periodic image m of nucleus , 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
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 , so that equation (5) becomes
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 -grid with Nk 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 in the solid, where is the electric displacement [25] and the polarization in the u.c. 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 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 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 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 is ill-defined in PBCs, since 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 . With that, one gets the well-known expression
that can be used to evaluate the non-diagonal matrix elements ()
Using equations (23) and (24) we obtain the representation
with
We can then recast the expectation value appearing in equation (18):
Here, the matrix elements can be directly evaluated as done in equation (25)
since they only depend on lattice periodic operators. Now, the matrix introduced in equation (21) is computed as
Similarly, the polarizability tensor components appearing in equation (10) can be rewritten as
using the matrix elements defined in equation (34). As explicitly highlighted in the notation, the matrix associated with a perturbation along the Cartesian axis δ has to be used in this case, whereas the matrix is associated with a perturbation along the Cartesian axis γ. Throughout the remainder of this work, the more general formulation in terms of Bloch-functions and wave vectors is used, since a simplification to finite systems is straightforward.
4. Details of the implementation
Our implementation closely follows the flowchart shown in figure 1: After a ground-state DFT calculation (see [13]) is completed, the matrix is computed. If 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 in a density matrix formalism (see section 4.1). As detailed in section 4.2, we then use to compute the remaining, individual ingredients that enter , i.e., the matrix elements defined in equation (38). The Sternheimer equation then provides a new matrix , 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 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.
Both the ground-state density n(0)(r) and the response of the density n(1)(r) are periodic, i.e., invariant against translations
by a lattice vector , 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.
Download figure:
Standard image High-resolution image4.1. Response of the electronic density
To numerically compute the electronic density in ground-state calculations [13], we use a density matrix formalism
which is obtained by inserting equations (24) and (23) into equation (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 in practice. Similarly, the response of the electronic density can thus be expressed as
using the response of the density matrix given by
In the sum over states o, we express in terms of via
whereby we split the sum over q into two separate sums over and u. In practice, these two sums can then be later evaluated by restricting the sum over 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 and one over . For the first one, we obtain
which vanishes due to , see equation (39). For the second double sum, we obtain
by switching the summation indices u, o in the second term and using , as done already for equation (48). By this means, the response of the density matrix can be written as
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 is split into different steps: The matrices and , which are defined in equations (34), (35) and which are required to calculate , 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 , which is defined in equation (38) and which is required to compute , explicitly depends on the response of the density 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 and as well as the response of the exchange-correlation potential , as discussed below. The matrix elements 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 is decomposed into two terms
The first term describes the density associated with a superposition of 'free', i.e., completely isolated, spherically symmetric atoms located at the positions of the nuclei and of their periodic images . The potentials of and are computed independently and then reassembled to get the full electrostatic potential that enters the Kohn–Sham equations. For this purpose, 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 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 entering the Kohn–Sham Hamiltonian in equation (2) is given by
Accordingly, its response can be obtained via
by integrating over the exchange-correlation kernel , i.e., the second functional derivative of the exchange-correlation energy , and the density response . 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 . 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
needs to be added to the entries of the Hamiltonian response matrix. Here, is the two-electron, four-index Coulomb integral defined and discussed in [15, 37, 38] and is the first-order density matrix defined in equation (51).
4.3. Stable evaluation of the expansion matrix
To compute , one can in principle just evaluate equation (39) as discussed in the beginning of section 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 , as shown and discussed for equation (51). Obviously, equation (56) becomes numerically unstable when quasi-degenerate eigenvalues are present close to the Fermi energy , since the denominator 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 σ
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 and re-write the problematic prefactor in equation (56) as
as detailed in [39]. This has virtually no influence in the regime . For , we replace and evaluate the rewritten problematic factor by its analytic limit for :
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):
In the first step, the use of reduces the summands to a real part , 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 appearing in equation (62) is associated with a perturbation along the Cartesian axis δ, while the matrix is associated to a perturbation along the axis γ.
5. Validation and performance
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 (C2H4), 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 [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 (11%) at a tier 2 level. Only at the tier 3 level we get a maximum absolute (relative) error of 0.23 (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%).
Download figure:
Standard image High-resolution imageThe 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 C2H4 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.
Download figure:
Standard image High-resolution imageTable 1. Influence of H- and C-like 'ghost' atoms on the diagonal elements of the polarizability of C2H4, 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 -mesh in the primitive Brillouin zone. We observe a maximum absolute (relative) error of 0.12 (0.15%) when using 16 × 16 × 16 -points with respect to the converged result. This convergence behavior is comparable or slightly slower than what is observed for the total energy.
Download figure:
Standard image High-resolution image5.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.
MAE () | MAPE | |
---|---|---|
Dimers | 0.0004 | 0.0007% |
Trimers | 0.0002 | 0.001% |
Molecules | 0.0002 | 0.0008% |
5.3. 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 LDA [42] and the GGA (PBE [34, 35]) and compared it with experimental and theoretical data compiled from literature [1, 22], see table 3. All calculations have been performed at the theoretical equilibrium lattice constant using 16 × 16 × 16 k-points in the primitive u.c. 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) [1, 22] 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 [22, 45].
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.
Exp. [46–50] | This work (all-electron) | NCPP | NCPP | PAW | PAW | ||
---|---|---|---|---|---|---|---|
1991 | 1996 | 2006 | 2016 | ||||
[22] | [1] | [43] | [44] | ||||
LDA | PBE | LDA | PBE | ||||
Si | 12.1 | 13.2 | 12.9 | 13.6 | — | 13.3 | 13.1 |
AlP | 7.5 | 8.4 | 8.2 | — | 8.2 | 8.3 | 8.1 |
AlAs | 8.2 | 9.5 | 9.5 | 9.2 | 9.3 | — | 9.5 |
AlSb | 10.24 | 11.7 | 11.9 | 12.2 | 11.4 | — | 12.1 |
GaP | 9.0 | 10.6 | 10.6 | — | 10.0 | — | 10.6 |
GaSb | 14.44 | 16.0 | 15.5 | 18.1 | 16.7 | — | — |
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 NCPPs: the agreement is generally better with literature results obtained using nonlinear core corrections [51, 52]. For instance, this can be observed for GaSb: the work of Dal Corso et al [1] made use of nonlinear core corrections, but not the earlier one of Giannozzi et al [22]. For the latter work, the use of a smaller -point grid may also be partially responsible for the observed deviations.
5.4. Performance and scaling of the implementation
To demonstrate the performance and scaling of our implementation, we show timings for the H(C2H4)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 (C2)32 to (C2)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 figure 6 (molecules), we find that the integration of the Hamiltonian response matrix determines the computational time for small 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 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 becomes dominant, since it scales with O(N2.5) in this regime. As discussed in section 4, the computation of requires matrix multiplication operations, which traditionally scale O(N3). 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.
Download figure:
Standard image High-resolution imageTable 4. Fitted CPU time exponents α for the H(C2H4)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 | |
---|---|---|
1.1 | 1.4 | |
1.4 | 1.5 | |
2.5 | 2.6 | |
Total | 1.3 | 1.4 |
In summary, we find an overall scaling behavior that is always smaller than O(N2) 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.
Download figure:
Standard image High-resolution image6. 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 [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 can be divided into an isotropic , and an anisotropic component ,
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
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 equation (64) by , where [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.
Download figure:
Standard image High-resolution imageConversely, 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−1, which means that these vibrational modes become more stiff [65, 66].
Download figure:
Standard image High-resolution imageIn 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−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 [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−1 solely due to nuclear quantum effects [68, 70].
Download figure:
Standard image High-resolution imageIn 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−1 for frequencies above 3000 cm−1 (see figure 9), while considering the quantum nature of the nuclei would redshift 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 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.
Download figure:
Standard image High-resolution imageIt 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.
7. Conclusions
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(N2) for both non-periodic (O(N1.3)) and periodic systems (O(N1.4)). In order to reduce the total time to O(N), more advanced algorithms [11, 12] for the evaluation of the density matrix response 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.
Acknowledgments
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).
: Appendix. Validation of the polarizability tensor for molecules
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 Å−1. All calculations were performed for fully relaxed geometries (remaining maximum force components smaller than 10−4 eV Å−1) 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 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 | 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 () | 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 () | AE · 103 () | 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 |