An understanding of density fluctuations in bulk water has made significant contributions to our understanding of the hydration and interactions of idealized, purely repulsive hydrophobic solutes. To similarly inform the hydration of realistic hydrophobic solutes that have dispersive interactions with water, here we characterize water density fluctuations in the presence of attractive fields that correspond to solute-water attractions. We find that when the attractive field acts only in the solute hydration shell, but not in the solute core, it does not significantly alter water density fluctuations in the solute core region. We further find that for a wide range of solute sizes and attraction strengths, the free energetics of turning on the attractive fields in bulk water are accurately captured by linear response theory. Our results also suggest strategies for more efficiently estimating hydration free energies of realistic solutes in bulk water and at interfaces.

Hydrophobic effects and solvent-mediated phenomena, in general, are important in a wide variety of contexts1–4 ranging from protein folding5–7 and aggregation,8–10 to colloidal assembly11–13 and detergency.14,15 When a hydrophobic solute is solvated by water, it excludes water molecules from the region that it occupies, thereby perturbing the inherent interactions between the nearby water molecules.1–4,16–18 Because water molecules have strong hydrogen bond interactions, it is the penalty for this disruption of the water structure that dominates hydrophobic hydration free energies. Hard sphere solutes that simply exclude water molecules from the region they occupy have thus served as idealized hydrophobic solutes, and their hydration and assembly have been extensively studied using both molecular simulations2,19–29 and theory.3,17,20,30–36

The use of hard solutes as ideal hydrophobes is particularly judicious because their excess hydration free energy, ΔGHS, is intimately tied to fluctuations in water density through32,37,38

β Δ G HS = ln P v ( N = 0 ) ,
(1)

where Pv(N) is the probability of observing N waters in the volume v corresponding to the size and shape of the hard solute and β = 1/kBT (kB is Boltzmann’s constant and T is the temperature). Thus, an understanding of density fluctuations in bulk water has the potential to inform free energies of hydrophobic hydration and association. Indeed, the notion that small fluctuations in water density obey Gaussian statistics lies at the heart of the Pratt-Chandler theory of hydrophobic hydration.30,31 Using molecular simulations, Hummer et al. explicitly showed that density fluctuations in small observation volumes (v ≲ 1 nm3) are Gaussian, allowing them to relate ΔGHS to the moments of Pv(N).32,33 This simplification has afforded molecular insights into various phenomena, including entropy convergence, wherein protein unfolding entropies converge at a common temperature,28,33 as well as the denaturation of proteins39–41 and the formation of clathrate hydrates41 at high pressures.

Noting that water at ambient conditions is in close proximity to liquid-vapor coexistence, the theory of Lum, Chandler, and Weeks (LCW) predicted34,36,42 that low-N fluctuations should be enhanced (relative to Gaussian statistics) in large volumes (v ≳ 1 nm3). Indeed, simulations subsequently verified that while Pv(N) remains Gaussian near its mean, the likelihood of low-N fluctuations in large v is enhanced significantly in bulk water and even more so near hydrophobic surfaces; that is, Pv(N) develops fat low-N tails.4,26,27,42–48 This perspective clarified that water near hydrophobic surfaces sits at the edge of a dewetting transition that can be readily triggered by small perturbations.48–53 It also led to the prediction that the assembly of small hydrophobic solutes in the vicinity of extended hydrophobic surfaces would be barrierless.47,54 Thus, an understanding of density fluctuations in small and large volumes as well as in the vicinity of interfaces has made significant contributions to our understanding of hydrophobic hydration and assembly.

In addition to excluding water, realistic solutes also possess favorable attractive (van der Waals) interactions with water. Here, our goal is to establish a connection between water density fluctuations and hydration free energies of realistic solutes of various sizes in bulk water and at interfaces. To accomplish this, we first turn on solute-water attractive interactions, and then characterize water density fluctuations corresponding to the emptying of the solute repulsive core in the presence of these attractions. The hydration free energy of realistic van der Waals solutes, ΔGvdW, is then given by

β Δ G vdW = β Δ G att ln P v ( att ) N = 0 ,
(2)

where ΔGatt is the free energy for turning on the solute-water attractions in bulk water, and P v ( att ) ( N ) is now the probability of observing N waters in v in the presence of the attractive field. We show that ΔGatt is accurately described by linear response theory (LRT), so that an understanding of P v ( att ) ( N ) informs ΔGvdW in much the same way as Pv(N) has informed ΔGHS.

We quantify P v ( att ) ( N ) for spherical volumes of different sizes and a range of attractive strengths. We find that the presence of attractions in both the solute core, v, and its hydration shell, following the Weeks-Chandler-Andersen (WCA) prescription,55 significantly alters water density fluctuations; it becomes progressively harder to empty v as the strength of attractions is increased. However, as far as the estimation of ΔGvdW is concerned, this relative difficulty in emptying v in the presence of attractions is largely offset by the favorable free energy, ΔGatt, of turning on those attractions in the first place.

To minimize cancellation between the favorable βΔGatt and the unfavorable ln P v ( att ) ( N = 0 ) terms in Eq. (2), we consider an alternative prescription for hydrating the same solute; one that involves turning on attractions in the hydration shell, but not in the solute core, v. The overall value of ΔGvdW does not depend on which prescription is used and, in particular, whether attractions are turned on in the solute core or not; attractions in the core simply increase the magnitude of the components of ΔGvdW. In the presence of attractions in the hydration shell alone, we find that the water density fluctuations in the core are remarkably unaltered; attractions effect only a subtle change in the mean density. We find this to be true for solutes of all sizes and reasonable attraction strengths. Thus, our primary finding is that water density fluctuations that are relevant to hydrophobic hydration are largely unaffected by the presence of attractions.

Our results also suggest strategies for circumventing the breakdown of perturbation theories of hydrophobic hydration, which occurs for large (≳ 1 nm3) solutes.51,56,57 ΔGvdW is typically estimated by first creating a cavity and subsequently turning on attractions. The free energy for turning on attractions can be readily estimated if accurate approximations are available for the response of the hydration shell water density to the attractions. This approach works well for small solutes because water structure around the cavity is not significantly altered when attractions are turned on, so that water density responds linearly to attractions. In contrast, a soft vapor-liquid interface is nucleated around large cavities and is readily displaced towards the solute when attractions are turned on. As a result, water density near the cavity increases in a sigmoidal fashion (and not linearly) as the strength of attractions is increased, thereby violating a key assumption that underpins perturbation theory. Because we turn on attractions in bulk water prior to the formation of a cavity, water density responds linearly to the strength of attractions, and ΔGatt is readily estimated using linear response theory. Water density fluctuations are largely unaltered when attractions are turned on in the hydration shell alone, so that ln P v ( att ) ( N = 0 ) , and therefore ΔGvdW, can be estimated from a knowledge of water density fluctuations in the absence of attractions.

Furthermore, recognizing that the presence of attractions in the core of a solute has no bearing on its hydration free energy, also allows for efficient estimation of hydration free energies in interfacial environments. When hydrating a solute in such inhomogeneous environments, attractions from neighboring solutes or interfaces acting on v make it harder to empty the solute core. We recommend that such attractions should be turned off (at least partially) in a first step, with the corresponding free energy readily estimated using linear response theory. The solute core can then be emptied more easily in a second step because the free energy for doing so is substantially reduced.

To obtain a qualitative understanding of the influence of attractive potentials on density fluctuations, we focus on the solvation of a square well particle in the extended simple point charge (SPC/E) model of water.58 While we focus on square-well solutes, our findings are general and hold for more realistic Lennard-Jones (LJ) solutes (see supplementary material59). We chose the SPC/E model for water because it reasonably mimics the bulk density, isothermal compressibility, surface tension, and liquid-vapor coexistence properties of real water; these properties primarily influence the behavior of water density fluctuations.18 

Bulk Hydration: The solutes we study interact with water oxygens through the pair potential shown in Figure 1(a). The potential has a hard core of radius RC and an attractive region of width RSRC and depth ϵ. We consider three core sizes spanning the small to large solute size regimes of hydrophobic solvation,34RC = 0.336 nm, RC = 0.6 nm, and RC = 0.9 nm. The smallest solute corresponds to the effective hard sphere radius23,60 of a united atom methane with Lennard-Jones parameters, σLJ = 0.3448 nm and ϵLJ = 0.8956 kJ/mol.61 For all solutes, we choose a spherical shell region 0.3 nm in width, that is, RS = RC + 0.3 nm, and set the well-depth, βϵ = 1.

FIG. 1.

(a) The square well potential, u(r), with a hard core repulsion at r = RC and a well of depth, ϵ, and width, RSRC. (b) This solute-water potential can readily be split into its repulsive, u0(r), and attractive, u1(r), components using a WCA-like prescription.55 (c) Alternatively, u(r) can be split into u0(r), and the attractive potential, u2(r), which is non-zero only in the shell region, RC < r < RS.

FIG. 1.

(a) The square well potential, u(r), with a hard core repulsion at r = RC and a well of depth, ϵ, and width, RSRC. (b) This solute-water potential can readily be split into its repulsive, u0(r), and attractive, u1(r), components using a WCA-like prescription.55 (c) Alternatively, u(r) can be split into u0(r), and the attractive potential, u2(r), which is non-zero only in the shell region, RC < r < RS.

Close modal

While the repulsive part of the square-well potential is unambiguously given by the hard sphere pair potential, u0, the attractive component can be defined in one of two ways (Figure 1). The first solute-water attractive potential, u1(r), arises from a WCA-like separation of a standard square-well potential55 shown in Figure 1(b) and acts on both the core and shell regions. The contribution of such solute-water attractions to the Hamiltonian of the system is

U 1 ( R ¯ ) = ϵ ( N v ( R ¯ ) + N v sh ( R ¯ ) )
(3)

and depends only on the number of molecules in the core and shell regions, N v ( R ¯ ) and N v sh ( R ¯ ) respectively, where R ¯ denotes a configuration vector containing the positions of all the water oxygens. The results obtained from this potential can be qualitatively compared to the WCA attractive portion of the LJ potential; such a comparison is shown in the supplementary material.59 Alternatively, an attractive potential that acts on the shell alone, u2(r), can be considered (Figure 1(c)), and its contribution to the system Hamiltonian is

U 2 ( R ¯ ) = ϵ N v sh ( R ¯ ) .
(4)

To vary the strength of attractions, we employ a linear coupling parameter in both cases, that is, we apply attractive potentials, λiUi, and vary λi to change the interaction strength.

Hydration at interfaces: We also consider the hydration of a hard solute in the vicinity of a square-well attractive surface. The surface excludes waters from a cuboid-shaped volume, vs = 2.5 × 2.5 × 1.0 nm3, and has attractive interactions of strength λ1ϵ with water in an adjoining volume, v1 = 2.5 × 2.5 × 0.3 nm3. The cuboid-shaped hard solute, v = 1.5 × 1.5 × 0.3 nm3, is then hydrated at the center of v1.

All simulation data presented here were obtained using version 4.5.3 of the GROMACS molecular dynamics (MD) simulation package,62 modified to include the various external and biasing potentials used in this work. MD simulations were performed in the isothermal-isobaric (NPT) ensemble, where the canonical velocity-rescaling thermostat of Bussi et al.63 was used to maintain a constant temperature of 300 K, and a Parinello-Rahman barostat64 was used to maintain a pressure of 1 bar. Electrostatic interactions were handled with the particle mesh Ewald method65 with a real space cutoff of 1 nm and grid spacing of 0.12 nm. Similarly, van der Waals interactions were cutoff at a distance of 1 nm, and standard energy and pressure corrections were used to account for the influence of the truncated interactions.66 

Bulk hydration: We denote the free energy of turning on the attractive field, λiUi (i = 1, 2), by ΔGi, and the fluctuations in the presence of the field by P v ( i ) ( N ) . For the specific case of square-well type potentials, the fields, U1 and U2, depend only on Nv and Nvsh, the number of water molecules in the solute core and hydration shell, respectively; that is, U i ( R ¯ ) = U i N v ( R ¯ ) , N v sh ( R ¯ ) . Thus, both ΔGi and P v ( i ) ( N ) can be readily obtained for a range of attractive strengths, λi, if the probability, P v , v sh N , N sh , of observing N and Nsh waters in v and vsh, respectively, in the absence of any solute-water attractions is known. In particular, P v ( i ) ( N ) can be obtained from the exact expression

P v ( i ) ( N ) = C N sh = 0 P v , v sh N , N sh e β λ i U i ( N , N sh ) ,
(5)

where C is a normalization constant. Similarly, ΔGi can be determined exactly from

β Δ G i = ln N N sh P v , v sh N , N sh e β λ i U i ( N , N sh ) .
(6)

Derivations of both Eqs. (5) and (6) are included in the supplementary material.59 To obtain P v , v sh N , N sh , we employed indirect umbrella sampling (INDUS)27 with harmonic biasing potentials of the form

U bias ( R ¯ ) = κ 1 2 ( N ̃ v ( R ¯ ) N ̃ ) 2 + κ 2 2 ( N ̃ v sh ( R ¯ ) N ̃ sh ) 2 .
(7)

The INDUS approach dictates coarse-graining the discrete variables Nv and Nvsh to the continuous variables N ̃ v and N ̃ v sh , respectively, through convolution with a truncated and shifted Gaussian; here, we use a width of 0.01 nm and a cutoff of 0.03 nm for this smoothing function. The force constants, κi, and the biased potential minima, N i ̃ , were tuned to achieve sufficient overlap between neighboring windows. The corresponding P v , v sh N , N sh distributions were then obtained from the biased simulations using the unbinned weighted-histogram analysis method (UWHAM).67,68

Hydration at interfaces: Because simulating the square-well surface using MD would result in impulsive forces, we instead employ closely related continuous potentials to mimic the square-well potential and make use of standard reweighting techniques69 to estimate the behavior of a true square-well surface. The exclusion of water molecules from the square-well surface is accomplished by applying the potential U ̃ s ( R ¯ ) = ϕ s N ̃ v s ( R ¯ ) , where N ̃ v s ( R ¯ ) is the coarse-grained number of waters in vs, and βϕs = 12 was chosen to yield a substantial number of configurations with the volume empty. Analogously, the attractive potential in the adjoining hydration shell volume, v1, is U ̃ 1 ( R ¯ ) = ϕ 1 N ̃ v 1 ( R ¯ ) , where N ̃ v 1 is the coarse-grained number of water molecules in v1, and ϕ1 ≡ − λ1ϵ. Both N ̃ v s and N ̃ v 1 were defined according to the INDUS prescription described above. To estimate the probability, P v ( 1 ) ( N ) , of observing N water molecules in v, a series of INDUS simulations were then performed using a biasing potential similar to the one in Eq. (7), but with κ2 = 0. Averages in a system with a true square-well surface, ( ) , can be related to averages in the mimic square-well system, ( ) (simulated using the potentials U ̃ s and U ̃ 1 ), through reweighting as

( ) = ( ) δ 0 , N v s ( R ¯ ) e β ϕ s N ̃ v s ( R ¯ ) β ϕ 1 [ N v 1 ( R ¯ ) N ̃ v 1 ( R ¯ ) ] δ 0 , N v s ( R ¯ ) e β ϕ s N ̃ v s ( R ¯ ) β ϕ 1 [ N v 1 ( R ¯ ) N ̃ v 1 ( R ¯ ) ] ,
(8)

where δi,j is the Kronecker delta function, and Nvs and Nv1 are the number of waters in vs and v1, respectively. This average, in turn, is evaluated using UWHAM.67,68 The probability, P v ( 2 ) ( N ) , of finding N waters in v, in the absence of favorable surface-water interactions in v, is obtained through another reweighting

P v ( 2 ) ( N ) = C P v ( 1 ) ( N ) e β λ 1 ϵ N ,
(9)

where C′ is a normalization constant.

The attractive field, λ1U1, which acts on both the core and shell regions, significantly affects water density fluctuations in the core, as shown in Figure 2(a) for RC = 0.6 nm. As expected, the average number of molecules in the core region, N v ( R ¯ ) λ 1 , increases with increasing field strength λ1. Here, 〈(⋯)〉λi represents an ensemble average in the presence of the attractive field λiUi. Importantly, the non-Gaussian nature, or the fatness of the low-N tails is also diminished, making it increasingly difficult to empty the core region as λ1 is increased. To better understand how water density fluctuations are affected by λ1U1, we extend the model for density fluctuations proposed by Huang and Chandler42 (see Appendix). The Huang-Chandler model assumes that reducing Nv below its average value results in the formation of a bubble, and that the low-N tail in the water density fluctuations can be described by the energetics of growing the bubble. As shown in the Appendix, in the presence of λ1U1, additional work has to be done against the attractive field to expand the bubble, which results in the low-N tails of the P v ( 1 ) ( N ) distributions being diminished. This qualitative behavior is also observed for both smaller and larger solutes (see supplementary material59).

FIG. 2.

(a) The presence of an attractive field that acts on both the core and shell regions, λ1U1, drastically alters density fluctuations in the core region, as shown here for a spherical solute with RC = 0.6 nm and RS = 0.9 nm. On increasing the attractive strength λ1, the tails at low N become steeper, indicating that it is more difficult to evacuate the core volume. (b) In contrast, fluctuations in the core region are not significantly perturbed by the attractive field, λ2U2, which acts on the hydration shell alone. As the strength of the attractive interactions, λ2, is increased, only a small decrease in the average number of molecules in the core is observed.

FIG. 2.

(a) The presence of an attractive field that acts on both the core and shell regions, λ1U1, drastically alters density fluctuations in the core region, as shown here for a spherical solute with RC = 0.6 nm and RS = 0.9 nm. On increasing the attractive strength λ1, the tails at low N become steeper, indicating that it is more difficult to evacuate the core volume. (b) In contrast, fluctuations in the core region are not significantly perturbed by the attractive field, λ2U2, which acts on the hydration shell alone. As the strength of the attractive interactions, λ2, is increased, only a small decrease in the average number of molecules in the core is observed.

Close modal

In contrast, in the presence of an attractive field such as λ2U2, which acts only in the hydration shell, the nature of the density fluctuations in the core region is not changed; see Figure 2(b). Only a small decrease in the mean, N v ( R ¯ ) λ 2 , is observed when the attractive strength, λ2, is increased. In Figure 3, the fluctuation spectra, P v ( 2 ) , are shifted horizontally so that their means are aligned. For all the RC-values studied, the shifted spectra, P v ( 2 ) ( N ) = P v ( 2 ) ( N Δ N ) , where Δ N N v ( R ¯ ) λ 2 N v ( R ¯ ) 0 , corresponding to various attractions, λ2, collapse onto a universal curve. Remarkably, this invariance of water density fluctuations around the mean is independent of whether the fluctuations are nearly Gaussian (small RC) or have pronounced fat tails (large RC). This is the central result of this work.

FIG. 3.

Water number distributions, P v ( 2 ) , in the presence of attractive fields λ2U2, shifted horizontally to align their means, N v ( R ¯ ) λ 2 , are shown for probe volumes with radii (a) equivalent to that of a methane, (b) 0.6 nm, and (c) 0.9 nm. Fluctuation spectra obtained over a range of attractions collapse onto a single curve when plotted as a function of N N Δ N , where Δ N N v ( R ¯ ) λ 2 N v ( R ¯ ) 0 , highlighting that the presence of an attractive potential in the hydration shell simply decreases the average number of water molecules in the core and not the fluctuations around the mean. This is true independent of solute size and whether the fluctuations follow Gaussian statistics (solid gray lines).

FIG. 3.

Water number distributions, P v ( 2 ) , in the presence of attractive fields λ2U2, shifted horizontally to align their means, N v ( R ¯ ) λ 2 , are shown for probe volumes with radii (a) equivalent to that of a methane, (b) 0.6 nm, and (c) 0.9 nm. Fluctuation spectra obtained over a range of attractions collapse onto a single curve when plotted as a function of N N Δ N , where Δ N N v ( R ¯ ) λ 2 N v ( R ¯ ) 0 , highlighting that the presence of an attractive potential in the hydration shell simply decreases the average number of water molecules in the core and not the fluctuations around the mean. This is true independent of solute size and whether the fluctuations follow Gaussian statistics (solid gray lines).

Close modal

To understand the basis of this invariance, we recognize that the addition of a linear potential (such as λ2Nvsh) to a perfectly harmonic basin (such as one arising from Pv,vsh being Gaussian) simply translates the basin (and hence the entire distribution) horizontally.48 Attractions in the shell favor configurations with more waters than the average, coupling to the high-Nsh region of P v , v sh N , N sh . Such high-N regions of water number distributions in bulk water have previously been shown to follow Gaussian statistics26,27 and P v , v sh N , N sh similarly follows a bivariate normal distribution with nonzero correlation between N and Nsh near its mean. Thus, the λ2U2 potential, which is linear in N v sh ( R ¯ ) and couples to the Gaussian high-Nsh region of P v , v sh N , N sh is expected to shift the entire P v , v sh N , N sh distribution towards higher Nsh-values. As λ2 is increased, the correlations between N v sh ( R ¯ ) and N v ( R ¯ ) then alter N v ( R ¯ ) , but leave the remainder of P v ( 2 ) ( N ) unaltered.

While the application of U2 increases water density in the shell region, packing effects can lead to a concomitant decrease in the average number of waters in the core region. The layering of water density at the core-shell interface, shown in Figure 4(a), highlights that attractions in the shell indeed decrease the average density in the core. Interestingly, as shown in Figure 4(b), this decrease in the average number of molecules in the core region upon application of the potential U2 can be accurately estimated using linear response theory

Δ N λ 2 β ϵ δ N v ( R ¯ ) δ N v sh ( R ¯ ) 0
(10)

given the correlation, δ N v ( R ¯ ) δ N v sh ( R ¯ ) 0 , between fluctuations in N v ( R ¯ ) and N v sh ( R ¯ ) in the absence of attractions. As discussed in Sec. V, this agreement will facilitate the development of approximations for estimating hydration free energies of realistic solutes in bulk water. For now, having quantified the ease with which the core can be emptied in the presence of the attractive fields λ1U1 and λ2U2, we turn to the free energies of turning on these potentials in bulk water.

FIG. 4.

(a) Attractions in the shell region, indicated by vertical gray lines, increase water density in the shell, but result in a decrease in the average water density in the core region, as shown here for RC = 0.6 nm. The decrease in the core density on increasing λ2 is a direct result of layering at the interface between the core and shell regions. (b) This decrease in the average number of water molecules in the core region, ΔN, as the strength of the potential λ2 is increased, can be accurately predicted with LRT, Eq. (10).

FIG. 4.

(a) Attractions in the shell region, indicated by vertical gray lines, increase water density in the shell, but result in a decrease in the average water density in the core region, as shown here for RC = 0.6 nm. The decrease in the core density on increasing λ2 is a direct result of layering at the interface between the core and shell regions. (b) This decrease in the average number of water molecules in the core region, ΔN, as the strength of the potential λ2 is increased, can be accurately predicted with LRT, Eq. (10).

Close modal

Two alternative paths for hydrating the square-well solute of Figure 1(a) are shown in Figure 5(a). In the upper path, the attractive potential λ1U1 is first turned on in the core and the shell, whereas in the lower path, the potential, λ2U2, is turned on in the shell alone. We denote the free energy for turning on the attractive field, λiUi, in bulk water by ΔGi (i = 1, 2). Then, in the second step, a cavity is created in the core region in the presence of the corresponding attractive potential. The free energy for creating such a cavity is given by k B T ln P v ( i ) ( N = 0 ) and is informed by the statistics of density fluctuations, P v ( i ) ( N ) , shown in Figure 2.

FIG. 5.

(a) The hydration of a realistic hydrophobic solute is carried out in two steps: an attractive field, λiUi, corresponding to the solute-water attractions is first turned on in bulk water, and water is emptied from the solute core in a subsequent step. The corresponding free energies are denoted by βΔGi and ln P v ( i ) ( 0 ) , respectively, where i = 1 represents attractions in both the core and shell regions (top) and i = 2 represents attractions in only the shell region (bottom). ((b) and (c)) The free energy of turning on the attractive potential U1 is much larger than that for turning on U2. (d) This results in a larger cancellation between the favorable βΔGi and the unfavorable ln P v ( i ) ( 0 ) free energies for i = 1, because the total hydration free energy Δ G vdW = Δ G i ln P v ( i ) ( 0 ) is the same for both paths, as shown here for RC = 0.9 nm and λiβϵ = 0.5.

FIG. 5.

(a) The hydration of a realistic hydrophobic solute is carried out in two steps: an attractive field, λiUi, corresponding to the solute-water attractions is first turned on in bulk water, and water is emptied from the solute core in a subsequent step. The corresponding free energies are denoted by βΔGi and ln P v ( i ) ( 0 ) , respectively, where i = 1 represents attractions in both the core and shell regions (top) and i = 2 represents attractions in only the shell region (bottom). ((b) and (c)) The free energy of turning on the attractive potential U1 is much larger than that for turning on U2. (d) This results in a larger cancellation between the favorable βΔGi and the unfavorable ln P v ( i ) ( 0 ) free energies for i = 1, because the total hydration free energy Δ G vdW = Δ G i ln P v ( i ) ( 0 ) is the same for both paths, as shown here for RC = 0.9 nm and λiβϵ = 0.5.

Close modal

As shown in Figures 5(b) and 5(c), turning on attractions in the both the core and shell regions results in a large and negative free energetic gain ΔG1, whereas the corresponding free energy, ΔG2, for turning on attractions in the shell alone is significantly smaller in magnitude. However, the total hydration free energy is given by the sum of the favorable ΔGi and unfavorable k B T ln P v ( i ) ( N = 0 ) terms. The diminished fluctuations in the presence of U1 (Figure 2(a)) mean that the unfavorable free energy to form a cavity is also larger in the presence of U1. The total hydration free energy, ΔGvdW, estimated through the upper path thus results in a significant cancellation between the favorable ΔG1 and the unfavorable k B T ln P v ( 1 ) ( N = 0 ) terms, as shown in Figure 5(d). Such large cancellations can lead to significant numerical uncertainty in the estimation of ΔGvdW.

In contrast, the lower path in Figure 5(a) reduces the extent of such a cancellation by recognizing that the presence or absence of attractions in the solute core plays no role in determining the overall solute hydration free energy and minimizing the volume over which the attractions act. In principle, the cancellation between the favorable βΔGi and unfavorable ln P v ( i ) ( N = 0 ) terms can be further reduced or even eliminated by turning on repulsions in the core in conjunction with attractions in the shell in the first step. Even more complex, low-variance pathways to hydration may also be chosen to cleverly minimize the computational overhead.70,71 However, the advantage of the paths considered here is that they allow for particularly simple analytical approximations for ΔGi, as discussed in Sec. V.

Because attractive fields couple to the Gaussian high-N tails of water density distributions, we expect linear response theory to provide an accurate estimate of the free energy for turning on such fields. In particular, we consider the following approximate forms of Eq. (6):

β Δ G i β λ i 2 U i ( R ¯ ) 0 + U i ( R ¯ ) λ i
(11)
β λ i U i ( R ¯ ) 0 β 2 λ i 2 2 ( δ U i ( R ¯ ) ) 2 0
(12)

referred to as LR1 and LR2, respectively. In the context of thermodynamic integration, LR1 represents the application of a trapezoid rule to integrate over the entire range of λi, while LR2 uses the initial slope of the thermodynamic force at λi = 0 and extrapolates over the range of λi to integrate. For the square-well potentials considered here, the free energy of turning on the attractive external field can further be approximated analytically, yielding

β Δ G i v i ρ B ( λ i β ϵ ) ρ B 2 κ T k B T 2 ( λ i β ϵ ) 2 ,
(13)

where vi is the volume over which the potential is applied, and ρB and κT are the density and isothermal compressibility of bulk water, respectively. To derive Eq. (13) (referred to as LR3) from LR2, we recognize that fluctuations in Ui are related to fluctuations in water density, which in turn are related to the isothermal compressibility. The exact derivations of LR3 and its generic form applicable to slowly varying potentials are provided in the supplementary material.59 

LR3 suggests that ΔGi can be estimated simply from a knowledge of the strength of the attractions, the volume over which they act, and bulk properties of the solvent. It predicts that ΔGi is proportional to the volume, vi, over which the attractive interactions act and is quadratic in the strength of the attractions, λi, consistent with previous findings.73 We find LR3 to be true to a very good approximation, as illustrated in Figure 6; it begins to break down only for small volumes and large λi (see supplementary material59). The discrepancy between the predicted and simulation results for typical attraction strengths72 is less than five percent. The approximations, LR1 and LR2, provide even more accurate estimates of ΔGi, with errors of less than a percent.

FIG. 6.

(a) The free energy, ΔGi, of turning on the attractive potential, normalized by the volume, vi, on which the attractions act, is largely independent of solute size and the range over which the potential acts. ΔGi is well described by the linear response prediction LR3 of Eq. (13) (solid black line). (b) The three linear response theory estimates for ΔGi are found to be highly accurate, as shown here for the RC = 0.9 nm solute. The dashed line indicates a typical magnitude of attractions for realistic solutes.72 

FIG. 6.

(a) The free energy, ΔGi, of turning on the attractive potential, normalized by the volume, vi, on which the attractions act, is largely independent of solute size and the range over which the potential acts. ΔGi is well described by the linear response prediction LR3 of Eq. (13) (solid black line). (b) The three linear response theory estimates for ΔGi are found to be highly accurate, as shown here for the RC = 0.9 nm solute. The dashed line indicates a typical magnitude of attractions for realistic solutes.72 

Close modal

In addition to being able to obtain accurate estimates for ΔGi from linear response theory, our essential finding that attractions (in the hydration shell of a solute) do not affect water density fluctuations, also permits accurate estimates of the free energy for emptying the solute core, ln P v ( 2 ) ( N = 0 ) . Turning on attractive interactions in the hydration shell only results in a small shift in the mean, ΔN (Figure 3), which itself can be readily estimated using linear response theory (Eq. (10)), so that fluctuations in the presence and absence of attractions can be related through

P v ( 2 ) ( N ) P v ( N + Δ N ) .
(14)

Evaluating the components of the total bulk hydration free energy (Eq. (2)) using Eqs. (13) and (14) (for N = 0) leads to accurate predictions for ΔGvdW, as shown in Figure 7 for the RC = 0.9 nm solute. This approach thus presents an efficient way to accurately estimate bulk hydration free energies of realistic attractive solutes, ΔGvdW, starting from information on the statistics of water density fluctuations in bulk water, Pv(N).

FIG. 7.

(a) The total hydration free energy of a realistic hydrophobic solute, Δ G vdW = Δ G 2 k B T ln P v ( 2 ) ( N = 0 ) , can be described with good accuracy using the approximations in Eqs. (13) and (14), as shown here for the solute with RC = 0.9 nm. (b) Indeed, for realistic values72 of λ2 (dashed line), the error in ΔGvdW is less than 10%.

FIG. 7.

(a) The total hydration free energy of a realistic hydrophobic solute, Δ G vdW = Δ G 2 k B T ln P v ( 2 ) ( N = 0 ) , can be described with good accuracy using the approximations in Eqs. (13) and (14), as shown here for the solute with RC = 0.9 nm. (b) Indeed, for realistic values72 of λ2 (dashed line), the error in ΔGvdW is less than 10%.

Close modal

Figure 5(d) highlights that it is efficient to turn off any attractive fields in the solute core before estimating the free energy for emptying the core. Such attractive fields (acting on water) can originate not only from the solute itself but also from neighboring solutes or interfaces. To demonstrate this, in Figure 8(a), we illustrate the hydration of a hard solute (shown in green) in the vicinity of a surface (shown in purple) that interacts with water through an attractive square well potential (dark blue gradient). The attractive interactions from the square well act on the core region of the solute, making it harder to empty (upper path). Alternatively, the attractive potential can first be turned off in the solute core region, as shown in the lower path of Figure 8(a), with the corresponding free energy accurately estimated by linear response theory. Turning off the attractions also ensures that emptying the solute cavity is easier, making it the more efficient alternative (Figure 8(c)).

FIG. 8.

(a) Top: The hydration of a 1.5 × 1.5 × 0.3 nm3 cuboid-shaped hard solute (green) adjacent to a 2.5 × 2.5 × 1.0 nm3 attractive square-well surface (purple, well depth βϵ = 1 and width 0.3 nm), entails the creation of a cavity in the presence of attractive surface-water interactions. Bottom: Alternatively, the solute can be hydrated in two steps; turning off the attractions in the core of the solute in the first step, followed by creating a cavity in the second step. (b) Emptying the solute core requires a significantly smaller free energy in the absence of attractions in the core. (c) While the solute hydration free energy is the same in both cases, the free energy for turning off the attractions, −βΔG2, can be accurately estimated by using linear response theory, and the free energy for subsequently emptying the solute core, ln P v ( 2 ) ( 0 ) , is smaller than ln P v ( 1 ) ( 0 ) , and therefore easier to estimate.

FIG. 8.

(a) Top: The hydration of a 1.5 × 1.5 × 0.3 nm3 cuboid-shaped hard solute (green) adjacent to a 2.5 × 2.5 × 1.0 nm3 attractive square-well surface (purple, well depth βϵ = 1 and width 0.3 nm), entails the creation of a cavity in the presence of attractive surface-water interactions. Bottom: Alternatively, the solute can be hydrated in two steps; turning off the attractions in the core of the solute in the first step, followed by creating a cavity in the second step. (b) Emptying the solute core requires a significantly smaller free energy in the absence of attractions in the core. (c) While the solute hydration free energy is the same in both cases, the free energy for turning off the attractions, −βΔG2, can be accurately estimated by using linear response theory, and the free energy for subsequently emptying the solute core, ln P v ( 2 ) ( 0 ) , is smaller than ln P v ( 1 ) ( 0 ) , and therefore easier to estimate.

Close modal

Indeed, the free energy needed to create a cavity in the presence of the attractive field is significantly higher (roughly 30 kBT versus 10 kBT), as seen in Figure 8(b). When the attractive field is turned off in the solute core, the average number of water molecules in the core decreases and a significant non-Gaussian tail emerges in P v ( 2 ) ( N ) at low N; both effects act in concert to lower the free energy needed to form a cavity. Such a lowering of the cavity formation free energy allows for more efficient estimation of hydration free energies in heterogeneous environments, because fewer biased simulations are needed to probe the entire range of density fluctuations in the solute core.

Caveat: While the free energy of turning off the attractive interactions inside the solute core was accurately described using linear response theory (LR1) here, this may not always be the case. In particular, turning off attractions completely in a sufficiently large volume adjacent to an extended hydrophobic surface can dewet the volume, causing linear response to break down.56 However, only weak attractions are needed to wet a hydrophobic surface; recent work suggests that an attractive strength of roughly 0.1 kBT74,75 is sufficient to prevent dewetting. Thus, attractive may be turned off partially without triggering dewetting, and the corresponding free energy safely estimated using linear response. In such cases, our revised recommendation is to reduce the surface-water attractions down to 0.1 kBT in a first step, followed by cavity creation in the presence of reduced attractions.

Water number fluctuations in small volumes follow Gaussian statistics, while those in large volumes have non-Gaussian tails at low density. An understanding of these fluctuations has made seminal contributions to our understanding of the hydration and interactions of idealized, purely repulsive hydrophobic solutes. To similarly inform the hydration of realistic hydrophobic solutes that have attractive interactions with water, here we have provided a description of these fluctuations in the presence of attractive fields. We find that WCA-like attractive potentials,55 which are non-zero both inside and outside the solute core, alter the nature of water density fluctuations significantly, making it more difficult to create a cavity. In contrast, when attractions are turned on in the hydration shell alone, water density fluctuations are largely unaltered in both small and large volumes.

We also find that the favorable free energy for turning on the various attractive fields is readily approximated by linear response theory and is proportional to the volume on which the fields act. Consequently, the free energy of turning on attractions in both the core and the shell is much larger than the free energy of turning on attractions in the hydration shell alone. When attractions are turned on in both the core and the shell, the favorable free energy of turning on those attractions and the unfavorable free energetics of emptying the solute core in the presence of those attractions largely cancel in the estimation of the total solute hydration free energy. This cancellation is minimized when the solute-solvent attractions are turned on in the hydration shell alone, making it the more efficient route for estimating hydration free energies of realistic solutes.

In addition to informing hydrophobic hydration and interactions, an understanding of density fluctuations has also facilitated the development of novel simulation methods. By leveraging an understanding of water density fluctuations, and the associated response of water density to perturbations, Patel and Garde recently introduced a method for estimating cavity hydration free energies that is two orders of magnitude more efficient as compared with umbrella sampling.29 Our characterization of density fluctuations in the presence of slowly varying attractive interactions has similarly allowed us to suggest strategies for more efficiently estimating hydration free energies of realistic solutes, both in bulk water and in the vicinity of interfaces.

Our approach involves turning on attractions in bulk water, followed by emptying the solute core, and is in contrast with traditional methods that first create a cavity in bulk water, and subsequently turn on attractions. Because attractions have little effect on the structure of water adjacent to small cavities (small enough for the hydrogen bond network of water to be maintained around them),23,28,73,76 traditional methods can readily estimate the free energy for turning on attractions using linear response theory. However, larger cavities induce dewetting in their vicinity; attractions rewet the solute, so that water density and consequently the solute-water interaction energy do not vary linearly with the strength of attractions.1 Indeed, recent work from Underwood and Ben-Amotz56 has shown a transition from linear to non-linear response as the solute size is increased,56 occurring at roughly 1 nm; the length scale corresponding to the crossover in hydrophobic hydration.34 Because our approach involves turning on attractions in bulk water before creating a solute core, it circumvents the breakdown of linear response theory that plagues the traditional perturbative treatment of attractive interactions.

The approach to solvation presented here is similar in spirit to the two-step method considered by Weeks and coworkers35,77 in their development of the molecular-scale van der Waals (MVDW) theory for inhomogeneous systems, which was subsequently combined with the Gaussian field model of Chandler31 to yield the LCW theory of hydrophobicity.34 The MVDW theory first considers turning on all slowly varying interactions in the bulk fluid, which include not only the solute-solvent attractions that we consider here but also long-ranged solvent-solvent interactions. Near large hydrophobic solutes, the long-ranged solvent-solvent attractive interactions are unbalanced and result in a net repulsion away from the solute.35 Such repulsions couple to the non-Gaussian low N tails in water number distributions and result in dewetting; thus their effect cannot be captured using linear response. Recognizing this, the LCW34 theory employed a Landau-Ginzburg free energy functional that enables repulsive potentials to induce dewetting. In contrast, because our focus here has been on slowly varying attractive potentials that couple to the Gaussian high-N region of water number distributions, we are able to make judicious use of linear response theory.

The results presented here rely on two essential properties of attractive interactions: (i) they are slowly varying and thereby minimally perturb the structure of water, and (ii) they couple to Gaussian tails of water number distributions. Thus, the lessons learned from this work will also be applicable to other interactions that possess these two characteristic features. In particular, recent work by Weeks and co-workers28,78–82 has shown that it is instructive to resolve electrostatic interactions into short- and long-ranged components, in analogy with the WCA prescription, which divides the Lennard-Jones potential into a short-ranged repulsive and a slowly varying attractive component.55 While the long-ranged electrostatic component follows linear response theory, underpinning dielectric continuum theories,83 the short-ranged interactions are more complex, leading to directional hydrogen bonds and specific ion effects.84,85 Thus, in hydrating ions, for example, it would be instructive to turn on the long-ranged component of electrostatic interactions first, and investigate how it affects water density fluctuations; such calculations will be the focus of future work.

The authors acknowledge partial financial support from the National Science Foundation through a Seed Grant from the University of Pennsylvania Materials Research Science and Engineering Center (NSF UPENN MRSEC DMR 11-20901).

Here, we describe the response of P v ( 1 ) ( N ) to an attractive potential by extending the model developed by Huang and Chandler (HC).42 The model assumes that reducing N v ( R ¯ ) below its average value results in the formation of a spherical bubble of radius rb, and approximates the density outside the bubble to be the bulk density, ρB. The free energetics of density fluctuations (in the absence of an external field) are then given by the work that must be performed to grow the bubble against the surface tension γ and the external pressure Δ P = P P sat , where P and P sat are the system pressure and the saturation pressure, respectively. Thus,

k B T ln P v ( N ) = G 0 r b ( N ) 4 π 3 r b 3 ( N ) Δ P + 4 π r b 2 ( N ) γ .
(A1)

For the sake of simplicity, the term corresponding to the translational entropy of the bubble has been omitted from the above expression. The radius of the bubble rb is related to N by

r b ( N ) = 3 4 π ρ B N v ( R ¯ ) 0 N 1 / 3 .
(A2)

In the presence of attractions in v, work must also be done against the attractive potential, such that

ln P v ( 1 ) ( N ) = ln P v ( N ) + 4 π 3 r b 3 ( N ) ρ B λ 1 β ϵ
(A3)
= 4 π 3 r b 3 ( N ) β Δ P ̃ + 4 π r b 2 ( N ) β γ ,
(A4)

where we have defined

Δ P ̃ Δ P + λ 1 ϵ ρ B
(A5)

as an effective pressure in the probe volume, to illustrate that an attractive potential that couples to N v ( R ¯ ) affects Pv(N) in the same manner as the pressure. Higher Δ P ̃ results in an increase in magnitude of the slope at low N, consistent with the results in Figure 2(a). Conversely, repulsive potentials, corresponding to negative values of λ1, should decrease Δ P ̃ and make it easier to create a cavity in the liquid.

3.
B. J.
Berne
,
J. D.
Weeks
, and
R.
Zhou
,
Annu. Rev. Phys. Chem.
60
,
85
(
2009
).
4.
S. N.
Jamadagni
,
R.
Godawat
, and
S.
Garde
,
Annu. Rev. Chem. Biomol. Eng.
2
,
147
(
2011
).
5.
N. T.
Southall
,
K. A.
Dill
, and
A. D. J.
Haymet
,
J. Phys. Chem. B
106
,
521
(
2002
).
6.
Y.
Levy
and
J. N.
Onuchic
,
Annu. Rev. Biophys. Biomol. Struct.
35
,
389
(
2006
).
7.
D.
Thirumalai
,
E. P.
O’Brien
,
G.
Morrison
, and
C.
Hyeon
,
Annu. Rev. Biophys.
39
,
159
(
2010
).
9.
M. G.
Krone
,
L.
Hua
,
P.
Soto
,
R.
Zhou
,
B. J.
Berne
, and
J.-E.
Shea
,
J. Am. Chem. Soc.
130
,
11066
(
2008
).
10.
D.
Thirumalai
,
G.
Reddy
, and
J. E.
Straub
,
Acc. Chem. Res.
45
,
83
(
2012
).
11.
G. M.
Whitesides
and
B.
Grzybowski
,
Science
295
,
2418
(
2002
).
12.
E.
Rabani
,
D. R.
Reichman
,
P. L.
Geissler
, and
L. E.
Brus
,
Nature
426
,
271
(
2003
).
13.
J. A.
Morrone
,
J.
Li
, and
B. J.
Berne
,
J. Phys. Chem. B
116
,
378
(
2012
).
14.
C.
Tanford
,
The Hydrophobic Effect: Formation of Micelles and Biological Membranes
(
John Wiley & Sons
,
1973
).
15.
L.
Maibaum
,
A. R.
Dinner
, and
D.
Chandler
,
J. Phys. Chem. B
108
,
6778
(
2004
).
16.
G.
Hummer
,
S.
Garde
,
A. E.
García
,
M. E.
Paulaitis
, and
L. R.
Pratt
,
J. Phys. Chem. B
102
,
10469
(
1998
).
17.
H. S.
Ashbaugh
and
L. R.
Pratt
,
Rev. Mod. Phys.
78
,
159
(
2006
).
18.
D.
Chandler
and
P.
Varilly
, in Proceeding of the International School of Physics “Enrico Fermi,” edited by F. Mallamace and H. E. Stanley (IOS, SIF, Amsterdam, Bologna, 2012), Vol. 176, pp. 75–111.
20.
F. H.
Stillinger
,
J. Solution Chem.
2
,
141
(
1973
).
21.
G. E.
Crooks
and
D.
Chandler
,
Phys. Rev. E
56
,
4217
(
1997
).
22.
D. M.
Huang
,
P. L.
Geissler
, and
D.
Chandler
,
J. Phys. Chem. B
105
,
6704
(
2001
).
23.
D. M.
Huang
and
D.
Chandler
,
J. Phys. Chem. B
106
,
2047
(
2002
).
24.
M. V.
Athawale
,
G.
Goel
,
T.
Ghosh
,
T. M.
Truskett
, and
S.
Garde
,
Proc. Natl. Acad. Sci. U. S. A.
104
,
733
(
2007
).
25.
G.
Goel
,
M. V.
Athawale
,
S.
Garde
, and
T. M.
Truskett
,
J. Phys. Chem. B
112
,
13193
(
2008
).
26.
A. J.
Patel
,
P.
Varilly
, and
D.
Chandler
,
J. Phys. Chem. B
114
,
1632
(
2010
).
27.
A. J.
Patel
,
P.
Varilly
,
D.
Chandler
, and
S.
Garde
,
J. Stat. Phys.
145
,
265
(
2011
).
28.
R. C.
Remsing
and
J. D.
Weeks
,
J. Phys. Chem. B
117
,
15479
(
2013
).
29.
A. J.
Patel
and
S.
Garde
,
J. Phys. Chem. B
118
,
1564
(
2014
).
30.
L. R.
Pratt
and
D.
Chandler
,
J. Chem. Phys.
67
,
3683
(
1977
).
31.
D.
Chandler
,
Phys. Rev. E
48
,
2898
(
1993
).
32.
G.
Hummer
,
S.
Garde
,
A. E.
Garcia
,
A.
Pohorille
, and
L. R.
Pratt
,
Proc. Natl. Acad. Sci. U. S. A.
93
,
8951
(
1996
).
33.
S.
Garde
,
G.
Hummer
,
A. E.
Garcia
,
M. E.
Paulaitis
, and
L. R.
Pratt
,
Phys. Rev. Lett.
77
,
4966
(
1996
).
34.
K.
Lum
,
D.
Chandler
, and
J. D.
Weeks
,
J. Phys. Chem. B
103
,
4570
(
1999
).
36.
P.
Varilly
,
A. J.
Patel
, and
D.
Chandler
,
J. Chem. Phys.
134
,
074109
(
2011
).
37.
B.
Widom
,
J. Chem. Phys.
39
,
2808
(
1963
).
38.
T. L.
Beck
,
M. E.
Paulaitis
, and
L. R.
Pratt
,
The Potential Distirbution Theorem and Models of Molecular Solutions
(
Cambridge University Press
,
2006
).
39.
P. R.
ten Wolde
,
J. Phys.: Condens. Matter
14
,
9445
(
2002
).
41.
G.
Hummer
,
S.
Garde
,
A.
Garcia
,
M.
Paulaitis
, and
L.
Pratt
,
Proc. Natl. Acad. Sci. U. S. A.
95
,
1552
(
1998
).
42.
D. M.
Huang
and
D.
Chandler
,
Phys. Rev. E
61
,
1501
(
2000
).
43.
J.
Mittal
and
G.
Hummer
,
Proc. Natl. Acad. Sci. U. S. A.
105
,
20130
(
2008
).
44.
S.
Sarupria
and
S.
Garde
,
Phys. Rev. Lett.
103
,
037803
(
2009
).
45.
R.
Godawat
,
S. N.
Jamadagni
, and
S.
Garde
,
Proc. Natl. Acad. Sci. U. S. A.
106
,
15119
(
2009
).
46.
H.
Acharya
,
S.
Vembanur
,
S. N.
Jamadagni
, and
S.
Garde
,
Faraday Discuss.
146
,
353
(
2010
).
47.
A. J.
Patel
,
P.
Varilly
,
S. N.
Jamadagni
,
H.
Acharya
,
S.
Garde
, and
D.
Chandler
,
Proc. Natl. Acad. Sci. U. S. A.
108
,
17678
(
2011
).
48.
A. J.
Patel
,
P.
Varilly
,
S. N.
Jamadagni
,
M. F.
Hagan
,
D.
Chandler
, and
S.
Garde
,
J. Phys. Chem. B
116
,
2498
(
2012
).
49.
R.
Zhou
,
X.
Huang
,
C. J.
Margulis
, and
B. J.
Berne
,
Science
305
,
1605
(
2004
).
50.
P.
Liu
,
X.
Huang
,
R.
Zhou
, and
B. J.
Berne
,
Nature
437
,
159
(
2005
).
51.
N.
Choudhury
and
B. M.
Pettitt
,
J. Am. Chem. Soc.
129
,
4847
(
2007
).
52.
N.
Choudhury
and
B. M.
Pettitt
,
Mol. Simul.
31
,
457
(
2005
).
53.
J. C.
Rasaiah
,
S.
Garde
, and
G.
Hummer
,
Annu. Rev. Phys. Chem.
59
,
713
(
2008
).
54.
S.
Vembanur
,
A. J.
Patel
,
S.
Sarupria
, and
S.
Garde
,
J. Phys. Chem. B
117
,
10261
(
2013
).
55.
J. D.
Weeks
,
D.
Chandler
, and
H. C.
Andersen
,
J. Chem. Phys.
54
,
5237
(
1971
).
56.
R.
Underwood
and
D.
Ben-Amotz
,
J. Chem. Phys.
135
,
201102
(
2011
).
57.
N.
Choudhury
and
B. M.
Pettitt
,
J. Am. Chem. Soc.
127
,
3556
(
2005
).
58.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
,
J. Phys. Chem.
91
,
6269
(
1987
).
59.
See supplementary material at http://dx.doi.org/10.1063/1.4905009 for a derivation of Eqs. (5), (6), (13), and a generalization of Eq. (13) that is applicable to a variety of slowly-varying attractive potentials. Detailed results for RC = 0.3 nm and RC = 0.9 nm, as well as comparisons to Lennard-Jones-like potentials are also described.
60.
H. C.
Andersen
,
J. D.
Weeks
, and
D.
Chandler
,
Phys. Rev. A
4
,
1597
(
1971
).
61.
W. D.
Cornell
,
P.
Cieplak
,
C. I.
Bayly
,
I. R.
Gould
,
K. M.
Merz
,
D. M.
Ferguson
,
D. C.
Spellmeyer
,
T.
Fox
,
J. W.
Caldwell
, and
P. A.
Kollman
,
J. Am. Chem. Soc.
117
,
5179
(
1995
).
62.
B.
Hess
,
C.
Kutzner
,
D.
van der Spoel
, and
E.
Lindahl
,
J. Chem. Theory Comput.
4
,
435
(
2008
).
63.
G.
Bussi
,
D.
Donadio
, and
M.
Parrinello
,
J. Chem. Phys.
126
,
014101
(
2007
).
64.
M.
Parrinello
and
A.
Rahman
,
J. Appl. Phys.
52
,
7182
(
1981
).
65.
U.
Essmann
,
L.
Perera
,
M. L.
Berkowitz
,
T.
Darden
,
H.
Lee
, and
L. G.
Pedersen
,
J. Chem. Phys.
103
,
8577
(
1995
).
66.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Oxford
,
New York
,
1987
).
67.
Z.
Tan
,
E.
Gallichio
,
M.
Lapelosa
, and
R. M.
Levy
,
J. Chem. Phys.
136
,
144102
(
2012
).
68.
M. R.
Shirts
and
J. D.
Chodera
,
J. Chem. Phys.
129
,
124105
(
2008
).
69.
G. M.
Torrie
and
J. P.
Valleau
,
J. Comput. Phys.
23
,
187
(
1977
).
70.
T. T.
Pham
and
M. R.
Shirts
,
J. Chem. Phys.
136
,
124120
(
2012
).
71.
L. N.
Naden
,
T. T.
Pham
, and
M. R.
Shirts
,
J. Chem. Theory Comput.
10
,
1128
(
2014
).
72.

We estimate an upper bound for the magnitude of solute-water attractions to be λiβϵ ≈ 0.7. By comparison, a square-well potential for the RC = 0.6 nm solute, with a well-depth, λiβϵ ≈ 0.71, has the same second virial co-efficient as the highly attractive C60-water interaction potential.

73.
L.
Maibaum
and
D.
Chandler
,
J. Phys. Chem. B
111
,
9025
(
2007
).
74.
R.
Godawat
,
S. N.
Jamadagni
,
V.
Venkateshwaran
, and
S.
Garde
, “
Connecting water correlations, fluctuations, and wetting phenomena at hydrophobic and hydrophilic surfaces
,” e-print arXiv:1409.2570 (
2014
).
75.
A. P.
Willard
and
D.
Chandler
,
J. Chem. Phys.
141
,
18C519
(
2014
).
76.
N.
Galamba
,
J. Phys. Chem. B
117
,
2153
(
2013
).
77.
J. D.
Weeks
,
K.
Katsov
, and
K.
Vollmayr
,
Phys. Rev. Lett.
81
,
4400
(
1998
).
78.
Y. G.
Chen
,
C.
Kaur
, and
J. D.
Weeks
,
J. Phys. Chem. B
108
,
19874
(
2004
).
79.
Y. G.
Chen
and
J. D.
Weeks
,
Proc. Natl. Acad. Sci. U. S. A.
103
,
7560
(
2006
).
80.
J. M.
Rodgers
and
J. D.
Weeks
,
J. Phys.: Condens. Matter
20
,
494206
(
2008
).
81.
J. M.
Rodgers
and
J. D.
Weeks
,
Proc. Natl. Acad. Sci. U. S. A.
105
,
19136
(
2008
).
82.
R. C.
Remsing
,
J. M.
Rodgers
, and
J. D.
Weeks
,
J. Stat. Phys.
145
,
313
(
2011
).
83.
R. C.
Remsing
, “
From structure to thermodynamics with local molecular field theory
,” Ph.D. thesis (
University of Maryland
,
2013
).
84.
Y.
Shi
and
T. L.
Beck
,
J. Chem. Phys.
139
,
044504
(
2013
).
85.

Supplementary Material