Random phase approximation up to the melting point: Impact of anharmonicity and nonlocal many-body effects on the thermodynamics of Au

Application of the generalized gradient corrected functional within standard density-functional theory results in a dramatic failure for Au, leading to divergent thermodynamic properties well below the melting point. By combining the upsampled thermodynamic integration using Langevin dynamics technique with the random phase approximation, we show that inclusion of nonlocal many-body effects leads to a stabilization and to an excellent agreement with experiment.

Recent years have seen significant progress in the development of ab initio-based databases (e.g., Ref. [1]).Such databases are valuable tools in modern materials design, given that the stored data are sufficiently accurate.Only then is full predictive power guaranteed.
Most present-day ab initio databases contain results from T = 0 K density-functional-theory (DFT) calculations using standard exchange-correlation functionals.This is a limitation in two respects: (a) First, materials development and application takes place at finite temperatures at which entropic contributions become important.In particular, anharmonic phonon-phonon excitations have been revealed to strongly modify materials properties, e.g., for dynamically unstable systems [2] or for defect formation [3].(b) Second, standard functionals are known to have intrinsic deficiencies, e.g., for elements with nearly full electron shells.Hybrid functionals provide an improvement by introducing a fraction of exact exchange [4][5][6], but the corresponding mixing parameter is not always well defined, especially for composite systems.Further improvement requires also the correlation energy to include nonlocal many-body effects.A possible route utilizes the random phase approximation (RPA) within the adiabatic-connection-fluctuation-dissipation theorem (ACFDT) [7].RPA has been successfully applied to various systems [8][9][10] and recent developments aim at improving the computational efficiency [11,12].However, despite the various developments and applications, RPA has been employed only at T = 0 K, whereas the impact of nonlocal many-body interactions at finite temperatures has remained unknown.
In this Rapid Communication, we compute the material properties of Au-a prototype closed shell element-using the RPA within the ACFDT formalism up to the melting point.This is made possible by advancing our previously introduced upsampled thermodynamic integration using Langevin dynamics (UP-TILD) technique [13] towards a combination of molecular dynamics simulations performed at the standard DFT level in conjunction with upsampled snapshots computed within RPA.We shall demonstrate that the inclusion of nonlocal many-body effects remedies the deficiencies introduced by standard DFT.
We start by investigating the thermodynamics of Au within the Perdew-Burke-Ernzerhof generalized gradient approximation (GGA-PBE [14]) in conjunction with the projectoraugmented-wave (PAW) formalism [15] as implemented in VASP [16,17] using the newly designed GW-PAW potentials [10].We concentrate on the heat capacity, which is a representative thermodynamic quantity.The state-of-the-art approach to calculate its temperature dependence includes electronic excitations in conjunction with the quasiharmonic approximation.We find that electronic excitations give a negligible contribution [18] due to the low electronic density of states of Au at the Fermi level.In contrast, the quasiharmonic approximation yields an unexpectedly large contribution, as revealed by the gray dashed line in Fig. 1 [19].In fact, the contribution is so large that the heat capacity diverges at 800 K, i.e., about 500 K below the experimental melting temperature.(Color online) Ambient pressure heat capacity of Au (k B = Boltzmann constant) utilizing various theory levels (qh = quasiharmonic, ah = anharmonic) and compared to experimental [20] and CALPHAD data [Scientific Group Thermodata Europe (SGTE) [21]] (T melt = 1337 K).
1098-0121/2015/91(20)/201103 (5) 201103-1 Published by the American Physical Society FIG. 2. (Color online) Linear thermal expansion calculated within the local-density approximation (LDA) and GGA-PBE (blue and orange lines; including anharmonicity), fully within RPA (green line; Au only), and using a mixed approach (red lines; discussed further below) employing the T = 0 K RPA curve in conjunction with GGA-PBE finite temperature excitations.For Au, the blue LDA curve lies fully behind the red curve.For Pt, the electronic entropy computed with GGA-PBE is included.Experimental values (solid circles) are from Ref. [23].
Recent work [22] showed that explicit anharmonicity can partly compensate the quasiharmonic heat capacity.As shown in Fig. 1 (orange line), anharmonicity has indeed a strong effect, but the divergence persists, shifting only to a higher temperature of 1100 K. Considering that relativistic effects are covered by the employed PAW potentials, the only remaining source for this unphysical divergence is the GGA-PBE approximation.
The failure of the GGA-PBE functional to describe properly the heat capacity of Au is highly unsatisfactory.We traced back this failure to an unphysical softness of bulk Au leading to an exponentially increasing thermal expansion well below the experimental melting point.This divergence is not observed in other metals (cf.orange lines in Fig. 2) where it is "hidden" since melting occurs already at temperatures below the onset of the divergent expansion.We can therefore conclude that GGA-PBE leads to an unphysical instability of fcc Au at temperatures that are more than 200 K below the melting point, in clear contrast to experimental observations (black circles in Figs. 1 and 2).
Knowing that RPA improves the T = 0 K properties of Au [10], we anticipate RPA to have an impact on finite temperature properties as well.Thus, an explicit finite temperature calculation would be desirable but is at present unfeasible due to the prohibitive computational workload.The strategy we follow here is to extend the UP-TILD method that was previously shown to largely reduce the number of configurations needed to get accurate thermodynamic averages [13].The original UP-TILD method was designed as a perturbation theory in the space of DFT convergence parameters (plane-wave cutoff, k-mesh sampling).Here, we extend the method in the following way.We utilize a molecular dynamics simulation performed within standard DFT, i.e., using GGA-PBE in the present case, to obtain a phase-space sampling at the desired volume and temperature.We then extract N uncorrelated snapshots (each corresponding to a 32 atomic supercell in this study) and calculate the UP-TILD energy E UP λ as an averaged difference between the RPA and GGA-PBE energies, E RPA λ,i and E GGA λ,i , as where i runs over the snapshots and where each energy is referenced with respect to the T = 0 K energies, E RPA 0K,sc and E GGA 0K,sc , calculated in the same supercell and with the same convergence parameters.Integrating the UP-TILD correction E UP λ over the coupling parameter λ, and adding the GGA-PBE quasiharmonic and anharmonic free energies, F GGA qh and F GGA ah , we obtain the RPA vibrational free energy as ( We observe that F UP is small [gray diamonds in Fig. 3(a)], that E UP λ quickly converges with the number of snapshots N , and that it is nearly independent of λ.The final RPA free energy is obtained as where E RPA 0K corresponds to a highly converged T = 0 K RPA energy-volume curve of a primitive cell as calculated previously [10].Further finite temperature excitations (e.g., electronic and magnetic) can be computed on the standard DFT level and added to Eq. ( 3).The introduced technique is not restricted to a combination of standard DFT with RPA.Any other higher level approach (e.g., quantum Monte Carlo [24]) can be coupled to standard DFT and used to obtain the accurate energies for the UP-TILD in Eq. (1) instead of the RPA.
In comparison to the original UP-TILD method, the present RPA UP-TILD technique requires additional convergence checks.It is known from previous work [9] that RPA is sensitive to the number of included unoccupied electronic states.To guarantee that our UP-TILD energy is converged, we have performed careful tests.Figure 3(b) shows that for the employed 32 atomic supercell the UP-TILD energy converges at 8000 electronic energy bands to about 1 meV/atom, which is sufficient for the present purpose.We have also checked the influence of the k-mesh sampling and found that a 2 × 2 × 2 k mesh leads to converged results.The energy cutoff was set to 450 eV for the wave function and 300 eV for the response 201103-2 function, and the frequency integration was done up to 800 eV at 16 frequencies.
Utilizing the RPA UP-TILD technique, we have computed a set of points at two temperatures and various volumes to obtain sufficient data for fitting an anharmonic free energy surface [cf.green squares in Fig. 3(a)].The resulting RPA heat capacity-implicitly including quasiharmonic and anharmonic excitations-is shown in Fig. 1 by the solid green curve.Remarkably, the divergence disappears and excellent agreement with experimental data is achieved.Consistently, the unphysical softening and exponential thermal expansion are removed within the shown temperature range (green line for Au in Fig. 2), i.e., they are shifted to temperatures above the experimental melting point.
A strength of the developed methodology is that the different contributions can be investigated separately.A corresponding analysis reveals that the difference between the GGA-PBE and RPA finite temperature free energies is rather small.This is quantified by the small values observed for the UP-TILD energy in Fig. 3(a) (gray diamonds).We stress that the absolute anharmonicity is significant, and it is only the difference between the GGA-PBE and RPA anharmonicity that is small [orange versus green solid lines in Fig. 3(a)].Considering this small difference, the question remains why RPA strongly impacts the thermodynamic properties at all.
Our analysis shows that this strong impact originates from differences in the T = 0 K energy-volume curves and corresponding values for the equilibrium lattice constant and bulk modulus.The GGA-PBE values (4.15 Å and 140 GPa) are substantially modified by the RPA (4.10 Å and 175 GPa) giving rise to a considerably stiffer bulk.This stiffness is the factor responsible for shifting the divergence in the thermodynamic quantities beyond the experimental melting point.
Utilizing this insight, we propose an efficient approach to compute accurate finite temperature properties beyond standard DFT.The approach combines T = 0 K RPA energies with entropic contributions from standard exchange-correlation functionals as GGA.Applying this approach, we obtained, with little computational effort, RPA finite temperature results for Ag and Pt, the elements showing the largest discrepancies within GGA-PBE except for Au [25].As shown by the red lines in Fig. 2, the mixed approach provides excellent agreement with experiment for all studied elements and with full RPA for Au.This indicates that nonlocal many-body effects as captured by RPA affect predominantly the absolute binding energies at T = 0 K, while thermal vibrations are unaffected to a good approximation, as evidenced by the small difference between GGA-PBE and RPA anharmonicity for Au [Fig.3(a)].
The insensitivity of high-temperature vibrations to nonlocal many-body effects is remarkable since recent studies [22] showed that close to melting the displacement of atoms from their equilibrium positions is significant and that the first-neighbor distance distribution is strongly affected by anharmonicity.Knowing that the GGA-PBE deficiency is not related to thermal vibrations but rather to the T = 0 K energy surface indicates a possible route for the development of accurate ab initio databases.T = 0 K computations are generally efficient since small supercells for single (or at most a few) atomic configurations are sufficient.It becomes therefore possible to employ computationally elaborate alternatives, x between LDA and RPA (blue lines) and between GGA-PBE and RPA (orange lines) for x = Pt (dotted lines), Ag (dashed lines), and Au (solid lines and light orange/blue shading) as a function of the lattice constant a x scaled with respect to the T = 0 K equilibrium lattice constant of RPA, a RPA x .Each curve is shifted to zero at a x /a RPA x = 1.Equilibrium lattice constants of LDA and GGA-PBE for each element are indicated at the bottom by the vertical lines and the arrows.(b), (c) T = 0 K binding pressure for GGA-PBE (orange solid line) and RPA (green solid line) and thermal phonon pressure (black solid lines) at four temperatures (300, 600, 900, 1200 K) for (b) Ag (melting point: 1235 K) and (c) Au (1337 K).For Au, the gray shaded area indicates the instability region (see main text).such as the RPA technique employed here, to supplement the databases.
Being able to assign the GGA-PBE error to a computationally inexpensive energy-volume curve allows us to further break down the responsible mechanism.To this end, we have analyzed the T = 0 K exchange-correlation energy of RPA (i.e., the sum of the exact exchange energy on the Kohn-Sham level and the RPA correlation energy) and the GGA-PBE exchange-correlation energy.We were interested in revealing differences among the investigated elements (Pt, Ag, Au) and in explaining why Au shows the largest deficiencies within GGA-PBE (orange curves in Fig. 2).Surprisingly, for all elements we find similar energy differences, E xc x , between the RPA and GGA-PBE exchange-correlation energy as quantified by the orange lines in Fig. 4(a) (dotted for Pt, dashed for Ag, solid for Au), showing E xc x as a function of the lattice constant.We observe E xc x to strongly decrease with increasing volume, resulting in a shift of the GGA-PBE equilibrium lattice constants towards larger values [cf. the bottom of Fig. 4(a)].Although E xc x is similar among the elements, the shifts in the lattice constant depend on the element with Ag showing the largest shift.This dependence can be understood by correlating it with the values for the bulk moduli of these elements (106, 175, and 267 GPa for Ag, Au, and Pt within RPA), i.e., the smaller the bulk modulus, the larger the shift in the lattice constant.
At this point, there seems to be a contradiction in the obtained results.As just reasoned, the T = 0 K equilibrium lattice constant is affected strongest for Ag due to the small bulk modulus.In contrast, finite temperature properties as the thermal expansion are affected strongest for Au when going from RPA to GGA-PBE.
The reason for this contradiction lies in the specific volume dependence of the thermal pressure of Au.In Figs.4(b) and 4(c) we use an unconventional, but for the present discussion convenient, representation to quantify the situation.The green and orange solid lines show a quantity that may be regarded as a (negative) T = 0 K binding pressure −P 0K that is obtained from the T = 0 K total energy as P 0K = −∂E 0K (V )/∂V .The T = 0 K binding pressure is the inner pressure due to the chemical bonds keeping the crystal together.The black solid lines show the thermal pressure P th (T ) = −∂[F (V ,T ) − E 0K (V )]/∂V , where F (V ,T ) is the free energy surface as a function of volume V and temperature T .The thermal pressure is the pressure built up by the vibrations and is opposite in sign to the binding pressure.It increases with temperature and leads to the expansion of the crystal.The crossing points between P 0K and P th (T ) give the equilibrium volume/lattice constant at a certain temperature, as shown for Ag and Au in Figs.4(b) and 4(c) for different temperatures.
The volume dependence of the thermal pressure has a similar qualitative dependence for Ag and Au (increasing with volume), but quantitatively the dependence is different in the relevant temperature interval (i.e., below melting).For Ag [Fig.4(b)], the upward curvature is modest and thus there is a well-defined crossing with the binding pressure curve −P 0K , regardless of whether the RPA (green) or the shifted GGA-PBE (orange) curve is considered.The situation is different for Au [Fig.4(c)], where P th strongly increases with volume, in particular, for higher temperatures.The increase is so strong that at temperatures above ∼1200 K, no crossing occurs.This is the critical temperature region where the thermal properties of Au start diverging, i.e., where the GGA-PBE Au bulk becomes unstable.This instability shows up also in the molecular dynamics simulations [gray shaded region in Fig. 4(c)] by frequent deviations from the perfect fcc structure (e.g., formation of defects), a phenomenon not occurring for Ag.We traced back these differences in the finite temperature behavior between Ag and Au to the Grüneisen parameter γ which measures the sensitivity of the phonon frequencies to volume changes.We find γ to be much larger for Au (γ = 3.3) than for Ag (γ = 2.4).Since a larger Grüneisen parameter indicates that the restoring forces-which realize the bonding pressure-quickly vanish with increasing volume, the crystal becomes prone to instabilities.
Finally, we study the performance of a second, popular exchange-correlation functional, the local-density approximation (LDA).Focusing on the thermodynamic properties in Figs. 1 and 2, we find a surprising agreement between the fully LDA derived data (blue lines) and RPA or experiment for all studied elements.Based on this agreement, one is tempted to conclude that LDA outperforms GGA-PBE.Such a conclusion is, however, too simplistic, as the following analysis shows.Similarly as for GGA-PBE, we plot in Fig. 4(a) the difference between the LDA and RPA exchange-correlation energy for Pt, Ag, and Au (blue lines).We observe that the differences are similar for the different elements and that they are even larger than for GGA-PBE (notice the steeper slope with volume for LDA).Consistently, the errors in the lattice constants are also larger than for GGA-PBE.
Considering this significant error in the T = 0 K energy, why does LDA show such remarkable finite temperature properties?The reason is that the T = 0 K error shifts the lattice constants towards smaller values.At these lattice constants the thermal pressure has a negligible volume dependence, and therefore the effect of the LDA error is small.Note that the absolute volume dependence of the thermal pressure is rather similar for LDA and GGA-PBE, and, in particular, LDA shows also a strongly increasing thermal pressure at larger volumes.These volumes are, however, not relevant for the thermodynamics of LDA due to the well-known overbinding of this functional.Thus, the agreement with the measured thermodynamic quantities is only due to a fortuitous error cancellation and cannot be expected in general.It would be interesting to extend our work in the future to an evaluation of other standard DFT functionals (e.g., PBEsol [26]).
In summary, we have extended the UP-TILD methodology, which now allows one to compute accurate finite temperature properties, including nonlocal many-body effects within RPA.The method is general and can be used also in combination with other higher level approaches.We have shown that the discrepancies introduced by standard functionals are removed and an excellent agreement with experiment is achieved.Identification of the T = 0 K energy as the main source of the error within standard DFT opens promising routes to a systematic increase of accuracy.

FIG. 3 .
FIG. 3. (Color online) (a) UP-TILD energy, F UP (gray diamonds), and anharmonic free energies computed using GGA-PBE, F GGA ah (orange circles), and RPA, F RPA ah = F GGA ah + F UP (green squares) at T melt = 1337 K.The vertical dashed (dotted) lines indicate the theoretical T = 0 K (experimental melting temperature) equilibrium lattice constant of RPA (green) and GGA-PBE (orange).For GGA-PBE only the T = 0 K value is shown.(b) Convergence of the UP-TILD energy as a function of the number of electronic bands included in the calculation.

FIG. 4 .
FIG. 4. (Color Difference in the T = 0 K exchangecorrelation energy E xcx between LDA and RPA (blue lines) and between GGA-PBE and RPA (orange lines) for x = Pt (dotted lines), Ag (dashed lines), and Au (solid lines and light orange/blue shading) as a function of the lattice constant a x scaled with respect to the T = 0 K equilibrium lattice constant of RPA, a RPA

Funding
by the European Research Council under the EU's 7th Framework Programme (FP7/2007-2013)/ERC Grant Agreement No. 290998 and by the German Federal Ministry for Education and Research under the BMBF NanoMatFutur Programme Grant No. 13N12972 is gratefully acknowledged.