# Simulating plastic ice VII with the data-driven many-body MB-pol potential

Francesco Guidarelli Mattioli, Xuanyu Zhu, Francesco Paesani, Livia E Bove,

John Russo, Francesco Sciortino

► To cite this version:

Francesco Guidarelli Mattioli, Xuanyu Zhu, Francesco Paesani, Livia E Bove, John Russo, et al. Simulating plastic ice VII with the data-driven many-body MB-pol potential. The Journal of Chemical Physics, 2025, 163 (20), 10.1063/5.0296428. [hal-05404161](https://hal-sorbonne-universite.fr/hal-05404161v1)

HAL Id: hal-05404161

https://hal.sorbonne-universite.fr/hal-05404161v1

Submitted on 8 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Simulating Plastic Ice VII with the Data-Driven Many-Body MB-pol Potential

Francesco Guidarelli Mattioli

*Department of Physics, Sapienza University of Rome, Piazzale Aldo Moro 5, I-00185 Roma, Italy*

Xuanyu Zhu, Francesco Paesani

*Department of Chemistry and Biochemistry, University of California, San Diego, La Jolla, California 92093, United States*

Livia Bove

*Department of Physics, Sapienza University of Rome, Piazzale Aldo Moro 5, I-00185 Roma, Italy and Sorbonne Université, UMR CNRS 7590, Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie (IMPMC), 5 Place Jussieu, 75005 Paris, France*

John Russo and Francesco Sciortino

*Physics Department, Sapienza University of Rome, Piazzale Aldo Moro 5, I-00185 Roma, Italy*

*(Dated: December 8, 2025)*

Computational studies using effective potentials that treat water as a rigid molecule were the first to predict the existence of a plastic phase of ice VII. Fifteen years later, this prediction was confirmed experimentally: in the plastic phase, water molecules retain the translational order of the body-centered cubic lattice while undergoing jump-like rotational motion typical of the liquid. Here, we show that this plastic phase also emerges in simulations employing the data-driven many-body MB-pol potential. The phase appears at high pressure, sandwiched between the liquid and ice VII phases. Our results confirm the first-order nature of both the ice VII–plastic and plastic–liquid transitions. Moreover, we find that ice VII exhibits rotational motion of hydrogen atoms on the nanosecond timescale near the melting line, although this motion is qualitatively distinct from that observed in the plastic phase. These results support the existence of the plastic phase beyond rigid-molecule models and provide new insight into the rotational dynamics of ice under extreme conditions.

# I. INTRODUCTION

Water continues to challenge our understanding of condensed matter with its rich and often counterintuitive phase behavior [1, 2], both in the liquid and in the crystal regions [3]. Over the past decades, the exploration of water's phase diagram under extreme thermodynamic conditions has revealed a multitude of crystalline and amorphous forms, each characterized by distinct structural, dynamical, and thermodynamic features. Among these, ice VII [4], a high-pressure polymorph of ice stable above 2 GPa, occupies a central place due to its relevance in planetary interiors, particularly for large icy moons and giant icy planets [5–8]. Structurally, ice VII consists of a body-centred cubic lattice of oxygen atoms, with each water molecule hydrogen-bonded to four neighbors in a disordered hydrogen configuration that distinguishes it from ordered counterparts such as ice VIII. The positions of the oxygens can be visualized as two interpenetrating diamond cubic latices [9].

Computational studies based on classical rigid water models — from the group of Hideki Tanaka [10–12] and of Carlos Vega [13, 14] — have been first in predicting the possible emergence of a plastic phase of ice VII, in which water molecules retain the translational order of the crystalline lattice while undergoing quasi-free rotational motion. In this phase, the hydrogen-bond network is topologically distinct from ice VII [?]. Free

energy calculations showed the thermodynamic stability of the plastic phase above a model dependent triple point (among ice VII, plastic and liquid phase) located below 10 GPa. Further studies examined the structure and dynamics of water plastic crystals [15] as well as their dielectric properties [16]. The potential existence of a water bcc plastic ice has also recently been supported by ab initio [17] and reactive force fields [18].

A different prediction was later on proposed in a numerically study based on a deep neural network potential trained on SCAN density functional calculations. This study [19] did not find evidence of a plastic phase, suggesting instead that above 10 GPa ice VII transforms into ice VII' and VII", ionic phases characterized by partial molecular dissociation and breaking of the ice rules. Experimental dielectric measurements have also suggested that at 10 GPa hydrogen atoms undergo a hopping translational dominant dynamics [20]. A plastic phase was also not found in simulations performed with a neural network potential [21]. Until recently the plastic phase had eluded experimental confirmation, largely due to the challenges in probing molecular orientations under the required high-pressure and high-temperature conditions. However, recent experimental observations [22] based on quasi-elastic neutron scattering [23] (which probe the dynamics of hydrogen atoms) have now confirmed the existence of a plastic state in the 5-6 GPa and 480-550 K region, in which molecules exhibits jump rotational
dynamics on the picosecond time scale. While provid-
ing direct evidence of a plastic phase, the experimental
study focused on dynamical observables, and thus could
not provide conclusive results on the thermodynamic na-
ture of the transition, e.g. whether it is continuous or
discontinuous.

Here, we revisit the problem of the existence of a ther-
modynamically stable plastic phase using the data-driven
many-body MB-pol potential. MB-pol accurately pre-
dicts the properties of liquid water across a wide range
of thermodynamic conditions [24, 25]. Below 1 GPa, it
is able to reproduce the stability field of all known crys-
talline forms of water [26]. It has also been shown that
MB-pol predicts a liquid-liquid critical point around 200
K and 1200 bar [27], in agreement with extrapolations
based on experimental data for the volume of liquid water
measured over a wide pressure-temperature range.[28]
The MB-pol potential represents a significant advance
in the modeling of water, being based on high-level elec-
tronic structure calculations, which do not suffer from the
approximations intrinsic in the density functional theory,
and its many-body expansion of the energy is rigorously
derived.[29–31] By construction, MB-pol neglects molec-
ular dissociation, limiting its applicability to non-ionic
phases. Unlike the previously studied TIP5P [32] and
TIP4P-2005 [33] models, MB-pol retains the flexibility
of water molecules. Recently, MB-pol has been refor-
mulated using permutationally invariant Fourier series
(PIFSs) as its functional basis.[34] This new formulation,
referred to as MB-pol(PIFS), improves the numerical ro-
bustness of the potential energy function and mitigates
the occurrence of catastrophically negative energies (so-
called “holes”), which may become relevant at high pres-
sures.

In the present study, we used MB-pol(PIFS) to per-
form simulations over a wide range of temperatures and
pressures, encompassing the conditions under which the
transition from fully hydrogen-disordered ice VII to the
plastic phase has been observed experimentally. We focus
on characterizing the equation of state and the dynam-
ics in the different phases and evaluating the pressure-
temperature coexistence lines between liquid, plastic and
ice VII phases. We find that MB-pol(PIFS) predicts the
existence of a plastic phase with a density intermedi-
ate between that of the crystal and the liquid, support-
ing the first-order nature of the ice VII-plastic transi-
tion. We observe that ice VII also exhibits rotational
dynamics near its melting line, but with important dif-
ferences compared to the plastic phase. First, rotational
motion in ice VII occurs on the nanosecond timescale,
significantly slower than the picosecond timescale ob-
served in the plastic phase and it is limited to 180° rota-
tions around the dipole vector. Second, in ice VII, rota-
tions are activated events between the stable configura-
tions where each molecule is bonded always to the same
diamond-cubic lattice. In contrast, the plastic phase fea-
tures stable configurations in which molecules can tran-
siently form hydrogen bonds across both interpenetrating

sublattices. Finally, we discuss the accuracy of the cur-
rent version of the MB-pol(PIFS) potential and examine
the implications of the commonly accepted definition of
“chemical accuracy” (e.g., 1 kcal/mol), particularly in
cases where relevant free-energy differences are on the
order of a few tenths of a kcal/mol. We also outline pos-
sible future improvements of MB-pol(PIFS).

II. METHODS

A. The MB-pol(PIFS/A1) potential

In this work we use the MB-pol(PIFS/A1) poten-
tial, referring to the PEF/A1 version presented in the
Supporting Information of Ref. 34. MB-pol(PIFS/A1)
is a variant of the data-driven many-body MB-pol
potential[29–31]. This potential is based on machine-
learned short-range *n*-body terms (with *n* ≤ 3) learned
from *ab initio* data calculated at the coupled cluster level
of theory, which are integrated with a physics-based term
representing 2-body dispersion and many-body electro-
statics at all molecular separations. MB-pol(PIFS/A1)
improves upon MB-pol by ensuring stability at high pres-
sures while maintaining the same robustness under am-
bient conditions.[34] For a detailed description of MB-
pol(PIFS/A1) and an assessment of its accuracy, the in-
terested reader is referred to Ref. [34].

B. Computational Details

We investigated the high-pressure phase diagram of
water by performing MD simulations in the isothermal-
isobaric (NPT: constant number of particles, pressure,
and temperature) ensemble at different pressures and
temperatures using the LAMMPS package[35] patched
with MBX[36, 37]. We started by generating an ice
VII crystal configuration of 864 water molecule with
GenIce [9] (replicating the unit box by 6-3-3 factors, pro-
ducing a box of approximately side 4-2-2 nm). We simu-
lated for at least 40 ps with a time step of 0.2 fs, running
with the openMP version of MBX on 40 cores. We ob-
served the spontaneous melting of ice VII to liquid water
at 5.5 GPa and 7 GPa by increasing the temperature. At
higher pressures (9 GPa and 11 GPa), we first observed
the spontaneous transition from ice VII to the plastic
phase, followed by melting of the plastic phase into the
liquid at higher temperatures. We performed direct coex-
istence simulations for liquid-ice VII and liquid-plastic ice
phase in the NPT ensemble to unveil the phase diagram.
To prepare the coexistence box, we first equilibrated the
liquid and ice phases using NPT simulations (each phase
being composed by 864 water molecules) and calculate
the corresponding equilibrium densities. We then scaled
the liquid box such as the *zy* side perfectly overlaps with
the ice *zy* side, and scaled the *z* axis of the liquid box
to preserve the correct liquid density. The interface be-
tween the liquid and ice was further equilibrated with a
Nosé-Hoover barostat applied to the *x* axis, keeping the
number of ice molecules fixed. Throughout these simu-
lations, we monitored ice growth by calculating a local
bond order parameter [38]. We also computed the coexis-
tence line by integrating the Clausis-Clapeyron equation
(following Ref. 39) with a pressure step of 0.25 GPa and
using the direct coexistence pressure and temperature
values as initial conditions.

C. Assessment of MB-pol(PIFS/A1)

Although MB-pol(PIFS/A1) preserves, and in some
cases improves, the accuracy of MB-pol for gas-phase
clusters and liquid water at ambient pressure, its per-
formance for high-pressure phases of water has not been
fully assessed. The simulations presented in this work re-
veal discrepancies with experimental observations, which
are discussed in this section. One notable issue is the MB-
pol(PIFS/A1)'s tendency to favor proton ordering at ele-
vated pressures. As a result, extended classical MD simu-
lations, particularly those near the melting temperature,
often evolve from the disordered ice VII phase toward
the fully proton-ordered ice VIII phase. This behavior
indicates that, under such conditions, classical MD sim-
ulations with MB-pol(PIFS/A1) predict ice VIII to be
the thermodynamically stable phase.

Our analysis suggests that this discrepancy likely arises
from *n*-body contributions with *n* > 3, which are not
quantitatively reproduced by many-body polarization.
Specifically, like MB-pol, the MB-pol(PIFS/A1) poten-
tial includes machine-learned corrections up to the 3-
body term, while 4-body and higher-order contributions
are approximated using a classical polarizable model. We
found that, despite remaining within “chemical accu-
racy,” the errors in the 4-body term become compara-
ble to the relevant free-energy differences at high pres-
sure, leading to a slight overestimation of both the lat-
tice energy (by ~0.5 kcal/mol) and the unit cell vol-
ume of ice VII. At high pressure, the Gibbs free energy
difference between ice VII and ice VIII (Δ*G*(*P*,*T*) =
Δ*E* + *P*Δ*V* − *T*Δ*S*) is highly sensitive to small vari-
ations (below ~0.5 kcal/mol) in all three contributing
terms, due to the elevated pressure and the relatively
small entropic differences between the two phases. As a
result, even deviations significantly smaller than “chemi-
cal accuracy” in MB-pol(PIFS/A1)’s predictions can in-
troduce noticeable uncertainty in the computed location
of the phase coexistence line. These observations sug-
gest that incorporating an explicit 4-body PIFS correc-
tion into MB-pol(PIFS/A1) could be a valuable direction
for future development. Additionally, it may be worth in-
vestigating nuclear quantum effects, as their impact on
enthalpy (and thus volume) and entropy could remain
non-negligible due to the sensitivity of free-energy differ-
ences, as previously observed in water clusters [40–44].

Numerical analyses of the free-energy differences associ-
ated with the ice VII–ice VIII equilibrium are provided
in the Supporting Information.

Nevertheless, for the purposes of this study, we have
focused exclusively on the metastable, disordered ice
VII phase. Owing to its long lifetime on the simula-
tion timescales considered here (nanoseconds), ice VII
remains well-defined and enables a meaningful analysis
of its phase behavior. We explicitly verified that all sim-
ulations, both bulk ice VII and direct coexistence setups
with the liquid, maintain proton disorder throughout
their duration, confirming that the observed transitions
are not artifacts of proton ordering.

It should also be noted that we encountered numerical
instabilities in some MD simulations, particularly in the
high-density plastic and liquid phases. These instabili-
ties suggest that MB-pol(PIFS/A1) occasionally predicts
dissociation events, which violate the assumptions un-
derlying the many-body framework on which the model
is based. It remains unclear whether these dissociation
events reflect physical phenomena or stem from numer-
ical artifacts. Further investigation is needed to assess
their sensitivity to possible sources of inaccuracy, includ-
ing the omission of 4-body energy corrections, nuclear
quantum effects, errors introduced by the numerical fit-
ting in MB-pol(PIFS/A1), and possible intrinsic inac-
curacies in the coupled cluster *n*-body energies that are
used for training.

III. RESULTS

1. *Thermodynamic differences*

Fig. 1(a) shows the temperature dependence of the
density at four different pressures as calculated from
*NPT* simulations. For each pressure, the results are
shown for simulations starting from an initial configu-
ration in the ice VII phase, from a liquid phase configu-
ration and (when existing) from a plastic configuration.
At pressures of ~9 GPa and above a temperature of ~600
K, ice VII transformed spontaneously (see Fig. 1(b)) into
a plastic phase, which remained stable on the timescale
of our simulations. As for the TIP4P/2005 model, the
transition was accompanied by a finite change in density
(~ 4% change with respect to the density of ice VII),
supporting its first order nature. At higher temperatures,
the plastic phase spontaneously transformed into a liquid
phase. A similar phenomenon was observed at 11 GPa,
with again a discontinuous change of density at the two
transitions. *NPT* simulations, which were initiated in
the plastic phase, remained in the plastic phase in a fi-
nite temperature interval, during the probed simulation
time scale.

Bulk simulations with periodic boundary conditions
tend to overestimate the stability range of any thermo-
dynamic phase, as they lack interfaces, heterogeneous
nucleation sites, and are affected by intrinsic barriers
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 1. (a) Temperature dependence of the density $\rho$ for the four studied pressures. Empty symbols represent the liquid density while the filled symbols follow the melting procedure and the dashed lines represent cooling runs, initialized from plastic phase. At P 11 GPa (blue) the melting curve has been calculated up to plastic transition while following increase in temperature made the potential instable. The liquid density at 11 GPa (blue dashed line) has been computed started from liquid bulk and reported to visualize the density drop between plastic and liquid. (b) Time evolution of the density during NPT simulations starting from an ice VII initial configuration at P 9 GPa. At low temperature (blue and green lines) the ice VII is preserved while increasing the temperature (brown and violet lines) the ice melts in the plastic phase. At higher temperature (cyan, magenta and orange lines) the plastic melts in the liquid phase. Initial time shift indicates that the starting configuration of the run has been selected from the equilibration of the previous run at lower temperature.</figcaption></figure>

to homogeneous nucleation. Consequently, although the abrupt density jump shown in Fig.1 suggests the existence of an intermediate phase between ice VII and the liquid, the temperature at which this change occurs is not thermodynamically meaningful, nor does it confirm the phase's thermodynamic stability. To accurately determine the transition temperature, we performed direct coexistence simulations [45] between the liquid and crystalline phases. At pressures where a plastic phase was found to emerge, we also simulated the coexistence between the plastic phase and the liquid at various temperatures, aiming to identify the point where the number of crystalline molecules fluctuates around 50%.

We defined a molecule as crystalline (either ice VII or plastic ice) evaluating $Q_6$, the rotational invariant first introduced by Steinhardt[46] and optimized for different crystal phase by Lechner and Dellago.[47] Fig. 2 shows the ice VII–liquid direct coexistence results for $P = 7$ GPa, suggesting that the coexistence temperature at this pressure is between 495 K ±5 K. Similar results were obtained at all studied temperature providing a first estimate of the phase diagram.

<figure><img src="image_4.png" /><figcaption>FIG. 2. Evolution of the fraction of crystalline molecules during a direct coexistence simulation of Ice VII and liquid water. A molecule is classified as crystalline based on the value of $l = 6$ rotational invariant [47].</figcaption></figure>

At $P = 9$ GPa, direct coexistence simulations between the liquid and the plastic phase always evolved toward the liquid, indicating that the plastic phase is not thermodynamically stable at this pressure. Although the plastic phase could be generated by heating ice VII at 9 GPa, it remained metastable with respect to the liquid across all temperatures, i.e. more stable than ice VII, but not more stable than the liquid.

In contrast, at $P = 11$ GPa, we found clear evidence that the plastic phase is thermodynamically more stable than the liquid. Specifically, in direct-coexistence simulations between liquid water and ice VII at $T = 635$ K, the crystal transformed into a plastic phase, despite the lat-
tice spacing being fixed in the orthogonal direction (see Fig. 3). At the same temperature, simulations of plastic-liquid coexistence showed that the plastic phase remained stable. Together, these observations demonstrated that at this temperature, the plastic phase is the most stable among the three phases (ice VII, liquid, and plastic) confirming its thermodynamic stability at $P = 11$ GPa and $T = 635$ K.

<figure><img src="image_5.png" /><figcaption>FIG. 3. Time evolution of the potential energy, density and fraction of crystalline molecules (as detected by $Q_6$) during a direct coexistence between liquid and crystal at $P=11$ GPa and $T=635$ K. Note that while the number of molecules in a crystalline environment remains similar, the density constantly decreases and the potential energy increases, signaling that Ice VII is turning into a plastic phase.</figcaption></figure>

To accurately trace the coexistence lines, we employed thermodynamic integration by numerically integrating the Clausius–Clapeyron equation, starting from the most reliable direct coexistence point at $P = 7$ GPa for the ice VII–liquid transition, as described in Ref. [39]. Similarly, we determined the plastic–liquid coexistence curve, initiating the integration from a direct coexistence result at $P = 11$ GPa. These two coexistence curves, illustrated in Fig. 4, intersect at a point in the $P-T$ plane. From this intersection, we initiated the Kofke integration [39] to map the coexistence boundary between plastic and ice VII. For this particular transition, direct coexistence methods are not applicable due to a mismatch in the lattice parameters of the two phases.

Fig. 4 shows that the plastic phase emerges at pressures above the triple point, around $P = 11$ GPa, where all three phases (ice VII, plastic, and liquid) coexist. Beyond this point, the plastic phase occupies the region between ice VII and the liquid in the phase diagram. The resulting coexistence lines confirm the first-order character of both transitions and clearly delineate the stability region of the plastic phase at high pressure. Experimentally the plastic phase is observed above $P = 5$ GPa [22]. The shift in pressure and temperature (Fig.S5) for the

location of plastic ice is consistent with the shift in the ice VII melting line, with respect to the experimental measurement [18].

<figure><img src="image_6.png" /><figcaption>FIG. 4. Phase diagram of the MB-pol(PIFS/A1) model. The pink symbols indicate a bracketing of the coexistence region as calculated via direct coexistence methods, where both phases appear metastable (during the simulation time) within the temperature range of the error bar. The lines are the results of Clausius Clapeyron integration. Dashed lines represent the upper bound for plastic stability region. Star black symbols are the lower and higher bound for the triple point estimation, respectively at $T = 570$ K, $P = 9.4$ GPa and $T = 590$ K, $P = 10.1$ GPa.</figcaption></figure>

## 2. Structural and dynamical differences between the phases

Having established the thermodynamic stability region of the VII plastic phase as predicted by the MB-pol(PIFS/A1) potential, we examined the structural and dynamic properties of the different phases.

Fig. 5 shows the radial distribution functions $g(r)$, separated in components (OO, HH, OH). Ice VII shows well-defined peaks in all partial $g(r)$. The plastic phase retains weaker and less defined correlations in $g_{OO}(r)$ signaling that the bcc lattice of the oxygen is distorted going from the ice VII to the plastic phase. Except for the first peak of the $g_{OH}(r)$ (that stands for the molecule integrity), the spatial distribution between oxygen and hydrogen of different water molecules is similar between plastic and liquid phase and differs from crystal phase. Together with observation done on $g_{OO}(r)$ this shows that the spatial order of hydrogen in plastic phase is more liquid-like than crystal-like. In agreement with simulations using rigid water models [15], very marked changes are seen in
$g_{HH}(r)$, which is quite similar for both plastic ice and the liquid phase.

The mean square displacement of oxygen atoms in the ice VII phase (on a bcc lattice) is primarily governed by vibrational motion around their equilibrium lattice positions. Similarly, the mean square displacement of the hydrogen atoms in ice VII is also expected to be governed by their vibrational motion. Interestingly, at 5.5 GPa we observed a net diffusion of the H atoms, significantly larger than the vibrational component. This diffusion process, which must be connected to changes in the molecular orientation, being MB-pol(PIFS/A1) potential non dissociable, indicates the presence of rotational jumps that rearrange the positions of the H atoms within the hydrogen-bond network. This rotational motion is not negligible even when the crystal is the thermodynamically stable phase. An analysis of the trajectories shows that these jump rotations take place around the C<sub>2v</sub> axis of the molecules (the dipole vector direction). In this ~180° rotation, the two hydrogen atoms of the molecule swap their positions, preserving the oxygen-oxygen hydrogen-bond connectivity. These slow hydrogen motions were observed, within the time scale of our simulations, also at higher pressures. Fig. 6 shows the mean square displacement for both hydrogen and oxygen atoms at P = 7 GPa for different temperatures in the stable ice VII phase. The data in Fig. 6 suggest that, with a switching timescale governed by temperature, the molecules undergo rotational motion and the hydrogen atoms switch their acceptor oxygens. Scaling along the time axis the different mean-square displacements (not shown) suggests that the slowing down of the dynamics is controlled by an Arrhenius process with an activation energy of about 1200 K at this pressure (Fig.S6).

Although ice VII can exhibit rotational dynamics, this behavior is fundamentally distinct from that observed in the plastic phase. In ice VII, molecular reorientation is thermally activated (Fig.S6) and corresponds to a flip of the molecule around the direction of the dipole which keeps the two interpenetrated diamond sub-lattices distinct. On the other hand, in the plastic phase the orientational motion is described by a jump-rotor mechanism nearly independent of temperature and pressure on picosecond-timescale [22], which cancels the distinction between the two sub-lattices. To better visualize the very shorter relaxation time of hydrogen rotation in plastic phase with respect to ice VII phase, we report in Fig.7 the mean square displacement of liquid, plastic and ice VII at P 8 GPa. On ps timescale the hydrogen in the plastic phase reaches a plateau whereas the hydrogen in the ice VII phase does not diffuse. The plastic and ice VII phase thus differ both for relaxation time of hydrogen and sub-lattice bonding pattern. To illustrate the difference in sub-lattices bonding, we compare the inherent structures of ice VII under two sets of conditions: one in which no molecular reorientations are observed on the simulation timescale (P = 7GPa, T = 400K), and one in which they are (P = 7GPa, T = 475K). This

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 5. Radial distribution functions at 11 GPa and different temperatures. Thermodynamically stable and metastable states are shown when possible, to allow comparison between different phases at the same temperature.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 6. Mean square displacement for H and O atoms at 7 GPa and different temperatures in the ice VII phase. Despite the thermodynamic stable phase is Ice VII, the H atoms undergo a significant diffusional process, while the O atoms remains anchored around the lattice positions.</figcaption></figure>

comparison is presented in the top and middle panels of Fig. 8, which show the probability distribution of the azimuthal angle between the hydrogen-oxygen bond vector and the [001] crystallographic direction. In a perfect bcc lattice, hydrogen atoms can orient along four equivalent diagonal directions within the [001] plane. In both panels, the red curves represent the angular distribution in the thermal (finite-temperature) simulations, showing broadened peaks around the four crystallographic directions, indicative of thermally activated reorientations. The black curves, obtained from energy-minimized configurations (i.e., inherent structures), reveal sharp peaks located precisely at the bcc orientations. In this way, all vibrational component are effectively suppressed, revealing the structure of the system at its local energy minimum. This analysis confirms that, despite thermal motion, the system remains confined to well-defined local

<figure><img src="image_11.png" /><figcaption>FIG. 7. Mean Square Displacement of oxygen (circles) and hydrogen (star) at P 8 GPa. The three different phases (Ice VII, plastic ice and liquid) are annotated and depicted in violet, blue and turquoise respectively. The dynamical behaviour of plastic ice is signaled by a non diffusive regime for oxygen (blue circles) and a diffusive regime for hydrogen (blue stars) on ps scale. In the Ice VII phase both oxygen and hydrogen (violet circles and stars) do not diffuse. The linearity of liquid curves for both oxygen and hydrogen (turquoise circles and stars) indicate a over ps diffusion regime.</figcaption></figure>

minima corresponding to the ordered ice VII phase, i.e., no stable intermediate orientations emerge.

By contrast, the bottom panel of Fig. 8 displays the same analysis for the plastic phase. Here, both the thermal and inherent state distributions exhibit a broad range of orientations, without sharp localization. This indicates that in the plastic phase, orientational disorder is intrinsic and persists even after energy minimization. Furthermore, we find that the inherent structures include configurations in which the two interpenetrating sublattices of the bcc lattice become indistinguishable, namely, molecules can form hydrogen bonds with neighbors from both sublattices, highlighting the breakdown of the orientational order characteristic of ice VII.

## IV. CONCLUSIONS

In this study, we have shown that water at pressures above 9 GPa exhibits a plastic phase by simulating water molecules with the high-level description of MB-pol (PIFS/A1) potential. We find the phase diagram region where the water plastic phase is thermodynamically more stable than the ice VII and the liquid phases. We find that the rotational dynamics of hydrogen atoms in plastic ice occur on the picosecond timescale, in agreement
<figure><img src="image_12.png" /><figcaption>FIG. 8. Probability distribution of the azimuthal angle between hydrogen-oxygen distance and the [001] direction for ice VII at *P* = 7GPa and *T* = 400K (top), ice VII with high rotational diffusion at *P* = 7 GPa and *T* = 475 K (middle), and plastic phase at *P* = 11 GPa and *T* = 625 K (bottom). Red lines represent the probability distribution in the thermal state, while black lines in the inherent states, showing that Ice VII is the same state with and without rotational diffusion, differently from the plastic phase.</figcaption></figure>

with recent inelastic neutron scattering experiments and the analysis of inherent structure configurations confirms that, as in the TIP4P/2005 model, hydrogen atoms in the plastic phase can transiently form bonds with both interpenetrating sublattices. Additionally, we observe a much slower rotational process in ice VII which, while effectively decorrelating hydrogen positions on nanosecond time scale, but consistently maintains hydrogen bonding within the original sublattice. We show that both the transitions from Ice VII to plastic and from plastic to liquid phase are signaled by a jump in density, suggesting a first order nature of the transitions. The first order nature of the phase transitions has been found also by the tracing of the phase diagram where we observe a coexistence line kinks near the triple point. Hence we demonstrate the existence of the plastic phase using an accurate water potential, taking into account flexibility, rigorous many-body derivation and polarizability. We further disclose structural, dynamical and thermodynamical peculiarity of water in Ice VII and plastic phase, highlighting differences in the molecular mechanisms. At high den-

sities, particularly in the plastic and liquid phases, the MB-pol(PIFS/A1) potential exhibits some numerical instabilities. This may indicate the presence of “holes” in the potential energy surface due to incomplete training data or to the lacking of the dissociation modelling, even with the improvements introduced in the PIFS variant. We also note that MB-pol(PIFS) predicts a shift with respect experimental phase diagram possibly reflecting limitations in the representation of the high pressure configurational space, limitations in the treatment of higher-body interactions (n > 3), the neglect of nuclear quantum effects and possible inherent errors in the raw n-body energies computed at the coupled cluster level of theory. Future work will focus on understanding the necessary physical modelling advancement to reach a more accurate description of water at high pressure, ultimately aiming to unravel the intricate behaviour of water at extreme pressure conditions.

V. ACKNOWLEDGMENTS

FGM, JR and FS acknowledge support from MIUR-PRIN 2022JWAF7Y, Cineca ISCRAB initiative and ICSC-Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by the European Union “NextGenerationEU”). X.Z. and F.P. acknowledge support from the Air Force Office of Scientific Research under grant no. FA9550-20-1-0351, and allocations of computational resources from the Department of Defense High Performance Computing Modernization Program (HPCMP) and the Advanced Cyber infrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants nos. 2138259, 2138286, 2138307, 2137603, and 2138296. L.E.B. acknowledges the financial support by the European Union - NextGenerationEU (PRIN N. 2022NRBLPT), the ANR-23-CE30-0034 EXOTIC-ICE, and the Swiss National Fund (FNS) under Grant No. 212889.

VI. DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available within the article. The public domain LAMMPS code has been used to generate the trajectories for all different state points. Code for the integration of Clausius-Clapeyron, input file and representative trajectory have been provided on Zenodo at https://doi.org/10.5281/zenodo.17350531.

<footnote>[1] J. L. Finney, *Water: a very short introduction*, Vol. 440 (Oxford University Press, 2015).</footnote>

<footnote>[2] P. Gallo, K. Amann-Winkel, C. A. Angell, M. A. Anisimov, F. Caupin, C. Chakravarty, E. Lascaris, T. Loerting, A. Z. Panagiotopoulos, J. Russo, et al., *Chem. Rev.*</footnote>
116, 7463 (2016).

[3] C. G. Salzmann, J. Chem. Phys. **150** (2019).

[4] B. Kamb and B. L. Davis, Proceedings of the National Academy of Sciences **52**, 1433 (1964).

[5] B. Journaux, K. Kalousová, C. Sotin, G. Tobie, S. Vance, J. Saur, O. Bollengier, L. Noack, T. Rückriemen-Bez, T. Van Hoolst, et al., Space Science Reviews **216**, 1 (2020).

[6] M. Lingam and A. Loeb, Rev. Mod. Phys. **91**, 021002 (2019).

[7] M. C. Nixon and N. Madhusudhan, Monthly Notices of the Royal Astronomical Society **505**, 3414 (2021).

[8] T. Van Hoolst, G. Tobie, C. Vallat, N. Altobelli, L. Bruzzone, H. Cao, D. Dirkx, A. Genova, H. Hussmann, L. Iess, J. Kimura, K. Khurana, A. Lucchetti, G. Mitri, W. Moore, J. Saur, A. Stark, A. Vorburger, M. Wieczorek, A. Aboudan, J. Bergman, F. Bovoło, D. Breuer, P. Cappuccio, L. Carrer, B. Cecconi, G. Choblet, F. De Marchi, M. Fayolle, A. Fienga, Y. Futaana, E. Hauber, W. Kofman, A. Kumamoto, V. Lainey, P. Molyneux, O. Mousis, J. Plaut, W. Puccio, K. Retherford, L. Roth, B. Seignovert, G. Steinbrügge, S. Thakur, P. Tortora, F. Tosi, M. Zannoni, S. Barabash, M. Dougherty, R. Gladstone, L. I. Gurvits, P. Hartogh, P. Palumbo, F. Poulet, J.-E. Wahlund, O. Grasset, and O. Witasse, Space Science Reviews **220**, 54 (2024).

[9] M. Matsumoto, T. Yagasaki, and H. Tanaka, J. Chem. Phys. **160**, 094101 (2024).

[10] Y. Takii, K. Koga, and H. Tanaka, J. Chem. Phys. **128** (2008).

[11] K. Himoto, M. Matsumoto, and H. Tanaka, Phys. Chem. Chem. Phys. **13**, 19876 (2011).

[12] K. Himoto, M. Matsumoto, and H. Tanaka, J. Phys. Soc. Japan **81**, SA023 (2012).

[13] J. Aragones, M. Conde, E. Noya, and C. Vega, Phys. Chem. Chem. Phys. **11**, 543 (2009).

[14] J. Aragones and C. Vega, J. Chem. Phys. **130** (2009).

[15] A. Henao, J. M. Salazar-Rios, E. Guardia, and L. C. Pardo, J. Chem. Phys. **154** (2021).

[16] I. Skarmoutsos, S. Mossa, and E. Guardia, J. Chem. Phys. **150** (2019).

[17] J.-A. Hernandez and R. Caracas, J. Chem. Phys. **148** (2018).

[18] Y. Adachi and K. Koga, J. Chem. Phys. **153** (2020).

[19] L. Zhang, H. Wang, R. Car, and W. E, Phys. Rev. Lett. **126**, 236001 (2021).

[20] R. Yamane, K. Komatsu, and H. Kagi, Physical Review B **104**, 214304 (2021).

[21] A. Reinhardt, M. Bethkenhagen, F. Coppari, M. Millot, S. Hamel, and B. Cheng, Nature Communications **13**, 4707 (2022).

[22] M. Rescigno, A. Toffano, U. Ranieri, L. Andriambariarijaona, R. Gaal, S. Klotz, M. M. Koza, J. Ollivier, F. Martelli, J. Russo, et al., Nature, **1**, (2025).

[23] S. Klotz, T. Strässle, and L. E. Bove, Applied Physics Letters **103**, 193504 (2013).

[24] S. K. Reddy, S. C. Straight, P. Bajaj, C. Huy Pham, M. Riera, D. R. Moberg, M. A. Morales, C. Knight, A. W. Götz, and F. Paesani, J. Chem. Phys. **145** (2016).

[25] E. Palos, E. F. Bull-Vulpe, X. Zhu, H. Agnew, S. Gupta, S. Saha, and F. Paesani, J. Chem. Theory Comput. **20**, 9269 (2024).

[26] S. L. Bore and F. Paesani, Nature Communications **14**, 3349 (2023).

[27] F. Sciortino, Y. Zhai, S. Bore, and F. Paesani, Nature Physics **21**, 480 (2025).

[28] O. Mishima and T. Sumita, J. Phys. Chem. B **127**, 1414 (2023).

[29] V. Babin, C. Leforestier, and F. Paesani, J. Chem. Theory Comput. **9**, 5395 (2013).

[30] V. Babin, G. R. Medders, and F. Paesani, J. Chem. Theory Comput. **10**, 1599 (2014).

[31] G. R. Medders, V. Babin, and F. Paesani, J. Chem. Theory Comput. **10**, 2906 (2014).

[32] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. **79**, 926 (1983).

[33] J. L. Abascal and C. Vega, J. Chem. Phys. **123** (2005).

[34] X. Zhu and F. Paesani, J. Chem. Theory Comput. **21**, 6950–6963 (2025).

[35] 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., Computer Physics Communications **271**, 108171 (2022).

[36] M. Riera, C. Knight, E. F. Bull-Vulpe, X. Zhu, H. Agnew, D. G. A. Smith, A. C. Simmonett, and F. Paesani, J. Chem. Phys. **159**, 054802 (2023).

[37] S. Gupta, E. F. Bull-Vulpe, H. Agnew, S. Iyer, X. Zhu, R. Zhou, C. Knight, and F. Paesani, J. Chem. Theory Comput. **21**, 1838–1849 (2025).

[38] J. Russo, F. Romano, and H. Tanaka, Nature materials **13**, 733 (2014).

[39] D. A. Kofke, J. Chem. Phys. **98**, 4149 (1993).

[40] C. Pérez, M. T. Muckle, D. P. Zaleski, N. A. Seifert, B. Temelso, G. C. Shields, Z. Kisiel, and B. H. Pate, Science **336**, 897 (2012).

[41] Y. Wang, V. Babin, J. M. Bowman, and F. Paesani, J. Am. Chem. Soc. **134**, 11116 (2012).

[42] V. Babin and F. Paesani, Chem. Phys. Lett. **580**, 1 (2013).

[43] L. Evangelisti, C. Perez, S. Lobsiger, N. A. Seifert, D. P. Zaleski, B. Pate, Z. Kisiel, B. Temelso, and G. C. Shields, in *69th International Symposium on Molecular Spectroscopy* (2014) p. FD12.

[44] S. E. Brown, A. W. Götz, X. Cheng, R. P. Steele, V. A. Mandelshtam, and F. Paesani, J. Am. Chem. Soc. **139**, 7082 (2017).

[45] M. Conde, M. Rovere, and P. Gallo, J. Chem. Phys. **147** (2017).

[46] P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B **28**, 784 (1983).

[47] W. Lechner and C. Dellago, J. Chem. Phys. **129** (2008).