The hydrophobic interaction, often combined with the hydrophilic or ionic interactions, makes the behavior of aqueous solutions very rich and plays an important role in biological systems. Theoretical and computer simulation studies have shown that the water-mediated force depends strongly on the size and other chemical properties of the solute, but how it changes with these factors remains unclear. We report here a computer simulation study that illustrates how the hydrophobic pair interaction and the entropic and enthalpic terms change with the solute size when the solute–solvent weak attractive interaction is unchanged with the solute size. The nature of the hydrophobic interaction changes qualitatively as the solute size increases from that of methane to that of fullerene. The potential of mean force between small solutes has several well-defined extrema, including the third minimum, whereas the potential of mean force between large solutes has the deep contact minimum and the large free-energy barrier between the contact and the water-bilayer separated configurations. The difference in the potential of mean force is related to the differences in the water density, energy, and hydrogen bond number distributions in the vicinity of the pairs of hydrophobic solutes.

The effective interactions between solute molecules in water or in aqueous solutions in many cases differ qualitatively from the corresponding direct interactions in vacuum due to the solvent-induced part of the interactions. More important than the mere difference from the direct interaction, however, is that the effective interaction depends in a complex way on the temperature, pressure, and composition of an aqueous solution, as well as the size of the solute and the solute–solvent interaction.

The hydrophobic interaction is one of the most studied effective interactions as it is the driving force for molecular self-assembly, such as the formation of micelles, membranes, and vesicles in aqueous solutions.1–6 Small hydrophobic molecules such as methane and inert gases are less attracted (or even repelled) in water than in a vacuum at low temperatures, but more attracted than in a vacuum at high temperatures.7 The pressure and salt-concentration dependences of the hydrophobic interaction are also similar in trend to the temperature dependence.8 The character of the hydrophobic hydration changes with the size of the solute9 and so would the nature of the hydrophobic interaction, including its temperature, pressure, and salt concentration dependencies.

In fact, the effect of the solute size on the strength of the hydrophobic interaction has been one of the issues that has attracted attention.10–28 Lum, Chandler, and Weeks showed that the solvation free energy and the effective interaction of hard-sphere solutes in water exhibit the crossover occurring on nanometer length scales.11 Zangi quantified the strength of hydrophobic interaction between graphene sheets with different sizes, showing that the association of small graphene sheets has an entropic origin while that of large graphene sheets is mainly due to the solute–solute direct pair interaction.23 Scheraga et al. calculated the potential w(r) of the mean force for hydrocarbons of different sizes in water,16,17 the results of which indicate that w(r) at the contact distance becomes higher than the potential in vacuum as the size of the solute exceeds that of neopentane. In a broader context, earlier studies based on the theory of liquids have examined the solute-size dependence of the solvent-mediated interactions in simple liquids.15,29–39 The numerical results from the integral equation method29 and molecular dynamics (MD) simulations27,28 indicate that the effective pair interactions between Lennard-Jones (LJ) particles in LJ liquids become less attractive with the size of the solute, provided that the LJ energy parameter for the solute–solvent pair is fixed to a value greater than or equal to the LJ energy parameter for the solvent–solvent pair.27–29 

The strength of the effective solute–solute pair interactions in a solvent can be quantified by the osmotic second virial coefficient B in the virial expansion of the osmotic pressure at constant chemical potential of the solvent species. The effective interaction is attractive overall when B is negative and repulsive otherwise. B is related to w(r) or the solute–solute radial distribution function g(r) via
(1)
where ρ is the solute density, k is Boltzmann’s constant, T is the temperature, and is an infinitesimal volume element.40 Recent simulation studies have found that B, which is negative, is proportional to some power α of the solute diameter, and that α ≃ 6 or 7, indicating an extremely strong size dependence for LJ solutes in water.27,28

Aqueous solutions of water-soluble polymers such as PNIPAM and PEO are homogeneous at room temperature but undergo phase separation at higher temperatures, indicating that the effective interactions between hydrophobic moieties in those polymers are entropic forces. It has been well-verified by simulations and theoretical methods that the strength of the hydrophobic interaction increases with temperature.7,8,10,18,20,22,27,41–58 Although the hydrophobic interaction is overall an entropic force, but there is a significant enthalpic (energetic) contribution, which is opposite in sign to the entropic one.14,18,20,22,23,41–52,56,59–63 To fully understand the solute size dependence of w(r), it is essential to investigate how the fraction of the entropic contribution to w(r) changes with solute size.

The aim of the present work is to quantify the enthalpic and entropic contributions to w(r) for solutes of different sizes and thereby to reveal how the underlying molecular mechanisms in the hydrophobic interaction change with solute size. First, we calculate, with high accuracy, w(r) and B for the LJ particles with different diameters σ in water at several temperatures. The enthalpic and entropic contributions to w(r) are obtained from the temperature derivative of the solvent-induced part w*(r). Second, we analyze the microscopic structures of water molecules around two solutes when the solute–solute distance r is fixed to some characteristic distances.

Isobaric–isothermal molecular dynamics (MD) simulations were performed for the model aqueous solutions consisting of TIP4P/2005 water molecules64 and spherical solute particles interacting via the Lennard-Jones (LJ) potential,
(2)
where ϵ is the depth of the pair potential well and σ is the particle diameter. The reference solute is chosen to be methane, whose LJ parameters are σm = 0.373 nm and ϵm = 1.23 kJ mol−1 in the TraPPE-UA force field.65 The reduced LJ diameters for solutes are σ* ≡ σ/σm = 1, 1.5, 2, 2.5, and 3; the solute–solute LJ energy parameter ϵ is fixed to ϵm. The solute–water pair potential is the LJ potential whose size and energy parameters are given by the Lorentz–Berthelot combining rules.

We used GROMACS 2018 software66 to perform MD simulations for the aqueous solutions under three-dimensional periodic boundary conditions. The pressure is set to 1 bar using the Parrinello–Rahman method, and the temperature is fixed at 270, 300, 330, or 360 K by using the Nosé–Hoover method. The equilibrium trajectories for T = 300 K are those obtained earlier.28 The time step interval is 1 fs, and the configurations of solute particles were recorded every 0.05 ps.

The model aqueous solutions consist of water molecules and solute particles. For the solutes with smaller diameters σ* = 1 and 1.5, the number N of solute particles is 40 (for σ* = 1) or 20 (for σ* = 1.5) and the number Nw of water molecules is 4000. The potential w(r) of mean force for each system was computed via the solute–solute g(r) obtained from the MD trajectories. The duration time t for the production run is 100 ns for σ* = 1 and 200 ns for σ* = 1.5. For the solutes with larger diameters σ* = 2, 2.5, and 3, N = 2 and Nw = 8000. For these systems, we performed the umbrella sampling simulations67,68 to compute w(r). The simulation time t for each window is 20 ns.

In the MD simulations for the solutes with σ* = 1 and 1.5, the solute–solute LJ pair potential ϕ(r) was replaced with a repulsive potential ϕrep(r) in order to suppress multi-body aggregations and minimize the finite-concentration effect on g(r). This procedure is particularly important to evaluate the osmotic second virial coefficient B from Eq. (1). The resulting radial distribution function grep(r) is then converted to g(r) by
(3)
with ϕatt(r) = ϕ(r) − ϕrep(r). For σ* = 1, ϕrep(r) is the repulsive part of the Weeks–Chandler–Andersen (WCA) potential.69 The solute particles with σ* = 1.5 tend to aggregate in water more strongly than those with σ* = 1. Therefore, ϕrep(r) is taken to be 4ϵ(σ/r).12 

The cutoff distance rcut for the LJ potentials depends on the diameter of the solute particle in each aqueous solution: rcut = 1.3 nm for σ* = 1 and 1.5, rcut = 2 nm for σ* = 2 and 2.5, and rcut = 2.4 nm for σ* = 3. The Coulomb potentials were treated using the particle mesh Ewald method, with the real-space cutoff distance being the same as rcut for the LJ potentials.

In the umbrella sampling simulations, the solute–solute separations were constrained via the harmonic potential with the spring constant set to 1000 kJ mol−1 nm−2. The constraint distance ranges from 0.6 to 2.9 nm in 0.1 nm increments for the solute with σ* = 2, from 0.7 to 2.8 nm for σ* = 2.5, and from 0.8 to 3.1 nm for σ* = 3. The potentials w(r) of mean force were obtained from the umbrella sampling simulations for the windows using the weighted histogram analysis method.70,71 The corresponding radial distribution functions g(r) are given by exp[−w(r)/kT].

The osmotic second virial coefficient B is evaluated from g(r) via Eq. (1). However, one cannot use g(r) directly obtained from a simulation since in a closed system, it does not converge to 1 and the volume integral of g(r) − 1 would tend to diverge. To overcome this issue, we employed the following computational procedure.27,28 First, the entire function g(r) was scaled by a factor adjusting the average value of g(r) over a certain range of large r to 1.7 Second, the Kirkwood–Buff (KB) integral G is evaluated from the following equations:72,73
(4)
(5)
where L is the upper limit of the integral with respect to r and D is a constant. G(L) is obtained as a function of L using Eq. (4) with corrected g(r). G is then determined as a slope in the plot of G(L)L vs L over a certain range where G(L)L is linear to L.
The solute–solute effective potential w(r) consists of the direct part ϕ(r) and the indirect, solvent-induced part w*(r) = w(r) − ϕ(r). The latter is the difference between the solvation free energy μpair*(r) of a pair of solutes distance r apart and that of a pair infinitely far apart. The temperature dependence of w(r) comes from that of w*(r) alone as ϕ(r) is independent of T. The temperature derivative of w*(r) at fixed pressure gives the enthalpic and entropic contributions to w(r),
(6)
where Δ means the constant-pressure process of changing the distance between two solute particles from infinity to r; h*(r) and s*(r) are the solvation enthalpy and entropy of the pair at fixed pressure, respectively. Here, the solvation enthalpy and entropy at fixed pressure mean the enthalpy and entropy changes upon inserting a molecule (now a pair of solute particles) in the solvent at a given pressure with its position and orientation fixed.74 
In the present study, we evaluated w*(r) at four temperatures (270, 300, 330, and 360 K) in order to fit the data to a quadratic function of T,
(7)
Then, we evaluated Δs*(r) and Δh*(r) at T = 300 K via
(8)

We performed additional isobaric–isothermal MD simulations for the aqueous solutions containing a pair of LJ solutes with fixed separation distance to examine the microscopic structures of water molecules around the pair. The solute diameters are σ* = 1, 2, and 3. In those simulations, the solute–solute distance R was fixed to the values corresponding to the first minimum, the first maximum, the second minimum, the second maximum, and the third minimum of w(r). Table I presents the values of R for σ* = 1, 2, and 3. The pressure and temperature are maintained at 1 bar and 300 K, respectively. The duration time of the production run is 20 ns for each simulation, and the configurations of solute particles and water molecules were recorded every 1 ps. We also investigate the solute with σ* = 3 interacting with water molecules via the repulsive WCA potential.

TABLE I.

Fixed solute–solute distances R corresponding to the first minimum, the first maximum, the second minimum, the second maximum, and the third minimum of the potentials w(r) of mean force for pairs of LJ solutes with σ* = 1, 2, and 3 at T = 300 K.

R (nm)
σ*First minimumFirst maximumSecond minimumSecond maximumThird minimum
0.388 0.562 0.706 0.874 1.026 
0.724 0.992 1.114 1.230 1.414 
1.032 1.400 1.506 1.570 1.764 
R (nm)
σ*First minimumFirst maximumSecond minimumSecond maximumThird minimum
0.388 0.562 0.706 0.874 1.026 
0.724 0.992 1.114 1.230 1.414 
1.032 1.400 1.506 1.570 1.764 

Figure 1 is the schematic illustration of the cylindrical coordinate system (d, θ, h). Based on this system, we calculated the distribution gcyl(h, d) of water molecules (oxygen), which is normalized to unity in the bulk region, the water–water pair interaction energy Eww(h, d), and the number NHB(h, d) of hydrogen bonds, with all the quantities being averaged over θ and being depicted on the dh plane. A pair of water molecules is taken into account for computing Eww when the distance rOO between their oxygen atoms is less than 0.35 nm. Two water molecules form one hydrogen bond if rOO < 0.35 nm and the H–O⋯O angle is less than or equal to 30°.75–77 We note that the cylindrical distributions of water molecules around two solute particles and other quantities were evaluated in earlier studies.14,16,17,19,49,51 Our focus here is on the question of how these distributions changes with the solute particle size.

FIG. 1.

Schematic illustration of water molecules around two solute particles in the cylindrical coordinate system (d, θ, h). The axis h passes through two centers of solutes (orange circles), d is the distance from the axis h, and θ is the azimuth. The origin of the coordinate system is set to the midpoint of the two centers of solutes. The red circles represent oxygen atoms of water molecules. The equilibrium distributions of oxygen atoms averaged over θ are depicted on the hd plane.

FIG. 1.

Schematic illustration of water molecules around two solute particles in the cylindrical coordinate system (d, θ, h). The axis h passes through two centers of solutes (orange circles), d is the distance from the axis h, and θ is the azimuth. The origin of the coordinate system is set to the midpoint of the two centers of solutes. The red circles represent oxygen atoms of water molecules. The equilibrium distributions of oxygen atoms averaged over θ are depicted on the hd plane.

Close modal

First, we examine the size and temperature effects on the solute–solute effective interactions in water. Figure 2(a) shows the potentials w(r) of mean force for pairs of LJ particles with different diameters in water at T = 300 and 360 K. The reduced solute diameter σ* ranges from 1 (0.373 nm, methane size) to 3 (1.119 nm, C60 size). At both temperatures, the first minimum of w(r) descends with increasing σ*. In contrast, both the first maximum and the second minimum of w(r) ascend with σ*, and the latter disappears at σ* greater than 2, being a shoulder. These features at ambient temperature have been reported earlier28 and now are confirmed at 360 K. One also notices that the curve of w(r) shifts downward with increasing T.

FIG. 2.

(a) Potentials w(r) of mean force for pairs of solutes with various sizes in water at different temperatures T. The LJ diameters of the solute particles are σ* ≡ σ/σm = 1, 1.5, 2, 2.5, and 3 with σm = 0.373 nm. The solute–solute LJ energy parameter ϵ is fixed to ϵm = 1.23 kJ mol−1 for all-size solutes. The solid and dotted curves are w(r) at T = 300 and 360 K, respectively. (b) The log–log plot of the osmotic second virial coefficients B at T = 300 K against σ*. The dotted and solid blue lines are the best fits to the data with the sixth and seventh power law, respectively. The second virial coefficient Bgas for the LJ gas is also plotted as a function of σ*. (c) B and Bgas at T = 360 K. The best fits of the sixth (dotted red line) and seventh (solid red line) power law are also plotted.

FIG. 2.

(a) Potentials w(r) of mean force for pairs of solutes with various sizes in water at different temperatures T. The LJ diameters of the solute particles are σ* ≡ σ/σm = 1, 1.5, 2, 2.5, and 3 with σm = 0.373 nm. The solute–solute LJ energy parameter ϵ is fixed to ϵm = 1.23 kJ mol−1 for all-size solutes. The solid and dotted curves are w(r) at T = 300 and 360 K, respectively. (b) The log–log plot of the osmotic second virial coefficients B at T = 300 K against σ*. The dotted and solid blue lines are the best fits to the data with the sixth and seventh power law, respectively. The second virial coefficient Bgas for the LJ gas is also plotted as a function of σ*. (c) B and Bgas at T = 360 K. The best fits of the sixth (dotted red line) and seventh (solid red line) power law are also plotted.

Close modal

It should be noted that w(r) for the C60-sized solutes is minimal at the distance smaller than the LJ diameter σ (Table I). This is because the water-induced attractive force between the C60-sized solutes is so strong that it balances with the equally strong repulsive force due to the direct solute–solute LJ pair potential at the short distance.

Figure 2(b) shows the plot of the osmotic second virial coefficients B at T = 300 K as a function of σ* in a log–log scale. The osmotic B is negative and that magnitude increases with increasing σ. This result indicates that the effective interaction between solute particles in water is attractive and becomes stronger as the particle size increases. The second virial coefficients Bgas for the LJ gas, which quantify the strength of the direct pair interaction in a vacuum, are also shown in Fig. 2(b). Both B and Bgas decrease with σ, but the magnitude of B is much greater than that of Bgas in the σ* range over 1.5.

The log–log plot of B vs σ* exhibits a linear relationship, indicating that
(9)
The best fits of the data to Eq. (9) with α = 6 and 7 are plotted by the dotted and solid blue lines, respectively. One may see that α = 7 fits the data better than 6. The gas virial coefficient Bgas is proportional to the cubic of the particle diameter, so the contribution from the water-induced part w*(r) = w(r) − ϕ(r) to B has a stronger size dependence than that from the direct part ϕ(r). Earlier studies27,28 showed that the power law between B and σ is partially understood based on the thermodynamic identity78,79 for B.

Figure 2(c) shows B and Bgas at T = 360 K as a function of σ*. At this temperature, the magnitude of B is greater than that of Bgas for any size of solute particles. The log–log plot of B vs σ* again suggests the power-law dependence of B on σ*. The best estimate of α at 360 K is 6 rather than 7.

Comparison between the plots shown in Fig. 2(b) and those shown in Fig. 2(c) indicates that B decreases (becomes more negative) with increasing T for any given solute, i.e., the effective solute–solute pair attractive force between hydrophobic particles in water becomes stronger at higher temperatures, which is characteristic of the hydrophobic interaction. On the other hand, Bgas increases with T, which is known for real gases and for model pair potentials including the LJ potential.

Now, we examine the distribution gcyl(h, d) of water molecules around a pair of solute particles with the inter-particle distance r fixed to particular values. Figure 3(a) shows gcyl(h, d) for r of the third minimum of w(r) (see Table I for the values of r). For all the diameters of solute particles, the first solvation shells (bright-colored rings) around the two particles have little effect on each other even in between the two particles. The local density of water at the solvation shell is much higher than the bulk density: gcyl(h, d) is ∼1.9, 1.8, and 1.6 for σ* = 1, 2, and 3, respectively. At this separation distance, gcyl(h, d) has two peaks in between the two solutes, indicating the bilayer structure of water.

FIG. 3.

Normalized two-dimensional distributions gcyl(h, d) of water molecules around a pair of solutes at T = 300 K. Two solute particles with LJ diameters σ* = 1(left), 2(middle), and 3(right) are indicated in black. The solute–solute separation distance R is fixed to that of (a) the third minimum, (b) the second maximum, (c) the second minimum, (d) the first maximum, and (e) the first minimum of w(r) for each solute (see Table I).

FIG. 3.

Normalized two-dimensional distributions gcyl(h, d) of water molecules around a pair of solutes at T = 300 K. Two solute particles with LJ diameters σ* = 1(left), 2(middle), and 3(right) are indicated in black. The solute–solute separation distance R is fixed to that of (a) the third minimum, (b) the second maximum, (c) the second minimum, (d) the first maximum, and (e) the first minimum of w(r) for each solute (see Table I).

Close modal

The bilayer structure disappears at r set to the second maximum of w(r): the two solvation shells partly overlap with each other in the gap between the two spherical particles and gcyl in the gap is lower than outside for σ* = 1 and 2 [Fig. 3(b)]. However, gcyl in the gap is higher than outside when the separation distance is fixed to that of the second minimum of w(r) [Fig. 3(c)]: gcyl(0, 0) is ∼2.4 and 2.1 for σ* = 1 and 2.

When the inter-particle distance is of the first maximum of w(r) [Fig. 3(d)], the distribution gcyl(0, 0) at the midpoint of the pair is ∼0, 0.4, and 0.6 for σ* = 1, 2 and 3, i.e., water molecules are scarce in the gap between the two particles. Finally, when the two solute particles are in contact [Fig. 3(e)], the distribution at the intersection of the two first solvation shells is greater than at the other part of the same solvation shells: see the brighter spots in the case of σ* = 1 and 2.

Next, we examine how the enthalpic and entropic contributions to w*(r) change with the solute size. Figure 4 shows curves of w*(r), Δh*(r), and Δs*(r) for σ* = 1, 1.5, 2, 2.5, and 3 at T = 300 K. The circle, square, diamond, and triangle on each curve indicates the distances of the first minimum (contact minimum), the first maximum, the second minimum (monolayer-separated minimum), and the third minimum (bilayer-separated minimum) of the potential w(r) of mean force.

FIG. 4.

Water-mediated part w*(r) of the hydrophobic interaction and the enthalpic and entropic contributions, Δh*(r) and −TΔs*(r), to w*(r) for the solutes with the diameters σ* = 1, 1.5, 2, 2.5, and 3 at T = 300 K: (a) w*(r) (solid curve) and the direct pair potential ϕ(r) (dotted curve); (b) Δh*(r); and (c) −TΔs*(r). The open circles, squares, diamonds, and triangles indicate the separation distances of the first minimum, the first maximum, the second minimum (monolayer-separated minimum), and the third minimum (bilayer-separated minimum) of w(r).

FIG. 4.

Water-mediated part w*(r) of the hydrophobic interaction and the enthalpic and entropic contributions, Δh*(r) and −TΔs*(r), to w*(r) for the solutes with the diameters σ* = 1, 1.5, 2, 2.5, and 3 at T = 300 K: (a) w*(r) (solid curve) and the direct pair potential ϕ(r) (dotted curve); (b) Δh*(r); and (c) −TΔs*(r). The open circles, squares, diamonds, and triangles indicate the separation distances of the first minimum, the first maximum, the second minimum (monolayer-separated minimum), and the third minimum (bilayer-separated minimum) of w(r).

Close modal

One can see in Fig. 4(a) that with increasing σ* from 1 to 3, the free-energy barrier of w*(r), which corresponds to the first maximum of w(r), rises sharply from 1 to 7 kJ/mol. At the same time, the local minimum of w*(r), which corresponds to the second minimum of w(r), also goes up with σ* and disappears at σ* = 2.5. Comparing Fig. 4(b) with Fig. 4(a), one notices that the major peak of Δh*(r) corresponds to that in w*(r) at a distance close to the first maximum of w(r). Thus, the remarkable evolution of the free-energy barrier with increasing solute size is due to the increase in the unfavorable solvation enthalpy. At the solute–solute contact distance, Δh*(r) is positive for σ* up to 2.5 but nearly zero for σ* = 3 [open circles in Fig. 4(b)], while −TΔs*(r) is always negative and the magnitude is much greater than Δh*(r) [open circles in Fig. 4(c)]. The minimum of w(r) at the solute–solute contact distance goes down with increasing solute diameter because the entropic contribution [−TΔs*(r)] becomes more negative. This statement is valid for σ* ≤2.5; we will discuss a different cause for σ* = 3 and beyond.

Plotted in Fig. 5 are values of w*(r), Δh*(r), and −TΔs*(r) at fixed distances r as functions of σ*. At the solute–solute contact distance [Fig. 5(a)], w*(r)(<0) decreases monotonically with σ*; −TΔs*(r) decreases up to σ* = 2.5 and turns to increase slightly; and Δh*(r) is maximal at around σ* = 2.

FIG. 5.

Solute size dependence of w*(r), Δh*(r), and −TΔs*(r) at r of (a) the first minimum, (b) the first maximum, (c) the second minimum, and (d) the third minimum of w(r). The solute diameter σ* ranges from 1 to 3. The blue circles, red diamonds, and green squares represent w*(r), Δh*(r), and −TΔs*(r), respectively.

FIG. 5.

Solute size dependence of w*(r), Δh*(r), and −TΔs*(r) at r of (a) the first minimum, (b) the first maximum, (c) the second minimum, and (d) the third minimum of w(r). The solute diameter σ* ranges from 1 to 3. The blue circles, red diamonds, and green squares represent w*(r), Δh*(r), and −TΔs*(r), respectively.

Close modal

For the methane-sized LJ particles in contact with each other, the magnitude of the favorable entropy change Δs*(r)/k (≃2.72) exceeds that of the unfavorable entropy change Δh*(r)/kT (≃1.95) to such an extent that the water-induced pair potential w*(r)/kT(≃−0.77) is slightly negative. On the other hand, for the C60-sized LJ particles in the contact configuration, Δh*(r)/kT (≃−0.06) is close to zero and thus w*(r)/kT (≃−6.69) is essentially equal to −Δs*(r)/k.

At distances corresponding to the first maximum and second minimum of w(r) [Figs. 5(b) and 5(c)], for the most part, w*(r) and Δh*(r), both positive, increase with σ* while −TΔs*(r) (<0) decreases with σ*. At the distance of the third minimum of w(r) or the bilayer-separated distance [Fig. 5(d)], w*(r) is near zero in the whole range of σ*. However, −TΔs*(r) (>0) increases with σ* while Δh*(r) (<0) decreases. For the largest solutes of σ* = 3, therefore, there is a large cancellation of −TΔs*(r) = 7.14 kJ mol−1 and Δh*(r) = −6.42 kJ mol−1, resulting in w*(r) = 0.72 kJ mol−1.

Having found the solvation-enthalpy and solvation-entropy contributions to w*(r), we now extract microscopic information from the spatial distributions of the potential energy Eww and the number NHB of hydrogen bonds, both per water molecule, around the pair of the hydrophobic solute particles in each system.

The distribution Eww(h, d) is the spatial distribution of the potential energy for a single water molecule due to water–water pair interactions within a cutoff distance of 0.35 nm. Figures 6(a)6(e) show 3 × 5 distributions: the three sized solutes (σ* = 1, 2, and 3) and the five distances r. When r is of the third minimum of w(r) [Fig. 6(a)], the energy distribution in the vicinity of one spherical particle is unaffected by the other particle. For r equal to or smaller than that of the second maximum of w(r) [Figs. 6(b)6(e)], Eww in between the solute particles and in the neighborhood of the intersection circle of the two first-solvation shells around the two particles is higher than Eww in the other region.

FIG. 6.

Spatial distributions of the potential energy Eww of a water molecule due to the water–water interaction around the fixed pair of solute particles at T = 300 K. Two solutes with LJ diameters σ* = 1(left), 2(middle), and 3(right) are represented in white. The solute–solute separation distances are fixed to those of (a) the third minimum, (b) the second maximum, (c) the second minimum, (d) the first maximum, and (e) the first minimum of w(r) for each solute (see Table I). Eww of the bulk water is −73.7 kJ mol−1.

FIG. 6.

Spatial distributions of the potential energy Eww of a water molecule due to the water–water interaction around the fixed pair of solute particles at T = 300 K. Two solutes with LJ diameters σ* = 1(left), 2(middle), and 3(right) are represented in white. The solute–solute separation distances are fixed to those of (a) the third minimum, (b) the second maximum, (c) the second minimum, (d) the first maximum, and (e) the first minimum of w(r) for each solute (see Table I). Eww of the bulk water is −73.7 kJ mol−1.

Close modal

Comparing three heat maps in each row of Fig. 6, one finds that the larger the solute particles, Eww in the vicinity of each particle is higher than Eww in the bulk. In Fig. 6(a), for example, Eww in the first solvation shells of the solute particles of σ* = 1 is almost the same as the bulk value −73.7 kJ mol−1, but the corresponding Eww is −65 and −60 kJ mol−1 for the larger solute particles of diameters σ* = 2 and 3, respectively. Water molecules in the first solvation shell are little perturbed by the methane-sized solute but are more frustrated by the larger solutes. We will see the equivalent picture in terms of the hydrogen bonding.

Except for σ* = 1, Eww between the two spherical solutes is higher than in bulk water when the solute–solute distance r is of the second maximum of w(r) or shorter: in Fig. 6(b), Eww(0, 0) at the midpoint is higher by 9 and 14 kJ mol−1 than bulk water for σ* = 2 and 3, respectively, and Eww(0, 0) increases as r decreases further [Figs. 6(c) and 6(d)].

The spatial distribution of the average number NHB of hydrogen bonds per water molecule are shown in Fig. 7. At the bilayer-separated distance, i.e., r of the third minimum of w(r) [Fig. 7(a)], NHB in the solvation shell of the small solute (σ* = 1) is almost the same as NHB = 3.66 in the bulk water, while NHB ≈ 3.4 and 3.0 in the solvation shells of the solutes with σ* = 2 and 3, respectively. In accordance with what we saw for Eww, we notice that the hydrogen-bond network is hardly disturbed by the methane-sized solutes, whereas it is strongly destabilized by the large particles. This result illustrates the difference in the hydrophobic hydration of small and large solute molecules—the idea expressed by Stillinger80 and theoretically confirmed by Chandler et al.9,11

FIG. 7.

Numbers NHB of hydrogen bonds between water molecules per one molecule at T = 300 K. Two solutes with LJ diameters σ* = 1(left), 2(middle), and 3(right) are represented in white. The solute–solute separations R are fixed to r of (a) the third minimum, (b) the second maximum, (c) the second minimum, (d) the first maximum, and (e) the first minimum of w(r) for each solute (see Table I). NHB of the bulk water is 3.66.

FIG. 7.

Numbers NHB of hydrogen bonds between water molecules per one molecule at T = 300 K. Two solutes with LJ diameters σ* = 1(left), 2(middle), and 3(right) are represented in white. The solute–solute separations R are fixed to r of (a) the third minimum, (b) the second maximum, (c) the second minimum, (d) the first maximum, and (e) the first minimum of w(r) for each solute (see Table I). NHB of the bulk water is 3.66.

Close modal

When the distance between two particles is of the second maximum of w(r) [Fig. 7(b)], there appears the effect of confinement on water: the effect is nearly absent for σ* = 1, slightly visible for σ* = 2 and clearly notable for σ* = 3. As the distance reduces to that of the second minimum of w(r) [Fig. 7(c)] and further down to that of the first maximum of w(r) [Fig. 7(d)], NHB in the region between two particles becomes lower and lower. When the two particles of σ* = 3 are in contact with each other [Fig. 7(e)], NHB at the intersection of the two solvation shells (the black regions at the cusps in the figure) is close to 2.4 compared to the bulk value 3.66. The set of the results for NHB is all consistent with that for Eww.

Earlier studies showed that the structure of water confined between hydrophobic planar walls changes from bilayer-like to monolayer-like structure as two walls approach each other.81–84  Figures 3, 6, and 7 show similar trends: when the solute–solute separation is of the third minimum of w(r), the bilayer-like structure is present in between the two particles of σ* = 3. At the separation distance of the second minimum of w(r), confined water assumes the monolayer-like structure. It is noteworthy that NHB ≃ 3.6 in the bilayer-like water is only slightly less than that in bulk water [the right panel shown in Fig. 7(a)] while NHB ≃ 2.9 in the monolayer-like structure is much less that the bulk value [the right panel shown in Fig. 7(c)]. This result is consistent with the earlier MD simulation result that the melting temperature of bilayer ice is much higher than that of monolayer ice:82 water molecules can form four-coordinated network structure without substantial energetic penalty if there are two molecular layers or more in confining geometry.

Here, we discuss, at the molecular level, the difference in hydrophobic interactions of small and large particles in the contact-minimum configuration.

First, let us consider the positive entropy change Δs* in the contact-minimum configuration [Fig. 5(a)]. The positive Δs* cannot be explained by an increase in the rotational entropy of water molecules because the hydrogen bond number per molecule is lower in the first hydration shell than in the bulk water. This is due to an increase in the translational entropy of solvent molecules or, equivalently, the excluded volume effect, as shown schematically in Fig. 8(a). Since the excluded volume effect is present in all liquids, one observes that Δs* > 0 for hard sphere solutes in a hard sphere solvent, LJ solutes in an LJ solvent, and hydrophobic solutes in water. A subtle issue is the solute size dependence of Δs*, as shown in Fig. 5(a): Δs* increases with the solute diameter and peaks out at about σ* = 2.5. The fact that Δs* does not increase monotonically with solute diameter may be due to a unique property of hydrophobic solutes in water.

FIG. 8.

Changes in the solvation shells and the solvation enthalpy when two solute particles are brought into contact. (a) Schematic picture of the first solvation shells around a pair of large solute particles in the separated and contact configuration. The water molecules in the solvation shells have Eww higher than in bulk (Fig. 6) and NHB less than in bulk (Fig. 7), but their solute–water interaction energy is lower than in bulk. (b) The solvation enthalpy change Δh*(r) and its components ΔUww(r) and ΔUws(r).

FIG. 8.

Changes in the solvation shells and the solvation enthalpy when two solute particles are brought into contact. (a) Schematic picture of the first solvation shells around a pair of large solute particles in the separated and contact configuration. The water molecules in the solvation shells have Eww higher than in bulk (Fig. 6) and NHB less than in bulk (Fig. 7), but their solute–water interaction energy is lower than in bulk. (b) The solvation enthalpy change Δh*(r) and its components ΔUww(r) and ΔUws(r).

Close modal
Second, we recall that Δh*(r) = h*(r) − h*(), where h*(r) is the solvation enthalpy of the pair of solute particles and is essentially the excess of the configurational energy U of the system with the pair over that without, for the pressure-volume term is negligible. The configurational energy U consists of Uww the sum of the pair interactions between water molecules and Uws the sum of the pair interactions between water molecules and the solute pair. It should be noted that
(10)
where gcyl(h, dr) and Eww(h, dr) are the normalized density distribution and the local water-water interaction energy, respectively, shown in Figs. 3 and 6; and the integral with the volume element is over the whole volume. The solvation enthalpy contribution to w*(r) is then expressed as
(11)
We obtained ΔUws(r) from MD simulations for the contact-minimum configuration and the infinitely separated configuration, and then evaluated ΔUww(r) as Δh*(r) − ΔUws(r). The results for the three sized LJ solutes in water are shown in Fig. 8(b).

For the LJ particles in the contact configuration, ΔUws(r) is positive because the number of water molecules in the first solvation shells is smaller when the two solute particles are in contact than when they are far apart [Fig. 8(a)] and increases monotonically with the particle diameter σ* because the larger the particle, the greater the decrease in the number of water molecules in the solvation shells. This behavior of ΔUws(r) should be observed for any solvent species as long as there is an attractive force between solute and solvent molecules. On the other hand, the particle size dependence of ΔUww(r) reflects the hydrogen-bonding property of water. For the methane-sized solutes, ΔUww(r) is slightly positive because water molecules in the vicinity of solute particles can form stronger hydrogen bonds when the solutes are far apart than they are in contact. For the C60-sized solutes, however, ΔUww(r) is large and negative because water molecules form weaker hydrogen bonds in the first solvation shells than in bulk. The opposite solute size dependencies of ΔUww(r) and ΔUws(r) give rise to the non-monotonic solute size dependence of Δh*(r), as shown in Fig. 5(a). For the fullerene-sized LJ solutes in the contact configuration, Δh*(r) is close to zero because the positive ΔUws(r) and the negative ΔUww(r) largely cancel each other out.

In summary, for the fullerene-sized LJ solutes (σ* = 3), the water-induced potential w*(r) is essentially the entropy change due to the excluded volume effect. The enthalpy contribution is close to zero for the reason described above. On the other hand, for methane-sized LJ solutes (σ* = 1) in water near room temperature, the magnitude of the favorable entropy and that of the unfavorable enthalpy are nearly balanced, so that w*(r) is slightly negative.

Up until now, we examined the solute-size effect on w(r), Δh*(r), −TΔs*(r), Eww, and NHB under the condition that the LJ energy parameter for the solute–water interaction is fixed to the value for methane–water pairs. We shall now examine the effect of the weak solute–water attractive interaction on the solvent-induced potential w*(r). If the weak attraction is completely turned off, the potential of mean force curve changes significantly. Figure 9 shows a comparison of the PMF curve obtained for solutes interacting with water molecules via the repulsive part of the Weeks–Chandler–Andersen (WCA) potential (“water-repellent” solutes) with the other for solutes interacting via the LJ potential, both types of solutes having the same diameter of σ* = 3. First of all, the PMF at the solute–solute contact distance is significantly lower for the water-repellent solutes: the potential-well depth exceeds 40 kJ/mol. Second, there exists no potential barrier between the contact distance and the monolayer- or bilayer-separated distance for the WCA solutes. The reason why the C60-sized WCA solutes attract each other much more strongly than the LJ solutes is that both ΔUws and ΔUww are negative and so Δh*(r) is large and negative; Δs*(r) is large and positive as it is for the LJ solutes. In summary, the large pairwise attractive interaction is driven by both the favorable entropy change and the favorable enthalpy change.

FIG. 9.

Potentials w(r) of mean force for pairs of solutes with σ* = 3 in water at T = 300 K, reported in the previous paper.28 The blue and red curves are w(r) for the solutes interacting with water molecules via the LJ potential and the repulsive WCA potential, respectively. The black dotted lines emphasize the distances r of the first minimum and maximum, the second minimum and maximum, and the third minimum of w(r) for the LJ particle with σ* = 3.

FIG. 9.

Potentials w(r) of mean force for pairs of solutes with σ* = 3 in water at T = 300 K, reported in the previous paper.28 The blue and red curves are w(r) for the solutes interacting with water molecules via the LJ potential and the repulsive WCA potential, respectively. The black dotted lines emphasize the distances r of the first minimum and maximum, the second minimum and maximum, and the third minimum of w(r) for the LJ particle with σ* = 3.

Close modal

It should be noted that for the all-atom model for C60 molecules in water, the water-induced potential w*(r) at the contact-minimum distance is positive17,85,86 in contrast to w*(r) being negative for the C60 sized LJ and WCA particles. This is due to the strong water–C60 attractive interaction.86 

Figure 10 shows gcyl(h, d) and NHB(h, d) for the pair of WCA solutes with σ* = 3. When r is of the third minimum of w(r), or equivalently the bilayer-separated distance, gcyl in the immediate neighborhood of the solutes is ∼1.4, which is slightly lower than that for the LJ solutes; however, NHB is not significantly affected by the absence of the water–solute attractive interaction [compare Fig. 10(a) with Fig. 7(a)]. When r is of the second maximum or of the second minimum of w(r), water molecules between the WCA solutes have less hydrogen bonds than those between the LJ solutes. See Figs. 10(b), 10(c), 7(b), and 7(c). At the same time, the local density between the WCA solutes is smaller than that between the LJ solutes. When r is of the first maximum, water molecules in between the WCA particles are completely depleted [Fig. 10(d)] while those between LJ solutes are scarce but present [Fig. 3(d)].

FIG. 10.

Microscopic structures of water molecules near the two solute particles interacting with water molecules via the repulsive WCA potential at T = 300 K. The solute LJ diameter is σ* = 3. Left panel: the normalized two-dimensional distributions gcyl(h, d) of water molecules. Right panel: the numbers NHB of hydrogen bonds per one molecule. The solute–solute separations R are r of (a) the third minimum, (b) the second maximum, (c) the second minimum, (d) the first maximum, and (e) the first minimum of w(r) for the LJ particle with σ* = 3 (see Table I).

FIG. 10.

Microscopic structures of water molecules near the two solute particles interacting with water molecules via the repulsive WCA potential at T = 300 K. The solute LJ diameter is σ* = 3. Left panel: the normalized two-dimensional distributions gcyl(h, d) of water molecules. Right panel: the numbers NHB of hydrogen bonds per one molecule. The solute–solute separations R are r of (a) the third minimum, (b) the second maximum, (c) the second minimum, (d) the first maximum, and (e) the first minimum of w(r) for the LJ particle with σ* = 3 (see Table I).

Close modal

Figure 9 shows that the free-energy barrier between the contact and solvent-separated configurations of two LJ particles in water disappears when the water–solute weak attractive interaction is completely turned off. The microscopic mechanism is now clearly shown by Fig. 10(d): water molecules are so scarce or absent in the gap between the WCA particles that the two particles can approach each other without squeezing out water.

We have observed that the nature of the hydrophobic interaction changes as the diameter of the solute increases from that of methane to that of C60. For the methane-sized solutes in the contact minimum configuration, the favorable entropy change Δs*(r)/k (>0) and the unfavorable enthalpy change Δh*(r)/kT(>0) largely cancel each other out to make the water-induced potential w*(r)/kT slightly negative. On the other hand, for the C60-sized solutes, the water-induced potential w*(r)/kT is large and negative, entirely due to the favorable entropy change [Fig. 5(a)]. To focus on the solute-size effect, the LJ energy parameters for the solute–solute and solute–water pair potentials have been fixed to those for methane–methane and methane–water pairs. We demonstrated that w*(r) is very sensitive to the strength of the solute–solvent attractive interaction (Fig. 9).

For small hydrophobic solutes, such as methane, the potential w(r) of mean force has the well-defined minima corresponding to the contact, monolayer-separated, bilayer-separated configurations. As the size of solute particles increases, the first minimum decreases monotonically and the PMF curve in the range between the contact minimum and the bilayer-separated minimum goes up and eventually bears no local extrema in that range. The consequence is the single free-energy barrier between the contact and bilayer-separated configurations.

The decrement of the contact minimum in w(r) with increasing solute diameter is accompanied by the increasing favorable solvation entropy term TΔs*(r), which exceeds the increasing unfavorable solvation enthalpy term Δh*(r). The increase in the excess solvation entropy of the contact pair over that of the infinitely separated pair has been confirmed for both aqueous solutions and simple liquids based on simulation and the theory of liquids.11,13,15,20–24,27,28 The essential feature of the solute size dependence of TΔs*(r) has already been captured at the level of the Asakura–Oosawa excluded-volume effect28,87 and should, therefore, be considered as a common nature of solvent-induced pair interactions. Whether the entropic term is the major factor in determining the solute-size dependence of the contact minimum in w(r), however, depends on details of the solute–solvent and the solvent–solvent interactions.24 

The development of the free-energy barrier between the contact pair and the bilayer-separated pair with increasing size of the solute is due to the increasing unfavorable solvation enthalpy contribution. The free energy barrier disappears when the solute–water attractive interaction is replaced by a purely repulsive force (Fig. 9). We also know that in simple liquids, the solvent-separated minima do not disappear with increasing solute size.28 A deeper insight into the free energy barrier would be gained by evaluating the enthalpic and entropic contributions to w(r) for different sized solutes in simple liquids.

Two notable features of the solvation of pairs of hydrophobic solutes are derived (Figs. 3, 6, and 7). First, methane molecules hardly perturb the hydrogen bond network of water while C60-sized hydrophobic particles cause the surrounding water molecules to have fewer hydrogen bonds than in the bulk (Fig. 7). For the C60-sized solutes in the contact configuration, the water–water potential energy change ΔUww(r) is large and negative, the water–solute potential energy change ΔUws(r) is large and positive, and consequently, the enthalpy change Δh*(r) is close to zero. That is, the hydrophobic attraction between the C60-sized solutes is driven by the solvation entropy change TΔs*(r). Second, the hydrogen bond network between two large hydrophobic solutes is much less perturbed when the two solutes are in the bilayer-separated configuration. That is, the signature of the solvation force between two macroscopic hydrophobic surfaces81–84,86 is already anticipated in the case of C60-sized solutes.

Finally, the power-law behavior [Eq. (9)] of the osmotic second virial coefficient B has been found at 360 K, affirming the robustness of the previously reported results.27,28

This work was supported by KAKENHI (Grant Nos. JP18KK0151 and JP20H02696) and by JST, the establishment of university fellowships toward the creations of science technology innovation (Grant No. JPMJFS2128). Part of the computation was performed using the Research Center for Computational Science, Okazaki, Japan (Project Nos. 22-IMS-C124, 23-IMS-C112, and 24-IMS-C106).

The authors have no conflicts to disclose.

Hidefumi Naito: Data curation (equal); Formal analysis (equal); Investigation (lead); Methodology (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Tomonari Sumi: Formal analysis (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). Kenichiro Koga: Conceptualization (lead); Formal analysis (equal); Investigation (equal); Project administration (lead); Supervision (lead); Validation (equal); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
W.
Kauzmann
, “
Some factors in the interpretation of protein denaturation
,”
Adv. Protein Chem.
14
,
1
63
(
1959
).
2.
C.
Tanford
,
The Hydrophobic Effect: Formation of Micelles and Biological Membranes
, 2nd ed. (
Wiley
,
New York
,
1980
).
3.
A.
Ben-Naim
,
Hydrophobic Interactions
(
Plenum
,
Oxford
,
1980
).
4.
L. R.
Pratt
and
A.
Pohorille
,
Chem. Rev.
102
,
2671
(
2002
).
5.
N. T.
Southall
,
K. A.
Dill
, and
A.
Haymet
,
J. Phys. Chem. B
106
,
521
(
2002
).
6.
B.
Widom
,
P.
Bhimalapuram
, and
K.
Koga
,
Phys. Chem. Chem. Phys.
5
,
3085
(
2003
).
7.
K.
Koga
,
J. Phys. Chem. B
117
,
12619
(
2013
).
8.
K.
Koga
and
N.
Yamamoto
,
J. Phys. Chem. B
122
,
3655
(
2018
).
10.
L. R.
Pratt
and
D.
Chandler
,
J. Chem. Phys.
67
,
3683
(
1977
).
11.
K.
Lum
,
D.
Chandler
, and
J. D.
Weeks
,
J. Phys. Chem. B
103
,
4570
(
1999
).
12.
H.
Tanaka
,
J. Chem. Phys.
86
,
1512
(
1987
).
13.
G.
Hummer
,
S.
Garde
,
A. E.
Garcia
,
A.
Pohorille
, and
L. R.
Pratt
,
Proc. Natl. Acad. Sci. U. S. A.
93
,
8951
(
1996
).
14.
N. T.
Southall
and
K. A.
Dill
,
Biophys. Chem.
101
,
295
(
2002
).
15.
T.
Sumi
and
H.
Sekino
,
J. Chem. Phys.
126
,
144508
(
2007
).
16.
E.
Sobolewski
,
M.
Makowski
,
C.
Czaplewski
,
A.
Liwo
,
S.
Ołdziej
, and
H. A.
Scheraga
,
J. Phys. Chem. B
111
,
10765
(
2007
).
17.
M.
Makowski
,
C.
Czaplewski
,
A.
Liwo
, and
H. A.
Scheraga
,
J. Phys. Chem. B
114
,
993
(
2010
).
18.
A.
Bartosik
,
M.
Wiśniewska
, and
M.
Makowski
,
J. Phys. Org. Chem.
28
,
10
(
2015
).
19.
M.
Bogunia
and
M.
Makowski
,
J. Phys. Chem. B
124
,
10326
(
2020
).
20.
G.
Graziano
,
J. Phys. Chem. B
113
,
11232
(
2009
).
21.
22.
23.
R.
Zangi
,
J. Phys. Chem. B
115
,
2303
(
2011
).
24.
D.
Ben-Amotz
,
J. Phys. Chem. Lett.
6
,
1696
(
2015
).
26.
Q.
Sun
,
X.
Su
, and
C.
Cheng
,
Chem. Phys.
516
,
199
(
2019
).
27.
H.
Naito
,
R.
Okamoto
,
T.
Sumi
, and
K.
Koga
,
J. Chem. Phys.
156
,
221104
(
2022
).
28.
H.
Naito
,
T.
Sumi
, and
K.
Koga
,
Faraday Discuss.
249
,
440
(
2024
).
29.
Y.
Kimura
and
Y.
Yoshimura
,
Mol. Phys.
72
,
279
(
1991
).
30.
P.
Attard
,
J. Chem. Phys.
91
,
3083
(
1989
).
31.
P.
Attard
and
G.
Patey
,
J. Chem. Phys.
92
,
4970
(
1990
).
32.
R.
Dickman
,
P.
Attard
, and
V.
Simonian
,
J. Chem. Phys.
107
,
205
(
1997
).
33.
T.
Biben
,
P.
Bladon
, and
D.
Frenkel
,
J. Phys.: Condens. Matter
8
,
10799
(
1996
).
34.
M.
Dijkstra
,
R.
van Roij
, and
R.
Evans
,
Phys. Rev. E
59
,
5744
(
1999
).
35.
R.
Roth
,
R.
Evans
, and
S.
Dietrich
,
Phys. Rev. E
62
,
5360
(
2000
).
36.
R.
Roth
and
M.
Kinoshita
,
J. Chem. Phys.
125
,
084910
(
2006
).
37.
38.
R.
Akiyama
,
Y.
Karino
,
Y.
Hagiwara
, and
M.
Kinoshita
,
J. Phys. Soc. Jpn.
75
,
064804
(
2006
).
39.
J.-P.
Simonin
,
J. Phys. Chem. B
105
,
5262
(
2001
).
40.
W. G.
McMillan
, Jr.
and
J. E.
Mayer
,
J. Chem. Phys.
13
,
276
(
1945
).
41.
D. E.
Smith
,
L.
Zhang
, and
A. D. J.
Haymet
,
J. Am. Chem. Soc.
114
,
5875
(
1992
).
42.
D. E.
Smith
and
A.
Haymet
,
J. Chem. Phys.
98
,
6445
(
1993
).
43.
P.
Jungwirth
and
R.
Zahradník
,
Chem. Phys. Lett.
217
,
319
(
1994
).
44.
S. W.
Rick
and
B.
Berne
,
J. Phys. Chem. B
101
,
10488
(
1997
).
45.
S. W.
Rick
,
J. Phys. Chem. B
107
,
9853
(
2003
).
46.
S.
Shimizu
and
H. S.
Chan
,
J. Chem. Phys.
113
,
4683
(
2000
).
47.
S.
Shimizu
and
H. S.
Chan
,
J. Am. Chem. Soc.
123
,
2083
(
2001
).
48.
D.
Paschek
,
J. Chem. Phys.
120
,
6674
(
2004
).
49.
D.
Paschek
,
J. Chem. Phys.
120
,
10605
(
2004
).
50.
G.
Graziano
,
J. Phys. Soc. Jpn.
85
,
024801
(
2016
).
51.
K.
Zieba
,
C.
Czaplewski
,
A.
Liwo
, and
G.
Graziano
,
Phys. Chem. Chem. Phys.
22
,
4758
(
2020
).
52.
M.
Bogunia
,
A.
Liwo
,
C.
Czaplewski
,
J.
Makowska
,
A.
Giełdoń
, and
M.
Makowski
,
J. Phys. Chem. B
126
,
634
(
2022
).
53.
M. I.
Chaudhari
,
S. A.
Holleran
,
H. S.
Ashbaugh
, and
L. R.
Pratt
,
Proc. Natl. Acad. Sci. U. S. A.
110
,
20557
(
2013
).
54.
L. R.
Pratt
,
M. I.
Chaudhari
, and
S. B.
Rempe
,
J. Phys. Chem. B
120
,
6455
(
2016
).
55.
M. I.
Chaudhari
,
S. B.
Rempe
,
D.
Asthagiri
,
L.
Tan
, and
L.
Pratt
,
J. Phys. Chem. B
120
,
1864
(
2016
).
56.
H. S.
Ashbaugh
,
K.
Weiss
,
S. M.
Williams
,
B.
Meng
, and
L. N.
Surampudi
,
J. Phys. Chem. B
119
,
6280
(
2015
).
57.
D.
Tang
,
C.
Delpo
,
O.
Blackmon
, and
H. S.
Ashbaugh
,
J. Chem. Phys.
148
,
016101
(
2018
).
58.
C. A.
Cerdeiriña
and
B.
Widom
,
J. Phys. Chem. B
120
,
13144
(
2016
).
59.
S.
Lüdemann
,
R.
Abseher
,
H.
Schreiber
, and
O.
Steinhauser
,
J. Am. Chem. Soc.
119
,
4206
(
1997
).
60.
S. W.
Rick
,
J. Phys. Chem. B
104
,
6884
(
2000
).
61.
N.
Islam
,
M.
Flint
, and
S. W.
Rick
,
J. Chem. Phys.
150
,
014502
(
2019
).
62.
R.
Zangi
and
B.
Berne
,
J. Phys. Chem. B
112
,
8634
(
2008
).
63.
E.
Sobolewski
,
S.
Ołdziej
,
M.
Wisniewska
,
A.
Liwo
, and
M.
Makowski
,
J. Phys. Chem. B
116
,
6844
(
2012
).
64.
J. L. F.
Abascal
and
C.
Vega
,
J. Chem. Phys.
123
,
234505
(
2005
).
65.
M. G.
Martin
and
J. I.
Siepmann
,
J. Phys. Chem. B
102
,
2569
(
1998
).
66.
M. J.
Abraham
,
T.
Murtola
,
R.
Schulz
,
S.
Páll
,
J. C.
Smith
,
B.
Hess
, and
E.
Lindahl
,
SoftwareX
1–2
,
19
(
2015
).
67.
G. M.
Torrie
and
J. P.
Valleau
,
Chem. Phys. Lett.
28
,
578
(
1974
).
68.
G. M.
Torrie
and
J. P.
Valleau
,
J. Comput. Phys.
23
,
187
(
1977
).
69.
J. D.
Weeks
,
D.
Chandler
, and
H. C.
Andersen
,
J. Chem. Phys.
54
,
5237
(
1971
).
70.
S.
Kumar
,
J. M.
Rosenberg
,
D.
Bouzida
,
R. H.
Swendsen
, and
P. A.
Kollman
,
J. Comput. Chem.
13
,
1011
(
1992
).
71.
M.
Souaille
and
B.
Roux
,
Comput. Phys. Commun.
135
,
40
(
2001
).
72.
P.
Krüger
,
S. K.
Schnell
,
D.
Bedeaux
,
S.
Kjelstrup
,
T. J.
Vlugt
, and
J.-M.
Simon
,
J. Phys. Chem. Lett.
4
,
235
(
2013
).
73.
N.
Dawass
,
P.
Krüger
,
S. K.
Schnell
,
O. A.
Moultos
,
I. G.
Economou
,
T. J.
Vlugt
, and
J.-M.
Simon
,
Nanomaterials
10
,
771
(
2020
).
74.
K.
Koga
,
Phys. Chem. Chem. Phys.
13
,
19749
(
2011
).
75.
A.
Luzar
and
D.
Chandler
,
Phys. Rev. Lett.
76
,
928
(
1996
).
76.
A. K.
Soper
and
M. G.
Phillips
,
Chem. Phys.
107
,
47
(
1986
).
77.
J.
Teixeira
,
M.-C.
Bellissent-Funel
, and
S.-H.
Chen
,
J. Phys.: Condens. Matter
2
,
SA105
(
1990
).
78.
B.
Widom
and
R. C.
Underwood
,
J. Phys. Chem. B
116
,
9492
(
2012
).
79.
K.
Koga
,
V.
Holten
, and
B.
Widom
,
J. Phys. Chem. B
119
,
13391
(
2015
).
80.
F. H.
Stillinger
,
J. Solution Chem.
2
,
141
(
1973
).
81.
K.
Koga
,
J. Chem. Phys.
116
,
10882
(
2002
).
82.
K.
Koga
and
H.
Tanaka
,
J. Chem. Phys.
122
,
104711
(
2005
).
83.
J.
Engstler
and
N.
Giovambattista
,
J. Phys. Chem. B
122
,
8908
(
2018
).
84.
F.
Leoni
,
C.
Calero
, and
G.
Franzese
,
ACS Nano
15
,
19864
(
2021
).
85.
L.
Li
,
D.
Bedrov
, and
G. D.
Smith
,
Phys. Rev. E
71
,
011502
(
2005
).
86.
L.
Li
,
D.
Bedrov
, and
G. D.
Smith
,
J. Chem. Phys.
123
,
204504
(
2005
).
87.
F.
Oosawa
and
S.
Asakura
,
J. Chem. Phys.
22
,
1255
(
1954
).