⇦Chinese Physics Letters, 2020, Vol. 37, No. 4, Article code 043101Express Letter Discriminating High-Pressure Water Phases Using Rare-Event Determined Ionic Dynamical Properties * Lin Zhuang (庄琳)1,2, Qijun Ye (叶麒俊)1,2**, Ding Pan (潘鼎)3,4, Xin-Zheng Li (李新征)1,2,5** Affiliations 1Institute of Condensed Matter and Material Physics, School of Physics, Peking University, Beijing 100871 2State Key Laboratory for Artificial Microstructure and Mesoscopic Physics, Frontier Science Center for Nano-optoelectronics and School of Physics, Peking University, Beijing 100871 3Department of Physics and Department of Chemistry, Hong Kong University of Science and Technology, Hong Kong 4HKUST Fok Ying Tung Research Institute, Guangzhou 511458 5Collaborative Innovation Center of Quantum Matter, Peking University, Beijing 100871 Received 7 March 2020, online 25 March 2020 *Supported by the National Basic Research Program of China under Grant Nos. 2016YFA0300900 and 2017YFA0205003, the National Science Foundation of China under Grant Nos. 11774003, 11634001, 11934003, and 11774072. D. Pan also acknowledges the support from Hong Kong Research Grands Council (Nos. ECS-26305017 and GRF-16307618), the Alfred R. Sloan Foundation through the Deep Carbon Observatory (DCO), and the Croucher Foundation through the Croucher Innovation Award.
**Corresponding authors. Email: xzli@pku.edu.cn; qjye@pku.edu.cn Citation Text: Zhuang L, Ye Q J, Pan D and Li X Z 2020 Chin. Phys. Lett. 37 043101    Abstract Recent discoveries of dynamic ice VII and superionic ice highlight the importance of ionic diffusions in discriminating high-pressure ($P$) water phases. The rare event nature and the chemical bond breaking associated with these diffusions, however, make extensive simulations of these processes unpractical to ab initio and inappropriate for force field based methods. Using a first-principles neural network potential, we performed a theoretical study of water at 5–70 GPa and 300–3000 K. Long-time dynamics of protons and oxygens were found indispensable in discriminating several subtle states of water, characterized by proton's and oxygen ion's diffusion coefficients and the distribution of proton's displacements. Within dynamic ice VII, two types of proton transfer mechanisms, i.e., translational and rotational transfers, were identified to discriminate this region further into dynamic ice VII T and dynamic ice VII R. The triple point between ice VII, superionic ice (SI), and liquid exists because the loosening of the bcc oxygen skeleton is prevented by the decrease of interatomic distances at high $P$'s. The melting of ice VII above $\sim$40 GPa can be understood as a process of two individual steps: the melting of protons and the retarded melting of oxygens, responsible for the forming of SI. The boundary of the dynamic ice VII and SI lies on the continuation line ice VII's melting curve at low $P$'s. Based on these, a detailed phase diagram is given, which may shed light on studies of water under $P$'s in a wide range of interdisciplinary sciences. DOI:10.1088/0256-307X/37/4/043101 PACS:31.15.xv, 66.30.jp, 81.30.Dz, 07.05.Mh © 2020 Chinese Physics Society Article Text Water is vital for all known forms of life in the universe. In water systems, the interacting electrons behave as a quantum glue for nuclei which interplay with the high dimensionality of the hydrogen bond (HB) network, resulting in a rich phase diagram over a wide range of temperatures ($T$'s) and pressures ($P$'s).[1–11] At different $T$'s, a variety of phases exist from the ambient pressures to $\sim $2 GPa.[1,4,12–15] Above this, body-centered cubic (bcc) or bcc-like structures of oxygens (O) dominate the phase diagram below the liquid phases till $\sim $400 GPa.[2,3,5,8,16] Among the different relatively low-$T$ phases, ice VII has a melting line intersecting a significant number of thermal profiles of the cold subducting slabs, meaning that it might exist in these cold subduction zones of the Earth.[17,18] This is critical to water transport in the Earth.[17,18] At higher $T$'s and $P$'s, a novel superionic phase exists, with its $P$–$T$ stability region penetrated by the isentropes of Neptune and Uranus.[19,20] The fast diffusing protons result in strong magnetic field, inducing a crucial role played by this superionic ice in depicting the features of these planets.[19,20] Therefore, understanding the properties of the bcc-based water phases and the liquid phases above them is a task of central importance, not only in physics and chemistry, but also in the Earth and planetary sciences. For all these bcc-based solids, ice VIII is the low-$T$ and relatively low-$P$ (i.e., $\sim $2 to $\sim $62 GPa) phase, upon which most of them are built.[1,21] In this structure, the O atoms site on a bcc-like lattice, with the H atoms arranged in order along two interpenetrating but not connected HB networks (inset in Fig. 1),[1,21,22] obeying the Bernal–Fowler–Pauling ice rules.[23] Upon increasing $T$'s, ice VIII transits into ice VII at $\sim $260 K from $\sim$2 to $\sim$10 GPa, with the O-atom bcc skeleton remaining but the H atoms getting disordered along the H-bonds.[1,22,24–26] Above $\sim $10 GPa, the transition $T$ deceases gradually to $\sim $100 K at $\sim $62 GPa, which is a triple point between these two phases and a new H-bond symmetrized phase named ice X.[16,22,26] Below (above) $\sim $100 K, water ice transits into phase X from phase VIII (VII) upon increasing $P$, as shown in Fig. 1. At $T$ and $P$ higher than $\sim $1000 K and $\sim $40 GPa, a superionic ice (SI) phase exists,[7,25,27–29] which has liquid-like protons and solid-like O ions staying on the bcc lattice. This SI was predicted by first-principles simulations in Ref. [30], characterized by an abrupt increase of diffusion coefficient of protons.[27,31–34] Using an improved diamond anvil cell technique in 2005, Goncharov et al. confirmed a triple point between ice VII, SI, and the liquid phase at $\sim $47 GPa and $\sim $1000 K,[27] which represents a seminal advance of our understanding of water under high $P$'s. Besides this, a dynamic ice VII was also suggested in Ref. [27], with more localized protons than SI and more disordered protons than the conventional ice VII. The rare event nature and the chemical bond breaking associated with the proton diffusions in this dynamic ice VII, however, made simulations of this part of the phase diagram unpractical to ab initio and inappropriate for force field based methods.[25,29,31,34–36] A systematic theoretical study of this part of the phase diagram is highly desired.
cpl-37-4-043101-fig1.png
Fig. 1. The melting curve of ice VII as a function of pressure given by previous experimental and theoretical studies: Goncharov et al.,[27] blue curve (Exp); Schwager et al.,[33] green curve (Exp); Redmer et al.,[37] faint yellow curve (Theory); Schwegler et al.,[25] red square (Theory). The black solid lines are the boundaries between ice X, ice VIII and ice VII. The insets are the structures of ice VIII, ice VII, and ice X. In ice VIII, the protons are ordered in the HB network. In ice VII, the protons become disordered in one of the two sites between oxygen pairs. In ice X, the protons are symmetrized between oxygen pairs.
In this Letter, we investigate the phase diagram of these bcc based solids and the liquid phases above them using a neural network (NN) potential,[38,39] which can enable simulations with accuracy comparable to density-functional theory (DFT) and computational cost similar to force fields. Beyond the reach of previous studies, extensive long-time & large scale molecular dynamic (MD) simulations at $T$ ranging from 300 to 3000 K and $P$ ranging from 5 to 70 GPa were carried out. We found that the long-time dynamics of these ions are indispensable in discriminating several subtle and non-trivial states of water, characterized by ionic diffusion coefficients and distributions of proton's displacements. This is beyond the reach of conventional structural properties. Two different mechanisms of protons transfer reflected by the distribution of proton's displacements, i.e., translational and rotational transfer, were found to be very important in the MD simulations, which can further discriminate dynamic ice VII into dynamic ice VII T (with T meaning that only translational proton transfer exists) and dynamic ice VII R (with R meaning that rotational proton transfer also exists). The triple point between ice VII, superionic ice (SI), and liquid exists because the decrease of interatomic distances can prevent the loosening of the bcc oxygen skeleton at high $P$'s. Therefore, the melting of ice VII above $\sim $40 GPa can be understood as a process of two individual steps, i.e., the melting of protons and the retarded melting of oxygens, responsible for the formation of SI. The boundary between the dynamic ice VII and SI lies on the continuation line of ice VII's melting curve at low $P$'s. Based on these, a detailed phase diagram is given, which will be useful to understand water under high $P$'s in a wide range of interdisciplinary sciences in future studies. A promising choice that can achieve both the accuracy of DFT and the efficiency of empirical potentials is offered by the machine learning techniques[40]. The ability of the machine learning methods is that complex high dimensional functions can be fitted to describe the potential energy surface using atomic coordinates as the basic variables $E=E(R_1,R_2,\dots,R_N)$. A large set of configurations generated via the MD or Monte-Carlo (MC) simulations are needed for the system of interest. The energies and forces of these configurations, used later as training and testing sets, are computed with the predictive first-principles calculations. The potential is then optimized using the machine learning algorithm. Among different options, NNs and Gaussian Process Regression have been extensively applied to simulations of a series of realistic system, in generating force fields without any empirical presets.[41,42] In this study, we have chosen the NN method.[38] The training of the NN potential is implemented via the DeePMD-kit package,[38] where the training set is generated from the DFT results implemented via VASP.[43,44] The DeePMD-kit package features to preserve all the natural symmetries by building for every atom a local coordinate frame.[38] In DeePMD, the Adam method is used to optimize the parameters of each layer of NN, as $$\begin{alignat}{1} \!\!\!\!\!\!L(p_\epsilon, p_f, p_\xi) = p_\epsilon \Delta \epsilon^2 \!+\! \frac{p_f}{3N} \sum_i|\Delta F_i|^2 \!+\! \frac{p_\xi}{9} \left\|\Delta\xi\right\|^2.~~ \tag {1} \end{alignat} $$ $N$ is the number of atoms. $\Delta \epsilon$, $\Delta F_i$, $\Delta \xi$ denote the differences between the DeePMD prediction and the training data of the energy per atom, the force on atom $i$, and the virial tensor divided by $N$, respectively. Iterations were performed to reduce the loss function until the aimed criterion was met. The MD simulations with NN potential are then preformed by the Large-scale Atomic/Molecular Massively Parallel Simulator (Lammps) with DeePMD module.[38,45] The input configurations for training this potential are summarized in Fig. 2(a). Considering that the atomic interacting environments, like O–O distances, are significantly affected by $P$, similar amount of training sets were used for all $P$'s studied. This can eliminate the potential preconception by inputs and convince the performance in the whole region of interest. In Fig. 2(b), we show the convergence of NN potential towards the DFT one by the reduced root-mean-square-error (RMSE) of energies and averaged forces. Above 200000 configurations, the RMSEs of trained energies (average forces) have reached the level of 2.4 meV/H$_2$O (120.7 meV/Å). For one specific $P$, the RMSEs of energy and averaged forces are even smaller, as shown in Fig. 2(c). The accuracy of trained potential can be perceived by the almost merged radial distribution functions (RDFs). For example, the RDFs of NN potential are in high coincidence with the ones of ab initio MD at 60 GPa and 2800 K shown in Fig. 2(d). More details are also given in our the Supplemental Materials (SM).
cpl-37-4-043101-fig2.png
Fig. 2. (a) The sampled frames of the training set at different pressures. (b) The RMSEs of energies (black square line) and forces (red circle line) with respect to trained frames. (c) The converged RMSEs of energies (green bar) and forces (blue bar) at different pressures. (d) The RDFs at 60 GPa and 2000 K. The solid (dashed) lines are the results from DFT (DeePMD).
In statistical physics, the phases are conventionally distinguished by the abrupt changes of the thermodynamic quantities. Accordingly, we start our discussions by looking at the densities (black up triangle line) and the potential energies (red down triangle line) at varying $T$'s, as shown in Fig. 3(a) for 60 GPa. With $T$ increased from 400 to 3000 K, the density (potential energy) decreases (increases) monotonously. There are two discontinuities at $\sim $1700 and $\sim $2600 K (blue shaded region in Fig. 3(a)), which split the water system into three known phases, i.e., conventional ice VII, superionic ice, and liquid. The transitions between them are of first order, as evidenced by the abrupt changes. These results are in alignment with the previous studies.[25,27,33] Then, motivated by the contemporary progress in understanding the concept of phases in spacetime,[46–48] we pose our system in the temporal degrees of freedom. Dynamic quantities are used to characterize the changes of system at different $T$'s, starting with the diffusion coefficients of protons ($D_{\rm H}$) and oxygens ($D_{\rm O}$). These diffusion coefficients are derived from the long-time behavior of mean square displacement (MSD), using $$ D_i = \lim_{\Delta t\to\infty}\frac{1}{6N \Delta t}\sum_{n=1}^{N}\left({\boldsymbol r}_{i}(t) - {\boldsymbol r}_{i}(t + \Delta t)\right)^2.~~ \tag {2} $$ Both $D_{\rm H}$ and $D_{\rm O}$ increase monotonously with $T$, and five zones with very different $D_{\rm H}$ and $D_{\rm O}$ can be seen in Fig. 3(b). Below 900 K, both values are negligibly small (less than $10^{-5}$ Å$^2$/ps, region A). Between 900 K and 1800 K, $D_{\rm O}$ remains negligibly small but $D_{\rm H}$ becomes finite and increases gradually. This qualitative change in $D_{\rm H}$ results in a division of region B, composed of B$^{\prime}$ and B$^{\prime\prime}$ whose details will be discussed later, out of the others. Above 1800 K, $D_{\rm H}$ is saturated and the regions are discriminated by $D_{\rm O}$ into C, D, and E. In region C ($T$ ranging from 1800 to 2100 K), the oxygens remain vibrating around the bcc sites while the protons become delocalized as in liquid, corresponding to the SI phase. When $T$ is higher than 2100 K, the oxygen atoms begin to diffuse and struggle to break the confinement of crystal structure, forming region D. At $\sim $2800 K, $D_{\rm H}$ and $D_{\rm O}$ saturate to different values (region E), meaning an ionic liquid (IL). The oxygens in region D are in the intermediate state between the SI and the diffusive IL phases, labeled as SI$^\prime$. It is clear that compared with Fig. 3(a), Fig. 3(b) presents much richer details of water under 60 GPa at different $T$'s. One deficiency of the approach of discriminating phases using diffusion coefficients is that these quantities might not give enough information of the dynamics at low $T$'s, when diffusions are rare events. Instead of the regular linear lineshape of MSD in the long time limit for liquid (with finite slope) and solid (horizontal), there are several plateaus arising randomly (see Fig. S4 in the SM). Taking the plateaus of the protons as an example, they mean that the protons are bounded in the local minima of the crystal sites for a long time before they escape, and it is only till such events happen very frequently when the plateaus merge into linear lineshape with finite slope. Otherwise, $D_{\rm H}$ is too small to be derived with reliable numerical accuracy (the gray shadow region in Fig. 3(b)). In these cases, we design another dynamic quantity, i.e., the distribution of proton's displacement, to analyze the low-$T$ phases. This distribution is a function of the proton's displacement $\Delta r$ in a certain time interval $\Delta t$, defined as $$ P_{\Delta t}(\Delta r) =\sum_i P\left(\Delta r_i\right),~~ \tag {3} $$ where $\Delta r_i(t) = |{\boldsymbol r}_i(t) - {\boldsymbol r}_i(t + \Delta t)|$, for $\forall t, \forall i=1,\ldots,N_{\rm H}$. Displacements $\Delta r_i(t)$'s for all protons in their trajectories are accumulated to derive the distribution. When $\Delta t$ is vanishingly short, $\Delta r$ is trivial and linear in velocity, and $P_{\Delta t}(\Delta r)$ follows the Maxwell velocity distribution (see Fig. S5 of the SM). With increased $\Delta t$, $\Delta r$ and $P_{\Delta t}(\Delta r)$ vary depending on the motion of the protons. This motion can be attributed to three types: the local vibration around one local minimum, translational proton transfer, and rotational proton transfer. In high-$P$ ice VII, the first two and their combinations happen much more frequently than the third one. At low $T$'s, absence of type-3 motion results in localized protons, either in one site or between the neighboring oxygen pairs. $P_{\Delta t}(\Delta r)$ converges with increasing $\Delta t$, during which the protons can sufficiently explore the local surroundings. Introducing the rotational transfer, however, protons can travel along the HB network. $\Delta r$ can exceed the O–O distance when large $\Delta t$ is accumulated, resulting in broadened $P_{\Delta t}(\Delta r)$. Therefore, one can rely on the convergence behaviors of $P_{\Delta t}(\Delta r)$ on $\Delta t$ and the distribution of $P_{\Delta t}(\Delta r)$ as a function of $\Delta r$ to quantitatively describe the motion of protons.
cpl-37-4-043101-fig3.png
Fig. 3. (a) The densities (black up triangle) and potential energies (red down triangle) as a function of temperature at 60 GPa. The blue shaded regions show the two sudden jump in the density and potential energy lines and the region is divided into three regions: conventional ice VII, superionic ice, and liquid, sequentially. (b) The diffusion coefficient of oxygen ($D_{\rm O}$, red oblique cross line) and protons ($D_{\rm H}$, black vertical cross line). The points in gray shaded region are shown for schematic reason to denote the situation when $D_{\rm O}$ and $D_{\rm H}$ are negligibly small to give numerical accurate values. (c) The distribution of proton's displacements $P_{\Delta t}(\Delta r)$ as a function of fractional coordinate at 500 K with different $P$'s. For convenience of comparison, the $\Delta r$ is in units of $l_{\rm a}$, the simulated unit cell lattice length at each ($P$, $T$), which for instance is 2.93 Å at 30 GPa (more complete information of $l_{\rm a}$ see Table S1 in the SM). (d) The same $P_{\Delta t}(\Delta r)$ at 60 GPa with different $T$'s.
An example of the 500 K simulations is shown in Fig. 3(c). $P_{\Delta t}(\Delta r)$ for all pressures simulated is converged using $\Delta t = 0.4$ ps (see Fig. S5 of the SM). When $P$ is lower than 30 GPa, there is a single peak in $P_{\Delta t}(\Delta r)$, which can be fitted to the Burr distribution. It corresponds to the case when neither translational nor rotational proton transfer happens, i.e., static ice VII. At higher $P$'s, a shoulder appears at a larger $\Delta r$ ($\sim$$0.20 l_{\rm a}$). This $\Delta r$, which is $\sim $0.59 Å at 30 GPa, equals the distance of translational proton transfer between the neighboring oxygens. It contributes to non-negligible $D_{\rm H}$, although the proton is confined within each HB. Therefore, this state of matter with a shoulder in $P_{\Delta t}(\Delta r)$ corresponds to a kind of dynamic ice VII, when translational proton transfer allows proton to jump within each HB, but the absence of rotational proton transfer prohibits it from traveling to further distances. We label it as dynamic ice VII T. With increasing $T$'s and $P$'s, the rotational proton transfer becomes important. To demonstrate how this happens, one example is given for the 60 GPa simulations in Fig. 3(d). A large $\Delta t = 0.4$ ps was taken for all $T$'s to ensure non-negligible $P_{\Delta t}(\Delta r)$ at large $\Delta r$ to be sampled, induced by the rotational plus translational proton transfer at higher $T$'s, as well as the convergence of $P_{\Delta t}(\Delta r)$ for the localized protons at lower $T$'s. For $T$'s ranging from 400 to 1200 K, the distribution is concentrated but not single-peaked, corresponding to the dynamic ice VII-T introduced in the former paragraph. This $T$ range matches perfectly with regions A plus B$^{\prime}$ of $D_{\rm H}$ in Fig. 3(b). Above 1200 K, part of the displacements falls into the region with $\Delta r$ larger than $\sqrt3/2 l_{\rm a}$, the length of the nearest O–O distance. Therefore, the additional peak is a clear evidence that protons have escaped from the restraint of the oxygen pair and begin to diffuse along the HB network. Based on these, we discriminate dynamic ice VII into the dynamic ice VII T (region A and B$^{\prime}$, without rotational proton transfer) and dynamic ice VII R (region B$^{\prime\prime}$, with rotational proton transfer). The inherent proton transfer mechanism could also explain the slope change of $D_{\rm H}$ at 1200 K in Fig. 3(b). Above 1800 K, the system becomes SI, with peaks being inconspicuous and protons diffusing freely. Using the dynamical quantities as discussed above, and exploiting the feature of the NN potential that it is much cheaper than the ab initio calculations, we now come to a thorough determination of the water phase diagram at 5–70 GPa and 300–3000 K, as shown in Fig. 4. Each symbol denotes two simulations composed by an $NPT$ one to average the volume and a $NVT$ one afterwards to derive the dynamical quantities. The color of the symbol shows the ratio between $D_{\rm H}$ and $D_{\rm O}$. In the liquid region, values close to 1 (red symbols) mean that they are molecular liquid (ML). In the solid region, $D_{\rm O}$ is negligible. Therefore, the symbols are empty and we mainly resort to the distribution of proton's displacement to discriminate different phases. The dashed lines are the boundary between phases determined by these criteria. To estimate the impact of hysteresis effects on the melting $T$'s,[25] we performed two-phase simulations as in Ref. [49], and found an error of $\sim $50 K (up and down triangles in Fig. 4, see Figs. S10 and S11 of the SM for details). This value is small compared to the distance of the phase boundary curves as shown in Fig. 4. It is clear that this phase diagram is much richer than that given by the conventional thermodynamically averaged structural properties, which by definition can tell only four regions labeled by (a)–(d) in Fig. 4. The proton's displacement distribution as introduced above can discriminate two subtle phases, i.e., dynamic ice VII T and dynamic ice VII R, within dynamical ice VII. Besides this, the intermediate state between the SI and the IL above it, i.e., SI$^\prime$, can also be characterized by the diffusion coefficients based on long time MD simulations. We also notice that the boundary between dynamic ice VII R and superionic ice above $\sim $40 GPa lies on the continuation line of the boundary between dynamic ice VII R and the liquid phase at lower $P$'s, indicating that they are of similar nature. In the liquid region, different from ML, ionic liquid (IL) exists at higher $P$'s. More details concerning the characteristics of these phases please see Table S2.
cpl-37-4-043101-fig4.png
Fig. 4. The phase diagram of H$_2$O. The symbols indicate the simulated $(T, P)$ conditions. The grey circles belong to solid phases and colored squares belong to liquid phases. Combining the aforementioned dynamic quantities, we can distinguish seven phases: static ice VII, dynamic ice VII T, dynamic ice VII R, SI, SI$^\prime$, molecular liquid (ML) and ionic liquid (IL). The blue and black dashed lines are determined by $D_{\rm O}$, $D_{\rm H}$ and the red dashed lines are determined by the $P_{\Delta t}(\Delta r)$. The upper panels show the accumulated atomic trajectories of 5 ps for four atomic patterns (from left to right): static ice VII (10 GPa and 400 K), dynamic ice VII (20 GPa and 700 K), SI (60 GPa and 2000 K), and ML (10 GPa and 1200 K). For liquids, we use the color map of the ratio ($D_{\rm H}$/$D_{\rm O}$) to indicate its tendency to ML (ratio close to 1, red colored) or IL (ratio significantly larger than 1, the other colored). The melting data determined using the two-phase simulation method are shown by up and down purple triangles, with the former meaning the highest crystalizing and the latter meaning the lowest melting $T$'s. The hysteresis effects change the melting line by $\sim $50 K.
One key difference between the melting of ice here and that in other atomic or molecular systems is that in the latter cases it is induced by the instability of acoustic modes.[50] Above melting point, all the identical least units (either the atoms or the molecules) equally escape from the crystal sites and the crystal structure breaks down. The melting at high-$P$ water, however, poses a clear two-step nature: the melting of protons, and a retarded melting of oxygens. The retarded melting suggests that the primary instability roots in a mode related to protons and a "sublattice melting" exists.[50] It is the same picture that previous studies interpreted the melting of alkali halide and other superionic conductors:[51–53] the heavier ions experience rattling motions in the cage formed by its neighboring heavier atoms, while the lighter ions fast diffuse and fill as lube to hold the cages. Mass effects were regarded prevail,[52] however the interatomic distances are also important. The retarded melting of oxygens appears only at high $P$'s, where interatomic distances are shortened and the Coulomb repulsion is strengthened. Oxygens with relative large electron shell are sensitive to this change. After the triple point, they are stuck and can only escape at even higher $T$'s. As exposed protons, the hydrogen ions are little affected and melt as usual. Accordingly, the melting line of protons (the black line in Fig. 4) is continued throughout the $P$ range. The different melting process could also be viewed from the binding tendency of oxygens and protons in low $T$ solid correspondence. At low $P$'s, the interatomic potential poses a robust molecular form, implied by the only-allowed local motion and rotational transfer. While at high $P$'s, the reduced interatomic distances drive the HBs to be symmetric and the motions of protons and oxygens tend to be departed. When $T$ is high, the former one leads to a simultaneous melting and the latter to a departed one. Another question of general interest concerns the dissociation of water at high $P$'s. Two scenarios were proposed, with the first one like that at ambient $P$'s, i.e., a bimolecular reaction as 2H$_2$O $\rightarrow$ OH$^-$ + H$_3$O$^+$ suggested by the measured lower enthalpy for the hydration than the ionization,[54] and the second one very different from that, i.e., unimolecular reaction as H$_2$O $\rightarrow$ OH$^-$ + H$^+$, as suggested by shock compression experiments.[55,56] The behavior of the protons in our simulations is in alignment with the second scenario, suggested by the majority of OH$^-$ than H$_3$O$^+$ species and the much shorter lifetime of the latter (see Fig. S9 of the SM). This is also consistent with the earlier theoretical study by Goldman et al.[57] Making use of the long-time and supercell nature of the MD simulations, we can further verify this scenario by comparing the computed conductivity with available experimental and theoretical results.[57,58] The Nernst–Einstein relation is used and the conductivity $\sigma (T)$ is expressed as $$\begin{alignat}{1} \!\!\!\!\!\!\sigma(T) = \left[ n_{-}(T)Z^{2}_{-}e^{2}D_{-} + n_{+}(T)Z^{2}_{+}e^{2}D_{+} \right] / k_{_{\rm B}}T,~~ \tag {4} \end{alignat} $$ $n(T)$ is the density of ionic carriers at a given $T$, $Ze$ is their charge, $D$ is their diffusion coefficient, and the $k_{_{\rm B}}$ is Boltzmann's constant. The charge $Ze$ is approximated by the average Bader charges over 1000 equilibrium configurations for each $(T, P)$ point. Our results agree on the order of magnitude with earlier theory and experiments (Fig. 5). Besides accuracy of the MD samplings, the remaining differences may also originate from three other sources. The first one is the fact that the NN potential does not allow the charges of ions to be calculated on the fly. The averaged Bader charges over sampled equilibrium configurations for each $(T, P)$ point may induce additional error for the estimation of ionic charges. The second one is the fact that the experiments were carried out by shock waves, whose thermodynamic condition is different from the pure equilibrium state simulated here. The third one is related to the fact that it is the total conductivity which was measured in experiments, whilst our simulations address ionic conductivity. Concerning the trend, $\sigma (T)$ increases rapidly with $T$, followed by a slower increase as the system enteres the SI region, which allows us to assign its superionic nature by making analogy to other known superionic conductors. Comparison with the experimental data of different types of superionic conductors in Ref. [59] shows that it is a type-II superionic conductor. However, different from PbF$_2$, a prototype of this class where the F$^-$ anions diffuse in the Pb$^{2+}$ sublattice, cations like protons diffuse in the O$^{2-}$ sublattice. The large conductivity also supports the explanation that Neptune's wonky magnetic field originates from the interfacial SI layer.[60]
cpl-37-4-043101-fig5.png
Fig. 5. The conductivity as a function of temperature at 40 GPa (black circles) and 60 GPa (red squares). The experimental data (empty square) were acquired from Yakushev et al.,[58] whose estimated $P$ is $\sim $62.6 GPa for $T=1600$ K and $\sim $94.2 GPa for $T=2300$ K. One earlier theoretical result (empty circle) from Goldman et al.[57] was also shown, whose estimated $P$ is $\sim $42.3 GPa. Inset: the conductivity of AgI (dashed curve), PbF$_2$ (solid curve) and Na-$\beta$-Al$_2$O$_3$ (dotted curve), which show type-I, type-II and type-III superionic solids, respectively, from Ref. [61].
Overall we present in this manuscript an extensive theoretical study on high-$P$ water, allowed by using an NN potential. The importance of ionic long-time dynamics in discriminating its phase diagram is highlighted. Given the crucial role played by this matter in the universe and the potential generalization of the idea, i.e., dynamic feature determined phases, to broad systems even those regarded well-understood, we believe that this work may stimulate many forthcoming studies in a wide range of physical chemistry systems and interdisciplinary sciences. The authors sincerely acknowledge the help from L. F. Zhang and H. Wang, the authors of the DeePMD code. The computational resources were provided by the supercomputer center in Peking University, China. References New High-Pressure Phase of IceInfrared absorption study of the hydrogen-bond symmetrization in ice to 110 GPaThe infrared spectra of amorphous solid water and ice Ic between 10 and 140 KTunnelling and zero-point motion in high-pressure iceIce XV: A New Thermodynamically Stable Phase of IceNew Phases of Water Ice Predicted at Megabar PressuresHigh pressure partially ionic phase of water iceThe phase diagram of high-pressure superionic iceStability of hydrous silicate at high pressures and water transport to the deep lower mantleSuperionic-Superionic Phase Transitions in Body-Centered Cubic H 2 O IceOptical Spectra of Orientationally Disordered Crystals. II. Infrared Spectrum of Ice Ih and Ice Ic from 360 to 50 cm −1Absorptivity of Ice I in the Range 4000–30 cm −1Covalency of the Hydrogen Bond in Ice: A Direct X-Ray MeasurementA determination of the crystal structure of ice XIDynamical Instabilities of Ice XPossible presence of high-pressure ice in cold subducting slabsMagnetism in cold subducting slabs at mantle transition zone depthsSuperionic to Superionic Phase Change in Water: Consequences for the Interiors of Uranus and NeptuneThe unusual magnetic fields of Uranus and NeptuneHigh pressure phase transitions and hydrogen‐bond symmetry in ice polymorphsStructure of Ice-VII and Ice-VIII: A Quantum Mechanical Study †The Structure and Entropy of Ice and of Other Crystals with Some Randomness of Atomic ArrangementReassigning Hydrogen-Bond Centering in Dense IceMelting of ice under pressureDynamic Ionization of Water under Extreme ConditionsDiffusion and electrical conductivity in water at ultrahigh pressuresDynamical Screening and Ionic Conductivity in Water from Ab Initio SimulationsSuperionic and Metallic States of Water and Ammonia at Giant Planet ConditionsHigh pressure-temperature Raman measurements of H[sub 2]O melting to 22 GPa and 900 KInfrared absorption study of the hydrogen-bond symmetrization in ice to 110 GPaMelting curve of H 2 O to 90 GPa measured in a laser-heated diamond cellExtended and accurate determination of the melting curves of argon, helium, ice ( H 2 O ) , and hydrogen ( H 2 ) Anharmonic Raman Spectra in High-Pressure Ice from Ab Initio SimulationsThe phase diagram of water and the magnetic fields of Uranus and NeptuneDeep Potential Molecular Dynamics: A Scalable Model with the Accuracy of Quantum MechanicsReinforced dynamics for enhanced sampling in large atomic and molecular systemsPerspective: Machine learning potentials for atomistic simulationsGeneralized Neural-Network Representation of High-Dimensional Potential-Energy SurfacesGaussian Approximation Potentials: The Accuracy of Quantum Mechanics, without the ElectronsEfficient iterative schemes for ab initio total-energy calculations using a plane-wave basis setFrom ultrasoft pseudopotentials to the projector augmented-wave methodFast Parallel Algorithms for Short-Range Molecular DynamicsSpace-time thermodynamics of the glass transitionSpectral Statistics and Many-Body Quantum Chaos with Conserved ChargeQuantum Time Crystals from Hamiltonians with Long-Range InteractionsQuantum simulation of low-temperature metallic liquid hydrogenTheory of melting based on lattice instabilityTransport properties of molten alkali halidesPotentials and correlation functions for the copper halide and silver iodide melts. II. Time correlation functions and ionic transport propertiesStatic structure and ionic transport in molten AgBr and AgClThe great fallacy of the H+ ion: And the true nature of H3O+Electrical conductivity of water compressed dynamically to pressures of 70–180 GPa (0.7–1.8 Mbar)Spontaneous Raman Scattering from Shocked WaterAb initio simulation of the equation of state and kinetics of shocked waterHigh-spin-low-spin transition in magnesiowüstite (Mg0.75,Fe0.25)O at high pressures under hydrostatic conditionsSuperionics: crystal structures and conduction processesExperimental evidence for superionic water ice using shock compressionSuperionic conductors: Transitions, structures, dynamics
[1]Whitworth R and Petrenko V 1999 Physics of Ice (Oxford: Oxford University Press)
[2] Benoit M, Bernasconi M, Focher P and Parrinello M 1996 Phys. Rev. Lett. 76 2934
[3] Aoki K, Yamawaki H, Sakashita M and Fujihisa H 1996 Phys. Rev. B 54 15673
[4] Hagen W, Tielens A and Greenberg J M 1981 Chem. Phys. 56 367
[5] Benoit M, Marx D and Parrinello M 1998 Nature 392 258
[6] Salzmann C G, Radaelli P G, Mayer E and Finney J L 2009 Phys. Rev. Lett. 103 105701
[7] Militzer B and Wilson H F 2010 Phys. Rev. Lett. 105 195701
[8] Wang Y, Liu H, Lv J, Zhu L, Wang H and Ma Y 2011 Nat. Commun. 2 563
[9] Sun J, Clark B K, Torquato S and Car R 2015 Nat. Commun. 6 8156
[10] Rozsa V, Pan D, Giberti F and Galli G 2018 Proc. Natl. Acad. Sci. USA 115 6952
[11] Hernandez J A and Caracas R 2016 Phys. Rev. Lett. 117 135503
[12] Bertie J E and Whalley E 1967 J. Chem. Phys. 46 1271
[13] Bertie J E, Labbé H J and Whalley E 1969 J. Chem. Phys. 50 4501
[14] Isaacs E D, Shukla A, Platzman P M, Hamann D R, Barbiellini B and Tulk C A 1999 Phys. Rev. Lett. 82 600
[15] Howe R and Whitworth R W 1989 J. Chem. Phys. 90 4450
[16] Caracas R 2008 Phys. Rev. Lett. 101 085502
[17] Bina C R and Navrotsky A 2000 Nature 408 844
[18] Kupenko I, Aprilis G, Vasiukov D M, McCammon C, Chariton S, Cerantola V, Kantor I, Chumakov A I, Rüffer R, Dubrovinsky L and Sanchez-Valle C 2019 Nature 570 102
[19] Wilson H F, Wong M L and Militzer B 2013 Phys. Rev. Lett. 110 151102
[20] Nellis W J 2015 Mod. Phys. Lett. B 29 1430018
[21] Schweizer K S and Stillinger F H 1984 J. Chem. Phys. 80 1230
[22] Kuo J L and Klein M L 2004 J. Phys. Chem. B 108 19634
[23] Pauling L 1935 J. Am. Chem. Soc. 57 2680
[24] Benoit M, Romero A H and Marx D 2002 Phys. Rev. Lett. 89 145501
[25] Schwegler E, Sharma M, Gygi F and Galli G 2008 Proc. Natl. Acad. Sci. USA 105 14779
[26]Water phase diagram http://www1.lsbu.ac.uk/ (2019) accessed 14 November 2019
[27] Goncharov A F, Goldman N, Fried L E, Crowhurst J C, Kuo I F W, Mundy C J and Zaug J M 2005 Phys. Rev. Lett. 94 125508
[28] French M, Mattsson T R and Redmer R 2010 Phys. Rev. B 82 174108
[29] French M, Hamel S and Redmer R 2011 Phys. Rev. Lett. 107 185901
[30] Cavazzoni C, Chiarotti G L, Scandolo S, Tosatti E, Bernasconi M and Parrinello M 1999 Science 283 44
[31] Lin J F, Militzer B, Struzhkin V V, Gregoryanz E, Hemley R J and Mao H K 2004 J. Chem. Phys. 121 8423
[32] Lin J F 2005 Geophys. Res. Lett. 32 L11306
[33] Schwager B, Chudinovskikh L, Gavriliuk A and Boehler R 2004 J. Phys.: Condens. Matter 16 S1177
[34] Datchi F, Loubeyre P and LeToullec R 2000 Phys. Rev. B 61 6535
[35]Ohtani E 2007 Advances in High-Pressure Mineralogy (Boulder: Geological Society of America)
[36] Putrino A and Parrinello M 2002 Phys. Rev. Lett. 88 176401
[37] Redmer R, Mattsson T R, Nettelmann N and French M 2011 Icarus 211 798
[38] Zhang L, Han J, Wang H, Car R and W E 2018 Phys. Rev. Lett. 120 143001
[39] Zhang L, Wang H and W E 2018 J. Chem. Phys. 148 124113
[40] Behler J 2016 J. Chem. Phys. 145 170901
[41] Behler J and Parrinello M 2007 Phys. Rev. Lett. 98 146401
[42] Bartók A P, Payne M C, Kondor R and Csányi G 2010 Phys. Rev. Lett. 104 136403
[43] Kresse G and Furthmüller J 1996 Phys. Rev. B 54 11169
[44] Kresse G and Joubert D 1999 Phys. Rev. B 59 1758
[45] Plimpton S 1995 J. Comput. Phys. 117 1
[46] Merolle M, Garrahan J P and Chandler D 2005 Proc. Natl. Acad. Sci. USA 102 10837
[47] Friedman A J, Chan A, De Luca A and Chalker J T 2019 Phys. Rev. Lett. 123 210603
[48] Kozin V K and Kyriienko O 2019 Phys. Rev. Lett. 123 210602
[49] Chen J, Li X Z, Zhang Q, Probert M I J, Pickard C J, Needs R J, Michaelides A and Wang E 2013 Nat. Commun. 4 2064
[50] Boyer L L 1985 Phase Transit. 5 1
[51] Ciccotti G, Jacucci G and McDonald I R 1976 Phys. Rev. A 13 426
[52] Trullas J and Giro A 1990 J. Phys.: Condens. Matter 2 6643
[53] Tasseven Ç, Trullàs J, Alcaraz O, Silbert M and Giró A 1997 J. Chem. Phys. 106 7286
[54] Giguere P A 1979 J. Chem. Educ. 56 571
[55] Chau R, Mitchell A C, Minich R W and Nellis W J 2001 J. Chem. Phys. 114 1361
[56] Holmes N C, Nellis W J, Graham W B and Walrafen G E 1985 Phys. Rev. Lett. 55 2433
[57] Goldman N, Reed E J, Kuo I F W, Fried L E, Mundy C J and Curioni A 2009 J. Chem. Phys. 130 124517
[58] Yakushev V V, Postnov V I, Fortov V E and Yakysheva T I 2010 J. Exp. Theor. Phys. 90 617
[59] Hull S 2004 Rep. Prog. Phys. 67 1233
[60] Millot M, Hamel S, Rygg J R, Celliers P M, Collins G W, Coppari F, Fratanduono D E, Jeanloz R, Swift D C and Eggert J H 2018 Nat. Phys. 14 297
[61] Boyce J and Huberman B 1979 Phys. Rep. 51 189