Efficient wide-range calculation of free energies in solids and liquids using reversible-scaling molecular dynamics

We elaborate a novel and efficient method to obtain multiphase Helmholtz free energies from molecular dynamics (MD) simulations over wide ranges of volume and temperature in materials that can be described by temperature-independent ion forces, with both higher accuracy and order-of-magnitude cost savings compared to direct thermodynamic-integration techniques. Our method leverages and significantly extends the technique of reversible-scaling molecular dynamics (RSMD) proposed by de Koning et al. [Phys. Rev. Lett. 83, 3973 (1999)], which allows a free-energy difference in a given phase at constant volume to be calculated as a function of temperature from a single MD simulation. In mechanically stable solid phases, our approach carefully combines quasiharmonic lattice dynamics at low temperatures with an accurate and fully isolated RSMD simulation of the anharmonic vibrational free energy at high temperatures to produce a seamless free energy from zero temperature to above melt along constant-volume isochores. In the liquid, we combine a unique calculation of the free energy along a high-temperature reference isotherm with isochoric RSMD simulations from that temperature to below melt. In metastable solid phases that are mechanically unstable at low temperature, we use two-phase MD melt simulations together with the liquid free energy to obtain the solid free energy along the solidus melt line and then perform isochoric RSMD simulations to temperatures above and below that point. While our free-energy method is general, we have specifically adapted it here to the case of metals in which the ion forces are well described by model generalized pseudopotential theory (MGPT) multi-ion interatomic potentials, and additive electron-thermal free-energy contributions can be included. Then using refined Ta6.8$x$ MGPT potentials, we have converged total free energies and their components to very high and unprecedented sub-milli-Rydberg (mRy) numerical accuracy in the stable-bcc, liquid, and metastable-fcc phases of tantalum for volumes ranging from up to 26% expansion to nearly twofold compression and for temperatures to 25 000 K. In turn, we have successfully used the free energies so obtained to calculate physically accurate thermodynamic properties and gain new insight into their behavior, including sensitive thermodynamic derivatives, bcc and fcc melt curves, and a multiphase equation of state for tantalum (Ta) over the same temperature range and for pressures as high as 600 GPa. We show that the anharmonic free-energy component in the bcc solid, although only 1–5 mRy in magnitude for Ta, can have a significant (15%–20%) effect on thermal expansivity, the Grüneisen parameter, and melt temperatures. We further show that the electron-thermal free-energy component can similarly impact the specific heat and thermal expansivity in both the solid and the liquid, while only minimally affecting (to$\le 3$%) the bcc and fcc melt curves.

Figure 1

Volume dependence of the MGPT potential parameters and pair potential for Ta6.8x. (a) Radial function exponent $p$ for $f\left(r\right)$ (right vertical scale) and $d$-band potential coefficients $v_{\alpha }$ (left vertical scale). (b) Two-ion pair potential $v_2$ for selected volumes, in x increments of 0.04 with $x={(\Omega /{\Omega }_0)}^{1/3}$.

Figure 2

Angular and volume dependence of the Ta6.8x multi-ion MGPT potentials $v_3$ (a) and $v_4$ (b) for specific near-neighbor configurations at selected volumes, in x increments of 0.04 with $x={(\Omega /{\Omega }_0)}^{1/3}$.

Figure 3

Effective pair potential $v_2^{\mathrm{eff}}$ in high-temperature liquid Ta at 25 000 K for $\Omega ={\Omega }_0$, as calculated from the present Ta 6.8x MGPT potentials $v_2$, $v_3$, and $v_4$ via Eq. (44) and as used to obtain an optimized $r^{-12}$ reference potential $v_{\mathrm{ref}}$ for the same conditions.

Figure 4

Volume and pressure dependence of the Debye temperature ${\Theta }_{\mathrm{D}}$ (a) and the quasiharmonic ion-Grüneisen parameter ${\gamma }_{\mathrm{ion}}^{\mathrm{qh}}$ (b) in bcc Ta, as calculated from the present MGPT Ta6.8x potentials, and as compared with corresponding results obtained from previous Ta3 [9] and Ta4 [9, 10, 11, 33] potentials and from experiment at ambient conditions [39, 40].

Figure 5

Anharmonic free energy in bcc Ta, as obtained from MGPT-RSMD simulations with the present Ta6.8x potentials. (a) Five selected isochores of the eight total calculated, with $x={(\Omega /{\Omega }_0)}^{1/3}$. Shown are both the raw simulation data (short dashed lines) and smooth fits (solid lines) to the data via Eq. (45). (b) Four selected high-temperature isotherms obtained via interpolation on the eight fitted isochores and displayed for $T<T_{\mathrm{L}}\left(\Omega \right)$.

Figure 6

Thermodynamic derivatives in bcc Ta at zero pressure, as calculated in the quasiharmonic (qh), quasiharmonic plus anharmonic $\left(\mathrm{qh}+\mathrm{ah}\right)$, and quasiharmonic plus anharmonic plus electron-thermal $\left(\mathrm{qh}+\mathrm{ah}+\mathrm{el}\right)$ limits of the present theory, and as compared to experiment. (a) Constant-pressure specific heat $C_{\mathrm{P}}$, with the measured data from Ref. [43] (solid circles) and Refs. [44–46] (solid square). (b) Thermal expansion coefficient or expansivity $\beta$, with the measured data from Ref. [47] (solid circles).

Figure 7

Temperature dependence of the Grüneisen parameter $\gamma ({\Omega }_0,T)$ in bcc Ta at ambient volume, as calculated in the quasiharmonic (qh), quasiharmonic plus anharmonic $\left(\mathrm{qh}+\mathrm{ah}\right)$, and quasiharmonic plus anharmonic plus electron-thermal $\left(\mathrm{qh}+\mathrm{ah}+\mathrm{el}\right)$ limits of the present theory, and as compared to experiment at room temperature [40].

Figure 8

Ion-thermal free energy $A_{\mathrm{ion}}^{\mathrm{liq}}$ in liquid Ta, as obtained with the present Ta6.8x MGPT potentials. (a) Reference isotherm $A_{\mathrm{ion}}^{\mathrm{liq}}(\Omega ,T_{\mathrm{ref}})$ for $T_{\mathrm{ref}}=25000$ K, as calculated both approximately from variational perturbation theory (VPT), using Eq. (30), and more accurately from thermodynamic integration (TI), using Eq. (21). (b) Four selected isochores of the eight total calculated via MGPT-RSMD simulation, in x increments of 0.08 with $x={(\Omega /{\Omega }_0)}^{1/3}$. Shown are both the raw simulation data (short dashed lines) and smooth fits (solid lines) to the data using Eq. (48).

Figure 9

Thermodynamic derivatives in liquid Ta at zero pressure, as calculated both without electron-thermal contributions (ion only) and with electron-thermal contributions (ion + el) in the present theory, and as compared to experimental isobaric expansion data from Ref. [45] (iso83) and Ref. [46] (iso86). (a) Constant-pressure specific heat $C_{\mathrm{P}}$, with the measured data at melt. (b) Thermal expansion coefficient or expansivity $\beta$.

Figure 10

High-pressure melting in Ta, as calculated from the present bcc-solid and liquid free energies. (a) Sensitivity of the melt curve to electron-thermal $\left(A_{\mathrm{el}}\right)$ and anharmonic $\left(A_{\mathrm{ah}}\right)$ free-energy contributions. (b) Full-theory melt curve and principal Hugoniot, with comparison to experimental isobaric [44, 45], DAC [49], and shock [50, 51] melt data, and to the calculated first-principles melt points of Taioli et al. [15].

Figure 11

High-pressure EOS properties in Ta, as obtained from the present bcc and liquid equations of state and compared to experiment. (a) Room-temperature (300 K) isotherm and principal Hugoniot calculated in the bcc phase (solid lines), and compared, respectively, to static DAC data (solid circles from Ref. [52] and open diamonds from Ref. [53]) and dynamic shock data (solid squares from Ref. [54]). (b) Principal Hugoniot calculated in both the bcc solid and the liquid, and compared to shock data (solid squares from Ref. [54] and solid diamonds from Ref. [55]).

Figure 12

High-pressure melt curve for the metastable-fcc structure in Ta, as calculated without electron-thermal free energy contributions $\left(A_{\mathrm{el}}=0\right)$. Solid square points denote input melt data from two-phase MGPT-MD simulations [12], while the solid line is the calculated melt curve from the output fcc free energy and the liquid free energy. The dashed lines represent approximate upper $\left(T_{\mathrm{c}}\right)$ and lower $\left(T_{\mathrm{ms}}\right)$ bounds on fcc mechanical stability.