Perspective: Atomistic simulations of water and aqueous systems with machine learning potentials

A large fraction of the surface of the Earth is covered by water and, still, some ice, giving our planet its distinctive blue color when viewed from space. Water is carried down deep into the Earth’s crust at subduction zones, influencing volcanism and plate tectonics, and in the atmosphere, in the form of vapor, liquid, or ice, water is a key climate factor from the troposphere up to the stratosphere and mesosphere. Down at the Earth’s surface, water shapes landscapes, provides the basis for life, and is central to many technologies that sustain humanity. Given its significance and abundance, it is no surprise that over the centuries, much research has been undertaken to understand the properties of water and their physical origin.

One of the central scientific questions addressed in water research is how the complex behavior of water, exhibiting many anomalies and a rich phase diagram, arises from the interactions of the chemically rather simple H₂O molecules. Due to the limited temporal and spatial resolution of many experimental probes, much of what we know about water has been learned from computer simulations. Specifically, atomistic simulations have provided detailed insights into the directed network of hydrogen bonds between molecules that governs the structure and dynamics of water and its interaction with solutes and surfaces.^1–4 Moreover, computer simulations have made it possible to investigate water under extreme conditions that are not accessible in experiments. For instance, simulations have been used to study water and ice under pressure and temperature conditions prevailing in the deep Earth⁵ and in the interiors of the giant planets Uranus and Neptune,⁶ as well as in the deeply supercooled state, the so-called no-man’s land, where crystallization occurs extremely quickly.^7,8

Following the pioneering Monte Carlo (MC) simulations of Barker and Watts⁹ and molecular dynamics (MD) studies of Rahman and Stillinger¹ in the late 1960s and early 1970s, respectively, many computer simulations of water and aqueous systems were carried out. Initially, these simulations were based on empirical potentials,¹⁰ but, later, they relied increasingly on forces and energies obtained from electronic structure calculations.¹¹ In the empirical potential—or force field—approach, the functional form of the interaction potential is constructed to capture the main physical interactions between molecules, with parameters adapted to reproduce some experimentally known quantities and/or quantum mechanical reference data.

Since the first water models for MD simulations were proposed more than half a century ago,^12–14 a vast number of empirical potentials were developed for water and ice, ranging from simple forms based on pair interactions to sophisticated many-body potentials, including polarization and charge transfer.^10,15–21 Despite their often simple functional form, these models have been remarkably successful in capturing the key properties of water across the phase diagram.^22,23 Modeling chemical reactivity, however, has proven difficult using empirical potentials, and with a few exceptions,^24–27 empirical water models are usually non-reactive. That is, they lack the capability to represent the dissociation and recombination of water molecules, which is, for example, essential for describing the famous Grotthuss mechanism of proton transport.²⁸ The inclusion of proton transfer events is not only central for simulations of acids and bases but is also of utmost importance in countless chemical reactions with water involved as a reactant, from the hydrolysis of biomolecules to electrochemical water splitting. Nevertheless, to this day, empirical force fields have been used for the vast majority of simulations involving water, particularly in studies of biological macromolecules, where aqueous solvation effects are of crucial importance.¹⁸ 

A more fundamental route to the computer simulation of water relies on determining interactions from first principles, i.e., ab initio by solving the electronic Schrödinger equation. The first ab initio molecular dynamics (AIMD) simulations of liquid water were carried out 30 years ago based on density functional theory (DFT). Since then, ab initio methods have been used extensively to study water^11,29 and, more generally, aqueous systems.³⁰ While, in principle, ab initio approaches have the potential for truly predictive simulations and also provide access to electronic properties, currently available approximate methods are still somewhat limited.³¹ For instance, predicting the equilibrium density of liquid water has been proven difficult within DFT based on standard exchange correlation functionals, and only the inclusion of dispersion forces produces satisfactory results.^30,32 Moreover, some attempts to study water using methods beyond DFT have been made, e.g., using MP2-based simulations³³ or employing the random phase approximation (RPA).³⁴ 

Compared to empirical force fields, ab initio approaches are computationally more expensive by many orders of magnitude, severely limiting the accessible system size and simulation times. Hence, many processes of interest occurring in aqueous systems, for instance freezing, glassy dynamics, the solvation of complex interfaces and biopolymers, as well as complex chemical reactions, are far beyond the capabilities of current ab initio simulations. Therefore, in order to transfer the reliability of ab initio methods to more complex systems, in the past two decades, considerable effort has been made to develop efficient but accurate potential energy surfaces (PESs) based on systematic and flexible functional forms.^35–38 In particular, the growing use of machine learning techniques, such as artificial neural networks (ANNs), and Gaussian processes for the construction of highly accurate and efficient machine learning potentials (MLPs) has revolutionized—and continues to revolutionize—the field of atomistic computer simulations.^39–53 These methodological advances in modern MLPs make it now possible to predict even complex properties of condensed systems from first principles, opening up exciting new possibilities in chemistry, materials science, and related disciplines.

Consequently, in recent years, atomistic modeling based on MLPs has also been increasingly applied to study water and, more generally, aqueous systems. The field has been evolving rapidly, in terms of both the underlying methodology and the complexity of systems that can be addressed, from small water clusters in vacuum via bulk water and its phase diagram to solution chemistry and processes at solid–liquid interfaces (cf. Fig. 1).

In this Perspective, we first provide a brief overview of the methodological basis and the current status of MLPs. In doing so, we focus on those types of MLPs that have been most frequently employed in simulations of systems involving water. In the subsequent sections, we discuss different applications of MLPs to aqueous systems. While it is not the goal of this Perspective to provide a comprehensive review covering all MLP-based studies of water and related systems exhaustively, we discuss a broad range of representative applications and point out future research directions to demonstrate the power and versatility of MLPs for the study of water and aqueous systems.

In recent years, machine learning potentials have become an increasingly important tool for atomistic simulations of complex systems in chemistry, physics, and materials science. As a consequence, the development of MLPs is a very active topic of research, and here, we will restrict our discussion to a concise overview about the current status of the field. Readers interested in further details are referred to a very large number of reviews covering all aspects of the methodology of MLPs.^39–54 

MLPs offer many advantages, like an excellent numerical agreement with the underlying electronic structure reference method, resulting in typical energy errors of only about 1 meV/atom and force errors in the order of 100 meV/Å. These errors are significantly smaller than the uncertainty, e.g., due to the choice of the exchange correlation functional in DFT, and thus, replacing electronic structure calculations by MLPs does only marginally affect the accuracy of the simulations. Moreover, MLPs can describe the making and breaking of chemical bonds and provide a rather high computational efficiency enabling simulations of systems containing many thousands of atoms. Still, they are usually about one to two orders of magnitude slower to compute than simple classical force fields. A main advantage of MLPs is their ability to transfer the predictive capabilities of ab initio approaches to large systems, which are required for aqueous systems to (1) achieve proper sampling of the liquid, (2) ensure bulk-like water properties at greater distances from solid surfaces, and (3) prevent artificially high electrolyte concentrations and/or artificial periodicity that would be present in relatively small simulation cells directly accessible by electronic structure calculations.

First, MLPs have been introduced more than a quarter of a century ago by Blank et al.,⁵⁵ who suggested to use a feed-forward neural network (NN) to represent the interactions between diatomic molecules and solid surfaces. A general limitation of this first generation of MLPs, which has been further explored by numerous groups for about a decade for different types of systems, has been the limitation to a few atomic degrees of freedom only, restricting the applicability to small molecules in vacuum or small molecules interacting with frozen surfaces, as summarized in some early reviews.^39,40 The major challenge for extending this method to high-dimensional condensed systems, such as liquid water, at that time has been the lack of suitable structural input descriptors for the machine learning algorithms, which ensure the imperative translational, rotational, and permutational invariances of the potential energy surface (PES). Only for some applications, system-specific approximate solutions could be derived when neglecting less important degrees of freedom of the system.^56–58 In parallel to these efforts in the development of early neural network-based MLPs, in pioneering work, Braams and Bowman^35,59 introduced permutation invariant polynomials (PIPs), which are closely related and enable the construction of very accurate PESs by linear regression based on symmetrized polynomials as basis functions. While PIPs do not employ traditional machine learning algorithms, such as neural networks, they show a similar flexibility and include all invariances exactly, which has been a frustrating challenge for early MLPs employing non-linear models. Still, PIPs share with early MLPs the restriction to small systems with a very limited number of degrees of freedom. Over the years, PIPs have enabled the construction of very accurate potentials and have been applied successfully, for example, to vibrations and reaction dynamics of small molecules in vacuum as well as a variety of water clusters.³⁵ 

The energetically relevant local environment determining the E_i is defined by a cutoff radius R_c such that all interactions beyond this radius, which is typically chosen between 5 and 10 Å, are not explicitly included. While the ansatz of Eq. (1) has been used in many empirical potentials for a long time, the introduction of second-generation MLPs has only become possible by the development of many-body descriptors with full translational, rotational, and permutational invariance. In the case of HDNNPs, atom-centered symmetry functions (ACSFs)⁶⁴ are most frequently used for this purpose, but at present, a wide range of alternative descriptors is available and employed in different types of MLPs.^65–68 All these descriptors represent structural fingerprints of the local atomic environments and serve as an input for the machine learning algorithms, which then construct the functional relation between the atomic environments and the atomic energies.

In the case of HDNNPs, which are often used for simulations of water and aqueous systems, there is one feed-forward NN to be parameterized per element, which is then evaluated as many times as atoms of the respective element are present in the system. Closely related is ANI (Accurate neural network engine for molecular energies ANAKIN-ME, ANI for short), which is a HDNNP with modified angular ACSFs⁶⁹ aiming for transferability across a wide range of organic molecules. Another NN-based MLP that is very frequently used for simulations of water is Deep Potential Molecular Dynamics (DeePMD),^70–72 which employs a local atomic coordinate system and descriptors in this reference frame as an input for the atomic NNs. Many other second-generation MLPs containing NNs have been proposed,^73,74 and in recent years also NN potentials learning descriptors as part of the training process employing message passing⁷⁵ have been put forward.^76–82 Beyond neural networks, Gaussian approximation potentials (GAPs)⁸³ combined with the SOAP (Smooth Overlap of Atomic Positions) descriptor⁶⁵ are among the most frequently used MLPs, with a few applications to aqueous systems reported to date. Many other second-generation MLPs are available in the literature, which can be expected to be used for systems containing water in the future.^68,84,85

An obvious limitation of second-generation MLPs is the truncation of atomic interactions at the cutoff radius. However, in many aqueous systems, long-range electrostatic interactions play an important role.⁸⁶ These are explicitly considered in third-generation MLPs, which include electrostatic interactions employing environment-dependent charges represented by machine learning models. Already in 2007, Popelier and co-workers showed that it is possible to construct electrostatic multipoles using neural networks⁸⁷ and Gaussian processes⁸⁸ to improve the description of electrostatics in classical force fields, and, in addition, applications to water clusters have been reported.⁸⁹ 

In 2011, HDNNPs of the third generation were proposed by introducing a second set of atomic neural networks providing atomic partial charges trained to DFT reference data.^90,91 From these charges, the electrostatic energy can be computed and combined with the short-range expression of Eq. (1) to yield the total energy of the system. By training the short-range part to represent only the energy component not covered by electrostatics, double counting of energy contributions can be avoided. Further MLPs, including long-range electrostatics, are, for example, the HDNNP TensorMol,⁹² the message passing network PhysNet,⁹³ and many others.^94–98 At present, the machine learning representation of atomic partial charges and electrostatic multipoles is a very active field of research, opening many routes to the construction of third-generation MLPs.

A remaining limitation of third-generation MLPs is the locality of the atomic charges, which does not allow us to describe systems exhibiting long-range charge transfer and other non-local dependencies between the geometric and electronic structures.⁵² These phenomena can be considered in fourth-generation MLPs. The first MLP of this generation has been the charge equilibration neural network (CENT) technique proposed by Ghasemi et al. in 2015.⁹⁹ Since the introduction of CENT, which employs a global charge equilibration step¹⁰⁰ and is intended for applications to ionic materials, several other fourth-generation MLPs have emerged, such as Becke population neural networks (BpopNNs),¹⁰¹ fourth-generation HDNNPs (4G-HDNNPs),¹⁰² and charge recursive NNs (QRNNs).¹⁰³ To date, fourth-generation MLPs have not been extensively applied to water but offer new interesting possibilities for studies of complex systems.

A key aspect in the training of MLPs⁵⁴ is the construction of suitable datasets covering the structures that are visited in the intended simulations, as MLPs often show a strongly reduced accuracy when extrapolating beyond the known part of configuration space. The size and composition of these datasets depend on the systems of interest, but due to the high flexibility of MLPs, often energies and forces of 10 000 or more electronic structure calculations are required for training reliable potentials. For a systematic and unbiased determination of these structures, often various forms of active learning are employed.^104–107 

Water clusters have received considerable attention already during the advent of MLPs as important benchmark systems. Even for small water clusters, there is a large structural variety with many energetically close local minima, posing a significant challenge for potential development. At the same time, their moderate size allows us to perform accurate high-level electronic structure reference calculations.

Early first-generation MLPs for water clusters include an MP2-based six-dimensional PES for the water dimer with frozen monomer geometries reported in 1997 that made use of a single feed-forward NN.¹⁰⁸ In 2006, a very accurate three-dimensional NN potential for the water monomer with a root mean squared error (RMSE) of only 1 cm⁻¹ (about 0.1 meV for the molecular potential energy) was published,¹⁰⁹ and, in addition, a NN potential for the same systems focusing on the permutation invariance of the PES was reported in 2012.¹¹⁰ From 2005 onward, water clusters have also been investigated in great detail using permutation invariant polynomials reaching MP2 and coupled cluster accuracies.^111–122 

To address larger water clusters, second-generation HDNNPs have been developed based on DFT data for a series of neutral clusters up to the decamer^123,124 as well as for several protonated clusters.¹²⁵ HDNNPs have also been applied to very large clusters containing hundreds of molecules.¹²⁶ Moreover, nuclear quantum effects (NQEs) in neutral and protonated water clusters have been studied in recent years using very accurate HDNNPs trained to coupled cluster data.^{104,127–131}

Starting in 2009, NNs were employed using water clusters as a test bed to find ways of improving the description of Coulomb interactions in classical force fields by learning environment-dependent electrostatic multipole moments.¹³² A comparison of NNs and Gaussian process regression (GPR) for the representation of multipoles found these methods to be similar in terms of accuracy and costs,⁸⁹ resulting in an extension of the work to a water molecule embedded in a water decamer.¹³³ Finally, this approach has been developed further to the FFLUX water model and applied to larger water clusters.¹³⁴ In related work, GPR has been employed in a similar way to express charges in water molecules.⁹⁵ 

A first MLP of the third generation expressing the full energy, i.e., electrostatics and short-range bonding, of the system by machine learning has been a third-generation HDNNP for the water dimer reported in 2012,⁹¹ which includes electrostatic interactions based on environment-dependent atomic partial charges expressed by a second set of atomic NNs. Another example in which electrostatics, as well as van der Waals interactions, were implemented on top of a short-range HDNNP for water clusters is TensorMol published in 2018.⁹² Moreover, small clusters have been used as a benchmark for the inclusion of long-range interactions in DeePMD.⁹⁷ 

Finally, machine learning was used to learn tensorial properties using water clusters as a benchmark,^135,136 and, in addition, the effect of noise on training HDNNPs for possible future applications on quantum computers has been tested for this system.¹³⁷ Beyond NNs and GPR, also Support Vector Regression (SVR) has been explored for water clusters.¹³⁸ Although not directly used in the development of high-dimensional PESs, the SVR method as well as random forest and Gaussian regression has also been used in specific applications, such as the prediction of the electron correlation energy in the water monomer and dimer.¹³⁹ Beyond atomistic potentials, machine learning methods have also been used in various forms in combination with electronic structure concepts and methods, such as coupled cluster theory with single, double and perturbative triple excitations [CCSD(T)] and variational quantum Monte Carlo calculations with benchmarks for monomers and small clusters.^140–146 

Of all aqueous systems, the construction of MLPs for bulk liquid water has received most attention not only due to its crucial importance as a solvent for a wide range of (bio)chemical reactions, but also because of its remarkable properties making water one of the most challenging benchmarks for the construction of interatomic potentials for molecular systems.

Due to the limited dimensionality that could be dealt with at that time, the first application of machine learning for the simulation of pure liquid water in 2002¹⁴⁷ aimed for the inclusion of polarization in the TIP4P water model.¹⁴⁸ For this purpose, a feed-forward neural network was used to represent the many-body interactions in dimers of rigid water molecules. This water model, named T4NN, was trained with MP2 reference data and was then used in Monte Carlo simulations to determine a range of properties of water, such as its density, heat capacity, and radial distribution functions. Overall, the agreement of many properties with experiment under standard conditions could be significantly improved with respect to the underlying TIP4P model, while the transferability to different temperatures posed a challenge.

In 2013, GAPs were used to enhance the accuracy of DFT-based ab initio MD simulations by introducing one- and two-body corrections trained on data of water monomers and dimers at the CCSD(T) level.¹⁴⁹ Their application to water clusters and MD simulations of liquid water showed an improved overall potential energy surface, but additional corrections beyond two-body terms were necessary for further improvements, as also demonstrated by comparison with quantum Monte Carlo data.¹⁵⁰ 

The first full-dimensional MLP for bulk liquid water and ice, which did not rely on a force-field or a DFT baseline potential, was reported in 2016.¹⁵¹ In this work, a series of HDNNPs were trained using DFT reference data obtained for different generalized gradient approximation (GGA) functionals, with and without dispersion corrections. This allowed us, for the first time, to benchmark the quality of common DFT functionals in the description of computationally demanding properties of water, such as the density anomaly, the melting temperature, the viscosity, and the dielectric constant. The calculation of such quantities requires extensive simulations of large systems, which are prohibitively expensive using ab initio MD directly, but become affordable with HDNNPs. In this work, particular attention was devoted to study the effect of van der Waals interactions, which were shown to govern the flexibility of the hydrogen bond network and, hence, play a crucial role in determining the properties of water and ice. In fact, if van der Waals forces are neglected, the density maximum of water disappears and ice becomes denser than liquid water. In follow-up work, HDNNPs based on BLYP-D3 and RPBE-D3 data were used to study the density anomaly of water at negative pressures (for both functionals)¹⁵² and the kinetics of the ice–water interface (only for RPBE-D3).¹⁵³ 

Shortly after the first full-dimensional MLP for water, a HDNNP trained with B3LYP + D3 data was used in conjunction with path integral MD simulations to study nuclear quantum effects (NQEs) of liquid water close to the triple point.¹⁵⁴ Since then, several other studies of nuclear quantum fluctuations in bulk liquid water and ice have followed,^155–158 including a thermodynamic stability analysis of liquid water as well as hexagonal and cubic ice employing hybrid DFT data,¹⁵⁹ and, more recently, a study of NQEs of liquid water based on the random phase approximation.¹⁶⁰ 

The past few years have witnessed a significant expansion in the use of MLPs for bulk water and ice. This growth encompasses many applications but also the development of methods, tools, and extensive benchmarking. For instance, the vibrational spectroscopy features of liquid water were extensively studied over the full frequency spectrum taking into account the effect of temperature and overcoordinated hydrogen-bond environments employing a HDNNP based on revPBE-D3.^161,162 Embedded atom neural networks (EANNs) have been used to represent tensorial properties in water describing vibrational features with the revPBE0-D3 functional.¹⁶³ Moreover, using water as a test example, DeePMD was introduced in 2018,⁷¹ which, alongside HDNNPs, emerged as one of the principal methods for modeling water using MLPs. DeePMD has also been coupled with empirical force fields¹⁶⁴ and employed to develop coarse-grained water models.¹⁶⁵ Some applications of DeePMD include the analysis of hydrogen bond dynamics in supercritical water,¹⁶⁶ the comparison of light and heavy water to assess isotope effects,^167,168 and the calculation of vibrational densities of states.¹⁶⁹ 

Crucial for the computation of vibrational features is the ability to determine the electronic polarizability tensor. In recent work, DeePMD was combined with an additional deep neural network to learn the environmental dependence of the polarizability tensor.¹⁷⁰ As demonstrated using the SCAN (Strongly-Constrained and Appropriately-Normed) functional as a reference, this approach yields accurate Raman spectra of liquid water. Furthermore, DeePMD has been integrated with a deep neural network trained to predict Wannier centers based on local environments. This approach allowed us to compute infrared spectra^171,172 and to determine the temperature dependence of the dielectric constant.¹⁷³ The methods can also be extended to account for quadrupole moments.¹⁷⁴ Recently, the description of tensorial properties, such as the polarizability tensor, has been fitted to MD simulations a posteriori using equivariant neural networks to describe infrared spectra.¹⁷⁵ Moreover, a combination of HDNNP and GPR has allowed us to model the hyper-Raman spectra of water, which helped us understand the differences in the OH stretch mode between infrared and Raman spectra.¹⁷⁶ 

The accuracy and computational efficiency of MLPs have made it possible to accurately determine the thermodynamic properties of water, including its phase diagram. For instance, the thermodynamic properties of water have been investigated with DeePMD based on SCAN^177,178 and with a HDNNP based on revPBE0-D3.^159,179 The phase behavior of water under extreme conditions expected in a planetary environment was also studied employing a HDNNP.¹⁸⁰ Other aspects addressed with MLPs include the study of heat transport^181,182 and the viscosity.¹⁸³ More recently, the thermodynamics of water has been investigated with a neuroevolution potential.¹⁸⁴ 

Studying the mechanism and kinetics of phase transitions is computationally very demanding and thus completely out of reach for ab initio simulations. MLPs, however, can be used to simulate systems of millions of water molecules with ab initio accuracy¹⁸⁵ such that the simulation of phase transitions is now possible. One example is a recent investigation of the homogeneous nucleation of ice in supercooled water studied using DeePMD, in combination with the seeding methods,^186,187 in a system of hundreds of thousands of water molecules.¹⁸⁸ These calculations yielded nucleation rates consistent with experimental measurements. Very recently, advanced sampling techniques covering 36 µs of total simulation time have been used to probe the atomic structure of the critical nucleus.¹⁸⁹ In addition, the liquid–liquid transition in supercooled water has been investigated. An initial attempt using DeePMD based on the SCAN functional found indications of this transition through anomalies in thermodynamic response functions.¹⁹⁰ Two years later, the existence of this transition was conclusively demonstrated¹⁹¹ and its relation with the melting curves of ice polymorphs has been investigated.¹⁹² Finally, building on previous studies, the transition between ice Ih and its proton-ordered counterpart ice XI, mediated by ionic defects, has been studied based on the DeePMD model.¹⁹³ 

Since modern MLPs can capture reactions, important processes such as proton transfer and autoionization are accessible. In this regard, HDNNPs have allowed us to describe the transport of hydronium and hydroxide ions, including nuclear quantum effects¹⁹⁴ and the free energetics and mechanics of water dissociation,¹⁹⁵ allowing us to compute the equilibrium pK_w of water.¹⁹⁶ 

The application of MLPs has been facilitated by careful benchmarking and transferability studies and the development of new ML-based methods and workflows. Some of these ML methods have been used beyond the fitting of potential energy surfaces as is the case for ML classifiers of phases^197,198 and dynamical processes.¹⁹⁹ One particular interesting finding regarding the transferability of MLPs is that liquid structures already contain the relevant information required to reproduce ice phases,²⁰⁰ including even ice–water interfaces. In particular, this was observed in studies of homogeneous nucleation, in which an empirical potential was benchmarked against its MLP representation containing only liquid structures.²⁰¹ The performance of different density functionals [Perdew–Burke–Ernzerhof (PBE), SCAN, vdW-cx, and optB88-vdW] for modeling water and ice has been compared.²⁰² Moreover, HDNNP and GPR trained on the same dataset have been shown to be equivalent when compared over different thermodynamic properties of liquid water although HDNNPs seem to be more demanding in terms of the required training data²⁰³ but are computationally more efficient. In fact, the role of the training data has been the focus in other studies.²⁰⁴ Furthermore, graph neural networks (GNNs), which do not require predefined structural descriptors, have been applied to accelerate molecular dynamics simulations.²⁰⁵ The selection of descriptors has also been automatized for HDNNPs,^206,207 whereas GPR-based potentials have been employed in on-the-fly learning workflows.^204,208

Empirical potentials have also benefited from the development of MLPs. Coarse-grained MLPs emerged,^{165,209–213} including an approach based on equivariant neural networks,²¹⁴ and empirical force fields were parameterized using ML algorithms.^215–218 Moreover, the addition of polarization to empirical force fields has been revisited, including charge transfer.²¹⁹ DeePMD has been used to fit an accurate but costly many-body potential,²²⁰ reducing its computational cost by one order of magnitude.²²¹ Moreover, a GNN has been applied to estimate Bayesian uncertainty in molecular dynamics simulations based on an empirical potential.²²² 

Further progress has come again from water clusters. HDNNPs, PIP-based potentials, and GAPs have been shown to be equivalent in representing many-body interactions in water clusters,²²³ which can be employed in the construction of improved water potentials for bulk water. In fact, recent advances suggest that reference data obtained exclusively for water clusters could be sufficient to train accurate MLPs even for the bulk liquid phase,^224–228 including results from Gaussian-moment neural networks (GMNNs).²²⁸ Recently, even gold-standard CCSD(T)-level accuracy for bulk water potentials has been reached by training to large clusters or periodic structures.^121,229,230

An important current topic of research is the inclusion of long-range interactions,⁸⁶ which are not explicitly considered in many MLPs. This problem has been addressed by introducing non-local representations of the system remapped as local and equivariant feature vectors, capturing non-local and non-additive effects.²³¹ Another approach to treat this issue is to learn the long-range response with a self-consistent field neural network, which has been shown to produce correct long-range polarization correlations in liquid water, as well as the correct response of liquid water to external electrostatic fields.^232,234

Architectures like equivariant neural networks have been combined with empirical electrostatics and dispersion.²³⁵ Such models are highly accurate in learning reference datasets,^{80,81,235–238} and their adoption is growing rapidly.

Neural network potentials have also been used to investigate the structure, thermodynamics, and spectroscopic properties of the liquid/vapor interface. As the local environments close to the interface are highly anisotropic and thus very different from the bulk, it is important that the training set explicitly includes data for interface configurations.²³⁹ The structure of such configurations has been analyzed in detail using SOAP descriptors and local order parameters.²⁴⁰ Investigating the structure of the interface reveals the prevalence of orientations with the dipole moment roughly parallel to the surface with one OH bond pointing out of it,²³⁹ corroborating insights gained from sum frequency generation (SFG) measurements.²⁴¹ By using the surface-sensitive velocity autocorrelation function,²⁴² such SFG spectra of the liquid/vapor interface were calculated from path integral molecular dynamics based on a HDNNP trained at the revPBE0+D3 level.²⁴³ Recently, SFG spectra have been computed fully from first principles using a HDNNP combined with GPR.²⁴⁴ 

In another study, it was found that a DeePMD potential relying only on local atomic energies can be applied to the liquid/vapor interface.²⁴⁵ However, the explicit inclusion of long-range interactions was shown to be beneficial, confirming the results of previous studies carried out for empirical potentials.²⁴⁶ The effect of long-range interactions was tested for a water molecule moving away from the liquid/vapor interface using an extension of DeePMD, including long-range electrostatics.⁹⁷ The case of curved liquid/vapor interfaces has been addressed as well. For instance, DeePMD has been employed to investigate the formation of bubbles in metastable water.²⁴⁷ Furthermore, it was shown that the free energy of water dissociation at the liquid/vapor interface of droplets and films deviates from the bulk, leading to an enrichment of hydronium cations at the interface and a depletion of hydroxide anions.²⁴⁸ 

Beyond pure water, MLPs have been used in numerous simulations of electrolyte solutions.²⁴⁹ Already in 1998, a feed-forward neural network was employed to represent the three-body interaction energies in H₂O–$Al3+$–H₂O clusters with the aim to improve the force field description of $Al3+$ ions dissolved in bulk water.⁵⁶ This work represents an important milestone in the incorporation of permutation symmetry in structural descriptors. Later, HDNNPs have enabled the construction of full-dimensional DFT-quality PESs for aqueous NaOH solutions over the entire solubility range.²⁵⁰ In this work, it has been found that as the NaOH concentration increases, the primary mechanism for proton transfer shifts from being acceptor-driven, influenced by the pre-solvation of hydroxide ions, to donor-driven, controlled by the pre-solvation of water molecules. In addition, with increasing concentration, octahedral coordination geometries become less favored, in contrast to trigonal prism geometries.²⁵¹ A novel water exchange mechanism has been identified around Na⁺(aq) ions in basic (high pH) solutions.²⁵² Studies comparing classical and ring-polymer molecular dynamics based on the HDNNP revealed that nuclear quantum effects significantly reduce proton transfer barriers, thus increasing proton transfer rates. This leads to an enhanced diffusion coefficient, especially for OH⁻, and a shorter mean residence time of molecules in the first hydration shell around Na⁺ at high NaOH concentrations.²⁵³ Moreover, elevated temperatures in concentrated NaOH solutions amplify both the contributions of proton transfer to ionic conductivity and deviations from the Nernst–Einstein relation.²⁵⁴ Further applications of HDNNPs include investigations of fluoride and sulfate ions in the solution.¹⁰⁷ Employing similar methodologies and training on revPBE + D3 data, the dissolution mechanisms of NaCl in water have also been addressed.²⁵⁵ Another example is the use of HDNNPs to study zinc ion hydration in water,²⁵⁶ with molecular dynamics simulations matching both the experimentally observed zinc–water radial distribution function and the x-ray absorption near edge structure spectrum. Moreover, HDNNP-based studies reveal the impact of surface stratification on the interfacial water structure in electrolyte solutions.²⁵⁷ Equivariant neural network potentials have also been employed to study various electrolyte solutions,²⁵⁸ including aqueous lithium chloride²⁵⁹ and aqueous sodium chloride.²⁶⁰ A genetic algorithm has also been utilized to study hydrated zinc(II) ion clusters.²⁶¹ In addition, microhydrated sodium ions with a few water molecules have been studied for both the potential energy and the dipole moment employing PIPs.²⁶² 

DeePMD potentials have been used to study sodium chloride, potassium chloride, and sodium bromide at various concentrations.²⁶³ These studies revealed that the structural changes due to the ions are confined to the immediate vicinity of the ions, where they disrupt the network of hydrogen bonds. Beyond these regions, the distribution of oxygen atoms relative to one another remains largely unchanged compared to pure water. In a related study, the dielectric permittivity of sodium chloride solutions has also been investigated.²⁶⁴ Using DeePMD potentials, the uptake of N₂O₅ into aqueous aerosols has been examined, a process that is challenging to study experimentally due to the fast reaction kinetics of N₂O₅.²⁶⁵ Furthermore, the diffusivity of water in aqueous cesium iodide and sodium chloride solutions has been examined using a DeePMD framework trained on DFT data using the revPBE-D3 functional.²⁶⁶ Such simulations addressing the characteristic behavior of different ions are not readily accessible through traditional force field-based molecular dynamics simulations due to less ion-specific description of ion–water interactions.

Solid–liquid interfaces are of high interest for catalysis and electrochemistry. Due to the very different bonding in liquid water and in crystalline surfaces, such as metals or oxides, constructing unified atomistic potentials that can describe all subsystems of solid–liquid interfaces with balanced high accuracy presents a substantial challenge for empirical potentials. Moreover, in many cases, water is not only in contact with the surface but can also dissociate and recombine at a much higher rate than in the bulk liquid. Consequently, the use of reactive potentials, which can describe the making and breaking of bonds, is mandatory. MLPs are ideally suited for this purpose.

In 2014, a HDNNP for a thin water layer on top of 55-atom CuAu alloy clusters with varying stoichiometries and a slab model were reported to study the effect of water on the stability of different interface compositions by Monte Carlo simulations.²⁶⁷ While this work exhibited a still rather large error of about 12 meV/atom, the accuracy of HDNNPs for solid–liquid interfaces has significantly improved in the following years. For instance, in studies of water at various surfaces of copper^268,269 and zinc oxide^270,271 energy, RMSEs of less than 1 meV/atom could be reached. Detailed convergence tests with respect to the required system size were carried out,²⁶⁸ showing that the diameter of the liquid phase between the slab surfaces needed to decouple the two surfaces by bulk-like water is at least about 35–40 Å.²⁶⁸ Such a size is beyond reach in ab initio molecular dynamics but necessary to ensure that the central water molecules have bulk-like environments in their local vicinity. This consideration is crucial to ensure that these molecules do not experience any significant influence from the altered water structure near the interfaces or from the surfaces directly.

While it has been found using HDNNPs that water does not spontaneously react with defect-free surfaces of certain metals, such as copper, on nanosecond time scales,²⁶⁸ fast dissociation and recombination processes leading to the formation of surface hydroxides have been observed on zinc oxide surfaces.^270,271 These processes are often governed by the surrounding hydrogen bond networks, significantly influencing the free energy barriers of proton transfer processes.²⁷⁰ This phenomenon has been confirmed for water at TiO₂ surfaces using DeePMD potentials,^272,273 which have also been employed to explore the impact of slab thickness.²⁷⁴ Furthermore, the use of HDNNPs for computing anharmonic frequencies has been suggested as a method to elucidate the role of hydrogen bonds in surface processes.²⁷⁵ 

Depending on the specific surface geometry, proton transfer events can lead to various topologies of proton transport networks along the surface, which can be either one-dimensional or two-dimensional. This has been demonstrated for several surfaces of zinc oxide²⁷⁶ and the lithium intercalation compound LiMn₂O₄²⁷⁷ using HDNNPs. Surface defects, often stabilized by solvation compared to the vacuum interface, have also been a subject of study. The mobility of adatoms has been found to significantly vary across different low-index surfaces of copper.²⁶⁹ Investigations into the defective Zr₇O₈N₄/H₂O and pristine ZrO₂/H₂O interfaces using neural network potentials²⁷⁸ revealed a bilayer water structure for Zr₇O₈N₄ and a monolayer structure for ZrO₂. Oxygen vacancies on the Zr₇O₈N₄ surface have been suggested as active sites for the oxygen reduction reaction. Furthermore, neural networks have been used to identify different oxidation states of transition metal ions at oxide-water–interfaces,²⁷⁷ which enables the characterization of electronic structures relevant for catalytic applications.

Due to its importance in catalysis, the TiO₂–water interface has so far been the most intensely studied interface using MLPs.^{107,272–274,279–282} Investigations into the water coverage on the anatase (101) TiO₂ surface using the DeePMD potential²⁷⁹ have shown that higher water coverage prompts significant reorganization of the water monolayer at O_2c sites, leading to the formation of a two-dimensional hydrogen bond network with closely linked pairs of water molecules on neighboring TiO_5c and O_2c sites. Other DeePMD-based studies have examined the impact of water dissociation on thermal transport at the TiO₂–water interface.²⁸⁰ While previous research on TiO₂–water interfaces mainly focused on the anatase (101) and rutile (110) TiO₂ surfaces, recent MLP-based studies²⁸¹ have explored seven different TiO₂ surfaces using three distinct functionals: SCAN, PBE, and optB88-vdW. These studies found that water dissociation is more likely on the anatase (100), anatase (110), rutile (001), and rutile (011) surfaces, while molecular adsorption is the primary process on the anatase (101) and rutile (100) surfaces. Moreover, simulations for rutile (110) showed that the slab thickness significantly influences the results, with thicker slabs favoring molecular adsorption. DeePMD has also been used to study amorphous TiO₂ (a-TiO₂) to compare its behavior with well-studied crystalline TiO₂ at aqueous interfaces.²⁸² These studies demonstrated that water molecules on the a-TiO₂ surface do not exhibit the distinct layering typical of the aqueous interface of crystalline TiO₂. This difference results in an approximately tenfold increase in water diffusion speed at the interface.

Other cases of employing MLPs to study solid–water interfaces include the use of HDNNPs for the Pt(111)–water interface to investigate the interaction between water and hydroxylated metal surfaces^283,284 and for the hematite–water interface,²⁸⁵ revealing solvation dynamics at various time scales. DeePMD has been utilized for studying the TiS₂/water interface²⁸⁶ to examine the influence of TiS₂ surface termination on the structure of interfacial water. Moreover, DeePMD has been applied to the IrO₂–water interface, exploring the hydration structure, proton transfer mechanisms, and acid–base characteristics,²⁸⁷ as well as to the GaP(110)–water interface,²⁸⁸ which has been shown to require about 12 ns to reach equilibrium, a duration not achievable with traditional AIMD simulations. DeePMD has also been used for the construction of a potential aimed at studying ice nucleation at the microcline feldspar surface using the SCAN functional²⁸⁹ and for investigating the impact of water dissociation on thermal transport at the Cu–water interface.²⁹⁰ Apart from HDNNPs and DeePMD, an equivariant graph neural network has been employed to study the oxygen reduction reaction at the Au(100)–water interface,²⁹¹ and on-the-fly learning kernel-based regression has been applied to investigate water adsorption on MgO and Fe₃O₄ surfaces, including surface reconstructions.²⁹² 

Confined water, which exhibits properties notably differing from bulk water, has also been studied using MLPs. Some examples include HDNNPs for water confined between two-dimensional boron nitride sheets²⁹³ and MoS₂¹⁰⁷ and water between graphite layers using DeePMD^294,295 and committees of HDNNPs;²⁹⁶ the latter method was also applied to confinement within a graphene-like material.²⁹⁷ Furthermore, MD simulations have been used to investigate water in single-walled carbon and boron nitride nanotubes,^107,298 finding a fivefold reduction in friction in carbon tubes compared to boron nitride, attributed to strong hydrogen–nitrogen interactions.²⁹⁸ Ion concentration profiles under nanoconfinement²⁹⁹ have also been studied using neural network potentials, focusing on the effects of channel widths, ion molarity, and ion types.

Apart from studies of pure water, electrolytes, and solid–liquid interfaces, which have been in the focus of MLP-based atomistic simulations for several years, the use of MLPs for aqueous systems has increasingly diversified and now covers essentially all fields of simulations involving water. An exhaustive coverage of the related literature is beyond the scope of this Perspective, and in this section, we just point the interested readers to several typical applications of MLPs in these rapidly growing fields.

A prominent use of atomistic simulations is to study chemical reactions of organic molecules in solution. Examples for the application of MLPs are corrections to quantum mechanical/molecular mechanical (QM/MM) simulations of S_N2 reactions in water,³⁰⁰ the solvation of protein fragments,⁹³ the decomposition of urea in water,³⁰¹ the computation of free energy profiles of reactions of organic molecules,³⁰² and enzyme reactions.³⁰³ Further studies include the quantum dynamics of an electron solvated in water³⁰⁴ and the excited state of CH₃NNCH₃ surrounded by several water molecules.³⁰⁵ 

Still, the all-atom description of chemical processes in solution can be demanding, and, consequently, simplified MLPs have been proposed as well. For example, the solvated alanine dipeptide^306,307 and the folding/unfolding of chignolin^306,308 have been studied using the coarse-grained CGNet potential. Moreover, a DeepPot-SE model has been used to describe molecules under the influence of an implicit solvent.³⁰⁹ 

Additional applications of MLPs for aqueous systems include the study of diffusion in hydrogen hydrates,³¹⁰ the determination of vibrational frequency shifts in formic acid C=O stretching and C=N stretching of MeCN in water,³¹¹ retinoic acid in water,³¹² graph-convolutional neural networks for benchmarking sets of solutes and chemical reactions in water,³¹³ hydration dynamics and IR spectroscopy of 4-fluorophenol,³¹⁴ zinc protein studies,³¹⁵ conformational shifts of stacked heteroaromatics,³¹⁶ and solvation free energy prediction of organic molecules in redox flow batteries.³¹⁷ 

In recent years, MLPs have reached a high level of maturity, and, currently, a transition from proof-of-concept and benchmark studies to practical simulations of a wide range of complex systems is taking place. Therefore, it can be anticipated that MLPs will allow us to overcome the limitations of conventional methods, such as empirical potentials in terms of accuracy and ab initio molecular dynamics in terms of efficiency, paving the way for simulations of extended systems with unprecedented accuracy. MLPs have demonstrated this capability already for a broad spectrum of aqueous systems, ranging from neutral and protonated water clusters to bulk liquid water and ice, liquid/vapor interfaces and from electrolyte solutions to complex solid-water interfaces (see Fig. 1). In all these studies, MLPs have enabled simulations with first-principles accuracy that previously have been prohibitively demanding, as evidenced by a rapidly growing number of publications in the field shown in Fig. 2.

MLPs bridge the gap between two traditional approaches in atomistic simulation, ab initio and force-field-based MD simulations, offering advantages over both. MLPs excel by enabling significantly longer length and time scales compared to AIMD simulations, thanks to their computational efficiency. This extension is crucial for accurately computing properties, such as equations of state across a broad range of parameters, kinetics of phase transitions, and interface mobility, and to study supercooled water as well as solvation dynamics on long time scales. In comparison with classical force fields, MLPs can potentially match the accuracy of the underlying electronic structure calculations, which usually surpass the accuracy of traditional force fields. Furthermore, MLPs are reactive since they do not rely on predefined chemical bonds, unlike most traditional force fields. This reactivity is crucial for correctly modeling systems with complex interactions, such as solid–liquid interfaces and aqueous solutions, where water dissociation, recombination, and proton transfer occur. It is also essential for accurately representing the entire range of acid/base chemistry, which fundamentally relies on proton transfer processes.

Over the past two decades, progress in MLPs for aqueous systems has focused on different frontiers. While right from the start the highly flexible functional form of ML algorithms has enabled a numerically very accurate representation of the electronic structure reference data, a severe challenge in the early years has been the very limited number of degrees of freedom that could be considered. Only the development of modern descriptors for the atomic environments allowed us to extend MLPs to condensed systems such as liquid water with all their associated degrees of freedom. Recently, message passing neural networks⁷⁵ have become a promising alternative to the use of predefined descriptors, opening many new exciting possibilities for the construction of MLPs with higher accuracy, based on less data, in particular if equivariant features are used. Another frontier that has increasingly received attention in recent years is the incorporation of physical concepts into hitherto purely mathematical machine learning potentials, with many developments specifically aiming for improved descriptions of long-range electrostatic interactions, van der Waals forces, long-range charge transfer,⁵² and electron densities.³¹⁸ This inclusion of physically meaningful terms, not only in the total energy expression but also in the form of novel descriptors,^231,319 will further increase the accuracy and transferability of the potentials.

To effectively simulate chemical processes in the aqueous phase, the underlying potential function must rely on accurate electronic structure calculations, since the chosen reference method represents a natural limit for the accuracy of MLPs. While coupled cluster accuracy has already been achieved in MLPs for liquid water,^121,229,230 which would have been unthinkable with conventional empirical potentials, reaching this gold-standard for more complex systems, such as solid–liquid interfaces, is very challenging. Therefore, DFT will likely remain the dominant method for the reference electronic structure calculations of many systems in the foreseeable future, in particular for those systems involving solids. Although the generalized gradient approximation (GGA) has been the most commonly employed functional in the study of aqueous systems as they offer a good compromise between computational cost and accuracy, more and more powerful computing resources increasingly enable the use of more advanced and computationally expensive functionals, such as meta-GGA (particularly SCAN) and hybrid functionals, in the investigation of several aqueous systems, and even this level of accuracy has remained essentially inaccessible by on-the-fly ab initio MD to date. Consequently, complex aqueous systems can now be investigated with a previously unattainable level of accuracy enabling predictive simulations. Moreover, MLPs can be employed to evaluate the accuracy of the underlying level of theory with respect to experimental values in a broad range of scenarios. This extends beyond the traditional comparison at standard temperature and pressure, allowing for a comprehensive evaluation that may provide insightful guidance for the development of theoretical methods.

The functional flexibility is a key property of MLPs, but it is a double-edged sword: on the positive side, this flexibility enables the accurate approximation of the PES based on the reference data. On the negative side, however, such flexibility severely limits the extrapolation capabilities of MLPs to chemical environments not adequately sampled during the training process. In fact, extrapolation to unfamiliar environments can lead to unphysical structures and completely wrong simulation results. Therefore, the construction and validation of MLPs have to be done with great care to ensure that all relevant local environments are included in the training set. Moreover, it is much more challenging than in case of simpler empirical potential to provide “boxed” MLPs for general usage, since not only the underlying parameters but also information about the range of validity is crucial information for successful applications. For instance, when studying solid–water interfaces, it becomes necessary to train the MLP not only on the bulk material and bulk water separately but also on systems that include all relevant interface configurations. Thus, increasing the complexity of the system also increases the number of local environments that must be included in the training data. Although AIMD simulations are sometimes used to generate initial training sets, they are usually insufficient since they often fail to capture the less frequently visited structures. Fortunately, this challenge has been largely overcome in recent years by incorporating active learning for the generation of the reference data. An alternative approach, recently emerging, focuses on establishing foundation models to be used as a starting point for MLP development.²³⁷ This might improve the transferability and contribute to the development of more accurate MLPs.

While, in this Perspective, we have focused on the construction and application of accurate and efficient interatomic potentials, machine learning approaches can be useful in several other ways for the atomistic simulation of aqueous systems, and, more generally, of materials and biomolecular systems.⁴⁷ For instance, neural networks have been employed to classify local structures in liquid water and various forms of ice.^197,198 The accurate identification of molecular structures with high spatial resolution is important, for instance, in the study of crystallization, melting, crystal growth, and the formation and migration of defect. Another challenge in molecular simulations lies in reducing or even completely removing correlations in the statistical sampling of configuration space as it occurs in sequential sampling methods, such as molecular dynamics and Markov chain Monte Carlo simulations. Here, normalizing flows^320,321 or other generative methods can play an important role for sampling the equilibrium distribution and also for creating new stable crystal structures.³²² Moreover, machine learning approaches have been suggested as a way to discover reaction coordinates and enhance the sampling of rare transitions occurring in complex molecular systems.^323,324 There is little doubt that in the years to come further, new ML and AI tools will be applied to the computational investigation of matter at the atomistic level, creating new opportunities for studying complex aqueous systems.

Machine learning and, more generally, artificial intelligence are currently revolutionizing our way to do science and are providing new opportunities not achievable with traditional approaches. In the field of computational materials science, the advent of accurate, flexible, and efficient MLPs has dramatically increased the time and length scales accessible by atomistic simulations of materials with ab initio accuracy. In this Perspective, we have provided an overview of the key concepts of such machine learned potentials with a focus on their application to water and aqueous systems. The broad spectrum of systems successfully explored with these techniques, including water clusters, bulk liquid water and ice, the liquid/vapor interface, electrolyte solutions, and solid–liquid interfaces, underscores the flexibility, efficiency, and the high level of maturity they have reached since recent years.

Since machined learned potentials accurately reproduce the underlying reference data obtained with electronic structure methods but in applications require a much lower computational effort, they now make it possible to compute complex materials properties, such as phase diagrams. In this way, they do not only allow us to gain new insights into a variety of systems, but they also provide a way to truly test the theoretical description underlying the reference data and to reveal their possible limitations. This facilitates the generation of high quality reference data in the future as a basis for truly predictive ML-based computer simulations of complex materials.

In summary, modern MLPs have created new opportunities for the investigation of aqueous systems that would have been unimaginable with conventional methods for the foreseeable future. Carefully trained and validated MLPs can be employed to study complex reactive aqueous systems accurately across large time and length scales, without imposing ad hoc empirical constraints. While predicting the future of this rapidly evolving field is challenging, the remarkable progress made to date suggests that we can expect exciting new developments and some surprising breakthroughs in the years to come.

J.B. and A.O. are grateful for the funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Grant No. TRR/CRC 247 (A10, Project No. 388390466) and under Germany’s Excellence Strategy—Grant No. EXC 2033 RESOLV (Project No. 390677874). C.D. and P.M.D.H. acknowledge the funding from the Austrian Science Foundation FWF through the Projects Doctoral College Advanced Functional Materials (Grant No. DOC 85-N) and the SFB TACO (Grant No. F-81).

The authors have no conflicts to disclose.

All authors contributed equally to this work.

Amir Omranpour: Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Pablo Montero de Hijes: Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Jörg Behler: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Christoph Dellago: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.