Melting of $\mathrm{MgSi}{\mathrm{O}}_3$ determined by machine learning potentials

Melting in the deep rocky portions of planets is important for understanding the thermal evolution of these bodies and the possible generation of magnetic fields in their underlying metallic cores. But the melting temperature of silicates is poorly constrained at the pressures expected in super-Earth exoplanets, the most abundant type of planets in the galaxy. Here, we propose an iterative learning scheme that combines enhanced sampling, feature selection, and deep learning, and develop a unified machine learning potential of ab initio quality valid over a wide pressure-temperature range to determine the melting temperature of $\mathrm{MgSi}{\mathrm{O}}_3$. The melting temperature of the high-pressure, post-perovskite phase, important for super-Earths, increases more rapidly with increasing pressure than that of the lower pressure perovskite phase, stable at the base of Earth's mantle. The volume of the liquid closely approaches that of the solid phases at the highest pressure of our study. Our computed triple point constrains the Clapeyron slope of the perovskite to post-perovskite transition, which we compare with observations of seismic reflectivity at the base of Earth's mantle to calibrate Earth's core heat flux.

Figure 1

Simulated structure factors (i.e., diffraction intensity) of bridgmanite (a) and liquid. (a) Schematic illustration of the crystal structure of bridgmanite, where the orthorhombic unit-cell is indicated by solid box. (b) and (c) Projected crystal structures of bridgmanite in the $x\text{−}y$ and $x\text{−}z$ planes, respectively. (d) and (e) display the simulated structure factors of bridgmanite and liquid from the three-dimensional and two-dimensional perspectives, respectively. The components of the collective variable, in which the descriptors, components of the collective variable $s_x$, are highlighted with black arrows. The corresponding planes of the two-dimensional descriptors are denoted in (b) and (c). The subscript $i$ (e.g., Mg, Si) indicates that only element $i$ is taken into accounts in the calculations.

Figure 2

Flowchart of the iterative training scheme. (1) We first build an initial dataset for liquid phase only and train a preliminary MLP as reported in our recent study [33]. Liquid, as a disordered phase, may encompass some of the local environments of the solid phases, and thus may serve as a good starting point for generating a unified MLP for both solid and liquid phases. (2) Multithermal-multibaric (MTMP) simulations are performed with the MLP using lammps interfaced with plumed 2. Here, plumed 2 is used to calculate the CVs and implement the enhanced sampling method. The target pressure and temperature ranges are very large in this study, making it difficult to cover in one MTMP simulation. We find that MTMP simulations with $\sim 20$ GPa and $\sim 2000$ K intervals yield good convergence and can sample the phase transition sufficiently. As a result, the target P/T ranges are divided into 20 GPa and 2000 K bins along the melting curves of bridgmanite and post-perovskite [4, 34]. We gradually update the P/T intervals with the iteration. (3) The resulting trajectories are saved every 500 time steps. The saved frames are converted to design matrices based on the smooth overlap of atomic positions descriptor [35]. We then perform principal component analysis on these design matrices and select the candidate configurations using the farthest point sampling technique [36, 37]. The size of candidate configurations N is large at first for a few iterations and gradually decreases at later iterations. Specifically, $\mathrm{N}=1000–2000$ for the first two iterations, $\mathrm{N}=50–100$ for the rest of iterations. (4) The selected frames were recalculated with DFT. The resulting energies and forces, are compared with the MLP predicted ones. For the first three iterations, the configurations with both energy difference $>15\mathrm{meV}/\mathrm{atom}$ and atomic forces difference $>0.5\mathrm{eV}/\AA$ are selected. For the rest of iterations, we relax the selection criteria to energy difference $>5\mathrm{meV}/\mathrm{atom}$ and atomic forces difference $>0.25\mathrm{eV}/\AA$. The size of the selected configurations in this step is M. We emphasize the selection criteria here is unlikely to be universal for all other systems but the principle that relaxing the selection criteria with iterations should apply. (5) The selected configurations will be combined with the initial dataset to train a new MLP. We re-iterate above steps until we cannot select frames in step (4) (i.e., $\mathrm{M}=0$) and all the target pressures and temperatures are covered by MTMP simulations in step 2).

Figure 3

Convergence tests of total energy (a) and pressure (b) with varying energy cutoff (ENCUT flag in VASP) for 32 $\mathrm{MgSi}{\mathrm{O}}_3$ bridgmanite at static condition. An energy cutoff of 800 eV is sufficient to obtain converged results for both energy and pressure.

Figure 4

Configurations explored by multithermal multibaric simulations at 40–60 GPa and 3000–5000 K. The energy of the system (160 atoms) is plotted against volume and color-coded by the value of collective variable (CV) defined in Eq. (4). Large and small CVs indicate a perovskitelike or liquidlike state, respectively. Snapshots of configurations are shown in the circles with atoms color-coded based on the orientationally targeted order parameters [32] building on the smooth overlap of atomic positions (SOAP) [35], with red indicating perovskitelike and blue indicating liquidlike local atomic environments. The yellow and green ellipses show the regions sampled by standard molecular dynamics simulations at the same P/T range for liquid and solid states, respectively.

Figure 5

Comparisons of energies (a), atomic forces (b), and stresses (c) between DFT and the machine learning potential (MLP) for all the test data at 2000 to 8000 K and pressures from $\sim 0$ to 220 GPa. 35 585 energies, 17 080 800 force components, and 21 3510 stress components are included in these comparisons. The red dashed lines are guides for perfect matches.

Figure 6

Comparisons of the total energy changes along molecular dynamics trajectories between the DFT (thick colored lines) and MLP potential (thin black lines) for $\mathrm{MgSi}{\mathrm{O}}_3$ bridgmanite (blue), post-perovskite (red), and liquid (green) at 140 GPa and 5000 K. The models used in this simulation contain 320 atoms, and none of the structures in the trajectories were included in the training set. The root-mean-square error of MLP is 4.2, 2.4, 7.1 meV/atom for perovskite, post-perovskite, and liquid, respectively.

Figure 7

Machine learning molecular dynamics simulations of the coexistence of $\mathrm{MgSi}{\mathrm{O}}_3$ bridgmanite and liquid at 140 GPa and 6100 K (upper panel) and 6000 K (lower panel). The simulation cell contains 600 $\mathrm{MgSi}{\mathrm{O}}_3$ formula units (3000 atoms). The simulation time step and the corresponding cell shape are also shown.

Figure 8

Machine learning molecular dynamics simulations of the coexistence of $\mathrm{MgSi}{\mathrm{O}}_3$ post-perovksite and liquid at 140 GPa and 5900 K (upper panel) and 6000 K (lower panel). The simulation cell contains 600 $\mathrm{MgSi}{\mathrm{O}}_3$ formula units (3000 atoms). The simulation time step and the corresponding cell shape are also shown.

Figure 9

Melting of $\mathrm{MgSi}{\mathrm{O}}_3$ bridgmanite (a), post-perovskite (b), and the phase boundary between bridgmanite and post-perovskite (c). Results from this study are shown solid blue circles for bridgmanite, red circles for post-perovskite, green circles for triple point and zero-K transition point. The uncertainties for melting temperatures are 50 K. The solid colored lines in (a) and (b) represent the Simon fit. The green dashed line in (a), (b), and (c) is the bridgmanite-post-perovskite phase transition boundary. Blue, green, and red shadings cover the stability fields of liquid, bridgmanite, and post-perovskite, respectively. (a) Previous experimental results on melting of $\mathrm{MgSi}{\mathrm{O}}_3$ bridgmanite are denoted by upward triangles [8], squares [51], leftward triangles [52], rightward triangles [53]. Experimental results of bridgmanite containing $\sim 10$ mol. % Fe are shown in solid dark upward triangles [8], downward triangles [9], and stars [54]. Prediction based on Lindemann law is shown dotted line [8]. Estimates based on atomistic modeling includes two-phase simulations based on classical potential with corrections [15] (dashed line), molecular dynamics simulations with empirical overheating correction [16] (loosely dashed line), and [17] (thin diamond), the integration of Clausius-Clapeyron equation by [34] (dashed-dotted line). (b) Previous results for the melting of post-perovskite include two shock compression experiments [11] (upward triangle) and [10] (downward triangle), the inferred melting curve using the Lindemann law [4] (dotted line), and two-phase simulations based on classical potential with corrections [15] (dashed line). (c) The results of subsolidus experiments of $\mathrm{MgSi}{\mathrm{O}}_3$ [55, 56] are shown with upward triangles, downward triangles, and squares denoting bridgmanite-only, post-perovskite-only, and bridgmanite-post-perovskite coexistence, respectively.

Figure 10

The volume (top) and the entropy (bottom) of melting from the bridgmanite (blue) and the post-perovskite (red) phases. For comparison we also show experimental values for melting from the low-pressure enstatite structure [57, 58], and theoretical results for the melting of monatomic systems interacting with inverse-power repulsion with the value of the power indicated ranging from 1 (one-component plasma) to infinite (hard-spheres) [59].