Understanding Anharmonicity in fcc Materials : From its Origin to ab initio Strategies beyond the Quasiharmonic Approximation

We derive the Gibbs energy including the anharmonic contribution due to phonon-phonon interactions for an extensive set of unary fcc metals (Al, Ag, Au, Cu, Ir, Ni, Pb, Pd, Pt, Rh) by combining densityfunctional-theory (DFT) calculations with efficient statistical sampling approaches. We show that the anharmonicity of the macroscopic system can be traced back to the anharmonicity in local pairwise interactions. Using this insight, we derive and benchmark a highly efficient approach which allows the computation of anharmonic contributions using a few T 1⁄4 0 K DFT calculations only.

For predicting real-world materials properties based on quantum mechanical ab initio approaches, it is crucial to take finite temperature effects accurately into account.This can be done by combining thermodynamic with statistical mechanics concepts employing specifically the partition function with ab initio computed total energies [1,2].A major challenge in doing so is the large number of configurations necessary to converge thermodynamic averages.Simple estimates show that brute force sampling requires on the order of 10 7 ab initio computed configurations to get the desired accuracy [3].Therefore, from the advent of ab initio guided materials design, sophisticated approaches have been suggested and employed to approximate free energies at finite temperatures with a minimum of sampling effort.
The presently most popular approach to compute thermodynamic properties from ab initio is the quasiharmonic (QH) approximation [4,5], for which the total energy of any atomic configuration is given by Here, N is the number of atoms, each Φ IJ is a 3 × 3 volume-dependent force constant matrix between atom I and J, and ũI is the displacement vector of atom I from its equilibrium position.Evaluating force constant matrices is nowadays almost routine within density-functional theory (DFT).Knowing them, the free energy can be straightforwardly obtained [6].
A drawback of the QH approach is that its accuracy is difficult to assess.Traditionally, anharmonic contributions beyond quasiharmonicity have been assumed to be small [7], most likely because accurate tools to compute them were missing.Only recent methodological advances, mainly based on thermodynamic integration techniques [8][9][10], allowed us to reduce the number of sampling configurations by 3-4 orders of magnitude providing, thereby, access to the numerically exact free-energy surface.Studies based on these approaches showed that anharmonic contributions may not only quantitatively but even qualitatively change thermodynamic quantities [11].Yet, the computations are still too expensive to be performed on a regular basis.
To avoid the computational cost of the numerically exact solution via statistical sampling, approaches have been developed to get an approximate description of the potential energy surfaces which go beyond Eq.(1).A main route in this direction is the treatment of anharmonic phononphonon interactions using second quantization [12][13][14], where the sum in Eq. (1) extends not only over pairwise interactions but includes higher-order terms.Other approaches use machine learning [15] or cluster expansion [16] concepts for approximating the potential energy surface or replace the T ¼ 0 K computed force constant matrix by one that optimally describes the system at a given temperature [17].In contrast to statistical sampling, a systematic convergence of these approaches proves challenging, making it difficult to assess their accuracy.
As a consequence of the difficulties arising when trying to go beyond Eq.( 1), systematic studies of anharmonicity for a wide range of materials are missing.In particular, it is unclear whether the huge configuration space can be compressed to a few physically intuitive collective coordinates.In this Letter we address this issue by systematically studying anharmonicity for a large set of fcc metals.Having the numerically exact free energy at hand allows us to confine the involved approximations to the choice of the DFT exchange-correlation functional and to unveil the importance of anharmonicity.Analyzing our ab initio molecular dynamics (MD) data, we gain a detailed understanding about the fundamental mechanism causing anharmonicity in these systems.This allows us to derive and test a highly efficient novel approach that enables accurate free-energy calculations over the whole temperature range to meV precision using a few sampling points derived at T ¼ 0 K.
We have calculated the total energies entering the electronic, quasiharmonic, and anharmonic free energies with the VASP code [18,19], projector augmented wave (PAW) potentials [20], and the local density approximation (LDA) as well as the generalized gradient approximation (GGA) [21][22][23].Spin polarization for Ni was included.Details regarding the electronic and quasiharmonic contributions were discussed previously [24].The anharmonic calculations were performed using the upsampled thermodynamic integration using Langevin dynamics (UP-TILD) method [3].Convergence errors due to DFT related parameters (e.g., k points, energy cutoff, supercell size) and statistical sampling have been carefully checked and are ensured to be below 1 meV=atom [25].
We focus first on the heat capacity which is a sensitive measure of the Gibbs energy surface and an experimentally extensively studied quantity.Figure 1(a) shows for the example of Ag the strong impact of anharmonicity: Employing the QH approximation (dashed curves for GGA and LDA) overestimates the final heat capacity by 1.4k B and 0.4k B (i.e.,37% and 10%) at the melting point.Including anharmonicity (solid curves) reduces the discrepancy between theory and experiment and, more importantly, the difference between LDA and GGA which has been identified in previous studies [3,26] as a reliable ab initio error measure (confidence interval) by more than a factor of 2. Figure 1(b) summarizes the GGA heat capacity contributions (electronic, quasiharmonic, anharmonic) for all studied elements and allows us, thus, to extract generic trends.All three contributions are roughly of the same order of magnitude, implying that neglecting anharmonicity (up to 1.2k B for LDA) or any of the other contributions would result in severe modifications.For Au, strong anharmonic together with many-body effects have been recently revealed to be crucial for an accurate thermodynamic description [26].Thus, even for high-symmetry closed packed fcc metals where anharmonic effects have been generally expected to be small, they are seen to be essential.The corresponding contributions to the Gibbs energy are shown in Fig. 1(c).We systematically observe anharmonic free energies to be about half the magnitude of the quasiharmonic free energies.The absolute anharmonic values of up to 25 meV (18 meV) for GGA (LDA) can easily modify phase transition temperatures by several hundred kelvin.
To understand the physical origin of the unexpectedly large anharmonicity, we carefully analyzed our MD data.We found that a particularly useful quantity for this purpose is the MD distribution of the first nearest-neighbor vectors d1NN .The corresponding QH and fully DFT-based distributions projected onto the (001) plane are shown in Fig. 2 direction, ẽL and ẽT .For the fully DFT-based MD, only the transversal direction symmetric, whereas a clearly nonsymmetric (anharmonic) behavior becomes evident in the longitudinal direction.Going towards the atom located at the origin gives rise to a sharp, almost planar edge beyond which the probability to find a neighboring atom drops essentially to zero, whereas on the opposite site a significant probability exists.
The observed distribution indicates that the functional dependence of the first nearest-neighbor distance is a good measure to represent anharmonicity.We therefore compute the Gauss-broadened first nearest-neighbor distribution function ρ 1NN ðdÞ ¼ P j δðj dj 1NN j − dÞ with j running over all first nearest-neighbor vectors of all MD time steps.In contrast to the fully symmetric QH distribution shown in Fig. 2(b) (black dashed line), the DFT distribution (red dashed line) shows a pronounced decrease (increase) of configurations smaller (larger) than the equilibrium distance.The redistribution of probabilities observed in the DFT MD case shown in Fig. 2 is fundamentally inaccessible in any harmonic description.It provides direct insight into the mechanism by which atomic correlation is able to lower the energy of the system and, thus, is the origin of anharmonicity: By shifting shorter bond configurations to on average longer bonds, the system very effectively avoids the strong Pauli repulsion while having to pay only the price for stretched bonds where the corresponding interaction is softer [33].This picture becomes even clearer when the atomic distribution function is used to construct the corresponding effective potential v eff ¼ −k B T ln ρ 1NN ðdÞ [see solid black and red lines in Fig. 2(b)].As expected, the full DFT effective potential (red line) is harder (softer) than the QH one for distances shorter (longer) than the equilibrium bond length.In fact, the thus constructed temperature-dependent nearestneighbor effective potential closely resembles the wellknown features of a Morse potential that is often used to describe the strength and anharmonicity of a T ¼ 0 K chemical bond.
Being able to relate finite temperature anharmonicity to the physically intuitive concept of local pairwise interactions, the question arises whether the corresponding interatomic potential can be directly derived from a few, ideally T ¼ 0 K calculations.Inspired by the distribution shown in Fig. 2(a), we use a local unitary transformation, U loc , which is an N × N block matrix with elements where Φ IJ corresponds to the quasiharmonic force constant matrix between atom I and J as in Eq. ( 1), ν IJ is the corresponding matrix of eigenvectors, and ẼIJ the vector of eigenvalues.Applying the 3 × 3 matrix U loc IJ on a Cartesian basis, A xyz ¼ ðẽ x ; ẽy ; ẽz Þ, results in a local basis with a longitudinal direction ẽL and two transversal vectors ẽT1 and ẽT2 [orthogonal to ẽL , compare Fig. 2(a)] for every atom pair.Using this new basis, the asymmetry in the DFT distribution function can be captured already from T ¼ 0 K calculations.For that purpose, we start with a perfect crystal with all atoms in equilibrium positions.We then displace atom I by a displacement u along each of the three principal vectors of A loc IJ and map the T ¼ 0 K force F0 K J and corresponding energy potential on the neighbor atom J as with a ¼ L; T1; T2.Crystal symmetries are employed to reduce the number of T ¼ 0 K calculations required to obtain the force and energy parametrizations to the irreducible ones.Having the parametrization of the forces and potentials available, the forces and energy for any atomic configuration, e.g., during a MD simulation, can be approximated.For example, the energy reads with n the number of nearest neighbors, u IJ ¼ j dIJ j − d eq IJ , u T1 ¼ dIJ • ẽT1 , u T2 ¼ dIJ • ẽT2 , dIJ the vector between atom I and J, and with d eq IJ the equilibrium distance of these atoms at T ¼ 0 K. Since the proposed method is based on probing the local anharmonic potential, we refer to formalism as the local anharmonic (LA) approximation.
A crucial difference between the LA approximation and the conventional QH direct-force-constant approach, where the force constants Φ IJ are formally restricted to the limit of infinitesimally small displacements (typical QH displacements are ≈0.01Å), is that significantly larger displacements are sampled during the LA parametrization [u in Eq. ( 4)] to accurately reproduce the anharmonic character.The magnitude of the displacements is dictated by the distribution function at the given temperature [see red dashed line in Figs.2(b) and 3(a)] and can reach values of > 1 Å, i.e., 2 orders of magnitude larger than those of the QH method.The resulting distance dependence, in particular, of the longitudinal component is no longer a linear function but highly anharmonic.
Figure 3(a) shows an example of the LA approximation.The longitudinal anharmonic forces (open red circles) derived from T ¼ 0 K DFT displacements are highly nonlinear.Trying to fit a linear function, as inherently done within the QH approximation (black solid line), is bound to introduce a large error.Rather, we generally find that a Morse potential provides an accurate description of the longitudinal DFT T ¼ 0 K forces (red line).A direct assessment of the quality of our LA approximation can be obtained by comparing finite-temperature forces extracted during a full DFT-MD run with forces computed for the identical configurations but using the LA approximation.Such a comparison is shown in Fig. 3(b) for Ag at its melting temperature (red dots).A dramatically improved correlation for the LA as compared to the QH approximation is observed.
To compute anharmonic free energies, the LA approximation can be used following two principle routes: (i) direct computation of LA anharmonic free energies using Eq. ( 5) by performing a thermodynamic integration from the quasiharmonic reference to the LA potential at negligible computational cost or (ii) calculation of the numerically exact DFT free energy using the LA potential as a reference system for thermodynamic integration to the DFT energy surface.Following the first route, we have computed for all studied elements the LA potential by a few (typically 4) DFT displacements at T ¼ 0 K , i.e., without any computationally expensive DFT MD.The free energy in the LA approximation is given by

where Fah
LA is the free energy obtained by thermodynamic integration from QH to the derived LA potential and ΔE LA DFT is an averaged difference between LA and DFT energies obtained from a few uncorrelated snapshots taken from a computationally inexpensive LA MD.We find ΔE LA DFT to converge for all studied elements within a few (< 5) uncorrelated structures to better than 1 meV.The accuracy achievable using route (i) in comparison to the fully ab initio computed anharmonic free energy [cf.Fig. 1(c)] is shown in row 3 of Table I, ΔF ah LAðLþT1Þ .Errors are on the order of a few meV for all elements highlighting the accuracy that can be achieved by applying the LA approach with a few (here < 10) T ¼ 0 K DFT calculations.
Following the second route (ii), we have computed numerically exact DFT free energies by thermodynamic integration with the LA approximation as reference.Since the LA forces faithfully reproduce the DFT forces [Fig.3(b)], the LA reference and the original system span very similar configuration spaces allowing us to obtain a very fast statistical convergence.The speed-up factors with respect to the previously applied QH reference given in Table I reveal that the computational effort can be reduced by about 2 orders of magnitude.For the free energy of Ag this resulted in a reduction from 4504 to less than 9 CPU hours [34] at the melting temperature and is now comparable to a QH calculation (∼6 CPU hours).This example illustrates that the LA method opens the path towards routine, numerically exact ab initio free energies.TABLE I. (Rows 1-3): Error in the free energy (in meV=atom and at T melt ) for the various approximations with respect to the numerically exact DFT free energy.(Row 1): QH approximation.(Row 2): LA approximation using a first nearest-neighbor longitudinal Morse parametrization only.(Row 3): As row 2, but augmented with a transversal function in ẽT1 direction corresponding to the second largest eigenvalue in Eq. ( 2).(Rows 4 and 5): Speed-up factors for the thermodynamic integration when replacing the QH by the LA reference.The speed-up is defined by n QH =n LA , where n QH and n LA are the MD steps necessary to reach a given standard error (typically 1 meV=atom) with a QH and LA reference.In conclusion, we computed statistically well-converged anharmonic free energies with DFT accuracy for a wide range of fcc metals.Using this computationally expensive data set we were able to identify the magnitude and origin anharmonicity in these materials.Based on this insight we extended the well-established quasiharmonic by the local anharmonic approximation, which allows us to obtain fully converged anharmonic energies at a fraction of the cost needed by previous approaches.We note that the LA approach may also be employed to study spectrally resolved anharmonic properties such as phonon renormalization and linewidths or thermal conductivity.First results on phonon linewidths are very promising and will be extended in future studies.The LA approach is not restricted to crystalline bulk systems but can be equally well used to compute the anharmonicity of alloys, defects, surfaces, etc. and promises a computational effort that is similar to quasiharmonic calculations.

FIG. 3 (
FIG. 3 (color online).Application of the local anharmonic (LA) approximation to Ag.(a) Longitudinal T ¼ 0 K DFT forces (open red circles) obtained by displacing an atom from its T ¼ 0 K equilibrium position, d eq;0 K1NN , towards and away from its nextnearest neighbor.The red solid line shows a fit according to the LA approximation and the black solid line using the QH approach.The necessary displacement range at T ¼ 0 K is dictated by the shown distance distribution at the melting temperature (red dashed line and shaded area).(b) Comparison between full DFT and LA forces (red dots), and between full DFT and QH forces (black dots) for configurations obtained from a fully DFT-based MD run at the melting temperature.