Computational Methods for Molecular Dynamics of Water at Low Temperatures

1. Introduction

Water is arguably the most studied liquid in science, yet its behavior continues to reveal surprising complexity, particularly at low temperatures. The thermodynamic and transport properties of water deviate increasingly from simple-liquid behavior on cooling, with the deviations already apparent below the temperature of maximum density (∼277 K) and becoming increasingly dramatic upon supercooling below the equilibrium melting point (273 K) [1,2]. These anomalies (i.e.: density, isothermal compressibility and heat capacity changes) have motivated decades of theoretical and experimental investigation [3,4,5].

Understanding water behavior at low temperatures is not merely of academic interest. For example: supercooled water plays critical roles in atmospheric chemistry, where it exists metastably in clouds at temperatures as low as 235 K [6].In biomedical sciences, supercooling preservation enables extended storage of organs, tissues, and cells for transplantation [7,8,9]. In food technology and material processing, controlling supercooled water states is essential for cryopreservation without ice crystal damage [10,11]. Amorphous solid water (ASW), formed by vapor deposition at cryogenic temperatures, is widely regarded as the dominant form of water in the universe, coating interstellar dust grains where it serves as a substrate for complex organic molecule (COM) formation [12,13,14,15]. In planetary science, water ice in various amorphous and crystalline forms is found on icy moons, comets, and Mars [16]. These diverse applications demand accurate theoretical descriptions, aimed to facilitate in silico experimentation, which has important adventages over wet lab experimentation as cost, time and reproducibility, among others.

Classical molecular dynamics (MD) simulations have provided valuable insights into water behavior, but they fundamentally cannot capture several phenomena critical to low-temperature systems. Chemical reactions, including the radical chemistry relevant to astrochemical environments, require explicit treatment of electronic structure. Moreover, the rigid non-polarizable classical water force fields neglect both electronic polarizability and nuclear quantum effects, which limits their accuracy for structural and thermodynamic properties of water at low temperatures [17].

Ab initio molecular dynamics (AIMD) simulations at low temperatures, face their own substantial challenges. First, the dynamics slow dramatically as temperature decreases, making adequate sampling increasingly difficult within accessible simulation timescales. Second, nuclear quantum effects (NQE), as zero-point energy, tunneling, and proton delocalization, become progressively more important as the thermal energy $k_{\mathrm{B}}T$ becomes comparable to vibrational zero-point energies. At 300 K, the ratio of hydrogen zero-point energy to thermal energy is about an order of magnitud, which increases when the temperature drops, indicating that classical treatment of nuclei introduces significant systematic errors [18]. Third, the accuracy of different density functionals for water remains a subject of ongoing investigation, with different approximations yielding qualitatively different predictions for structure and dynamics [19].

This review addresses these challenges systematically. We aim to provide practitioners with a comprehensive guide to performing reliable AIMD simulations of water in the 200–273 K temperature range. Section 2 briefly reviews the structural characteristics of water at low temperatures, including supercooled liquid water and amorphous ices, establishing the physical context and validation targets for simulations. Section 3 presents the principal AIMD approaches, from Born–Oppenheimer MD (BOMD) to path-integral MD (PIMD), and discusses how each addresses the specific challenges that water simulations face at ultra-low temperatures: slow dynamics and degraded sampling, the growing importance of nuclear quantum effects, and the sensitivity of structural and dynamical properties to the underlying potential energy surface. Section 4 addresses density functional selection for water. Section 5 addresses the central methodological challenge: nuclear quantum effects and their interplay with density functional selection.

2. Water Structure at Low Temperatures

First, it is essential to establish the physical phenomena that simulations must capture. This section reviews the structural evolution of water from ambient conditions through the supercooled regime to the glassy and crystalline states.

2.1. Supercooled Liquid Water

This metastable supercooled state has been accessed experimentally down to its homogeneous freezing limit of approximately 235 K in micrometre-sized samples [20], and to even lower temperatures using rapidly cooled nanodroplets [21] and ultrafast X-ray probing of evaporatively cooled microdroplets [22], the latter revealing metastable bulk liquid water down to 227 K.

The structure of supercooled water evolves continuously from ambient conditions, with an accelerating increase in ordering as temperature decreases. X-ray diffraction studies reveal that upon cooling, the first oxygen–oxygen radial distribution function (RDF) peak near 2.8 Å grows in height, while the running O–O coordination number exhibits a temperature-independent isosbestic point near 3.3 Å, and the second peak near 4.5 Å becomes increasingly pronounced below the compressibility-minimum temperature, reflecting enhanced tetrahedral coordination [23]. Measurements extending to 234.8 K show that intermediate-range correlations continue to strengthen upon supercooling, with the 4th and 5th peaks near 9 and 11 Å increasing in height and the 4th-peak position shifting to shorter distances [24]. At ambient conditions, water already exhibits nanometer-scale density fluctuations representing tetrahedral-like and hydrogen-bond-distorted structures [25]. These structural changes are often interpreted within a two-state framework, where water is viewed as a fluctuating mixture of locally tetrahedral, low-density liquid (LDL) and more disordered, high-density liquid (HDL) structural motifs [2]. More recent experiments reported observation of a liquid-liquid transition under pressure at 205 ± 10 K [26]. The temperature range between approximately 160 K and 232 K, often called the “no-man’s land”, has historically been inaccessible to bulk experiments due to rapid crystallization [22]. Below approximately 140 K, solid water has been obtained in different distinct glassy forms, from low-density amorphous ice (LDA) to high-density amorphous ice (HDA) [27,28]. However, our understanding of water’s structure below 235 K, whether it represents continuous evolution or discrete phase transitions, continues evolving, and computational simulations provide much of the structural insight for this experimentally challenging region.

2.2. Amorphous Ice

Amorphous ice phases are key benchmarks for low-temperature simulations. The principal forms relevant to this review are LDA, MDA, HDA, VHDA, and HGW.

2.2.1. Low-Density Amorphous Ice (LDA)

Formation temperature and conditions. LDA can be formed by several routes: (i) vapor deposition onto cold substrates below ∼120–130 K [29,30]; (ii) decompression of HDA below ∼150 K [28]; and (iii) hyperquenching of liquid aerosol droplets at cooling rates exceeding ${10}^5$ K/s to temperatures below ∼140 K [31]. In the interstellar medium, LDA forms via vapor deposition of water molecules onto cold dust grains (<130 K) under ultra-high vacuum conditions, making it the most abundant form of water in the universe [14,15].

Structure. The first peak of the oxygen–oxygen radial distribution function $g_{\mathrm{OO}}\left(r\right)$ is located at $r_{\mathrm{OO}}\approx 2.75$ Å [32,33]. LDA exhibits a locally tetrahedral hydrogen-bond network resembling ice I_h but lacking long-range crystallographic order [32].

Key remarks. LDA is considered the most common form of water in the universe, constituting the amorphous solid water coatings on interstellar dust grains [14,15]. The glass transition temperature is $T_g\approx 136$ K, as determined calorimetrically [34]. LDA is metastable and crystallizes to stacking-disordered ice I_sd upon heating at about ∼200 K at ambient pressure [35,36]. LDA has a density of approximately 0.94 g cm⁻³ at 77–80 K and ambient pressure. This value is derived from reported atomic number densities of 0.094 atoms/ų at 80 K [32] and 0.0937 atoms/ų [33]; for H2O, 0.094 atoms/ų converts to ∼0.937 g cm⁻³.

2.2.2. Medium-Density Amorphous Ice (MDA)

Density. MDA has a density of $1.06\pm 0.06$ g cm⁻³ at 77 K and ambient pressure [37], falling within the density gap previously thought to be empty between LDA (∼0.94 g cm⁻³) and HDA (∼1.17 g cm⁻³).

Formation temperature and conditions. MDA is produced by ball milling hexagonal ice (ice I_h) at 77 K using stainless-steel balls under cryogenic conditions [37]. Molecular dynamics simulations suggest that isobaric cooling at $P\approx 125$ MPa can produce structurally similar amorphous states [37]. More recent simulations identify MDA as a shear-driven amorphous ice that can also be formed by shearing LDA or HDA at controlled shear rates [38].

Structure. The O–O radial distribution function of MDA shows peaks intermediate between those of LDA and HDA [37]. Importantly, MDA is not a simple mixture of LDA and HDA; its structure factor exhibits distinct features that cannot be reproduced by linear combination of LDA and HDA patterns [37].

Key remarks. MDA was first reported in 2023 by Rosu-Finsen et al. [37]. Its density near that of liquid water (∼1.00 g cm⁻³) raises the possibility that MDA may represent the true glassy state of liquid water, though this remains debated [37,38]. The discovery suggests a continuum of amorphous ice states rather than discrete phases separated by sharp density gaps.

2.2.3. High-Density Amorphous Ice (HDA)

Density. HDA has a density of approximately 1.17 g cm⁻³ at 77 K and ambient pressure [39], though the originally reported value upon pressure-amorphization was ∼1.19 g cm⁻³ [28]. Values in the range 1.15–1.17 g cm⁻³ are reported depending on the specific preparation protocol and degree of structural relaxation [39,40].

Formation temperature and conditions. HDA is formed by: (i) pressure-induced amorphization of ice I_h at 77 K under ∼1 GPa [28]; or (ii) compression of LDA at temperatures below ∼150 K, with the transition occurring at ∼0.5–0.6 GPa [28,41]. These transitions occur below the crystallization temperature of amorphous ice (∼150 K at ambient pressure) [39]. Importantly, unannealed HDA (uHDA) and equilibrated (expanded) HDA (eHDA) differ in their structural and thermal properties, reflecting the sensitivity of the amorphous state to thermal history [42,43].

Structure. The first O–O coordination shell is located at $r_{\mathrm{OO}}\approx 2.8$ Å, with approximately 4 nearest oxygen neighbors in the first shell [33,39]. The structure features interpenetrating hydrogen-bond networks in which the second coordination shell collapses inward compared to the open tetrahedral arrangement of LDA [33,39]. Neutron diffraction studies of the five amorphous ices at ambient pressure show that HDA is structurally distinct from both LDA and VHDA in its intermediate-range order [32].

Key remarks. HDA is approximately 20–25% denser than LDA. Once formed, it can be recovered to ambient pressure at 77 K and stored indefinitely [28,39]. The apparently first-order LDA–HDA transition is widely interpreted as evidence supporting the liquid–liquid critical point hypothesis for water [1,44].

2.2.4. Very High-Density Amorphous Ice (VHDA)

Density. VHDA has a density of $1.25\pm 0.01$ g cm⁻³ at 77 K and ambient pressure [45], approximately 9% higher than that of HDA. Values up to 1.26 g cm⁻³ have been reported depending on preparation pressure [39].

Formation temperature and conditions. VHDA is prepared by isobaric annealing of HDA at pressures of 1.1–1.9 GPa and temperatures of 160–175 K, followed by quenching to 77 K [39,45,46]. VHDA forms at pressures >0.8 GPa [39]. A stepwise transformation sequence LDA → HDA (at ∼0.5 GPa) → VHDA (at ∼0.95 GPa) has been observed at 125 K [41].

Structure. The first O–O coordination shell remains near ∼2.8 Å [45,47]. The Raman spectrum of VHDA differs clearly from HDA, with the hydrogen-bonded O–O distance increasing from 2.82 Å in HDA to 2.85 Å in VHDA [45]. Under pressure (2 GPa), the second shell collapses to ∼3.0 Å with 5 neighbors [48], reflecting different interstitial occupancy compared to HDA.

Key remarks. VHDA prepared at high pressures and temperatures (1.9 GPa, 175 K) represents the amorphous state least affected by nanocrystalline remnants [46]. The activation barrier for the VHDA → HDA back-conversion is higher than for HDA → LDA [39]. The crystallization temperature of VHDA increases by up to 4 K when prepared at 1.9 GPa compared to 1.1 GPa, reflecting improved structural homogeneity [46].

2.2.5. Hyperquenched Glassy Water (HGW)

Density. HGW has a density of approximately 0.94 g cm⁻³, similar to LDA [31].

Formation temperature and conditions. HGW is produced by: (i) rapid spraying of micrometer-sized (∼5 $\mu$m) water droplets into a cryogenic medium (e.g., liquid propane) at ∼80 K [31]; (ii) deposition of fine droplets onto a cold substrate at ≤80 K; or (iii) cooling of water in capillary tubes (∼100 $\mu$m diameter) with liquid helium at 4.2 K. The required cooling rates are ${10}^5$–${10}^7$ K/s [31,49].

Structure. HGW is structurally very similar to LDA and ASW, particularly after annealing [34]. It possesses a fully hydrogen-bonded, open tetrahedral network characteristic of low-density amorphous forms [50].

Key remarks. The glass transition temperature of HGW is $T_g\approx 136$ K according to Johari, Hallbrucker, and Mayer [34,50]. However, this value is debated; Velikov, Borick, and Angell argued based on enthalpy relaxation scaling that $T_g=165\pm 5$ K [49], while Kohl et al. presented DSC evidence supporting the 136 K assignment [51]. HGW is kinetically stable at 77 K and can be stored for years. Above $T_g$, HGW behaves as an ultraviscous liquid before crystallizing near ∼150–154 K [50,51].

2.2.6. Summary of Amorphous Ice Properties

Table 1 provides a consolidated comparison of the principal amorphous ice phases.

2.3. Crystalline Ice Phases

For completeness and to provide validation targets for simulations, Table 2 summarizes the key crystalline ice phases relevant to low-temperature conditions at atmospheric and moderate pressures.

Ice I_h (hexagonal ice) is the thermodynamically stable form of ice I at ambient pressure, whereas ice I_c is metastable cubic ice [57]. In practice, what has often been reported as “cubic ice” is typically stacking-disordered ice I_sd, containing cubic sequences interlaced with hexagonal sequences [57], which results in variable density and $r_{\mathrm{OO}}$.

X-ray diffraction measurements on freezing micron-sized water droplets show a droplet-size dependence, where smaller droplets freeze dominantly to cubic ice with stacking faults and larger droplets favour hexagonal ice [64].

Ice XI is the proton-ordered form of ice I_h and is thermodynamically stable below ∼72 K (near ambient pressure) [65].

3. Ab Initio Molecular Dynamics Methods for Water

The theoretical description of atomic and molecular systems based on the Schrödinger equation emerged from early theoretical developments of electron and nuclear behavior. The Schrödinger equation molecular Hamiltonian comprises terms for nuclear kinetic energy, electronic kinetic energy, electron-nuclear attraction, electron-electron repulsion, and nuclear-nuclear repulsion. Solving this equation exactly for all but the simplest systems presents formidable computational challenges, given the exact mathematical form of the wavefunction is unknown. Thus, there are some common approximation methods that address the tradeoff between resource consumption and accuracy. This section presents the principal AIMD methods used for water simulations in cold regimes, emphasizing the fundamental origin of their limitations or adventages.

3.1. Born-Oppenheimer Molecular Dynamics (BOMD)

In BOMD, the electronic structure is converged self-consistently at each nuclear configuration, yielding forces from the ground-state energy surface [66]. Nuclei evolve according to Newton’s equations:

$$M_I{\ddot{\mathbf{R}}}_I=-{\nabla }_{{\mathbf{R}}_I}{E}_0\left[\left\{\mathbf{R}\right\}\right]$$

(1)

where $M_I$ and ${\mathbf{R}}_I$ are the mass and position of nucleus I, and $E_0$ is the ground-state electronic energy.

The importance of these nuclear quantum effects can be quantified by the ratio of zero-point energy to thermal energy, $\hslash \omega /2k_{\mathrm{B}}T$. For the O–H stretching mode in water ($\omega \approx 3600$ cm⁻¹ [67]), the zero-point energy $E_0={\displaystyle \frac{1}{2}}\hslash \omega \approx 5.1$ kcal mol⁻¹, which already exceeds the thermal energy $k_{\mathrm{B}}T$ by an order of magnitude at room temperature [18]. As temperature decreases, this ratio grows further: at 200 K, $k_{\mathrm{B}}T\approx 0.4$ kcal mol⁻¹, making the ZPE approximately 13 times larger than the thermal energy available to the O–H coordinate.

BOMD performs a clean separation of electronic and nuclear timescales, preserving a rigorous energy conservation trhough the simulation, in contrast to CPMD.

However, BOMD is unable to compute zero-point energy or nuclear tunneling; errors grow systematically as T decreases and the de Broglie thermal wavelenght $\Lambda$ increases.

BOMD is implemented in CP2K [68], VASP, and Quantum ESPRESSO. The Gaussian and augmented plane waves (GPW) method in CP2K offers favorable scaling for 100–500 molecule systems [69].

3.2. Car-Parrinello Molecular Dynamics (CPMD)

CPMD avoids explicit SCF convergence by treating electronic wavefunctions as dynamical variables with fictitious mass $\mu$ [70]:

$$\begin{array}{cc}\hfill M_I{\ddot{\mathbf{R}}}_I& =-{\nabla }_{{\mathbf{R}}_I}E[\left\{\psi \right\},\left\{\mathbf{R}\right\}]\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \mu {\ddot{\psi }}_i& =-{\displaystyle \frac{\delta E}{\delta {\psi }_i^{*}}}+\sum _j{\Lambda }_{ij}{\psi }_j\hfill \end{array}$$

(3)

where ${\Lambda }_{ij}$ enforces orthonormality.

The method’s validity rests on adiabatic decoupling: the fictitious electronic frequency ${\omega }_e\propto \sqrt{E_{\mathrm{gap}}/\mu }$ must greatly exceed nuclear frequencies ${\omega }_n$, keeping electrons near the ground state. For water, two problems arise at low temperatures. First, the O-H stretching frequency (∼3600 cm⁻¹) sets a high ${\omega }_n$ that constrains the choice of $\mu$, requiring small timesteps of ∼0.12–0.24 fs [71,72]. Second, as nuclear motion slows with decreasing T, maintaining the ratio ${\omega }_e/{\omega }_n$ becomes more delicate—small perturbations can transfer energy from electrons to nuclei, causing drift from the Born-Oppenheimer surface [72].

Advantages: Faster per step than BOMD; no SCF convergence at each step.

Limitations: Requires small timesteps; sensitive to the fictitious electronic mass $\mu$ tuning; potential energy drift for light atoms; shares BOMD’s classical-nuclei limitation.

Standard CPMD with classical nuclei is generally not recommended for water when nuclear quantum effects are required: at room temperature it already shows over-structuring relative to experiment, an effect that is partially corrected by including NQEs via PI-CPMD [73]. The underlying fictitious-mass sensitivity is also exacerbated in water specifically because of the high-frequency O–H stretching mode, which couples to the electronic dynamics and can cause drift from the Born–Oppenheimer surface [72]. For low-temperature water, both limitations remain: classical nuclei become a progressively worse approximation as the de Broglie wavelength of the proton grows — making methods with explicit NQEs (PIMD, PIGLET) the more appropriate choice.

3.3. Path-Integral Molecular Dynamics (PIMD)

Classical treatment of nuclei fails for light atoms at low temperatures. PIMD addresses this by exploiting the isomorphism between quantum statistical mechanics and classical mechanics of ring polymers [74,75]: each quantum particle of mass M is represented by P classical “beads” connected by harmonic springs, with spring constant $k_P=MP{k}_{\mathrm{B}}^2{T}^2/{\hslash }^2$. The temperature dependence of $k_P$ is the central feature: at low T the springs weaken (quadratically), allowing the ring polymer to “open” and directly represent the growing thermal de Broglie wavelength of the quantum particle.

Convergence to the quantum limit requires $P\gtrsim \beta \hslash {\omega }_{\max }$, where ${\omega }_{\max }$ is the highest physical frequency in the system [18]. For water, ${\omega }_{\max }\approx {\omega }_{\mathrm{OH}}\approx 3600$ cm⁻¹.

For water at room temperature, $P\approx 32$ replicas are typically required to converge structural properties [18]; the requirement grows roughly as $1/T$ for fixed ${\omega }_{max}$, making low-temperature PIMD considerably more expensive.

Since standard PIMD evaluates forces independently for each bead, computational cost scales linearly with P [18]. Note that converging O–H structural properties requires at least $P\approx 32$ at 300 K [18]; the $\hslash {\omega }_{max}/k_{B}T$ scaling implies $P\approx 48$–64 at 150–200 K, placing standard ab initio PIMD in the “tens to hundreds of times” more expensive regime relative to classical BOMD [18], which is often prohibitive for production runs.

3.4. Accelerated Path-Integral Methods

Several approaches reduce PIMD’s computational burden while preserving essential quantum effects:

Ring-polymer contraction (RPC) exploits the observation that the physical potential $V\left(\mathbf{R}\right)$ varies more slowly than the inter-bead springs. Expensive DFT forces can be evaluated on a contracted ring ($P^{\prime }<P$ beads), with projection back to the full ring [76]. This is particularly effective when a reference potential can absorb the rapidly varying forces, leaving a smoothly varying remainder for the contracted ring [18].

PI-GLE (Ceriotti, Manolopoulos, Parrinello) [77] employs a generalized Langevin equation thermostat with colored noise to accelerate the convergence of path integral molecular dynamics in anharmonic systems exhibiting both zero-point energy and tunnelling effects. Unlike white noise (uncorrelated in time), colored noise mimics the quantum fluctuation spectrum, designed so that a P-replica PIMD simulation converges structural observables in the harmonic limit for any P. PIGLET [78] extends this scheme with additional constraints on the noise spectrum, so that the quantum kinetic energy estimator also converges in the harmonic limit at small P. Together, these approaches achieve converged NQEs with as few as ∼6 beads at 300 K for liquid water [18].

The key insight underlying both methods is that quantum statistics require energy distributed according to ${\displaystyle \frac{1}{2}}\hslash \omega coth(\beta \hslash \omega /2)$ per mode rather than the classical equipartition value $k_{\mathrm{B}}T$; the GLE thermostat enforces this distribution while reducing the number of beads required for convergence. Both methods are implemented in i-PI [79], which interfaces with CP2K and other electronic-structure codes.

3.5. Machine-Learning Interatomic Potentials

Machine-learning interatomic potentials (MLPs), including neural-network potentials and kernel-based methods such as Gaussian Approximation Potentials, are trained on a database of ab initio (typically DFT) energies and forces, and then used as a drop-in replacement for the underlying potential energy surface during MD [80,81]. Once trained, force evaluations cost a small fraction of a DFT step. They can be combined with BOMD, PIMD, PIGLET, RPMD, or CMD, and the speedup applies in each case. In the path-integral context in particular, the P-fold overhead of bead sampling becomes affordable, and nanosecond-scale path-integral trajectories of condensed-phase water at near-DFT accuracy are now feasible [82].

A comprehensive review of MLPs for water is beyond the scope of this work; we refer the reader to recent dedicated reviews on neural-network potentials [83] and on data-driven potentials for water specifically [84]. We note one caveat directly relevant to the present review: the great majority of MLPs for water reported to date have been trained and validated on reference data sampled at or near ambient conditions. Their transferability to the deeply supercooled regime (200–250 K), and especially to the amorphous ice phases discussed in Section 2.2, is not guaranteed and remains an active research question. Studies that use MLPs for low-temperature water should therefore include explicit validation against ab initio reference points at the temperatures and densities of interest, rather than relying on extrapolation from ambient-condition training sets.

3.6. Comparison of AIMD Methods

Table 3 provides a comprehensive comparison of AIMD methods for water simulations.

General remarks for method selection:

Method selection depends on the target observable, the required accuracy, and the temperature regime, rather than on a single criterion. Table 3 summarizes the trade-offs. Broadly: classical-nuclei AIMD (BOMD, or CPMD with a small fictitious mass) is appropriate when NQE corrections to the target property are smaller than other modeling errors; path-integral methods (PIMD, PIGLET) are needed when quantum statistics of the nuclei matter, with the number of beads P growing as T decreases; and approximate quantum-dynamics methods (RPMD, CMD/QCMD, instanton) are needed when real-time correlation functions, vibrational spectra, or rate constants are required. The boundaries between these regimes are not sharp and should be assessed against benchmarks for the specific property of interest.

4. Density Functional Selection for Water

The functional choice is perhaps the most consequential decision in DFT-based AIMD of water. Water presents a stringent test for density functionals because accurate simulation requires simultaneous description of:

• Intramolecular covalent bonding (O-H stretching, H-O-H bending)
• Intermolecular hydrogen bonding (electrostatic + charge transfer)
• London dispersion interactions
• Many-body polarization effects

Different classes of functionals handle these interactions with varying fidelity, and errors in any component propagate into structural and dynamical predictions [19]. Below we discuss these nuances family by family, presenting representative members of each class, including both widely-used functionals and those that have performed comparatively well in benchmark studies of liquid water and ice. We note an important caveat: most systematic benchmarks of these functionals have been carried out at or near ambient conditions, including the broad surveys of Gillan et al. [19] and the more recent force-based benchmark of 50+ functionals by Cui [90]. Direct DFT-AIMD validation of any functional in the supercooled regime (200–250 K) and for the amorphous ice phases of Section 2.2 remains sparse. Most of what is known about low-temperature water at the DFT level comes indirectly through machine-learning potentials trained on DFT reference data [82,91,92,93], which then enable the much longer trajectories needed to probe supercooled and amorphous states.

4.1. Generalized Gradient Approximation (GGA) Functionals

GGA functionals are computationally efficient but suffer from well-documented limitations for water:

BLYP [94,95]: Without dispersion correction, BLYP over-structures liquid water and underestimates its diffusivity (i.e., the dynamics becomes too sluggish) [96,97]; the density of the liquid is also significantly under-predicted [97]. Adding an empirical Grimme dispersion correction (BLYP-D in the original Schmidt et al. study, equivalent to modern BLYP-D3 in current practice) significantly improves the predicted density and brings the O–O radial distribution function into close agreement with experiment [97]. The resulting BLYP-D3 potential energy surface is comparatively “soft,” and several authors have argued that this softness partially compensates for the missing nuclear quantum delocalization in classical-nuclei simulations [19]. Performance of BLYP-D3 in the supercooled regime (200–250 K) is sparse. However some structural features (ie. g(r)O–O) are covered in the 225–350 K range in Tatini2015.

Advantages: Fast computation; soft PES partially compensates for missing NQE; good for equilibration.

Disadvantages: Bare BLYP underbinds the water dimer (4.18 kcal mol⁻¹ vs 5.44 kcal mol⁻¹ experimental) [98] and produces an overstructured liquid with too slow diffusion ($\sim 4\times$ below experiment) [98,99]. Adding D3 brings density and the O–O RDF close to experiment [97], but barrier heights remain unreliable: the GMTKN55 BH76 WTMAD-1 is 5.19 kcal mol⁻¹ for BLYP and 5.53 kcal mol⁻¹ for BLYP-D3(BJ), reflecting an underlying self-interaction error that dispersion correction does not address [100].

PBE [101]: PBE tends to over-structure radial distribution functions of water even at ambient conditions, predicting too many hydrogen bonds and underestimating diffusivity, with respect to experiments [98]. PBE beneffits from adding dispersion (PBE-D3) improving structural parameters. However, dispersion correction has limited effect on bulk water structure because PBE already overbinds molecules [102].

Advantages: Well-validated; widely available; reasonable lattice constants for ice I_h [98].

Disadvantages: Over-structures liquid water, an artifact not fully removed by adding dispersion corrections [96,98]; poor dynamics in the supercooled regime[103]; severely overestimates the melting point of ice [104], equivalent to over-stabilizing ice relative to liquid.

revPBE-D3 [105,106]: revPBE without dispersion correction performs poorly for condensed-phase water, but combined with Grimme’s D3 dispersion correction it has been reported to describe ambient liquid water (298 K, 1 atm) remarkably well in classical AIMD, with equilibrium density, $g_{\mathrm{OO}}\left(r\right)$, self-diffusivity, IR spectrum, and dipole moment in reasonable agreement with experiment [107]. This performance partly reflects a cancellation between intrinsic functional error and the neglect of nuclear quantum effects in classical-nuclei trajectories [107]. Direct validation in the supercooled regime is sparse: in the ab initio benchmark of Herrero et al. down to ∼248 K, the SCAN functional was found to perform better, while revPBE-D3 was identified as a plausible alternative that describe the temperature evolution of the shear viscosity and self-diffusion coefficient of liquid waterwith good performance, but was not itself tested across the full supercooled range [103]. Consequently, applicability of revPBE-D3 to the 200–273 K regime should be regarded as inferred from ambient-temperature benchmarks rather than directly validated.

Table 4. Summary of GGA functional performance for water at ambient conditions (∼300 K, 1 atm). Density and $r_{\mathrm{OO}}^{\left(1\right)}$ errors are deviations from experimental liquid-water values, computed as (functional value − experimental value) using raw simulation data tabulated in the cited primary references. The experimental reference values $\rho =1.00$ g cm⁻³ (room temperature, 1 atm) and $r_{\mathrm{OO}}^{\left(1\right)}=2.80$Å [108] are also taken from the same primary references. The barrier error is the mean absolute deviation from CCSD(T)/CBS reference values on the gas-phase BH76 subset of GMTKN55 [100], taken directly from the supporting information of that work, and is reported as a generic indicator of transition-state energetics, not a water-specific benchmark. revPBE values are from simulations performed at 330 K (rather than 300 K), as reported in Ref. [109] to compensate for revPBE’s over-structuring tendency, and is therefore not directly comparable to the other rows. No directly comparable benchmark data was found for any of these functionals in the supercooled regime (200–273 K).

Functional	Density Error	${\mathit{r}}_{\mathbf{OO}}^{\left(1\right)}$ Error	Barrier Error (BH76)
	(%)	(Å)	(kcal mol−1)
BLYP	$-20.3$ [110]	$+0.03$ [110]	$5.19$ [100]
BLYP-D3	$+6.6$ [110]	$-0.02$ [110]	$5.53$ [100]
PBE	$-12$ [97]	$-0.03$ [98]	$9.15$ [100]
PBE-D3	$+5.5$ [110]	$-0.07$ [110]	$9.62$ [100]
revPBEa	$-31$ [109]	$+0.08$ [109]	$7.11$ [100]
revPBE-D3 [105,106]	$-6$ [111]	$+0.01$ [111]	$8.32$ [100]

^a

Table 4. Summary of GGA functional performance for water at ambient conditions (∼300 K, 1 atm). Density and $r_{\mathrm{OO}}^{\left(1\right)}$ errors are deviations from experimental liquid-water values, computed as (functional value − experimental value) using raw simulation data tabulated in the cited primary references. The experimental reference values $\rho =1.00$ g cm⁻³ (room temperature, 1 atm) and $r_{\mathrm{OO}}^{\left(1\right)}=2.80$Å [108] are also taken from the same primary references. The barrier error is the mean absolute deviation from CCSD(T)/CBS reference values on the gas-phase BH76 subset of GMTKN55 [100], taken directly from the supporting information of that work, and is reported as a generic indicator of transition-state energetics, not a water-specific benchmark. revPBE values are from simulations performed at 330 K (rather than 300 K), as reported in Ref. [109] to compensate for revPBE’s over-structuring tendency, and is therefore not directly comparable to the other rows. No directly comparable benchmark data was found for any of these functionals in the supercooled regime (200–273 K).

Functional	Density Error	${\mathit{r}}_{\mathbf{OO}}^{\left(1\right)}$ Error	Barrier Error (BH76)
	(%)	(Å)	(kcal mol−1)
BLYP	$-20.3$ [110]	$+0.03$ [110]	$5.19$ [100]
BLYP-D3	$+6.6$ [110]	$-0.02$ [110]	$5.53$ [100]
PBE	$-12$ [97]	$-0.03$ [98]	$9.15$ [100]
PBE-D3	$+5.5$ [110]	$-0.07$ [110]	$9.62$ [100]
revPBEa	$-31$ [109]	$+0.08$ [109]	$7.11$ [100]
revPBE-D3 [105,106]	$-6$ [111]	$+0.01$ [111]	$8.32$ [100]

^a

4.2. Meta-GGA Functionals

Meta-GGA functionals incorporate the kinetic energy density $\tau$, providing additional flexibility:

SCAN [112]: The strongly constrained and appropriately normed (SCAN) functional satisfies all 17 known exact constraints for a semi-local functional. For water, SCAN captures the density difference between liquid water and ice I_h at ambient conditions and reproduces many structural, electronic, and dynamic properties of liquid water without empirical dispersion corrections [113], and correctly predicts ice I_h as more stable than ice II_c along with reasonable melting properties of both polymorphs [114]. However, SCAN overestimates the liquid water density by approximately 5% relative to experiment [113] and can exhibit numerical grid sensitivity in some implementations [115,116].

Advantages: No empirical parameters; good ice energetics including correct I_h/II_c stability ordering; captures intermediate-range dispersion.

Disadvantages: Over-dense liquid (∼5%); grid sensitivity; computational overhead (1.5–2× GGA); can be unstable in some codes.

r²SCAN [116]: A regularized version of SCAN with improved numerical stability.

For water, r²SCAN-driven AIMD reproduces the structure of bulk liquid water at ambient density when simulations are run at an elevated temperature of 350 K, similar to the compensation required for SCAN, which exhibits excessive structuring at 300 K and an augmented first solvation shell at 330 K [117]. The first O–O peak position of r²SCAN (2.865 Å) sits slightly above the experimental value (2.801 Å), indicating mild residual over-structuring [117]. Furness et al. report that r²SCAN converges with fewer iterations than SCAN [116].

Advantages: Improved numerical convergence over SCAN [116]; comparable structural accuracy for bulk liquid water at ambient conditions [117].

Disadvantages: Inherits SCAN’s need for an elevated simulation temperature (∼350 K) to match the 300 K experimental RDF [117]; mild over-structuring of the first O–O peak relative to experiment [117].

DC-SCAN (Density-corrected SCAN): Applies a Hartree–Fock density correction scheme to SCAN, minimizing density-driven errors and elevating the accuracy of the SCAN functional to that of coupled-cluster theory for water clusters and interaction energies [118].

Advantages: Reaches near-coupled-cluster accuracy for water cluster energetics; reduces SCAN’s overbinding/over-structuring tendencies; provides a reliable reference for parameterizing many-body and machine-learned potentials.

Disadvantages: Not standard in most codes; requires a Hartree–Fock density evaluation in addition to the SCAN energy [118]; additional computational overhead; condensed-phase MD applications rely on derived many-body potentials (e.g., MB-SCAN(DC)) rather than direct AIMD [118].

4.3. Hybrid Functionals

Hybrid functionals include a fraction of exact (Hartree-Fock) exchange, improving descriptions of charge transfer and reducing self-interaction error [119].

In the context of liquid water, hybrid functionals combined with nuclear quantum effects (NQEs) have been shown to provide a markedly improved description of hydrogen bonding, with revPBE0-D3 + PIMD identified as one of the few approaches that closely agrees with experiment [120].

PBE0 [121]: With 25% exact exchange and, when paired with D3 dispersion correction (PBE0-D3), this functional provides excellent barrier heights for reactions [122,123] and reliable magnetic properties [124]. Under classical nuclei, PBE0-D3 tends to over-structure liquid water [125].

The computational cost is substantially higher than GGA due to the evaluation of exact exchange (commonly an order of magnitude or more in conventional Gaussian or plane-wave AIMD implementations), although schemes like ADMM can reduce this overhead by roughly a factor of 20 [126].

Advantages: Excellent barrier heights (MAD ∼1–2 kcal mol⁻¹ on representative benchmarks) [123]; reliable magnetic properties [124]; reduced self-interaction error [119]; improved hydrogen-bond description, particularly when combined with NQEs [127].

Disadvantages: Substantial cost increase over GGAs [126]; tends to over-structure water under classical nuclei [125]; requires dispersion correction (e.g., D3 [106]) for condensed phases [125].

revPBE0: Combines the exchange of revPBE with 25% exact exchange. This functional has been shown to accurately reproduce the structure, dynamics, and vibrational spectroscopy of bulk liquid water when combined with path-integral simulations and D3 dispersion. [82,120,128], and is one of the few approaches that yields quantitative agreement with experiment for the O–O radial distribution function, self-diffusion coefficient, and full IR and Raman spectra at ambient conditions [120,128].

Advantages: Among the best-performing functionals for bulk liquid water when combined with NQEs, capturing the O–O RDF, diffusion, and IR/Raman spectra in close agreement with experiment [82,120,128]; balanced description of covalent bonding and hydrogen bonding such that competing nuclear quantum effects on the hydrogen bond cancel correctly [120].

Disadvantages: Substantial cost overhead from the evaluation of exact exchange in periodic AIMD, although acceleration schemes such as the auxiliary density matrix method (ADMM) significantly reduce this cost [126]; NQE treatment essential for quantitative agreement of dynamics and vibrational spectra [120,128].

$\omega$B97X-V [129]: A range-separated hybrid GGA with non-local correlation functional VV10, designed against a broad training set covering thermochemistry, kinetics, and non-covalent interactions. Its meta-GGA counterpart $\omega$B97M-V [130] extends this design to the meta-GGA rung. On benchmark sets of small water clusters with coupled-cluster reference binding energies, $\omega$B97X-V and $\omega$B97M-V have been reported to give the smallest errors among the DFT functionals tested in that set [131]. A closely related meta-GGA, B97M-rV (B97M-V with the rVV10 non-local correlation [132]), gives bulk liquid water structure, diffusion, and IR spectroscopy of comparable accuracy to the much more expensive hybrid revPBE0-D3 when combined with path-integral simulations [128], and is therefore an attractive meta-GGA-level alternative to hybrid functional for bulk water.

Advantages: Smallest binding-energy errors among DFT functionals tested on water-cluster benchmarks against coupled-cluster references [131]; built-in VV10 non-local correlation, so no separate dispersion correction is required [129,133]; the related meta-GGA B97M-rV reaches near-hybrid accuracy for bulk water at meta-GGA cost [128].

Disadvantages: Substantial cost relative to GGA functionals, from range-separated exact exchange combined with VV10 non-local correlation evaluation; less commonly used in condensed-phase AIMD of liquid water than global hybrids such as PBE0 and revPBE0.

4.4. Dispersion Corrections

Standard DFT functionals miss long-range dispersion interactions, which contribute significantly to ice cohesion (see [134,135] ). Common corrections include:

• D3(BJ) [106,136]: Damped pairwise $C_6/R^6$ correction with Becke-Johnson damping. Well-validated, widely available.
• D4 [137] Improved atom-pairwise dispersion with charge-dependent coefficients.
• VV10 [133]: Vydrov and Van Voorhis non-local correlation functional that captures dispersion from the electron density alone. More physically motivated than atom-pairwise schemes, but computationally more demanding due to the six-dimensional double integral over the density [138].

For water simulations, D3(BJ) is chosen as the default method due to its extensive validation, low computational overhead, and availability in all major codes.

5. Nuclear Quantum Effects: The Critical Challenge Below 250 K

5.1. Physical Origin of Nuclear Quantum Effects

Because electrons move much faster than nuclei due to their significantly smaller mass, the Schrödinger equation is commonly simplified using the Born-Oppenheimer approximation. This approximation underlies the most widely used methods for molecular dynamics, including BOMD, CPMD, and standard DFT-based approaches. According to Born-Oppenheimer, the electronic wavefunction is solved quantum mechanically at fixed nuclear positions, separating the electronic and nuclear degrees of freedom and the nuclei movement. This generates a potential energy surface on which nuclei evolve. A further approximation (treating nuclei as classical point particles obeying Newton’s equations) is applied in standard MD methods. This classical treatment of nuclei provides accurate results with reasonable computational cost under most conditions but inherently neglects nuclear quantum effects (NQE) such as zero-point energy (the minimum energy a particle possesses even at 0 K) and quantum tunneling, where a nucleus can traverse a potential energy barrier even when its kinetic energy is classically insufficient. For most systems at ambient conditions, neglecting NQE introduces minimal error and is therefore an acceptable approximation [66].

The fundamental limitation for low-temperature water lies in treating nuclei as classical point particles obeying Newton’s equations. This approximation neglects zero-point energy, quantization of vibrational energy levels, and tunneling. The importance of these nuclear quantum effects can be quantified by the ratio of zero-point energy to thermal energy, $\hslash \omega /2k_{\mathrm{B}}T$. For the O–H stretching mode in water ($\omega \approx 3600$ cm⁻¹ [67]), the zero-point energy $E_0={\displaystyle \frac{1}{2}}\hslash \omega \approx 5.1$ kcal mol⁻¹, which already exceeds the thermal energy $k_{\mathrm{B}}T$ by 9 times at room temperature [18]. As temperature decreases, this ratio grows further: at 200 K, $k_{\mathrm{B}}T\approx 0.4$ kcal mol⁻¹, making the ZPE approximately 13 times larger than the thermal energy available to the O–H coordinate.

Consequently, Eq. (1) introduces systematic errors that grow progressively with decreasing temperature, as the classical treatment cannot account for the proton delocalization, enhanced hydrogen-bond fluctuations, and tunneling contributions that become increasingly dominant in the supercooled regime [18,127].

The three principal phenomena of NQE are:

1. Zero-point energy (ZPE): Quantum mechanics mandates that a particle confined in a potential well possesses a minimum energy $E_0={\displaystyle \frac{1}{2}}\hslash \omega$ even at absolute zero temperature, where $\omega$ is the vibrational frequency. For the O-H stretching mode in water, the estimated ZPE 5.1 kcal mol⁻¹ is comparable to the hydrogen bond strength itself (∼5 kcal mol⁻¹). This ZPE causes hydrogen atoms to be more delocalized than classical mechanics predicts at all temperatures, effectively “swelling” the proton distribution perpendicular to the bond direction. This delocalization persists at ambient conditions but becomes increasingly dominant as temperature decreases and the ratio ZPE/$k_{\mathrm{B}}T$ grows.

2. Quantum tunneling: Nuclei can traverse potential energy barriers even when their kinetic energy is classically insufficient to surmount them. For light nuclei such as hydrogen, this effect is particularly pronounced: the tunneling probability depends exponentially on both the barrier width and the particle mass, and the low mass of the proton makes hydrogen transfer reactions especially susceptible.

3. Exchange effects: In principle, identical nuclei should obey quantum exchange statistics (Bose-Einstein or Fermi-Dirac). However, in molecular water and ice, protons are localized in deep covalent potential wells (∼1 eV), which exponentially suppresses the inter-molecular wavefunction overlap required for exchange. All path integral simulations of water treat nuclei as distinguishable particles [18,139], and this approximation is considered valid at all temperatures relevant to condensed-phase water, in contrast to weakly bound systems such as superfluid helium where exchange effects dominate [140]. Exchange effects are therefore not considered further in this review.

Competing quantum effects in hydrogen bonding. Importantly, NQE do not simply weaken or strengthen hydrogen bonds—they produce competing effects [127,141]:

• Intramolecular ZPE (O-H stretch): Causes proton delocalization perpendicular to the bond, effectively weakening the hydrogen bond by increasing the average O…H distance.
• Intermolecular ZPE (O…O stretch): Causes delocalization along the hydrogen bond axis, which can strengthen hydrogen bonds by allowing closer O-O approach in the delocalized wavefunction.

The balance between these competing effects is functional-dependent and temperature-dependent, contributing to the challenge of accurate low-temperature simulations.

5.2. Magnitude of NQE Versus Temperature

The ratio of zero-point energy (ZPE) to thermal energy provides a quantitative measure of NQE importance:

Table 5. Comparison of thermal energy and H-atom zero-point energy in water. The ZPE is computed from the O-H stretching frequency (∼3600 cm⁻¹, corresponding to ∼10.3 kcal mol⁻¹ for two O-H bonds, or ∼5.1 kcal mol⁻¹ per bond). The “NQE Regime” classification follows Ref. [18].

T (K)	${\mathit{k}}_{\mathbf{B}}\mathit{T}$ (kcal mol−1)	ZPE(OH)/(${\mathit{k}}_{\mathbf{B}}\mathit{T}$)	${\mathit{\Lambda }}_{\mathbf{H}}$ (Å)	NQE Regime
300	0.60	1.5–2	1.0	Moderate
250	0.50	2–2.5	1.1	Significant
200	0.40	2.5–3	1.2	Strong
150	0.30	3–4	1.4	Very strong
77	0.15	6–8	2.0	Extreme

Table 5. Comparison of thermal energy and H-atom zero-point energy in water. The ZPE is computed from the O-H stretching frequency (∼3600 cm⁻¹, corresponding to ∼10.3 kcal mol⁻¹ for two O-H bonds, or ∼5.1 kcal mol⁻¹ per bond). The “NQE Regime” classification follows Ref. [18].

T (K)	${\mathit{k}}_{\mathbf{B}}\mathit{T}$ (kcal mol−1)	ZPE(OH)/(${\mathit{k}}_{\mathbf{B}}\mathit{T}$)	${\mathit{\Lambda }}_{\mathbf{H}}$ (Å)	NQE Regime
300	0.60	1.5–2	1.0	Moderate
250	0.50	2–2.5	1.1	Significant
200	0.40	2.5–3	1.2	Strong
150	0.30	3–4	1.4	Very strong
77	0.15	6–8	2.0	Extreme

When ZPE/$k_{\mathrm{B}}T>2$, classical treatment of nuclei introduces qualitative errors in properties sensitive to hydrogen motion. The thermal de Broglie wavelength ${\Lambda }_{\mathrm{H}}$ exceeding the O-H bond length signals that protons cannot be treated as localized classical particles.

5.3. Structural Consequences of NQE

Classical nuclei predictions diverge from quantum reality in several ways at low temperature [127,142]:

Density trends: Classical MD predicts monotonically increasing density upon cooling, whereas PIMD reveals density maxima in both liquid water and LDA [142]. This qualitative difference arises because ZPE-induced volume expansion counteracts thermal contraction at low T.

Hydrogen bond geometry: Classical simulations produce overly linear, overly short O-H…O hydrogen bonds. NQE-induced proton delocalization (perpendicular to the O-H bond direction) results in larger O-O separations and more bent hydrogen bonds [73,143].

Radial distribution functions: First-neighbor peaks in $g_{\mathrm{OO}}\left(r\right)$ and $g_{\mathrm{OH}}\left(r\right)$ are sharper in classical simulations than with NQE. Including quantum effects “smears” these peaks toward experimental data.

Vibrational spectra: The O-H stretching frequency is red-shifted in condensed phases due to hydrogen bonding. NQE further modulate these shifts through anharmonic coupling and Fermi resonances.

6. Conclusions

Simulating water in the 200–273 K range requires confronting three coupled challenges: degraded sampling as dynamics slow, the growing magnitude of nuclear quantum effects (NQE), and the residual sensitivity of structural and dynamical observables to the underlying potential energy surface. The picture that emerges from the methods and benchmarks surveyed here is that none of these issues can be treated in isolation, and that the most reliable strategies are those that match the rigor of the nuclear treatment to the rigor of the electronic structure.

At the methodological level, the hierarchy is clear. BOMD on a converged Born–Oppenheimer surface remains the baseline when NQE corrections are smaller than other modeling errors [18,66]. CPMD [70] is competitive per step but its fictitious-mass dependence and the high O–H stretching frequency make it a poor choice when low-temperature accuracy is required [72,73]. Path-integral methods (PIMD [74,75], PIGLET [77,78], RPC [76]) are the appropriate tool when quantum statistics of the nuclei matter, with the bead count P scaling with $\beta \hslash {\omega }_{\max }$. Real-time approximations (RPMD [85,86], CMD/QCMD [88,89]) remain the route to vibrational spectra and rate constants, but each carries known artefacts that are amplified at low T [87].

At the electronic-structure level, density functional selection cannot be decoupled from the nuclear treatment. GGAs with dispersion correction (BLYP-D3, revPBE-D3) benefit from a known cancellation between intrinsic functional softness and the neglect of NQE in classical-nuclei trajectories [107]; their good performance at ambient temperature is not, in itself, evidence of accuracy at 200 K. Meta-GGAs (SCAN, r²SCAN) capture intermediate-range dispersion and ice polymorph energetics well [112,113,114,116], but require elevated simulation temperatures to match ambient experimental structure [117]. Hybrid functionals combined with NQE (revPBE0-D3 + PIMD) are among the few approaches that quantitatively reproduce bulk water structure and spectroscopy [120,127,128], at a computational cost that has historically restricted them to short trajectories on small cells.

Below approximately 250 K, NQE shift from a correction to an organizing principle. The ratio of zero-point energy to thermal energy for the O–H stretch grows past order unity, and properties as basic as the temperature of maximum density of low-density amorphous ice cannot be recovered without quantum nuclei [18,142]. Within this regime, choices that look conservative at ambient conditions—classical nuclei, GGA functionals selected on room-temperature benchmarks—become uncontrolled. The literature consistently treats protons as distinguishable particles [139], and the deep covalent localization of the O–H bond makes this approximation robust for condensed water, unlike systems where exchange dominates [140].

The most consequential recent development is the maturation of machine-learning interatomic potentials trained on high-level ab initio reference data [80,81,83,84]. By reducing the per-step cost of force evaluation by orders of magnitude, MLPs make nanosecond-scale path-integral trajectories with near-coupled-cluster accuracy feasible on system sizes appropriate for nucleation, phase-equilibrium, and amorphous-ice studies [82,91,92,93]. This is the route by which most quantitative information about low-temperature water at the DFT level has reached the literature in recent years. The principal caveat is that MLPs trained on ambient-condition data are not guaranteed to transfer to 200–250 K or to amorphous phases; explicit ab initio validation at the target conditions remains essential.

Several methodological gaps remain. Direct DFT-AIMD benchmarks in the supercooled regime are sparse, and most existing functional rankings for water rest on ambient-temperature data. Thermodynamic reweighting between functional ensembles is statistically unreliable at low temperature, so functional comparisons require independent simulations. Approximate quantum-dynamics methods need further development for vibrational spectroscopy below 200 K. We expect the next generation of MLP-driven, path-integral simulations, trained on hybrid or density-corrected reference data and validated against the experimental structural targets summarised in Section 2.2, to constitute the practical workflow for the next decade of low-temperature water simulation.

Beyond methodological refinement, the value of these advances lies in the questions they make accessible. Supercooled and amorphous water are not abstract benchmarks: they are the relevant phases for ice nucleation in mixed-phase clouds and the freezing of micrometric atmospheric droplets [64,144], for the survival of vitrified organs and tissues under sub-zero preservation protocols [7,8], for the control of nucleation in cryopreservation and food processing [10,11], and for the chemistry that occurs on amorphous solid water grain mantles in dense interstellar clouds, where it provides the substrate for the formation of complex organic molecules and prebiotic precursors [12,13,14,15]. Amorphous and crystalline ices are also dominant surface phases on icy moons, comets, and Mars [16]. In each of these settings, predictive simulation requires not just qualitative agreement with experiment but quantitative reproduction of structure, dynamics, and phase behaviour at temperatures where NQE are non-negligible and where direct experimental access is fragmented. The methodological convergence outlined above—path-integral sampling driven by machine-learning potentials trained on hybrid-functional reference data—is what brings these problems within reach of routine in silico investigation, at a cost that finally permits the system sizes and trajectory lengths that the underlying physics demands.