Thermoelastic properties of bridgmanite using Deep Potential Molecular Dynamics

To generate the ab initio benchmark datasets for MgPv, we performed ab initio MD simulations employing the Vienna ab initio simulation package (VASP)[28]. We employed multiple exchange-correlation functionals to model the system, including the local density approximation (LDA)[29], the Perdew-Burke-Ernzerhof parametrization (PBE) generalized gradient approximation[30], the revised PBE for solids (PBEsol)[31], and the SCAN functional[26]. The projector augmented wave pseudopotential[32, 33] was employed, along with a plane-wave cutoff energy of up to 100 Ry. In AIMD simulations, we used 160-atom supercells of MgSiO₃ with $\Gamma$ k-point sampling. These DFT settings warranted good convergence of the calculated results[34].

Using the ab initio datasets of forces and energies, we developed deep potential (DP) models for BM employing the DeePMD-kit package[35, 36]. Two-body embedding with coordinates of the neighboring atoms (se_e2_a) was used for the descriptor. The embedding network was designed with a shape of (25, 50, 100), while the fitting network had a shape of (240, 240, 240). We used a cutoff radius of 6 Å and a smoothing parameter of 0.5 Å. The model was trained using the Adam optimizer[37] for $1\times 10^{6}$ training steps, with the learning rate exponentially decaying from $1\times 10^{-3}$ to $3.51\times 10^{-8}$ during the training process. The loss function $\mathcal{L}(p_{e},p_{f})$ is given by[35]:

where $p_{e}$ linearly decays from 1.00 to 0.02, while $p_{f}$ linearly increases from $1\times 10^{0}$ to $1\times 10^{3}$ throughout the training process.

We employed the DP-GEN concurrent learning scheme to create the reference dataset and generate the potential[38]. Initially, we randomly extracted 100 labeled configurations from 25 MD runs spanning various P-T ranges, covering 0-160 GPa and 300-4000 K to generate the initial potentials, shown in Fig. S1. We performed four DP-GEN iterations to explore the configuration space and ultimately achieve a potential that meets the desired accuracy for MD simulations. We trained four candidate DP potentials initialized with different random seeds in each iteration. The error estimator (model deviation) $\epsilon_{t}$ is determined based on the force disagreement between the candidate DPs[38, 39]. The expression for $\epsilon_{t}$ is:

where $F_{\omega,i}(\mathcal{R}_{t})$ represents the force on the $i$-th atom predicted by the $\omega$-th potential for the configuration $\mathcal{R}_{t}$. After the DP-GEN iterations, the final training dataset comprises 3,000 configurations annotated with ab initio force and energy information. We trained four potentials to reproduce the results of four exchange-correlation functionals: LDA, PBE, PBEsol, and SCAN. Consequently, we generated DP-LDA, DP-PBE, DP-PBEsol, and DP-SCAN potentials.

We used the LAMMPS package[35] to perform MD simulations. We obtained phonon dispersions using the finite-displacement method implemented in Phonopy[40] and PhonoLammps codes[41] using a supercell of $2\times 2\times 2$ with 160 atoms. The elastic properties were determined with DPMD using stress-strain relation[42, 43, 44] with 1,280 atoms and infinitesimal deformation.

We first assess the accuracy of our machine learning potential. A comprehensive comparison between the DP and DFT calculations is undertaken for energies and atomic forces. We quantified the root-mean-square error (RMSE) for the DP’s predictions on the validation sets. Fig. 1 shows the outcomes of this analysis for the DP-LDA calculation. The RMSE for potential energies amounted to approximately 0.82 $meV/atom$, while the RMSE for atomic force prediction errors was around 0.07 $eV/\AA$. Fig. S2-Fig. S4 validate the results of force and energy for DP-PBE, DP-PBEsol, and DP-SCAN, respectively. All these potentials show small RMSE compared to the DFT data, indicating the present deep potentials are well trained.

Fig. 2 shows the radial distribution functions of Bm at 4000 K and 120 GPa, obtained from both DFT-LDA and DP-LDA simulations. The DP potentials faithfully reproduce the $g(r)$ data derived from DFT calculations, indicating the DP’s capability to capture Bm’s structural and bonding properties accurately.

Since phonon dispersions are crucial in determining high-temperature properties, we also compare in Fig. 3 these dispersions at 0 GPa obtained with DP-LDA and DFT-LDA calculations. We see excellent agreement across almost all phonon branches, with only minor deviations in a few high-frequency optical branches along $\Gamma-Z$ path. The agreement between DPMD and DFT predictions, encompassing force, potential energy, $g(r)$, and phonon dispersion, attest to the robustness of our DP potential. Therefore, we can confidently proceed and compute the EOS and elastic coefficients, $\textbf{c}_{ij}$, at high-pressure and high-temperature conditions.

We compute the EOS across a broad range of P-T conditions. Fig. 4 compares EOSs obtained from our DP and previous LDA calculations[34]. Notably, the third-order Birch-Murnaghan EOS derived from DPMD results agrees exceedingly well with the DFT calculations. However, it is worth noting that both DPMD and DFT calculations tend to underestimate the pressure at a given volume when compared to experimental values[45, 46, 47, 48, 49, 50, 51, 52, 53, 54]. This discrepancy can be attributed mainly to the LDA functional employed in the present calculations. We also note that the DP calculations are based on classical MD simulations without taking into account the quantum effect of zero-point motion (ZPM). This ZPM effect was included in the previous LDA calculation[34], which manifests in the slight underestimation of DP-300 K compared with the DFT-300 K results. This ZPM effect is clearly reduced at 2000 K as the system approaches the classical limit.

Fig. 5 compares the EOS obtained with various functionals and experimental data at 300 K and 2000 K. The commonly used DP-PBE and DP-PBEsol functionals tend to overestimate the volume; DP-LDA slightly underestimates it, while DP-SCAN performs significantly better.

We first compare in Fig. 6 the elastic coefficients calculated using DPMD with those from previous PBE/GGA (DFT1, DFT2)[16, 55], LDA (DFT3) calculations[9], as well as experimental measurement[4, 5] at 300 K and 0 GPa; a comprehensive list of these coefficients is provided in Table S1. One can see DP-SCAN yields results more consistent with experimental data, particularly $\textbf{c}_{11}$, $\textbf{c}_{22}$, and $\textbf{c}_{33}$. This underscores DP-SCAN’s superior predictive capability for elastic properties. DP-LDA also agrees well with the experimental data, with minor discrepancies observed in $\textbf{c}_{44}$ and $\textbf{c}_{55}$. The DP-PBE and DP-PBEsol methods cannot accurately describe the elastic properties. Their deviation from the experimental data is approximately 12.5%. It is noteworthy that previous DFT calculations yield different elastic properties at room temperature, particularly for $\textbf{c}_{11}$ by Wehinger et al.[55]. Acoording to Wehinger et al.[55], this discrepancy was primarily attributed to previous calculations directly deriving results from acoustic phonon dispersions, bypassing numerical challenges arising from finite grid sampling in reciprocal space and a limited number of plane waves. In contrast, the present study with DP calculations enables large-scale MD simulations to directly investigate the stress-strain relationship and obtain accurate elastic behavior, which should ensure the reliability and robustness of our results.

Fig. 7 depicts the dependence of the elastic coefficients on P and T, revealing an almost linear increase with P and a linear decrease with T. These results are consistent with previous measurements[4] performed at ambient conditions. Furthermore, both the DP-SCAN and DP-LDA functionals agree well with experimental data at room temperature. Actually, for higher temperatures, our results closely match previous DFT calculations at 1500K[16]. At 3500K, our results remain consistent with earlier ones for $\textbf{c}_{22}$[16], while showing minor discrepancies in $\textbf{c}_{11}$, $\textbf{c}_{44}$, and $\textbf{c}_{66}$. Since the present results fully address anharmonic effects and overcome the simulation size limitation, it is expected to provide more accurate results than previous ones. In principle, it requires high P-T experimental data to validate these predictions, but such measurements are unavailable. We note that the comparison with experimental data is based on the isothermal elastic coefficients, as the difference between isothermal and adiabatic elastic coefficients is relatively small. The inclusion of adiabatic correction may increse the non-shear elastic coefficients, particularly at high temperatures.[13, 56, 57]

Fig. 8 illustrates the pressure-temperature dependence of the Voigt-Reuss-Hill (VRH)-averaged[58, 59] bulk modulus ($K_{S}$), shear modulus ($G$), as well as the compressive ($V_{p}$) and shear ($V_{S}$) wave velocities from DP-LDA and DP-SCAN. These properties are derived from the elastic coefficients ($\textbf{c}_{ij}$) and exhibit a consistently positive pressure dependence and negative temperature dependence. In Fig. 8 (a) and (b), we compare the predicted velocities and density of Bm with seismic values from the Preliminary Reference Earth Model (PREM)[60]. The density ($\rho$) is found to be relatively insensitive to temperatures, changing by only 1% with a 1000 K variation at the conditions present at the bottom of the lower mantle. Both $V_{S}^{PREM}$ and $V_{p}^{PREM}$ closely match the 2500 K to 3500 K isothermal shear and compressive wave velocities, in the mantle. Moreover, $V_{p}^{PREM}$ intersects two $V_{p}$ isotherms, at 2500 K and 3500 K, suggesting a temperature difference of less than 1000 K from top to bottom in the lower mantle if it were entirely composed of BM. Considering the shear and bulk moduli depicted in Fig. 8 (c) and (d), our results agree well with experimental data at 300K. Both $K_{S}^{PREM}$ and $G^{PREM}$ closely follow the corresponding 2500 K isotherms. Also, our shear modulus displays favorable agreement with previous DFT calculations[16], even at higher temperatures. Comparing (a) and (b), or (c) and (d), both DP-LDA and DP-SCAN can well describe these properties.