Probing the state of hydrogen in $\delta$-AlOOH at mantle conditions with machine learning potential

We employ the SCAN functional [54] to describe $\delta$’s potential energy (Born-Oppenheimer) surface and interatomic forces. Compared to LDA’s underestimation and PBE’s overestimation of the pressure, SCAN’s static compression curve is the closest to the 300 K measured one (see Fig. 1). The slight underestimation of the SCAN volume is eliminated when thermal effects at 300 K are included, as shown in the subsequent discussion. Quantum nuclear effects disregarded in the present study might spoil this good agreement [57], and this effect must be further inspected. We are interested in the temperature range of typical subducting tectonic plates (slabs), i.e., along “slab geotherms” [58, 59], where this phase is stable in the mantle, i.e., $\sim\!\mbox{1,000~{}K}<T<\mbox{3,000~{}K}$ [60, 21]. At these temperatures, classical MD should adequately describe the ionic dynamics in $\delta$.

Our SCAN-DP interatomic potential trained on just a few thousand reference configurations can reproduce SCAN-DFT’s forces and energies for configurations sampled from SCAN-BOMD simulations with 128 atoms, as illustrated in Fig. 2. The root-mean-square error (RMSE) [61] of the SCAN-DP predictions for these tests at various $(V,T)$ states are listed in Table SI. Configurations generated at higher temperatures produce larger RMSEs in general. At 3,000 K, the RMSEs for potential energy and force predictions are $\sim$2 meV/atom, and $\sim$0.12 eV/Å, respectively. Such accuracy is similar to previous benchmarks (e.g., [49]) in similar DP studies.

SCAN-DPMD also reproduces SCAN-BOMD’s structure and bonding properties at high-$P,T$, as illustrated by the comparison of SCAN-BOMD’s and SCAN-DPMD’s radial distribution functions, $g(r)$, in Fig. 3. The results are for 54 Å³ volume, corresponding to 9.5–28.3 GPa at 600–3,000 K obtained in 128-atom simulations. Plots detailing the contribution from each type of atomic pair at other $P,T$’s are shown in Fig. S1 ¹¹1See Supplemental Material online for Tables SI and SII, Figs. S1–S8..

The $g(r)$ at $V=$ 54 Å³ shown in Fig. 3 is particularly significant because this volume approximately corresponds to the critical volume of the experimentally observed H-bond disorder transition at 300 K [22]. At 600 K, the first peak at $\sim$1.1 Å corresponds to the ionic OH bond length. This peak is asymmetric at this volume and at all sampled temperatures due to asymmetric proton motion. At 600 K, $g(r)$ shows a near double-peak distribution of OH ionic bonds and H-bonds; the broad shoulder centered at $\sim$1.4 Å corresponds to the distribution of H-bond lengths. The increased overlap of the two peaks with increasing temperature suggests the onset of a disordered state, a precursor to H-bond symmetrization [23]. A similar OH bond-length distribution was observed in 300 K classical BOMD simulations with a quantum thermostat [63]. Some features on $g(r)$ disappear at elevated temperatures. For example, the 2.4 Å peak, a combination of Al-O and O-H contributions, and the 2.7 Å peak associated with O-O bonds (see Fig. S1) are split at 600 K but merge above 1,200 K. This phenomenon results from increased atomic vibration amplitudes, proton dynamic disorder, and diffusion onset at elevated temperatures. We will elaborate further on the diffusion process in Sec. II.3. Under pressure, the first peak at $\sim$1.1 Å and the shoulder at 1.4 Å evolve into a single symmetric peak due to H-bond disorder and finally symmetrization [23, 27] (see Fig. S1).

SCAN-DPMD with ab initio-level accuracy and improved efficiency will help predict the thermoelastic properties of hydrous solids, a property of first-order importance in geophysics [53, 64, 18]. Here, we compare its predictions for $\delta$’s high-temperature compressive behavior with ab initio SCAN-BOMD predictions and measurements in the 300–3,000 K and 20–115 GPa range. Fig. 4(a) shows excellent agreement between SCAN-DPMD’s high-temperature compression curves fitted to third-order Birch-Murnaghan equations of state (EoS) and the SCAN-BOMD’s predicted $P,V$ data points. Comparison with measurements [21, 22] in Fig. 4(b) shows a promising 3 GPa maximum error in pressure prediction at the same $V,T$s. The next section will present a more detailed analysis of the 300 K compression curve at lower pressure and compare it with measurements.

These good agreements between the SCAN-DFT and SCAN-DP predictions of forces, potential energies, $g(r)$, and high-$P,T$ EoSs suggest that the deep-learning potentials have reached satisfactory accuracy and are predictive at least in the 20–120 GPa pressure range and up to 3,000 K.

The 300 K compression curve is one of $\delta$’s best-determined properties. Experimentally, $\delta$ exhibits an anomalous change in compressive behavior at $\sim$10 GPa, which is usually attributed to a change in H-bond structure, most often with H-bond symmetrization [67, 68, 69, 24], but also with H-bond disorder [22, 23, 70]. However, it has been impossible to reproduce this compressive behavior using QHA-based DFT calculations. This is because (a) static ab initio calculations without multiconfiguration analysis do not account for H-bond disorder; (b) the HC phase atomic configuration is dynamically stable above $\sim$35 GPa only in static GGA/PBE calculations [25]. Strong anharmonicity produces unstable phonon modes at lower pressures, hindering the application of the QHA for high-temperature calculations [27, 25]. Using SCAN-DPMD, we optimize and equilibrate the cell shape at several pressures to obtain the static and the 300 K compression curves shown in Figs. 5(a,b). With these, we can analyze the effect of atomic motion on the H-bond state, its effect on the volume, and compare the latter with measurements [65, 22, 21]. As shown in Fig. 5(b), extrapolations of a Birch-Murnaghan EoS fit to results below and above $\sim$10 GPa produce divergent curves (dotted blue and red dotted lines), indicating a change in the OH-bond state at $\sim$10 GPa.

The 300 K isothermal bulk modulus, $K_{T}=V\,(\partial P/\partial V)_{T}$, shown in Fig. 5(c), amplifies the subtle change in compressibility across this critical point. Across this state change in the 0–20 GPa range, the OH bond-length distribution, $g_{\mathrm{OH}}(r)$, evolves from a double peak to a single modal distribution, though not symmetric (Fig. 5(d)). In this pressure range, the H-bond compresses quickly while the ionic OH bond stretches, a well-known behavior of these bonds (e.g., [71, 57]). $\delta$’s compressive behavior obtained from SCAN-DP calculation confirms previous PBE/GGA results [26] that this change in compressive behavior occurs in the 30–40 GPa range in static calculations (dashed line in Fig. 5(c)).

SCAN-DPMD results at 300 K account for the thermal expansion and reproduce well the 300 K experimental compression curve [65, 22, 21]: in the 0–15 GPa range, the difference between results and measurements is less than the differences between measurements; at higher pressures, the predicted SCAN-DPMD volume is slightly overestimated, with an error smaller than 0.25 Å³/f.u. Although predicted [24], it is striking to see that the inclusion of classical ionic motion lowers the transition pressure by $\sim$27 GPa, giving a transition pressure that agrees quite well with the experimentally measured one at $\sim$9 GPa (e.g., [22, 69]). This 300 K result is surprising and puzzling. At this low temperature, quantum proton motion, particularly tunneling [72], should impact the H-bond behavior (e.g., [73]).

Above 10 GPa, $V$ and $K_{T}$ display a smooth monotonic dependence on pressure. Table 1 summarizes the Birch-Murnaghan EoS parameters resulting from the fitting of measured and calculated 300 K $P,V$ data above 10 GPa. Fitting the SCAN-DPMD results within the 10–64 GPa pressure range with $K^{\prime}_{0}=4$ gives $V_{0}=55.5$ Å³, $K_{0}=225\pm 2$ GPa (red dotted line in Fig. 5(b)). The predicted $V_{0}$ is in excellent agreement with the measured $V_{0}$ [22] fit to the same Birch-Murnaghan EoS, 55.47 Å³. The predicted bulk modulus differs by $\sim$3% from the experimental one, 219 GPa [22], an impressive agreement for a result extrapolated to 0 GPa.

Our SCAN-DPMD results are significantly more accurate than previous PBE-BOMD results, which reported a $V_{0}=58.5$ Å³ and $K_{0}=183.4$ GPa using the Vinet EoS fit to results above 10 GPa [63]. This improvement can be attributed to (a) SCAN’s more accurate description of the $P$-$V$ relation, (b) the denser sampling of $(P,V)$ states over a broader pressure range, and (c) larger, longer, and better-converged simulation runs enabled by DPMD’s performance leap compared to BOMD. Nevertheless, this level of agreement between classical MD results and measurements at 300 K is unexpected.

This good agreement between SCAN-DPMD results and experimental data above 10 GPa suggests that our simulations describe well the H-bond state above this pressure. As indicated in Fig. 5(d), not even at 20 GPa, the first two peaks in the $g_{\mathrm{OH}}(r)$ pair distribution function have merged, indicating that H-bond symmetrization has not yet been achieved in the simulations. Therefore, the change in compressive behavior at $\sim$10 GPa cannot be attributed to H-bond symmetrization in this classical simulation. Fig. 6 shows in detail the evolution of $g_{\mathrm{OH}}(r)$ and $g_{\mathrm{HH}}(r)$ pair distribution functions from 0 to 10 GPa. Both show two clear peaks at 0 GPa and tend to merge with increasing pressure. However, only the $g_{\mathrm{HH}}(r)$ peaks are fully merged at 10 GPa. $g_{\mathrm{OH}}(r)$ still displays two superposing broad peaks. This indicates that H-bonds are not symmetric in these classical simulations at 10 GPa (or 20 GPa). The change in the H-bonds’ state at 10 GPa seems related to a change in H-ordering, likely into a disordered state. Our previous multi-configuration PBE-QHA study [27] suggested several short- to medium-range ordered H-bonded configurations remain dynamically stable beyond the pressure where the compressive behavior changes ($\sim$10 GPa in experiments, $\sim$12 GPa in PBE-QHA calculations). This result supported the notion of a more disordered H-bond state. The current results showing an asymmetric first peak in $g_{\mathrm{OH}}(r)$ and the fully merged broad peak in $g_{\mathrm{HH}}(r)$, and the previously identified multi-configuration equilibrium state [27] suggest that the change in compressive behavior at $\sim$10 GPa is associated with the emergence of a disordered state in the simulations. Static or dynamic, disorder involves proton hopping and is a precursor to proton diffusion at higher temperatures. This 300 K behavior might change significantly if the proton dynamics is addressed quantum mechanically using path integrals [73].

The DP potential allows us to systematically perform large-scale DPMD simulations to investigate proton diffusion in $\delta$ over a wide $P,T$ range. We perform 96 $NVT$ DPMD simulations covering the pressure range of 35–140 GPa and temperature range of 1,500–3,000 K (see Fig. S2). They contain 8,192 atoms and run for 2 ns. We do not observe diffusion, structural change, or melting of the Al and O sub-lattices in this $P,T$ range.

Proton diffusion is quantitatively characterized by protons’ mean-square displacement (MSD). Protons’ MSD over a simulation run time $t$, $L_{\mathrm{H}}^{2}(t)$, is defined by [74, 75]

$$L_{\mathrm{H}}^{2}(t)=\frac{1}{N_{\mathrm{H}}}\bigg{\langle}\sum_{i=1}^{N_{\mathrm{H}}}\big{\lvert}R_{\mathrm{H},i}(t+\tau)-R_{\mathrm{H},i}(\tau)\big{\rvert}^{2}\bigg{\rangle},$$

(1)

where $N_{\mathrm{H}}$ is the total number of protons, $R_{\mathrm{H},i}(\tau)$ denotes the position of $i$-th proton at moment $\tau$, and “$\langle\,\cdot\,\rangle$” denotes the ensemble average over the start time $\tau$. A linear (non-linear) dependence of the MSD on the simulation run time, $t$, reflects a continuous (irregular) diffusive behavior. The self-diffusion coefficient, $D_{\mathrm{H}}$, or the diffusivity, is related to the slope of the MSD ($L^{2}$) vs. $t$ line,

$$D_{\mathrm{H}}=\lim_{t\to\infty}\frac{L_{\mathrm{H}}^{2}}{6t}\,.$$

(2)

Figs. 7(a,b) show the protons’ MSD over a time span of 1 ns for $V=47$ Å³ and $T=$ 1,500 K (60.9 GPa) and 2,400 K (67.5 GPa), respectively. MSD at several other $V,T$ conditions are available in Fig. S3. Increasing the temperature significantly increases protons’ diffusivity (see Fig. S3). Diffusion is not observed at 1,500 K. It starts between 1,800 K and 2,100 K and becomes steady at 2,400 K. Figs. 7(c,d) show the corresponding proton trajectories. Fig. 7(b), shows that the protons’ MSD along the $c$ direction is approximately twice as large as those along the $a$ and $b$ axes, indicating they move more freely through the interstitial channels along the $c$ direction (see Figs. 3(a,b)). A similar phenomenon has also been reported in the NAM phase hydrous Al-bearing stishovite, which has a similar CaCl₂-type Si-O framework [76].

At 1,800–3,000 K, the relationship between $\log D_{\mathrm{H}}$ and $1/T$ is linear at all volumes and can be fitted using the Arrhenius equation (see Fig. S4),

$$D_{\mathrm{H}}(T)=D_{0}\,\exp\bigg{(}\frac{-E_{\mathrm{a}}}{k_{\mathrm{B}}T}\bigg{)},$$

(3)

where $k_{\mathrm{B}}$ is the Boltzmann constant and $E_{\mathrm{a}}$ denotes the activation energy. This behavior indicates that proton diffusion in $\delta$ is a common thermally activated process. The volume dependence of the activation energy, $E_{\mathrm{a}}$, is nearly linear (Fig. S5). Using these relationships, we interpolate $D_{\mathrm{H}}$ vs. $V$ and $T$ using the finite temperature EoS displayed in Fig. 4. The magnitude of protons’ diffusivity is shown as the background color in Fig. 8. For diffusivity corresponding to MSD $>1~{}\mbox{\AA}^{2}/\mbox{ns}$, (e.g., Fig. 7), diffusion is steady, stable, and characterized by a linear dependence of the MSD vs. time (reddish background area); for diffusivity of $<0.1~{}\mbox{\AA}^{2}/\mbox{ns}$ (bluish background area), diffusion can be intermittent with non-linear dependence of the MSD vs. time and the system is considered to be solid [78]. Depending on the pressure, for $\mbox{1,800}~{}\mathrm{K}<T<\mbox{2,100}~{}\mathrm{K}$, $D_{\mathrm{H}}$ approaches 1 Å²/ns and starts deviating from the high-$T$ Arrhenius behavior (see Fig. S4). This region of the diffusivity diagram is transiting from a blueish to a whitish background. At these $P,T$ conditions, one can recognize in Fig. S3 a change in the MSD’s behavior from non-linear to linear. We designate this regime as the onset of the fully superionic behavior in the $P,T$ phase space and represent this $D_{\mathrm{H}}=1$ Å²/ns diffusion boundary in Fig. 8 with a solid black curve.

The experimental dehydration boundary reported in Ref. [21] resembles our diffusion boundary. Similar measurements in Ref. [42] pinpointed two $\delta$ dehydration temperatures in two distinct pressure ranges: 1,850–1,900 K at 28–68 GPa and $\sim$1,959 K at 100–110 GPa. Two divergent boundaries were reported in the mid-lower mantle (68–100 GPa) due to uncertainties related to unidentified X-ray diffraction peaks and uncertain experimental conditions (Fig. 8). Their lower boundary directly connects the low-pressure and high-pressure boundaries and runs parallel to our boundary with a 200 K downward shift. It has been argued [42] that the disagreement between these studies could be due to uncertainties in granularity or the preferred orientation of samples in the high-pressure cell, which could affect the dehydration product detection sensitivity. The proximity of our diffusion boundary and the experimental dehydration boundary, particularly their similar pressure gradients, suggests that continuous proton diffusion controls the dehydration process as in antigorite [33] or diopside with H defects [34].

The diffusion boundary previously obtained with PBE-BOMD [29] largely overestimates the superionic transition temperature. It is known that supercell size could affect the superionic diffusivity calculations, often leading to overestimation of the transition temperature [79, 80]. The short simulation time and the PBE’s inaccurate description of the H-bond could also contribute to the overestimation of the diffusion boundary.

It is also possible to investigate the change in protons’ diffusivity by inspecting the constant-volume heat capacity, $C_{V}$, through the $P,T$ range shown in Fig. S6. The constant-volume heat capacity, $C_{V}$, can be calculated from MD ensemble averages [81] as

$$NC_{V}=\frac{\mathrm{var}(E)}{k_{\mathrm{B}}T^{2}}=\frac{1}{k_{\mathrm{B}}T^{2}}\!\left[\big{\langle}E^{2}\big{\rangle}-\big{\langle}E\big{\rangle}^{2}\right],$$

(4)

where $k_{\mathrm{B}}$ denotes the Boltzmann’s constant, $N$ denotes the number of atoms, $\langle E^{2}\rangle$ and $\langle E\rangle^{2}$ denotes the run average of energy $E$ and its square $E^{2}$. $C_{V}$ computed using this method does not exhibit a size effect for simulation cell sizes varying from 128–27,648 atoms in both the solid and diffusive states (Fig. S6). $C_{V}$ of an approximately harmonic solid is expected to plateau at the Dulong-Petit limit, $3NK_{\mathrm{B}}$. Deviation from this value suggests strongly anharmonic behavior, in the present case, entering the diffusive regime.

Fig. 9 shows $\delta$’s $C_{V}$ vs. $P,T$ as background color. For $T<\mbox{1,500}$ K (blueish background), $C_{V}$ generally follows the Dulong-Petit limit with at most 5% excess. For $T>$ 1,800 K, $C_{V}$ increases quickly. Depending on the pressure, the background color turns from blue to white then red between 1,800 K and 2,100 K, indicating $>10$% deviation from $3Nk_{B}$ within this color change range. The $D_{\mathrm{H}}=1$ Å²/ns diffusion boundary corresponds to $C_{V}$ $\sim$6% above the Dulong-Petit limit. Measurements of $C_{V}$ could be useful to determine the onset of $\delta$’s superionic state.

The NN potential for $\delta$ was developed based on the Deep Potential Smooth Edition (DeepPot-SE) model [48] implemented in DeePMD-kit v2.1 [88, 89]. Two-body embedding with coordinates of the neighboring atoms (se_e2_a) was used for the descriptor. The embedding network shape is (25, 50, 100). The fitting network shape is (128, 128, 128). The cut-off radius is 6 Å, and the smoothing parameter is 0.5 Å.

The model was trained using the Adam optimizer [90] for $1\times 10^{6}$ training steps, with the learning rate exponential decaying from $1\times{10}^{-3}$ to $3.51\times{10}^{-8}$ throughout the training process. The loss function $\mathcal{L}(p_{e},p_{f})$ is [88]

$$\mathcal{L}(p_{e},p_{f})=p_{e}\lvert\Delta e\rvert^{2}+\frac{p_{f}}{3N}\lvert\Delta f_{i}\rvert^{2},$$

(5)

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

The DP-GEN concurrent learning scheme [91] was employed to create the training data set and to generate the potential. We randomly extracted 118 labeled configurations from 59 BOMD runs at various $V,T$s to generate the initial potentials and to kickstart the DP-GEN training process. 6 DP-GEN iterations were performed to explore the configuration space and to eventually generate a potential reaching satisfactory accuracy requirement for DPMD for a temperature range of $300<T<3,000$ K. 4 candidate DP potentials initialized with different random seeds were trained in each iteration. They were used to perform $NVT$ DPMD simulations at various $V,T$s for a few thousand timesteps. After the simulations, the error estimator (model deviation), $\epsilon_{t}$, is calculated every 50 MD steps based on the force disagreement between the candidate DPs [51, 91]:

$$\epsilon_{t}=\max_{i}\sqrt{\left<\|F_{w,i}(\mathcal{R}_{t})-\left<F_{w,i}(\mathcal{R}_{t})\right>\|^{2}\right>}\,,$$

(6)

where $F_{w,i}(\mathcal{R}_{t})$ denotes the force on the $i$-th atom predicted by the $w$-th potential for the $\mathcal{R}_{t}$ configuration. For a particular configuration, if $\epsilon_{t}$ satisfies $\epsilon_{\mathrm{min}}\leq\epsilon_{t}\leq\epsilon_{\mathrm{max}}$, the corresponding configuration is collected, then labeled with DFT forces and total energy then added into the training dataset; if $\epsilon_{t}<\epsilon_{\mathrm{min}}$, these configurations are considered “covered” by the current training dataset; if $\epsilon_{t}>\epsilon_{\mathrm{max}}$ are considered failed and were discarded. After a few iterations, almost no new configurations are collected according to this standard ($>99$% are “accurate” for a few iterations), the DP-GEN process is then complete. The parameters for these DP-GEN iterations are listed in Table SII in detail.

After these DP-GEN iterations, our training dataset consists of 3,487 configurations labeled with ab initio force and energy. DPMD simulations were performed using the LAMMPS code [92] with 0.5 fs timesteps.

Training and testing data sets were based on ab initio calculations performed with the Vienna Ab initio Simulation Package (VASP) v6.3 [93]. The strongly constrained and appropriately normed (SCAN) [94] meta-GGA functional with PAW basis sets were adopted. The cutoff energy for the plane-wave-basis set was set to 520 eV. The Brillouin zone sampling for the $2\times 2\times 4$ supercells (with 16 f.u. of $\delta$, or 128 atoms) used a shifted $2\times 2\times 2$ Monkhorst-Pack $k$-point mesh. BOMD simulations were performed with a timestep of 0.5 fs.

A dataset for validating the DP potential is created by performing both BOMD and DPMD $NVT$ canonical ensemble simulations on $2\times 2\times 4$ supercells (128 atoms) with the Nóse-Hoover (NH) thermostat/barostat [95]. These simulations are performed at 6 $T$s ($T=$ 300, 600, 1,200, 1,800, 2,400, and 3,000 K) and 4 $V$s (the conventional cell volume $V=$ 40, 44, 48, and 54 Å³). This mesh covers the pressure range of $\sim$20–150 GPa. These simulations start from configurations produced after $10^{4}$ DPMD equilibration timesteps at each given $(V,T)$ and run for a minimum of 1 ns at each $(V,T)$. The DPMD simulation time scales linearly with the number of atoms (see Fig. S7), similarly with previous benchmarks [96].

The investigation of the 300 K compression curve involves DPMD simulations with 1,024-atom (or $4\times 4\times 8$) supercells. We perform $NPT$ simulations for 50 ps to equilibrate the structure. Then, perform $NVT$ simulations for another 50 ps to obtain the pressure at the given volume. The investigation of high $P,T$s, involved systematic DPMD simulations with 8,192-atom (or $8\times 8\times 16$) supercells. We perform $NVT$ simulations. 96 DPMD simulations at different $(P,T)$’s cover the $P,T$ range from 1,500 K to 3,000 K and 10 GPa to 180 GPa (see Fig. S2). Each MD simulation at equilibrated $V,T$ conditions runs for 2 ns, or $4\times 10^{6}$ timesteps. These simulations were performed concurrently on 96 GPUs.

Our workflows, including the DP-GEN active learning iterations and subsequent MD simulations, were implemented and managed using the Snakemake [97, 98] workflow management system.