Aided by a neural network representation of the density functional theory potential energy landscape of water in the Revised Perdew–Burke–Ernzerhof approximation corrected for dispersion, we calculate several structural and thermodynamic properties of its liquid/vapor interface. The neural network speed allows us to bridge the size and time scale gaps required to sample the properties of water along its liquid/vapor coexistence line with unprecedented precision.

The strong, directed network of hydrogen bonds1 that confers water its rich phase diagram, and its numerous anomalous properties,2 is responsible for the peculiar structure of its liquid/vapor interface.3 Thanks to its hydrogen bond network, water is one of the most cohesive liquids known with exceptionally high surface tension.

The anomalous surface tension is, by far, not the only surprising property of water interfaces. Unlike simple liquids, water undergoes substantial structural changes in the interface’s first molecular layer. In its preferential configuration, a surface molecule at the liquid/vapor interface reorients its dipole vector along the liquid surface, with one OH bond and one lone pair both directed toward the vapor phase.4–6 The rearrangement of water molecules at the interface with other phases, with the corresponding entropic change, lies at the very heart of phenomena such as the hydrophobic effect,7 one of the main actors of self-organization in living organisms, and of many physicochemical based applications, from detergents to novel materials based on microemulsions.

Understanding the relationship between the microscopic structure of the water/vapor interface and its thermodynamic properties is then key to obtaining better control over a vast array of processes. Due to its molecular-scale extension, interfacial water can be investigated experimentally by a small set of techniques, most importantly x-ray reflectivity8 and sum-frequency generation spectroscopy.9 In this sense, computer simulation techniques, providing direct access to the molecular configurations, are a unique asset. However, the accuracy of the calculations and their computational cost are two significant challenges for the simulation of water, in general, and its interfaces, in particular. Ideally, one would aim at an accurate parameter-free description. While much progress has been made to achieve better ab initio simulations of liquid water,10,11 the computational resources required to simulate liquid interfaces using state-of-the-art methods are still daunting. Artificial neural networks have proven to be a valuable tool12 to reproduce the features of the potential energy landscapes of water.13 Using an artificial neural network to reproduce the forces acting between nuclei as computed with an ab initio method of choice can decrease the computational cost of the problem by several orders of magnitude.14 However, the neural network-predicted forces associated with a specific configuration can be unreliable if the network has not encountered sufficiently similar configurations during its training phase.

Here, we extend previously parameterized neural networks13,15 by including explicitly interfacial configurations and use the optimized potentials to study the phase diagram of water along its liquid/vapor coexistence line. The need for relatively large systems is particularly pressing for interfacial systems as they suffer more than bulk systems from finite-size effects.16 The computational advantage warranted by the neural network allowed us to simulate for many nanoseconds simulation boxes containing 1024 water molecules based on the Revised Perdew–Burke–Ernzerhof (RPBE) generalized gradient approximation for the exchange–correlation functional17 supplemented by Grimme’s D3 dispersion corrections.18 

As a first step, we simulated a water/vapor interface using the neural network for the RPBE exchange–correlation functional with Grimme’s D3 corrections presented in Ref. 13. This network was trained using configurations taken from the liquid, solid, and liquid/solid coexistence phases.15 As expected, the network could not recognize a relatively high fraction of configurations at the liquid/vapor interface, resulting in unphysically high vapor densities and large interfacial widths. A preliminary test with the Flexible Simple Point-Charge (SPC/F)19 empirical model showed that using a training dataset of about 400 frames generated from equilibrium trajectories of a liquid/vapor interface (216 molecules) yields excellent convergence and physically sound results. The density in the liquid phase obtained by the neural network differed from explicit SPC/F simulation from a minimum of 0.6% at 300 K up to a maximum of 3% at 460 K, while the surface tension differences were always below 1 mN/m.

We used these 400 frames to augment the dataset presented in Ref. 13, taking care of perturbing the nuclear coordinates (by a uniform distribution up to ±0.003 nm) to enhance the sampling of the neighborhood of the free energy minima. The density functional theory (DFT) energies and forces at the RPBE level were calculated for the whole training dataset using Fritz Haber Institute ab initio molecular simulations package (FHI-aims),20 and D3 corrections were added. We then used this initial extension to train the neural network and generated a trajectory with a temperature ramp from 300 K to 600 K over 1 ns of integration. Next, we took about 500 (unperturbed) frames from this trajectory to extend the training set further, performing in this way a self-consistent refinement step. We report additional methodological details in Sec. IV.

Here, we would like to stress that without the extension of the training dataset, no meaningful simulation of the liquid/vapor equilibrium could be performed. For example, the interface would become unstable already at about 550 K (to be compared with a critical point of 632 K), and unphysical configurations in the first liquid layer, with the majority of molecular dipoles pointing toward the gas phase, instead of being parallel to the interface. Note that the new training set with interfacial configurations also includes the bulk configurations used in Ref. 13, and the thermodynamic properties of the liquid state are compatible with those reported there, at comparable thermodynamic points.

Simulating water using the SPC/F empirical potential also allowed us to test the effect of using a cutoff in the neural network, as opposed to the inclusion of long range forces with mesh Ewald methods. We performed additional simulations with the SPC/F model at 300 K, employing the same cutoff used for the neural network (0.635 nm), which yielded a similar (1% difference) density of the liquid phase, but a remarkably different surface tension, with a difference of about 9 mN/m. This confirms, albeit indirectly that the neural network encodes the effects of long-range forces through the local arrangement of atoms and their forces, even if a simple cutoff is used. With the present simulations, it is not possible to tell, however, to which extent the inaccuracy of the results obtained with the neural network can be ascribed to the lack of explicit information about further neighbors, by the choice of symmetry functions, or by intrinsic limitations of the neural network itself.

We ran molecular dynamics simulations in the canonical ensemble using the neural network potential for systems of 1024 water molecules in slab configurations within simulation boxes of size 3 × 3 × 10 nm3 under periodic boundary conditions at 15 different temperature values ranging from 300 K to 620 K (Fig. 1). Each simulation started from a pre-equilibrated configuration of the SPC/F model at the corresponding temperature using a time step of 0.5 fs. After 100 ps, we observed no drift in the potential energy and in the density profile of each trajectory, and we deemed to have reached equilibrium. We saved configurations to disk during the subsequent 1 ns of trajectory at an interval of 1 ps for further analysis. The values of energy and pressure were dumped every 10 fs.

FIG. 1.

Simulation snapshot of the simulation box of 1024 water molecules at 550 K. The oxygen atoms are colored in red (liquid phase), blue (interfacial layer), and orange (vapor phase).

FIG. 1.

Simulation snapshot of the simulation box of 1024 water molecules at 550 K. The oxygen atoms are colored in red (liquid phase), blue (interfacial layer), and orange (vapor phase).

Close modal

From the mass density profiles reported in Fig. 2, one can appreciate the progressive broadening of the slab width, accompanied by the density increase in the vapor phase and corresponding decrease in the liquid one. In Fig. 3, we report the density values of the 1 nm-wide regions located in the middle of the liquid and vapor slab. To extrapolate the critical point location, we performed the best fit of the simulation results to the expression proposed by Wilding,23 

ρ(T)=ρc+aTTc±b(TTc)β,
(1)

where ρc and Tc are the critical density and temperature, and the sign plus and minus applies to the liquid and vapor branch, respectively. We used the scaling exponent β as a fitting parameter, in addition to a and b, obtaining as best fit estimates ρc = 310 ± 3 kg/m3, Tc = 632 ± 2 K, β = 0.27 ± 0.01, a = 0.45 ± 0.01, and b = 94 ± 5.

FIG. 2.

Mass density profiles across the whole simulation box for the whole range of temperature investigated (see Table I) from 300 K (highest density in the liquid phase) to 620 K (lowest density in the liquid phase).

FIG. 2.

Mass density profiles across the whole simulation box for the whole range of temperature investigated (see Table I) from 300 K (highest density in the liquid phase) to 620 K (lowest density in the liquid phase).

Close modal
FIG. 3.

Water liquid/vapor coexistence from simulations using the RPBE-D3 approximation (circles, this work and crosses, Ref. 21) and using the TIP4-2005 empirical potential (triangles, from Ref. 22). Best fit of the data from the present work to Eq. (1) (dotted-dashed line), the corresponding estimate of the critical point (star), and the IAPWS95 equation of state (dashed line).

FIG. 3.

Water liquid/vapor coexistence from simulations using the RPBE-D3 approximation (circles, this work and crosses, Ref. 21) and using the TIP4-2005 empirical potential (triangles, from Ref. 22). Best fit of the data from the present work to Eq. (1) (dotted-dashed line), the corresponding estimate of the critical point (star), and the IAPWS95 equation of state (dashed line).

Close modal
TABLE I.

Liquid (ρl) and vapor (ρg) densities and surface tension (γ) as a function of the temperature. The estimated critical point is reported in the last line. Values in parentheses represent one standard deviation in the least significant digits.

T/Kρl (kg/m3)ρg (kg/m3)γ (mN/m)
300 899(1) 0.03(1) 68(2) 
320 889(1) 0.02(1) 64(2) 
350 876(1) 0.20(2) 57(2) 
370 859(1) 0.37(5) 51(2) 
400 836(1) 1.41(5) 46(2) 
420 817(1) 2.57(8) 42(2) 
440 797(1) 3.9(2) 40(2) 
450 788(1) 4.7(1) 36(2) 
480 752(1) 9.4(3) 33(2) 
500 728(1) 12.7(2) 28(2) 
520 702(1) 15.8(2) 21(2) 
550 657(1) 31.2(3) 15(2) 
580 597(1) 53.3(6) 10(2) 
600 550(1) 86(1) 8(2) 
620 509(1) 132(2) 3(2) 
632(2) 310(2) 
T/Kρl (kg/m3)ρg (kg/m3)γ (mN/m)
300 899(1) 0.03(1) 68(2) 
320 889(1) 0.02(1) 64(2) 
350 876(1) 0.20(2) 57(2) 
370 859(1) 0.37(5) 51(2) 
400 836(1) 1.41(5) 46(2) 
420 817(1) 2.57(8) 42(2) 
440 797(1) 3.9(2) 40(2) 
450 788(1) 4.7(1) 36(2) 
480 752(1) 9.4(3) 33(2) 
500 728(1) 12.7(2) 28(2) 
520 702(1) 15.8(2) 21(2) 
550 657(1) 31.2(3) 15(2) 
580 597(1) 53.3(6) 10(2) 
600 550(1) 86(1) 8(2) 
620 509(1) 132(2) 3(2) 
632(2) 310(2) 

In Fig. 3, we also report the IAPWS95 curve,24 which matches the experimental data, and the points from the only other estimate of the coexistence line of RPBE-D3 water we are aware of, taken from the work of Schienbein and Marx.21 Notice the steeper trend of Schienbein and Marx’s data, obtained via ab initio Gibbs ensemble Monte Carlo25,26 of 128 water molecules, which is arguably a finite-size effect yielding an effective critical temperature Tc(L) that shifts toward higher values with decreasing system (linear) size L as27 

Tc(L)Tc1L1/ν,
(2)

where ν is the critical exponent governing the scaling of the correlation length.

From the average values of the pressure tensor elements, pij, it is straightforward to compute the surface tension γ using the mechanical route as

γ=Lz2pzzpxx+pyy/2,
(3)

where z identifies the direction normal to the macroscopic interface. In Fig. 4, we report the surface tension as a function of temperature, the result of the best fit to Eq. (1) from Ref. 28, and we compare them with the interpolated experimental values. Even though the present results do not match the experimental curve and the best empirical potential models such as TIP4P-2005,29 the agreement is still very good and superior to several other mainstream empirical models.22 

FIG. 4.

Surface tension of RPBE-D3 as computed in this work (circles) and of TIP4P-2005 (triangles, taken from Ref. 22). The dotted-dashed line is the result of a best fit to Eq. (1) from Ref. 28. The dashed line is the best-fit to experimental data reported in Ref. 28.

FIG. 4.

Surface tension of RPBE-D3 as computed in this work (circles) and of TIP4P-2005 (triangles, taken from Ref. 22). The dotted-dashed line is the result of a best fit to Eq. (1) from Ref. 28. The dashed line is the best-fit to experimental data reported in Ref. 28.

Close modal

Density and surface tension values along the coexistence line are all reported in Table I, together with the estimated critical point.

With plenty of configurations at hand, it is now possible to investigate, in a statistically meaningful way, how water’s structural features are dependent on whether the molecules are located right at the interface or below it. Molecular dynamics simulations with empirical potentials show that at the liquid/vapor interface, water exhibits a much weaker order than van der Waals liquids, and many structural and dynamical properties reach their bulk value roughly after the second interfacial layer. In our investigation, we witnessed the same behavior, and here, we only concentrate on the properties of the first molecular layer, detected on a per-frame basis using the Pytim analysis package,30 as described in Sec. IV. For brevity, we will refer to “bulk properties” extracted from the trajectories presented in this work, as those properties computed from the molecules beyond the third surface layer, where the observables already converged to position-independent values.

We also calculated the bond length and angle distributions in the first and subsequent layers and found that the molecules in the first layer are less stretched, with a difference, at 300 K, of 0.3 pm and 0.3° for the bond length and angle, respectively. These differences become less pronounced when the temperature is increased.

Next, we characterize the orientation of the water molecules at the surface. To this purpose, we employ two angles. The first one, θ, is the angle that the molecular axis vector (pointing from the oxygen atom to the midpoint between the hydrogen ones) is making with the macroscopic surface normal (pointing from the liquid to the vapor phase). The second one, ϕ, is the angle that the molecule would need to rotate along its molecular axis to have the H–H vector aligned parallel to the surface plane.4,31 In Fig. 5, we report the free energy of molecules in the first layer as functions of cos θ and ϕ; for a randomly distributed molecular orientation, the histograms would be homogeneous.

FIG. 5.

Free energy map at 300 K of the orientation of water molecules in the first molecular layer (interpolated using cubic splines). The overlaid water molecules, represented using the position of the nuclei and their electronic charge density (85% isodensity surface), show (approximately) their orientation in the physical space with respect to the macroscopic surface plane for some selected pairs (ϕ, cos θ).

FIG. 5.

Free energy map at 300 K of the orientation of water molecules in the first molecular layer (interpolated using cubic splines). The overlaid water molecules, represented using the position of the nuclei and their electronic charge density (85% isodensity surface), show (approximately) their orientation in the physical space with respect to the macroscopic surface plane for some selected pairs (ϕ, cos θ).

Close modal

The free energy plot shows a prevalence of molecular orientations (the minima) when the dipole moment is aligned parallel to the surface plane or pointing slightly below it (cos θ ≃ −0.25) and when the OH bonds either lying in the surface plane (ϕ ≃ 0) or pointing out of it (ϕπ/2). This result agrees qualitatively with the previous experimental findings and simulations with empirical potentials. The orientations with the dipole moment pointing toward the vapor (cos θ = 1) or the liquid (cos θ = −1) are less likely to be found than the parallel orientation, although all orientations are accessible within an energy of kBT, where kB is the Boltzmann constant. When the temperature increases, the free energy map becomes gradually flatter, and although the preferential orientation of the dipole vector is always parallel to the surface, above 550 K, all orientations are accessible within an energy of 0.2 kBT (see the supplementary material).

Having access to the set of surface molecules, it is possible to build a local reference frame (x, y, ζ(x, y)) on the corrugated interface and use it to compute the intrinsic density profile,32 as reported in Fig. 6. The usual density profile expresses the correlation between molecular positions in the system and the location of its center of mass. In contrast, the intrinsic density profile, ρI(z′), represents the correlation between molecular positions and the local position of the interface,

ρI(z)=mAiNδ(zzi+ζ(xi,yi)),
(4)

where A is the simulation box cross-sectional area, N is the number of molecules, m is their mass, and angular brackets stand for the canonical average. The molecules in the surface layer are located at z′ = 0, giving rise to a delta peak, which is not shown in Fig. 6. By convention, the positive values of z′ are located in the vapor phase. At a relatively low temperature, where the vapor density is small and the liquid’s cohesive strength is the highest, the structure of the local packing emerges with a clear peak and a small shoulder, similar to the results from simulations with empirical potentials. With increasing temperature, the local structure at the interface becomes less pronounced and almost disappears at 620 K. On the vapor side, one observes a similar behavior, with a peak of the vapor density next to the interface that vanishes at high temperature. One can expect condensing vapor at the interface because the molecules in the vapor phase feel both the attractive dispersion forces and residual dipolar interactions. Upon a closer look at the vapor phase region in a semi-logarithmic scale [Fig. 6(b)], we note that the density, beyond the local maximum, decays exponentially toward the bulk vapor density like exp(−z/ξ), as it is expected for a dilute vapor next to a weakly attractive surface (the liquid phase) in the mean field approximation.33,34 The closer to the critical temperature, the more the profile tends toward homogeneity (with vapor density approaching the liquid one), except for an excluded volume region right at the interface that resembles the square-well one would expect from a hard-sphere.

FIG. 6.

(a) Intrinsic mass density profiles of water molecules next to the interface for the whole range of temperature investigated (see Table I) from 300 K (highest density in the liquid phase) to 620 K (lowest density in the liquid phase). The vapor phase is located by convention at positions on the right (positive), relative to the surface. The delta-like contribution at the origin is not shown.(b) Details of the vapor side close to the interface, in the semilogarithmic scale. Thin dashed lines are the result of best fits to exponential decays to the bulk vapor density.

FIG. 6.

(a) Intrinsic mass density profiles of water molecules next to the interface for the whole range of temperature investigated (see Table I) from 300 K (highest density in the liquid phase) to 620 K (lowest density in the liquid phase). The vapor phase is located by convention at positions on the right (positive), relative to the surface. The delta-like contribution at the origin is not shown.(b) Details of the vapor side close to the interface, in the semilogarithmic scale. Thin dashed lines are the result of best fits to exponential decays to the bulk vapor density.

Close modal

The neural network potential has been set up using the approach by Behler and Parrinello,14 using the implementation n2p2 of Singraber and co-workers.35 We used the same selection of symmetry functions and cutoff values reported in Ref. 13. For the sake of consistency and simplicity, this network is not trained to reproduce charge distributions (although it is, in principle, possible). Hence, no electrostatic quantities can be directly computed from the present simulation results. The network was able to reproduce the forces acting on atoms within 1.5 meV/Å (1 standard deviation, ≃68%), 2.2 meV/Å (2 standard deviations, ≃95%), and 4 meV/Å (3 standard deviations, ≃99.7%), with a typical rms error of 3% (see the supplementary material). DFT calculations were performed using the FHI-AIMS package,20 with the second tier level of basis functions for oxygen and hydrogen atoms, reproducing without any significant difference the forces and energies of the old dataset of Ref. 13. Neural network molecular dynamics simulations were performed using Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) (http://lammps.sandia.gov),36 linked to the neural network library of Singraber and co-workers (https://github.com/CompPhysVienna/n2p2).35 Each system consisted of 1024 water molecules simulated in the canonical ensemble using the Nosé–Hoover thermostat (damping constant 0.5 ps) and an integration time step of 0.5 fs. The center of mass was prevented from drifting by subtracting the center of mass velocity at every step, and the liquid slab was kept in this way in the middle of the simulation box along the surface normal direction. The surface layer was determined using the MDAnalysis37-based pytim package (https://github.com/Marcello-Sega/pytim)30 via a combination of the Identification of Truly Interfacial Molecules (ITIM)38 method (probe sphere radius 2 Å) and DBSCAN39-based prefiltering of the vapor phase,40 using the automatic determination of the density threshold for the clustering procedure.

We have performed what is arguably the most accurate calculation, to date, of the liquid/vapor coexistence of water described by the RPBE exchange–correlation functional, supplemented by dispersion corrections. This result was possible thanks to the use of a neural network-based fit of the DFT potential energy surface. We reported the coexistence curve of the system and estimated the critical temperature of the model (Tc = 632 ± 2 K), the surface tension curve as a function of temperature, and two order parameters, namely, the density profile and the orientation of water molecules in the surface layer. While the most refined empirical potential models are still superior in describing some aspects of the thermodynamics of water at interfaces, ab initio calculations are becoming increasingly more accurate, albeit still very expensive computationally. Neural network-based approaches like the present one, alone or by exploiting the promotion to the DFT level of choice,41 open up the possibility to explore with superior statistical accuracy systems that were, until now, almost exclusively in the realm of empirical potential-based simulations.

See the supplementary material for the neural network training force histogram and the orientation free energy maps.

The computational results presented have been achieved using the Vienna Scientific Cluster (VSC).

We thank Andreas Singraber for providing help with the n2p2 package.

The datasets of the neural network training configurations with forces, the optimal weights, as well as the input files are openly available in Zenodo at http://doi.org/10.5281/zenodo.3944892, reference number [3944892].

1.
R.
Kumar
,
J. R.
Schmidt
, and
J. L.
Skinner
,
J. Chem. Phys.
126
,
204107
(
2007
).
2.
J. R.
Errington
and
P. G.
Debenedetti
,
Nature
409
,
318
(
2001
).
3.
M.
Sovago
,
R.
Kramer Campen
,
H. J.
Bakker
, and
M.
Bonn
,
Chem. Phys. Lett.
470
,
7
(
2009
).
4.
A.
Morita
and
J. T.
Hynes
,
Chem. Phys.
258
,
371
(
2000
).
5.
L. F.
Scatena
,
Science
292
,
908
(
2001
).
6.
T. D.
Kühne
,
T. A.
Pascal
,
E.
Kaxiras
, and
Y.
Jung
,
J. Phys. Chem. Lett.
2
,
105
(
2011
).
7.
D.
Chandler
,
Nature
437
,
640
(
2005
).
8.
A.
Braslau
,
M.
Deutsch
,
P. S.
Pershan
,
A. H.
Weiss
,
J.
Als-Nielsen
, and
J.
Bohr
,
Phys. Rev. Lett.
54
,
114
(
1985
).
9.
Y. R.
Shen
,
Nature
337
,
519
(
1989
).
10.
M. J.
Gillan
,
D.
Alfè
, and
A.
Michaelides
,
J. Chem. Phys.
144
,
130901
(
2016
).
11.
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
,
Proc. Natl. Acad. Sci. U. S. A.
114
,
10846
(
2017
).
12.
A. P.
Bartók
,
S.
De
,
C.
Poelking
,
N.
Bernstein
,
J. R.
Kermode
,
G.
Csányi
, and
M.
Ceriotti
,
Sci. Adv.
3
,
e1701816
(
2017
).
13.
T.
Morawietz
,
A.
Singraber
,
C.
Dellago
, and
J.
Behler
,
Proc. Natl. Acad. Sci. U. S. A.
113
,
8368
(
2016
).
14.
J.
Behler
and
M.
Parrinello
,
Phys. Rev. Lett.
98
,
146401
(
2007
).
15.
T.
Morawietz
and
J.
Behler
, HDNNP training data set for H2O,
2019
, type: dataset.
16.
M. P.
Gelfand
and
M. E.
Fisher
,
Physica A
166
,
1
(
1990
).
17.
B.
Hammer
,
L. B.
Hansen
, and
J. K.
Nørskov
,
Phys. Rev. B
59
,
7413
(
1999
).
18.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
,
J. Chem. Phys.
132
,
154104
(
2010
).
19.
K.
Toukan
and
A.
Rahman
,
Phys. Rev. B
31
,
2643
(
1985
).
20.
V.
Blum
,
R.
Gehrke
,
F.
Hanke
,
P.
Havu
,
V.
Havu
,
X.
Ren
,
K.
Reuter
, and
M.
Scheffler
,
Comput. Phys. Commun.
180
,
2175
(
2009
).
21.
P.
Schienbein
and
D.
Marx
,
J. Phys. Chem. B
122
,
3318
(
2018
).
22.
M.
Sega
and
C.
Dellago
,
J. Phys. Chem. B
121
,
3798
(
2017
).
23.
N. B.
Wilding
,
Phys. Rev. E
52
,
602
(
1995
).
24.
W.
Wagner
and
A.
Pruß
,
J. Phys. Chem. Ref. Data
31
,
387
(
2002
).
25.
A. Z.
Panagiotopoulos
,
Mol. Phys.
61
,
813
(
1987
).
26.
M. J.
McGrath
,
J. I.
Siepmann
,
I.-F. W.
Kuo
,
C. J.
Mundy
,
J.
VandeVondele
,
J.
Hutter
,
F.
Mohamed
, and
M.
Krack
,
J. Phys. Chem. A
110
,
640
(
2006
).
27.
K. K.
Mon
and
K.
Binder
,
J. Chem. Phys.
96
,
6989
(
1992
).
28.
N. B.
Vargaftik
,
B. N.
Volkov
, and
L. D.
Voljak
,
J. Phys. Chem. Ref. Data
12
,
817
(
1983
).
29.
J. L. F.
Abascal
and
C.
Vega
,
J. Chem. Phys.
123
,
234505
(
2005
).
30.
M.
Sega
,
G.
Hantal
,
B.
Fábián
, and
P.
Jedlovszky
,
J. Comput. Chem.
39
,
2118
(
2018
).
31.
P.
Jedlovszky
,
Á.
Vincze
, and
G.
Horvai
,
Phys. Chem. Chem. Phys.
6
,
1874
(
2004
).
32.
E.
Chacón
and
P.
Tarazona
,
Phys. Rev. Lett.
91
,
166103
(
2003
).
33.
Fundamentals of Inhomogeneous Fluids
, edited by
D.
Henderson
(
Dekker
,
New York
,
2009
).
34.
S.
Fisk
and
B.
Widom
,
J. Chem. Phys.
50
,
3219
(
1969
).
35.
A.
Singraber
,
J.
Behler
, and
C.
Dellago
,
J. Chem. Theory Comput.
15
,
1827
(
2019
).
36.
S.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
).
37.
N.
Michaud-Agrawal
,
E. J.
Denning
,
T. B.
Woolf
, and
O.
Beckstein
,
J. Comput. Chem.
32
,
2319
(
2011
).
38.
L. B.
Pártay
,
G.
Hantal
,
P.
Jedlovszky
,
Á.
Vincze
, and
G.
Horvai
,
J. Comput. Chem.
29
,
945
(
2008
).
39.
M.
Ester
,
H.-P.
Kriegel
,
J.
Sander
,
X.
Xu
 et al, in
Proceedings of the 2nd International Conference on Knowledge Discovery and Data Mining
(
AAAI Press
,
Portland, OR
,
1996
), Vol. 96, pp.
226
231
.
40.
M.
Sega
and
G.
Hantal
,
Phys. Chem. Chem. Phys.
19
,
18968
(
2017
).
41.
B.
Cheng
,
E. A.
Engel
,
J.
Behler
,
C.
Dellago
, and
M.
Ceriotti
,
Proc. Natl. Acad. Sci. U. S. A.
116
,
1110
(
2019
).

Supplementary Material