Computational studies of liquid water and its phase transition into vapor have traditionally been performed using classical water models. Here, we utilize the Deep Potential methodology—a machine learning approach—to study this ubiquitous phase transition, starting from the phase diagram in the liquid–vapor coexistence regime. The machine learning model is trained on ab initio energies and forces based on the SCAN density functional, which has been previously shown to reproduce solid phases and other properties of water. Here, we compute the surface tension, saturation pressure, and enthalpy of vaporization for a range of temperatures spanning from 300 to 600 K and evaluate the Deep Potential model performance against experimental results and the semiempirical TIP4P/2005 classical model. Moreover, by employing the seeding technique, we evaluate the free energy barrier and nucleation rate at negative pressures for the isotherm of 296.4 K. We find that the nucleation rates obtained from the Deep Potential model deviate from those computed for the TIP4P/2005 water model due to an underestimation in the surface tension from the Deep Potential model. From analysis of the seeding simulations, we also evaluate the Tolman length for the Deep Potential water model, which is (0.091 ± 0.008) nm at 296.4 K. Finally, we identify that water molecules display a preferential orientation in the liquid–vapor interface, in which H atoms tend to point toward the vapor phase to maximize the enthalpic gain of interfacial molecules. We find that this behavior is more pronounced for planar interfaces than for the curved interfaces in bubbles. This work represents the first application of Deep Potential models to the study of liquid–vapor coexistence and water cavitation.

Water is a fundamental substance, crucial for life and relevant in many environmental, engineering, and biological processes.1–9 Due to this, the past decades have seen a significant effort devoted to the development of models to reproduce the behavior of water in computer simulations.10–22 Although some of these models account for flexibility and polarizability,18–21,23 the most widely employed models for water are rigid and nonpolarizable. These include the TIP4P/2005,10 TIP4P/ICE,11 SPC/E,15 TIP3P,13 and TIP4P-Ew24 among others.16,17 These empirical models have parameters obtained by fitting to experimentally measured properties, such as coexistence lines between thermodynamic phases or critical points. They are frequently used to describe ionic solutions25–30 and for biomolecular simulations31,32 as well as other applications.33–36 

In contrast to the classic semiempirical approach to water modeling, ab initio models are determined from first principles and therefore do not require fitting to experimental data.37 Traditionally, this approach has not been applied to large systems due to its computational cost.38,39 Nonetheless, recent advances in machine learning (ML) have allowed the development of Deep Potential generators40,41 capable of constructing Molecular Dynamics (MD) potentials based on ab initio models. In this work, we use a ML-based model that has successfully recapitulated the different solid phases of water.42 This model has been constructed based on the SCAN quantum mechanical density functional, which successfully reproduces several properties of water.43,44 This approach to ab initio-based models is efficient enough to carry out simulations with tens of thousands of water molecules in a computationally affordable way.45,46

Making use of the machinery provided by ML ab initio-based models, here we study the liquid–vapor coexistence of water by means of MD simulations. The liquid–vapor properties of water are well known from experiments.47 From the computational side, the TIP4P/2005 non-polarizable rigid water model10 has been highly successful at describing the liquid–vapor coexistence properties, reproducing the experimental phase diagram48 and surface tension49–51 reasonably well, as well as transport properties.10,48,52 The TIP4P/2005 model has also been extensively benchmarked in the study of liquid-to-vapor and vapor-to-liquid phase transitions.53–60 Therefore, we compare the Deep Potential Molecular Dynamics (DPMD) water model performance with that of TIP4P/2005 in addition to experimental data, when available.

To test the DPMD model, we first evaluate the liquid–vapor phase diagram and measure properties at equilibrium, such as the surface tension or the enthalpy of vaporization. We also focus on the liquid-to-vapor phase transition at negative pressure, a phenomenon known as cavitation. Under negative pressures, water can remain metastable with respect to the vapor (the most stable phase under these conditions) for a finite amount of time before undergoing a phase transition.61 This results from the fact that a phase transition is an activated process and it requires the formation of a critical bubble, i.e., one that has surmounted a free energy barrier and can continue growing irreversibly without a free energy penalty.34,62 This happens because the formation of a bubble intrinsically requires the formation of a liquid–vapor interface, which comes associated with an energetic cost, the surface tension. When the phase transition is initiated by the formation of water bubbles within the liquid bulk and in the absence of any surface or external agent, the process is termed homogeneous cavitation.

It has been experimentally determined that, at ambient temperature, water can sustain negative pressures of up to −120 MPa before transitioning into the vapor phase.63–67 While experiments have determined the cavitation pressure, which is the pressure at which the phase transition is observed, computational studies using the TIP4P/2005 model have been able to compute the nucleation rate, a crucial quantity to characterize the cavitation process. The nucleation rate is defined as the number of critical clusters formed per unit of time and volume. The nucleation rate obtained in previous studies for the TIP4P/2005 model53–55,58,59 will be used as a reference for our DPMD calculations since there are no reliable experimental data for this quantity.

Aside from the nucleation rate, we also compare in this work the nucleation free energy barrier and the Tolman length,68 a parameter employed to describe the change in surface tension with curvature. Finally, we characterize the orientational distribution of water molecules at the interface. We find that the DPMD model can reproduce well the phase diagram of water but displays a lower surface tension than experimental results or the TIP4P/2005 model. The nucleation rate is consequently greater for the DPMD model compared to the TIP4P/2005 model.

We use the recently developed DPMD model for water42 to perform simulations in the liquid–vapor coexistence regime. The model was generated using an iterative concurrent learning scheme, Deep Potential generator,41 to construct a potential energy landscape based on SCAN,43 a nonempirical functional that recapitulates several properties of water,44 such as molecular geometry and solid structures. The final training set used to construct this model included ice and liquid phases snapshots. In Ref. 42, the phase diagram for this model was calculated for the different ice phases, reaching a reasonable agreement with experimental data.69–71 For computational purposes, we employ the compressed version of this potential, making use of the scheme developed in Ref. 72. Thanks to this approach, we are capable of reaching a computational performance of 5.2 nanoseconds of simulation time per wall clock day for a system of about 10 000 water molecules running with a single 2.8 GHz Intel Ice Lake node using four NVIDIA A100 80 GB GPUs and 28 CPU-cores. Example files of simulations employing this potential are available in the Princeton DataSpace repository https://doi.org/10.34770/ms7d-wm45.

Simulations of the DPMD water model were performed using the LAMMPS package,73 built with the DeePMD-kit.74 Seeding and Direct Coexistence (DC) simulations were performed in the NVT ensemble, keeping the number of particles N, system volume V, and temperature T constant with the Nose–Hoover thermostat.75–77 Additionally, to compute equations of state and to observe crystallization directly at high supersaturations, we performed simulations in the NPT ensemble using the Nose–Hoover barostat.75–77 The cutoff for this potential is 6 Å. The equations of motion were integrated using the velocity-Verlet integrator. The simulation time step was 0.5 fs and the thermostat relaxation time 0.1 ps. In NPT simulations, the barostat relaxation time was 1 ps.

For the DC simulations, a system size of 1024 molecules was used and the density profiles were obtained with at least 10 ns of sampling. Coexistence densities were obtained by fitting the density profile to the following expression:
ρ(z)=ρl+ρv2ρlρv2tanhzz0d,
(1)
where ρl and ρv are the coexistence liquid and vapor densities, respectively, z0 the position of the interface, and d its thickness.
The surface tension was calculated from DC simulations at each temperature according to the Kirkwood–Buff equation,78 
γ=Lz2[Pzz0.5(Pxx+Pyy)],
(2)
where Pii are the diagonal components of the pressure tensor and Lz are the box length in the elongated dimension, perpendicular to the slab interfaces.

We also performed simulations with the TIP4P/2005 water model10 using the GROMACS 4.6.7 Molecular Dynamics package79 in the NPT and NVT ensembles, keeping temperature constant with the velocity-rescale thermostat80 and pressure constant (for NPT simulations) with the Parrinello–Rahman barostat.81 In GROMACS, we integrated the equations of motion using the Leap-Frog integrator.82 The simulation time step was 2 fs, and the thermostat and barostat relaxation times were 0.75 and 2 ps, respectively. We set the cutoff of both dispersion interactions and the real part of the electrostatic interactions at 12 Å. Long-range Coulombic interactions were treated with the Particle-Mesh Ewald (PME) solver.83,84 We kept the O–H bond length (0.9572 Å) and H–O–H angle (104.52°) values constant with the LINCS algorithm.85 With this model, we reached a computational performance of 40 ns of simulation time per wall clock day for a system of about 10 000 water molecules running with Intel(R) Xeon(R) Platinum 8368Q CPU @ 2.60 GHz, using 32 CPUs in parallel.

Seeding is a method that consists of using classical nucleation theory (CNT) in combination with MD simulations.34,35 More specifically, we use the NVT seeding approach,86 in which a cluster (in this case a bubble) close in size to the critical one is artificially inserted into the system, then spontaneously equilibrated into the critical size and tracked along time. With this approach, a critical bubble can be characterized for long timescales because the maximum in free energy barrier in a nucleation process represents a minimum in the Helmholtz free energy landscape in the canonical ensemble.54,87 Therefore, more precise measurements can be made compared to seeding in the NPT ensemble, where the bubble will rapidly either shrink or grow.86 This method is suitable to measure nucleation rates along isotherms since the pressure at which the cluster is critical cannot be known a priori and is obtained from the simulations.

CNT88,89 is a theoretical framework that describes nucleation processes under saturation conditions. It can be used to obtain the free energy barrier, interfacial free energy, and nucleation rate. The limitations of quantitatively characterizing nucleation rates using CNT are due to assumptions inherent in the theory.90–95 Despite these potential limitations, multiple studies position CNT as a powerful tool to estimate free energy barriers and nucleation rates for phase transitions,34,35,96–106 including cavitation.53,86,107,108 According to CNT,88,109 the nucleation rate (J) can be computed as
J=ρl2γπmexpΔG*kBT,
(3)
where ρl represents the density of the liquid phase, γ the liquid–vapor surface tension, m the mass of water, ΔG* the free energy barrier for nucleation, kB the Boltzmann’s constant, and T the temperature. Within the CNT framework, the free energy barrier is obtained as
ΔG*=43γπRc,
(4)
where Rc is the critical bubble radius. Additionally, we obtain the interfacial free energy from Laplace’s equation as
γ=RcΔP2,
(5)
where ΔP is the pressure difference between the vapor and liquid phases. This approach provides more reliable estimations than assuming the capillarity approximation (i.e., inserting the surface tension at planar interface and coexistence conditions into the CNT).86,108,110,111 Combining Eqs. (3)(5), we reach the final equation for J,
J=ρlRcΔPπmexp4πRc2ΔP3kBT.
(6)
To summarize, in order to compute J, we require the pressure difference between the liquid and the vapor phases and the critical radius of the bubble at the corresponding thermodynamic conditions of P and T. Although the difference in pressure can be computed in principle,86,110 in this study the pressure inside the bubbles is ∼0;48 therefore, ΔP can be easily estimated as −Pliq, which is directly obtained through the virial expression. We additionally confirmed that the virial pressure obtained in the system containing a bubble matches that of the bulk liquid (Fig. S1), as recently shown.112 
Finally, the critical radius, Rc, was obtained employing a local order parameter. Although multiple parameters have been proposed to track the size of a simulated bubble,53,55,58,59,113,114 here, we adopted the “equidensity” criterion,115 which was shown to provide the surface of tension radius (i.e., the radius that, when inserted into Laplace’s equation, provides a consistent value of γ) for the Lennard-Jones system.86,108 As illustrated in Fig. 1(a), the center of the bubble is first calculated through the minimum in the density profile along the three cartesian directions. Afterward, a radial density profile from the bubble center is computed, in which the critical radius (Rc) corresponds to the point in which the density equals the average of the liquid and vapor densities [Fig. 1(b)]. This point is found via nonlinear fitting to the equation
ρ(r)=ρl+ρv2+ρlρv2tanhrRcα,
(7)
where r is the distance from the bubble center and α is a fitting parameter. We confirmed that, as assumed by CNT, the bubbles have a close-to-spherical shape (Fig. S2).
FIG. 1.

(a) Density profiles along the three cartesian directions. Vertical dashed lines depict the location of the minimum density, which corresponds to the center of the bubble in each direction. In these density profiles, we averaged the density of each point with two other neighboring points in order to make the curves smoother. (b) Radial density profile calculated from the center of the bubble. The blue curve indicates the density calculated at each distance while the black dashed curve is the fit of the blue curve to Eq. (7), from which we obtained the critical radius Rc.

FIG. 1.

(a) Density profiles along the three cartesian directions. Vertical dashed lines depict the location of the minimum density, which corresponds to the center of the bubble in each direction. In these density profiles, we averaged the density of each point with two other neighboring points in order to make the curves smoother. (b) Radial density profile calculated from the center of the bubble. The blue curve indicates the density calculated at each distance while the black dashed curve is the fit of the blue curve to Eq. (7), from which we obtained the critical radius Rc.

Close modal

It is important to note that all DPMD data shown in this work are shifted by 40 K, so that the simulations for a given temperature have been performed at 40 K higher than the reported one in the presented figures. Similar shifts in temperature have been performed for AIMD simulations using SCAN116 and in other works using SCAN-based ML models.117–119 The rationale for the shift was originally attributed to nuclear quantum effects, but it is likely mainly due to the limitations of the density functional itself. SCAN is known to overestimate the strength of the hydrogen bond.117 In this work, the shift in temperature was adjusted by calculating the mean square error between the coexistence vapor densities predicted by the model and the experimental ones, considering different values for the shift. The value for the shift was iteratively modified until we obtained the one that gave the minimum mean square error, which was 40 K. In all plots and tables that follow, the results of the DPMD model have this shift already applied. We note that this 40 K shift is applicable only to the liquid–vapor equilibrium regime, and different model adjustments should be made for a precise study of the different solid phases.

To test the DPMD model in the liquid–vapor regime, we begin by computing the phase diagram. We do so using simulations in the canonical ensemble, in which both the liquid and vapor phases coexist. In Fig. 2(a) (inset), we show a snapshot of a typical DC simulation box. In Fig. 2(a), we plot the temperature against density phase diagram of the DPMD model (green points), where the filled points represent the densities directly obtained from DC simulations (2 at each temperature). We include results for the TIP4P/2005 model10,48 (blue circles), as well as experimental data47 (black line). We can observe that for the DPMD model, the density of the liquid branch is slightly higher than the experimental results at low temperatures (<500 K) but matches well at higher ones. The critical point is estimated through the universal scaling law of coexistence densities near a critical point121 and the law of rectilinear diameters as follows:122 
ρl(T)ρv(T)3.06=d1TTc
(8)
and
(ρl(T)+ρv(T))/2=ρc+s2(TcT),
(9)
where ρl and ρv refer to the coexisting densities of the liquid and vapor phases, respectively, ρc is the critical density, Tc is the critical temperature, and d and s2 are fitting parameters. The critical temperature obtained for the DPMD model is of 632.6.6 K, which is lower than the experimental one by 14.5 K.
FIG. 2.

(a) Phase diagram in the Tρ plane of the DPMD model (green), TIP4P/2005 model (blue),48 and experimental measurements of water (black).47 Filled circles represent the vapor and liquid densities, estimated from the averaged bulk densities from DC simulations. Inset: Snapshot of a DC simulation performed at 385 K with the DPMD model, rendered making use of Ovito software.120 (b) Liquid–vapor interfacial free energy (γ) as a function of temperature for the DPMD model (green) and comparison with the TIP4P/2005 model (blue)49 and experimental values (black).47 (c) Vapor saturation pressure (Psat) as a function of temperature for the DPMD model (green), TIP4P/2005 (blue)48 and comparison with experimental values (black).47 Inset: Enthalpy of vaporization (ΔHvap) as a function of temperature for DPMD, TIP4P/200548 and experimental values47 as indicated in the legend.

FIG. 2.

(a) Phase diagram in the Tρ plane of the DPMD model (green), TIP4P/2005 model (blue),48 and experimental measurements of water (black).47 Filled circles represent the vapor and liquid densities, estimated from the averaged bulk densities from DC simulations. Inset: Snapshot of a DC simulation performed at 385 K with the DPMD model, rendered making use of Ovito software.120 (b) Liquid–vapor interfacial free energy (γ) as a function of temperature for the DPMD model (green) and comparison with the TIP4P/2005 model (blue)49 and experimental values (black).47 (c) Vapor saturation pressure (Psat) as a function of temperature for the DPMD model (green), TIP4P/2005 (blue)48 and comparison with experimental values (black).47 Inset: Enthalpy of vaporization (ΔHvap) as a function of temperature for DPMD, TIP4P/200548 and experimental values47 as indicated in the legend.

Close modal

Moreover, we compute the liquid–vapor surface tension for the DPMD model at different temperatures from the DC simulations. This quantity can be directly estimated according to Eq. (2). As can be seen in Fig. 2(b), the DPMD model provides lower values of γ than both the TIP4P/2005 model49 and experimental measurements.47 From DC simulations, we also calculate the saturation pressure (Psat) as a function of temperature. Psat is obtained as the component of the pressure tensor normal to the interface. We plot it in Fig. 2(c), compared to the values obtained from TIP4P/2005 (blue)48 and experiments (black).47 This quantity closely matches with experimental measurements, while the TIP4P/2005 model underestimates it, which is a natural consequence of the way the DPMD model was shifted to match the vapor densities. We note that the temperature shift was applied in order to obtain a better match of the vapor phase behavior; nonetheless, this shift affects negatively on the surface tension prediction. With no temperature shift, the surface tension of the DPMD model matches the experimental one at temperatures above 450 K and only underestimates it by 5% at T < 450 K (see Fig. S4). By analogy to the “short-blanket” dilemma recently discussed,123 here, we choose to prioritize the description of the Tρ phase diagram and Psat, disfavoring the γ estimates.

We estimated the enthalpy of vaporization (ΔHvap) from independent bulk simulations of both liquid and vapor phases. For this, we perform canonical simulations at the equilibrium density of the given phase, which was previously obtained from DC simulations. From each simulation, the enthalpy is directly obtained as H = UPV, where U is the internal energy. The enthalpy of vaporization is simply calculated as ΔHvap = HvaporHliquid for every temperature. The values of ΔHvap are plotted against temperature in Fig. 2(c) (inset), along with values from TIP4P/200548 and experiments.47 This quantity slightly deviates from the experimental values at low temperatures (<550 K) but matches at higher temperatures. Similarly, the TIP4P/2005 model matches the experimental trend at T > 500 K and overestimates it at lower T.

In summary, the DPMD model describes the liquid–vapor coexistence properties after applying the temperature shift of 40 K reasonably well. Some discrepancies may arise from the fact that this model has been trained on solid and liquid data only,42 but the main source of differences from experimental data are limitations in accuracy for the SCAN density functions used to train the model. The biggest difference with experimental values is found in the surface tension, which is underestimated by ∼20%. Nonetheless, since the training was not focused on liquid–vapor configurations, we do not expect the SCAN functional to necessarily provide an accurate description of γ. Other than this discrepancy, the phase diagram, saturation pressure, and enthalpy of vaporization are well described using the DPMD model. In addition to the plots in Fig. 2, we provide the equilibrium data of the DPMD model in Table I.

TABLE I.

Data for the liquid (ρl) and vapor (ρv) densities, saturation pressure (Psat), and surface tension (γ) as a function of temperature for the DPMD model. The numbers in parenthesis depict the uncertainty of our measurements and apply to the numeral left of themselves, for instance, 41(3) stands for 41 ± 3.

T (K)ρl (g cm−3)ρv (g cm−3)Psat (bar)γ (mN m−1)
385 0.9697(5) 0.0011(3) 2.3(3) 41(4) 
410 0.9484(8) 0.0022(4) 4.9(5) 35(5) 
435 0.9242(6) 0.0037(6) 8.3(9) 32(3) 
460 0.8969(8) 0.0067(9) 15(1) 27(3) 
485 0.865(3) 0.012(4) 25(3) 23(4) 
510 0.831(2) 0.016(2) 35(4) 16(4) 
535 0.791(1) 0.024(1) 48(7) 15(5) 
560 0.740(6) 0.039(5) 71(6) 13(4) 
585 0.684(8) 0.058(4) 98(15) 7(2) 
597 0.62(1) 0.079(5) 114(5) 3(3) 
T (K)ρl (g cm−3)ρv (g cm−3)Psat (bar)γ (mN m−1)
385 0.9697(5) 0.0011(3) 2.3(3) 41(4) 
410 0.9484(8) 0.0022(4) 4.9(5) 35(5) 
435 0.9242(6) 0.0037(6) 8.3(9) 32(3) 
460 0.8969(8) 0.0067(9) 15(1) 27(3) 
485 0.865(3) 0.012(4) 25(3) 23(4) 
510 0.831(2) 0.016(2) 35(4) 16(4) 
535 0.791(1) 0.024(1) 48(7) 15(5) 
560 0.740(6) 0.039(5) 71(6) 13(4) 
585 0.684(8) 0.058(4) 98(15) 7(2) 
597 0.62(1) 0.079(5) 114(5) 3(3) 

After establishing the equilibrium properties of the DPMD water model, we proceeded to investigate its cavitation. Although some experimental studies of water cavitation have been conducted,63–67 it is difficult to establish a direct comparison due to the lack of measurements of the nucleation rate. Menzl et al.53 performed a nucleation study utilizing Umbrella Sampling (US) calculations124 for the TIP4P/2005 model, in which the nucleation free energy barrier and the nucleation rate were reported, without comparisons to experimental data. Here, we employ the NVT seeding technique at the same temperature (296.4 K) as in Ref. 53 for both the TIP4P/2005 (to establish the validity of our methods) and DPMD models (to provide new data for this ab initio-based model).

We prepared various systems in which we artificially generated a cavity of a given size, starting from a bulk liquid configuration. As detailed in Sec. II C, the system spontaneously evolves and equilibrates into a state in which there is a critical bubble that remains stable over time, due to the fact that in the canonical ensemble, a critical bubble represents a local minimum in the Helmholtz free energy landscape.54,87 Once equilibrated, we measured the system pressure by means of the virial expression,125 which corresponds to the liquid phase pressure.86 To track the critical radius, we made use of our order parameter (see Sec. II C). We repeated this process for each configuration and then averaged over more than 500 independent radial density profiles for the calculation of the radius.

We used our data for ΔP and Rc along with Eq. (6) to compute J, which is plotted in Fig. 3(a) (blue and green squares for TIP4P/2005 and DPMD respectively) against P. It can be seen that there is agreement within the simulation uncertainties between the US and seeding simulations for the TIP4P/2005 model. We also include a continuous line for each model, which represents a fit to the CNT equation, in which we linearly fit γ against P and insert values from such fit to solve Eq. (6). The uncertainty is estimated from the standard deviation of the radius between different independent configurations.

FIG. 3.

(a) Nucleation rate (J) as a function of pressure (P) for TIP4P/2005 (blue) and DPMD (green) cavitation at 296.4 K, including data from Ref. 53. Continuous lines are obtained by linearly fitting the surface tension (γ) as a function of pressure and then inserting such γ into Eqs. (5) and (6). The shaded region is obtained in the same way but making use of the upper and lower bounds of the surface tension error and, therefore, represents the error limits in J. (b) Free energy barrier for bubble nucleation as a function of pressure at 296.4 K for TIP4P/2005 (blue) and DPMD (green) models, estimated from CNT [Eq. (4)].

FIG. 3.

(a) Nucleation rate (J) as a function of pressure (P) for TIP4P/2005 (blue) and DPMD (green) cavitation at 296.4 K, including data from Ref. 53. Continuous lines are obtained by linearly fitting the surface tension (γ) as a function of pressure and then inserting such γ into Eqs. (5) and (6). The shaded region is obtained in the same way but making use of the upper and lower bounds of the surface tension error and, therefore, represents the error limits in J. (b) Free energy barrier for bubble nucleation as a function of pressure at 296.4 K for TIP4P/2005 (blue) and DPMD (green) models, estimated from CNT [Eq. (4)].

Close modal
We computed J at higher superstreching conditions (green and blue diamonds in Fig. 3) through “brute force” simulations. In these simulations, we observed the metastable bulk liquid under high superstreching conditions in the NPT ensemble for a sufficient time before spontaneous cavitation takes place. Then, J is calculated as
J=1tV,
(10)
where ⟨t⟩ is the average time required for cavitation to occur and V is the volume of the metastable liquid phase. The onset of cavitation can easily be identified with a sudden and sharp change in properties such as the simulation box volume or the potential energy. From these simulations, we obtain J = 2.52 · 10−6 ps−1 nm−3 at P = −150 MPa for the DPMD model and J = 2.44 · 10−5 ps−1 nm−3 and P = −200 MPa for the TIP4P/2005 potential. These results are also shown in Fig. 3(a) and match with the trend of US and seeding simulations. This result, in addition to the agreement with the US calculations from Menzl et al., provides confidence in the validity of the results obtained using CNT.

In addition to the nucleation rate, we also obtained the free energy barrier (ΔG*), which can be estimated from Eq. (4). This quantity is the main output from US simulations.53 In Fig. 3(b), we compare the calculated free energy barriers for the DPMD (green) and TIP4P/2005 (blue) models, also including data from Ref. 53. As expected, we find good agreement between seeding and US calculations as for the nucleation rates.

It can be seen in Fig. 3(a) that the DPMD model returns nucleation rates many orders of magnitude greater than the TIP4P/2005 potential, outside the uncertainty bounds, despite the close resemblance of the phase diagrams from the two models [Fig. 2(a)]. This difference is also present for the free energy barrier [Fig. 3(b)], with the TIP4P/2005 model possessing a higher free energy barrier. This is likely the crucial factor behind its lower nucleation rate. In order to understand the differences between the two models, we also compared the change of the surface tension with curvature for both models: In Fig. 4, we plot γ as a function of P, where filled points correspond to the value obtained at the coexistence pressure from DC simulations, while the empty points depict the value of γ obtained through Laplace’s equation [Eq. (5)] from our seeding simulations. From this analysis, we observe that the surface tension is significantly lower for the DPMD model not only under coexistence conditions but also in the cavitation regime. This directly points toward the surface tension being the decisive factor behind the quantitative difference in J and ΔG* between the DPMD and TIP4P/2005 models. In Table II, we detail the different quantities playing a role in Eqs. (2)(4). Since the kinetic prefactor in the calculation of J is of the same order of magnitude in all cases, we can conclude that the different nucleation rates between the TIP4P/2005 and DPMD models arises from a quantitative difference in the surface tension, which is lower for the DPMD.

FIG. 4.

Liquid–vapor surface tension (γ) as a function of pressure at 296.4 K for TIP4P/2005 (blue) and DPMD (green). Empty circles are obtained from cavitation simulations (and therefore with a curved interface) by means of Laplace’s equation (see Sec. II C), while filled points were estimated from DC simulations as in Fig. 2.

FIG. 4.

Liquid–vapor surface tension (γ) as a function of pressure at 296.4 K for TIP4P/2005 (blue) and DPMD (green). Empty circles are obtained from cavitation simulations (and therefore with a curved interface) by means of Laplace’s equation (see Sec. II C), while filled points were estimated from DC simulations as in Fig. 2.

Close modal
TABLE II.

NVT seeding data for the DPMD and TIP4P/2005 models, including the nucleation pressure (P), the critical radius (Rc), the liquid density (ρl), the total number of water molecules in the system (NT), the surface tension (γ), the free energy barrier (ΔG*), and the logarithm of the nucleation rate (log10J).

P (MPa)Rc (nm)ρl (g cm−3)NTγ (mN m−1)ΔG*/kBTlog10 (J/(ps−1 nm−3))
DPMD 
−81.9 1.25 0.965 7 527 51.0 71.5 −29.5 
−69.1 1.50 0.972 7 386 51.9 105.7 −44.4 
−60.3 1.74 0.977 7 161 52.6 144.4 −61.2 
−54.4 1.96 0.980 6 888 53.2 183.5 −78.2 
−47.1 2.29 0.984 10 324 53.8 253.7 −108.6 
−43.7 2.50 0.986 9 825 54.7 309.2 −132.7 
TIP4P/2005 
−96.9 1.24 0.953 6 855 59.9 93.6 −39.1 
−87.6 1.38 0.956 6 796 60.3 116.8 −49.2 
−77.0 1.59 0.960 23 814 61.2 158.5 −67.3 
−73.0 1.68 0.961 11 385 61.5 178.2 −75.8 
−61.9 2.01 0.965 11 006 62.0 255.3 −109.3 
−54.9 2.29 0.968 10 570 62.8 337.3 −144.9 
P (MPa)Rc (nm)ρl (g cm−3)NTγ (mN m−1)ΔG*/kBTlog10 (J/(ps−1 nm−3))
DPMD 
−81.9 1.25 0.965 7 527 51.0 71.5 −29.5 
−69.1 1.50 0.972 7 386 51.9 105.7 −44.4 
−60.3 1.74 0.977 7 161 52.6 144.4 −61.2 
−54.4 1.96 0.980 6 888 53.2 183.5 −78.2 
−47.1 2.29 0.984 10 324 53.8 253.7 −108.6 
−43.7 2.50 0.986 9 825 54.7 309.2 −132.7 
TIP4P/2005 
−96.9 1.24 0.953 6 855 59.9 93.6 −39.1 
−87.6 1.38 0.956 6 796 60.3 116.8 −49.2 
−77.0 1.59 0.960 23 814 61.2 158.5 −67.3 
−73.0 1.68 0.961 11 385 61.5 178.2 −75.8 
−61.9 2.01 0.965 11 006 62.0 255.3 −109.3 
−54.9 2.29 0.968 10 570 62.8 337.3 −144.9 
Another quantity we can extract from our simulations is the Tolman length, which describes the deviation of the surface tension with respect to its value at the planar interface and, therefore, the coexistence conditions. The Tolman length can also be defined as the deviation of the surface of tension from the equimolar dividing surface. In 1949, Tolman showed that the change in surface tension with curvature follows the equation68,
γ=γ01+2δRc,
(11)
where γ0 is the value of the surface tension under coexistence conditions and δ is the Tolman length. This quantity has been extensively studied for Lennard-Jones particles,62,126–130 hard spheres,87 and other systems.126,127,131,132 Here, we make use of Tolman’s expression and compute δ for both the DPMD and TIP4P/2005 models at 296.4 K. We fit our surface tension and critical radius data (obtained from the NVT seeding simulations) to Eq. (11), performing a nonlinear regression. We choose to have γ0 and δ as fitting parameters, despite having estimated the former from the DC simulations. We took this approach in order to corroborate the value of γ obtained from both approaches.

In Fig. 5, we plot the surface tension against the inverse of the critical radius for the DPMD (green filled points) and TIP4P/2005 (blue filled points) models. From nonlinear regression, we obtain δ = (0.091 ± 0.008) nm for the DPMD model and δ = (0.070 ± 0.004) nm for the TIP4P/2005 model. Both values are positive as expected since the surface tension decreases for smaller bubbles. In Fig. 5, we also show with dashed lines how the surface tension changes against the inverse of the critical radius according to Eq. (11). From the fit, we also obtain values for the surface tension at planar interface of 58 ± 2 and 67 ± 3 mN m−1 for the DPMD and TIP4P/2005 models, respectively. These values are in agreement, within statistical uncertainties, with those obtained from DC simulations (54 and 7049 mN m−1 for DPMD and TIP4P/2005, respectively). The uncertainty of δ is simply determined by the error in the nonlinear fit, while the uncertainty in γ0 is the sum of the errors coming from NVT seeding simulations (see Sec. III B) and the error in the fit to Eq. (11).

FIG. 5.

Surface tension γ against the inverse of the critical radius (Rc1). Filled points represent the surface tension obtained from Laplace’s equation [Eq. (5)]. The dashed lines indicate the change in surface tension according to Tolman’s equation [Eq. (11)], where the parameters δ and γ0 were obtained through nonlinear regression. We plot the values of surface tension at planar interface (Rc1 = 0) obtained from the fit with empty points. The surface tension values obtained from DC are also represented here with filled diamonds, also at P = 0.

FIG. 5.

Surface tension γ against the inverse of the critical radius (Rc1). Filled points represent the surface tension obtained from Laplace’s equation [Eq. (5)]. The dashed lines indicate the change in surface tension according to Tolman’s equation [Eq. (11)], where the parameters δ and γ0 were obtained through nonlinear regression. We plot the values of surface tension at planar interface (Rc1 = 0) obtained from the fit with empty points. The surface tension values obtained from DC are also represented here with filled diamonds, also at P = 0.

Close modal

Our results are also in good agreement with previous simulation results54,60,133 (at different temperatures: 0.199 nm at 300 K,54 0.09 nm at 250 K,133 0.18 nm at 350 K,133 and 0.1 nm at 450 and 550 K60), which indicate that the Tolman length is positive for water bubbles, in contrast to the negative sign in water droplets (i.e., condensation),56,57 where the surface tension increases with curvature. We can establish direct comparison with the calculations from Menzl et al.,53 which were performed with TIP4P/2005 at the same temperature (296.4 K). They obtained values of δ = 0.195 nm and γ0 = 82.79 mN m−1, which while having the same order of magnitude, moderately disagree with our calculations of γ0 and δ. In the case of Ref. 53, the employed local order parameter does not necessarily identify an accurate value of the radius but is instead used to bias the sampling of the configurational space for US simulations and may therefore not represent accurate values of δ and γ0. As mentioned before, the good agreement in the nucleation free energy barrier between US and seeding calculations is a good sign. In our case, a good indicator for the calculation of the Tolman length, although not definitive, is the fact that the surface tension obtained from the nonlinear regression to Eq. (11) provides a value of γ0 that matches within the uncertainty the surface tension obtained from DC simulations, as aforementioned.

Finally, we examined the organization of the liquid–vapor interface in our simulations for both planar and curved interfaces. The water–air planar interface has been widely studied in the past134 with IR vibrational spectra experiments135–138 and ab initio139–142 and MD simulations,143–145 all of which generally agree regarding the orientation of water molecules near the interface. Interfacial molecules closer to the vapor phase tend to orient the O–H bond parallel to the surface normal vector and expose the H atom, resembling the (1000) crystal face of ice Ih (only in the dimension perpendicular to the surface).145 The (1000) ice Ih–liquid is the direction with lower interfacial free energy of the different solid–liquid interfaces for this polymorph;146 therefore, in the liquid–vapor interface, a similar arrangement may also reduce the surface tension, apart from maximizing the enthalpic gain of the more exposed interfacial molecules to the vapor phase. To perform this analysis, we follow the same approach as in Refs. 142 and 145 where, once the interfacial region has been identified, the angle formed by the O–H bonds and the vector normal to the interface (θH) is computed. Fan et al.145 and Vassilev et al.142 identified two distinct layers in non-curved interfaces: an external layer, which comprehends the region in which the density is between 5% and 50% of the bulk liquid, and an internal layer, where the density ranges between 50% and 95% of the bulk liquid density. Fan et al.145 measured the orientation distributions for the TIP3P, TIP4P-EW, TIP5P, and SPC/E water models at 300 K, and for comparison purposes, we perform the same analysis with the DPMD and TIP4P/2005 models and obtain results also for curved interfaces.

In Fig. 6(a), we depict how the different interfacial regions and θH angle are identified for both planar (top) and curved (bottom) interfaces. There, the dashed black line indicates the liquid–vapor interface as defined by our order parameter (see Sec. II C and the supplementary material). In the zoomed image, the internal and external layers of the interface are labeled and highlighted in red and black as well as the bulk liquid and vapor phases in green and white, respectively. The two arrows indicate the normal vector to the interface and the O–H bond vector and the angle between these two vectors defines θH. In Figs. 6(b) and 6(c), we plot the θH distributions normalized by the random distribution sin(θH). Consistent with previous calculations,141,142,144,145 the molecules expose one hydrogen atom to the vapor phase in the external layer, as can be seen in Fig. 6(b) (top) for the TIP4P/2005 model, where we show the probability distribution of the angle θH for the different interfacial regions. Smaller angles are more probable in the external layer than in the bulk, where all directions are equally probable (flat green line). As discussed in Ref. 145, the internal layer of the interface also displays preferential orientations, in order to maximize the interactions with the structure created in the external layer, in an arrangement that leads to a maximum near θH ∼ 80°. It must be noted that the preferential molecular orientations are not fixed but rather transient. The distributions presented in Fig. 6(b) (top) are also consistent with the results presented in the work of Fan et al.145 for other rigid semiempirical models.

FIG. 6.

(a) Top: Slab snapshot (including only a reduced amount of molecules for better visualization) of a planar interface. Bottom: Snapshot for a curved interface. (b) Normalized histograms of θH at the external region of the interface, the internal one, and bulk water for the TIP4P/2005 model for planar (top) and curved (bottom) interfaces. (c) Same as in (b) but for the DPMD model.

FIG. 6.

(a) Top: Slab snapshot (including only a reduced amount of molecules for better visualization) of a planar interface. Bottom: Snapshot for a curved interface. (b) Normalized histograms of θH at the external region of the interface, the internal one, and bulk water for the TIP4P/2005 model for planar (top) and curved (bottom) interfaces. (c) Same as in (b) but for the DPMD model.

Close modal

We extended this analysis to curved interfaces using our NVT seeding simulations. In curved interfaces, the vector normal to the interface points at the bubble center. In Fig. 6(b) (bottom), we show the probability distribution of θH for the three regions and observe how these resemble those of planar interface. A significant loss of ordering is observed for the curved interface since the peaks at ∼15° and ∼115° in the external layer and at ∼80° in the internal layer are less prominent compared to the corresponding peaks for the flat interface. This could be a consequence of the molecular steric hindrance caused by the curved geometry of the interfaces in the bubbles.

We also evaluated the θH distribution for the DPMD model. In Fig. 6(c) (top), we present the distribution of the three interfacial regions for the planar interface at 296.4 K. There is an obvious preferential orientation toward small angles in the external layer, once again indicating the preference of molecules to expose one hydrogen atom to the vapor phase. The distribution is remarkably more prominent for the DPMD model relative to TIP4P/2005 although both models result in similar molecular arrangements. In Ref. 142, this same structure was found in ab initio calculations using PW91 functional, confirming that the DPMD model based on the SCAN functional yields similar results to those obtained from ab initio MD simulations using different density functionals.

Looking at the DPMD model distributions for curved interfaces [Fig. 6(c) (bottom)], we still find a preferential ordering toward small angles, although it is much less pronounced than in the case of planar interfaces, in agreement with TIP4P/2005 simulations. The bubbles generated by both models shown here (TIP4P/2005 and DPMD) have comparable sizes; however, the effect of curvature is more dramatic in the case of the DPMD model. Other bubble size distributions are reported in the supplementary material, where we observe that there is a slight change in the distributions, with bigger bubbles resulting in curves more similar to those found at coexistence (Fig. S3).

In this work, we explored the liquid–vapor phase behavior of a Deep Potential model based on ab initio energies and forces. This model was derived from the SCAN approximation of density functional theory.42 The model has been shown to closely follow the experimentally determined phase diagram for the different ice phases (excepting ices III and XV)42 and to have a liquid–liquid phase transition in the supercooled regime.118 We computed the phase diagram via DC simulations and adjusted our vapor equilibrium densities to the experimental values, resulting in a shift of 40 K in the DPMD model. Once the model was tested and shifted, we compare its equilibrium properties with experimental data and the TIP4P/2005 model,10 one of the most benchmarked classical models for water.48,49,53,143 The surface tension of the DPMD model is lower than the one obtained from TIP4P/2005 and experiments. Nonetheless, by construction the DPMD model provides accurate results for other properties such as the vapor saturation pressure for a wide range of temperatures, between 300 and 600 K (Fig. 2). Overall, once the temperature is shifted, the model reproduces most of the properties of liquid–vapor equilibrium, despite not having this regime included in its training. This result is especially remarkable since this is the first Deep Potential-based model to provide liquid–vapor properties in a computationally affordable time, with the primary drawback being ∼20% deviation in the surface tension.

Moreover, we study bubble nucleation in the cavitation regime. We make use of the NVT seeding method, a rare event technique already shown to be successful in determining the nucleation free energy barrier and the nucleation rate for simpler systems86,108 in cavitation events. We first performed NVT seeding calculations at 296.4 K, a temperature at which previous data from Umbrella Sampling are available for the TIP4P/2005 model, and we confirm that the seeding technique provides consistent results with those from the work of Menzl et al.53 The DPMD water model provides higher nucleation rates than the TIP4P/2005 model under the same stretching conditions. We show that this quantitative difference can be explained by the difference in surface tension between models, which persists for curved interfaces (Fig. 4). Our results highlight once more the relevance of surface tension and its change with curvature to critically control nucleation events.53,133 We could have obtained closer agreement between the TIP4P/2005 nucleation rates with those from the DPMD model if we did not apply the temperature shift; nonetheless, we prioritized adjusting the model to obtain better equilibrium densities, in line with prior studies of water properties using the SCAN-derived DPMD model.117–119 

Furthermore, we provide an estimate of the Tolman length by performing a nonlinear regression to the relevant expression,68 which describes how the surface tension changes with curvature [Eq. (11)]. Using data from our NVT seeding simulations, we obtain estimates of δ = (0.091 ± 0.008) nm for the DPMD model and δ = (0.070 ± 0.004) nm for the TIP4P/2005 model, both at 296.4 K. This confirms that a Deep Potential-based model also predicts a positive sign of the Tolman length, confirming previous results showing a decrease of the surface tension with curvature for the case of water bubbles.53,54,133

Finally, we studied the orientation of water molecules in the interface, corroborating previous studies that have indicated that the molecules closer to the vapor phase have a preference so as to expose a hydrogen atom facing the vapor.134–140,143–145 We quantify this behavior by measuring the angle formed between the normal vector to the interface and the O–H bond. We find that a preferential molecular orientation appears for both the TIP4P/2005 and DPMD models. We also confirm that this phenomenon also takes place in the curved interface of water bubbles, although the possibility to orient more O–H bonds toward the vapor is diminished for curved interfaces due to the increasing curvature of the bubbles.

Overall, this study confirms that machine-learning ab initio-based models are viable for prediction of equilibrium as well as dynamic properties that require large system sizes and long sampling timescales. The computational cost of ab initio-based models in long-scale Molecular Dynamics is now affordable, being only an order of magnitude greater than that for empirical potentials, thanks to recent improvements such as the compressed Deep Potential modeling scheme.72 DPMD’s higher computational cost can be explained by the requirement of a smaller time step to take into account hydrogen vibrations and its many-body nature and more complex functional form when compared to pairwise classical potentials. Further work and models trained with more liquid and vapor data, as well as liquid–vapor interfaces,147 will only improve the already existing models, which will gradually be better in describing the real behavior of water.

See the supplementary material for additional information on the ρ against P equation of state at 296.4 K for the DPMD and TIP4P/2005 models; the bubble sphericity analysis; and the orientational analysis of bubbles of different sizes.

I.S.-B. acknowledges funding support provided by Derek Brewer scholarship of Emmanuel College and EPSRC Doctoral Training Program studentship, Grant No. EP/T517847/1. JR.E. acknowledges funding support provided by the Roger Ekins Research Fellowship of Emmanuel College, Oppenheimer Research Fellowship of the University of Cambridge and the Ramon y Cajal fellowship (Grant No. RYC2021-030937-I). M.C.M and A.Z.P acknowledge the support from the “Chemistry in Solution and at Interfaces” (CSI) Center funded by the U.S. Department of Energy Award Nos. DE-SC001934 and DE-SC0002128. M.C.M also acknowledges funding from the High Meadows Environmental Institute at Princeton University through the Mary and Randall Hack ’69 Graduate Award for Water and the Environment. Computational resources were provided by the Princeton Research Computing at Princeton University, which is a consortium of groups led by the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology’s Research Computing.

The authors have no conflicts to disclose.

Ignacio Sanchez-Burgos: Conceptualization (equal); Investigation (equal); Validation (equal); Writing – original draft (lead); Writing – review & editing (lead). Maria Carolina Muniz: Conceptualization (equal); Investigation (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Jorge R. Espinosa: Writing – review & editing (equal). Athanassios Z. Panagiotopoulos: Conceptualization (lead); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material. Example files of simulations employing the DPMD potential are available in the Princeton DataSpace repository https://doi.org/10.34770/ms7d-wm45.

1.
A.
Chanda
and
V. V.
Fokin
, “
Organic synthesis ‘on water
,’”
Chem. Rev.
109
(
2
),
725
748
(
2009
).
2.
M.
Chaplin
, “
Do we underestimate the importance of water in cell biology?
,”
Nat. Rev. Mol. Cell Biol.
7
(
11
),
861
866
(
2006
).
3.
E.
Lester
,
P.
Blood
,
J.
Denyer
,
D.
Giddings
,
B.
Azzopardi
, and
M.
Poliakoff
, “
Reaction engineering: The supercritical water hydrothermal synthesis of nano-particles
,”
J. Supercrit. Fluids
37
(
2
),
209
214
(
2006
).
4.
A.
Ponomarenko
,
O.
Vincent
,
A.
Pietriga
,
H.
Cochard
,
É.
Badel
, and
P.
Marmottant
, “
Ultrasonic emissions reveal individual cavitation bubbles in water-stressed wood
,”
J. R. Soc. Interface
11
(
99
),
20140480
(
2014
).
5.
T. D.
Wheeler
and
A. D.
Stroock
, “
The transpiration of water at negative pressures in a synthetic tree
,”
Nature
455
(
7210
),
208
212
(
2008
).
6.
U.
Adhikari
,
A.
Goliaei
, and
M. L.
Berkowitz
, “
Mechanism of membrane poration by shock wave induced nanobubble collapse: A molecular dynamics study
,”
J. Phys. Chem. B
119
(
20
),
6225
6234
(
2015
).
7.
D.
Yu
,
B.
Liu
, and
B.
Wang
, “
The effect of ultrasonic waves on the nucleation of pure water and degassed water
,”
Ultrason. Sonochem.
19
(
3
),
459
463
(
2012
).
8.
P.
Kumar
and
R. P.
Saini
, “
Study of cavitation in hydro turbines—A review
,”
Renewable Sustainable Energy Rev.
14
(
1
),
374
383
(
2010
).
9.
J.
Ma
,
A.
Michaelides
,
D.
Alfe
,
L.
Schimka
,
G.
Kresse
, and
E.
Wang
, “
Adsorption and diffusion of water on graphene from first principles
,”
Phys. Rev. B
84
(
3
),
033402
(
2011
).
10.
J. L. F.
Abascal
and
C.
Vega
, “
A general purpose model for the condensed phases of water: TIP4P/2005
,”
J. Chem. Phys.
123
(
23
),
234505
(
2005
).
11.
J. L. F.
Abascal
,
E.
Sanz
,
R.
García Fernández
, and
C.
Vega
, “
A potential model for the study of ices and amorphous water: TIP4P/ice
,”
J. Chem. Phys.
122
(
23
),
234511
(
2005
).
12.
M. W.
Mahoney
and
W. L.
Jorgensen
, “
Quantum, intramolecular flexibility, and polarizability effects on the reproduction of the density anomaly of liquid water by simple potential functions
,”
J. Chem. Phys.
115
(
23
),
10758
10768
(
2001
).
13.
W. L.
Jorgensen
,
J.
Chandrasekhar
,
J. D.
Madura
,
R. W.
Impey
, and
M. L.
Klein
, “
Comparison of simple potential functions for simulating liquid water
,”
J. Chem. Phys.
79
(
2
),
926
935
(
1983
).
14.
H. W.
Horn
,
W. C.
Swope
, and
J. W.
Pitera
, “
Characterization of the TIP4P-Ew water model: Vapor pressure and boiling point
,”
J. Chem. Phys.
123
(
19
),
194504
(
2005
).
15.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
, “
The missing term in effective pair potentials
,”
J. Phys. Chem.
91
(
24
),
6269
6271
(
1987
).
16.
M. W.
Mahoney
and
W. L.
Jorgensen
, “
A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions
,”
J. Chem. Phys.
112
(
20
),
8910
8922
(
2000
).
17.
V.
Molinero
and
E. B.
Moore
, “
Water modeled as an intermediate element between carbon and silicon
,”
J. Phys. Chem. B
113
(
13
),
4008
4016
(
2009
).
18.
T.
Hasegawa
and
Y.
Tanimura
, “
A polarizable water model for intramolecular and intermolecular vibrational spectroscopies
,”
J. Phys. Chem. B
115
(
18
),
5545
5553
(
2011
).
19.
P.
Ahlström
,
A.
Wallqvist
,
S.
Engström
, and
B.
Jönsson
, “
A molecular dynamics study of polarizable water
,”
Mol. Phys.
68
(
3
),
563
581
(
1989
).
20.
L.-P.
Wang
,
T.
Head-Gordon
,
J. W.
Ponder
,
P.
Ren
,
J. D.
Chodera
,
P. K.
Eastman
,
T. J.
Martinez
, and
V. S.
Pande
, “
Systematic improvement of a classical molecular model of water
,”
J. Phys. Chem. B
117
(
34
),
9956
9972
(
2013
).
21.
V.
Babin
,
C.
Leforestier
, and
F.
Paesani
, “
Development of a ‘first principles’ water potential with flexible monomers: Dimer potential energy surface, vrt spectrum, and second virial coefficient
,”
J. Chem. Theory Comput.
9
(
12
),
5395
5403
(
2013
).
22.
B.
Santra
,
A.
Michaelides
, and
M.
Scheffler
, “
On the accuracy of density-functional theory exchange-correlation functionals for H bonds in small water clusters: Benchmarks approaching the complete basis set limit
,”
J. Chem. Phys.
127
(
18
),
184104
(
2007
).
23.
E.
Lambros
and
F.
Paesani
, “
How good are polarizable and flexible models for water: Insights from a many-body perspective
,”
J. Chem. Phys.
153
(
6
),
060901
(
2020
).
24.
H. W.
Horn
,
W. C.
Swope
,
J. W.
Pitera
,
J. D.
Madura
,
T. J.
Dick
,
G. L.
Hura
, and
T.
Head-Gordon
, “
Development of an improved four-site water model for biomolecular simulations: TIP4P-Ew
,”
J. Chem. Phys.
120
(
20
),
9665
9678
(
2004
).
25.
A. L.
Benavides
,
J. L.
Aragones
, and
C.
Vega
, “
Consensus on the solubility of nacl in water from computer simulations using the chemical potential route
,”
J. Chem. Phys.
144
(
12
),
124504
(
2016
).
26.
J. R.
Espinosa
,
J. M.
Young
,
H.
Jiang
,
D.
Gupta
,
C.
Vega
,
E.
Sanz
,
P. G.
Debenedetti
, and
A. Z.
Panagiotopoulos
, “
On the calculation of solubilities via direct coexistence simulations: Investigation of nacl aqueous solutions and Lennard-Jones binary mixtures
,”
J. Chem. Phys.
145
(
15
),
154111
(
2016
).
27.
H.
Jiang
,
A.
Haji-Akbari
,
P. G.
Debenedetti
, and
A. Z.
Panagiotopoulos
, “
Forward flux sampling calculation of homogeneous nucleation rates from aqueous nacl solutions
,”
J. Chem. Phys.
148
(
4
),
044505
(
2018
).
28.
I.
Sanchez-Burgos
and
J. R.
Espinosa
, “
Direct calculation of the planar nacl-aqueous solution interfacial free energy at the solubility limit
,” arXiv:2208.08322 (
2022
).
29.
C. P.
Lamas
,
C.
Vega
, and
E. G.
Noya
, “
Freezing point depression of salt aqueous solutions using the Madrid-2019 model
,”
J. Chem. Phys.
156
(
13
),
134503
(
2022
).
30.
S.
Blazquez
,
M. M.
Conde
,
J. L. F.
Abascal
, and
C.
Vega
, “
The Madrid-2019 force field for electrolytes in water using TIP4P/2005 and scaled charges: Extension to the ions F, Br, I, Rb+, and Cs+
,”
J. Chem. Phys.
156
(
4
),
044505
(
2022
).
31.
C.
Bergonzo
and
T. E.
Cheatham
III
, “
Improved force field parameters lead to a better description of RNA structure
,”
J. Chem. Theory Comput.
11
(
9
),
3969
3972
(
2015
).
32.
M. D.
Smith
,
J. S.
Rao
,
E.
Segelken
, and
L.
Cruz
, “
Force-field induced bias in the structure of Aβ21–30: A comparison of OPLS, AMBER, CHARMM, and GROMOS force fields
,”
J. Chem. Inf. Model.
55
(
12
),
2587
2595
(
2015
).
33.
P. M.
Piaggi
and
R.
Car
, “
Phase equilibrium of liquid water and hexagonal ice from enhanced sampling molecular dynamics simulations
,”
J. Chem. Phys.
152
(
20
),
204116
(
2020
).
34.
E.
Sanz
,
C.
Vega
,
J. R.
Espinosa
,
R.
Caballero-Bernal
,
J. L. F.
Abascal
, and
C.
Valeriani
, “
Homogeneous ice nucleation at moderate supercooling from molecular simulation
,”
J. Am. Chem. Soc.
135
(
40
),
15008
15017
(
2013
).
35.
J. R.
Espinosa
,
C.
Vega
,
C.
Valeriani
, and
E.
Sanz
, “
Seeding approach to crystal nucleation
,”
J. Chem. Phys.
144
(
3
),
034501
(
2016
).
36.
H.
Niu
,
Y. I.
Yang
, and
M.
Parrinello
, “
Temperature dependence of homogeneous nucleation in ice
,”
Phys. Rev. Lett.
122
(
24
),
245501
(
2019
).
37.
R. A.
Friesner
, “
Ab initio quantum chemistry: Methodology and applications
,”
Proc. Natl. Acad. Sci. U. S. A.
102
(
19
),
6648
6653
(
2005
).
38.
V.
Mlýnský
,
P.
Banáš
,
J.
Šponer
,
M. W.
van der Kamp
,
A. J.
Mulholland
, and
M.
Otyepka
, “
Comparison of ab initio, DFT, and semiempirical QM/MM approaches for description of catalytic mechanism of hairpin ribozyme
,”
J. Chem. Theory Comput.
10
(
4
),
1608
1622
(
2014
).
39.
R. B.
Gerber
,
D.
Shemesh
,
M. E.
Varner
,
J.
Kalinowski
, and
B.
Hirshberg
, “
Ab initio and semi-empirical molecular dynamics simulations of chemical reactions in isolated molecules and in clusters
,”
Phys. Chem. Chem. Phys.
16
(
21
),
9760
9775
(
2014
).
40.
L.
Zhang
,
J.
Han
,
H.
Wang
,
R.
Car
, and
W.
E
, “
Deep potential molecular dynamics: A scalable model with the accuracy of quantum mechanics
,”
Phys. Rev. Lett.
120
(
14
),
143001
(
2018
).
41.
Y.
Zhang
,
H.
Wang
,
W.
Chen
,
J.
Zeng
,
L.
Zhang
,
H.
Wang
, and
W.
E
, “
DP-GEN: A concurrent learning platform for the generation of reliable deep learning based potential energy models
,”
Comput. Phys. Commun.
253
,
107206
(
2020
).
42.
L.
Zhang
,
H.
Wang
,
R.
Car
, and
W.
E
, “
Phase diagram of a deep potential water model
,”
Phys. Rev. Lett.
126
(
23
),
236001
(
2021
).
43.
J.
Sun
,
A.
Ruzsinszky
, and
J. P.
Perdew
, “
Strongly constrained and appropriately normed semilocal density functional
,”
Phys. Rev. Lett.
115
(
3
),
036402
(
2015
).
44.
J.
Sun
,
R. C.
Remsing
,
Y.
Zhang
,
Z.
Sun
,
A.
Ruzsinszky
,
H.
Peng
,
Z.
Yang
,
A.
Paul
,
U.
Waghmare
,
X.
Wu
et al, “
Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional
,”
Nat. Chem.
8
(
9
),
831
836
(
2016
).
45.
P. M.
Piaggi
,
J.
Weis
,
A. Z.
Panagiotopoulos
,
P. G.
Debenedetti
, and
R.
Car
, “
Homogeneous ice nucleation in an ab initio machine-learning model of water
,”
Proc. Natl. Acad. Sci. U. S. A.
119
(
33
),
e2207294119
(
2022
).
46.
J.
Han
,
L.
Zhang
,
R.
Car
et al, “
Deep potential: A general representation of a many-body potential energy surface
,” arXiv:1707.01478 (
2017
).
47.
P.
Linstrom
and
W.
Mallard
, NIST Chemistry WebBook, NIST Standard Reference Database No. 69, Gaithersburg MD, National Institute of Standards and Technology.
48.
C.
Vega
,
J. L. F.
Abascal
, and
I.
Nezbeda
, “
Vapor-liquid equilibria from the triple point up to the critical point for the new generation of TIP4P-like models: TIP4P/Ew, Tip4p/2005, and TIP4P/Ice
,”
J. Chem. Phys.
125
(
3
),
034503
(
2006
).
49.
C.
Vega
and
E.
de Miguel
, “
Surface tension of the most popular models of water by using the test-area simulation method
,”
J. Chem. Phys.
126
(
15
),
154707
(
2007
).
50.
R. D.
Mountain
, “
An internally consistent method for the molecular dynamics simulation of the surface tension: Application to some TIP4P-type models of water
,”
J. Phys. Chem. B
113
(
2
),
482
486
(
2009
).
51.
J.
Alejandre
and
G. A.
Chapela
, “
The surface tension of TIP4P/2005 water model using the Ewald sums for the dispersion interactions
,”
J. Chem. Phys.
132
(
1
),
014701
(
2010
).
52.
H. L.
Pi
,
J. L.
Aragones
,
C.
Vega
,
E. G.
Noya
,
J. L.
Abascal
,
M. A.
Gonzalez
, and
C.
McBride
, “
Anomalies in water as obtained from computer simulations of the TIP4P/2005 model: Density maxima, and density, isothermal compressibility and heat capacity minima
,”
Mol. Phys.
107
(
4–6
),
365
374
(
2009
).
53.
G.
Menzl
,
M. A.
Gonzalez
,
P.
Geiger
,
F.
Caupin
,
J. L. F.
Abascal
,
C.
Valeriani
, and
C.
Dellago
, “
Molecular mechanism for cavitation in water under tension
,”
Proc. Natl. Acad. Sci. U. S. A.
113
(
48
),
13582
13587
(
2016
).
54.
S. H.
Min
and
M. L.
Berkowitz
, “
Bubbles in water under stretch-induced cavitation
,”
J. Chem. Phys.
150
(
5
),
054501
(
2019
).
55.
J. L. F.
Abascal
,
M. A.
Gonzalez
,
J. L.
Aragones
, and
C.
Valeriani
, “
Homogeneous bubble nucleation in water at negative pressure: A Voronoi polyhedra analysis
,”
J. Chem. Phys.
138
(
8
),
084508
(
2013
).
56.
M. N.
Joswiak
,
N.
Duff
,
M. F.
Doherty
, and
B.
Peters
, “
Size-dependent surface free energy and Tolman-corrected droplet nucleation of TIP4P/2005 water
,”
J. Phys. Chem. Lett.
4
(
24
),
4267
4272
(
2013
).
57.
M. N.
Joswiak
,
R.
Do
,
M. F.
Doherty
, and
B.
Peters
, “
Energetic and entropic components of the Tolman length for mw and TIP4P/2005 water nanodroplets
,”
J. Chem. Phys.
145
(
20
),
204703
(
2016
).
58.
M. A.
Gonzalez
,
J. L. F.
Abascal
,
C.
Valeriani
, and
F.
Bresme
, “
Bubble nucleation in simple and molecular liquids via the largest spherical cavity method
,”
J. Chem. Phys.
142
(
15
),
154903
(
2015
).
59.
M. A.
González
,
G.
Menzl
,
J. L.
Aragones
,
P.
Geiger
,
F.
Caupin
,
J. L. F.
Abascal
,
C.
Dellago
, and
C.
Valeriani
, “
Detecting vapour bubbles in simulations of metastable water
,”
J. Chem. Phys.
141
(
18
),
18C511
(
2014
).
60.
C. P.
Lamas
,
C.
Vega
,
E. G.
Noya
, and
E.
Sanz
, “
The water cavitation line as predicted by the TIP4P/2005 model
,”
J. Chem. Phys.
158
,
124504
(
2023
).
61.
P. G.
Debenedetti
,
Metastable Liquids
(
Princeton university press
,
2021
).
62.
D.
Kashchiev
,
Nucleation
(
Elsevier
,
2000
).
63.
J. L.
Green
,
D. J.
Durben
,
G. H.
Wolf
, and
C. A.
Angell
, “
Water and solutions at negative pressure: Raman spectroscopic study to -80 megapascals
,”
Science
249
(
4969
),
649
652
(
1990
).
64.
Q.
Zheng
,
D. J.
Durben
,
G. H.
Wolf
, and
C. A.
Angell
, “
Liquids at large negative pressures: Water at the homogeneous nucleation limit
,”
Science
254
(
5033
),
829
832
(
1991
).
65.
A. D.
Alvarenga
,
M.
Grimsditch
, and
R. J.
Bodnar
, “
Elastic properties of water under negative pressures
,”
J. Chem. Phys.
98
(
11
),
8392
8396
(
1993
).
66.
M. E. M.
Azouzi
,
C.
Ramboz
,
J.-F.
Lenain
, and
F.
Caupin
, “
A coherent picture of water at extreme negative pressure
,”
Nat. Phys.
9
(
1
),
38
41
(
2013
).
67.
G.
Pallares
,
M.
El Mekki Azouzi
,
M. A.
González
,
J. L.
Aragones
,
J. L. F.
Abascal
,
C.
Valeriani
, and
F.
Caupin
, “
Anomalies in bulk supercooled water at negative pressure
,”
Proc. Natl. Acad. Sci. U. S. A.
111
(
22
),
7936
7941
(
2014
).
68.
R. C.
Tolman
, “
The effect of droplet size on surface tension
,”
J. Chem. Phys.
17
(
3
),
333
337
(
1949
).
69.
C. G.
Salzmann
,
P. G.
Radaelli
,
B.
Slater
, and
J. L.
Finney
, “
The polymorphism of ice: Five unresolved questions
,”
Phys. Chem. Chem. Phys.
13
(
41
),
18468
18480
(
2011
).
70.
W.
Wagner
,
T.
Riethmann
,
R.
Feistel
, and
A. H.
Harvey
, “
New equations for the sublimation pressure and melting pressure of H2O ice Ih
,”
J. Phys. Chem. Ref. Data
40
(
4
),
043103
(
2011
).
71.
A. J.
Brown
and
E.
Whalley
, “
Preliminary investigation of the phase boundaries between ice VI and VII and ice VI and VIII
,”
J. Chem. Phys.
45
(
11
),
4360
4361
(
1966
).
72.
D.
Lu
,
W.
Jiang
,
Y.
Chen
,
L.
Zhang
,
W.
Jia
,
H.
Wang
, and
M.
Chen
, “
DP compress: A model compression scheme for generating efficient deep potential models
,”
J. Chem. Theory Comput.
18
(
9
),
5559
5567
(
2022
).
73.
S.
Plimpton
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
74.
H.
Wang
,
L.
Zhang
,
J.
Han
, and
W.
E
, “
DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics
,”
Comput. Phys. Commun.
228
,
178
184
(
2018
).
75.
S.
Nosé
, “
A unified formulation of the constant temperature molecular dynamics methods
,”
J. Chem. Phys.
81
(
1
),
511
519
(
1984
).
76.
W. G.
Hoover
, “
Canonical dynamics: Equilibrium phase-space distributions
,”
Phys. Rev. A
31
,
1695
1697
(
1985
).
77.
W. G.
Hoover
, “
Constant-pressure equations of motion
,”
Phys. Rev. A
34
(
3
),
2499
(
1986
).
78.
J. G.
Kirkwood
and
F. P.
Buff
, “
The statistical mechanical theory of surface tension
,”
J. Chem. Phys.
17
(
3
),
338
343
(
1949
).
79.
H.
Bekker
,
H.
Berendsen
,
E.
Dijkstra
,
S.
Achterop
,
R.
Vondrumen
,
D.
Vanderspoel
,
A.
Sijbers
,
H.
Keegstra
, and
M.
Renardus
, “
GROMACS—A parallel computer for molecular-dynamics simulations
,” in 4th International Conference on Computational Physics (PC 92) (
World Scientific Publishing
,
1993
), pp.
252
256
.
80.
G.
Bussi
,
D.
Donadio
, and
M.
Parrinello
, “
Canonical sampling through velocity rescaling
,”
J. Chem. Phys.
126
(
1
),
014101
(
2007
).
81.
M.
Parrinello
and
A.
Rahman
, “
Polymorphic transitions in single crystals: A new molecular dynamics method
,”
J. Appl. Phys.
52
(
12
),
7182
7190
(
1981
).
82.
R. W.
Hockney
,
S. P.
Goel
, and
J. W.
Eastwood
, “
Quiet high-resolution computer models of a plasma
,”
J. Comput. Phys.
14
(
2
),
148
158
(
1974
).
83.
T.
Darden
,
D.
York
, and
L.
Pedersen
, “
Particle mesh Ewald: An Nlog(N) method for Ewald sums in large systems
,”
J. Chem. Phys.
98
(
12
),
10089
10092
(
1993
).
84.
U.
Essmann
,
L.
Perera
,
M. L.
Berkowitz
,
T.
Darden
,
H.
Lee
, and
L. G.
Pedersen
, “
A smooth particle mesh Ewald method
,”
J. Chem. Phys.
103
(
19
),
8577
8593
(
1995
).
85.
B.
Hess
,
H.
Bekker
,
H. J. C.
Berendsen
, and
J. G. E. M.
Fraaije
, “
Lincs: A linear constraint solver for molecular simulations
,”
J. Comput. Chem.
18
(
12
),
1463
1472
(
1997
).
86.
P.
Rosales-Pelaez
,
I.
Sanchez-Burgos
,
C.
Valeriani
,
C.
Vega
, and
E.
Sanz
, “
Seeding approach to nucleation in the n v t ensemble: The case of bubble cavitation in overstretched Lennard Jones fluids
,”
Phys. Rev. E
101
(
2
),
022611
(
2020
).
87.
P.
Montero de Hijes
,
J. R.
Espinosa
,
V.
Bianco
,
E.
Sanz
, and
C.
Vega
, “
Interfacial free energy and Tolman length of curved liquid–solid interfaces from equilibrium studies
,”
J. Phys. Chem. C
124
(
16
),
8795
8805
(
2020
).
88.
M.
Volmer
and
A.
Weber
, “
Keimbildung in übersättigten gebilden
,”
Z. Phys. Chem.
119U
(
1
),
277
301
(
1926
).
89.
R.
Becker
and
W.
Döring
, “
Kinetische behandlung der keimbildung in übersättigten dämpfen
,”
Ann. Phys.
416
(
8
),
719
752
(
1935
).
90.
R. P.
Sear
, “
The non-classical nucleation of crystals: Microscopic mechanisms and applications to molecular crystals, ice and calcium carbonate
,”
Int. Mater. Rev.
57
(
6
),
328
356
(
2012
).
91.
J.
Merikanto
,
E.
Zapadinsky
,
A.
Lauri
, and
H.
Vehkamäki
, “
Origin of the failure of classical nucleation theory: Incorrect description of the smallest clusters
,”
Phys. Rev. Lett.
98
(
14
),
145702
(
2007
).
92.
A.
Cacciuto
,
S.
Auer
, and
D.
Frenkel
, “
Breakdown of classical nucleation theory near isostructural phase transitions
,”
Phys. Rev. Lett.
93
(
16
),
166105
(
2004
).
93.
M.
Horsch
,
J.
Vrabec
, and
H.
Hasse
, “
Modification of the classical nucleation theory based on molecular simulation data for surface tension, critical nucleus size, and nucleation rate
,”
Phys. Rev. E
78
(
1
),
011603
(
2008
).
94.
D.
Moroni
,
P. R.
Ten Wolde
, and
P. G.
Bolhuis
, “
Interplay between structure and size in a critical crystal nucleus
,”
Phys. Rev. Lett.
94
(
23
),
235703
(
2005
).
95.
F.
Leoni
and
J.
Russo
, “
Nonclassical nucleation pathways in stacking-disordered crystals
,”
Phys. Rev. X
11
(
3
),
031006
(
2021
).
96.
J. R.
Espinosa
,
C.
Vega
, and
E.
Sanz
, “
Homogeneous ice nucleation rate in water droplets
,”
J. Phys. Chem. C
122
(
40
),
22892
22896
(
2018
).
97.
B. C.
Knott
,
V.
Molinero
,
M. F.
Doherty
, and
B.
Peters
, “
Homogeneous nucleation of methane hydrates: Unrealistic under realistic conditions
,”
J. Am. Chem. Soc.
134
(
48
),
19544
19547
(
2012
).
98.
C. P.
Lamas
,
J. R.
Espinosa
,
M. M.
Conde
,
J.
Ramírez
,
P.
Montero de Hijes
,
E. G.
Noya
,
C.
Vega
, and
E.
Sanz
, “
Homogeneous nucleation of nacl in supersaturated solutions
,”
Phys. Chem. Chem. Phys.
23
(
47
),
26843
26852
(
2021
).
99.
I.
Saika-Voivod
,
F.
Romano
, and
F.
Sciortino
, “
Nucleation barriers in tetrahedral liquids spanning glassy and crystallizing regimes
,”
J. Chem. Phys.
135
(
12
),
124506
(
2011
).
100.
D.
Richard
and
T.
Speck
, “
Crystallization of hard spheres revisited. ii. thermodynamic modeling, nucleation work, and the surface of tension
,”
J. Chem. Phys.
148
(
22
),
224102
(
2018
).
101.
L.
Separdar
,
J. P.
Rino
, and
E. D.
Zanotto
, “
Molecular dynamics simulations of spontaneous and seeded nucleation and theoretical calculations for zinc selenide
,”
Comput. Mater. Sci.
187
,
110124
(
2021
).
102.
I.
Sanchez-Burgos
,
A.
Garaizar
,
C.
Vega
,
E.
Sanz
, and
J. R.
Espinosa
, “
Parasitic crystallization of colloidal electrolytes: Growing a metastable crystal from the nucleus of a stable phase
,”
Soft Matter
17
(
3
),
489
505
(
2021
).
103.
I.
Sanchez-Burgos
,
E.
Sanz
,
C.
Vega
, and
J. R.
Espinosa
, “
Fcc vs. hcp competition in colloidal hard-sphere nucleation: On their relative stability, interfacial free energy and nucleation rate
,”
Phys. Chem. Chem. Phys.
23
(
35
),
19611
19626
(
2021
).
104.
L.
Filion
,
M.
Hermes
,
R.
Ni
, and
M.
Dijkstra
, “
Crystal nucleation of hard spheres using molecular dynamics, umbrella sampling, and forward flux sampling: A comparison of simulation techniques
,”
J. Chem. Phys.
133
(
24
),
244115
(
2010
).
105.
G.
Fiorucci
,
G. M.
Coli
,
J. T.
Padding
, and
M.
Dijkstra
, “
The effect of hydrodynamics on the crystal nucleation of nearly hard spheres
,”
J. Chem. Phys.
152
(
6
),
064903
(
2020
).
106.
I.
Saika-Voivod
,
P. H.
Poole
, and
R. K.
Bowles
, “
Test of classical nucleation theory on deeply supercooled high-pressure simulated silica
,”
J. Chem. Phys.
124
(
22
),
224709
(
2006
).
107.
V. G.
Baidakov
, “
Stability of metastable phases and kinetics of nucleation in a simple single-component system (molecular dynamics simulation) (a review)
,”
Russ. J. Gen. Chem.
92
(
4
),
611
628
(
2022
).
108.
I.
Sanchez-Burgos
,
P. M.
de Hijes
,
P.
Rosales-Pelaez
,
C.
Vega
, and
E.
Sanz
, “
Equivalence between condensation and boiling in a Lennard-Jones fluid
,”
Phys. Rev. E
102
(
6
),
062609
(
2020
).
109.
M.
Blander
and
J. L.
Katz
, “
Bubble nucleation in liquids
,”
AIChE J.
21
(
5
),
833
848
(
1975
).
110.
K. M.
Bal
and
E. C.
Neyts
, “
Extending and validating bubble nucleation rate predictions in a Lennard-Jones fluid with enhanced sampling methods and transition state theory
,”
J. Chem. Phys.
157
(
18
),
184113
(
2022
).
111.
J.
Diemand
,
R.
Angélil
,
K. K.
Tanaka
, and
H.
Tanaka
, “
Direct simulations of homogeneous bubble nucleation: Agreement with classical nucleation theory and no local hot spots
,”
Phys. Rev. E
90
(
5
),
052407
(
2014
).
112.
P.
Montero de Hijes
,
K.
Shi
,
E. G.
Noya
,
E. E.
Santiso
,
K. E.
Gubbins
,
E.
Sanz
, and
C.
Vega
, “
The Young–Laplace equation for a solid–liquid interface
,”
J. Chem. Phys.
153
(
19
),
191102
(
2020
).
113.
S. L.
Meadley
and
F. A.
Escobedo
, “
Thermodynamics and kinetics of bubble nucleation: Simulation methodology
,”
J. Chem. Phys.
137
(
7
),
074109
(
2012
).
114.
I.
Kusaka
and
D. W.
Oxtoby
, “
Identifying physical clusters in bubble nucleation
,”
J. Chem. Phys.
111
(
3
),
1104
1108
(
1999
).
115.
R.
Angélil
,
J.
Diemand
,
K. K.
Tanaka
, and
H.
Tanaka
, “
Bubble evolution and properties in homogeneous nucleation simulations
,”
Phys. Rev. E
90
(
6
),
063301
(
2014
).
116.
M.
Chen
,
H.-Y.
Ko
,
R. C.
Remsing
,
M. F.
Calegari Andrade
,
B.
Santra
,
Z.
Sun
,
A.
Selloni
,
R.
Car
,
M. L.
Klein
,
J. P.
Perdew
, and
X.
Wu
, “
Ab initio theory and modeling of water
,”
Proc. Natl. Acad. Sci. U. S. A.
114
(
41
),
10846
10851
(
2017
).
117.
Y.
Yao
and
Y.
Kanai
, “
Temperature dependence of nuclear quantum effects on liquid water via artificial neural network model based on scan meta-GGA functional
,”
J. Chem. Phys.
153
(
4
),
044114
(
2020
).
118.
T. E.
Gartner
,
L.
Zhang
,
P. M.
Piaggi
,
R.
Car
,
A. Z.
Panagiotopoulos
, and
P. G.
Debenedetti
, “
Signatures of a liquid–liquid transition in an ab initio deep neural network model for water
,”
Proc. Natl. Acad. Sci. U. S. A.
117
(
42
),
26040
26046
(
2020
).
119.
P. M.
Piaggi
,
A. Z.
Panagiotopoulos
,
P. G.
Debenedetti
, and
R.
Car
, “
Phase equilibrium of water with hexagonal and cubic ice using the scan functional
,”
J. Chem. Theory Comput.
17
(
5
),
3065
3077
(
2021
).
120.
A.
Stukowski
, “
Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool
,”
Modell. Simul. Mater. Sci. Eng.
18
(
1
),
015012
(
2009
).
121.
J. S.
Rowlinson
and
B.
Widom
,
Molecular Theory of Capillarity
(
Courier Corporation
,
2013
).
122.
J. A.
Zollweg
and
G. W.
Mulholland
, “
On the law of the rectilinear diameter
,”
J. Chem. Phys.
57
(
3
),
1021
1025
(
1972
).
123.
Y.
Zhai
,
A.
Caruso
,
S. L.
Bore
,
Z.
Luo
, and
F.
Paesani
, “
A ‘short blanket’ dilemma for a state-of-the-art neural network potential for water: Reproducing experimental properties or the physics of the underlying many-body interactions?
,”
J. Chem. Phys.
158
(
8
),
084111
(
2023
).
124.
G. M.
Torrie
and
J. P.
Valleau
, “
Monte Carlo free energy estimates using non-Boltzmann sampling: Application to the sub-critical Lennard-Jones fluid
,”
Chem. Phys. Lett.
28
(
4
),
578
581
(
1974
).
125.
A. P.
Thompson
,
S. J.
Plimpton
, and
W.
Mattson
, “
General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions
,”
J. Chem. Phys.
131
(
15
),
154107
(
2009
).
126.
A. E.
van Giessen
and
E. M.
Blokhuis
, “
Direct determination of the Tolman length from the bulk pressures of liquid drops via molecular dynamics simulations
,”
J. Chem. Phys.
131
(
16
),
164705
(
2009
).
127.
T. V.
Bykov
and
X. C.
Zeng
, “
A patching model for surface tension and the Tolman length
,”
J. Chem. Phys.
111
(
8
),
3705
3713
(
1999
).
128.
Y. A.
Lei
,
T.
Bykov
,
S.
Yoo
, and
X. C.
Zeng
, “
The Tolman length: Is it positive or negative?
,”
J. Am. Chem. Soc.
127
(
44
),
15346
15347
(
2005
).
129.
J.
Julin
,
I.
Napari
,
J.
Merikanto
, and
H.
Vehkamäki
, “
A thermodynamically consistent determination of surface tension of small Lennard-Jones clusters from simulation and theory
,”
J. Chem. Phys.
133
(
4
),
044704
(
2010
).
130.
D.
Kashchiev
, “
Nucleation work, surface tension, and Gibbs–Tolman length for nucleus of any size
,”
J. Chem. Phys.
153
(
12
),
124509
(
2020
).
131.
A.
Tröster
and
K.
Binder
, “
Positive Tolman length in a lattice gas with three-body interactions
,”
Phys. Rev. Lett.
107
(
26
),
265701
(
2011
).
132.
M.
Iwamatsu
and
K.
Horii
, “
Temperature dependence of liquid-vapor nucleation rate for the Yukawa model fluid
,”
Aerosol Sci. Technol.
27
(
5
),
563
574
(
1997
).
133.
F.
Magaletti
,
M.
Gallo
, and
C. M.
Casciola
, “
Water cavitation from ambient to high temperatures
,”
Sci. Rep.
11
(
1
),
20801
(
2021
).
134.
F.
Tang
,
T.
Ohto
,
S.
Sun
,
J. R.
Rouxel
,
S.
Imoto
,
E. H. G.
Backus
,
S.
Mukamel
,
M.
Bonn
, and
Y.
Nagata
, “
Molecular structure and modeling of water–air and ice–air interfaces monitored by sum-frequency generation
,”
Chem. Rev.
120
(
8
),
3633
3667
(
2020
).
135.
Y.
Nagata
,
T.
Hasegawa
,
E. H. G.
Backus
,
K.
Usui
,
S.
Yoshimune
,
T.
Ohto
, and
M.
Bonn
, “
The surface roughness, but not the water molecular orientation varies with temperature at the water–air interface
,”
Phys. Chem. Chem. Phys.
17
(
36
),
23559
23564
(
2015
).
136.
M.
Bonn
,
Y.
Nagata
, and
E. H. G.
Backus
, “
Molecular structure and dynamics of water at the water–air interface studied with surface-specific vibrational spectroscopy
,”
Angew. Chem., Int. Ed.
54
(
19
),
5560
5576
(
2015
).
137.
S.
Nihonyanagi
,
J. A.
Mondal
,
S.
Yamaguchi
, and
T.
Tahara
, “
Structure and dynamics of interfacial water studied by heterodyne-detected vibrational sum-frequency generation
,”
Annu. Rev. Phys. Chem.
64
,
579
603
(
2013
).
138.
S.
Nihonyanagi
,
T.
Ishiyama
,
T.-k.
Lee
,
S.
Yamaguchi
,
M.
Bonn
,
A.
Morita
, and
T.
Tahara
, “
Unified molecular view of the air/water interface based on experimental and theoretical χ(2) spectra of an isotopically diluted water surface
,”
J. Am. Chem. Soc.
133
(
42
),
16875
16880
(
2011
).
139.
C.
Liang
,
J.
Jeon
, and
M.
Cho
, “
Ab initio modeling of the vibrational sum-frequency generation spectrum of interfacial water
,”
J. Phys. Chem. Lett.
10
(
5
),
1153
1158
(
2019
).
140.
T.
Ohto
,
K.
Usui
,
T.
Hasegawa
,
M.
Bonn
, and
Y.
Nagata
, “
Toward ab initio molecular dynamics modeling for sum-frequency generation spectra; an efficient algorithm based on surface-specific velocity-velocity correlation function
,”
J. Chem. Phys.
143
(
12
),
124702
(
2015
).
141.
J.
Kessler
,
H.
Elgabarty
,
T.
Spura
,
K.
Karhan
,
P.
Partovi-Azar
,
A. A.
Hassanali
, and
T. D.
Kühne
, “
Structure and dynamics of the instantaneous water/vapor interface revisited by path-integral and ab initio molecular dynamics simulations
,”
J. Phys. Chem. B
119
(
31
),
10079
10086
(
2015
).
142.
P.
Vassilev
,
C.
Hartnig
,
M. T. M.
Koper
,
F.
Frechard
, and
R. A.
van Santen
, “
Ab initio molecular dynamics simulation of liquid water and water–vapor interface
,”
J. Chem. Phys.
115
(
21
),
9815
9820
(
2001
).
143.
R.
Khatib
,
T.
Hasegawa
,
M.
Sulpizi
,
E. H. G.
Backus
,
M.
Bonn
, and
Y.
Nagata
, “
Molecular dynamics simulations of SFG librational modes spectra of water at the water–air interface
,”
J. Phys. Chem. C
120
(
33
),
18665
18673
(
2016
).
144.
C. Y.
Lee
,
J. A.
McCammon
, and
P. J.
Rossky
, “
The structure of liquid water at an extended hydrophobic surface
,”
J. Chem. Phys.
80
(
9
),
4448
4455
(
1984
).
145.
Y.
Fan
,
X.
Chen
,
L.
Yang
,
P. S.
Cremer
, and
Y. Q.
Gao
, “
On the structure of water at the aqueous/air interface
,”
J. Phys. Chem. B
113
(
34
),
11672
11679
(
2009
).
146.
J. R.
Espinosa
,
C.
Vega
, and
E.
Sanz
, “
Ice–water interfacial free energy for the TIP4P, TIP4P/2005, TIP4P/Ice, and mw models as obtained from the mold integration technique
,”
J. Phys. Chem. C
120
(
15
),
8068
8075
(
2016
).
147.
S. P.
Niblett
,
M.
Galib
, and
D. T.
Limmer
, “
Learning intermolecular forces at liquid–vapor interfaces
,”
J. Chem. Phys.
155
(
16
),
164101
(
2021
).

Supplementary Material