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 $\ensuremath{\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\ifmmode\times\else\texttimes\fi{}$ 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 $\ifmmode\pm\else\textpm\fi{}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.