Entropy of Hydrophobic Hydration: Extension to Hydrophobic Chains

1900

J. Phys. Chem. 1996, 100, 1900-1913

Entropy of Hydrophobic Hydration: Extension to Hydrophobic Chains Henry S. Ashbaugh and Michael E. Paulaitis* Center for Molecular and Engineering Thermodynamics, Department of Chemical Engineering, UniVersity of Delaware, Newark, Delaware 19716 ReceiVed: August 16, 1995; In Final Form: October 25, 1995X

A statistical mechanical formulation for the entropy in terms of multiparticle correlation functions, used previously to calculate entropies of the hydration of simple hydrophobic solutes, has been generalized to molecular solutes of arbitrary shape by recasting the correlation function expansion as a summation over sites that define the solute molecule. The new formulation for the entropy is applied to a Monte Carlo simulation study of normal alkanes at infinite dilution in water to calculate contributions to the entropy of hydration from water-solute site pair correlations and to examine the role these contributions play in stabilizing different solute conformations. In this implemention, the water-solute site pair correlations are determined only for individual water molecules with their nearest solute site and are defined by water orientational and water oxygen radial distributions around the site, independent of solute orientation relative to the water molecule. We show that these distribution functions give an accurate representation of water structure around the individual n-alkane sites for methane through normal butane, account for the large negative entropies of hydration of these alkanes at 25 °C, and predict the stabilization of gauche-butane relative to trans-butane in water on the basis of an entropically favorable (energetically unfavorable) trans f gauche transition. Contributions to the entropy of hydration arising from solute-induced perturbations in water-water correlations (i.e., water structure enhancement) have also been examined, and we show energy-entropy compensation of these contributions within the framework of the correlation function expansion for the entropy.

Introduction The solubility of nonpolar compounds in aqueous solution at room temperature is dominated by a large negative entropy of solvation and an uncharacteristically large positive heat capacity increment. These properties, and a slightly favorable enthalpy of solvation at room temperature, are the distinguishing features of what is commonly called hydrophobic hydration.1-3 Much of the recent interest in hydrophobic hydration has been stimulated by the recognition that hydrophobic interactions play an important role in self-assembly processes in aqueous solutions. The hydrophobic interaction is viewed as a partial reversal of hydrophobic hydration and has been shown to be entropy driven at room temperature.4,5 Therefore, revealing the molecular origins of the large negative entropies of hydration is central to understanding the mechanisms underlying selfassembly phenomena such as micelle formation, molecular recognition, and protein folding. While a number of solution theories of hydrophobic hydration have met with some success in describing the thermodynamic properties of nonpolar compounds dissolved in water, many recent studies of hydration phenomena have turned to molecular simulations to gain insights into interactions between water and hydrophobic solutes that are not experimentally accessible. Recently, we presented a statistical mechanical formulation for the entropy of hydration of simple hydrophobic solutes based on an expansion6,7 for the entropy in terms of multiparticle correlation functions.8-11 Contributions to the entropy of hydration from water-solute pair correlations were evaluated by Monte Carlo simulations for methane, the inert gases, and model Lennard-Jones (LJ) particles at infinite dilution in water. It was found that solute-water translational and orientational pair correlations alone can account for the large negative * To whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, January 1, 1996.

0022-3654/96/20100-1900$12.00/0

entropies of hydration at room temperature and, for methane in water, the large positive heat capacity increment at room temperature.8 The primary factor determining entropies of hydration of the LJ solutes was found to be solute size, with the solute-water pair interaction energy playing a secondary, but still significant, role.10 Contributions to the hydration entropy from solute-induced perturbations in water-water correlations were not evaluated. However, for solutes the size of methane, these contributions appear to be small when the solute is added at constant temperature and pressure. In this work, we extend the correlation function expansion to aqueous solutions containing molecular solutes of arbitrary shape. The new formulation is applied to a simulation study of the hydration of normal alkanes in dilute aqueous solution to calculate water-solute pair correlation entropies of hydration and to examine the role these contributions play in stabilizing butane conformations. We also consider contributions arising from solute-induced perturbations in the water-water correlations and energy-entropy compensation of these contributions in the application of the correlation function expansion to calculate entropies of hydration of simple hydrophobic solutes. Correlation Function Expansion for the Entropy The entropy of a pure fluid can be expressed in terms of an ideal-gas entropy at the same temperature, volume, and number density and a residual or excess contribution due to intermolecular interactions. Thus,

S(T,V,N) ) SIG(T,V,N) + S*(T,V,N)

(1)

where SIG is the ideal-gas entropy and S* is the excess entropy. For a monatomic molecule,

SIG )

5NkB - NkB ln(FΛ3) 2

© 1996 American Chemical Society

(2)

Entropy of Hydrophobic Hydration

J. Phys. Chem., Vol. 100, No. 5, 1996 1901

where N is the number of molecules, kB is Boltzmann’s constant, T is temperature, F is the number density, and Λ is the thermal de Broglie wavelength. The excess entropy can be expanded in terms related to two-body, three-body, and higher order correlations. As shown by Baranyai and Evans,12 this expansion can be cast into the form

S * ) S (2) + S (3) + higher order terms

(3)

with

S (2) ) -

kBF2 ∫∫g (2)(r1,r2) ln[g (2)(r1,r2)] dr1 dr2 + 2 kBF2 ∫∫[g (2)(r1,r2) - 1] dr1 dr2 (4) 2

kBF3 ∫∫∫g (3)(r1,r2,r3) × 6 kBF3 ln[δg (3)(r1,r2,r3)] dr1 dr2 dr3 + ∫∫∫[g (3)(r1,r2,r3) 6 3g (2)(r1,r2) g (2)(r2,r3) + 3g (2)(r1,r2) - 1] dr1 dr2 dr3 (5)

S (3) ) -

and δg

(3)

5NSkB 5NAkB - NSkB ln(FSΛS3) + - NAkB ln(FAΛA3) 2 2 kBNSFS ∫gSS(r) ln[gSS(r)] dr - ∫[gSS(r) - 1] dr 2

S)

{ k N F { ∫g k N F ∫g 2 { B A S

B A A

δg (3)(r1,r2,r3) ≡

SA(r)

ln[gSA(r)] dr - ∫[gSA

AA(r)

ln[gAA(r)] dr - ∫[gAA

where Ni is the number of molecules of component i (i ) S for solvent; i ) A for solute), Fi is the number density of component i, and gSS, gSA, and gAA are the solvent-solvent, solvent-solute, and solute-solute pair correlation functions, respectively. The solute partial molar entropy, ShA, is obtained directly from eq 8 by differentiating with respect to NA at constant T, P, and NS. In the limit of infinite dilution,

ShA∞ )

3kB j A∞ - kB ln(FAΛA3) + kBFSυ 2

[ ∫g

k BF S

ln(gSA) dr - ∫(gSA - 1) dr -

[ ∫( )

SA

kBFS 2

∂gSS ∂xA

kBFS2υ j A∞

defined by

g (3)(r1,r2,r3) g (2)(r1,r2) g (2)(r1,r3) g (2)(r2,r3)

2 (6)

The functions g (2) and g (3) are the two- and three-body correlation functions, respectively, while S (2) and S (3) are the corresponding two- and three-body entropies. The higher order terms in eq 3 represent contributions to the entropy from correlations involving four or more molecules.13,14 The second term in the expressions for S (2) and S (3) arises in the grand canonical ensemble derivation for the entropy and makes the expansion local; i.e., the two terms in eq 4 and in eq 5 together converge more rapidly than does either term alone.12 It has also been shown that the second term in each expression accounts for the long-range behavior of the correlation functions in the canonical ensemble derivation of the expansion.10 For simple molecules, the pair correlation functions are easily evaluated from molecular simulations. However, triplet and higher order correlation functions are much more difficult to obtain because inordinately long simulation times are required or inordinately large system sizes must be considered to acquire reasonable statistical averages. Only recently have triplet correlations been derived from simulations of hard-sphere (HS) fluids,15 LJ fluids,16 and primitive-model electrolyte fluids.17 In addition, truncation of the correlation function expansion at the level of pair correlations; i.e.,

S ≈ S IG + S (2)

(7)

has been shown to provide accurate approximations for the entropy of liquid metals,18 noble gases,19 HS fluids,19,20 and LJ fluids12,16,20,21 at liquidlike densities. These results motivated our application of the correlation function expansion in previous studies of hydrophobic hydration, in which we calculated entropies of hydration on the basis of pair correlations alone. For the entropy of a binary mixture, the correlation function expansion truncated at the level of pair correlations is given by

} (r) - 1] dr} (r) - 1] dr} (8)

[ ∫g

SS

∞

]

]

ln(gSS) dr + T,P,NS

ln(gSS) dr - ∫(gSS - 1) dr (9)

]

where υ j A∞ is the infinite-dilution partial molar volume of the solute, and xA is the solute mole fraction. The arguments of the pair correlation functions have been dropped in this expression and in the remainder of this section to simplify notation. Thus, in addition to the ideal gas contributions, contributions to the partial molar entropy arise from solvent-solute pair correlations, perturbations in the solvent structure induced by the presence of the solute, and solvent expansion to accommodate the added solute at constant pressure. These last two contributions were neglected in previous studies,8-11 assuming for the hydration of simple solutes (e.g., methane) that they were nearly equal in magnitude and opposite in sign,9 and the infinitedilution partial molar entropy was calculated from

ShA∞ ≈

3kB j A∞ - kB ln(FAΛA3) + kBFSυ 2 kBFS

[ ∫g

SA

ln(gSA) dr - ∫(gSA - 1) dr (10)

]

Alternatively, the partial molar entropy can be derived from11

ShA∞ ) Sˆ A∞ + υ j A∞

(∂V∂S)

T,NS,NA)0

(11)

where the solute partial specific entropy, Sˆ A, is obtained by differentiating eq 8 at constant volume rather than constant pressure. In the limit of infinite dilution, this quantity is

Sˆ A∞ )

3kB - kB ln(FAΛA3) - kBFS ∫gSA ln(gSA) dr 2 kBFS2 ∂gSS ∞ (g 1) dr ∫ SA ∫ ∂FA T,V,N ln(gSS) dr (12) 2 S

[

]

( )

The first two terms in this expression are the ideal-gas contributions to the partial specific entropy, the third is the contribution from solvent-solute pair correlations, and the last term is the contribution due to solute-induced perturbations in

1902 J. Phys. Chem., Vol. 100, No. 5, 1996

Ashbaugh and Paulaitis

the solvent structure, which we define as the solvent reorganization entropy at the level of pair correlations. The solute chemical potential is calculated using either the partial molar or the partial specific entropy. At infinite dilution,

∆S* j A∞ A ) -kBTR° + kBFSυ

j A∞ + E h A∞ - TShA∞ µA∞ ) Pυ

The last term in this equation is the solvent-solute pair correlation entropy, hereafter denoted as S* SA.

) Pυ j A∞ + kBT ln(FAΛA3) - kBTFSυ j A∞ + SA

SA

SA

kBTFS2υ j A∞ 2

SS

SS

( )

kBTFS ∂gSS ∫ ∂xA 2

SS

∞

ln(ySS) dr (13) T,P,NS

or equivalently,

µA∞ ) Eˆ A∞ - TSˆ A∞

[ ∫g

) kBT ln(FAΛA3) + kBTFS

∫(gSA - 1) dr] +

SA

2

kBTFS 2

∫

ln(ySA) dr -

( ) ∂gSS ∂FA

∞

ln(ySS) dr (14) T,V,NS

where E h A∞ is the partial molar and Eˆ A∞ is the partial specific energy for the solute at infinite dilution, and yRβ is the cavity distribution function for species R and β, defined by22

yRβ ≡ exp(ΦRβ/kBT)gRβ

(15)

with ΦRβ the pair interaction energy. The equivalence of eqs 13 and 14 is established by calculating the pressure given the truncation of the entropy expansion in eq 7; i.e.,

P)T

(∂V∂S)

T,NS,NA)0

-

kBTFS2 ) kBTFS + 2

(∂V∂E)

T,NS,NA)0

[ ∫g

SS

ln(ySS) dr - ∫(gSS - 1) dr + 3

kBTFS 2

∫

( ) ∂gSS ∂FS

]

T,NS,NA)0

ln(ySS) dr (16)

It should be noted that the first term in this expression, when j A∞ in multiplied by υ j A∞, cancels identically with the term kBTFSυ ∞ eq 13. This cancellation does not occur if the Pυ j A term in eq 13 is neglected, as is frequently done for liquid solutions at low pressures. The connection between the partial molar entropy at infinite dilution and Ben Naim’s standard entropy of solvation,23

∆S*A ) -

∂µ* A ∂T

(17)

with 3 µ* A ) µA - kBT ln(FAΛA )

(18)

is given by

∆S*A ) ShA∞ + kB ln(FAΛA3) - kBTR° -

3kB 2

[ ∫g

SA

ln(gSA) dr - ∫(gSA - 1) dr (20)

]

Solvent Reorganization

[∫g ln(y ) dr - ∫(g - 1) dr] [∫g ln(y ) dr - ∫(g - 1) dr] +

kBTFS

kBFS

(19)

where R° the thermal expansion coefficient of the solvent. Substituting eq 10 into this expression, we obtain the desired result for the standard entropy of solvation in terms of solutesolvent pair correlations,

Solvent reorganization contributions to the entropy and the energy of solvation are difficult to evaluate using correlation functions derived from computer simulations because they arise from perturbations in solvent structure averaged over a large number of solvent molecules. In the case of hydrophobic hydration, these difficulties are magnified by the large number of degrees of freedom required to describe water-water orientational correlations. Using thermodynamic integration, Yu and Karplus24 derived coupling parameter expressions for the partial specific energy and partial specific entropy of solvation. The excess free energy of solvation was then expressed in terms of these properties, and contributions to the energy and the entropy that exactly cancel one another in this expression were identified as the solvent reorganization contributions. Solvent reorganization was defined in their analysis as solute-induced perturbations in the solvent radial distribution function, which give rise to a nonvanishing solvent-solvent interaction energy contribution24,25 and, similarly, the solventsolvent entropic contribution in eq 12. Since these contributions cancel in the free energy expression, they showed that only solvent-solute correlations are required to calculate solvation free energies. Guillot and Guissani26 also recognized that these solvent reorganization contributions to the solvation entropy and energy do not contribute to the free energy in their implementation of Widom’s test particle insertion method27 and, as such, focused on calculating only those contributions to the entropy due to solvent-solute correlations. Matubayasi et al.28 proposed an alternative decomposition of the correlation function expansion for the excess entropy, identifying a solvent reorganization contribution that is different from the one presented by Yu and Karplus, but incorrectly asserted that their dissemination of terms corresponds to that of Lazaridis and Paulaitis and likewise to that presented here. To some extent, the decomposition of the solvation free energy is semantic, as pointed out by Lee.29 Nonetheless, we believe that, to develop a mechanistic understanding of solvation thermodynamics, it is essential to identify those energetic and entropic contributions to the free energy that do compensate, no matter what they are called. To this end, we show below that those terms in the correlation function expansion for the solvation entropy identified as solvent reorganization contributions cancel exactly with similar solvent reorganization contributions to the solvation energy and, therefore, can be compared term by term with the free energy expressions of Yu and Karplus or Guillot and Guissani. In what follows, we consider the more general case of energy-entropy compensation starting from completely independent expressions for the partial specific energy and the partial specific entropy. This approach is based on the energy equation22 and the correlation function expansion for the entropy given above, for which solvent reorganization energy-entropy compensation has not been established. A general proof of compensation is complicated by the fact that the correlation function expansion is formally an infinite series, and exact compensation of those contributions identified as the solvent reorganization entropy with the solvent reorganization energy requires contributions from each term in this series. Nevertheless, it is possible to examine only leading order terms in the


J. Phys. Chem., Vol. 100, No. 5, 1996 1903

expansion for the entropy by considering the low-density limit of the correlation functions and infer from this analysis whether or not energy-entropy compensation holds at liquidlike densities. From eqs 3-5, the excess entropy associated with solventsolvent correlations in the pure fluid is given by

kBFS ∫∫gSS(2)(r1,r2) ln[gSS(2)(r1,r2)] dr1 dr2 + 2 kBFS2 ∫∫[gSS(2)(r1,r2) - 1] dr1 dr2 2 2

S*solvent ) -

Differentiating this expression with respect to NA at constant T, V, and NS yields the solvent reorganization entropy at infinite dilution (xA f 0),

( )

∞

) T,V,NS

(

)

FS2 ∂gSS(1,2) ∫∫ 2T ∂NA

dr1 dr2 -

(

∞

ΦSS(1,2) ×

)

[(

∞

[

(

) ] ∞

× T,V,NS

∞

kBFS3 ∂gSS(2,3) -6gSS(1,2) ∫∫∫ 6 ∂NA ∂gSS(1,2) 3 ∂NA

( ) ∂S* solvent ∂NA

×

T,V,NS ∞

) T,V,NS

(

)

FS2 ∂gSS(1,2) ∫∫ 2T ∂NA

∞

ΦSS(1,2) dr1 dr2 T,V,NS

kBFS3 ∫[(fSS(1,2) + 1)fSS(1,3) fSS(2,3) fSA(1,4) fSA(2,4)] × 2V kBFS3 ∫[(fSS(1,2) + 1)(fSS(1,3) + 1) × 2V (fSS(2,3) + 1)fSA(1,4) fSA(2,4)] dr1 dr2 dr3 dr4 -

dr1 dr2 dr3 dr4 +

dr1 dr2 dr3 + ... (22)

The first two terms in this equation are the solvent reorganization contributions at the level of pair correlations. It is easy to show that this contribution is equivalent to the solvent reorganization term in eq 12 for the partial specific entropy by applying the definition for the cavity distribution function in eq 15. The first term in eq 22 also corresponds to the change in the total solvent-solvent interaction energy upon dissolving the solute at constant volume. It is this term that compensates identically with the solvent reorganization energy. Thus, for complete energy-entropy compensation, all remaining terms in eq 22 must sum to zero. In the limit FS f 0, all terms in the expansion approach zero more rapidly than the first; hence, energyentropy compensation holds in the zero-density limit. At higher densities, the pair and triplet correlation functions can be expanded in terms of the density. For the pair correlation function, this expansion is30

gSS(1,2) ) (fSS(1,2) + 1) 1 + FS∫fSS(1,3) fSS(2,3) dr3 +

(

FA∫fSA(1,3) fSA(2,3) dr3 + ... (23)

)

)

(

)

FS2 ∂g (1,2) ∫∫ SS∂NA 2T

• ° +• ° • • • •

+

T,V,NS

T,V,NS

where fRβ(i,j) are Mayer f functions, defined by

∞

(26a)

∞

) ] ( )

∂ ln(gSS(1,2)) ln(δgSSS(1,2,3)) - 3gSSS(1,2,3) ∂NA dr1 dr2 dr3 +

gSSS(1,2,3) ) (fSS(1,2) + 1)(fSS(1,3) + 1)(fSS(2,3) + 1) × [1 + O(FS)] (25)

(fSS(1,2) + 1)fSA(1,4) fSA(2,4)] dr1 dr2 dr3 dr4 + O(FS4)

)

kBFS3 ∂gSSS(1,2,3) ln(ySS(1,2)) dr1 dr2 ∫∫∫ 6 ∂NA

(

The triplet correlation function is similarly expanded, giving

× T,V,NS

(24)

kBFS3 ∫[2(fSS(1,2) + 1)(fSS(2,3) + 1)fSA(2,4) fSA(3,4) 2V

T,V,NS

kBFS2 ∂gSS(1,2) ∫∫ 2 ∂NA

)

ΦRβ(i,j) -1 kBT

Substituting eqs 23 and 25 into eq 22, the solvent reorganization entropy to order FS3 becomes

kBFS3 (3) (3) (r1,r2,r3) ln[δgSSS (r1,r2,r3)] dr1 dr2 dr3 + ∫∫∫gSSS 6 kBFS3 (3) (2) (2) (r1,r2,r3) - 3gSS (r1,r2) gSS (r2,r3) + ∫∫∫[gSSS 6 (2) 3gSS (r1,r2) - 1] dr1 dr2 dr3 + ... (21)

∂S* solvent ∂NA

(

fRβ(i,j) ) exp -

∞


• °+ • ° + • • • • • • • ° +• ° +2 • • • • • •

kBFS3 × 2V

° + • ° +2 • ° +2 • ° + • • • • • • • ° +2 • ° + O(F 4) (26b) S • • •

–

where we have represented the integrals containing Mayer f functions diagrammatically13 in eq 26b. In these diagrams, the solid circles represent integrations over solvent coordinates, the open circles represent integrations over the solute coordinates, and the lines represent Mayer f bonds between circles acting on each other with the appropriate pair potential. Energyentropy compensation holds to order FS3 only if these diagrams sum to zero, which is obviously the case as diagrams from different integrals exactly cancel one another. Thus, the terms in the solvent reorganization entropy corresponding to pair, triplet, and higher order correlations are not zero by themselves, but every term in the expansion cancels with another term one order higher in the hierarchy such that, at least to order FS3,


∞

) T,V,NS

)

(

)

FS2 ∂gSS(1,2) ∫∫ 2T ∂NA

(

)

1 ∂E*solvent T ∂NA

∞


∞

(27) T,V,NS

Hence, energy-entropy compensation also holds at higher densities. We can extend this analysis to liquidlike densities by considering eq 22 in the context of Widom’s test particle insertion method. Our starting point is the Gibbs entropy formula,22

1904 J. Phys. Chem., Vol. 100, No. 5, 1996


S* NS ) -kB〈ln PNS〉NS

(28)

which relates the excess entropy to the probability density of observing a given configuration of NS solvent molecules. The brackets 〈 〉NS denote the canonical ensemble average over NSparticle configuration space. The canonical ensemble probability density is defined as

( )

exp P NS )

UNS

kBT ZNS

(29)

where UNS is the total potential energy of interactions between the NS molecules, and ZNS is the configurational partition function. The Gibbs entropy formula is the basis for the correlation function expansion and, therefore, is formally equivalent to eq 21. In the thermodynamic limit, NS f ∞, the derivative of eq 28 with respect to NA at constant T, V, and NS is given by the finite difference before and after the solute is added,

( ) ∂S* ∂NA

∞ T,V,NS

) -kβ〈ln PNS+1〉NS+1 + kB〈ln PNS〉NS (30)

The probability of observing a given configuration in the (NS + 1)-particle ensemble is

(

exp PNS+1 )

) ( ) ( ) 〈 ( )〉

UNS+1

kBT ZNS+1

exp -

)

∆UNS+1

exp -

kBT ∆UNS+1 exp kBT

UNS

kBT

)

Z NS

NS

P*NS+1PNS (31) where ∆UNS+1 is the potential energy of solute-solvent interactions, and P*NS+1 is the conditional probability of inserting the solute at a particular position in the solvent for a given configuration of the NS solvent molecules. Thus,

( ) ∂S* ∂NA

∞ T,V,NS

) -kB〈ln P* NS+1〉NS+1 - kB〈ln PNS〉NS+1 + kB〈ln PNS〉NS (32)

where the first term in this expression is the entropy due to solvent-solute correlations, and the last two terms represent the change in entropy due to perturbations in solvent structure, i.e., the solvent reorganization contribution. Equating these two terms to eq 22 gives, after simple algebraic manipulation,


∞ T,V,NS

) -kB〈ln PNS〉NS+1 + kB〈ln PNS〉NS 1 ) (〈UNS〉NS+1 - 〈UNS〉NS) T

(33)

Hence, energy-entropy compensation is shown to hold in general. An important implication of energy-entropy compensation in eqs 27 and 33 is that the solvent reorganization entropy can be evaluated indirectly by calculating the solvent reorganization energy. Thus, water-water orientational correlations required in the correlation function expansion for the entropy are not evaluated explicitly, but enter indirectly through the calculation of the orientation-dependent potential energy. Matubayasi et al.28 recently calculated the solvent reorganization contribution

to the partial molar energy for methane at infinite dilution in TIP4P water at 10 °C by computing the water-water interaction energy as a function of radial distance from the solute. The solvent reorganization energy was then obtained by integrating this quantity over one or more hydration shells. Their calculated value of -1.7 kcal/mol corresponds to a solvent reorganization partial molar entropy that is approximately 30% of the value we calculated for the standard entropy of hydration of methane in TIP4P water at 10 °C using eq 20.8 Thus, at least for methane hydration, this contribution is significant, but not the dominant one. Extension to Molecular Solutes of Arbitrary Shape For molecular solutes in water, the pair correlation function will depend on water-solute separation, relative orientation, and solute conformation. The additional number of degrees of freedom required to describe these orientations and solute conformations makes rigorous determination of the water-solute pair correlation function intractable, even for relatively simple molecular solutes. Thus, we propose an approximation to this pair correlation function to calculate pair correlation entropies of hydration for molecular solutes of arbitrary shape,

S* SA ) -kBFS[∫gSA ln(gSA) dr - ∫(gSA - 1) dr] (34) To do this, we note that solvent-solute pair correlation functions can be defined relative to the solute center of mass (COM) or, through a coordinate transformation, to any other site on the solute. It follows, therefore, that solute correlations with different solvent molecules can be defined using different solute sites. Ethane in the united atom approximation, for example, can be treated as two atomic sites corresponding to the centers of mass of the two carbon atoms. The correlation function integrals are also invariant to coordinate transformations, and as a consequence, the pair correlation entropy can be expressed equivalently as

{ ∫g

S*SA ) -kBFS

{ ∫g

) -kBFS

SA(rCOM)

SA(rR)

ln[gSA(rCOM)] drCOM -

∫[gSA(rCOM) - 1] drCOM}

ln[gSA(rR)] drR -

∫[gSA(rR) - 1] drR}

(35)

where gSA(rCOM) is the solvent-solute pair correlation function relative to the solute COM, rCOM is the position of a solvent molecule relative to the solute COM coordinates, gSA(rR) is the solvent-solute pair correlation function relative to an arbitrary solute site R, and rR is the position of a solvent molecule relative to site R. Furthermore, the correlation function integrals can be evaluated by dividing the domain of integration into a number of subdomains and performing the integration in each subdomain with respect to a different solute size. We define the subdomains by a Voronoi tesselation of space with nodes centered on the solute sites. Thus, the number of subdomains equals the number of solute sites.31 For example, the two sites on ethane define two subdomains divided by the plane midway between the carbon atoms and normal to the carbon-carbon bond vector (Figure 1a). This definition of subdomains can also be readily generalized to solutes of arbitrary shape. Schematic representations of the subdomains for propane and trans-butane are shown in Figure 1b. The solvent-solute pair correlation function in each subdomain is evaluated relative to the solute site defining that


J. Phys. Chem., Vol. 100, No. 5, 1996 1905

Figure 2. Definition of the angles θ and χ describing the orientation of water with respect to proximal solute site R. θ is the angle between vectors rR and rµ, while χ is the angle between vectors rHH and rRµ.

Figure 1. (a, top) Schematic representation of the two subdomains for ethane defined by a Voronoi tessellation of space in which the nodes are centered on the two carbon sites. Ethane in the united atom approximation is treated as two sites corresponding to the centers of mass of the two carbon atoms. (b, bottom) Schematic representations of the subdomains for propane and trans-butane defined by the same method and making the united atom approximation in treating the alkanes.

subdomain, such that the pair correlation entropy is given by Rmax

S*SA ) -kBFS ∑

R)1

{∫ g

(r ) VR SA R

ln[gSA(rR)] drR -

∫V [gSA(rR) - 1] drR} R

(36)

where Rmax is the total number of solute sites and VR is the subdomain volume of integration for site R. Equation 36 suggests that a site additivity approximation can be used to calculate S* SA if the solvent-solute site pair correlation functions for the individual solute sites are independent of one another; i.e., if these pair correlation functions are independent of the molecular solute. If the solute sites are not independent, a collection of sites can be defined such that this group is independent of other sites or groups of sites comprising the molecular solute, which naturally leads to a group contribution approximation for calculating the pair correlation entropy. The complexities of evaluating solvent-solute site pair correlation functions within subdomains have not been removed at this point since these pair correlations still depend on solventsolute site separations, relative orientations, and molecular conformations. These complexities are avoided, however, by assuming the pair correlation function within each subdomain is independent of solute orientation relative to the solvent molecule; i.e., that the solute site is “locally” spherically symmetic. Thus, if the solvent molecule is also spherically symmetric, the solvent-solute pair correlation function will be a function of separation alone, prox (rR) gSA(rR) ≈ 〈gSA(rR)〉rR ≡ gSR

prox main associated with site R, and gSR (rR) is defined as the proximal distribution function. The integrated sum of the proximal distribution functions over all solute sites is subject to the same normalization conditions as gSA(rR).32 The proximity approximation given by eq 37 is an extension of a similar approach to characterizing solvent structure around molecular solutes proposed by Ben-Naim31 and implemented by Mehrotra and Beveridge32 to calculate coordination numbers and hydration energies. However, in our implementation of this approximation to calculate pair correlation entropies of hydration, water orientations relative to the individual solute sites are still taken into account. The two angles needed to describe these orientations are shown in Figure 2. The two angles are evaluated over the range, 0-π. The range for χ has been reduced by a factor of 2 due to the symmetry of water. To analyze contributions to the entropy from water-solute site orientational correlations, we use the same factorization of prox (rR,θ,χ) applied previously to simple hydrophobic solgSR utes.8 That is, orientational correlations are independent of distance within the first hydration shell of the solute site and are negligible beyond the first shell. Thus,

(37)

where 〈gSA(rR)〉rR is averaged over solute orientations relative to a solvent molecule constrained to remain within the subdo-

prox gSR (rR,θ,χ)

)

{

prox prox (θ,χ) gSR (rR) pSR

rR e rsh R

prox gSR (rR)

rR > rsh R

(38)

prox (θ,χ) requires where the normalization condition for pSR

∫∫pSRprox(θ,χ) sin θ dθ dχ ) Ω ) 4π

(39)

A discussion of eq 38 is given elsewhere.10 When this factorization is coupled with eq 37, the following separation of translational and orientational contributions in the site-additive expression for the pair correlation entropy is obtained, Rmax

S*SA ) -kBFS ∑

R)1

∫

[gprox(rR) VR SR

{∫ g

prox (rR) VR SR

}

- 1] drR -

prox ln[gSR (rR)] drR -

kBFSRmax

prox Vsh ∑ { R ∫∫pSR (θ,χ) × Ω R)1

} (40)

prox (θ,χ)] sin θ dθ dχ ln[pSR sh

rR prox where FSVsh R ) FS∫0 gSR (rR) drR is the coordination number of site R. Equation 40 allows a qualitative discussion of the observed solute surface area dependence of thermodynamic properties of aqueous solutions containing nonpolar solutes. It has been shown experimentally that many of these thermodynamic properties (e.g., free energy, enthalpy, heat capacity) can be correlated with the solvent-accessible surface area of the

1906 J. Phys. Chem., Vol. 100, No. 5, 1996


solute.2,33 These correlations typically require as input the hardsphere radii of water and the atoms comprising the solute to calculate the solvent-solute contact surface area. Hard-sphere radii for water range from 1.35 to 1.5 Å and, for carbon, from 1.75 to 1.85 Å.33,34 Although surface areas determined over the range of these values are not that different, the uncertainties in hard-sphere radii lead to contact surface areas that are not uniquely defined. An alternative approach has been to use the coordination number of solvent molecules around the solute, rather than the surface area, in these correlations.35 However, despite the general success of such approaches, there is little theoretical justification for either the solvent-accessible surface area or the coordination number dependence of hydration thermodynamic properties. In fact, eq 40 suggests that not one, but several, surface areas (weighted by the proximal distribution functions) must be considered in such correlations. Since the integrals in eq 40 are dominated by contributions from the first hydration shells of the solute sites, it follows that contact surface areas corresponding to the first peak in the water-solute site radial distribution functions should correlate strongly with the entropy of hydration. Therefore, empirical surface area correlations are expected to be reasonably successful for estimating entropies of hydrophobic hydration as long as the entropy is not sensitive to subtle details of the chemical structure or molecular conformation of the solute. A practical description of the solvation thermodynamics of molecular solutes must address the effects of the connectivity of polyatomic solutes and the molecular conformation of the solute. The proximity approximation leading to eq 40 makes use of the reasonable expectation that the entropy of hydration of a molecular solute composed of equivalent sites is dominated by pair correlations between the individual water molecules and the nearest solute site. Therefore, proximal distribution functions for the same site on similar solute moleculessfor example, the methyl group in normal and branched alkane chainssare expected to be equivalent. Structural details, such as solute connectivity and molecular conformation, are accounted for in eq 40 by differences in the Voronoi topologies of the subdomains for different solutes or solute conformations and by the indirect effects of neighboring solute sites on the proximal distribution functions evaluated from molecular simulations. Thus, eq 40 will predict different hydration entropies for isomers based on differences in the connectivity of the solute sites. The proximity approximation assumes that the solvent-solute site correlation function at a given radial distance from a particular solute site can be represented by an average value for the entire subdomain corresponding to that site. If the actual pair correlation function at this radial distance deviates from the average value by δ, then the deviation in the entropy obtained from eq 40 will be on the order of δ2. Therefore, for small values of δ, the uncertainty in the calculated entropy will be much less than the uncertainties in the associated pair correlation functions calculated on the basis of the proximity approximation. To test whether deviations in the actual pair correlation functions are indeed small or not, site-site radial distribution functions for water oxygens around individual solute sites are predicted using the proximal distribution functions and then compared to those obtained directly from simulations. The site-site radial distribution function is defined as the angleaveraged site-site pair correlation function, rdf gSR (rR) )

1 ∫g (r ) dω 4π SA R

(41)

where ω represents the angular coordinates and the integration is performed over a spherical shell centered on site R. Unlike

TABLE 1: Interaction Parameters and Geometrical Constraints of Simulated Alkanesa and Waterb Lennard-Jones and Coulombic Interaction Parametersc σ

molecular site alkane carbon TIP4P water oxygen hydrogen M (auxiliary site)

q

3.73

0.294

0.0

3.15365 0.0 0.0

0.155 0.0 0.0

0.0 +0.52 -1.04

Molecular Bond Lengths bond

r, Å

alkanes C-C TIP4P water O-H O-M

1.53 0.9572 0.15

Molecular Bong Angles bond angle

angle, deg

alkanes C-C-C C-C-C-C dihedral angle gauche-butane trans-butane TIP4P water H-O-H a

b

110.0 60 180 104.52

c

See ref 38. See ref 39. The units are Å for σ and kcal/mol for , and q is expressed in terms of the fundamental electrostatic charge.

proximal distribution functions, solvent molecules in the sitesite radial distribution functions are not constrained to be in the subdomain for the solute site of interest. For a solute consisting of one site, the radial and proximal distribution functions are equivalent; however, for molecular solutes, the entire set of proximal distribution functions for the molecule is required to predict the site-site radial distribution functions. One merely substitutes the proximal distribution function representation of gSA(rR) into the integral in eq 41. Simulation Method and Computational Details Computer simultions were used to obtain the distribution prox prox functions, gSR (rR) and pSR (θ,χ), needed to calculate the pair correlation entropy via eq 40. Monte Carlo (MC) simulations36 of one solute at infinite dilution in water were performed in the isothermal/isobaric ensemble using Jorgensen’s BOSS program37 (version 2.8) modified to calculate the required distribution functions. A series of normal alkanes from methane to normal butane were studied. Alkanes were simulated as a chain of identical LJ particles modeled using the methane OPLS potential parameters.38 The TIP4P potential model39 was used for water. All site-site interaction potential energies were calculated with a spherical cutoff of 8.5 Å. The potential parameters and structural parameters are given in Table 1. The starting point for each simulation was a box of 267 water molecules preequilibrated at 25 °C and 1 atm (216 water molecules were used for methane). After inserting the solute, one to three water molecules were deleted from the simulation box to eliminate overlaps in the initial calculation of intermolecular potential energies. Each MC run consisted of 30 million configurations for equilibration, followed by an additional 100 million configurations for analysis. Block averages of 20 million configurations were used to perform statistical analyses on the results. All simulations were performed at a constant temperature of 25 °C and a constant pressure of 1 atm. The simulations were run on an IBM RISC/6000 Model 590 computer. On this machine, 15 million configurations took


J. Phys. Chem., Vol. 100, No. 5, 1996 1907

a

d

b

e

c

Figure 3. Calculated water oxygen-solute site proximal distribution functions for (a) methane, (b) ethane, (c) propane, (d) gauche-butane, and (e) trans-butane. The open circles denote terminal carbon sites, and the solid circles denote interior sites. For methane, the proximal distribution function and the radial distribution function are equivalent.

approximately 24 CPU h. The preferential sampling facility of BOSS was employed.40 Proximal and site-site radial distribution functions were calculated with a radial increment of 0.1 Å. The proximal distribution functions were evaluated by counting the number of water molecules in the individual spherical shells around the proximal solute site, but only within the shell volume of the subdomain for that site. These populations were then divided by the number of water molecules in each subdomain shell at the bulk density of water. The integrals for the translational entropy were evaluated using rectangular midpoint integration. More refined numerical integration schemes were also investigated but did not lead to significant differences in the calculated value of the translational entropy. Although the integrals were evaluated to a radial distance of 8.5 Å, it was found that they converged quickly beyond 7 Å to within the statistical uncertainties of the simulation. As dictated by eq 38, orientational correlations were calculated only in the first hydration shell of

each solute site and were neglected beyond that distance. The first hydration shell for a solute site extends to the first minimum prox (rR) for that site: approximately 5.4 Å for all sites. The in gSR orientational entropy was evaluated using intervals of 3°, 6°, 9°, and 12° for both θ and χ and then optimized as described in our original study.8 As before, the interval of 6° was found to give an optimal value of the orientational entropy for all solute sites. Results and Discussion The calculated water oxygen-solute site proximal distribution functions for the n-alkanes, including both the trans and gauche conformations of butane, are shown in Figure 3. The first peak prox (rR) appears at a separation of 3.8 Å in all cases, and the in gSR peak height increases with increasing alkane chain length for methane, ethane, and propane but is essentially the same for the terminal carbon sites on the chain ends of propane and the

1908 J. Phys. Chem., Vol. 100, No. 5, 1996

a

Ashbaugh and Paulaitis TABLE 2: Orientational Correlation Entropy per Water Molecule for Individual Alkane Sites Calculated via Eq 45a alkane site

orient -SSR

methane ethane propane terminal carbons interior carbon gauche-butane terminal carbons interior carbons trans-butane terminal carbons interior carbons

0.302 ( 0.040 0.326 ( 0.032 0.326 ( 0.040 0.368 ( 0.034 0.332 ( 0.046 0.368 ( 0.040 0.348 ( 0.040 0.428 ( 0.022

a The units are eu/water (eu ) cal/(mol K)). The uncertainties are reported as one standard deviation.

b

c

Figure 4. Contour plots of the calculated orientational probability distribution p(θ,χ) of water molecules in the first hydration shell of (a) methane, (b) proximal to the terminal carbon sites of trans-butane, and (c) proximal to the interior carbon sites of trans-butane. Letters denote characteristic orientations: (A) hydrogen radial inward (θ ) 125°, χ ) 90°); (B) hydrogen radial outward (θ ) 55°, χ ) 90°); (C) lone pair radial outward (θ ) 125°, χ ) 0°); (D) lone pair radial inward (θ ) 55°, χ ) 0°).

two butane conformations. For trans- and gauche-butane, the interior carbons have significantly higher first peaks than those for the terminal carbons, indicating that these sites have a higher density of water molecules in their first hydration shells relative to the chain ends. There is, however, little difference in first peak heights for the interior and terminal carbons of propane. The density of water molecules is also higher around the interior carbons of butane in the trans conformation compared to the those carbons in the gauche conformation. Since methane potential parameters were used for all sites on the chains, rather

than optimized OPLS parameters for the terminal (CH3-) and interior (CH2-) sites, the observed increase in peak heights with increasing carbon number might be attributed to differences between the OPLS parameters for these groups and those for methane.38 To evaluate this possibility, we also performed an identical simulation for propane using instead the OPLS parameters recommended for CH3- and CH2- and found that the first peak heights from this simulation were only marginally different from those obtained using the methane potential parameters. This result was expected since water structure around a LJ solute is governed largely by the repulsive core of the solute,41 and the two parameter sets give essentially the same values for the repulsive cores of the different sites, as characterized by their calculated thermal radii.10 The higher first peak prox (rR) for the interior carbons of butane relative to the in gSR chain ends reflects a deeper potential energy minimum for water molecules in the vicinity of these sites due to additional van der Waals interactions with the two neighboring carbon sites. This conclusion is supported by the observation that first peak heights for sites on the chain ends do not change significantly from propane to butane, while they do increase substantially for the interior carbons. Another contributing factor may be differences in the structure or "local curvature” of the solute at the interior sites relative to the chain ends.42,43 The calculated water orientational distributions around the solute sites were found to be essentially the same for all the alkanes. Contour plots of these distributions for methane and the terminal and interior sites of trans-butane are given as examples in Figure 4. The dominant feature of these contour plots is the maxima at B and C, corresponding to preferred orientations in which one vertex of the tetrahedral hydrogenbonding structure of water points directly away from the carbon while the other three vertices “straddle” the site to optimize hydrogen bonding with other water molecules. Conversely, the minima at A and D correspond to the least favorable orientations in which water forfeits one hydrogen bond by pointing one vertex directly toward the carbon. A quantitative comparison of water orientational distributions around the different alkane sites is made by calculating the orientational entropy per water molecule in the first hydration shell of each site. From eq 40, this entropy for an arbitrary site R is orient SSR )-

kB ∫∫pSRprox(θ,χ) ln[pSRprox(θ,χ)] sin θ dθ dχ (42) Ω

The results are summarized in Table 2. Within the statistical uncertainties reported in the table, the orientational entropies per water molecule are equivalent for all sites, except the interior carbons of trans-butane, where it is more negative. Orientational entropies for the other interior carbons of propane and


J. Phys. Chem., Vol. 100, No. 5, 1996 1909

a

e

b

f

c

g

d

Figure 5. Comparison of calculated and predicted water oxygen-solute site radial distribution functions. The open circles denote radial distribution functions obtained from the simulations, and the solid circles denote radial distribution functions predicted from eq 41. The solute sites are (a) ethane, (b) the terminal sites of propane, (c) the interior site of propane, (d) the terminal sites of gauche-butane, (e) the interior sites of gauchebutane, (f) the terminal sites of trans-butane, and (g) the interior sites of trans-butane.

gauche-butane are also consistently more negative than those for the terminal carbons: however, these differences are not

statistically significant. It is interesting to note that waters of hydration around the interior carbons of trans-butane, which

1910 J. Phys. Chem., Vol. 100, No. 5, 1996


TABLE 3: Translational and Orientational Contributions to the Pair Correlation Entropies of Hydration Calculated via Eq 40 on a per Site Basis and for the Entire Solutea alkane site methane ethane site total propane terminal carbons interior carbon total gauche-butane terminal carbons interior carbons total trans-butane terminal carbons interior carbons total

translational entropy

orientational entropy

total entropy

10.9 ( 0.9

6.4 ( 0.6

17.3 ( 1.1

7.9 ( 0.4 15.8 ( 0.8

3.9 ( 0.4 7.8 ( 0.8

11.8 ( 0.6 23.6 ( 1.1

7.4 ( 0.1 4.8 ( 0.4 19.6 ( 0.5

3.5 ( 0.4 2.3 ( 0.2 9.4 ( 0.9

10.9 ( 0.5 7.1 ( 0.5 29.0 ( 1.0

6.8 ( 0.4 4.8 ( 0.3 23.2 ( 1.0

3.2 ( 0.4 2.1 ( 0.2 10.6 ( 1.0

10.0 ( 0.6 6.9 ( 0.4 33.8 ( 1.4

7.7 ( 0.6 4.8 ( 0.4 25.0 ( 1.4

3.8 ( 0.4 2.0 ( 0.2 11.7 ( 0.9

11.5 ( 0.7 6.8 ( 0.4 36.7 ( 1.7

a The values for the entropy are eu. Uncertainties are reported as one standard deviation.

have the least translational freedom, also have the least orientational freedom compared to water around the other alkane sites. Site-site radial distribution functions for water oxygens around the different alkane sites obtained directly from the simulations are compared to those predicted on the basis of proximal distribution functions in Figure 5. In general, the agreement is very good. The position of the first peak and the peak heights are nearly the same in all cases. However, the proximal distributions tend to underpredict slightly the magnitude of the first peak height with the largest discrepancy observed for the interior carbons of propane and butane. These underpredictions, which are systematic, occur because the proximal distribution function represents an average value of the actual pair distribution function for an entire subdomain and, as such, does not account for higher water densities that exist near the subdomain interfaces. The calculated first peak position for the interior carbons are also shifted to slightly larger separations relative to those for the terminal carbons. This shift reflects the influence of waters in the first hydration shells of neighboring carbons on the sitesite radial distribution function. The proximal distribution functions capture this shift, although the predicted first peak is somewhat more asymmetric in shape. This is also a consequence of neglecting higher water densities near the subdomain interfaces. The long-range behavior of the site-site radial distribution functions is also reasonably well represented in all cases, except water oxygens around the interior carbons of transbutane. The discrepancies in the predicted and calculated sitesite radial distribution functions noted above are, nonetheless, minor. Thus, we conclude that the proximal distributions give an accurate representation of water structure around the carbon sites of these alkane chains. Calculated translational and orientational pair correlation entropies of hydration for the individual n-alkane sites are given in Table 3. As expected, the total pair correlation entropy becomes increasingly more negative with increasing alkane chain length. The orientational contribution to the pair correlation entropy is also approximately 30% of the total value in each case. However, the most striking feature of these results is the nearly constant contributions for the different sites. That is, excluding methane, an average value of approximately -11 eu can be assigned to the pair correlation entropy of hydration for a terminal carbon site, and approximately -7 eu can be

TABLE 4: Standard Entropies of Hydrophobic Hydration for n-Alkanes through Butane Calculate via Eq 20a alkane species

b -∆S* A (experimental)

υ j A∞

-∆S* A (calculated)

methane ethane propane butanec

15.9 20.1 22.8 26.0

60 94 127 161

13.9 ( 1.1 18.0 ( 1.1 21.2 ( 1.1 25.4 ( 1.6

a The values for entropy are eu, and the values for υ j A∞ are Å3/molecule. Uncertainties are reported as one standard deviation.b See ref 44. c Butane is assumed to be a 56/44% trans/gauche mixture.46

assigned to the pair correlation entropy of hydration for an interior carbon site. The standard entropies of hydration of the n-alkanes can be calculated from the pair correlation entropies in Table 3 using eq 20. This calculation requires the infinite-dilution partial molar volume of the solute in water, which was estimated from the value for methane in TIP4P water (υ j ∞methane ) 60.1 Å3/ 45 molecule ) scaled by the ratio of the sum of hard-sphere volumes for the alkane sites divided by the hard-sphere volume for methane. A radius of 1.9 Å is used to calculate the hardsphere volume of a site. The thermal expansion coefficient of TIP4P water, also needed in eq 20, was obtained from the literature (R° ) 94 × 10-5 K-1).39 The results are compared with experimental values for the standard entropy of hydration at 25 °C in Table 4. The agreement is surprisingly good, although differences between the calculated and experimental entropies become more pronounced with increasing alkane chain length if these values are scaled by the values for methane. The agreement must result, in part, from a fortuitous cancellation of neglected terms. At the level of pair correlations, these terms are presumably the solvent reorganization and solvent expansion contributions to the partial molar entropies in eq 9. That these contributions do not appear to become more important with increasing alkane chain length, as is suggested by the results in Table 4, is noteworthy. The pair correlation entropies of hydration given in Table 3 for the interior carbons of butane are independent of solute conformation, even though waters around the interior carbon sites of trans-butane are more localized and have less orientational freedom. The observed conformational invariance of the entropy, including its translational and orientational components, occurs because the subdomains for the interior sites of transbutane are also smaller than those for the interior sites of gauchebutane. The two effects offset one another to such an extent that the pair correlation entropy of hydration of the interior sites remains essentially constant, independent of conformation. In contrast, pair correlation entropies of hydration for the terminal carbons (Table 3) do depend on conformation. This conformation dependence is also observed for both the translational and orientational entropies. In this case, however, the water-solute site pair correlations are approximately the same for both conformations, but the subdomains are smaller for the terminal carbons in the gauche conformation. The gauche conformation is, therefore, entropically favored over the trans conformation of butane as a consequence of smaller subdomains for the chain ends of gauche-butane. The trans-butane conformation is energetically more favorable than the gauche conformation in an ideal gas, which leads to a 68/32% population distribution of trans/gauche conformations. However, in TIP4P water, a stabilization of the gauche conformation occurs such that the trans/gauche population shifts to 56/44%.46 In contrast, stabilization of the gauche conformation is not observed in liquid butane.38 Thus, it appears that the conformational stabilization in aqueous solution is driven by hydrophobic interactions. As shown above, the entropic


Figure 6. Comparison of the calculated free energy of hydration of butane with the free energy contribution from the solvent-solute correlational entropy as a function of butane conformation. The molecular conformation of butane is defined by the C1-C2-C3-C4 dihedral angle, Φ, where Φ ) 60° and Φ ) 180° correspond to the gauche and trans conformations, respectively. The open circles denote calculated values of -T∆SSA, and the solid circles denote free energies calculated via statistical perturbation theory by Jorgensen and Buckner.46 These values are calculated relative a zero free energy for the trans conformation.

penalty for hydrating an interior site is essentially the same for both conformations, with an entropic benefit arising from smaller subdomains for the chain ends in the gauche conformation; i.e., the terminal sites are shielded from the aqueous environment to a greater extent in the gauche conformation. Unfortunately, the statistical noise in the entropy calculations (e.g., see Table 3) precludes obtaining a meaningful value for the entropy change of the trans f gauche transition in water. However, making the reasonable assumptions that the proximal and orientational distributions for water around the chain ends are independent of conformation and that the total entropy of hydration per interior site is also independent of conformation (from Table 3, it is -6.9 eu for gauche-butane and -6.8 eu for trans-butane), we can calculate the change in the entropy of hydration for butane as a function of its primary (C1-C2-C3-C4) dihedral angle. As shown in Figure 6, the gauche conformation is indeed entropically stabilized relative to the trans conformation in aqueous solution. When compared to the calculated conformational free energy of hydration,46 the larger negative change in the entropic contribution to the free energy, relative to the change in free energy itself, indicates that the trans f gauche transition is energetically unfavorable. However, this energetic penalty is overcome by the favorable entropy change, consistent with the notion of hydrophobic interactions. Similar observations were made in an earlier simulation study of the trans f gauche transition for normal butane in TIP3P water.47 Assuming the entropy of hydration per water molecule is independent of the solute site leads to a surface area-dependent correlation for the entropy of hydration of molecular solutes. This assumption also implies a simple linear relationship between the coordination number of water molecules around a solute and its entropy of hydration. We have generated this correlation for the n-alkanes studied here by calculating the coordination number of waters around the alkane chains from the sum of coordination numbers for the individual alkane sites. The site coordination numbers were obtained by integrating the proximal distribution functions for each site over the first hydration shell. As shown in Figure 7, this correlation is linear, suggesting that the entropic contributions of the solute are indeed additive. However, closer inspection of the calculated entropies of hydration for the individual solute sites reveals significant

J. Phys. Chem., Vol. 100, No. 5, 1996 1911

Figure 7. Correlation of solvent-solute pair correlation entropies (eq 40) with alkane coordination numbers calculated using the proximal distribution functions. The numbers denote (1) methane, (2) ethane, (3) propane, (4) gauche-butane, and (5) trans-butane.

TABLE 5: Coordination Numbers and Pair Correlation Hydration Entropy per Water Moleculea alkane site

coordination number

entropy

methane ethane propane terminal carbons interior carbon gauche-butane terminal carbons interior carbons trans-butane terminal carbons interior carbons

20.0 12.0

-0.86 -0.99

10.7 6.0

-1.01 -1.18

9.5 5.7

-1.05 -1.22

10.8 4.8

-1.06 -1.43

a

The values for entropy are eu/water.

differences, which become apparent when the entropy of hydration per water molecule is calculated for each site; i.e., the entropy of hydration for the individual sites is divided by the site coordination number. These results are given in Table 5. The entropy of hydration per water molecule for the terminal carbons, excluding methane, is nearly constant at roughly -1.0 eu per water molecule. The interior carbon sites of propane and gauche-butane contribute roughly -1.2 eu per water molecule, while the interior carbon sites of trans-butane contribute -1.43 eu per water molecule, consistent with the higher water-solute site orientational and translational correlations observed for that site. The large difference between the entropies of hydration of the interior sites of trans- and gauchebutane leads to the entropic stabilization of the gauche conformation in water, which cannot be predicted on the basis of a surface area-dependent correlation alone given an average value for the entropy of hydration per water molecule for an interior site of butane. A similar conclusion was reached in earlier work.35 Moreover, the results in Table 5 emphasize that surface area-dependent correlations for thermodynamic properties of hydration will not capture those features of water structure that are specific to the alkane chain conformation. Conclusions The extent to which the proximal distribution functions accurately predict site-site radical distribution functions obtained from simulation and the observed similarities in the proximal and orientation distribution functions for different alkane sites suggest that proximal distribution functions have general utility as the basis for predicting hydration thermodynamic properties of n-alkanes. These thermodynamic properties

1912 J. Phys. Chem., Vol. 100, No. 5, 1996


include the chemical potential, which can be calculated directly from the water-solute pair correlation entropy and energy. That is, the water oxygen-solute site radial distribution functions, predicted from the proximal distribution functions, can be used to calculate the water-solute interaction energy. Given this quantity and the water-solute pair correlation entropy, the solute chemical potential is easily obtained. It should be noted that this approach to deriving the chemical potential avoids the inefficiencies of Widom’s particle insertion method at liquidlike densities48 or the need to consider a large number intermediate hydration states in the application of statistical perturbation theory.49 Also, since the normalization of the proximal distribution functions is the same as that for the complete distribution function, the proximal distribution functions can be applied to calculating solute partial molar volumes at infinite dilution via the Kirkwood-Buff expression,49 i.e.,

υ j A∞ ) kBTκ° -

Rmax

∫V [gSRprox(rR) - 1] drR ∑ R)1 R

(43)

where κ° is the solvent isothermal compressibility. As an example, the water-solute pair interaction energy for methane in TIP4P water at 25 °C, determined using the radial distribution function obtained from our simulations, is -3.22 ( 0.04 kcal/ mol. Combining this value with the water-methane pair correlation entropy in Table 3 gives an excess chemical potential of 1.94 ( 0.33 kcal/mol, which compares favorably with the value calculated via statistical perturbation theory for the same system: 2.5 ( 0.4 kcal/mol.49 The experimental value is 2.00 kcal/mol.44 Moreover, agreement between the two calculated chemical potentials can be improved to give virtually the same prox (rR,θ,χ) in eq 38 is not value if the factorization of gSR invoked. Work is now in progress to calculate chemical potentials for the other normal alkanes in TIP4P water, which can be compared to those obtained from statistical perturbation theory, as well as to experimental values. The principal advantage proximal distribution functions have over, for example, solvent-accessible surface area as the basis for predicting thermodynamic properties is that they give detailed information about solvent structure in the vicinity of a solute as a function of the chemical structure and molecular conformation of the solute, and yet this information is embodied in a relatively small number of well-characterized distribution functions. In applications such as the one considered above, namely, the trans f gauche transition for butane in water, it is the conformational dependence of the proximal distribution functions that leads to the entropic stabilization of gauchebutane, which cannot be predicted on the basis of surface area correlations alone. Calculation of the entropy change (or the free energy change) for this transition can also be made much less computationally demanding by considering only the stable conformational states of interest, i.e., trans- and gauche-butane. Other methods capable of distinguishing solute chemical structures and molecular conformations, such as free energy perturbation, are computationally much more intensive because thermodynamic averages are required at a significant number of intermediate conformational states along the trajectory that defines the transition. The application of proximal distribution functions are also likely to have advantages when considering the effects of molecular roughness or curvature on the hydration of hydrophobic surfaces. In this case, the surface is characterized by how it is divided into subdomains, while hydration is characterized by the proximal distribution functions within the subdomains and how they vary with the surface features. While the solvent reorganization contributions to solvation are not required to determine the chemical potential, it has been

assumed that these contributions play only a small part in the determination of the partial molar entropy of hydration of simple hydrophobic solutes, such as methane. However, a more thorough study of the phenomenon needs to be undertaken, especially in light of the results for standard entropies of hydration of the n-alkanes in Table 4. First, a general proof of energy-entropy compensation within the context of the correlation function expansion for the entropy must be found. The method proposed by Matubayasi et al.28 to calculate the solvent reorganization contribution to the partial molar energy can also be generalized to molecular solutes, thereby providing an indirect evaluation of the solvent reorganization contribution to the partial molar entropy. This approach will allow calculation of the complete partial molar entropy of hydrophobic hydration of molecular solutes that will ultimately lead to relating the entropy of hydration to water structure around these molecular solutes. The present analysis has made use of proximal distribution functions to determine the water-solute pair correlation to the entropy of hydrophobic hydration for a series of n-alkanes. Compared to the general case where the solute consists of sites having a range of sizes and electrostatic charges, the sites of the n-alkanes are, to very good approximation, equivalent. In the case of solutes comprised of heterogeneous sites, the utility of proximal distribution functions in calculating entropies of hydrationsor, for the matter, energies or free energies of hydrationsis questionable, since the influence of neighboring sites is expected to be large, leading to wide variations in water structure around the individual solute sites. Thus, if this approach is to be applied to calculating thermodynamic properties for a wide range of solutes, it will be necessary to reconsider the underlying assumptions that lead from the site-additive expression for the entropy in eq 36 to eq 40, in which we have implemented proximal distributions functions, to determine how these assumptions can be modified for solutes having heterogeneous sites. Nonetheless, proximal distribution functions for heterogeneous solutes do have utility in applications other than the calculation of thermodynamic properties. These include, for example, applying proximal distribution functions to study the local solvation environment of an individual solute site surrounded by different groups of heterogeneous sites or to characterize the effects of solvation on transferring a single site from dilute solution, where the proximity analysis is exact, to a heterogeneous molecular solute or a surface. Acknowledgment. We are indebted to Dr. Lawrence Pratt and Prof. Eric Kaler for numerous helpful discussions. This work is supported the National Aeronautics and Space Administration (NAG3-1424), the National Science Foundation (Grant BCS9210401), and a National Science Foundation Fellowship for H.S.A. (Grant GER9253850). References and Notes (1) Frank, H. S.; Evans, M. W. J. Chem. Phys. 1945, 13, 507. (2) Privalov, P. L.; Gill, S. J. AdV. Protein Chem. 1988, 39, 191. (3) Frank, F.; Reid, D. S. In Water, a ComprehensiVe Treatise; Franks, F., Ed.; Plenum Press: New York, 1973; Vol. 2, Chapter 5. (4) Tanford, C. The Hydrophobic Effect, 2nd ed.; Wiley: New York, 1980. (5) Tucker, E. E.; Lane, E. H.; Christian, S. D. J. Solution Chem. 1981, 10, 1. (6) (a) Nettleton, R.; Green, M. S. J. Phys. Chem. 1958, 29, 1365. (b) Raveche´, H. J. Chem. Phys. 1971, 53, 2242. (7) (a) Wallace, D. C. Phys. ReV. A 1989, 39, 4843. (b) Green, H. S. The Molecular Theory of Fluids; North-Holland: Amsterdam, 1952. (8) Lazaridis, T.; Paulaitis, M. E. J. Phys. Chem. 1992, 96, 3847. (9) Lazaridis, T.; Paulaitis, M. E. J. Phys. Chem. 1993, 97, 5789. (10) Lazaridis, T.; Paulaitis, M. E. J. Phys. Chem. 1994, 98, 635.

Entropy of Hydrophobic Hydration (11) Paulaitis, M. E.; Ashbaugh, H. S.; Garde, S. Biophys. Chem. 1994, 51, 349. (12) Baranyai, A.; Evans, D. Phys. ReV. A 1989, 40, 3817. (13) Morita, T.; Hiroke, K. Prog. Theor. Phys. 1961, 25, 537. (14) De Dominicus, C. J. Math. Phys. 1962, 3, 983. (15) Bildstein, B.; Kahl, G. J. Chem. Phys. 1994, 100, 5882. (16) Baranyai, A.; Evans, D. Phys. ReV. A 1990, 42, 849. (17) Hummer, G.; Soumpasis, D. J. Chem. Phys. 1993, 98, 581. (18) Wallace, D. J. Chem. Phys. 1987, 87, 2282. (19) Mountain, R.; Raveche´, H. J. Chem. Phys. 1971, 53, 2250. (20) Laird, B.; Haymet, A. Phys. ReV. A 1992, 45, 5680. (21) Borzsa´k, I.; Baranyai, A. Chem. Phys. 1992, 165, 227. (22) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (23) Ben-Naim, A. J. Phys. Chem. 1978, 82, 792. (24) Yu, H.; Karplus, M. J. Chem. Phys. 1988, 89, 2366. (25) Garisto, F.; Kusalik, P. G.; Patey, G. N. J. Chem. Phys. 1983, 79, 6294. (26) Guillot, B.; Guissani, Y. J. Chem. Phys. 1993, 99, 8075. (27) Widom, B. J. Chem. Phys. B 1963, 39, 2808. (28) Matubayasi, N.; Reed, L.; Levy, R. J. Phys. Chem. 1994, 98, 10640. (29) Lee, B. Biophys. Chem. 1994, 51, 271. (30) Henderson, D. J. Chem. Phys. 1967, 46, 4306. (31) Ben-Naim, A. Water and Aqueous Solutions; Plenum Press: New York, 1974. (32) Mehrotra, P. K.; Beveridge, D. L. J. Am. Chem. Soc. 1980, 102, 4287. (33) Hermann, R. J. Phys. Chem. 1972, 76, 2754. (34) Pratt, L. R.; Chandler, D. J. Chem. Phys. 1977, 67, 3683. (35) Jorgensen, W. L.; Gao, J.; Ravimohan, C. J. J. Phys. Chem. 1985, 89, 3470.

J. Phys. Chem., Vol. 100, No. 5, 1996 1913 (36) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (37) Jorgensen, W. L. BOSS, version 2.8; Yale University: New Haven, CT, 1989. (38) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. Am. Chem. Soc. 1984, 106, 6638. (39) Jorgensen, W. L.; Chendrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. Jorgensen, W. L.; Madura, J. D. Mol. Phys. 1985, 56, 1381. (40) Owicki, J.; Scheraga, H. Chem. Phys. Lett. 1977, 47, 600. (41) Weeks, J. D.; Chandler, D.; Andersen, H. C. J. Chem. Phys. 1971, 54, 5237. (42) Stillinger, F. H. J. Solution Chem. 1973, 2, 141. (43) Pratt, L. R.; Pohrille, A. Proc. Natl. Acad. Sci. U.S.A. 1992, 89, 2995. (44) Ben-Naim, A. SolVation Thermodynamics; Plenum Press: New York, 1987. (45) This value was estimated on the basis of MC simulations in which 100 million configurations were used for sampling, as described in ref 11. A lower value for the infinite-dilution partial molar volume of methane in TIP4P water was reported in ref 11 on the basis only 30 million configurations. (46) Jorgensen, W. L.; Buckner, J. K. J. Phys. Chem. 1987, 91, 6083. (47) Tobias, D. J.; Brooks, C. L. J. Chem. Phys. 1990, 92, 2582. (48) Shing, K. S.; Gubbins, K. E. Mol. Phys. 1982, 46, 1109. (49) Jorgensen, W. L.; Blake, J. F.; Buckner, J. K. Chem. Phys. 1989, 129, 193. (50) Kirkwood, J. G.; Buff, F. P. J. Chem. Phys. 1951, 19, 774.

JP952387B

Entropy of Hydrophobic Hydration: Extension to Hydrophobic Chains

Recommend Documents