Hydration free energy of hard-sphere solute over a wide range of size studied by various types of solution theories

The hydration free energy of hard-sphere solute is evaluated over a wide range of size using the method of energy representation, information-theoretic approach, reference interaction site model, and scaled-particle theory. The former three are distribution function theories and the hydration free energy is formulated to reflect the solution structure through distribution functions. The presence of the volume-dependent term is pointed out for the distribution function theories, and the asymptotic behavior in the limit of large solute size is identified. It is indicated that the volume-dependent term is a key to the improvement of distribution function theories toward the application to large molecules.


Introduction
Hydrophobicity plays a central role in the formation of self-organizing structures such as protein, micelle, and membrane in aqueous solution.The key idea of hydrophobicity is that water tends to reduce its contact with nonpolar solutes and promotes their association.Among various components of solute-solvent (water) interactions, the repulsive part is the origin of hydrophobicity.It is then insightful to focus on the effect of the repulsive interaction in hydration thermodynamics since one of the purposes of solution theory is to clarify how molecular interactions manifest themselves in the observable.
To isolate the effect of the repulsive interaction, the hard core particle is an ideal target of investigation.Its interaction with the solvent is purely repulsive and infinitely hard at the contact distance.The solvation (hydration) of a hard particle is a process to exclude the solvent (water) molecules from the region which the hard particle occupies.In fact, the concept of exclusion is not restricted to a molecular object.On the macroscopic or hydrodynamic scale, the solvent exclusion is a common process which leads, for example, to buoyancy.Given that the exclusion effect is operative for both macroscopic and molecular objects, it is of interest to examine the dependence on size.When a macroscopic object is of concern, its geometrical properties such as volume and surface area describe the effect of solvent exclusion; since a macroscopic system has a sharp boundary in its length scale, the volume and surface area can be defined without ambiguity.As for a molecule, in contrast, some ambiguity is always present in the definitions of volume and surface area since the "boundary" of a molecule is vague.From the viewpoint of solvent exclusion over a wide range of size, the hard particle is a useful target of examination.Even on the molecular level, its boundary is sharp and the geometrical properties such as volume and surface area can be defined without ambiguity.
In the present work, we focus on the hydration of hard-sphere particle.We evaluate the free energy of hydration over a wide range of size, and discuss the volume dependence to connect macroscopic and microscopic behaviors.The hydration free energy is calculated using 4 approaches: (1) method of energy representation, [1][2][3][4] (2) information-theoretic approach, [5][6][7] (3) reference interaction site model (RISM), [8][9][10][11][12] and (4) scaled-particle theory [13,14].( 1), (2), and (3) are molecular approaches in the sense that the structural information of the solution needs to be incorporated, while (4) is not a molecular scheme in that sense.Although (1) and ( 3) are applicable to general types of solute-solvent interactions, their usage is restricted to hard-sphere solute in the present work.In their formulations, on the other hand, (2) and ( 4) are specialized to the treatment of hard particle.As shown in section 2, (1) and ( 2) adopt the structural information of pure water from molecular simulation, and (3) and ( 4) are processed relying only upon potential functions or geometrical parameters.

Theory and methodology
In this section, we briefly review the approaches employed in the present work.The system treated is the aqueous solution of hard-sphere solute at infinite dilution.A single solute molecule is dissolved in water.
The method of energy representation is a theory of distribution functions in solution with emphasis on the free energy of solvation (hydration) [1][2][3].In the energy representation, the solutesolvent interaction energy is adopted as the coordinate of the distribution functions, and a functional for the solvation free energy is constructed in terms of energy distribution functions in the solution and reference solvent (pure water) systems of interest.Let ρ be the bulk density of the pure solvent and v be the volume from which the solvent is excluded.When the solute is a hard particle, the currently used version of functional [2] for the solvation free energy ∆µ reduces to a simple form as derived in Appendix A and is given by where k B is the Boltzmann constant, T is the temperature, and n is the occupancy number of solvent molecules in the region which the hard solute is to occupy after insertion.Ω is defined by equation (1) from the first and second cumulants of n, where the average • • • is taken in the pure solvent.When the system is the solution and the hard solute is present at full coupling of the solute-solvent interaction, n is always zero.When the system is the pure solvent and the region of insertion is specified by treating the solute as a test particle, n is variable.n is the average number of solvent molecules in the region of insertion and is equal to ρv.The term n 2 − n 2 reflects the solvent-solvent correlation in the pure solvent.This correlation affects the solvation free energy ∆µ through Ω.As noted in the footnote 44 of [2], no input is necessary from the solution system when the solute is a hard particle.Molecular simulation of the pure solvent system is necessary only to obtain n 2 − n 2 .The factor written as corresponds to the number fluctuation of the solvent molecules contained in the region of exclusion.
In the limit of large size of the region, this factor approaches ρk B T χ T , where χ T is the isothermal compressibility of the pure solvent system.At liquid conditions, the factor is much smaller than unity.Actually, the functional described in [2] is a combined functional of the Percus-Yevick (PY)-and hypernetted-chain (HNC)-type approximations.Since Ω < 0 holds in the conditions studied in the present work, the PY-type form is adopted and the functional is simplified into equation ( 2) [15].The second approach examined in the present work is the information-theoretic approach [5][6][7].In this approach, a Gaussian-like approximation is adopted for the probability of the occupancy number n of solvent molecules contained within the region of solute insertion in the pure solvent system.In Hummer et al.'s version, [5] the probability p n is set to where λ 0 , λ 1 , and λ 2 are determined from a set of constraints expressed as In equation ( 5), • • • is the average in the pure solvent system, as is so for equation (1).n 2 is an input brought from molecular simulation of the pure solvent and carries the information of solvent-solvent correlation.The solvation free energy ∆µ of a hard particle is then given by The third approach employed in the present work is the reference interaction site model (RISM) [8][9][10][11][12].In this approach, the site-site form of Ornstein-Zernike equation is formulated to introduce the direct correlation function and a closure relationship is adopted to give a self-consistent integral equation for pair correlation functions.The site-site radial distribution functions are obtained from the integral equation and determine the solvation free energy.
The fourth approach is the scaled-particle theory.It provides an analytically convenient expression for the solvation free energy of hard-sphere solute [13,14].In this theory, the detailed interaction and structure of the solution are not taken into account explicitly and affect the solvation free energy of the hard-sphere solute only through the solvent density ρ, pressure P , and temperature T of the system.The solvent molecule is modeled as a hard sphere with a properly chosen diameter σ v .When the radius of solvent exclusion (cavity) is λ, the corresponding solute has a diameter of 2λ − σ v .Since we are concerned with cavities of realistic size, we restrict our attention to the solutes which involve positive diameters and satisfy λ σ v /2.The free energy ∆µ of solvation (cavity formation) is then expressed in the scaled-particle theory as [13,14] where y is the packing fraction of the solvent defined as and R is the ratio of the solute diameter to the solvent diameter given by When the solvent is water, σ v is actually an effective diameter of the solvent molecule.It is an adjustable parameter in equation (7), and a value of 2.8 Å is adopted in the present work.This value is the best-fit value obtained in a previous molecular simulation [16].Actually, an exact expression for ∆µ is available when the cavity radius λ is small enough or large enough; when λ is very large, ∆µ is simply the product of the pressure P and the cavity volume.Equation ( 7) was formulated to smoothly interpolate the limits of small λ and large λ [13,14].
In this paper, we express the size of the hard sphere by the radius λ of solvent exclusion.As noted above, it corresponds to the hard sphere diameter 2λ − σ v when the solvent diameter is set to σ v .We employ water as the solvent at 25 • C and 1.0 g/cm 3 , and identify the region of solvent exclusion in terms of the hard-sphere center and the position of the oxygen atom.In this case, the factor of equation ( 3) appearing in equations ( 1) and ( 5) is given by where g(r) is the oxygen-oxygen radial distribution function of pure water.In the present work, g(r) used to calculate ∆µ with the method of energy representation and the information-theoretic approach was obtained from a molecular dynamics simulation of pure water.The potential model employed was the SPC/E model [17], and the simulation was performed for 100 ps in the canonical ensemble at 25 • C and 1.0 g/cm 3 .The number of water molecules in the unit cell was 10000, and the edge length of the cubic cell was ∼67 Å.When we calculated equation (10) with g(r) as an input, we used the g(r) from the SPC/E simulation in the range of r < 30 Å.
In the RISM calculation, we actually employed the dielectrically consistent version referred to as the DRISM theory [12].The HNC approximation was adopted as the closure.The SPC/E model was used as the water model, and the temperature and density were the same as the above.A necessary modification in the DRISM calculation was that a repulsive core is incorporated at the H site.This is needed to avoid catastrophic overlap of the interaction sites [18,19].Accordingly, the H site was also treated as a hard sphere of diameter 0.4 Å when a hard-sphere solute is inserted into water.Further, a hard-sphere diameter of 2.8 Å was assigned to the O site.Thus, when the radius of solvent exclusion measured with respect to the hard-sphere center and the O site is λ, additional hard interaction of radius (λ − 1.2) Å is incorporated between the hard-sphere center and the H site in the DRISM calculation [20].The details of the DRISM methodology are described in [21,22].A supplementary discussion using the molecular Ornstein-Zernike approach is presented in Appendix B.

Results and discussion
In figure 1, we show the hydration free energy ∆µ as a function of the cavity radius λ.In the range of λ < ∼ 3 Å, the ∆µ values from the method of energy representation and the informationtheoretic approach are indistinguishable from those by the scaled-particle theory.In this range, since the scaled-particle theory provides essentially exact values, [16] the above two methods are shown to be superior in treating a hard particle of relatively small size.At λ > ∼ 10 Å, the information-theoretic approach was found not to give a converged result since the argument in the exponential function of equation ( 4) was too large and the calculation in double precision was terminated.
When the size of a molecule is large, ∆µ is commonly expressed to be proportional to the surface area [23][24][25].In figure 2, we show ∆µ/(4πλ 2 ) as a function of λ.It is seen that ∆µ calculated by the scaled-particle theory grows with the surface area 4πλ 2 when λ > ∼ 10 Å.In the other approaches, ∆µ/(4πλ 2 ) increases with λ over the whole range of λ examined in the present work.Actually, the increase is linear, and this shows the presence of a volume-dependent term in the molecular length scale.
In real systems, the dispersion attraction is operative between a hydrophobic solute and water.Since the dispersion interaction is short-ranged, its contribution to ∆µ is expected to be proportional to the surface area and to λ 2 when λ is large.In previous treatments, although the effect of dispersion attraction is incorporated as perturbation, [13,14,[24][25][26] it is not negligible in ∆µ.This shows that the λ 2 term is reduced from that given in figures 1 and 2 when the dispersion attraction is introduced in the calculation.The hydration free energy ∆µ as a function of the cavity radius λ calculated using the method of energy representation (denoted as ER in the figure), the information-theoretic approach (denoted as IT), the RISM-HNC approach (denoted as RISM), and the scaled-particle theory (denoted as SPT).
Figure 2. ∆µ/(4πλ 2 ), the hydration free energy scaled by the surface area of the cavity as a function of the cavity radius λ calculated using the method of energy representation (denoted as ER in the figure), the information-theoretic approach (denoted as IT), the RISM-HNC approach (denoted as RISM), and the scaled-particle theory (denoted as SPT).
In the limit of large λ, the exact ∆µ approaches P v, where P is the pressure, and v is the volume of the exclusion region and is equal to (4/3)πλ 3 .By construction, the scaled-particle theory gives the correct limit value.In the other approaches examined, in contrast, the large-λ behavior of ∆µ is not exact and is discussed in the following.
To determine the asymptotic value of ∆µ at large v in the method of energy representation and the information-theoretic approach, it is useful to note that when the region of solvent exclusion is large enough, the factor of equation ( 3) reduces to where χ T is the isothermal compressibility of the pure solvent system.The expression in the energy representation given by equation (2) then leads, by virtue of ρk B T χ T 1 at liquid conditions, to In the information-theoretic approach, it was shown that the asymptotic form is expressed as [6] ∆µ It should be noted that equations ( 12) and ( 13) are valid for a hard particle of any shape; the validity is not restricted to a hard-sphere solute.
In water at 25 • C and 1.0 g/cm 3 , with the χ T value of 4.8 • 10 −5 atm −1 for the SPC/E model (this value corresponds well to the experimental value of 4.6 • 10 −5 atm −1 ) [16].With the RISM approach, the inspection of figure 2 indicates that ∆µ/v is ∼10000 atm at large v, and indeed agrees with the analytical expression at large v.The analytical expression in the RISM-HNC approach is presented in Appendix B. The ∆µ/v calculated above are the asymptotic values of ∆µ/v at large enough v.The comparison of the values shows that the method of energy representation provides a better estimate of ∆µ for a large solute.The exact value of ∆µ/v is 1 atm, [27] however, and it is concluded that except for the scaled-particle theory, the theories of solutions examined in the present work overestimate the v dependent term.Actually, since the surface tension of water is ∼70 mN/m, the volume term is dominant only when λ > ∼ 1 µm in the exact treatment.The scaled-particle theory incorporates the correct limiting behavior of the P v term in its construction.According to the asymptotic analysis, in contrast, the other approaches may exhibit an evident volume dependence even in the nanometer-scale.This is caused by the fact that the large v behavior is determined by the isothermal compressibility χ T , rather than the pressure P .The solvation free energy needs to be nicely computed in the nanometer-scale when its quantitative analysis is necessary for self-organizing structures such as protein, micelle, and membrane.A key to improve the large v behavior of the solvation free energy in distribution function theories will be to accurately and smoothly incorporate the crossover of the volume-dependent term from the χ T -dominated form to the P -dominated one.

Appendix B
We supplement the integral equation treatment with the molecular Ornstein-Zernike (MOZ) theory, in which the pair correlation functions and the integral equation are explicitly dependent on angular coordinates as well as on the radial distance.The multipole model was adopted for the water molecule; it is a hard sphere of diameter 2.8 Å in which a point dipole and a point quadrupole of tetrahedral symmetry are embedded [28].The dipole moment was determined through the iterative procedure given in [28] to incorporate the polarization effect in the mean-field sense.The resulting value at 25 • C and 1.0 g/cm 3 was found to be enhanced by 40% from the value at isolation (in dilute gas); it is equal to 2.6 D and is in agreement with the commonly accepted value in liquid [29].The quadrupole moment was taken to be 2.6 • 10 −26 esu cm 2 , the gas-phase value.The HNC approximation was adopted as the closure relationship.In order to numerically solve the MOZ-HNC equation, the pair potential and the correlation function were expanded in terms of rotational invariants.The methodological details are presented in [30,31].
The results of the MOZ-HNC scheme combined with the multipole water model were found to be within ∼20% from those of the method of energy representation in the small radius range of λ < ∼ 5 Å.When λ is larger, ∆µ/(4πλ 2 ) was observed to involve an evident slope as a function of λ; the slope corresponds to the coefficient of the volume-dependent term.The asymptotic behavior of ∆µ for large λ can be obtained from the closed-form expression for ∆µ.In the MOZ-HNC approach, ∆µ is expressed as [32] where x collectively denotes the position and orientation of the solvent (water) relative to the solute (hard sphere), and g uv , h uv , and c uv are the pair correlation function, total correlation function, and direct correlation function, respectively, between solute and solvent.When the solute is at infinite dilution, the integrals of h uv and c uv over the whole configuration are given through the MOZ equation as where h vv is the total correlation function between solvent and solvent and its integral over the whole configuration gives the compressibility χ T [33].Equation ( 21) is then rewritten as Within the region of solvent exclusion, g uv = 0 and h uv = −1.We further assume that h uv and c uv decay on the molecular length scale for any λ.Equation ( 23) then leads to when the excluded volume v is large enough and ρk B T χ T 1.Note the correspondence of equation (24) with equation ( 13) [15].In the multipole water model employed with the MOZ-HNC equation, χ T was calculated to be 3.6 • 10 −4 atm −1 , and is larger than the experimental value by an order of magnitude.In turn, the χ T value leads to ∆µ/v ≈ 1000 atm for a large enough cavity.This shows that with respect to the large v behavior, the MOZ-HNC approach combined with the multiple water model agrees apparently with the scaled-particle theory better than the other approaches treated in section 3.
The asymptotic analysis in the integral equation with the HNC closure is also possible for the RISM approach.In the RISM-HNC approximation, ∆µ is expressed as [34] where r is the radial distance, and g ui,vj , h ui,vj , and c ui,vj are the radial distribution function, total correlation function, and direct correlation function, respectively, between the i-th site of the solute and the j-th site of the solvent.The second term can be written as [35][36][37] where ∆V is the excess partial molar volume of the solute.∆µ ≈ v/(2χ T ) then follows from the arguments similar to those deriving equation (24) from equation ( 21).In the RISM-HNC approach, χ T was calculated to be 5.7•10 −5 atm −1 .This leads to ∆µ/v ≈ 10000 atm for a large enough cavity.

Figure 1 .
Figure1.The hydration free energy ∆µ as a function of the cavity radius λ calculated using the method of energy representation (denoted as ER in the figure), the information-theoretic approach (denoted as IT), the RISM-HNC approach (denoted as RISM), and the scaled-particle theory (denoted as SPT).