Aqueous solution chemistry in silico and the role of data-driven approaches

Liquid water is one of the key ingredients for life.¹ It forms a central lubricant for biological materials and is perhaps the most ubiquitous solvent in physical, chemical, engineering, and technological applications.^2–4 Many of the unique properties of water arise from its hydrogen-bond network^5–8 and how it changes in different thermodynamic conditions. Dissecting the microscopic structure of water and aqueous solutions in terms of both the static and dynamical properties of hydrogen bonding has been the subject of numerous experimental and theoretical studies and rather lively controversies.^9–11 

Computer simulations of varying levels of complexity have played an important role in providing molecular-level insights into the structural, dynamical, and electronic properties of both bulk water and solutes in aqueous solutions.^11–13 Over the last five decades, the typical models that are used to simulate aqueous systems can be broadly separated into two categories namely, molecular mechanics (MM)¹⁴ empirical potentials and first-principles ab initio molecular dynamics (AIMD) approaches.¹⁵ While MM-based water models provide a powerful way to explore, for example, the phase diagram of water^6,9 and perform simulations of complex biological systems,^16–18 they typically do not allow for chemistry to occur. Although being more computationally prohibitive, electronic structure-based AIMD simulations overcome this limitation allowing for modeling chemical reactions where bond breaking or formation occurs. One of the most popular electronic structure approaches for modeling the properties of water has been Density Functional Theory (DFT),^12,13,19 although over the last decade more advanced approaches have also been employed, using, for example, Quantum Monte-Carlo²⁰ and quantum chemistry techniques such as Møller–Plesset perturbation theory (MP2).²¹ 

In both the MM and AIMD approaches to study aqueous systems, there are several challenges that have emerged if one is interested in producing meaningful simulations that can be interpreted and compared with experiments. The first is the quality of the electronic structure theory which ultimately controls the underlying potential energy surface (PES) associated with the hydrogen bonds needed to reproduce structural, dynamical, and spectroscopic properties of water.^19,22,23 Second, there is the problem of sampling arising from the fact that it takes a long time for the system to hop over free-energy barriers.^24,25 This is particularly true for modeling chemical reactions such as water ionization in different environments^26–28 with AIMD, as well as exploring the complex phase diagram of water, for example, under supercooled conditions.²⁹ Finally, in many particle systems involving the coupling of both the solute and its aqueous environment, identifying the relevant and important degrees of freedom (commonly referred to as order parameters, reaction coordinates or collective variables) along which physical, chemical, or biological processes occur, is far from trivial.^30,31

Herein, we will provide an overview of the most recent advances in the field aimed at overcoming these challenges. Figure 1 provides a schematic outline of our review. In the first part, we begin with a brief historical perspective on the challenges in modeling liquid water from first principles electronic structure simulations. Thereafter, we discuss the advances made in the development of reactive, many-body, and machine-learning (ML) based potentials which are opening new turf in studying the physical chemistry of liquid water [Figs. 1(a) and 1(b)]. We focus specifically on highlighting important advances that have been made in the applications of these techniques using them to also underscore the challenges in modeling aqueous chemistry in the bulk, at interfaces, and under confinement. Furthermore, a critical aspect associated with understanding aqueous chemistry is the determination and exploration of relevant order parameters or reaction coordinates. Due to the collective nature associated with processes involving chemistry in water, identifying these coordinates requires going beyond chemical imagination; we highlight in this context the importance of data-driven approaches [Fig. 1(c)]. Finally, we conclude with some perspectives on the future of data-driven approaches in empowering conceptual advances in the chemical physics and physical chemistry of aqueous systems.

The present section is essentially divided into three parts, as illustrated in Fig. 1. We begin by providing an overview of the main approaches that have been used to simulate aqueous chemistry focusing on DFT, many-body and reactive potentials. As indicated in the Introduction, our focus is to highlight key applications related to aqueous chemistry that underscore the challenges. This subsequently serves as motivation for the second part where we focus on recent developments in the use of machine-learning potentials for modeling water and its constituent ions in different contexts. In the last part of this section, we discuss advances and challenges in data-mining for water chemistry.

DFT has become one of the most popular electronic structure methods for simulating complex chemical and biological systems from first principles^32,33 due to its favorable compromise between computational efficiency and accuracy. The standard way in which one can couple DFT with finite temperature molecular dynamics is using Born–Oppenheimer (BO) or Car–Parrinello (CP) MD. For more details on the theory of DFT, BOMD, and CPMD, the interested reader is referred to specific books and other reviews on this subject and references therein.^13,34,35

DFT-based ab initio molecular dynamics simulations (AIMD) of bulk liquid water and its constituent ions, the proton and hydroxide, offered a first glimpse into the coupling of nuclear and electronic degrees of freedom in water.^36–39 In particular, such techniques have opened up a microscopic window into the celebrated Grotthuss mechanism⁴⁰ involving the inter-conversion of covalent and hydrogen bonds in the water network. These simulations have since guided the interpretation of many IR and Raman-based spectroscopy experiments of acidic and basic water.^41–44 

Although DFT is an exact theory, within the Kohn–Sham formalism, the functional form for the exchange-correlation (XC) energy is not known and various approximations need to be made to treat this term. As a result, this poses limits on the accuracy that DFT can achieve for the prediction of electronic structure properties. Numerous studies over the last decades have highlighted how the choice of the XC functional can, in fact, hugely impact the description of the structural and dynamical properties of water. For a detailed review of the effect of the DFT XC functional on water in different conditions, the reader is referred to a perspective by Michaelides and coworkers.¹⁹ Panel (a) in Fig. 2 shows the radial distribution functions (RDF) obtained for two different popular functionals with and without the inclusion of dispersion interaction corrections (e.g., Grimme's semiempirical D3⁴⁵) compared to experimental measurements.^46,47 In summary, it can be seen that accounting for dispersion interactions plays an important role in correcting for approximations made in standard generalized-gradient approximation (GGA) functionals, an aspect that has been reinforced in numerous studies.^36,48–50 Specifically, standard GGA functionals yield a much more glassy liquid. In addition to the inclusion of dispersion corrections,^19,51–53 the inclusion of exact Hartree–Fock exchange through the use of hybrid functionals has also been shown to give more accurate estimates of water polarizability.⁵⁴ However, despite the generally favorable scaling of DFT and especially in cases where sophisticated functionals and extra dispersion corrections are employed, the computational cost of the simulations poses severe limits on both the size of the system investigated [relatively small box sizes (<2 nm)] and the length of the dynamics [short simulation times (∼100 ps)].

In addition to the challenges of dealing with the quality of the electronic structure, another issue in simulating the physical and chemical behavior of water is the high zero-point energy (ZPE) of the O–H covalent bonds, which is almost about an order of magnitude larger than thermal energy at room temperature. This is reflected in the O–H pair correlation function of water, which tends to be overlocalized in simulations where the nuclei are treated classically. Car and coworkers were the first to elucidate the importance of nuclear quantum effects (NQEs) on the structural properties of the hydrogen bonds in water by using path-integral molecular dynamics (PIMD) together with CPMD simulations.⁵⁵ Examining the O–H RDF, one observes that in the PIMD simulations, the proton is much more delocalized, leading to a broader first peak consistent with the experiments. In this context, path-integral approaches coupled with BOMD or CPMD simulations have become somewhat routine, especially in combination with generalized Langevin-based thermostats.^56–58 These and other types of simulations that include NQEs have shown that these effects are important for understanding the structural, dynamical, and even electronic properties of bulk water as well as chemical and biologically relevant solutes.^59–67 

The broader O–H bond length distribution found in PIMD was associated with enhanced transient autoionization events where charged-pairs of waters' constituent ions apparently form due to extreme events involving protons delocalizing along hydrogen bonds.⁶⁸ The extent of these proton fluctuations are, however, very sensitive to the quality in the underlying electronic structure. Specifically, the details of the DFT functional along with dispersion corrections compete with NQEs in highly non-trivial ways, for example, in affecting spectroscopic properties of hydrogen bonds. Figure 2(b) nicely demonstrates this effect by illustrating how the high frequency modes of the O–H vibrations change as a function of combining hybrid DFT functionals with nuclear quantum dynamics. Marselek and Markland show that standard GGA functionals, even when including dispersion corrections, exaggerate the fraction of these transient autoionization events leading to larger red shifts in the O–H vibrational stretch frequencies.⁶⁹ Upon using hybrid functionals together with dispersion corrections, the inclusion of quantum effects almost perfectly reproduces the experimental spectra.

The transient autoionization events previously discussed are, of course, the initial seeds for water dissociation that leads to the creation of hydronium (H₃O)⁺ and hydroxide (OH) $−$ ions. The preceding issues regarding the choice of the quality in the electronic structure description within the framework of DFT have been shown to play a critical role in affecting structural and dynamical properties of these ionic topological defects. In Fig. 2(c), Voth and coworkers perform a systematic analysis on the factors affecting the diffusion constant of the excess proton in water.⁷⁰ The reported values of the diffusion constant are extremely sensitive to the inclusion of dispersion corrections, the choice of the basis set, the density functional employed and, finally, the initial conditions. Furthermore, Wu and Car recently studied the role of including exact exchange on the diffusion mechanism of the excess proton and hydroxide in liquid water.⁷¹ Previous studies from some of us had suggested that the Grotthuss mechanism involves concerted proton hopping events for both the proton and hydroxide ion.⁷² Specifically, it was found that protons and proton-holes (hydroxide ions) diffuse along water wires involving concerted double jumps. In Fig. 2(d), Car and coworkers demonstrate that while the use of GGA functionals display consistent results for the excess proton compared to the hybrids, the story is much more complicated for the hydroxide ion. In particular, they show that the migration of hydroxide is dominated by a stepwise mechanism rather than concerted hopping. The qualitative difference between the effect of the DFT electronic structure on the proton and hydroxide is rooted in the fact that the use of hybrid functionals along with ab initio dispersion corrections sensitively affects the solvent environment of the hydroxide, favoring a hyper-coordinated vs a three-coordinated solvation structure.

In addition to looking at the structural, thermodynamic, and dynamical properties of water, there have also been several studies examining its electronic properties, for example the bandgap, refractive index, and polarizability—the reader is referred to those references for more details.^74–77 In this respect, improvements over traditionally used GGA density functionals were obtained by employing more recent meta-GGA functionals,⁷⁸ even though the electronic density of states is far from being satisfactory over the whole relevant energy scales. To improve the DFT predictions, in some of these studies,⁷⁶ AIMD simulations are performed in the ground-state in order to obtain a sampling of nuclear geometries from which one then performs very accurate many-body perturbation theory in order to extract quantities such as the electron affinity or absorption spectra. Within these contexts, the sampling protocol for the nuclear geometries can sensitively affect the results. Although historically most electronic-structure schemes were unable to accurately reproduce the experimental electronic energy levels of liquid water, a recently proposed many-body perturbation theory with effective vertex corrections produced photoemission and absorption spectra in excellent agreement with experiments.⁷⁹ 

Although DFT-based MD simulations (discussed above) provide a very practical way to study aqueous solutions and have proven useful in obtaining new insights into mechanisms that cannot be observed experimentally, the quality of the underlying potential energy surfaces corresponding to widely used XC functionals does not reach the chemical accuracy required in several applications.^19,80,81 For example, even though a specific XC functional may work for bulk water at ambient conditions, studies have shown that there is no guarantee of transferability to other regions of the phase diagram.¹⁹ At the same time, due to the computational bottleneck of converging the electronic structure, one is often limited to rather small system sizes and simulation times (on the order of 100 atoms and picoseconds). For example, while the simulations shown in Fig. 2 are certainly commendable, they are limited to timescales on the order of hundreds of picoseconds. Alternatively, one can avoid the explicit calculation of the electronic structure of a system through the use of model potentials, often called molecular mechanics (MM) force fields. This approach is widely used due to its low cost compared to DFT, but is often hindered by limits in accuracy and transferability.

Over the last decade, the Paesani group has made significant strides to fill some of these gaps in DFT and MM approaches through the development of the MB-pol many-body potential for water.^22,82–85 It should be noted that there currently are several other many-body potentials available, such as MB-UCB⁸⁶ and q-AQUA;⁸⁷ however, we choose to focus this section on MB-pol since it has had 10 years of validation since its release.

The idea behind these many-body potentials follows from two major principles: (1) a large difference in the complexity of close-range interactions compared to long-range interactions; (2) the dominance of the first three terms in the many-body expansion (MBE) for aqueous systems. In brief, long-range interactions such as permanent electrostatics, classical polarization, and London dispersion interactions can be accurately accounted for through standard MM-based force field approaches.^{14,23,88–90} At close range, interactions become significantly more complicated, including repulsive interactions due to Pauli exclusion, charge penetration effects between overlapping orbitals of different molecules, and charge transfer effects.^{86,88,91–102} DFT and other electronic structure methods capture the full extent of all long-range and close-range intermolecular interactions (within the accuracy of their underlying approximations). However, obtaining simple analytic functions which provide an accurate representation is extremely challenging in general.

In the context of many-body interactions in water, $E 1 B$ corresponds to the sum of distortion (stretching and bending) energies of individual water molecules, $E 2 B$ corresponds to interactions between pairs of water molecules (e.g., hydrogen bonding), and each $E nB$ term for n > 2 corresponds to (excess) polarization effects (including charge transfer) between n water molecules.

Typical rigid, nonpolarizable MM force fields^90,104,105 such as the SPC^106,107 and TIP^108–110 families of models only include explicit 2B potentials, often in the form of the Lennard-Jones potential describing long-range dispersion and close-range repulsive interactions, and pairwise Coulomb potentials capturing permanent electrostatic interactions. Although classical induction can be added on top of pairwise potentials to include higher-body interactions, they are inadequate for capturing the reorganization of electronic structure that occurs at close-range. A critical improvement to classical potentials can be made by recognizing that the MBE for water–water interactions converges rapidly, with the first three terms in Eq. (1) containing the vast majority of the total binding energy.^111–114 Building on earlier ideas used in the development of, for example, the CC-pol^115–119 and WHBB^120–123 families of models for interactions within gas phase water clusters, the Paesani group combined a long-range potential (consisting of dispersion, permanent electrostatics, and classical induction) with highly accurate models for close-range 2B and 3B interactions fitted to coupled-cluster [CCSD(T)^124,125 and CCSD(T)-F12b^126,127] energies. The resulting model, MB-pol, is able to reproduce thermodynamic and dynamical properties of condensed phase water across much of the phase diagram as discussed below. There is, of course, a computational cost associated with the use of many-body potentials, lying around one order of magnitude above classical polarizable potentials, but several orders of magnitude below AIMD.²² 

For neutral bulk water, MB-pol has proven to be among the most accurate potentials currently available,^{22,23,128,129} and extensions such as MB-nrg potentials^{101,130–137} for modeling ions in solution also display comparable accuracy. Comparisons between MB-pol and CCSD(T) for two- and three-body energies of clusters taken from MD simulations demonstrates the quality of the model potential [Figs. 3(a) and 3(b)]. The RMSD obtained for these two energies are ten times smaller than thermal energy at room temperature. It is worth noting that errors obtained with classical polarizable models such as TTM and AMOEBA have been reported to range between up to 2-to-3 times larger than thermal energy.²³ This feature is manifested even at higher temperature and in the condensed phase of liquid water. Note that while low-order many-body energies can be seen as gas phase properties, evaluating the model accuracy on dimers and trimers extracted from MD simulations allows for validation of condensed phase energetics, given the rapid convergence of the MBE in water. As an illustration of this point, Fig. 3(b) shows the average (AD), average absolute (AAD) and maximum absolute (MAD) deviations of the total binding energy per water molecule calculated with a number of classical potentials relative to those calculated with MB-pol for configurations obtained from a molecular dynamics simulation of MB-pol at 300 K.²³ 

Figure 3(c) shows the pair-correlation function of MB-pol liquid water across a range of temperatures compared to experiments. Although there is a slight tendency for MB-pol to overestimate the height of the first peak, very likely due to the fact that these simulations were conducted without nuclear quantum effects, there is a consistency in the distributions from the supercooled regime to below the boiling temperature.²² Another key experimental probe for the structure of liquid water is x-ray absorption spectroscopy (XAS).^138,139 In Ref. 140, it was shown that Bethe–Salpeter electronic structure calculations on top of approximate quantum dynamics with the MB-pol potential correctly reproduces the experimentally measured XAS, further validating the liquid structure resulting from the model potential. In particular, the calculations reproduce the maximum amplitude of the pre-edge feature at 533 eV (although slightly overestimating its width) indicating that MB-pol accurately captures the experimentally measured disorder in the hydrogen bonding network [see Fig. 3(d)]. It is important to note that condensed phase structure and thermodynamics were not directly included in the development of MB-pol, so accuracy in reproducing these properties reflects the importance of close-range low-order many-body effects in aqueous systems.

As alluded to earlier in the review, water is characterized by different types of anomalies, which are not found in simple liquids. These include, for example, the density maximum at 4° C¹⁴¹ and the minimum in isothermal compressibility at 46.5° C.¹⁴² The MB-pol model has been successful at reproducing many of these anomalies where standard DFT-based methods at the GGA level of accuracy typically fail,¹⁴³ and effective MM potentials have experienced significant challenges. In part, this originates from the lack of transferability of effective MM potentials as they are often parameterized to work at a specific region of the phase diagram.¹⁴⁴ In Fig. 3(e), we show the isothermal compressibility of liquid water from ambient temperatures to the supercooled regime, comparing MB-pol to experimental measurements. The agreement is extremely impressive. In addition to the well-established minimum at 46.5° C, MB-pol also reproduces the value of the measured compressibility at the tentative maximum at around −44^° C (albeit with large error bars due to long relaxation times near the Widom line).^145,146 The quality of the potential energy surface sampled by MB-pol is clearly superior to the vast majority of MM models and commonly used XC functionals on the market.

Among the many interesting properties of water, perhaps the most challenging is determining spectroscopic properties both in the bulk and at interfaces. Specifically, vibrational spectra in liquid water are sensitive to both the details of microscopic structure, including properties of the hydrogen bonding network, as well as the potential energy surface along specific vibrational modes. Early studies validating the MB-pol potential demonstrated that the model successfully reproduces experimentally measured features of the bulk vibrational (infrared and Raman) spectra, within the limits of the approximate quantum dynamics method used, demonstrating that liquid water simulations with MB-pol produce accurate hydrogen bond dynamics.^147,148 Additionally, the ability of MB-pol to reproduce the features of the experimental SFG spectrum of the air–water interface (within the limits of the approximate quantum dynamics used) is a noteworthy achievement. In particular, the MB-pol calculations predicted the lack of a positive feature around 3000 cm⁻¹ prior to the correction of experimental artifacts between 2011 and 2015.¹⁴⁹ These examples further show that the bottom-up approach involving converging the MBE expansion ensures the transferability of the potential.

Finally, we would like to showcase recent extensions of MB-pol (MB-nrg) to go beyond bulk water and allow for modeling solutions such as ions in water.^{101,130–137} An instructive example is the dilute aqueous Na⁺ solution, for which pairwise classical potentials tend to over-structure the solvation shells around Na⁺, while classical polarizable models under-structure the Na⁺ solvation, predicting an almost continuous transition between the first and second shells.¹³⁷ The accurate description of close-range quantum mechanical effects, such as charge transfer and charge penetration, at the 2B and 3B levels provided by a full many-body potential¹⁰¹ allows MB-nrg to predict the correct solvation structure, and thus reproduce the experimentally measured extended x-ray absorption fine structure (EXAFS) spectrum with unparalleled accuracy.^137, Figure 3(f) compares the K-edge EXAFS spectra for water solvating the Na⁺ ion obtained from MB-nrg and experiments. Interestingly, in the case of Na⁺ solvation, it was found that the MB-nrg potential without 3B corrections produced a solvation structure similar to the classical polarizable potential. EXAFS spectra computed for solutions of halides and other alkali ions using the MB-nrg potentials are presented in several references, demonstrating the applicability of the approach to each system.^134–137 An MB-nrg model for treating solvated hydrocarbons, such as methane, has also been recently reported.^150,151 Inspired by both experiments¹⁵² and previous DFT-based AIMD simulations,¹⁵³ it was shown that under pressure, methane picks up a dipole moment which leads to its enhanced solubility. In this case, the use of MB-nrg potentials allows for running larger and longer simulations than the far more expensive AIMD approaches permit, offering the possibility of obtaining accurate predictions on the corresponding thermodynamics.¹⁵¹ 

Although many-body polarizable potentials, such as MB-pol, allow for an extremely accurate characterization of the structural, dynamical, and spectroscopic properties of bulk water, they do not allow for bond-breaking and bond-formation. One of the most fundamental processes in aqueous chemistry that requires chemistry to occur in water is the dissociation of water into its constituent ions, proton and hydroxide.^154,155 This equilibrium determines the pH of water and has numerous implications for the biochemistry of biomolecules in solution. A lot of our understanding of water dissociation and the diffusion of protons and hydroxide ions in solution has come from molecular simulations that allow for reactive chemistry to occur. In this regard, DFT-based AIMD simulations of liquid water and its ionic products have been instrumental over the last three decades in providing a molecular lens into the structure of the excess proton in water. These simulations have guided and motivated state-of-the-art spectroscopy experiments.^41,156 As alluded to earlier, for chemical reactions such as those involving the breaking of covalent bonds, the standard DFT simulations are hindered by short timescales and small system sizes.

In addition to the ML-based force fields designed to overcome these challenges that will be discussed in Sec. II B, there are also empirical potentials that have been developed to study water dissociation. Attempts to construct dissociative water potentials date back to work almost four decades ago by Stillinger and Rahman.^157–160 These initial studies laid the groundwork for one of the most popular dissociative water potentials, which describes both intramolecular and intermolecular interactions associated with water and its constituent ions, namely, the OSS family.¹⁶¹ The functional forms of these potentials include two-body radial, three-body angular contributions, and polarizable oxygen ions. The first form of the OSS potential was shown to be rather promising at characterizing potential energy surfaces in protonated clusters. Over the last two decades, it has been systematically improved to allow for modeling proton and hydroxide transfer in the bulk.^162–164 The left panel of Fig. 4(a) illustrates the mean square displacement (MSD) and, thus, the inferred diffusion constants for the excess proton, hydroxide, and neutral water, obtained from the OSS2 potential by Rasaiah and coworkers.¹⁶³ It is first important to note that the potential successfully predicts the relative differences in the diffusion constants for the three species, namely, that $D H +$ > $D O H −$ > $D H 2 O$⁠, consistent with experiments. Instead, the right panel of Fig. 4(a) shows the magnitude of the diffusion constants for the proton and hydroxide, at different temperatures again compared with experimental results. While the trends in the change of the diffusion constants as a function of temperature are consistently reproduced by the model, the differences appear to get more pronounced at larger temperatures. Being able to quantitatively predict these types of dynamical properties of water's constituent ions at different thermodynamic conditions remains an open challenge.

Another family of dissociative potentials that have caught some momentum in the literature are those developed by Garofalini and coworkers to deal with water dissociation in the bulk and near inorganic metal oxide interfaces such as silica.¹⁶⁵ Similar to the OSS potential, water dissociation is facilitated by incorporating intramolecular interactions, which encompass two and three-body terms. This approach has been employed in recent studies to investigate the mechanisms involved in water dissociation. Garofalini and coworkers show that liquid water is characterized by a larger fraction of transient autoionization events, similar to previous studies where this was observed upon inclusion of nuclear quantum effects.¹⁶⁶ They also point out the challenge in identifying the relevant collective variables that ultimately lead to separation of the hydronium and hydroxide, a topic that will be discussed at greater length later in the review. Over the last decade, Wiedemair and coworkers have extended the original potential proposed by Garofalini, in an effort to improve various dynamical properties by replacing the non-Coulombic Morse-like potential with either a Lennard-Jones or Buckhingham potential, allowing for it to be of greater practical use in other common MD softwares.¹⁶⁷ 

In the last two decades, the Voth group has pioneered the development of reactive potentials based on the multi-state empirical valence bond formalism (MS-EVB). The EVB-type potentials have been typically trained with different levels of DFT simulations and have allowed for examining problems such as proton transfer in bulk liquid water with impressive detail and accuracy.^168–173 Since this has already been tackled in several previous reviews,^171,174 here we focus on highlighting a specific problem that has been the subject of raging debate and regards the propensity of water's constituent ions for the air–water interface and the apparent negative charge near hydrophobic interfaces.^175,176 Figure 4(b) illustrates the potential of mean force (PMF) calculations performed by the Voth group, which show that the excess proton has a slight propensity for the surface of water while the hydroxide ion is repelled.^177,178 On the other hand, previous studies using DFT-based AIMD simulations reached different conclusions regarding the propensity of the hydronium and hydroxide ions for the surface of water—the latter, hydroxide ion, is instead found to be weakly attracted while the proton has no affinity for the air–water interface.^179,180 Such a dual—acidic or basic—behavior of the surface of water has also recently been the subject of a large-scale QM/MM study, attempting at resolving the controversial finding that the excess proton presents a higher affinity for the surface compared to the hydroxide.¹⁸¹ The origin of this behavior is attributed to differences in the local solvation structures which are challenging to sample with standard DFT-based simulations.

The Hertzfeld group have also developed LEWIS-based potentials, which include valence electrons as semi-classical particles interacting with each other through pairwise potentials.^182,183 The LEWIS model tends to predict a much higher affinity of hydroxide ions for the surface of water compared to both MS-EVB- and DFT-based simulations [see Fig. 4(c)].^177,184 In addition, potentials from the Netz group, by including several extra bond-stretching and angle potentials, as well as altering the point charges, have been interfaced with the SPC/E water model. While this does not allow for studying the Grotthuss mechanism, thermodynamic properties such as solvation energies and surface activities are accurately reproduced.¹⁸⁵ This is illustrated in Fig. 4(d), which compares the surface activities extracted from experiments, simulations and through PMF calculations, for HCl, NaOH, and NaCl. One observes that the proton is surface active, while the hydroxide is not, consistent with the predictions of the PMFs by Voth shown in Fig. 4(c). Needless to say, it is clear that this is an extremely challenging problem and, depending on the flavor of the underlying potential that is used, can lead to very different qualitative and quantitative conclusions.

Finally, another dissociative potential of relevance is the ReaxFF force-field developed by van Duin, Goddard III, and coworkers.^186–188 This is a bond-order-based force field, typically employed for molecular dynamics simulations involving chemical reactions. Conventional force fields face limitations in representing chemical reactions due to their explicit bond definition requirements. In contrast, ReaxFF employs bond orders rather than explicit bonds, enabling continuous bond formation and breaking, making it more suitable for modeling chemical reactions. The first-generation ReaxFF water force field (water-2010) underestimated the bulk density of liquid water by $∼ 8 %$ at ambient conditions¹⁸⁹ (an issue also occurring with many DFT XC functionals), but also incorrectly predicted the order of the diffusion coefficients as D(H₂O) < D(H₃O⁺) < D(OH⁻).¹⁹⁰ 

The second-generation of ReaxFF force field¹⁹¹(water-2017) was developed by fitting energies and structures to QM data (as in water-2010) and, in addition, by explicitly including both the experimental density of water and the correct diffusion coefficients of H₂O, H₃O⁺, OH⁻ in the training procedure. This resolved the aforementioned issues and resulted in a more accurate density of bulk water, as well as the correct order of the diffusion coefficients: D(H₂O) < D(OH⁻) < D(H₃O⁺). This force field has also been shown to describe water dissociation and the Grotthuss mechanism underlying the propagation of protons. However, in comparison with water-2010, the water-2017 force field slightly misjudges the location and intensity of the first peak in the O–O and H–O RDFs, even though they are comparable to what is predicted by common classical water models, such as SPC/E, TIP3P, and TIP4P-2005.¹⁹¹ It was hypothesized that this was caused by including the density of bulk water in the training data, which might have resulted in an overestimation of the intermolecular interactions. The reactive nature of this force field has allowed it to be used to simulate water dissociation at liquid–solid interfaces.¹⁹² 

One of the challenges with developing accurate potential energy surfaces (PESs) for use in empirical potentials like the ones described earlier is the construction of an appropriate functional form for the interacting particles. The reactive potentials discussed in the Reactive Empirical Potentials can often have very complicated functional forms. This has triggered the use of models for which the knowledge of the functional form is not needed, with Neural Networks (NN) being a natural choice, as they are universal approximators of any mathematical function.¹⁹³ Therefore, in principle, by properly optimizing the internal weights of the NN, one can fit any PES as a function of the atomic positions without having to explicitly separate out the reactive or bonded/non-bonded nature of the interactions.

The vast majority of NN potentials for water used in molecular dynamics applications builds on the seminal work by Behler and Parrinello,¹⁹⁴ who developed a generalized NN representation of high-dimensional PES that was trained on DFT data. The reader is referred to several other detailed reviews on the topic.^195–197 Here, we outline the essential principles of the method. The idea is to express the total energy of a system as a sum of individual atomic contributions E_i which depend on the local chemical environment. Each atom is dressed with its own NN so that the total energy of the system is given by $E = ∑ i = 1 N E i$⁠. In order to be physically meaningful, each E_i must be invariant under the permutation of identical atoms in the local environment and under the rotation/translation of their coordinates. However, learning these symmetries (although theoretically possible) would be extremely costly since the algorithm would need much data points during the training. Therefore, it is a common practice to employ some kind of symmetry functions (G_i) that encode the local atomic structure and satisfy these symmetry conditions, providing a direct link between the atomic energies E_i and the local atomic environments. Specifically, in the Behler–Parrinello NN (BPNN) framework, the atomic coordinates are mapped onto a set of two- and three-body symmetry functions.

The structure of the atomic NN shown in Fig. 1(b) is inspired by the architecture of artificial neurons which consist of nodes organized in various layers. The symmetry function coordinates described earlier are used as input for the first layer, while the last layer produces as output the atomic energy. The weights of the connecting nodes are optimized to reduce the error of these energies with respect to electronic structure calculations.¹⁹⁷ 

Several versions of NN potentials of water have been developed over the last decade which began by originally fitting the potential energy surface of the water dimer based on environment-dependent atomic energies and charges.¹⁹⁸ Over time, this potential was refined to handle water clusters highlighting the importance of dispersion interactions.¹⁹⁹ These NN potentials clearly allow for the simulation of systems over much longer timescales. Naturally, however, the quality of the potential depends on where exactly the training dataset comes from. For example, training an NN potential using DFT-based molecular dynamics without dispersion corrections leads to overstructured and glassy water dynamics.²⁰⁰ 

Figure 5(a) illustrates a powerful example of using different types of NN potentials to study the behavior of the density maximum of liquid water, as well as the thermodynamics associated with the ice-water equilibrium. Behler and Dellago developed a series of NN potentials with which they demonstrated that the inclusion of dispersion corrections with different DFT functionals reflects the correct existence of a density maximum.²⁰⁰ In this work, they also computed melting temperatures for systems like the one shown in Fig. 5(a), consisting of over 2000 water molecules and with numerous simulations on the nanosecond timescale. Table in Fig. 5(a) shows that the melting temperatures are significantly improved upon the inclusion of dispersion corrections, in comparison with the experiments. Again, these conclusions are based on the possibility of conducting large-scale and long-time (several nanoseconds) simulations with the NN potentials.

The use of NN-potentials for water has opened up exciting applications in which they are coupled with path integral molecular dynamics (PIMD). A nice example of this synergy is a recent application involving Ceriotti and Behler that uses a NN-potential of water combined with path integral techniques to predict the thermodynamic properties of liquid water, as well as hexagonal and cubic ice²⁰¹ By training a NN with a hybrid functional, these simulations rigorously take into account NQEs, the disorder of the protons, and anharmonic fluctuations. The upper panel of Fig. 5(b) shows the density isobars computed with the NN for liquid water, hexagonal, and cubic ice. The temperature of maximum density is in perfect agreement with the experiments, while the density isobars for all studied systems are within 3% of the experiments. The lower panel of Fig. 5(b), instead, shows the coupling of the NN with path integral simulations being used to determine the role of NQEs in controlling the structure of liquid water. The PIMD simulations coupled with the NN show excellent agreement with experiments for all three pair-correlation functions. Furthermore, as pointed out earlier, NQEs are essential for capturing the delocalization of the proton along the hydrogen bond as seen in the bottom two panels of the RDFs.

In recent years, Weinan and Car have extended the scope of standard NN potentials with the development of deep neural network potentials (DNN).^202–205 DNN potentials have emerged in order to tackle some critical limitations with the classical NN methodology, including the somewhat ad hoc choice of symmetry functions. Furthermore, DNNs have demonstrated a substantial enhancement in the efficiency of the learning phase for potential functions. This improvement is achieved through the incorporation of optimal loss functions and local decomposition methods, which enable training DNNs on relatively small systems while maintaining their applicability to much larger systems. It is noteworthy that the computational cost of employing DNNs scales linearly with system size, primarily because they are highly parallelizable due to the earlier mentioned local decomposition approach. Their adoption has expanded the scope of applications in simulating various aqueous systems. Notably, DNNs exhibit remarkable accuracy in reproducing average energy, density, and RDFs across a range of water and ice systems when compared to AIMD simulations.²⁰² 

Car and coworkers have recently also used the DNN framework to not only learn structural properties, but also the environment-dependent polarizability tensor of water molecules.²⁰⁶ This, in turn, has allowed determining the Raman spectrum of liquid water.^207, Figure 5(c) shows the Raman spectrum for bulk water obtained from a DNN which is used to sample the configurations during a molecular dynamics simulation, but also predicts the Wannier functions necessary for computing Raman intensities. The DNN approach was also used to develop a model that faithfully reproduces the potential energy surface of water compared to the SCAN (Strongly Constrained and Appropriately Normed) approximation of DFT over a wide range of temperature and pressure conditions (up to 50 GPa and 2400 K). The computational efficiency of this model enables the prediction of the phase diagram of water through MD simulations, in good agreement with experimental data.²⁰⁸ Another recent application of this methodology was investigation of the possibility of a second critical point in water under supercooling. Debendetti and Car trained a DNN to represent the ab initio potential energy surface coming from DFT calculations using the SCAN functional.²⁰⁹ Although the calculations are rather limited in terms of both system size and timescales, they strongly indicate the validity and the existence of the liquid–liquid critical point (LLCP), consistent with previous classical point-charge models.²⁹ Furthermore, Paesani and coworkers have significantly expanded the scope of MB-pol by coupling it with DNN, which has recently allowed the exploration of the phase diagram of water.²¹⁰ 

Finally, highlighting further opportunities and possibilities in the future, Car and Baroni present a very recent study on using a DNN of liquid water to predict the viscosities and thermal conductivities from first principles using the Green–Kubo theory.^211,212 Figure 5(d) shows the comparison of the behavior of thermal conductivity for bulk water across a range of temperatures with different functionals. While the temperature dependence is in moderately good agreement with the experiments and there is a shift in the right direction when going from standard GGAs to meta-GGAs, it is clear that the quality of the generated DNN ultimately depends on the electronic structure theory used in the training.

Fully dissociative NN and DNN allowing for water ionization and simulating the excess proton and hydroxide simultaneously and, in general, chemical reactions, have only recently been gathering momentum. Behler and coworkers have pioneered some of the early NNs stemming from neutral water that have displayed encouraging results on capturing the potential energy surface of gas-phase protonated clusters.^213,214 Interestingly, NN potentials have been developed to allow simulations of sodium hydroxide solutions (NaOH) at different concentrations, which have also been coupled with PIMD simulations to investigate the role of NQEs.^215,216 Recently, Markland and coworkers have demonstrated very promising steps in the development of an NN potential that allows for simultaneous modeling of the proton and hydroxide ion in water.²¹⁷ Laage and coworkers have also very recently trained NN potentials for simulating the excess proton in the bulk as well as at interfaces the latter of which is also used to construct SFG vibrational spectra.^218,219 In the last year, Car and coworkers have reported a water dissociable DNN potential which is used to compute the pKw of water.²²⁰ The findings are extremely encouraging and open up new possibilities for studying water chemistry at interfaces, under confinement, and in organic and salt solutions. For example, a recent AIMD study by Di Pino and coworkers showed that the pKw of water is rather insensitive to confinement until one gets to very small cluster sizes where the local solvation of the hydronium and hydroxide changes.²⁸ We expect that over the next decade, there will be new advancements on the methodological side that will allow for routine application of DNNs to study water ionization in the bulk as well as at interfaces and under confinement such as in this recent application.

Another interesting family of models not relying on NN are the general-purpose gradient-domain machine learning potentials (GDML), which have originally shown tremendous promise for medium to large-size single molecules.²²¹ In this approach, the coordinates are mapped into the eigenvalues of the Coulomb Matrix, which contains the inverse distances between all distinct pairs of atoms in the system. The functional relationship between atomic coordinates and interatomic forces is then obtained by employing a generalization of the Kernel Ridge regression technique. GDML models using small training sets have been shown to perform similarly to potentials requiring large training sets. However, GDML applicability is limited to the same system it was trained on.²²² A recent GDML-driven many-body expansion framework (mbGDML) has enabled size transferability, opening the doors for the simulation of bulk water. Predictions on static clusters of up to 16 monomers had energy errors of less than 0.251 kcal mol⁻¹ per monomer, and the simulated RDFs with mbGDML for bulk water agree remarkably well with the reference experimental data.²²² 

An important aspect of training NN potentials is generating accurate and sufficiently large enough datasets that adequately represent the environments that are of interest to the physical problem of study. The wider diversity of environments that are seen during the training step, the more likely it will be stable across different conditions. A commonly employed strategy involves generating configurations across a broader range of temperature and pressure than those being considered. In particular, adding to the dataset environments at temperatures and pressures that are higher than that under consideration reinforces the model and leads to more accurate predictions. Litman and coworkers recently trained a NN for different electrolytes at the surface of water performing ab initio calculations on slab geometries.²²³ A challenge within this context that remains an open problem is how transferable a potential that has been trained for the bulk is to anisotropic conditions such as an interface and vice versa. Another interesting and open area of research in the field is how one takes advantage of data that have been trained with different levels of electronic structure theory using techniques involving fine-tuning transfer learning approaches.²²⁴ 

On the more practical side of things, it is worth also comparing the computational performances of DNN compared to classical potentials. A critical game changer within this context has been the ability to run simulations on GPUs. As a benchmark, we compare the performance of conducting a simulation of approximately 300 water molecules under supercooled conditions close to the putative second critical point of water. Using the GROMACS code²²⁵ on a single GPU with 8 OpenMP threads, one can obtain about 1 microsecond a day while in the case of the DNN, using 1 GPU one obtains roughly 1 ns per day. For large enough systems, the scaling for DNN potentials grows linearly with the number of GPUs. Although the performance of the classical potentials is clearly better, the DNN should be compared to the situation only decade ago when one would need to wait several months to obtain 100s of picoseconds with AIMD approaches.

As alluded to earlier, for many problems that involve probing fluctuations of the hydrogen bond network—such as exploring the phase diagram in the supercooled regime or studying the breaking of bonds in the ionization of water—there is a computational challenge of dealing with slow timescales of the underlying activated processes. In the literature, this challenge is commonly termed as the sampling problem. At the same time, molecular simulations in the condensed phase often yield large amounts of data that are extremely challenging to interpret and understand using only chemical imagination. In this section, we will highlight efforts that have been made in recent years to overcome the sampling problem, with specific examples associated with water chemistry. Afterward, we will present recent developments in data-mining techniques that are providing new insights into the chemical physics and physical chemistry of aqueous systems.

In the last two decades, numerous methods have been developed to deal with the sampling problem. The reader is referred to other detailed reviews in the topic.^{25,31,226–228} Here, we briefly overview some of the main applications in this area that have been instrumental in driving our understanding of various aqueous-related phenomena. Specifically, the equilibrium between water and its constituent ions (proton and hydroxide), and more generally acids and their conjugated bases, requires fluctuations that are rather slow and cannot be achieved by brute force simulations. Through Car-Parrinello molecular dynamics (CPMD) simulations, Trout, Parrinello, and Sprik computed the free energy associated with water dissociation using the Blue-Moon ensemble technique.^229–231 The principle idea was to constrain the system along a pre-defined reaction coordinate, typically chosen by chemical intuition, with an external bias potential. Although these early studies were limited by short simulation times, they demonstrated that most of the free energy of ionization arises from breaking the O–H covalent bond. This observation is essentially reproduced using more accurate DFT functionals such as SCAN.²³² These constrained simulations, however, while effective at extracting thermodynamic quantities, do not say anything about dynamical mechanisms.

Chandler and Parrinello combined Transition Path Sampling²³³ together with CPMD simulations to study the pathways associated with water ionization.²³⁴ By generating reactive trajectories, they suggested that the breaking of the O–H covalent bond involves large electric field fluctuations followed by the reorganization of the hydrogen-bond network that separates the hydronium and hydroxide ion. Several other studies have used metadynamics simulations to investigate acid–base equilibria, for example, in acetic and carbonic acid.^235,236 These studies have shown that, due to the collective nature of the fluctuations of the hydrogen-bond network, coordinates on both short and intermediate-length scales are required. The time-reversed process of the ionization, namely, the recombination of the hydronium and hydroxide ions, was revisited by Hassanali and Parrinello.²³⁷ The AIMD simulations confirmed the importance of medium-range correlations in the water network, as revealed by fluctuations of picosecond compressions of water wires that connect the two ions as illustrated in Fig. 6(a). Though to a lesser extent, these water wires compressions are also observed during proton transfer events in crystalline ice phases.²³⁸ In addition, the amphiphilic nature of these ions^239,240 implies that there are rather specific changes in the solvation patterns required on the donating/accepting hydrogen-bond directions of the ions.

In the last decade, machine learning and information science techniques have also been integrated into atomistic simulations to automatically discover relevant reaction coordinates, minimize chemical bias, and, subsequently, characterize the free energy landscapes of complex systems.²⁴¹ In addition, identifying and constructing reaction coordinates involve a dimensionality reduction which can lead to uncontrolled information loss.²⁴² A recent example that attempts to overcome this problem for aqueous chemistry is the work by van Erp and coworkers²⁴³ where they revisit the ionization of water using a regression tree analysis, pictorially shown in Fig. 6(b). This approach allows the determination of the relative weights of different variables (starting from 138 collective variables) and a more rigorous distinction between the reactive and non-reactive trajectories in the ionization event. In addition to the role of the vibrations of the proton wire shown in Fig. 6(a), which the decision tree ranks as one of the most important, van Erp and coworkers also demonstrate a non-insignificant role of the variables such as the number of accepting hydrogen bonds along the wire and tetrahedrality.

Chemical reactions involving proton transfer, as is the case of water ionization, are particularly challenging due to the fact that since protons move through the hydrogen-bond network via the Grotthuss mechanism, their identity always changes.^11,40 Thus, a good reaction coordinate must automatically account for this feature. We highlight a relatively recent work by Pietrucci and coauthors that makes a significant stride in dealing with this problem, in which they develop chemical-topology-based coordinates.^244,245 For a more detailed review of the methodological aspects and developments of graph-based approaches for studying chemical reactions, the reader is directed to other papers.²⁴⁶ The idea is that chemical species in a reactant and product state can be defined by their bonding topology without a need to describe bonding interactions between specific atoms. This is depicted in the top panel of Fig. 6(c) which illustrates the transformation from formamide to ammonia and carbon-monoxide. These collective variables have been interfaced with techniques such as metadynamics and have been used to explore chemical reactions relevant to problems associated with the origin of life.^245,247 A key finding of these simulations is the importance of solvent reorganization which changes the thermodynamics and dynamics of chemical reactions quite significantly. Specifically, the two bottom panels of Fig. 6(c) show free energy surfaces coming from metadynamics simulations—formamide chemistry in the gas phase (left) vs in solution (right) resulting in completely different pathways, intermediates, and products.

Another very active and growing area of research currently is the use of unsupervised learning approaches to characterize and understand patterns that arise in molecular simulations without prior imposition of knowledge often introduced by chemical bias.²⁴¹ In this context, a wide range of local atomic descriptors are now being used to encode information about atomic environments because they preserve important symmetries. A particularly popular one, that is used in the study of liquid water for example, is the smooth overlap of atomic positions (SOAPs) which essentially expands the atomic neighbor density of a chemical species onto a basis of radial basis functions and spherical harmonics.^248,249 These types of atomic descriptors are high-dimensional vectors encoding details of the local environment. In order to extract useful and interpretable information from them, one needs to perform some form of dimensionality reduction and subsequently project along the relevant degrees of freedom using clustering. Figure 6(d) illustrates an example of this approach by Cheng and coworkers: SOAP descriptors of bulk liquid water are built and later used for a principal component analysis.²⁵⁰ These authors also try to relate the fluctuations in liquid water to milestone structures involving different phases of ice. In particular, they demonstrate that one can think about fluctuations in liquid water as transiently forming different local structures resembling phases of ice. Pavan and coworkers have also demonstrated similar ideas, again by taking advantage of the generality of the SOAP descriptors, comparing similarity measures of different empirical potentials of liquid water to phases of ice.²⁵¹ 

Recently, our group has also taken important steps in this direction in an effort to understand the thermodynamic landscape of liquid water and how it changes around model hydrophobic polymers, amino acids, peptide groups, and at interfaces.^252–257 Specifically, using the SOAP descriptors we have investigated the number of independent degrees of freedom [often referred to as the intrinsic dimension (ID) of a dataset] needed to characterize the fluctuations in liquid water at room temperature and, thereafter, we extracted the high dimensional free energies.²⁵⁸ One important observation is that even at the local structure level, for example, between the first and second solvation shells, the ID of the system is quite large. This implies that a correct description of the underlying thermodynamics and dynamics of the system requires looking at many orthogonal degrees of freedom concertedly. Contrary to current descriptions of water in terms of two-state liquid, the left panel of Fig. 7 shows that liquid water at room temperature is characterized by a single broad minimum. Furthermore, we do not find any evidence for two states or stable local structures corresponding to low-density (LDL) and high-density liquid (HDL). More recently, we have applied similar techniques to address the question regarding the structure of the proton in liquid water, typically discussed in the literature in terms of a competition between idealized limiting states, namely, the Eigen and Zundel.²⁵⁹ Our agnostic approach, instead, shows that the Eigen and Zundel are neither limiting nor stable thermodynamic states. Contrasting the two panels of Fig. 7 shows that the simulations with the excess proton (concentrated HCl to be specific) lead to the creation of two additional minima on the free energy landscape. The excess proton in water is best seen as a charged topological defect that strongly perturbs its local environment, leading to an enhancement of the concentration of neutral water defects.

Another topic that has generated large interest over the last decade is the relation between reorientational dynamics of water molecules and the reorganization of the hydrogen-bond network in liquid water.^260–263 While it is now well known that water molecules reorient mostly by sudden large-amplitude angular jumps rather than by slow rotation through small diffusive steps as was previously believed,^264–267 the collective and concerted nature of these jumps has only been indicated by the spectroscopy experiments of water dynamics.^{262,263,268–270} In order to overcome the challenge of identifying the appropriate collective variables that capture the underlying reorganization of the hydrogen-bond network, our group has developed an automatic angular-jump detection tool which enabled us to observe numerous small and large angular jumps happening simultaneously in the system.²⁷¹ This large data analysis that correctly characterizes the ongoing dynamics without a priori knowledge of the underlying chemistry, has proven to be a powerful tool in unraveling the microscopic mechanism that involves a cascade of hydrogen-bond fluctuations underlying angular jumps. Challenges lying ahead involve understanding how the collective reorientational dynamics can be used as a signal of various chemical processes in liquid water such as ion recombination or for characterizing transitions between different phases, as recently shown for the supercooled water.²⁷² 

The goal of this review is to give a broad overview of the historical context of studying aqueous chemistry in solution, both from the perspective of the use of DFT-based AIMD simulations as well as modern empirical potentials with a focus on dissociative schemes. This background motivated our discussion on the development of machine-learning potentials and, more generally, data-driven approaches to both model aqueous systems with greater accuracy and sophistication as well as to drive the learning of complex phenomena in these systems. Where possible, we have made an attempt to elucidate critical challenges that have arisen in different types of problems relevant to the physical chemistry and chemical physics of solvation.

The growth of machine-learning approaches to study aqueous chemistry, in particular the dissociation of water, offers enormous potential for future applications. Specifically, how these models can be extended to tackle other problems involving salt solutions,²⁷³ organic molecules such as amino acids, proteins, and DNA,¹⁸ and, finally, processes out of equilibrium such as in external electric fields,^66,274–281 remains an open question. The challenge here involves both the generation of accurate datasets for training as well as assessing the transferability of existing models. For example, the exchange of protons in solution is a critical step underlying isotope fractionation of biomolecules in water, which is thought to be catalyzed by negatively charged hydroxide ions.²⁸² However, the underlying mechanisms by which this happens still remain uncharted territory.

In the context of externally driven aqueous solutions, large external electric fields, for example, have been shown to enhance autoprotolysis in solution.^66,275,276 It is clear that the breadth and depth of computer simulations that have been achieved with standard methods such as AIMD will play a critical role in helping push forward the development and application of machine-learning-based approaches, as recently reported in the literature.²⁸³ Another very interesting and challenging area that will certainly benefit from data-driven approaches is excited-state chemistry in solution, which presents many other methodological issues at the moment that are typically absent in the ground-state such as the near-sightedness of the electronic degrees of freedom. The reader is referred to the following review and references therein for more details.²⁸⁴ 

In addition to the use of ML-based potentials to simulate aqueous chemistry with higher accuracy and on larger system sizes, another active area that we believe will continue to grow is the use of both supervised and unsupervised learning to unravel the complexity of chemical processes in aqueous solutions. While these approaches allow an agnostic inference of models based on the underlying structure of the data, chemical and physical interpretability is often a critical missing link. Therefore, attempts to provide physical interpretability, for example, to atomic descriptors used in neural networks (NN), is an area that requires more attention.²⁵⁷ This understanding has important implications on how to translate physics-based interaction models into the structure of NNs.

While the focus of this review has been on motivating data-driven potentials relevant for first-principles simulations, there are also other examples of ML being used to improve coarse-grained potentials. The mW model developed by the Molinero group,²⁸⁵ despite being a coarse-grained model, outperforms many of the other empirical potentials such as SPC and TIP in predicting thermodynamic properties such as surface-tension. While this model does not capture realistically the dynamics of water, it has been successful at studying phenomena such ice-nucleation and thermodynamics of water under confinement.^286–292 More recently, it has also been interfaced to combining it with ions, molecules and polymers to tackle more complex aqueous solutions as well as systems involving interfaces such as membranes.^293–297 ML approaches have the potential to also enable coarse-grained models such as mW. To highlight one specific case, the coarse-grained ML potentials developed by Sankaranarayanan and coworkers have devised an analytical force-field for mono-atomic water based on the Pauling bond-order concept, containing up to 4-body terms.²⁹⁸ In particular, they have used MD trajectories and experimental data to train their NN for the optimal parameters to feed their MB potential, named ML-BOP, and also the mW model. The possibility to combine ML approaches with coarse-grained potentials like mW has potential to open up other interesting applications of aqueous solutions.

D.B., C.K.E., G.D.M., M.M., and A.H. thank the European Commission for funding on the ERC Grant No. HyBOP 101043272. K.A. thanks the Center of International Science and Technology Cooperation (CISTC) of the Iranian Vice-Presidency of Science and Technology for the postdoctoral grant.

The authors have no conflicts to disclose.

Debarshi Banerjee: Conceptualization (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Khatereh Azizi: Writing – original draft (equal); Writing – review & editing (equal). Colin K. Egan: Writing – original draft (equal); Writing – review & editing (equal). Edward Danquah Donkor: Writing – original draft (equal); Writing – review & editing (equal). Cesare Malosso: Writing – original draft (equal); Writing – review & editing (equal). Solana Di Pino: Writing – original draft (equal); Writing – review & editing (equal). Gonzalo Diaz Miron: Writing – original draft (equal); Writing – review & editing (equal). Martina Stella: Writing – original draft (equal); Writing – review & editing (equal). Giulia Sormani: Writing – original draft (equal); Writing – review & editing (equal). Germaine Neza Hozana: Writing – original draft (equal); Writing – review & editing (equal). Marta Monti: Writing – original draft (equal); Writing – review & editing (equal). Uriel N. Morzan: Writing – original draft (equal); Writing – review & editing (equal). Alex Rodriguez: Writing – original draft (equal); Writing – review & editing (equal). Giuseppe Cassone: Writing – original draft (equal); Writing – review & editing (equal). Asja Jelic: Writing – original draft (equal); Writing – review & editing (equal). Damian Scherlis: Writing – original draft (equal); Writing – review & editing (equal). Ali Hassanali: Conceptualization (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.