An accurate calculation of the melting point of iron at various pressures in the Earth's core is important for understanding the core structure, geodynamo, and the Earth's history. Previous studies have assessed the melt line of iron at these extreme conditions using various experimental measurement techniques as well as both ab initio and classic molecular dynamics simulations. However, experimental measurements have uncertainties up to several hundred Kelvin, and inconsistencies remain among simulation results. In this work, we propose an iterative framework that couples density functional theory (DFT) calculations and molecular dynamics simulations performed using an ensemble of interatomic potentials to assess the effect of electronic temperature on the melting point. We systematically validate the potentials by comparing lattice constants and phonon dispersion curves at 0 K and enthalpy differences between liquid and HCP, FCC, BCC phases of iron close to the melt line at 300 GPa with DFT. Our results show that HCP iron melts at 6144 K (at 300 GPa), BCC phase is thermodynamically unstable, and FCC is metastable at this temperature. The melting points of FCC and BCC phases at 300 GPa are 5858 and 5647 K, respectively.
The Earth's core stores huge amounts of energy to power the Earth's geodynamo so that it creates a suitable living environment for all living things on its surface.^{1,2} Iron has consistently been regarded as the major component in the earth's core due to its abundant presence and its density.^{3} Earth's core consists of an inner solid core and an outer molten core, and the pressure at the boundary between the inner and outer cores is around 330 GPa. Thus, the melting point of iron at 330 GPa can be used to estimate the Earth's core temperature.
Unfortunately, studying the highpressure melt line and the behavior of the solid–liquid interface of iron under highpressure conditions is very challenging. In the past few decades, researchers have used multiple approaches to probe the behavior of iron in such extreme conditions, including theoretical analyses such as statistical moment method (SMM)^{4} and atominjellium model;^{5} atomistic simulations using free energy calculations,^{6} thermodynamic integration,^{7} and twophase simulations using density functional theory (DFT)^{8} or classical potentials;^{9,10} and a variety of experiments using shockwave,^{11,12} diamond anvil cell (DAC),^{13} and highenergy laser compression coupled with in situ xray diffraction.^{14} However, there are large discrepancies between the melting points obtained from these various methods. For instance, as illustrated in Refs. 7 and 15, melting points obtained from different measurement techniques differ by more than 1000 K.
On the simulation front, discrepancies in melting points obtained from different simulation methods stem from various factors. For example, the embedded atom method (EAM) potential employed in classical molecular dynamics (MD) simulations can model some ab initio properties well, but it may inadequately represent defect properties that are crucial for describing iron's behavior near its melting point^{16,17} and also predicts the bodycentered cubic (BCC) phase of iron to be stable in contradiction to experimental and other simulation results.^{18,19} Predicting the highpressure melt line using DFT also remains challenging, mainly due to two factors: First, large simulation boxes are needed for typical twophase simulations because the solid–liquid interface region can span seven or more atomic layers,^{20} causing the interface in a smaller supercell (accessible in DFT) to interact with periodic images. Second, highpressure conditions at the Earth's core necessitates the inclusion of semicore electron states in the valence electron system,^{6} significantly increasing the computational overhead. Another problem, which is often overlooked, is the effect of electronic temperature on the highpressure melting point of iron. The electronic temperature affects occupancies of different electronic states as well as total energies and forces on atoms. However, due to the high computational cost associated with the generation of a DFT dataset, the convergence (or iteration) with respect to electronic temperature is typically not analyzed. In addition, most simulations reported in the literature were performed using supercells containing only a few hundred atoms, and statistical (as well as systematic) uncertainties associated with calculated melting points are often not assessed.
To address these problems, we use an iterative computational framework to investigate the effect of electronic temperature on the equilibrium melt line of iron near the Earth's core pressure conditions and estimate statistical and systematic uncertainties associated with the predictions, converged for system size. This iterative framework involves developing an ensemble of interatomic potentials by using a DFT dataset generated at a fixed electronic temperature and estimating the mean melting temperature by using the twophase coexistence approach for each potential in the ensemble. In the next iteration, the mean melting point obtained in the previous iteration is used as the electronic temperature to recalculate the DFT energies and forces for the training structures and, then obtain another ensemble of potentials. This framework allows for the evaluation of the effect of electronic temperature on the converged equilibrium melting temperature and the uncertainty associated with it. We use the recently developed generalized embedded atom method (GEAM) potential framework, which has been used to predict melting temperatures of Ta, W, Nb, Mo, and V within 3% of the corresponding experimental results.^{16,17}
Our iterative approach to converge $ T e$ to $ T m$ involves on the following fivesteps:

For the first iteration, $ T e ( 1 ) = T m ( 0 )$ is set to an estimated temperature of 6500 K. For subsequent $ ( k + 1 )$th iterations, we regenerate the training dataset by performing DFT calculations such that the electronic temperature, $ T e ( k + 1 )$, is equal to the melting point obtained in the previous iteration, $ T m ( k )$. As we discuss below, two iterations are sufficient to obtain the converged results.

Next, we train an ensemble of GEAM potentials by varying the number of Gaussian basis functions used to represent the two and threebody interaction terms, the local and nonlocal atomic densities, and the polynomial order for the embedding energy term. Further details about the training set and DFT calculations are presented in the supplementary material.
At each iteration, we generate between 8 and 10 GEAM potentials and select the four potentials with the lowest mean absolute error (MAE) in total energies and forces (for structures in the training set) for subsequent melting point calculations. Representative hyperparameters, i.e., number of Gaussian functions and order of polynomials used in pair, threebody, embedding energy and nonlocal terms and the corresponding force and energy MAE values are shown in the supplementary material.

We use the potentials obtained in step (ii) to calculate the equilibrium melting points with MD, along with the standard deviations of instantaneous temperatures in each melting simulation.
To obtain the melting point, we use the twophase simulation approach within the isoenthalpicisobaric (NHP) ensemble at 300 GPa using the Largescale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software.^{21} For each simulation, we use an orthorhombic supercell containing 12 000 iron atoms initialized by replicating the conventional hexagonal closepacked (HCP) unit cell (containing four iron atoms) by $ 5 \xd7 12 \xd7 50$ along the three lattice vectors (solid–liquid interface is the (0001) plane of HCP). Similar system sizes and procedure (see the supplementary material) were used previously to compute melting points with close agreement with experimental measurements for W, Ta, Mo, Nb, and V.^{16,17} In addition, our analysis with GEAMP8 shows that the difference in melting points for $ ( 0 1 \xaf 01 )$ and (0001) interfaces at 300 GPa is less than 5 K, which is smaller than the statistical and systematic uncertainties discussed below.

Next, we accumulate data from the final 500 000 steps (1 ns) from each of the four twophase coexistence simulations performed in the current iteration to obtain the overall mean melting point, $ T m ( k + 1 )$, and the overall standard deviation, $ \Delta T m ( k + 1 )$. These values are summarized in the third column of Table I. The second column of Table I also shows melting points and corresponding standard deviations obtained from the four potentials.

If $  T m ( k + 1 ) \u2212 T e ( k + 1 )  < \Delta T m ( k + 1 )$, we terminate the computation, or else we set $ T e ( k + 1 ) = T m ( k + 1 )$ and repeat steps (i)–(v).
.  Melting point (K) .  Mean (K) . 

GEAMP1 ( $ T m , 1$)  6113 ± 44  $ T \xaf m 1 = 6128 \u2009 \xb1 \u2009 47$ 
GEAMP2 ( $ T m , 2$)  6146 ± 46  
GEAMP3 ( $ T m , 3$)  6130 ± 44  
GEAMP4 ( $ T m , 4$)  6124 ± 46  
GEAMP5 ( $ T m , 5$)  6175 ± 49  $ T \xaf m 2 = 6142 \u2009 \xb1 \u2009 55$ 
GEAMP6 ( $ T m , 6$)  6096 ± 49  
GEAMP7 ( $ T m , 7$)  6154 ± 44  
GEAMP8 ( $ T m , 8$)  6144 ± 45 
.  Melting point (K) .  Mean (K) . 

GEAMP1 ( $ T m , 1$)  6113 ± 44  $ T \xaf m 1 = 6128 \u2009 \xb1 \u2009 47$ 
GEAMP2 ( $ T m , 2$)  6146 ± 46  
GEAMP3 ( $ T m , 3$)  6130 ± 44  
GEAMP4 ( $ T m , 4$)  6124 ± 46  
GEAMP5 ( $ T m , 5$)  6175 ± 49  $ T \xaf m 2 = 6142 \u2009 \xb1 \u2009 55$ 
GEAMP6 ( $ T m , 6$)  6096 ± 49  
GEAMP7 ( $ T m , 7$)  6154 ± 44  
GEAMP8 ( $ T m , 8$)  6144 ± 45 
As shown in Table I, the statistical uncertainties shown in the second column for different GEAM potentials in the two iterations are similar, i.e., around 46 K. In addition, results in the third column show that the overall mean melting temperature and the standard deviation for iterations one and two are $ T \xaf m ( 1 ) = 6128 \xb1 47$ and $ T \xaf m ( 2 ) = 6142 \xb1 55$ K, respectively.
The standard deviations associated with $ T \xaf m 1$ and $ T \xaf m 2$ and the spread in the melting points obtained from different GEAM potentials in the two iterations quantify the systematic uncertainty arising due to the choice of hyperparameters, i.e., number of Gaussian basis functions used to model the different terms. Since the MAE for the predicted forces (for structures in the training set) for potentials GEAMP1 to GEAMP4 used in the first iteration are close (56.09–56.50 meV/Å), it is intuitive to expect the systematic errors arising due to these potentials to be small. In the second iteration, MAE for forces (57.89–60.09 meV/Å) and energies (2.49–2.84 meV/atom) for potentials GEAMP5 to GEAMP8 span a wider range than potentials in the first iteration, and correspondingly, the difference between the maximum and minimum melting points (79 K) is also higher. Therefore, variations due to different hyperparameters are similar to or smaller than statistical uncertainty in the predicted melting temperatures.
To compare the uncertainties discussed above with small supercells used in previous work in the literature,^{7,9,18,22} we perform MD simulations (using the NPT ensemble) using supercells containing 432 and 256 Fe atoms for BCC and FCC phases, respectively, at the respective melting points of these phases. The standard deviations of temperatures obtained from these simulations are 201 and 336 K for BCC and FCC phases, respectively, which are about 4–7 times larger than those listed in Table I demonstrating that large supercells used in our calculations significantly reduced statistical errors.
We note that the difference between $ T \xaf m ( 1 )$ and $ T \xaf m ( 2 )$ ( $ 6142 \u2212 6128 = 14$ K) is much less than statistical and systematic uncertainties, which suggests that the change in melting temperatures when the electronic temperature changes from 6500 to 6128 K is negligible, and it will not be significantly affected by using different hyperparameters of GEAM potentials. In addition, statistical uncertainties in the predicted melting points at $ T \xaf m ( 1 )$ (electronic temperature is $ T \xaf e ( 1 ) = 6500$ K) and $ T \xaf m ( 1 )$ ( $ T \xaf e ( 2 ) = 6128$ K) do not show any correlation with the electronic temperature. However, in other simulations where the temperature changes are larger, electronic temperature can have a larger effect and still should be considered.
Among the four potentials selected in the second iteration, i.e., at electronic temperature of 6128 K, the melting point of P8 is the closest to $ T \xaf m 2$, so we use P8 to extend the pressure range of our calculated melt line from 200 to 500 GPa, obtaining 5074 ± 35, 6144 ± 45, 7021 ± 48, and 7830 ± 58 K at 200, 300, 400, and 500 GPa, respectively. A polynomial fit of these melting points takes the following form $ T m = 2597.6 + 13.72 P \u2212 0.0065 P 2$. The predicted melting points are compared with recent experimental (static laser and in situ xray diffraction)^{14,15} and simulation results^{23} in Fig. 1. In addition, we find that BCC is unstable close to the melt line at 300 and 400 GPa, i.e., it spontaneously changes to HCP structure during NPT and NPH simulations performed at ∼6142 K, but the FCC phase is metastable close to the melt line. The stability of the HCP phase is consistent with recent experimental results.
Figure 2 shows the phonon dispersion curves for BCC and HCP phases at 300 GPa by using the Phonopy package^{25} from GEAMP8 and DFT using lattice constants shown in Table II. Note that structures with perturbations along high symmetry directions are not present in the training set. Therefore, the close agreement between GEAMP8 and DFT highlights the reliability of the potential.
.  BCC .  FCC .  HCP . 

a (GEAMP8)  2.401  3.013  2.139 
a (DFT)  2.401  3.013  2.141 
c (GEAMP8)  ⋯  ⋯  3.439 
c (DFT)  ⋯  ⋯  3.433 
T_{m} (300 GPa)  5647 K  5858 K  6144 K 
$ \Delta H l \u2212 s$(T_{m}, DFT)  0.36  0.48  0.50 
$ \Delta H l \u2212 s$(T_{m}, GEAMP8)  0.39  0.47  0.53 
.  BCC .  FCC .  HCP . 

a (GEAMP8)  2.401  3.013  2.139 
a (DFT)  2.401  3.013  2.141 
c (GEAMP8)  ⋯  ⋯  3.439 
c (DFT)  ⋯  ⋯  3.433 
T_{m} (300 GPa)  5647 K  5858 K  6144 K 
$ \Delta H l \u2212 s$(T_{m}, DFT)  0.36  0.48  0.50 
$ \Delta H l \u2212 s$(T_{m}, GEAMP8)  0.39  0.47  0.53 
Since, the difference in enthalpies of bulk solid and liquid phases, $ \Delta H l \u2212 s$, plays an important role in determining the melting point, we compare $ \Delta H l \u2212 s$ obtained from GEAMP8 for the HCP, FCC, and BCC phases of iron at different temperatures with calculations using DFT and the EAM potential developed by Sun et al.^{10} To this end, we randomly select 50 solid and 50 liquid configurations from two separate MD simulations (in NPT ensemble) at 6144 K and 300 GPa performed using bulk solid and bulk liquid structures. $ \Delta H l \u2212 s$ for the HCP phase obtained from GEAMP8 is 0.53 eV/atom compared to 0.50 eV/atom from DFT for the 50 selected bulk HCP and bulk liquid structures, which corresponds to only 8% deviation from the DFT value. The enthalpy differences for BCC and FCC iron at their respective melting points (at 300 GPa) are summarized in Table II. Due to the close agreement between GEAMP8 and DFT predicted $ \Delta H l \u2212 s$, it is reasonable to expect that the instability of BCC phase at 6144 K evident in GEAMP8 is also correctly captured. In contrast, $ \Delta H l \u2212 s$ calculated using the potential developed by Sun et al.^{10} for the same 50 HCP and 50 liquid structures is considerably lower at 0.17 eV/atom, and the enthalpy difference between 50 FCC and 50 liquid structures is only 0.19 eV/atom. This large discrepancy between the DFT value and the EAM potential points to the narrow range of applicability of that potential.
To conclude, we used an iterative computational approach to investigate the effect of electronic temperature and statistical uncertainty of twophase simulation methodologies on the melting point of iron at highpressure conditions similar to those in the Earth's core. Our analyses show that the statistical uncertainty (calculated from standard deviation of instantaneous temperatures) associated with melting points we calculated at 300 GPa is $ \Delta T m \u223c$ 46 K and variations in the melting points due to changes in GEAM hyperparameters are less than $ 2 \Delta T m$. These statistical uncertainties are much higher than the difference in melting temperatures (about 14 K at 300 GPa) due to change in the electronic temperature from 6500 to 6128 K.
Advanced interatomic potential models, such as Gaussian approximation potential (GAP)^{26} and moment tensor potential (MTP),^{27} developed in the past decade include explicit angular effects, but interatomic potentials developed for iron^{7,10,18,22} at highpressure conditions still rely on either pairinteractions or do not include explicit threebody interactions and, therefore, do not have the accuracy needed to probe equilibrium phase boundaries with high accuracy. The GEAM framework includes various higher order interactions, yet minimizes the number of free parameters and allows us to accurately model the energy surface from a minimal set of DFT calculated forces and energies. Our iterative framework provides an efficient and accurate avenue for future studies of equilibrium phase diagrams of metals and alloys at highpressure conditions, especially when good estimates of the melt line or large temperature and pressure ranges must be considered.
SUPPLEMENTARY MATERIAL
See the supplementary material for a description of the GEAM model, training dataset, DFT simulation settings, training procedure, and hyperparameters for the different GEAM potentials and enthalpies of solid and liquid phases.
This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DEAC5207NA27344. This work was funded by the Laboratory Directed Research and Development (LDRD) Program at LLNL under project tracking Codes 21ERD005 and 23SI006. The authors thank Babak Sadigh and Sebastien Hamel for stimulating discussions and Haoyuan Shi for help in setting up simulations. LLNL information management release number is LLNLJRNL855137.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Liming Zhao: Conceptualization (equal); Data curation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Vincenzo Lordi: Conceptualization (equal); Resources (equal); Writing – review & editing (equal). Amit Samanta: Conceptualization (equal); Methodology (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available within the article and its supplementary material.