Reliable lattice dynamics from an efficient density functional

First principles predictions of lattice dynamics are of vital importance for a broad range of topics in materials science and condensed matter physics. The large-scale nature of lattice dynamics calculations and the desire to design novel materials with distinct properties demands that first principles predictions are accurate, transferable, efficient, and reliable for a wide variety of materials. In this work, we demonstrate that the recently constructed r2SCAN density functional meets this need for general systems by demonstrating phonon dispersions for typical systems with distinct chemical characteristics. The functional's performance opens a door for phonon-mediated materials discovery from first principles calculations.


I. INTRODUCTION
Each new age of human technology has been enabled by the discovery of new materials. From the knapped stone and simple metallurgy of history to the semiconductor revolution of recent decades, a new understanding of materials has expanded the horizons of possibility. This trend for materials to drive progress has not gone unnoticed and ever increasing effort is being devoted to using computational models to search the vast materials space for desirable new compounds.
Phonons are the quanta of lattice waves driven by the elementary thermal excitation of the atoms or molecules that make up a condensed matter system. Intuitively, long-wavelength phonons are perceived as sound. Phonons can interact with electronic structure and have a profound impact on a wide range of observed material phenomena, from thermal and electrical conductivity through to more exotic charge density waves and superconductivity, alongside their decisive role controlling the dynamic stability of materials. This position at the center of materials property design has driven prediction of phonon spectra to become an important aspect of materials space searches.
The connection between the vibrational frequency of a phonon, ω(k), and the wavevector, k, is known as the phonon dispersion. It can be measured experimentally by inelastic neutron or x-ray scattering. The phonon dispersion can also be predicted from theory using force constants calculated with computational models, though the cost of such calculations is generally high. This results in a simultaneous requirement for phonon calculations to efficiently scale for high-throughput workflows while maintaining sufficient accuracy to usefully guide experiments. Density functional theory (DFT) [1] using a semi-local exchange-correlation (XC) functional offers an appealing balance of these considerations and has become the workhorse computational method for high throughput materials discovery. Efficient evaluation of phonon * jning1@tulane.edu † jsun@tulane.edu spectra can be obtained from density functional calculations using density functional perturbation theory [2], or through direct displacements of the atoms [3]. Indeed, density functional methods have already proved effective tools for identifying new phonon phenomena, with the discovery of high/room temperature hydrogen-based superconductors as a prominent example [4,5].
High throughput calculations of thermodynamic properties [6] is vital for phase diagram predictions [7,8], and discovery of new meta-stable materials [9] places a particularly high demand for simultaneous accuracy, transferability, efficiency, and reliability. While the accuracy of a DFT calculation is largely determined by the accuracy of the chosen XC functional, the high computational demand of phonon calculations largely excludes expensive nonlocal XC functionals, like hybrid density functionals [10]. The current choice for phonon calculations remains conventional density functionals, including local density approximation (LDA) and the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA). While efficient and reasonably accuracy, one problem with these conventional density functionals is their transferability in the materials space where different compounds can have very different chemical bonds. Recent progress has shown that semi-local meta-GGAs can maintain this efficiency while being accurate for a wide variety of materials [9,[11][12][13][14][15][16][17][18][19][20][21][22], exemplified by the strongly constrained and appropriately normed (SCAN) meta-GGA [11,12]. Unfortunately, extensive use has shown that SCAN suffers numerical problems that are exaggerated in phonon calculations, making reliably obtaining accurate phonon spectra from SCAN calculations a challenging task. Here, we show that a revised version of SCAN which solves the numerical problems, called r 2 SCAN [23], delivers accurate, transferable, and reliable lattice dynamics. This is demonstrated in a selected set of materials that have different bonding characteristics. We then further explain the origins of such excellent performance of r2SCAN for phonon calculations.  [24], 1972 [25], and 1994 [26] for Si, of 1990 [27] for GaAs, from Ref. [28] for Fe, and from Ref. [29] for NiO. The second 3 experimental acoustic band data points along Γ-K direction for NiO were directly taken from their figure of Ref. [29].

II. RESULTS
The selected small test set includes four solids. Two are industrially important semiconductors with covalent or mixed covalent-ionic interactions: covalent Si in the diamond structure, and GaAs with zinc blende structure. One magnetic metal: body-centered-cubic (bcc) Fe. We also include a magnetic oxide with covalent-ionic bonding NiO. These solids have been widely studied experimentally and theoretically and their phonon dispersions have been accurately determined from experiments. In all cases we are interested in establishing both the accuracy and numerical stability of the phonon spectra calculated by different density functionals.

A. Si and GaAs
Figures 1 a) and b) compare calculated phonon dispersions with experimental results for Si and GaAs respectively. The aim of these calculations is to establish the relative accuracy of the functionals under ideal conditions with high-accuracy computational settings tuned to ensure well converged results for all functionals. We find that LDA predicts a relatively accurate spectrum for GaAs but underestimates lower frequency phonon bands in Si, while PBE underestimates phonon frequen-cies across the board. The r 2 SCAN meta-GGA shows the most consistent accuracy across both materials and at all energy ranges, closely matching the experimental data.
To establish the relative efficiency of r 2 SCAN compared to its parent SCAN functional, we repeat the calculations of Figures 1 a) and b) using the default VASP computational settings that better reflect a highthroughput workflow. The resulting phonon spectra are shown in Figures 2 a) and b) for Si and GaAs respectively. For Si of Figure 2 a), SCAN and r 2 SCAN show similar accuracy across much of the spectrum, though SCAN's error is significant for low frequency bands between the L and Γ points. For GaAs in Figure 2 b) however, the numerical problems of the SCAN functional are immediately apparent. Here spurious imaginary frequencies occur across the SCAN spectrum and the higher frequency bands show generally poor accuracy. Conversely, the r 2 SCAN functional remains well behaved under these cheaper settings, predicting an accurate and well converged spectrum for high and low frequencies. Note that imaginary frequencies are predicted by SCAN despite using a fully relaxed ionic structure that should be stable along all wave vectors. This prediction of spurious imaginary frequencies can be attributed to incomplete sampling of sharp oscillations in the SCAN XC potential as the ionic positions are displaced [23,30]. With special FIG. 2. Phonon dispersions of (a) Si and (b) GaAs calculated by SCAN and r2SCAN with default settings, compared with available experimental data from 1963 [24], 1972 [25], and 1994 [26] for Si, and from 1990 [27] for GaAs.
tuning of the parameters like Fourier transform grid density and atomic displacement size, SCAN can deliver accurate phonon dispersions for these two solids, as shown in the supplementary materials. Such tuning tricks just highlight the serious numerical problems of SCAN due to potential surface oscillations however. It is not guaranteed that these tuning tricks can solve SCAN's numerical problems for all solids.  It has commonly been observed that SCAN tends to overestimate magnetic moments and the magnetisation energy of simple magnetic metals like Fe, Co, and Ni [43][44][45][46][47][48][49]. As we include bcc Fe we have calculated the local magnetic moments for the transition metal atoms and present the results in Table I. We see that r 2 SCAN maintains the over-magnetisation of SCAN, as should be expected from r 2 SCAN's construction as a regularisation of SCAN. Figure 1 c) shows that the r 2 SCAN calculated phonon dispersions are not unduly degraded by this overmagnetisation however.
C. NiO Figure 1 d), compares the calculated and experimental phonon dispersions for NiO. The r 2 SCAN meta-GGA shows significant improvements over LDA and PBE for this material. In particular, the high-frequency optical bands from LDA and PBE are qualitatively wrong, while r 2 SCAN is reasonably accurate. Note we allow the crystal structure to be fully relaxed from the ideal AFM FCC structure, resulting in a small shrinkage of the lattice in the direction perpendicular to the ferromagnetic Ni planes. This symmetry breaking then leads to three optical bands. For this system we expect the main source of error to be self-interaction error, which causes the d orbital electron density to become too diffuse and fractionally occupied. In comparison with LDA and PBE, SCAN and r 2 SCAN reduce self-interaction errors and localize the d electrons around the Ni ion to a greater degree [16,17], stabilizing the magnetic moment as shown Table I. The self-interaction error of semi-local functionals can be remedied by including a Hubbard U term. The ad hoc nature of its parameterization limits predictive power however.
Polar bonds, such as those found in systems like NiO and GaAs, can cause the longitudinal optical and transverse optical (LO-TO) splitting in the experimentally observed phonon dispersion. A non-analytical correction to the phonon dispersion [50][51][52] based on the high-frequency dielectric constants ∞ and Born effective charge Z* must therefore be considered. Table I shows that r 2 SCAN delivers slightly better Z* in comparison with LDA and PBE. As ∞ is related to the response of electrons to the external electric field and can be strongly affected by the self-interaction error, r 2 SCAN significantly improves ∞ over LDA and PBE, although the discrepancy in ∞ from the experimental values is notable. In order to better illustrate the comparison of calculated force constants, we use the r 2 SCAN ∞ and Z* for the non-analytical term corrections [50][51][52] of the phonon dispersions of NiO and GaAs for all functionals.

III. DISCUSSION
Viewing Figure 1 as a whole, we can see a broad trend of accuracy across the different materials: PBE < LDA < r 2 SCAN. It is perhaps surprising that the inclusion of gradient information into the PBE GGA results in worse accuracy than the simpler LDA functional for the phonon dispersions. This effect can also be viewed from the other direction: it is surprising that LDA is as successful as it is for the phonon dispersion, particularly given its well known tendency to underestimate bond lengths and lattice constants. The reason for this can be found in the force constant dynamic matrix [2], where R I is the position of nucleus I, n R (r) and V R (r) are the ground state electron density and nuclear potential respectively with nuclei in positions R, and E N (R) is the Coulomb repulsion between the nuclei at positions R.
As the LDA bond lengths are too short, the second order derivative of the nuclear repulsion energy is overestimated (the final term of Eq. 1). This error is compensated however, by an overestimation of the linear response of electron density n R (r) to the nuclear distortion in the first term of Eq. 1. Since the first and final terms of Eq. 1) have opposite signs, their errors are favorably cancelled. The overestimation of the linear response of electron density is a consequence of the self-interaction errors intrinsic to semilocal density functionals, including LDA, PBE, and r 2 SCAN. PBE tends to overestimate bond lengths without correcting the linear response, so the favorable cancellation is lost. Like its parent functional r 2 SCAN improves both these aspects, giving accurate lattice constants [23] while simultaneously improving linear response characteristics [12] as demonstrated in Table I for lattice constants and ∞ . This results in a more accurate phonon spectrum with greater transferability across different classes of materials.
As previously mentioned, Figure 2 shows how r 2 SCAN improves on the SCAN functional by avoiding the numerical sensitivities that necessitate the expensive tuning of the fast Fourier transform grid. A full analysis of the origin of the numerical issues in SCAN, and their solution in r 2 SCAN, is presented in Ref. [23]. When calculating lattice dynamics from finite atomic displacements the smooth exchange-correlation potential of r 2 SCAN is well sampled by a coarse grid while the sharp oscillations of the SCAN potential are not [23,30]. This poor sampling results in slow and unpredictable convergence of the SCAN phonon spectrum with grid density, and the appearance of spurious imaginary frequencies.

IV. CONCLUSIONS
We have tested the performance of r 2 SCAN for calculating the phonon dispersions of typical systems relative to experimental data and other commonly used functionals (LDA and PBE). Our results for these systems suggest that r 2 SCAN can calculate accurate lattice dynamics for general systems with good transferability between different bonding characteristics. Across all the materials tested we find r 2 SCAN is either the best choice, or competitive with the best choice in the case of magnetic metals. While we find that SCAN can be accurate when Fourier transform grid and atomic displacement settings are tuned, however its poor numerical stability makes identifying the ideal parameters burdensome. Additionally, the necessary use of expensive Fourier transform grids prevents the SCAN functional being truly useful to high throughput studies. When default low-cost computational settings are used we find that SCAN predicts spurious imaginary bands. These problems are avoided in the r 2 SCAN functional which predicts accurate phonon spectra even from low-cost default parameters. While we find that while LDA and PBE can be quite accurate for some systems, they do not show the same generally transferable accuracy as r 2 SCAN does. With these inspiring findings, we strongly recommend r 2 SCAN to the community as an effective computational tool for future phonon dispersion studies.

V. METHODS
DFT [1] calculations with the LSDA, PBE [53], SCAN [11], and r 2 SCAN [23] XC functionals were performed using the Vienna Ab-initio Simulation Package (VASP) [54]. The projector-augmented wave (PAW) method was used to treat the core ion-electron interaction [55,56]. An energy cutoff of 600 eV was used to truncate the plane wave basis. A Γ-centered mesh with a spacing threshold of 0.15Å −1 was used for k-space sampling for unit cell relaxations of semiconducting systems Si, GaAs and NiO, and 0.1Å −1 for metallic Fe. For supercell atomic force calculations, only a single Γ point is used for semiconducting systems Si, GaAs and NiO, and a 2 × 2 × 2 k-point mesh for metallic Fe. A Gaussian smearing with 0.02 eV is used for semiconducting systems Si, GaAs and NiO, and Methfessel-Paxton smearing with 0.2 eV for Fe. For atomic force calculations, 3 × 3 × 3, 3 × 3 × 3, 5 × 5 × 5, and 4 × 4 × 2 supercells of the conventional unit cells (as shown in Table 1) are used for Si, GaAs, Fe and NiO, respectively. The ionic positions of all systems were relaxed for all functionals until the maximum ionic forces were below 1 meVÅ −1 . We used the Phonopy code [57] to obtain the harmonic force constants from VASP atomic force calculations within finite displacement method (0.015Å). For Figure 1, PREC = High; ENAUG = 2000 is specified for r 2 SCAN, while for Figure 2, we used the VASP officially recommended accurate defaults (PREC = Accurate) together with special tuned fast Fourier transform grid density for comparison. The full set of comparison is referred to the supplementary material. Table 1 and Figures 1-2 contain the detailed information for the comparison of numerical stability between SCAN and r 2 SCAN with respect to specific input settings, including the basic accuracy setting (the PREC tag in INCAR), the FFT-grid and the size of small displacements used in atomic force calculations, being tested on Si and GaAs. Figure 3 shows the phonon dispersions of GaAs and NiO calculated by LDA, PBE and r 2 SCAN with the non-analytical correction using Born effective charge and high-frequency dielectric constant from the underlying density functional. In contrast, the Figures 1(b) and (d) in the main text use the r 2 SCAN Born effective charge and high-frequency dielectric constants for all LDA, PBE and r 2 SCAN results.