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.

## I. INTRODUCTION

Hydrophobic effects and solvent-mediated phenomena, in general, are important in a wide variety of contexts^{1–4} ranging from protein folding^{5–7} and aggregation,^{8–10} to colloidal assembly^{11–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 simulations^{2,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, Δ*G*_{HS}, is intimately tied to fluctuations in water density through^{32,37,38}

where *P _{v}*(

*N*) is the probability of observing

*N*waters in the volume

*v*corresponding to the size and shape of the hard solute and

*β*= 1/

*k*

_{B}

*T*(

*k*

_{B}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 nm

^{3}) are Gaussian, allowing them to relate Δ

*G*

_{HS}to the moments of

*P*(

_{v}*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 proteins

^{39–41}and the formation of clathrate hydrates

^{41}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) predicted^{34,36,42} that low-*N* fluctuations should be enhanced (relative to Gaussian statistics) in large volumes (*v* ≳ 1 nm^{3}). Indeed, simulations subsequently verified that while *P _{v}*(

*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,

*P*(

_{v}*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, Δ*G*_{vdW}, is then given by

where Δ*G*_{att} 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 Δ*G*_{att} is accurately described by linear response theory (LRT), so that an understanding of $ P v ( att ) ( N ) $ informs Δ*G*_{vdW} in much the same way as *P _{v}*(

*N*) has informed Δ

*G*

_{HS}.

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 Δ*G*_{vdW} is concerned, this relative difficulty in emptying *v* in the presence of attractions is largely offset by the favorable free energy, Δ*G*_{att}, of turning on those attractions in the first place.

To minimize cancellation between the favorable *β*Δ*G*_{att} and the unfavorable $\u2212ln 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 Δ*G*_{vdW} 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 Δ*G*_{vdW}. 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 nm^{3}) solutes.^{51,56,57} Δ*G*_{vdW} 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 Δ*G*_{att} 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 $\u2212ln P v ( att ) ( N = 0 ) $, and therefore Δ*G*_{vdW}, 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.

## II. METHODS

### A. Models

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 material^{59}). 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 *R*_{C} and an attractive region of width *R*_{S} − *R*_{C} and depth *ϵ*. We consider three core sizes spanning the small to large solute size regimes of hydrophobic solvation,^{34} *R*_{C} = 0.336 nm, *R*_{C} = 0.6 nm, and *R*_{C} = 0.9 nm. The smallest solute corresponds to the effective hard sphere radius^{23,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, *R*_{S} = *R*_{C} + 0.3 nm, and set the well-depth, *βϵ* = 1.

While the repulsive part of the square-well potential is unambiguously given by the hard sphere pair potential, *u*_{0}, the attractive component can be defined in one of two ways (Figure 1). The first solute-water attractive potential, *u*_{1}(*r*), arises from a WCA-like separation of a standard square-well potential^{55} 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

and depends only on the number of molecules in the core and shell regions, $ N v ( R \xaf ) $ and $ N v sh ( R \xaf ) $ respectively, where $ R \xaf $ 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, *u*_{2}(*r*), can be considered (Figure 1(c)), and its contribution to the system Hamiltonian is

To vary the strength of attractions, we employ a linear coupling parameter in both cases, that is, we apply attractive potentials, λ_{i}*U _{i}*, 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, *v*_{s} = 2.5 × 2.5 × 1.0 nm^{3}, and has attractive interactions of strength λ_{1}*ϵ* with water in an adjoining volume, *v*_{1} = 2.5 × 2.5 × 0.3 nm^{3}. The cuboid-shaped hard solute, *v* = 1.5 × 1.5 × 0.3 nm^{3}, is then hydrated at the center of *v*_{1}.

### B. Simulation details and methods

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 barostat^{64} was used to maintain a pressure of 1 bar. Electrostatic interactions were handled with the particle mesh Ewald method^{65} 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, λ_{i}*U _{i}* (

*i*= 1, 2), by Δ

*G*, 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,

_{i}*U*

_{1}and

*U*

_{2}, depend only on

*N*and

_{v}*N*

_{vsh}, the number of water molecules in the solute core and hydration shell, respectively; that is, $ U i ( R \xaf ) = U i N v ( R \xaf ) , N v sh ( R \xaf ) $. Thus, both Δ

*G*and $ P v ( i ) ( N ) $ can be readily obtained for a range of attractive strengths, λ

_{i}_{i}, if the probability, $ P v , v sh N , N sh $, of observing

*N*and

*N*

_{sh}waters in

*v*and

*v*

_{sh}, respectively, in the absence of any solute-water attractions is known. In particular, $ P v ( i ) ( N ) $ can be obtained from the exact expression

where *C* is a normalization constant. Similarly, Δ*G _{i}* can be determined exactly from

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

The INDUS approach dictates coarse-graining the discrete variables *N _{v}* and

*N*

_{vsh}to the continuous variables $ N \u0303 v $ and $ N \u0303 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,

*κ*, and the biased potential minima, $ N i \u0303 \u2217 $, 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).

_{i}^{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 techniques^{69} 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 \u0303 s ( R \xaf ) = \varphi s N \u0303 v s ( R \xaf ) $, where $ N \u0303 v s ( R \xaf ) $ is the coarse-grained number of waters in *v*_{s}, 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, *v*_{1}, is $ U \u0303 1 ( R \xaf ) = \varphi 1 N \u0303 v 1 ( R \xaf ) $, where $ N \u0303 v 1 $ is the coarse-grained number of water molecules in *v*_{1}, and *ϕ*_{1} ≡ − λ_{1}*ϵ*. Both $ N \u0303 v s $ and $ N \u0303 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, $ ( \cdots ) $, can be related to averages in the mimic square-well system, $ ( \cdots ) \u2032 $ (simulated using the potentials $ U \u0303 s $ and $ U \u0303 1 $), through reweighting as

where *δ*_{i,j} is the Kronecker delta function, and *N*_{vs} and *N*_{v1} are the number of waters in *v*_{s} and *v*_{1}, 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

where *C*′ is a normalization constant.

## III. THE INFLUENCE OF ATTRACTIONS ON WATER DENSITY FLUCTUATIONS

The attractive field, λ_{1}*U*_{1}, which acts on both the core and shell regions, significantly affects water density fluctuations in the core, as shown in Figure 2(a) for *R*_{C} = 0.6 nm. As expected, the average number of molecules in the core region, $ N v ( R \xaf ) \lambda 1 $, increases with increasing field strength λ_{1}. Here, 〈(⋯)〉_{λi} represents an ensemble average in the presence of the attractive field *λ _{i}U_{i}*. 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 λ

_{1}

*U*

_{1}, we extend the model for density fluctuations proposed by Huang and Chandler

^{42}(see Appendix). The Huang-Chandler model assumes that reducing

*N*below its average value results in the formation of a bubble, and that the low-

_{v}*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 λ

_{1}

*U*

_{1}, 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 material

^{59}).

In contrast, in the presence of an attractive field such as λ_{2}*U*_{2}, 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 \xaf ) \lambda 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 *R*_{C}-values studied, the shifted spectra, $ P v ( 2 ) ( N ) = P v ( 2 ) ( N \u2212 \Delta N ) $, where $\Delta N\u2261 N v ( R \xaf ) \lambda 2 \u2212 N v ( R \xaf ) 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 *R*_{C}) or have pronounced fat tails (large *R*_{C}). This is the central result of this work.

To understand the basis of this invariance, we recognize that the addition of a linear potential (such as λ_{2}*N*_{vsh}) to a perfectly harmonic basin (such as one arising from *P*_{v,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-*N*_{sh} 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 statistics^{26,27} and $ P v , v sh N , N sh $ similarly follows a bivariate normal distribution with nonzero correlation between *N* and *N*_{sh} near its mean. Thus, the λ_{2}*U*_{2} potential, which is linear in $ N v sh ( R \xaf ) $ and couples to the Gaussian high-*N*_{sh} 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 *N*_{sh}-values. As λ_{2} is increased, the correlations between $ N v sh ( R \xaf ) $ and $ N v ( R \xaf ) $ then alter $ N v ( R \xaf ) $, but leave the remainder of $ P v ( 2 ) ( N ) $ unaltered.

While the application of *U*_{2} 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 *U*_{2} can be accurately estimated using linear response theory

given the correlation, $ \delta N v ( R \xaf ) \delta N v sh ( R \xaf ) 0 $, between fluctuations in $ N v ( R \xaf ) $ and $ N v sh ( R \xaf ) $ 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 λ_{1}*U*_{1} and λ_{2}*U*_{2}, we turn to the free energies of turning on these potentials in bulk water.

## IV. ATTRACTIONS IN THE SOLUTE CORE DO NOT INFLUENCE ITS HYDRATION FREE ENERGY

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 λ_{1}*U*_{1} is first turned on in the core and the shell, whereas in the lower path, the potential, λ_{2}*U*_{2}, is turned on in the shell alone. We denote the free energy for turning on the attractive field, λ_{i}*U _{i}*, in bulk water by Δ

*G*(

_{i}*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 $\u2212 k B Tln P v ( i ) ( N = 0 ) $ and is informed by the statistics of density fluctuations, $ P v ( i ) ( N ) $, shown in Figure 2.

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 Δ*G*_{1}, whereas the corresponding free energy, Δ*G*_{2}, 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 Δ*G _{i}* and unfavorable $\u2212 k B Tln P v ( i ) ( N = 0 ) $ terms. The diminished fluctuations in the presence of

*U*

_{1}(Figure 2(a)) mean that the unfavorable free energy to form a cavity is also larger in the presence of

*U*

_{1}. The total hydration free energy, Δ

*G*

_{vdW}, estimated through the upper path thus results in a significant cancellation between the favorable Δ

*G*

_{1}and the unfavorable $\u2212 k B Tln 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 Δ

*G*

_{vdW}.

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 *β*Δ*G _{i}* and unfavorable $\u2212ln 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 Δ

*G*, as discussed in Sec. V.

_{i}## V. EFFICIENT ESTIMATION OF BULK HYDRATION FREE ENERGIES

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):

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

where *v _{i}* is the volume over which the potential is applied, and

*ρ*

_{B}and

*κ*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

_{T}*U*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.

_{i}^{59}

LR3 suggests that Δ*G _{i}* 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 Δ

*G*is proportional to the volume,

_{i}*v*, over which the attractive interactions act and is quadratic in the strength of the attractions, λ

_{i}_{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 material

^{59}). The discrepancy between the predicted and simulation results for typical attraction strengths

^{72}is less than five percent. The approximations, LR1 and LR2, provide even more accurate estimates of Δ

*G*, with errors of less than a percent.

_{i}In addition to being able to obtain accurate estimates for Δ*G _{i}* 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, $\u2212ln 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

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 Δ*G*_{vdW}, as shown in Figure 7 for the *R*_{C} = 0.9 nm solute. This approach thus presents an efficient way to accurately estimate bulk hydration free energies of realistic attractive solutes, Δ*G*_{vdW}, starting from information on the statistics of water density fluctuations in bulk water, *P _{v}*(

*N*).

## VI. EFFICIENT ESTIMATION OF HYDRATION FREE ENERGIES IN HETEROGENEOUS ENVIRONMENTS

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)).

Indeed, the free energy needed to create a cavity in the presence of the attractive field is significantly higher (roughly 30 *k*_{B}*T* versus 10 *k*_{B}*T*), 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 *k*_{B}*T*^{74,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 *k*_{B}*T* in a first step, followed by cavity creation in the presence of reduced attractions.

## VII. CONCLUSIONS

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-Amotz^{56} 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 coworkers^{35,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 Chandler^{31} 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 LCW^{34} 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-workers^{28,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.

## Acknowledgments

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).

### APPENDIX: EXTENSION OF THE HUANG-CHANDLER MODEL

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 \xaf ) $ below its average value results in the formation of a spherical bubble of radius *r _{b}*, 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 $\Delta P=P\u2212 P sat $, where $P$ and $ P sat $ are the system pressure and the saturation pressure, respectively. Thus,

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 *r _{b}* is related to

*N*by

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

where we have defined

as an effective pressure in the probe volume, to illustrate that an attractive potential that couples to $ N v ( R \xaf ) $ affects *P _{v}*(

*N*) in the same manner as the pressure. Higher $\Delta P \u0303 $ 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 $\Delta P \u0303 $ and make it easier to create a cavity in the liquid.

## REFERENCES

*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.

*R*

_{C}= 0.3 nm and

*R*

_{C}= 0.9 nm, as well as comparisons to Lennard-Jones-like potentials are also described.

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 *R*_{C} = 0.6 nm solute, with a well-depth, λ_{i}*βϵ* ≈ 0.71, has the same second virial co-efficient as the highly attractive C_{60}-water interaction potential.