Melting properties from ab initio free energy calculations: Iron at the Earth's inner-core boundary

We present a general scheme to accurately determine melting properties of materials from ab initio free energies. This scheme does not require prior fitting of system-specific interatomic potentials and is straightforward to implement. For the solid phase, ionic entropies are determined from the phonon quasiparticle spectra (PQS), which fully account for lattice anharmonicity in the thermodynamic limit. The resulting free energies are nearly identical (within 10 meV/atom) to those from the computationally more demanding thermodynamic integration (TI) approach. For the liquid phase, PQS are not directly applicable and free energies are determined via TI using the Weeks-Chandler-Andersen (WCA) gas as the reference system. The WCA is a simple, short-range, purely repulsive potential with established equation of states. As such, it is an ideal reference for ab initio TI of liquids. We apply this scheme to determine melting properties of hexagonal close-packed (hcp) iron at the Earth's inner core boundary $(P=330$ GPa), a subject of great significance in Earth sciences. The important influences of system size and pseudopotentials are carefully analyzed. The results (melting temperature equals $6170\pm 200$ K, latent heat $56\pm 2$ kJ/mol, Clapeyron slope $8.1\pm 0.2$ K/GPa) are consistent with experiments as well as previous calculations.

Figure 1

Phonon quasiparticles of hcp iron at $V=6.957{\AA }^3$/atom and $T=6400$ K. (a) Mode projected VDoS at $\mathrm{\Gamma }[\mathbf{q}$=(000)]. Vertical lines correspond to the frequencies of harmonic TO (left) and LO (right) phonons. Results from the 64 (180)-atom cell are denoted by red solid (dotted) lines, respectively. (b) Dispersions of phonon quasiparticles (red) and harmonic phonons (gray). Those from force constant matrices of 64 (180)-atom supercell are denoted by solid (dotted) lines, respectively. (c) Complete VDoS of phonon quasiparticles, where the effective harmonic force constant matrices are from MD involving 64 atoms (red solid line) and 180 atoms (red dotted line). Gray area denotes the VDoS of harmonic phonons, where the harmonic force constant matrix is from the 64-atom supercell. For both phonon quasiparticles and harmonic phonons, Fourier interpolations are performed on a dense $\mathbf{q}$-point mesh equivalent to a supercell containing 64000 atoms. Blue dashed line denotes the VDoS of MD involving 64 atoms without interpolation. (d) The Gibbs free energy per atom $\left(G\right)$ in the thermodynamic limit as a function of temperature $\left(T\right)$. Red circles (crosses) denote results from the PQS approach using simulation cells with $64\left(180\right)$ atoms. Blue solid (dotted) lines denote results from TI using $64\left(180\right)$ atoms. Inset: the differences with respect to $G$ from TI using 64 atoms. Uncertainties in $G$ from phonon quasiparticles are estimated from the variance of five independent simulations at 6400 K as 6 (4) meV/atom for the 64-atom (180-atom) simulations, respectively. Uncertainty associated with TI is estimated to be 3 meV/atom.

Figure 2

(a) The integrand of TI along the transition path at $V=6.957{\AA }^3$/atom and $T=6400$ K. (b)The $H/k_B{T}^2$ as a function of $T$. It is used in Eq. (17) to determine $G$ along the isobar. Inset: $H$ as a function of $T$. Colored (gray) symbols denote data associated with the solid (liquid) phase. Circles and solid lines (crosses and dotted lines) denote results from simulations involving $64\left(180\right)$ atoms.

Figure 3

Size dependence of free energies at $T=6400$ K.

Figure 4

TI of the liquid phase at $V=7.063{\AA }^3$/atom and $T=6400$ K. (a) Time evolution of $\mathrm{\Phi }$ (solid line) and ${\mathrm{\Phi }}_{\mathrm{WCA}}$ (dotted line) at ${\lambda }_2$ (red) and ${\lambda }_6$ (blue). As ${\mathrm{\Phi }}_{\mathrm{WCA}}$ is close to zero, for clarity it is multiplied by a factor of 10. (b) The integrand of TI after integral transformation. Inset: the original integrand as a function of $\lambda$.

Figure 5

Gibbs free energies of the hcp (circles) and liquid (crossed) phases from PAW-8 (dashed lines) and PAW-16 (solid lines) potentials.

Figure 6

Static equation of state of hcp iron determined from different PAW potentials. (a) $P-V$ curves, (b) pressure differences with respect to LAPW calculations.

Figure 7

Size dependence of melting temperature. Free energies are from PAW-8.