Accurate and precise ab initio anharmonic free-energy calculations for metallic crystals: Application to hcp Fe at high temperature and pressure

A framework for computing the anharmonic free energy (FE) of metallic crystals using Born-Oppenheimer ab initio molecular dynamics (AIMD) simulation, with unprecedented efficiency, is introduced and demonstrated for the hcp phase of iron at extreme conditions (up to $\approx 290$ GPa and 5000 K). The advances underlying this work are: (1) A recently introduced harmonically-mapped averaging temperature integration (HMA-TI) method that reduces the computational cost by order(s) of magnitude compared to the conventional TI approach. The TI path starts from zero Kelvin, where it assumes the behavior is given exactly by a harmonic treatment; this feature restricts application to systems that have no imaginary phonons in this limit. (2) A Langevin thermostat with the HMA-TI method that allows the use of a relatively large MD time step (4 fs, which is about eight times larger than the size needed for the Andersen thermostat) without loss of accuracy. (3) AIMD sampling is accelerated by using density functional theory (DFT) with a low-level parameter set, then the measured quantities of selected configurations are robustly reweighted to a higher level of DFT. This introduces a speedup of about 20–$30\times$ compared to directly simulating the accurate system. (4a) The temperature ($T$) dependence of the hcp equilibrium shape (i.e., $c/a$ axial ratio) is determined (including anharmonicity), with uncertainty less than $\pm 0.001$. (4b) Electronic excitation is included through Mermin's finite-temperature extension of the $T=0$ K DFT. A simple FE perturbation method is introduced to handle the difficulty associated with applying the TI method with a $T$-dependent geometry and (due to electronic excitation) potential-energy surface. (5) The FE in the thermodynamic limit is attained through extrapolation of only the (computationally inexpensive) quasiharmonic FE, because the anharmonic FE contribution has negligible finite-size effects. All methods introduced here do not affect the AIMD sampling—results are obtained through post-processing—so established AIMD codes can be employed without modification. Analytical formulas fitted to the results for the variation of the equilibrium $c/a$ ratio and FE components with $T$ are provided. Notably, effects of magnetic excitations are not included and may yet prove important to the overall FE; if so, it is plausible that such contributions can be added perturbatively to the FE values reported here. Notwithstanding these considerations, FE values are obtained with an estimated accuracy and precision of 2 meV/atom, suggesting that the capability to compute the phase diagram of iron at Earth's inner core conditions is within reach.

Figure 1

Convergence of the anharmonic energy (using the DFT-$0^{*}$ model) with respect to AIMD time-step size at 5000 K, using AT/LT thermostats and Conv/HMA averaging methods.

Figure 2

Variation of hcp $c/a$ axial ratio with temperature using lattice+quasiharmonic components (lat+qh) and including anharmonicity (full) using Eq. (16). The line is a second-order polynomial fit of the full $\beta \left[\alpha \left(T\right)-\alpha \left(0\right)\right]$ data, constrained to go to $\alpha \left(0\right)$ at 0 K (1.5925). The red points are experimental data for the indicated volumes per atom [76].

Figure 3

Variation of the Helmholtz FE contributions with temperature. (a) Lattice energy, relative to that at 0 K; the line is obtained from fitting $\beta \left[U^{\mathrm{lat}}\left(T\right)-U^{\mathrm{lat}}\left(0\right)\right]$ to force the data to go to zero at 0 K. (b) Quasiharmonic FE; the line is obtained from fitting $\beta \left[A^{\mathrm{qh}}\left(T;T_{\mathrm{el}}=T\right)-A^{\mathrm{qh}}\left(T;0\right)\right]$ to force the data to go to zero at 0 K. (c) Anharmonic FE of DFT-0 model (with $c/a=\widehat{\alpha }$) obtained from linear fitting of the HMA data (Fig. 4); dashed lines are the uncertainty limits. (d) Anharmonic FE correction between DFT-0 and DFT-el models; the line is obtained from fitting ${\beta }^2{A}_{0\to \mathrm{el}}^{\mathrm{ah}}$ to force both the data and slope to go to zero at 0 K. Uncertainties are based on $68\%$ confidence limits. See Supplemental Material for fitting details.

Figure 4

Anharmonic FE integrand of DFT-$0^{*}/0$ models, using the Conv/HMA averaging methods, with the temperature controlled using the Langevin thermostat. Values for $0^{*}$ are shifted slightly to the left to allow them to be distinguished from the DFT-0 values. Uncertainties are based on $68\%$ confidence limits.

Figure 5

Anharmonic finite-size effects of the DFT-0 model, using $c/a=\widehat{\alpha }$. The dashed lines delineate the $68\%$ confidence intervals.

Figure 6

Quasiharmonic finite-size effects of the DFT-el model at 5000 K, using $c/a={\alpha }_{\mathrm{eq}}\left(5000\mathrm{K}\right)=1.6164$. The $-k_{\mathrm{B}}TlnN/N$ term is added to reduce curvature in the data (see text for details).