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 high-pressure melt line and the behavior of the solid–liquid interface of iron under high-pressure 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 atom-in-jellium model;5 atomistic simulations using free energy calculations,6 thermodynamic integration,7 and two-phase simulations using density functional theory (DFT)8 or classical potentials;9,10 and a variety of experiments using shock-wave,11,12 diamond anvil cell (DAC),13 and high-energy laser compression coupled with in situ x-ray 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 point16,17 and also predicts the body-centered cubic (BCC) phase of iron to be stable in contradiction to experimental and other simulation results.18,19 Predicting the high-pressure melt line using DFT also remains challenging, mainly due to two factors: First, large simulation boxes are needed for typical two-phase 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, high-pressure conditions at the Earth's core necessitates the inclusion of semi-core 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 high-pressure 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 two-phase co-existence 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 re-calculate 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 five-steps:

  • 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 three-body interaction terms, the local and non-local 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, three-body, 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 two-phase simulation approach within the isoenthalpic-isobaric (NHP) ensemble at 300 GPa using the Large-scale 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 close-packed (HCP) unit cell (containing four iron atoms) by 5 × 12 × 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 GEAM-P8 shows that the difference in melting points for ( 0 1 ¯ 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 two-phase co-existence simulations performed in the current iteration to obtain the overall mean melting point, T m ( k + 1 ), and the overall standard deviation, Δ 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 ) T e ( k + 1 ) | < Δ 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).

TABLE I.

Melting temperatures for different GEAM potentials in the first and second iterations obtained from two-phase co-existence simulations in the NPH ensemble.

Melting point (K) Mean (K)
GEAM-P1 ( T m , 1 6113 ± 44  T ¯ m 1 = 6128 ± 47 
GEAM-P2 ( T m , 2 6146 ± 46   
GEAM-P3 ( T m , 3 6130 ± 44   
GEAM-P4 ( T m , 4 6124 ± 46   
GEAM-P5 ( T m , 5 6175 ± 49  T ¯ m 2 = 6142 ± 55 
GEAM-P6 ( T m , 6 6096 ± 49   
GEAM-P7 ( T m , 7 6154 ± 44   
GEAM-P8 ( T m , 8 6144 ± 45   
Melting point (K) Mean (K)
GEAM-P1 ( T m , 1 6113 ± 44  T ¯ m 1 = 6128 ± 47 
GEAM-P2 ( T m , 2 6146 ± 46   
GEAM-P3 ( T m , 3 6130 ± 44   
GEAM-P4 ( T m , 4 6124 ± 46   
GEAM-P5 ( T m , 5 6175 ± 49  T ¯ m 2 = 6142 ± 55 
GEAM-P6 ( T m , 6 6096 ± 49   
GEAM-P7 ( T m , 7 6154 ± 44   
GEAM-P8 ( 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 ¯ m ( 1 ) = 6128 ± 47 and T ¯ m ( 2 ) = 6142 ± 55 K, respectively.

The standard deviations associated with T ¯ m 1 and T ¯ 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 GEAM-P1 to GEAM-P4 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 GEAM-P5 to GEAM-P8 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 ¯ m ( 1 ) and T ¯ m ( 2 ) ( 6142 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 ¯ m ( 1 ) (electronic temperature is T ¯ e ( 1 ) = 6500 K) and T ¯ m ( 1 ) ( T ¯ 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 ¯ 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 0.0065 P 2. The predicted melting points are compared with recent experimental (static laser and in situ x-ray diffraction)14,15 and simulation results23 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.

FIG. 1.

Comparison of GEAM predicted melting points (triangle markers) with other works. The orange and the green dashed line denote the experimental measurement done by Anzellini et al.15 and Kraus et al.,14 respectively. The DFT calculation was done by Morard et al.8 The purple dashed line represents the classic MD simulation using the GAP potential.24 The brown and pink dashed lines show the predictions of the statistical moment method4 and the atom-in-jellium calculation.5 Statistical uncertainties associated with GEAM calculated melting points are also shown.

FIG. 1.

Comparison of GEAM predicted melting points (triangle markers) with other works. The orange and the green dashed line denote the experimental measurement done by Anzellini et al.15 and Kraus et al.,14 respectively. The DFT calculation was done by Morard et al.8 The purple dashed line represents the classic MD simulation using the GAP potential.24 The brown and pink dashed lines show the predictions of the statistical moment method4 and the atom-in-jellium calculation.5 Statistical uncertainties associated with GEAM calculated melting points are also shown.

Close modal

Figure 2 shows the phonon dispersion curves for BCC and HCP phases at 300 GPa by using the Phonopy package25 from GEAM-P8 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 GEAM-P8 and DFT highlights the reliability of the potential.

FIG. 2.

Phonon dispersion curves for BCC and HCP phases of iron at 300 GPa. The red and black curves correspond to results obtained from the GEAM-P8 potential and DFT calculations, respectively.

FIG. 2.

Phonon dispersion curves for BCC and HCP phases of iron at 300 GPa. The red and black curves correspond to results obtained from the GEAM-P8 potential and DFT calculations, respectively.

Close modal
TABLE II.

Lattice constants (Å), melting points, and enthalpy differences (eV/atom) between bulk solid and bulk liquid phases (for 50 selected structures at the respective melting points) of different phases at 300 GPa. The electronic smearing temperature for the DFT training set is 6128 K.

BCC FCC HCP
a (GEAM-P8)  2.401  3.013  2.139 
a (DFT)  2.401  3.013  2.141 
c (GEAM-P8)  ⋯  ⋯  3.439 
c (DFT)  ⋯  ⋯  3.433 
Tm (300 GPa)  5647 K  5858 K  6144 K 
Δ H l s(Tm, DFT)  0.36  0.48  0.50 
Δ H l s(Tm, GEAM-P8)  0.39  0.47  0.53 
BCC FCC HCP
a (GEAM-P8)  2.401  3.013  2.139 
a (DFT)  2.401  3.013  2.141 
c (GEAM-P8)  ⋯  ⋯  3.439 
c (DFT)  ⋯  ⋯  3.433 
Tm (300 GPa)  5647 K  5858 K  6144 K 
Δ H l s(Tm, DFT)  0.36  0.48  0.50 
Δ H l s(Tm, GEAM-P8)  0.39  0.47  0.53 

Since, the difference in enthalpies of bulk solid and liquid phases, Δ H l s, plays an important role in determining the melting point, we compare Δ H l s obtained from GEAM-P8 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. Δ H l s for the HCP phase obtained from GEAM-P8 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 GEAM-P8 and DFT predicted Δ H l s, it is reasonable to expect that the instability of BCC phase at 6144 K evident in GEAM-P8 is also correctly captured. In contrast, Δ H l 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 two-phase simulation methodologies on the melting point of iron at high-pressure 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 Δ T m  46 K and variations in the melting points due to changes in GEAM hyperparameters are less than 2 Δ 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 iron7,10,18,22 at high-pressure conditions still rely on either pair-interactions or do not include explicit three-body 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 high-pressure conditions, especially when good estimates of the melt line or large temperature and pressure ranges must be considered.

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. DE-AC52-07NA27344. This work was funded by the Laboratory Directed Research and Development (LDRD) Program at LLNL under project tracking Codes 21-ERD-005 and 23-SI-006. 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 LLNL-JRNL-855137.

The authors have no conflicts to disclose.

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).

The data that support the findings of this study are available within the article and its supplementary material.

1.
B. A.
Buffett
, “
Earth's core and the geodynamo
,”
Science
288
(
5473
),
2007
2012
(
2000
).
2.
J.
Badro
,
J.
Siebert
, and
F.
Nimmo
, “
An early geodynamo driven by exsolution of mantle components from Earth's core
,”
Nature
536
(
7616
),
326
328
(
2016
).
3.
J.
Li
,
Q.
Wu
,
J.
Li
,
T.
Xue
,
Y.
Tan
,
X.
Zhou
,
Y.
Zhang
,
Z.
Xiong
,
Z.
Gao
, and
T.
Sekine
, “
Shock melting curve of iron: A consensus on the temperature at the Earth's inner core boundary
,”
Geophys. Res. Lett.
47
(
15
),
e2020GL087758
, https://doi.org/10.1029/2020GL087758 (
2020
).
4.
T. D.
Cuong
,
N. Q.
Hoc
,
N. D.
Trung
,
N. T.
Thao
, and
A. D.
Phan
, “
Theoretical predictions of melting behaviors of hcp iron up to 4000 Gpa
,”
Phys. Rev. B
106
(
9
),
094103
(
2022
).
5.
D. C.
Swift
,
T.
Lockard
,
R. F.
Smith
,
C. J.
Wu
, and
L. X.
Benedict
, “
High pressure melt curve of iron from atom-in-jellium calculations
,”
Phys. Rev. Res.
2
(
2
),
023034
(
2020
).
6.
D.
Alfè
,
M.
Gillan
, and
G.
Price
, “
Complementary approaches to the ab initio calculation of melting properties
,”
J. Chem. Phys.
116
(
14
),
6170
6177
(
2002
).
7.
Y.
Sun
,
M. I.
Mendelev
,
F.
Zhang
,
X.
Liu
,
B.
Da
,
C.-Z.
Wang
,
R. M.
Wentzcovitch
, and
K.-M.
Ho
, “
Ab initio melting temperatures of bcc and hcp iron under the Earth's inner core condition
,”
Geophys. Res. Lett.
50
(
5
),
e2022GL102447
, https://doi.org/10.1029/2022GL102447 (
2023
).
8.
G.
Morard
,
J.
Bouchet
,
D.
Valencia
,
S.
Mazevet
, and
F.
Guyot
, “
The melting curve of iron at extreme pressures: Implications for planetary cores
,”
High Energy Density Phys.
7
(
3
),
141
144
(
2011
).
9.
A. B.
Belonoshko
,
R.
Ahuja
, and
B.
Johansson
, “
Quasi–Ab initio molecular dynamic study of Fe melting
,”
Phys. Rev. Lett.
84
(
16
),
3638
(
2000
).
10.
Y.
Sun
,
F.
Zhang
,
M. I.
Mendelev
,
R. M.
Wentzcovitch
, and
K.-M.
Ho
, “
Two-step nucleation of the Earth's inner core
,”
Proc. Natl. Acad. Sci. U. S. A.
119
(
2
),
e2113059119
(
2022
).
11.
J.
Brown
and
R.
McQueen
, “
Melting of iron under core conditions
,”
Geophys. Res. Lett.
7
(
7
),
533
536
, https://doi.org/10.1029/GL007i007p00533 (
1980
).
12.
C.
Yoo
,
N.
Holmes
,
M.
Ross
,
D.
Webb
, and
C.
Pike
, “
Shock temperatures and melting of iron at Earth core conditions
,”
Phys. Rev. Lett.
70
(
25
),
3931
(
1993
).
13.
R.
Sinmyo
,
K.
Hirose
, and
Y.
Ohishi
, “
Melting curve of iron to 290 GPa determined in a resistance-heated diamond-anvil cell
,”
Earth Planet. Sci. Lett.
510
,
45
52
(
2019
).
14.
R. G.
Kraus
,
R. J.
Hemley
,
S. J.
Ali
,
J. L.
Belof
,
L. X.
Benedict
,
J.
Bernier
,
D.
Braun
,
R.
Cohen
,
G. W.
Collins
,
F.
Coppari
et al, “
Measuring the melting curve of iron at super-Earth core conditions
,”
Science
375
(
6577
),
202
205
(
2022
).
15.
S.
Anzellini
,
A.
Dewaele
,
M.
Mezouar
,
P.
Loubeyre
, and
G.
Morard
, “
Melting of iron at Earth's inner core boundary based on fast x-ray diffraction
,”
Science
340
(
6131
),
464
466
(
2013
).
16.
H.
Sun
and
A.
Samanta
, “
Exploring structural transitions at grain boundaries in Nb using a generalized embedded atom interatomic potential
,”
Comput. Mater. Sci.
230
,
112497
(
2023
).
17.
B.
Sharma
,
Y. S.
Teh
,
B.
Sadigh
,
S.
Hamel
,
V.
Bulatov
, and
A.
Samanta
, “
Development of an interatomic potential for the W-Ta system
,”
Comput. Mater. Sci.
230
,
112486
(
2023
).
18.
A. B.
Belonoshko
,
T.
Lukinov
,
J.
Fu
,
J.
Zhao
,
S.
Davis
, and
S. I.
Simak
, “
Stabilization of body-centred cubic iron under inner-core conditions
,”
Nat. Geosci.
10
(
4
),
312
316
(
2017
).
19.
A. B.
Belonoshko
,
R.
Ahuja
, and
B.
Johansson
, “
Stability of the body-centred-cubic phase of iron in the Earth's inner core
,”
Nature
424
(
6952
),
1032
1034
(
2003
).
20.
A.
Samanta
and
J. L.
Belof
, “
The thermodynamics of a liquid-solid interface at extreme conditions: A model close-packed system up to 100 GPa
,”
J. Chem. Phys.
149
(
12
),
124703
(
2018
).
21.
A. P.
Thompson
,
H. M.
Aktulga
,
R.
Berger
,
D. S.
Bolintineanu
,
W. M.
Brown
,
P. S.
Crozier
,
P. J.
in't Veld
,
A.
Kohlmeyer
,
S. G.
Moore
,
T. D.
Nguyen
et al, “
LAMMPS—A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales
,”
Comput. Phys. Commun.
271
,
108171
(
2022
).
22.
F.
González-Cataldo
and
B.
Militzer
, “
Ab initio determination of iron melting at terapascal pressures and Super-Earths core crystallization
,”
Phys. Rev. Res.
5
,
033194
(
2023
).
23.
L.
Stixrude
, “
Melting in super-earths
,”
Philos. Trans. R. Soc., A
372
(
2014
),
20130076
(
2014
).
24.
Z.
Zhang
,
G.
Csányi
, and
D.
Alfè
, “
Partitioning of sulfur between solid and liquid iron under Earth's core conditions: Constraints from atomistic simulations with machine learning potentials
,”
Geochim. Cosmochim. Acta
291
,
5
18
(
2020
).
25.
A.
Togo
and
I.
Tanaka
, “
First principles phonon calculations in materials science
,”
Scr. Mater.
108
,
1
5
(
2015
).
26.
A. P.
Bartók
,
M. C.
Payne
,
R.
Kondor
, and
G.
Csányi
, “
Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons
,”
Phys. Rev. Lett.
104
(
13
),
136403
(
2010
).
27.
A. V.
Shapeev
, “
Moment tensor potentials: A class of systematically improvable interatomic potentials
,”
Multiscale Model. Simul.
14
(
3
),
1153
1173
(
2016
).

Supplementary Material