Pure Appl. Chem., Vol. 73, No. 12, pp. 1–11, 2001. DRAFT
© 2001 IUPAC
1
Using simulation to study solvation in water*
R. M. Lynden-Bell
1,‡
, J. C. Rasaiah
2
, and J. P. Noworyta
2
1
Atomistic Simulation Group, School of Mathematics and Physics, Queen’s
University, Belfast BT7 1NN, UK;
2
Department of Chemistry, University of Maine,
Orono, ME 04469, USA
Abstract: Simulations of simple solutes (charged and uncharged spheres) in model water have
been performed in order to elucidate aspects of solvation in water at ambient and supercriti-
cal states. The variation of solvation entropy as a function of solute charge has been used to
investigate hydrophobic and hydrophilic ordering and the structure-making and structure-
breaking effects of ions. Simulations with model solvents, which differ from water in certain
features, have been used to try to identify the particular properties of water that are associat-
ed with these phenomena.
INTRODUCTION
This contribution to the 27
th
International Conference on Solution Chemistry describes some of our
recent work on solvation of ions and small solvents in water, which is a rather unusual molecular liquid.
We shall discuss the role that simulation can play in assisting our understanding of solvation at the molec-
ular level. Our particular interest in this paper is the solvation of small atomic ions and small spherical
uncharged, hydrophobic solutes in water both at room temperature and at supercritical temperatures.
Water at ambient and at biological conditions is a common and important solvent. As is well
known, it is a good solvent for ions and polar molecules (hydrophilic solutes), and an exceptionally poor
solvent for nonpolar or hydrophobic solutes [1,2]. The phenomenon of hydrophobicity has been dis-
cussed for many years before and there are many relevant experimental measurements [3]. In recent
years, simulation has been added to the range of tools used to try to understand this phenomenon. In
fact, there are a range of phenomena included in the term hydrophobic effect—the low solubility of
methane in water compared to methanol or organic liquids; the immiscibility of oil and water; the wet-
ting or non-wetting of surfaces by water; hydrophobic hydration. Here, we shall be concerned with the
solvation of small solutes in very dilute solutions.
It is generally understood that hydrophobicity is a property of the water rather than the solute. It
is not so much that methane is water-hating but that water is methane-hating [4]. The properties of liq-
uid water depend on the intermolecular potential between water molecules and in particular on the for-
mation of hydrogen bonds. In ice, each water molecule forms four hydrogen bonds (two donated and
two accepted), which results in a rather open network structure. In the liquid under ambient conditions,
this structure persists to some extent, at least at short range, and it is this that is thought to give water
its unique properties. The interaction of ions with water can lead to breaking of this structure, which is
replaced by a different ordering of the molecules in the local field of the ion. The balance between the
structure-breaking and structure-making tendencies can be seen in the entropy of solvation, which we
shall examine in simulations. It has also been proposed that hydrophobic solutes are intrinsically struc-
ture-making, that is, they induce local structure around themselves in what has been described as ice-
berg formation. Simulations have shown that this probably does not occur to any significant extent.
*Plenary lecture presented at the 27
th
International Conference on Solution Chemistry, Vaals, The Netherlands, 26–31 August
2001. Other presentations are published in this issue, pp. XXXX–XXXX.
‡
Corresponding author
Supercritical water is a fluid above the critical point. The density can be varied continuously over
a wide range of values from gas-like to liquid-like. It is commercially important as it can be used as a
solvent for reactions (e.g., oxidation of waste) and then removed by reducing the density. It also has the
great advantage of being nontoxic. Hydrogen bonds are less important at these high temperatures, and
the solvation properties are more
normal or organic-like. The fact that the density can be varied wide-
ly makes it an interesting medium for theoretical studies of solvation, as the effects of density and tem-
perature changes can be separated.
SIMULATION
There are a number of types of atomistic simulation available for studies. While quantum simulations
are more accurate, they are very expensive computationally and it is only realistically possible to do
small simulations (50 molecules or less) for short times (a few ps). This is quite inadequate for explor-
ing the range of configurations that are typical of a solution, and hence classical simulations are appro-
priate. These can be carried out for hundreds (or even thousands) of molecules for times of a nanosec-
ond or longer.
As the aim is to study microscopic aspects of solvation and the relationship between microscop-
ic and macroscopic properties, it is not necessary to have complete quantitative accuracy (the aim is not
to supersede experimental work, but to complement and understand it). However, the model must con-
tain the correct physics. Thus, the model intermolecular potential can be approximate, but must not be
too crude. One real advantage of simulations is that there is the possibility of changing the model (or
the state point) to inaccessible or unphysical values. In this work, for example, we use ions with charges
of a fraction of an electron in order to understand the solvation thermodynamics. One should note that
even in classical simulations there are severe limitations of time scale, distance scale, and accessible
concentration range.
Intermolecular potentials for water
The basic physics that must be included in the model are the molecular size, the molecular shape, and
the intermolecular interactions. These are usually described by an intermolecular site–site potential with
several sites on each molecule. The inter-site interactions include repulsion at short distances which
describes the shape of the molecule; attractive dispersion (van der Waals’) interactions at intermediate
range; and electrostatic interactions. The latter may be attractive or repulsive, are both short and long
range, and are particularly important in determining the local arrangement of molecules in the liquid.
There are other terms that contribute to intermolecular forces, such as partial covalency and charge-
transfer, but, provided that no chemical bonds exist between molecules, it is sufficient to include repul-
sion, dispersion, and electrostatic terms even for hydrogen-bonded liquids.
The work described here used the SPC/E model for water [5]. This is one of a number of models
with fixed charges on a few sites and a Lennard–Jones center on the oxygen atom. In this particular
model there are 3 sites (O and two H) with a bond length of 1 Å and a tetrahedral bond angle. As the
molecule is neutral and symmetric, the electrostatics is described by a single parameter (the charge on
the oxygen atom, say). Two parameters describe the Lennard–Jones interaction and two the geometry.
The molecular size is determined by the value of the Lennard–Jones sigma parameter, while the direc-
tionalities of the intermolecular interactions including the hydrogen bonds are determined by the par-
tial charges and the positions of the H sites relative to the oxygen. This model is remarkably successful
for modeling liquid water under ambient conditions and is reasonably successful under other condi-
tions. We may conclude that it contains the essential physics of the intermolecular potential, although
we would not expect it to be quantitatively accurate. One missing ingredient is the effect of polariz-
ability. In real liquids, the instantaneous electron distribution in each molecule responds to the instan-
R. M. LYNDEN-BELL et al.
© 2001 IUPAC, Pure and Applied Chemistry 73, 0000–0000
2
taneous environment, while in this model the electrostatic properties of each molecule are fixed. The
dipole moment at 2.35 D is higher than in the gas phase (1.85 D), representing the average effect of the
polarization due to the environment, but does not fluctuate. This large enhancement is supported by ab
initio calculations which, when analyzed in different ways, give values of 2.45 D [6] or 2.95 D [7] in
the liquid state, however, the same calculations show significant fluctuations.
While the use of this model is reasonable for ambient conditions, we have also used it in the
supercritical region where the effects of fluctuations may be more important and the average polariza-
tion less. There are polarizable models [8,9] that are used in this region, but the aim of our work is to
try to isolate effects of temperature and density on solvation using the simplest realistic model avail-
able. An important point is that the critical point and static dielectric properties of SPC/E water are not
too different from real water [10,11].
Our calculations [12–18] have been carried out using standard methods of molecular dynamics
[19], with some extensions in the runs with variable size or charge [12]. It is important to treat the long-
range part of the electrostatics in a consistent way. We have used both the reaction field method and the
Ewald summation method. The reader is referred to the original papers for more details of the simula-
tions.
SOLVATION OF IONS AND UNCHARGED SPHERES AT AMBIENT CONDITIONS
We have investigated the solvation of spherical solutes in very dilute solutions, that is, with one solute mol-
ecule per simulation box. The solvation thermodynamics is defined by considering a box of volume V con-
taining N water molecules (this is the periodically repeated unit in the simulation). The solvation energy,
free energy, entropy, etc. are defined as the difference in these quantities for the box containing the water
molecules and the solute molecule and two separate isolated boxes of volume V, one with the water mol-
ecules and one with the solute molecule. One can also consider this as the difference between the ther-
modynamic properties of the box with the solute molecule and the water molecules when they interact and
when they do not interact. As the volume is held constant, the natural thermodynamic variables are the
Helmholz free energy, A, and the internal energy U, rather than the Gibbs free energy and the enthalpy.
The solvation free energy is related to the Ostwald solubility L
ost
of a gas by A
solv
= –kTln L
ost
and is equal
to the chemical potential of the solute. It is possible to measure A
solv
and U
solv
from simulations and
hence to determine the solvation entropy S
solv
from their difference. It should be noted that S
solv
is the
partial derivative of the chemical potential with respect to temperature at constant volume and differs
from the partial molar entropy of the solute, which is the partial derivative of the chemical potential with
respect to temperature at constant pressure. The difference is negligible at ambient conditions, but may
be significant in the supercritical region.
The solutes considered are simple Lennard–Jones spheres with a single charge. The solvation
thermodynamics depends on just three parameters, the oxygen charge, q, the size σ, and the well depth
ε of the Lennard–Jones interaction between the solute molecule and the water oxygen. We have kept
the Lennard–Jones well depth constant and equal to the value [20] recommended for a number of atom-
ic ions and determined the free energy as a function of charge and size. This is done by a sequence of
simulations:
• For zero charge and small values of solute size, the solvation free energy was found by the Widom
particle insertion method. This gives an absolute value for the free energy.
• For larger sizes, the variation of free energy with size was determined by thermodynamic inte-
gration of the average derivative of the energy with respect to the value of σ. This was done in an
extended ensemble simulation in which the size was allowed to vary.
• With a fixed value of σ (that for a Na
+
ion) the variation of the solvation free energy with charge
was found by similar methods.
• This was repeated for other values of the solute size.
© 2001 IUPAC, Pure and Applied Chemistry 73, 0000–0000
Using simulation to study solvation in water
3
Figure 1 shows the variation of the solvation free energies with charge for solutes of two differ-
ent sizes. The first point to note is that the chemical potential decreases (the solute is more stabilized
by solvation) as the magnitude of the charge increases. This is primarily due to electrostatic interactions.
and is what would be expected from the Born model. It can be partly explained by the lowering of the
free energy in solution compared to the gas phase by the polarization of the medium by the charge. This
is a particularly large effect due to the high dielectric constant of water. However there are features that
are peculiar to the molecular structure of water, in particular the considerable asymmetry between pos-
itive and negative ions. Negative ions act as hydrogen bond acceptors and so have a greater and asym-
metrical interaction with the nearest water molecules than the positive ions do. In the latter case, the
water molecules align with their dipole moments pointing (on average) along the radial directions from
the ion. A molecular description of water is necessary to describe this asymmetry as continuum mod-
els, such as the Born model, predict a symmetric quadratic dependence of the free energy on the charge.
Solvation entropies
Figure 2 shows the solvation entropies for solutes the size of a sodium and a chloride ion, respectively.
Further examples are given in our paper [12]. While the solvation free energy depends on both short-
and long-range interactions, the solvation entropy is a more local effect, and is more dependent on the
molecular nature of the solvent. The unusual double maximum in these curves is an interesting feature
of the curves in Fig. 2. At high values of the magnitude of the charge, the solvation entropy becomes
more negative as the solvent molecules become aligned in the electric field of the ion. One anticipates
this latter “structure-making effect” in any solvent formed from polar molecules. The double maximum,
on the other hand, is less common. The increase in solvation entropy with charge at small magnitudes
of the charge is the result of structure-breaking by the ion, which depends, of course, on the solvent hav-
ing structure to be broken. The double maximum occurs when the structure-making balances the struc-
ture-breaking. This is an old concept, which can be demonstrated unambiguously in these simulations
as the charge can be varied continuously and independently of other parameters in simulations.
The extent and pattern of structure-breaking depends on the size of the solute. For larger solutes,
the depth of the minimum and the separation between the two maxima increases, showing that the struc-
ture-making tendencies take longer to swamp the structure-breaking tendencies. We note that again
there is considerable asymmetry between positive and negative ions of the same size. The positive ion
Na
+
is a net structure maker in the sense that q = 1 lies well outside the two double maxima. Structure-
breaking and structure-making are about equal for an ion the size of Cs
+
while larger ions such as Rb
+
R. M. LYNDEN-BELL et al.
© 2001 IUPAC, Pure and Applied Chemistry 73, 0000–0000
4
Fig. 1 Variation of the solvation free energies (chemical potentials) of spheres the size of sodium and chloride ions
with ion charge at room temperature.
lie on the inside of the maximum and so can be described as structure-breakers. The negative ions are
all net structure-makers according to this definition with I
–
lying just beyond the maximum for a nega-
tively charged solute. However, there are other possible definitions of structure-breaking and structure-
making, and it has been more common to base these terms on the absolute values of the solvation
entropies.
Indeed another way of thinking about these curves is to concentrate on the difference between the
observations and those expected for a simple polar solvent. From this point of view, it is the region of
low charge that is anomalous and we can describe the abnormally low solvation entropy of uncharged
solutes to hydrophobic ordering. The low entropy associated with the solution of hydrophobic solutes
has also been known for many years from experiments and was ascribed to a local freezing of the water
induced by the solute, sometimes known as the iceberg effect. However, the word “freezing” implies a
decrease in enthalpy as well as entropy, which is not found. One can now ask the question “What is spe-
cial about water which gives rise to all the phenomena described as hydrophobic?” A supplementary
question is whether this double maximum can be used as a signature of the hydrophobic effect.
From the measurements of the solvation thermodynamics of uncharged spheres as a function of
size (see operations 1 and 2 above), we find that the energy of solvation is essentially zero and the free
energy of solvation is determined by the entropy. The fact that the energy is very small shows that the
entropy of solvation is determined by the probability of finding a cavity in the liquid that is large enough
to insert the sphere. Thus, it depends both on the free volume in the liquid and its distribution. It is not
surprising that the free energy of solvation is positive and increases with solute size, or alternatively that
the entropy of solvation is negative and becomes more negative as the solute size increases.
One important factor that contributes to the large negative entropy of solvation is simply that the
water molecules are smaller than the molecules of most organic liquids. For the same fraction of free
volume, the size of the voids scales with the size of the solvent molecules and much of the abnormally
low entropy of solvation of hydrophobic solutes can simply be attributed to the small size of water mol-
ecules. This is only indirectly a property of the water molecule in that, if it did not have hydrogen bonds,
it would not be a liquid under ambient conditions. However, this explanation is insufficient to account
for the more subtle effects of hydrophobicity such as the double maximum in the entropy versus charge
curves. The existence of a network structure might be expected to be important in this context. It has
been shown [22] that the distribution of cavities in the neat liquid is sharpened by the network struc-
ture, lowering the solvation entropy for larger solutes.
Another contribution to the solvation entropy of a nonpolar solute is from the formation of a
hydrogen-bonded cage enclosing the solute. This is confirmed in computer simulations by visual obser-
© 2001 IUPAC, Pure and Applied Chemistry 73, 0000–0000
Using simulation to study solvation in water
5
Fig. 2 Variation of the solvation entropies of spheres the size of sodium and chloride ions with ion charge at room
temperature.
