1 Introduction
Aqueous systems are ubiquitous and fundamental to fields ranging from biology [
1−
5] and environmental [
6] science to industrial chemistry. Their essential role is governed by the unique properties of water, most notably its high dielectric constant and complex hydrogen-bonding network, which dictate critical behaviors in processes such as cellular function, protein dynamics, electrolyte chemistry, and catalytic reactions [
7,
8]. The interplay between water molecules and solutes determines key macroscopic properties such as ionic conductivity, chemical reactivity, and phase behavior. Therefore, achieving a predictive, atomic-level understanding of these systems bridges fundamental science with critical technological and environmental applications. However, unraveling these properties is severely hindered by a computational bottleneck: achieving both electronic-structure accuracy (to describe complex interactions) and macroscopic-scale sampling. The complexity is further compounded in practically relevant systems involving nanoconfinement, interfaces, and chemical reactions.
Computational simulations have become indispensable for probing atomic-scale details. Classical molecular dynamics (MD), which uses empirical force fields, can simulate large systems over long timescales, but it often lacks accuracy because empirical force fields frequently fail to capture quantum effects. In aqueous systems, for instance, MD may poorly describe the structure and dynamics of the hydrogen-bond network [
9], or proton transfer processes [
10,
11], or electronic structure [
12], leading to systematic errors in predicting properties like ion hydration free energies and reaction pathways. In contrast,
ab initio molecular dynamics (AIMD) relies on quantum mechanical principles to describe electronic interactions with high accuracy. Based on AIMD, researchers have achieved accurate predictions and gained deeper insights into fundamental aqueous properties, such as the solvation structure of ion [
12−
16], self-diffusion of water [
17−
19], and proton transfer dynamics in aqueous solutions [
20,
21]. Despite its precision, computational cost of AIMD is extremely high, typically limiting simulations to small systems (~100 water molecules) and short timescales (picoseconds), which is insufficient for studying rare events or large-scale collective dynamics. This long-standing accuracy-efficiency trade-off has been a major obstacle to reliable simulations of complex aqueous systems.
The emergence of machine learning potentials (MLPs) has revolutionized molecular simulations by bridging the gap between quantum accuracy and classical efficiency. MLPs employ machine learning algorithms, primarily neural networks, to construct high-dimensional potential energy surfaces (PES) from
ab initio data. A fundamental requirement for MLPs is invariance under translation, rotation, and permutation of identical atoms, achieved via rotationally invariant descriptors (e.g., atom-centered symmetry functions or bispectrum components) [
22]. Another key is decomposing the energy into local atomic contributions, enabling near-linear scaling. The Behler-Parrinello neural network (BPNN) [
23] pioneered the use of atom-centered symmetry functions to achieve rotational invariance and locality, enabling linear-scaling simulations of materials such as bulk silicon. Building on this, the Gaussian Approximation Potential (GAP) [
24] introduced Gaussian process regression with Smooth Overlap of Atomic Positions (SOAP) descriptors, improving the faithfulness of the MLP [
25]. To further enhance accuracy and automation, the Deep Potential Molecular Dynamics (DPMD) model [
26] employed deep neural networks to automatically learn optimal features from DFT data, achieving near-linear scaling and DFT-level accuracy for large systems. Around the same time, SchNet [
27] introduced continuous-filter convolutional networks on molecular graphs, demonstrating that graph-based representations could learn molecular properties with high accuracy and transferability. More recently, equivariant architectures have improved data efficiency and accuracy by incorporating physical symmetries directly into the network structure: NequIP [
28] uses E(3)-equivariant graph neural networks, while MACE [
29] introduces higher-order message passing network to capture complex many-body interactions with fewer parameters. These models, together with advances in handling long-range electrostatics (e.g., DPLR), provide a versatile toolkit for simulating aqueous systems with
ab initio accuracy over extended scales.
The computational advantage of MLPs over AIMD is most evident in their scaling behavior. AIMD, typically based on density functional theory (DFT), exhibits a cubic scaling
with number of molecules
due to the electronic structure problem, while well-designed MLPs achieve near-linear scaling
(see Fig. 1) [
26]. This reduction is achieved because the MLP inference step involves evaluating local atomic environments within a cutoff radius, bypassing the expensive global electronic minimization. Once trained, MLPs can predict energies and forces with
ab initio accuracy at a computational cost that is orders of magnitude less expensive than direct AIMD, enabling microsecond-scale simulations of systems containing tens of thousands of atoms. Techniques like active learning further allow MLPs to autonomously explore configuration space [
30,
31], enhancing their robustness and generalizability. Due to the high accuracy and efficiency, MLPs have emerged as a powerful and efficient tool for simulating complex aqueous systems, effectively overcoming the longstanding accuracy-efficiency trade-off.
This review aims to provide a comprehensive overview of the applications of MLPs in complicated aqueous systems, highlighting recent breakthroughs, methodological innovations, and persistent challenges. We will begin by exploring how MLPs have advanced our understanding of structural and dynamical properties in bulk water and aqueous solutions, such as ion hydration and water diffusion. Then we will discuss insights gained from machine learning molecular dynamics (MLMD) on confined aqueous systems, including water in confined systems and at interfaces, where finite-size effects and surface interactions dominate. Finally, we will examine the role of MLPs in probing reactive processes, such as proton transport, acid-base reactions, and catalytic mechanisms at interfaces. Although DFT serves as the primary reference for most MLPs, its intrinsic accuracy limitations — discussed further in the Outlook — must be acknowledged, with higher-level methods like quantum Monte Carlo (QMC) and CCSD(T) providing benchmark references.
2 Structural and dynamical properties of bulk systems
Building on the motivation for MLPs in aqueous systems, this section delves into the fundamental structural and dynamical properties of bulk water and aqueous solutions, which serve as a baseline for understanding complex interactions. MLP methods like DPMD enable ab initio level insights into hydrogen bonding networks, reorientation dynamics, and solute-solvent effects over nanosecond timescales, overcoming limitations of traditional simulations.
2.1 Bulk water
Bulk water constitutes the fundamental system for understanding aqueous environments, where its structural and dynamical properties govern behaviors in more complex solutions. Prior to the advent of MLPs, AIMD established a foundational understanding of the microscopic structure of liquid water. AIMD studies demonstrated that advanced density functionals, including hybrid functionals with exact exchange and van der Waals corrections, as well as the meta-GGA, could accurately predict the structure, thermodynamics, and electronic properties of water [
17−
19].
However, constrained by computational cost, AIMD simulations are typically limited to small system sizes and short timescales. An insufficiently large system size may give rise to finite-size effects, while an insufficiently long simulation time results in poor statistical convergence of dynamical convergence. Early benchmark studies revealed that for basic structural properties such as the O-O radial distribution function (RDF), simulations with about 64 water molecules were sufficient to eliminate finite-size effect [
32,
32]. While modeling the X-ray absorption spectra (XAS) of pure water requires a system size of 128 water molecules to achieve numerically converged results [
33]. For dynamic properties, the limitation of finite-size effects is more severe. Classical MD studies quantified that the self-diffusion coefficient of water exhibits significant finite-size dependence, underestimating values by ~10% even for systems of ~2000 molecules and necessitating hydrodynamic corrections (see Fig. 2) [
34]. Properties such as dielectric response, solvation free energies, and collective dynamics in electrolyte solutions also require very large simulation cells for numerical convergence [
35,
36].
For dynamical and time-correlation properties, achieving statistical convergence requires not only the mitigation of finite-size effects but also extensive sampling over long simulation times. In a system of 64 water molecules, obtaining a statistically converged diffusion coefficient necessitates at least 20–30 ps of trajectory [
17−
19], while properties involving rotational motions, which require well-averaged time-correlation functions, demand even longer simulation durations [
37,
38]. These requirements for large-scale, long-timescale sampling exceeded the practical limits of conventional AIMD, thus highlighting the need for a simulation method that could simultaneously deliver both electronic-structure accuracy and high computational efficiency.
Recent years, MLPs have emerged as a transformative solution, bridging the accuracy-sampling gap by extending the scope and depth of first-principles simulations in two key dimensions: (i) overcoming the severe sampling limitation, (ii) mapping atomic configurations to electronic and spectroscopic observables. MLPs overcome the sampling limitations of AIMD by enabling large-scale, long-time simulations, thereby allowing for the rigorous convergence of properties sensitive to system size and trajectory length. State-of-the-art MLP-driven simulations have provided definitive evidence for a dynamic, “dancing” distorted Eigen cation (H
9O
4+) as the prevalent species, resolving a long-standing structural debate through computationally intensive sampling that would be prohibitive for direct AIMD [
39]. In the study of heat transport, a deep neural network potential trained on DFT data enabled the calculation of thermal conductivity using the nonempirical SCAN meta-GGA functional, which significantly improved agreement with experiment — a task exceedingly difficult with direct AIMD due to the analytical complexity of the required energy flux (see Fig. 3) [
40].
Beyond addressing sampling limitations, MLPs dramatically broaden the spectrum of accessible properties by directly mapping atomic configurations to electronic and spectroscopic observables, moving far beyond the traditional output of forces and energies. This capability manifests on two interconnected fronts. First, MLPs enable the efficient prediction of electronic response properties. For instance, deep neural networks trained to predict molecular multipole moments from electron density allow for the calculation of dielectric constants and other susceptibilities from long trajectories, creating a direct bridge from structure to electronic response [
41]. Second, MLPs provide a framework for
ab initio-quality spectroscopy at molecular dynamics scales. To predict such spectroscopic observables, MLPs are trained not only on energies and forces but also on response properties computed from DFT (e.g., via finite-field methods). These quantities are local functions of atomic configurations, making them learnable using the same descriptors. A dual-network framework is often employed: one network learns the potential energy, while a separate network learns the dipole or polarizability. Once trained, the MLP can compute time-correlation functions over microsecond-long trajectories — a task far beyond direct AIMD. This yields statistically converged spectral line shapes and relaxation times, enabling direct comparison with experiment. For instance, such a dual-network framework modeling both the potential energy surface and electronic polarizability has enabled the calculation of Raman spectra over nanoseconds [
42]. This spectroscopic power extends across multiple techniques: MLPs have been used to compute infrared (IR), Raman, and sum-frequency generation spectra for direct comparison with experiment [
43]; to reveal how isotopic substitution (D
2O vs. H
2O) blue-shifts the X-ray absorption (XAS) edge due to a more ordered hydrogen-bond network [
44]. Most recently, by analyzing shifts in the OH-stretching bands, researchers successfully linked these spectral changes to the degree of ordering within hydrogen-bond chains [
45]. Collectively, these advancements transform MLP simulations into a versatile spectroscopic platform, yielding direct structural and dynamical insights that both validate and refine our understanding of the hydrogen-bonding network.
The ability of MLPs to deliver first-principles accuracy across extended scales not only solves sampling problems but also creates a new pathway to efficiently incorporate nuclear quantum effects (NQEs), which are central to water’s behavior. NQEs, arising from the quantum nature of atomic nuclei, are central to the behavior of water and far from minor perturbations. The importance of NQEs is unambiguously evidenced by isotopic differences, such as the ~20% higher self-diffusion coefficient of H
2O versus D
2O [
46], and by the NQE-induced softening of the oxygen-oxygen radial distribution function observed in neutron scattering [
47]. Accurately capturing these effects requires both an extremely precise potential energy surface and computationally feasible access to quantum statistical mechanics. Prior to the advent of MLPs, path-integral AIMD (PI-AIMD) was widely used to provide foundational insights of NQEs [
48]. However, the computational cost of PI-AIMD is very expensive. This stems from the path-integral formalism, which maps each quantum nucleus onto a ring of P classical replicas (“beads”), discretizing its Feynman path. Since each bead requires an independent electronic structure calculation, the cost per MD step scales linearly with P, which is typically 16−32 for aqueous systems. Even though applying colored-noise thermostat could enhance sampling efficiency [
49−
51], when combined with the necessity for a smaller integration timestep, the total computational effort escalates by two to three orders of magnitude compared to an AIMD simulation of the same physical system and duration [
52].
The integration of machine learning potentials with path-integral molecular dynamics (MLP-PIMD) has emerged as a transformative strategy to overcome this barrier. A landmark demonstration of this power is the work of Cheng
et al. [
53], which employed an MLP trained on the SCAN functional to perform extensive PIMD simulations. This study achieved an
ab initio thermodynamic description of water, accurately predicting its phase diagram while quantitatively capturing the isotopic shift between H
2O and D
2O — a direct and stringent test of NQEs. This exemplifies how MLP-PIMD enables quantum-accurate simulations at scales previously inaccessible, fundamentally extending the scope of first-principles investigation into nuclear quantum phenomena.
This capability to reliably incorporate NQEs has opened new avenues for probing their influence on dynamical properties. For instance, a detailed analysis across a temperature range shows that NQEs slow down the self-diffusion and rotational motion of water molecules, with intermolecular effects diminishing as temperature rises while intramolecular effects persist [
54]. Further dynamical insights emerge from studies of the dielectric response of water under pressure, where a nonlinear increase of the dielectric constant is attributed to pressure-induced strengthening of the hydrogen bond network, enhancing collective dipole fluctuations [
55].
This systematic progression — from establishing structural benchmarks and hydrogen-bonding networks, through elucidating NQE dynamics and pressure responses, to spectroscopic validations — demonstrates how MLPs provide a unified framework for unraveling the complexity of bulk water. These foundational insights create the necessary basis for exploring more intricate aqueous solutions in the following subsection.
2.2 Aqueous solutions
MLPs extend their predictive power to aqueous solutions, unraveling how ions reshape structural and dynamical behavior of water. This section summarizes key insights across diverse systems, organized by solute-induced perturbations.
In this review we mainly focus on two fundamental and widely debated questions regarding aqueous solutions: (i) do ions influence water’s structure and dynamics through short-range or long-range effects? (ii) how does concentration effect influence the structural and dynamical properties? The first concerns the spatial range of ion-induced perturbations. Since the early 2000s, experimental findings have presented a stark contrast. Evidence supporting a long-range influence comes from neutron diffraction studies, dielectric spectroscopy, and femtosecond elastic second-harmonic scattering experiments [
56−
58]. Conversely, femtosecond infrared spectroscopy and X-ray absorption experiments have concluded that ion-induced perturbations are confined to the first solvation shell [
59,
60]. Classical MD simulations, whose outcomes are sensitive to the empirical force field employed [
61], fail to resolve this debate. Meanwhile, AIMD, while providing first-principles accuracy, lacks the computational throughput to simulate the large systems and long timescales necessary to gather statistically robust evidence for long-range correlations [
9]. Notably, Gaiduk
et al. [
62] highlighted that addressing this question may require moving beyond single-ion models to consider ion-pair effects. This context is ideally suited for MLPs, which uniquely combine AIMD-level accuracy with the computational efficiency needed to tackle this long-standing puzzle.
To address this question, Li
et al. [
63] investigated intermediate-range of Mg
2+ ions. By employing an MLP-based approach, they compared the properties of MgCl
2 and MgBr
2 solutions at varying concentrations. Their simulations revealed that the retarding effect of Mg
2+ on water dynamics extends beyond the first hydration shell. This point was further proved by the results of hydration number, which is consistent with the experiment (see Fig. 4) [
38]. This medium-range effect was attributed to a competition between ion-water and hydrogen-bond interactions. This finding provides a new perspective for understanding the dichotomy between short-range and long-range perturbation theories. Furthermore, their work demonstrated that the choice of the underlying density functional significantly influences the predicted dynamical properties [
63].
We now turn to the second question regarding concentration effects. Experimental observations dating back nearly a century have shown that the static dielectric constant of aqueous NaCl solutions decreases with increasing salt concentration [
64]. Subsequent quantitative measurements further revealed that this decrease follows a distinct nonlinear trend [
65−
67]. However, the molecular origin of this nonlinear behavior has remained elusive. To resolve this, Zhang and colleagues leveraged an MLP approach to simulate NaCl solutions across a wide concentration range. Predicting the dielectric response of aqueous solutions demands both large system sizes (hundreds of molecules) and nanosecond-scale trajectories to get numerically and statistically converged dipolar correlations, yet AIMD simulations of such scale are computationally prohibitive. Moreover, standard MLPs with local cutoffs cannot intrinsically capture the long-range
tail of dipole−dipole and ion-water interactions, which are critical for accurate dielectric constants. To overcome these limitations, the authors employed the Deep Potential Long-Range (DPLR) framework [
68], which augments a local deep potential with a physics-based electrostatic term that explicitly treats long-range interactions beyond the cutoff. By classifying water molecules into distinct subpopulations (e.g., bulk and hydration water) and quantitatively deconvoluting their contributions to the total system dipole, their work demonstrated the nonlinearity to the evolving coupling in dipole correlations between bulk water and hydration water (see Fig. 5) [
69].
In addition to probing concentration effects, the capability of MLPs to efficiently perform simulations across multiple concentrations enables the routine screening of solution properties, thereby facilitating the investigation of a broader spectrum of problems in aqueous chemistry. For NaCl solutions, MLP simulations has shown that the structure perturbation originates from concentration effects is not equivalent to that originates from applying pressure, demonstrating distinct hydrogen-bonding network disruptions that cannot be mapped to pressure-induced compression [
70].
In recent years, water-in-salt electrolytes (WISEs) have attracted significant research interest due to their promising applications. Wang
et al. [
71] employed MLP-based molecular dynamics to demonstrate that the high reductive stability stems from a fundamental switching of redox levels between anions and water. This critical inversion, driven by the isolation of water molecules, suppresses the hydrogen evolution reaction and promotes the formation of a protective, anion-derived solid-electrolyte interphase [
71]. Beyond thermodynamic stability, MLPs also elucidate the unique transport dynamics in WISEs. In concentrated Zn
2+ systems, a subsequent study revealed a sophisticated three-stage carrier transport mechanism, progressing from independent ion diffusion to strongly correlated motion, and finally to the diffusion of small positive charges through anionic polymeric clusters [
72]. Together, these works elucidate the power of MLPs in unraveling the coupled thermodynamic and kinetic phenomena that govern the performance of next-generation electrolytes.
Beyond these equilibrium properties, MLPs provide profound insights into the dynamical behavior and transport mechanisms governing ionic solutions. Utilizing DPMD, Avula
et al. [
73] successfully simulated the anomalous diffusion of water in CsI and NaCl solutions, obtaining results that are qualitatively consistent with experimental observations across a range of concentrations. The researchers employed DPMD to elucidate the competitive mechanism by which different ions exert opposing effects on water diffusion. Long-timescale MLP simulations further resolve NaCl dissociation kinetics, revealing a stochastic pathway involving dynamic intermediates [
74], while potential of mean force (PMF) calculations show minimal energy differences between contact ion pairs and solvent-shared ion pairs [
75]. Similar PMF-based analyses have been extended to divalent metal ions interacting with biologically relevant ligands. Using neural network potentials, recent studies systematically characterized the binding of Mg
2+, Ca
2+, and Zn
2+ to amino acid side-chain analogues and the Mg
2+-acetate ion pair in water, revealing ion-specific coordination structures and a ~5.8 kcal/mol barrier for monodentate-bidentate conversion, while also providing benchmarks for classical force fields [
76,
77].
3 Properties of confined and interfacial aqueous systems
The exploration of water and electrolytes under spatial confinement and at solid interfaces has long been hindered by the dual challenges of experimental characterization and computational cost. The advent of MLPs is now overcoming these barriers, enabling a precise and comprehensive dissection of the unique physical chemistry that emerges when water is taken out of its bulk environment.
Water in nanoscale confinement is ubiquitous in nature. Such extreme spatial confinement can profoundly alter the dynamic transport behavior of liquid water, with significant implications for a wide range of biological and industrial processes [
78,
79]. Over the past decade, confinement effects have been recognized as a key factor modulating liquid transport. Experiments have demonstrated that water can flow through nanochannels at ultrafast speeds. This enhanced flow surpasses the predictions of classical Navier-Stokes hydrodynamics by over an order of magnitude. This ultrafast transport is highly dependent on the confinement dimension and the material of the channel walls, being particularly pronounced within carbon-based nanostructures [
80,
81]. In contrast, it is markedly less evident in systems such as those confined by MoS
2 or hexagonal boron nitride (hBN) [
82,
83]. The prevailing experimental interpretation of this phenomenon has largely relied on macroscopic hydrodynamic arguments [
84], attributing the differences to wall hydrophobicity and atomic-scale smoothness, which collectively determine the interfacial friction and thus the transport properties [
85]. However, these continuum models fail to account for the anomalous water transport observed specifically in graphene-based channels. While classical MD simulations have been employed to probe the role of interfacial friction, their results are notably sensitive to the empirical force field chosen, often leading to significant discrepancies with experimental measurements [
86]. To systematically unravel the underlying confinement effects, an
ab initio accurate description of the water-confining material interaction is essential [
87]. Herein lies the key strength of the MLP approach, which provides a powerful and accurate means to investigate how nanoconfinement influences water transport properties.
Recent advancements leveraging MLPs have clarified the dependence of nanofluidic transport on channel material and confinement geometry. First, addressing the long-standing puzzle of material specificity, a seminal MLP study directly compared water flow in carbon nanotubes (CNTs) and boron nitride nanotubes (BNNTs) [
88]. It quantified
a ≈ 5 times lower friction on carbon surfaces, tracing the origin to atomic-scale mechanisms: facile oxygen motion in CNTs versus specific hydrogen-nitrogen interactions in BNNTs. This moves the explanation beyond macroscopic hydrophobicity. At a more fundamental level, MLPs are quantifying how confinement distorts the hydrogen-bond network itself. Simulations reveal that hydrogen bond orient preferentially parallel to the interface, imparting a “two-dimensional” character that significantly reduces the water layers needed for network percolation compared to bulk water [
89]. This structural insight directly accounts for the geometric dependence of transport. Another pivotal study on graphene nanoslits revealed that the interlayer distance (
H) dictates dynamics through a competition between parallel flow and perpendicular fluctuations, with a critical transition at
H ≈ 12.5 Å (see Fig. 6). Importantly, this mechanistic understanding yields a rational design principle: transport can be enhanced several-fold by optimizing friction at specific confinement sizes [
90].
MLPs have uniquely enabled the exploration of the novel physical properties of confined water. Zhao, Qiu, and Guo [
98] developed a deep neural network potential for water confined in graphene nanocapillaries, trained on DFT calculations, which achieves near-DFT accuracy at substantially reduced computational cost and accurately reproduces a wide range of structural and dynamical properties of nanoconfined water. A landmark study revealed that water confined to a monolayer exhibits a rich phase diagram, including a superionic phase with proton conductivity exceeding that of battery materials — a phenomenon dictated by nanoscale pressure and inaccessible to conventional sampling [
99]. Furthermore, MLPs have elucidated how such confinement fundamentally alters collective electronic responses. Simulations of water between graphene and hBN directly connected the spontaneous polarization of interfacial water to the system’s anomalously low dielectric constant, providing a mechanistic link between subnanometer structure and macroscopic dielectric properties [
100].
MLPs are uniquely suited to address the fundamental challenges of understanding aqueous solutions in confinement. For NaCl solutions in graphene slits, MLP simulations have demonstrated that increasing confinement not only reduces ionic conductivity but also triggers a fundamental shift in the ion transport mechanism from a collective process to vehicular motion [
101]. Furthermore, the effective interactions between ions are profoundly altered. PMF calculations, now feasible with
ab initio accuracy over relevant distances, reveal that ion-ion interactions in nanoconfinement deviate significantly from bulk solutions, explaining anomalous phase behavior and ion selectivity [
102].
Among interfacial aqueous systems, one of the most widely studied systems is the electric double layer (EDL). The EDL plays a central role in electrochemical energy storage, catalysis, and corrosion, yet its atomic-scale structure and potential distribution have long remained elusive to direct experimental probing. Conventional AIMD is severely limited by its computational cost when simulating large, realistic electrode-electrolyte interfaces over relevant timescales. Recent MLP developments, particularly those incorporating long-range electrostatics, have opened new avenues for EDL research. For instance, Li
et al. [
103] combined AIMD with MLP simulations to reveal how alkali metal cations in base alter the hydrogen-bond network continuity of interfacial water within the EDL, thereby modulating the proton transfer kinetics in alkaline hydrogen electrocatalysis [
103]. Zhang and coworkers employed a DFT-based MLP with explicit long-range electrostatic corrections to systematically investigate the anatase TiO
2-electrolyte interface under different pH conditions. Their large-scale simulations captured water dissociation/recombination and proton transport, showing that the higher affinity of cations for the oxide surface under basic conditions leads to an increased differential capacitance, in quantitative agreement with experiment [
104]. Zhu and Cheng further developed a hybrid MLP framework that simultaneously treats local electronic polarization in the electrolyte and non-local charge transfer in the metal electrode. This approach successfully reproduced the bell-shaped differential Helmholtz capacitance at the Pt(111)-electrolyte interface and provided new insights into the dielectric response of the EDL [
105]. Collectively, these studies demonstrate that MLPs can effectively unravel the complex behavior of EDLs at atomic resolution.
At solid−liquid interfaces, the subtle interplay between water and the atomic specifics of the substrate dictates structure and dynamics. While AIMD has been limited to simulating small patches for short durations, MLPs now enable direct comparison with experimental observables across large, realistic interfaces. For example, MLPs combined with experimental atomic force microscopy can differentiate the distinct hydrogen-bonding structures of water confined between graphene and hexagonal boron nitride surfaces — a task beyond the reach of conventional force-field molecular dynamics [
106]. Beyond static structure, MLPs reveal how mechanical coupling governs transport. Simulations of flexible channels show that the solid−liquid friction in materials like graphene scales inversely with lateral size — a key relationship absent in rigid models and one that opens new avenues for flow control [
107]. Similarly, MLP simulations have resolved the long-standing question of hydronium (H
3O
+) and hydroxide (OH
−) ion distribution at the graphene-water interface, linking it to interfacial polarization and charge rearrangement effects [
108]. These studies exemplify how MLPs bridge the gap between nanoscale simulation and macroscopic experiment, providing mechanistic insights that are both atomically precise and directly relevant to observable interfacial phenomena.
In summary, MLP-based simulations are not merely incremental improvements but represent a paradigm shift in the study of confined and interfacial aqueous systems. By bridging the gap between the accuracy of quantum mechanics and the scale of classical molecular dynamics, they are systematically unraveling the complex, atomic-level physics that govern these technologically critical environments, from energy storage and catalysis to nanofluidics and desalination.
4 Reactions in aqueous systems
The accurate simulation of chemical reactions in aqueous environments represents one of the most demanding challenges in computational chemistry. These processes, which include proton transfer (PT), acid dissociation, and interfacial catalysis, are characterized by concerted atomic motions, significant solvent reorganization, and the making and breaking of chemical bonds. Capturing these processes requires both the electronic-level accuracy necessary to describe bond rearrangements and the ability to sample rare events and achieve statistical convergence for reliable kinetics. MLPs enable long-timescale simulations at ab initio accuracy. This capability allows for the rigorous mapping of multi-step reaction pathways, the quantification of free energy landscapes, and the direct connection of simulated dynamics to experimental observables. This section showcases the transformative impact of MLPs across a spectrum of aqueous reactivity, from elucidating the intricate mechanism of proton transport to unraveling complex reaction dynamics at interfaces and in bulk solution.
AIMD has laid the foundation for understanding PT in aqueous systems, capturing key aspects of the Grotthuss mechanism through visualization of hydrogen-bond network rearrangements and the dynamics of hydrated proton complexes [
21]. These studies have provided insights into proton hopping mechanisms and the quantification of free-energy barriers [
20]. However
, AIMD’s severe computational constraints prevent adequate sampling of complete PT coordinates and reaction kinetics. A critical illustration of this limitation is the finding that reliable estimation of proton and hydroxide (OH
−) ion diffusion coefficients requires simulations spanning hundreds of picoseconds — a timescale far beyond routine AIMD capabilities [
109]. Therefore, there is a pressing need for efficient and highly accurate sampling methods, enabled by MLPs, to investigate PT events.
MLPs enable the precise quantification of PT mechanisms in bulk systems. A seminal study leveraged DPMD to unravel a three-step mechanism for proton transport in bulk water, which is governed by two consecutive hydrogen-bond exchanges (see Fig. 7) [
110]. This finding provides unprecedented statistical confidence in the detailed reaction coordinate of the Grotthuss mechanism. Furthermore, MLPs have proven crucial for mapping PT pathways at complex interfaces. For instance, DPMD simulations of the air-water interface enabled the computation of the SFG spectrum of hydronium (H
3O
+) ions, directly connecting simulated interfacial structures with experimental spectroscopic signatures [
111]. Researchers have further extended their investigations of proton transfer to aqueous solutions. Using DPMD, Zhao
et al. [
112] has found that the dissociation kinetics of H
2CO
3 are governed not by the inherent stability of carbonic acid, but by the hydrogen-bonding ring structure of specific conformational isomers, underscoring the critical role of solvation topology (see Fig. 8). Collectively, these studies have elucidated the PT mechanism.
While proton transfer remains foundational in bulk systems, it is the processes at surfaces and interfaces, most notably the air−water interface, that are often more critical in practical scenarios and represent a domain well-suited to the strengths of MLP simulations. DPMD simulations of the air-water interface revealed a decoupling between water and proton dynamics: while water molecular diffusion accelerates, proton diffusion slows down considerably [
113]. This contrast stems from reduced hydrogen bond coordination at the interface, which facilitates water mobility but impedes the stable proton hopping central to the Grotthuss mechanism. MLP-based reactive force fields have also resolved long-standing debates about interfacial acidity. Studies on nitric and formic acids at the air−water interface showed that the measured degree of dissociation is highly sensitive to the experimental probing depth and system size [
114]. The simulations revealed that acidity drops at the interface due to desolvation but can be enhanced just below it, offering a molecular explanation for contrasting experimental results. This concept extends to water droplets, where MLP simulations showed that although pH is uniform in the core, substantial concentration gradients of H
3O
+ and OH
− ions across the outermost molecular layers lead to increased acidity at the interface [
115].
Moving beyond proton hopping, MLPs have enabled a deep dive into the mechanisms of complex chemical reactions across diverse aqueous environments. At fluid interfaces, MLP simulations have uncovered unconventional reaction paradigms, benefiting from their high-efficiency sampling of molecular configurations [
116]. A prime example is the discovery, via MACE-based MD simulations, of an “in-and-out” mechanism for CO
2 hydration at the air-water interface, wherein CO
2 diffuses into the aqueous surface layer to react and the product is subsequently expelled; this confirms a significant contribution from surface-adsorbed CO
2 to the overall acidification rate [
117]. Similarly, the dissociation dynamics of N
2O
5 on water surfaces are dominated by interfacial hydrolysis and competitive re-evaporation rather than bulk-solution processes [
118]. In bulk solution, MLPs coupled with enhanced sampling provide quantitative insights into reaction free energies and mechanisms. This is exemplified by the study of weak carboxylic acids, where MLPs successfully captured the concerted proton dissociation of formic and acetic acid to yield accurate p
Ka estimates and characterized their dimerization processes [
119]. MLPs also illuminate reactions of fundamental importance in synthesis, as demonstrated by the exploration of the amide bond formation mechanism in solution [
120].
5 Outlooks
This review highlights the versatility of MLPs and their importance across multiple aqueous systems — including bulk water, aqueous solutions, confined systems, and interfaces — due to their ability to capture water−water, ion−water, and water−interface interactions with ab initio accuracy. Although the MLP approach has provided solutions to long-standing problems concerning accuracy and efficiency, fundamental questions persist regarding its applicability to the most complex aqueous scenarios.
A fundamental open question remains: do ions exert genuine long-range effects, and if so, how do these effects influence the structure and dynamics of water? Resolving this issue crucially depends on the development of MLP approaches capable of capturing long-range interatomic interactions, which still remains difficult. To date, most models utilize local descriptors that only capture the atomic environment within a cutoff radius, implicitly attributing all interactions to atoms within this sphere. Due to computational constraints, this cutoff cannot be extended arbitrarily, thus precluding the modeling of true long-range physical effects. Recent efforts fall into several categories. One strategy involves incorporating physics-based corrective terms for specific interactions beyond the cutoff, such as electrostatics (e.g., DPLR or CACE). However, this approach requires prior knowledge of the functional form of the long-range interactions. Another strategy employs message-passing neural networks to propagate information over longer ranges (e.g., SO3krates, MACE) [
29,
121]. Yet, the computational cost of these models escalates significantly with the number of message-passing layers. While the latest research explores capturing information in Fourier or spherical harmonic space [
122], a general and efficient model for long-range interactions remains elusive.
As mentioned in the introduction, the accuracy of MLPs is inherently bounded by the DFT reference data on which they are trained. While DFT remains the most practical electronic structure method for generating large training sets, its well-known limitations — including self-interaction error, incomplete treatment of van der Waals dispersion, and approximate exchange-correlation functionals — can propagate into MLP predictions for hydrogen-bonding networks, proton transfer barriers, and interfacial charge distributions. If the reference functional fails to accurately describe these key interactions, MLPs trained on such data will inevitably inherit the same deficiencies, regardless of their numerical efficiency. To overcome this limitation, the community is increasingly turning to higher-accuracy benchmarks such as quantum Monte Carlo (QMC) and CCSD(T) [
123−
128]. Both methods provide near-exact energies and forces for small water clusters and reaction intermediates, albeit at high computational cost. Although still too demanding for direct dynamics, QMC/CCSD(T) data can serve as a reference for training or correcting MLPs, as discussed below.
Beyond the quest for higher-accuracy reference data, another major challenge lies in the efficient generation of training datasets that comprehensively sample the complex configurational space of aqueous systems. Practical applications of aqueous solutions often involve multi-element systems, solid−liquid or gas−liquid interfaces, and long-range interactions, all of which make it prohibitively expensive to acquire sufficient training data using traditional DFT methods. To address this, current strategies focus on improving data efficiency. One powerful approach is active learning, which iteratively explores configuration space using ML-driven MD rather than
ab initio MD, identifies regions where the MLP has high uncertainty, and selectively adds new training points from those regions [
30]. This minimizes the amount of DFT data required while maximizing model robustness. Other methods include Δ-learning, where a low-cost MLP learns the correction between a cheaper baseline method and a higher-level theory (e.g., DFT vs. QMC/CCSD(T)) [
129], and multi-fidelity approaches that combine data from multiple levels of theory [
130]. Additionally, the emergence of universal foundation models pretrained on extensive datasets offers a promising route: such models can be fine-tuned for specific aqueous systems with limited additional data [
131,
132]. Despite these advances, a systematic and robust solution for generating training sets for highly complex solutions and chemical reactions remains an open challenge.
Once the aforementioned challenges are overcome, MLPs will evolve into a transformative tool, not only deepening our fundamental understanding of the mechanisms of aqueous systems but also directly contributing to the rational design of materials and processes to address critical challenges in environmental science and industrial chemistry [
133].