A thermodynamic analysis of solvation
HsiangAi Yu and Martin Karplus
Citation: J. Chem. Phys. 89, 2366 (1988); doi: 10.1063/1.455080
View online: http://dx.doi.org/10.1063/1.455080
View Table of Contents: http://jcp.aip.org/resource/1/JCPSA6/v89/i4
Published by the American Institute of Physics.
Additional information on J. Chem. Phys.
Journal Homepage: http://jcp.aip.org/
Journal Information: http://jcp.aip.org/about/about_the_journal
Top downloads: http://jcp.aip.org/features/most_downloaded
Information for Authors: http://jcp.aip.org/authors
Downloaded 03 May 2012 to 222.195.186.126. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
A thermodynamic analysis of solvation
HSiang-Ai Yu and Martin Karplus
Department o/Chemistry, Harvard University, Cambridge, Massachusetts 02138
(Received 20 January 1988; accepted 17 March 1988)
The free energy, energy, and entropy of solvation, relative to the pure liquid, are analyzed. Bya
coupling parameter integration it is shown that only averages over the solute-solvent
interaction energy contribute to the free energy and that the solvent-solvent interaction term,
which contributes the so-called cavity (solvent reorganization) term to the energy, is cancelled
exactly by a corresponding term in the entropy. These terms exist even in the infinite dilution
limit since they arise from the derivative of the free energy with respect to the solute density.
Following the approach of Garisto et 01. [J. Chem. Phys. 79, 6294 (1983)], the site-site
Ornstein-Zernike integral equations and HNC closures are used to determine the derivatives
of the distribution functions with respect to the density. This makes it possible to calculate the
energetic and entropic contributions to the solvation free energy in the infinite dilution limit.
The method is applied to pure solvent and to infinitely dilute aqueous solutions of cations,
anions and neutral Lennard-Jones particles. The results are in agreement with numerical
calculations of the thermodynamic quantities by use of finite difference values for the
temperature derivatives. A simple empirical relation for the charge dependence of the solvation
free energy is observed; it is shown for the case of an ion in a dipolar solvent, as typified by
aqueous electrolyte solutions, that the free energy of solvation varies quadratically with the
charge and is very nearly equal to one-halfthe solute-solvent portion of the solvation energy.
Some discussion of the relation of the present results to entropy-enthalpy compensation and to
computer simulations is given.
I. INTRODUCTION
The statistical thermodynamics of solution has evolved
considerably since the seminal developments of McMillian
and Mayerl and Kirkwood and Buff.2 Integral equation the-
ories, 3-11 and molecular dynamics, 12-18 or Monte Carlo sim-
ulations l 9-24 have made possible the calculation of thermo-
dynamic quantities for pure liquids and the solvation of
nonpolar and polar species. An essential issue in these stud-
ies is the contribution of the structural perturbations of the
solvent to the solution thermodynamics. To avoid complica-
tions from solute-solute interactions, most integral equation
studies have been made in the limit of infinite dilution; corre-
spondingly, simulations have tried to approach this limit by
using a single solute molecule in a box of solvent. It appears
intuitive that the change in solvent structure, relative to pure
solvent, should not playa role in the thermodynamic proper-
ties at infinite dilution. This is not correct, as has been point-
ed out by Garisto et 01. in their study of hard-sphere ions in
dipolar solvents.4 We use more general arguments to show
that the solvation energy and entropy do have contributions
from the change in solvent structure that are nonzero at infi-
nite dilution. However, these contributions cancel exactly in
the solvation free energy.
In Sec. II, we derive expressions for the solvation free
energy, energy, and entropy at finite concentration in terms
of canonical averages over the potentials associated with the
various interactions. The terms in the energy and entropy
that cancel in the free energy are isolated. These terms are
shown to arise from the change in the solvent structure even
in the infinite dilution limit. Focusing on the energy of solva-
tion, we show how to express the difference between the
average solute-solvent interaction energy < Uuv ) and the to-
tal solvation energy AEsol in terms of integrals over distribu-
tion functions. To calculate the magnitudes of the various
contributions, the formulation of the site-site Ornstein-Zer-
nike (SSOZ) integral equations and their closures3,s,2s-27 are
extended to determine the derivatives of the distribution
functions with respect to the density. Section III gives the
parameters for the specific solute and solvent models that are
used to illustrate Sec. II. Results of the calculations are pre-
sented and discussed· in Sec. IV. The conclusions are out-
lined in Sec. V. The Appendix shows explicitly that no effect
from solvent structural change contributes to the chemical
potential of solvation by use of the SSOZ equations with the
site-site HNC closure.
II. THEORY
In what follows, we consider a solution composed of a
single solute and a single solvent species for notational sim-
plicity; extensions of the formalism to more complex solu-
tions is straightforward. Both the solute and the solvent can
be flexible polyatomic molecules, although we finally spe-
cialize to simpler systems. Since we are concerned with ex-
cess quantities, all kinetic energy contributions cancel and
we can confine ourselves to the configurational contribution
to the thermodynamic properties.
A. Solvation thermodynamics at finite concentration
We employ a thermodynamic coupling parameter inte-
gration28 to express the excess Helmholtz free energy AAool
relative to pure solvent, for a system of Nv solvent molecules
and Nu solute molecules in the form
AA = t d)., (aUto" l ({R};A») . (1)
sol Jo aA..t
Here UtOtal ({i};A) is the total potential energy
2366 J. Chern. Phys. 89 (4),15 August 1988 0021-9606/88/162366·14$02.10 @) 1988 American Institute of PhysiCS
Downloaded 03 May 2012 to 222.195.186.126. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
H. Yu and M. Karplus: Thermodynamic analysis of solvation 2367
U total ( On;A. )
[
Nv • 1 Nv Nv • • ]
= i~1 Uv(lv) +"2 i~IJ~1 (1-6ij)Uvv (lv']v)
[
N.. 1 N. N. • •
+ i~IUu(lu;A.) +"2 i~lj~1 (l-6ij)Uuu (zu']u;A.)
N. Nv ]
+ i~lj~1 Uuv(iu,jv;A.)
= Uv,v({Rv}) + Uv,u({R};A.) (2)
for the solution in a fixed total volume of V. The quantity
U A (i A ;A.) represents the intramolecular potential for the ith
molecule of type A with configurational coordinates i A; for
rigid polyatomics, U A (i A;A.) will be singular, but this causes
no difficulties as the term can be eliminated. The quantity
UAB (iA,jB;A.) is the intertnolecuiar potential between the
ith molecule of type A and the j th molecule of type B. The
symbol on without a SUbscript represents all of the configu-
rational variables and the subscripted symbol {Rv} stands
for the solvent variables alone; the bar over a variable indi-
cates that it is a vector. The letter A. designates the coupling
parameter; A. = 0 corresponds to the pure solvent reference
state and A. = 1 to the solution under consideration; we have
written the equations in a general form that is not restricted
to a linear A.-coupling scheme. The bracket ( ) A denotes a
canonical configurational average at coupling strength A.;
that is, for the quantity P( on,A. ), we have
- Se -PUtotol({R};A) P( {R"};A.)dOn (P({R};,i»A= (3) S e - PU,otol({R};A) d{R}
with P = 1/ k BT, where k B is the Boltzmann constant.
To calculate the solvation energy and entropy per solute
molecule (aEool INu and TliSool INu)' we use the thermody-
namic relationships
aEsol = _1_ (ap!Msol ) ; TliSool = L ( a!Msol )
Nu Nu ap N,V Nu Nu ap N,V
(4)
together with Eqs. (1 )-( 3) to obtain
aEsol 1 [{-} i l d {- } ) ]
-- = - ( Uv,u ( R ;,i = 1» A = I + dA. d/l. (Uv,v ( Rv ) A
Nu Nu 0
1 - P i l [ - ( a -}) = - (Uv,u ({R};A. = 1» A= I + - dA. (Uv,v ({Rv}» A a
1
Uv,u ({R ;,i)
~ ~o ~ A
- a - ]
- (Uv,v ({Rv}) a,i Uv,u ({R};A) ) A
(5)
and
TliSool P t [ - ( a {-) ( {- a {-})
--;v:- = Nu Jo d,i (Uv,v ({Rv}» A aA. Uv,u ( R};A.) A - Uv,v ( Rv}) a,i Uv,u ( R ;A.) A
+ :u f d,i [ (Uv,u ( {R"};A) ) A (! Uv,u ( {R"};A.) t -( Uv,u ( {R"};A.) ! Uv,u ( {R"};,i) ) J . (6)
Thus, the excess free energy of solvation per solute molecule has the form
!Mool 1 1 - P i l [( {-} a {-})
--= - (aESOI - TliSsol ) = - (Uv,u ( {R};A. = 1) A = I + - dA. Uv,u ( R ;,i) --;-;- Uv,u ( R ;A) Nu Nu Nu Nu 0 aA A
- (Uv,u({R"};A.»A (! Uv,u({R"};A»)J. (7)
These results are exact and valid at finite solute concen-
trations, as well as at infinite dilution. In the usual discus-
sions of the solvation energy [Eq. (5)], the first term is de-
scribed as the average solute-solvent interaction energy and
the second as the cavity formation or solvent reorganization
term.8,9 As we show below, both terms contribute to the sol-
vation energy even at infinite dilution. However, whether
the system is at infinite dilution or at a finite concentration,
this cavity term does not contribute to the solvation free en-
ergy. It is cancelled exactly by an identical contribution from
the solvation entropy [Eq. (6)]. Moreover, averages over
the solvent-solvent potential energy do not contribute di-
rectly to the free energy when the pure solvent is chosen as
the reference state. Its contribution is indirect in that the
canonical averages [Eq. (3) ] have the total potential energy
in the exponential weighting factor.
The total free energy of solvation is the sum of an energy
contribution involving the canonical average of the solute-
solvent interaction potential [the first term in Eq. (7)] and
an entropic contribution involving more complex canonical
averages over the solute-solvent interaction potential. This
second term in Eq. (7), which is positive semidefinite (e.g.,
consider the case of a linear coupling scheme where it is in
the form offtuctuations of the Uv,u), represents an entropic
opposition to the solvation process. Thus, a negative free
energy of solvation will have to come entirely from a nega-
tive average solute-solvent interaction energy. In the case of
a hard-sphere solute the first term in Eq. (7) is zero and only
the second (entropic) term contributes to the free energy of
solvation,29 which is therefore always positive. We empha-
J. Chem. Phys., Vol. 89, No.4, 15 August 1988
Downloaded 03 May 2012 to 222.195.186.126. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
2368 H. Yu and M. Karplus: Thermodynamic analysis of solvation
size that our usage of the cavity formation term is restricted
to the solvent rearrangement defined by the change in the
solvent radial distribution functions [third term of Eq.
(9) ] .8,9 This contrasts with other usage in the literature
where the total free energy change upon introduction of a
repulsive solute is referred to as the cavity term. 29,30
To explore the infinite dilution limit and provide expres-
sions for calculations, we recast the expression for the solva-
tion energy in terms of distribution functions. We write Eq.
(5) in the integrated form
+ (Uv,v({1iV}»A=1 - (Uv,v({RV}»A=O]'
(8)
Explicitly evaluating each of the terms for a translation-
ally invariant system and specializing to a monatomic solute
and solvent so as to avoid cumbersome notation, we obtain
!::..Eso1 = Pv f d'7 Uuv Cr,A = 1 )guv (r,A. = 1)
Nu
+ P; f d'7 UuuCr,A = 1 )guu (r,A = 1)
+ p~ f ar Uvv ('7) 2pu
X [gvv(r,A = 1) -gvv(r,A = 0)]. (9)
Here Pv and Pu are the solvent and solute number density,
respectively, and the quantities g('7) are radial distribution
functions. A corresponding expression for the free energy
can be obtained; it is rather complicated and, since we do not
make use of it, we do not give it here. In general, evaluation
of the free energy requires an integration over the coupling
parameter A.; for the SSOZ equation with the HNC closure at
infinite dilution, the integration can be performed analyti-
cally to obtain a closed form expression in terms of distribu-
tion functions31 (see Sec. II B). For the energy itself, the A.
integration can be performed more generally as shown in Eq.
(9).
To take the infinite dilution limit of Eq. (9), we must
know the nature ofthe dependence on Pu for gvv (r,A = I).
Following H0ye and Stelln and Garisto et al.,4 we assume a
simple power series in Pu with the first term being linear in
Pu; that is, we write
gvv (r,pvpu,A. = 1) = gvv (r,pv,O,A = 1)
+ Pu~hvv (r,pv,O,A. = 1) + 8(p';- 1),
(10)
where ~hvv (=~gvv) is the correction to gvv ('7) that is first
order in Pu and 8(p';- 1) represents higher order terms; the
functional dependence of the distribution functions on the
solute and solvent densities is emphasized by including Pu
and Pv explicitly in the arguments. The total solvation ener-
gy per solute molecule at infinite dilution !::..E~?;o1 is obtained
from Eqs. (9) and (10) in the form
= Pv f ar Uuv (r,A. = 1 )g~~) (7;A. = 1)
+ P; f ar Uvv ('7)~h ~~) (r,A. = 1) (11a)
with the superscript (0) denoting the limit of infinite dilu-
tion. The first term on the right-hand side ofEq. (1Ia) gives
the average solute-solvent interaction energy per solute mol-
ecule
( ) (0) _ l' I ( ) Euv - 1m - Uuv
Pu-O Nu
(12)
in the presence of the unperturbed solvent at the limit of
infinite dilution; this corresponds to the first term of Eq. (5)
in this limit. The second term is the correction to the solva-
tion energy due to the perturbation of the solvent-solvent
correlation functions; this is the second term of Eq. (5) in
the infinite dilution limit and we denote this term as !::..E~~~ity .
It contributes even at this limit because the last term in Eq.
(9) varies inversely withpu, cancelling out thepu depend-
ence in Eq. (10). The essential point is that one is consider-
ing a partial molecular quantity (i.e., solvation energy per
solute molecule) and this represents a derivative with re-
spect to the solute number density to which there is a contri-
bution even at infinite dilution.
For the solvation free energy by contrast, the infinite
dilution result has the form31
Thus, the excess Helmholtz free energy per solute molecule
for an infinite dilute solution involves only g~~) ('7,A.). This
quantity is calculated directly from g~~) ('7,A.), the pure sol-
vent correlation function which is independent of A. [see Eq.
(20a)]. Consequently, in contrast to !::..E~?~1'~~?~1 does not
have contributions from changes in the solvent distribution
function. For the SSOZ equations with the HNC closure we
demonstrate this point explicitly in the Appendix.
In computer simulations of solvated species, it is com-
mon to treat a single solute molecule in a periodic box of
solvent molecules. The quantities that are calculated are the
average solute-solvent interaction energy and the change in
the solvent-solvent energy due to the presence of the solute.
If the simulation conditions corresponded to the infinite di-
lution limit (typical simulations consider I solute in 255 wa-
ter molecules, which leads to a concentration of about 0.2
M), the first quantity would be (Euv) (0) and the second
!::..E~~~ity = [!::..E~?~I - (Euv)(O)] . Although (Euv)(O) has a rel-
atively fast convergence, it is well-known that the second
quantity is difficult to calculate precisely by simulation
methods due to the large statistical fluctuations in the sol-
vent-solvent energy. An alternative to this approach, which
provides more precise results, is to determine !::...u~?~1 by a
free energy perturbation method28 and to evaluate!::..E~?~1
J. Chem. Phys., Vol. 89, No.4, 15 August 1988
Downloaded 03 May 2012 to 222.195.186.126. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions
H. Yu and M. Karplus: Thermodynamic analysis of solvation 2369
from the temperature derivative of ~~'!'I obtained by nu-
merical finite ditferencea;. i.e.,
aE(O) = ~(O) _ T(a~~'!'I) . (llb)
",001 "'001 aT N. v
This approach has been used in simulations33 and can also be
employed with the integral equation formulation.8•9•34 The
fact that ~~'!'I is independent of averages of U •• or guo (r,ii.)
is one source of the high precision in simulations for the free
energy of solvation (or more precisely differences in the free
energy of solvation of similar solutes) achieved by thermo-
dynamic perturbation methods.35
B. Site-slte Ornstein-Zernike equation
In this section we employ the renormalized SSOZ equa-
tion to obtain the change in the solvent-solvent distribution
function required for evaluating the solvation energy at infi-
nite dilution; most of the theory used here has been described
in the literature.5.25 The formalism is restricted to rigid
polyatomic solvent and solute molecules with site-site pair-
wise decomposable central force potentials. We briefly out-
line the main equations and emphasize the new features
needed to calculate 13h •• in Eq. (10). The renormalized
SSOZ integral equation for a solution composed of n. sol-
vent and n" solute interaction sites is given by5
Alii. A "","" A A
A=h - Q = [w + Qp]c*[w + ph], (14)
where p, w, h, c*, Q are n. + n" by n. + n .. dimensional site-
labeled matrices and the caret denotes functions in Fourier k
space. The symbol p represents the matrix of site densities
and w the block-diagonal intramolecular correlation matrix
with elements given as
A k sin (kLay ) (15) Way ( ) = ,
kLay
where Lay is the distance between sites a and r within a
given molecule; h = II - 1 is the matrix of intermolecular
total correlation functions with 1 the matrix where all ele-
ments are I; and c* is a modified direct correlation functions
defined by
A '" '" c* = c - + = c + (l(V - V*) (16)
'" with V, V* the site-labeled matrices of total and short-
ranged interaction potentials, respectively . We deal with po-
tentials which have a simple Coulomb interaction + as their
only long-range part; the quantity Q is the renormalized
Coulomb interaction defined, in Rossky and Oale's36 chain
sum notation, as
'" A '" Q = Ch[+lpw]. (17)
To solve Eq. (14), we adopt a site-site hypernetted chain
like closure in real space,
hay = exp( -/3U:y + hay - c:y ) - 1. (18)
Equations (14) and (18) can be solved iteratively for c· and
ii.5 The hypernetted chain-like closure is chosen because an
analytic expression is available for the excess chemical p0-
tential at infinite dilution aJl~?'!'1 .8.31 It has the form
A,,(O)
""""'tIOI
=~ f ar{Hh~~)(r,A.= 1)]2-c~~)(r,A.= 1)
- ~h ~~)(r,A. = 1)c~~)(r,A. = I)}. (19)
Some limitations of the site-site HNC closure have been dis-
cussed recently by Chandler et 01.37
To proceed to extract a first-order correction to the sol-
vent distribution functions within the spirit of Eq. (10), we
make use of the fact that Eq. (14) is equivalent to a set of
three coupled matrix equations when the elements are or-
dered by sites and species are explicitly identified as parts of
the solvent or solute; i.e., there are the solvent-solvent
equations (vv), the solute-solvent equations (uv), and the
solute-solute equations (uu). Expanding these sets of equa-
tions [Eq. (14)] around the limit p" = 0 and retaining
terms to first order inp", we obtain
A~~) = [w. + Q~~)P. ]c:.(O) [w. + p.ii~~)], (20a)
AA(O) = [-'* cA*(O) + Q"'(O)p c*(O)] [w + p ii(O)