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]. Amorphous solid water (ASW), formed by vapor deposition at cryogenic temperatures, is a widely form of water found in the universe, for instance, coating interstellar dust grains where it serves as a substrate for complex organic molecule (COM) formation [11,12,13,14]. In planetary science, water ice in various amorphous and crystalline forms is found on icy moons, comets, and Mars [15]. 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 [16].

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 magnitude, which increases when the temperature drops, [17] indicating that classical treatment of nuclei introduces significant systematic errors. 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 [18].

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 [19], and to even lower temperatures using rapidly cooled nanodroplets [20] and ultrafast X-ray probing of evaporatively cooled microdroplets [21], 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 [22]. 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 [23]. At ambient conditions, water already exhibits nanometer-scale density fluctuations representing tetrahedral-like and hydrogen-bond-distorted structures [24]. 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 [25]. 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 [21]. Solid water has been obtained at 77 K in different distinct glassy forms, from low-density amorphous ice (LDA) to high-density amorphous ice (HDA) [26,27]. However, our understanding of water’s structure below the homogeneous nucleation temperature, $T_{\mathrm{H}}\approx 235$ K, whether it reflects continuous structural evolution or a discrete liquid–liquid transition, continues to evolve, and computational simulations provide much of the structural insight into 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.

The low-density amorphous ices can be prepared by three routes that yield structurally similar forms of the same density (∼0.94 g cm⁻³): (i) vapour deposition of water onto a cold substrate, giving amorphous solid water (ASW) [28]; (ii) decompression (or controlled heating) of HDA, giving low-density amorphous ice (LDA), with the reversible HDA→LDA transition occurring at ∼135 K and ∼0.2 GPa [29]; and (iii) hyperquenching of micrometre-sized liquid droplets at cooling rates of order ${10}^7$ K/s, which vitrifies water to hyperquenched glassy water (structurally akin to LDA) [30,31].

In the interstellar medium, water is the most abundant ice species in grain mantles, built up by freeze-out of gas-phase species onto cold dust grains; this ice is predominantly amorphous [32].

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$ Å [33,34]. LDA exhibits a locally tetrahedral hydrogen-bond network resembling ice ${\mathrm{I}}_{\mathrm{h}}$ but lacking long-range crystallographic order [33].

Key remarks. The glass transition temperature is $T_g\approx 136$ K, as determined calorimetrically [35]. LDA is metastable and crystallizes to stacking-disordered ice ${\mathrm{I}}_{\mathrm{sd}}$ upon heating at ambient pressure, at about ∼150 K [36].

LDA has a density of approximately 0.94 g cm⁻³ at 77–80 K and ambient pressure. This value is consistent with reported atomic number densities of ∼0.094 atoms/ų for LDA [33,34]; 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⁻³) [26,38].

Formation temperature and conditions. MDA is produced by ball milling hexagonal ice (ice ${\mathrm{I}}_{\mathrm{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 [39]. 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 [40].

Structure. The oxygen–oxygen pair-distribution function of MDA reveals a progressive collapse of the second coordination shell along the sequence ice ${\mathrm{I}}_{\mathrm{h}}$ → LDA → MDA → liquid water, placing MDA structurally closer to liquid water than to LDA [37]. Importantly, MDA is not a simple mixture of LDA and HDA; comparison of its diffraction pattern with other ice forms identifies it as a structurally unique material rather than a transition state or a mixture of existing forms [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]. 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 [27]; the value reported for HDA obtained by compressing LDA at 77 K was ∼1.19 g cm⁻³ at zero pressure [27], while the originally reported HDA densities span ∼1.13–1.17 g cm⁻³ depending on the specific preparation protocol and degree of structural relaxation, with relaxed (expanded) HDA being ∼0.02 g cm⁻³ less dense than unannealed HDA [41,42].

Formation temperature and conditions. HDA is formed by: (i) pressure-induced amorphization of ice ${\mathrm{I}}_{\mathrm{h}}$ at 77 K under ∼1 GPa [26]; or (ii) compression of LDA at temperatures below ∼150 K, with the transition occurring at ∼0.6 GPa [27,43]. These transitions occur below the crystallization temperature of amorphous ice (∼150 K at ambient pressure) [44]. 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 [44,45].

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 [34]. The structure features interpenetrating hydrogen-bond networks in which the second coordination shell collapses inward compared to the open tetrahedral arrangement of LDA [34]. Neutron diffraction studies of the five amorphous ices at 80 K and ambient pressure show that HDA is structurally distinct from both LDA and VHDA in its intermediate-range order [33].

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

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 [38], approximately 9% higher than that of HDA. Values up to 1.26 g cm⁻³ have been reported depending on preparation pressure [48].

Formation temperature and conditions. VHDA is prepared by isobaric annealing of HDA at pressures of 1.1–1.9 GPa and temperatures of 165–177 K, followed by quenching to 77 K [38,48,49]. VHDA forms at pressures >0.8 GPa [48]. A stepwise transformation sequence LDA → HDA → VHDA (with HDA → VHDA at ∼0.95 GPa) has been observed at 125 K [43].

Structure. The first O–O coordination shell remains near ∼2.8 Å [38,50]. 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 [38]. Under pressure (∼2 GPa), the second coordination shell moves inward as interstitial neighbors are incorporated, distinguishing VHDA from HDA [51].

Key remarks. VHDA prepared at high pressures and temperatures (1.9 GPa, 175 K) represents the amorphous state of water at >0.7 GPa least affected by nanocrystalline remnants [49]. The activation barrier for the VHDA → HDA back-conversion is higher than for HDA → LDA [52]. 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 [49].

2.2.5. Hyperquenched Glassy Water (HGW)

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

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 on the order of ${10}^6$–${10}^7$ K/s [31].

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

Key remarks. The glass transition temperature of HGW is $T_g\approx 136$ K according to Johari, Hallbrucker, and Mayer [35,54]. However, this value is debated; Velikov, Borick, and Angell argued based on enthalpy relaxation scaling that $T_g=165\pm 5$ K [55], while Kohl et al. presented DSC evidence supporting the 136 K assignment [56]. 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 [56].

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 ${\mathrm{I}}_{\mathrm{h}}$ (hexagonal ice) is the thermodynamically stable form of ice I at ambient pressure, whereas ice ${\mathrm{I}}_{\mathrm{c}}$ is metastable cubic ice [63]. In practice, what has often been reported as “cubic ice” is typically stacking-disordered ice ${\mathrm{I}}_{\mathrm{sd}}$, containing cubic sequences interlaced with hexagonal sequences [63], 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 [69].

Ice XI is the proton-ordered form of ice ${\mathrm{I}}_{\mathrm{h}}$ and is thermodynamically stable below ∼72 K (near ambient pressure) [70].

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 [71]. 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⁻¹ [72]), 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 [17]. 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. With sufficiently tight self-consistent-field convergence at each step, it conserves the physical total energy, whereas CPMD conserves a different (extended-Lagrangian) constant of motion [71].

However, because the nuclei are propagated classically, BOMD is unable to compute zero-point energy or nuclear tunneling; errors grow as T decreases and the de Broglie thermal wavelength $\Lambda$ increases.

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

3.2. Car-Parrinello Molecular Dynamics (CPMD)

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

$$\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^{\ast }}}+\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. 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 [76,77]. As nuclear motion slows with decreasing T, maintaining the ratio ${\omega }_e/{\omega }_n$ becomes harder: small perturbations can transfer energy from electrons to nuclei, causing drift from the Born-Oppenheimer surface [77].

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. Classical first-principles water is over-structured relative to experiment at room temperature, and including nuclear quantum effects softens the radial distribution functions and systematically improves agreement with experiment [78]. 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 [77]. 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 [79,80]: each quantum particle of mass M is represented by P classical “beads” connected by harmonic springs, with spring constant $k_P=M{\omega }_P^2=M{P}^2{k}_{\mathrm{B}}^2{T}^2/{\hslash }^2$, where ${\omega }_P=P{k}_{\mathrm{B}}T/\hslash$. 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 [17]. 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 [17]; 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 [17]. Note that converging O–H structural properties requires at least $P\approx 32$ at 300 K [17]; the $\hslash {\omega }_{max}/k_{B}T$ scaling implies $P\approx 48$–64 at 150–200 K, placing standard ab initio PIMD tens to hundreds of times more expensive than classical BOMD, 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 [81]. This is particularly effective when a reference potential can absorb the rapidly varying forces, leaving a smoothly varying remainder for the contracted ring [17].

PI-GLE (Ceriotti, Manolopoulos, Parrinello) [82] 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, 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 [83] 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 [83].

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 [84], 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 [85,86]. 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 [87].

A comprehensive review of MLPs for water is beyond the scope of this work; we refer the reader to a recent dedicated review on neural-network potentials [88], and note that data-driven quantum simulations have recently reproduced the water phase diagram at near-first-principles accuracy [89]. 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 [18]. 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. [18] and the more recent force-based benchmark of 50+ functionals by Cui [95]. 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 [87,96,97,98], 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 [99,100]: Without dispersion correction, BLYP over-structures liquid water and underestimates its diffusivity (i.e., the dynamics becomes too sluggish) [101,102]; the density of the liquid is also significantly under-predicted [102]. 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 [102]. 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 [18]. 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) [103] and produces an overstructured liquid with too slow diffusion ($\sim 4\times$ below experiment) [103,104]. Adding D3 brings density and the O–O RDF close to experiment [102], 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 [105].

PBE [106]: 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 [103]. 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 [18].

Advantages: Well-validated; widely available; reasonable lattice constants for ice ${\mathrm{I}}_{\mathrm{h}}$ [103].

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

revPBE-D3 [109,110]: 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 [111]. This performance partly reflects a cancellation between intrinsic functional error and the neglect of nuclear quantum effects in classical-nuclei trajectories [111]. 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 [107]. 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$Å [112] 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 [105], 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. [113] 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$ [114]	$+0.03$ [114]	$5.19$ [105]
BLYP-D3	$+6.6$ [114]	$-0.02$ [114]	$5.53$ [105]
PBE	$-12$ [102]	$-0.03$ [103]	$9.15$ [105]
PBE-D3	$+5.5$ [114]	$-0.07$ [114]	$9.62$ [105]
revPBEa	$-31$ [113]	$+0.08$ [113]	$7.11$ [105]
revPBE-D3 [109,110]	$-6$ [115]	$+0.01$ [115]	$8.32$ [105]

^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$Å [112] 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 [105], 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. [113] 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$ [114]	$+0.03$ [114]	$5.19$ [105]
BLYP-D3	$+6.6$ [114]	$-0.02$ [114]	$5.53$ [105]
PBE	$-12$ [102]	$-0.03$ [103]	$9.15$ [105]
PBE-D3	$+5.5$ [114]	$-0.07$ [114]	$9.62$ [105]
revPBEa	$-31$ [113]	$+0.08$ [113]	$7.11$ [105]
revPBE-D3 [109,110]	$-6$ [115]	$+0.01$ [115]	$8.32$ [105]

^a

4.2. Meta-GGA Functionals

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

SCAN [116]: 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 ${\mathrm{I}}_{\mathrm{h}}$ at ambient conditions and reproduces many structural, electronic, and dynamic properties of liquid water without empirical dispersion corrections [117], and correctly predicts ice ${\mathrm{I}}_{\mathrm{h}}$ as more stable than ice ${\mathrm{I}}_{\mathrm{c}}$ along with reasonable melting properties of both polymorphs [118]. However, SCAN overestimates the liquid water density by approximately 5% relative to experiment [117] and can exhibit numerical grid sensitivity in some implementations [119,120].

Advantages: No empirical parameters; good ice energetics including correct ${\mathrm{I}}_{\mathrm{h}}$/${\mathrm{I}}_{\mathrm{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 [120]: 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 [121]. The first O–O peak position of r²SCAN (2.865 Å) sits slightly above the experimental value (2.801 Å), indicating mild residual over-structuring [121]. Furness et al. report that r²SCAN converges with fewer iterations than SCAN [120].

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

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

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 [122].

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 [122]; additional computational overhead; condensed-phase MD applications rely on derived many-body potentials (e.g., MB-SCAN(DC)) rather than direct AIMD [122].

4.3. Hybrid Functionals

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

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 [124].

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

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 [130].

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

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

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. [87,124,132], 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 [124,132].

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 [87,124,132]; balanced description of covalent bonding and hydrogen bonding such that competing nuclear quantum effects on the hydrogen bond cancel correctly [124].

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 [130]; NQE treatment essential for quantitative agreement of dynamics and vibrational spectra [124,132].

$\omega$B97X-V [133]: 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 [134] 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 [135]. A closely related meta-GGA, B97M-rV (B97M-V with the rVV10 non-local correlation [136]), 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 [132], 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 [135]; built-in VV10 non-local correlation, so no separate dispersion correction is required [133,137]; the related meta-GGA B97M-rV reaches near-hybrid accuracy for bulk water at meta-GGA cost [132].

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 [138,139] ). Common corrections include:

• D3(BJ) [110,140]: Damped pairwise $C_6/R^6$ correction with Becke-Johnson damping. Well-validated, widely available.
• D4 [141] Improved atom-pairwise dispersion with charge-dependent coefficients.
• VV10 [137]: 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.

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 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 [71].

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, $E_0/k_{\mathrm{B}}T=\hslash \omega /\left(2k_{\mathrm{B}}T\right)$. For the O–H stretching mode in water ($\omega \approx 3600$ cm⁻¹ [72]), 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 [17]. 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 [17,131].

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. In standard path integral simulations of water, nuclei are treated as distinguishable particles [17], 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 [142]. Bosonic-exchange path integral methods do exist [143] but are unnecessary in this regime. Exchange effects are therefore not considered further in this review.

Competing quantum effects in hydrogen bonding. NQE do not simply weaken or strengthen hydrogen bonds; they produce competing effects [131,144]:

• Intramolecular ZPE (O-H stretch): Causes proton delocalization perpendicular to the bond, 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.

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).

T (K)	${\mathit{k}}_{\mathbf{B}}\mathit{T}$ (kcal mol−1)	ZPE(OH)/(${\mathit{k}}_{\mathbf{B}}\mathit{T}$)	${\mathbf{\Lambda }}_{\mathbf{H}}$ (Å)
300	0.60	∼8.5	1.0
250	0.50	∼10	1.1
200	0.40	∼13	1.2
150	0.30	∼17	1.4
77	0.15	∼34	2.0

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).

T (K)	${\mathit{k}}_{\mathbf{B}}\mathit{T}$ (kcal mol−1)	ZPE(OH)/(${\mathit{k}}_{\mathbf{B}}\mathit{T}$)	${\mathbf{\Lambda }}_{\mathbf{H}}$ (Å)
300	0.60	∼8.5	1.0
250	0.50	∼10	1.1
200	0.40	∼13	1.2
150	0.30	∼17	1.4
77	0.15	∼34	2.0

Across the entire supercooled and glassy range the zero-point energy of the O–H stretch exceeds the thermal energy by roughly an order of magnitude or more, so classical treatment of nuclei introduces qualitative errors in properties sensitive to hydrogen motion. The thermal de Broglie wavelength ${\Lambda }_{\mathrm{H}}$ approaching or exceeding the O–H bond length (∼1 Å) signals that protons cannot be treated as localized classical particles.

5.3. Isotope Substitution as an Experimental Probe of NQE

Nuclear quantum effects manifest experimentally as isotope effects: replacing H by D (or ¹⁶O by ¹⁸O) alters the zero-point energy and hence the measured structure, dynamics, and thermodynamics, which makes the H₂O/D₂O comparison the principal experimental route for validating quantum-nuclei simulations [131]. The competing intra- and intermolecular contributions introduced above (Section 5.1) are central to interpreting these data. Although the underlying changes are sizeable, with the O–H covalent bond roughly 3% longer than O–D and the associated hydrogen bond about 4% shorter, the opposing shifts largely offset, so that the net geometric isotope effect on the O–H distance in liquid water and ice ${\mathrm{I}}_{\mathrm{h}}$ is only ∼0.001 Å: the competing quantum effects almost exactly cancel under these conditions [131]. Water’s isotope effects consequently span a wide range across properties, from very large to almost negligible depending on the property and the thermodynamic conditions [17], a dichotomy that the competing-quantum-effects principle rationalises. Because reproducing this balance requires an accurate description of a broad region of the potential energy surface, isotope-resolved observables (geometries, proton kinetic energies, and fractionation ratios) provide a stringent simultaneous test of both the electronic-structure method and the nuclear treatment [131].

5.4. Structural Consequences of NQE

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

Density trends: Classical MD predicts monotonically increasing density upon cooling, whereas PIMD reveals a density maximum in low-density amorphous ice that is absent in the classical treatment [145]. This qualitative difference arises because ZPE-induced volume expansion counteracts thermal contraction at low T.

Hydrogen bond geometry: Classical simulations tend to produce overly structured, overly linear O–H···O hydrogen bonds. NQE-induced proton delocalization perpendicular to the O–H bond direction softens the intermolecular structure and yields more bent hydrogen bonds, broadening the O–O radial distribution function relative to the classical result [78,146].

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 broadens these peaks toward experimental values [78,146].

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 [124].

6. Conclusions and Outlook

Simulating water between 200 and 273 K confronts three coupled challenges: sampling degrades as the dynamics slow, nuclear quantum effects (NQE) grow from a correction into a dominant contribution, and structural and dynamical observables remain sensitive to the underlying potential energy surface. These cannot be treated in isolation; a reliable low-temperature simulation matches the rigour of the nuclear treatment to the rigour of the electronic structure. For a newcomer, this reduces to a small number of practical decisions.

Choosing the propagation scheme. BOMD on a well-converged Born–Oppenheimer surface is the natural baseline whenever NQE are smaller than the other modelling errors [17,71]. CPMD [75] is competitive in cost per step, but its fictitious-mass dependence and the difficulty of maintaining adiabaticity in the presence of the high O–H stretching frequency make it the weaker choice when low-temperature accuracy is the goal [77]. When the quantum statistics of the nuclei must be captured, imaginary-time path-integral methods (PIMD [79,80], accelerated through PIGLET [82,83] or ring-polymer contraction [81]) are the appropriate tools, with the bead number P scaling as $\beta \hslash {\omega }_{\max }$ and therefore rising as temperature falls. Approximate real-time methods (RPMD [90,91]; CMD and its quasi-centroid variant [93,94]) remain the route to vibrational spectra and reaction rates, but each carries characteristic artefacts (spurious resonances in RPMD and curvature-induced frequency shifts in CMD) that are amplified as T decreases [91,92].

Choosing the functional. Electronic-structure accuracy cannot be decoupled from the nuclear treatment. Dispersion-corrected GGAs such as revPBE-D3 owe part of their good ambient performance to a cancellation between the functional’s intrinsic error and the neglect of NQE in classical-nuclei trajectories [111]; this is a fortuitous balance, not evidence of accuracy at 200 K. Meta-GGAs (SCAN, r²SCAN) describe intermediate-range dispersion and ice-polymorph energetics well [116,117,118,120], but reproduce ambient liquid structure only when simulated at elevated temperature [121], a workaround that mimics the missing NQE rather than supplying them. Hybrid functionals combined with explicit path integrals (e.g. revPBE0-D3 + PIMD) are among the few combinations that quantitatively reproduce both structure and spectroscopy [124,131,132], at a cost that has historically restricted them to short trajectories on small cells.

Why the regime is special. Below roughly 250 K, NQE shift from a correction to an organizing principle: the ratio of O–H zero-point energy to thermal energy grows well past unity, and properties as basic as the temperature of maximum density of low-density amorphous ice cannot be recovered with classical nuclei [17,145]. Choices that look conservative at ambient conditions, such as classical nuclei or a GGA selected on room-temperature benchmarks, therefore become uncontrolled here. Standard path-integral simulations treat protons as distinguishable particles [17], an approximation kept robust for condensed water by the deep covalent localization of the O–H bond, in contrast to systems such as superfluid helium where exchange dominates [142]; bosonic-exchange path-integral methods exist [143] but are unnecessary in this regime.

Where the field is heading. A consequential recent development is the maturation of machine-learning interatomic potentials (MLPs) trained on high-level ab initio reference data [85,86,88,89]. By lowering the per-step cost of a force evaluation by orders of magnitude, MLPs make nanosecond-scale path-integral trajectories feasible, at the accuracy of their underlying ab initio reference, on the system sizes required for nucleation, phase-equilibrium, and amorphous-ice studies [87,96,97,98], and they are now the route by which most quantitative low-temperature results reach the literature. The principal caveat is transferability: an MLP trained on ambient data is not guaranteed to hold at 200–250 K or in amorphous phases, so explicit ab initio validation at the target conditions remains essential. Several gaps follow directly: direct DFT-AIMD benchmarks in the supercooled regime are sparse, functional rankings still rest largely on ambient data, and approximate quantum-dynamics methods need further development below 200 K. We expect MLP-driven, path-integral simulations, trained on hybrid or density-corrected reference data and validated against the structural targets of Section 2.2, to become the practical workflow of the coming decade.

Why it matters. Supercooled and amorphous water are the relevant phases for ice nucleation in mixed-phase clouds and in freezing atmospheric droplets [69,147], for the survival of organs under sub-zero preservation [7,8], for nucleation control in cryopreservation and food processing [10], and for the chemistry on amorphous-solid-water grain mantles in dense interstellar clouds that builds complex organic and prebiotic molecules [11,12,13,14]; amorphous and crystalline ices likewise dominate the surfaces of icy moons, comets, and Mars [15]. In each setting, prediction demands quantitative agreement with experiment at temperatures where NQE are non-negligible and direct measurement is fragmented. The convergence outlined above, path-integral sampling driven by machine-learning potentials trained on hybrid-functional reference data, brings these problems within reach of routine simulation at the system sizes and trajectory lengths the underlying physics demands.