Table Of ContentFaraday Discussiosn, volume 141, 251-276, (2009)
What ice can teach us about water interactions: a critical
comparison of the performance of different water models
9
C. Vega*, J. L. F. Abascal, M. M. Conde and J. L. Aragones
0
0
Departamento de Qu´ımica F´ısica,
2
n
Facultad de Ciencias Qu´ımicas,
a
J
Universidad Complutense de Madrid,
3
1
28040 Madrid, Spain ,
]
i *) [email protected]
c
s
- (Dated: January 13, 2009)
l
r
t
m Abstract
.
t
a The performance of several popular water models (TIP3P, TIP4P, TIP5P and TIP4P/2005) is
m
analysed. For that purpose the predictions for ten different properties of water are investigated,
-
d
n namely: 1. vapour-liquidequilibria(VLE)andcritical temperature; 2. surfacetension; 3. densities
o
c of the different solid structures of water (ices); 4. phase diagram; 5. melting point properties; 6.
[
1 maximum in density at room pressure and thermal coefficients α and κT; 7. structure of liquid
v
3 water andice; 8. equation ofstateathighpressures;9. diffusioncoefficient; 10. dielectric constant.
0
8 For each property, the performance of each model is analysed in detail with a critical discussion
1
.
1 of the possible reason of the success or failure of the model. A final judgement on the quality of
0
9 these models is provided. TIP4P/2005 provides the best description of almost all properties of the
0
: list, with the only exception of the dielectric constant. In the second position, TIP5P and TIP4P
v
i
X yield an overall similar performance, and the last place with the poorest description of the water
r
a
properties is provided by TIP3P. The ideas leading to the proposal and design of the TIP4P/2005
are also discussed in detail. TIP4P/2005 is probably close to the best description of water that
can be achieved with a non polarizable model described by a single Lennard-Jones (LJ) site and
three charges.
Keywords:
1
I. INTRODUCTION
Water is probably the most important molecule in our relation to nature. It forms the
matrix of life,1 it is the most common solvent for chemical processes, it plays a major
role in the determination of the climate on earth, and also it appears on planets, moons
and comets.2 Water is interesting not only from a practical point of view, but also from a
fundamental point of view. In the liquid phase water presents a number of anomalies when
compared to other liquids.3,4,5,6,7 In the solid phase it exhibits one of the most complex phase
diagrams, having fifteen different solid structures.3 Dueto its importance andits complexity,
understanding the properties of water from a molecular point of view is of considerable
interest. The experimental study of the phase diagram of water has spanned the entire
20th century, starting with the pioneering works of Tammann8 and Bridgman9 up to the
recent discovery of ices XII, XIII and XIV.10,11 The existence of several types of amorphous
phases at low temperatures,12,13,14 the possible existence of a liquid-liquid phase transition
in water15,16 , the properties of ice at a free surface17,18,19,20,21,22,23 and the interaction with
hydrophobic molecules24,25 have also been the focus of much interest in the last two decades.
Computer simulations of water started their road with the pioneering papers by Watts
and Barker26 and by Rahman and Stillinger27 about forty years ago. A key issue when
performing simulations of water is the choice of the potential model used to describe the
interaction between molecules.28,29,30,31,32 A number of different potential models have been
proposed. An excellent survey of the predictions of the different models proposed up to 2002
was made by Guillot.33. Probably the general feeling is that no water potential model is
totally satisfactory and that there are no significant differences between water models.
Intherecent years the increase inthecomputing power hasallowed the calculationofnew
properties which can be used as new target quantities in the fit of the potential parameters.
More importantly, some of these properties have revealed as a stringent test of the water
models. In particular, recently we have determined the phase diagram for different water
models and we have found that the performance of the models is quite different.34,35,36 As
a consequence, new potentials have been proposed. In this paper we want to perform a
detailed analysis of the performance of several popular water models including the recently
proposed TIP4P/2005, and including in the comparison some properties of the solid phases,
thus extending the scope of the comparison performed by Guillot33. The importance of
2
including solid phase properties in the test of water models was advocated by Morse and
Rice37 and Whalley38 , among others, and Monson suggested that the same should be true
for other substances39,40.
It must be recognised that water is flexible and polarizable. That should be taken into
account in the next generation of water models.41,42 However, it is of interest at this stage
to analyse how far it is possible to go in the description of real water with simple (rigid,
non-polarizable) models. For this reason we will restrict our study to rigid non-polarizable
models. Besides, in case a certain model performs better than others, it would be of interest
to understand why isthis so since this informationcan beuseful inthedevelopment offuture
models. As commented above, the study of the phase diagram of water can be extremely
useful to discriminate among different water models. Thus, inthe comparison between water
models, we incorporate not only properties of liquid water but also properties of the solid
phases of water. The scheme of this paper is as follows. In Section II, the potential models
used in this paper are described. In Section III, we present the properties that will be used
in the comparison between the different water models. Section IV gives some details of the
calculations of this work. Section V presents the results of this paper. Finally, in Section
VI, the main conclusions will be discussed.
II. WATER POTENTIAL MODELS
In this work we shall focus on rigid, non-polarizable models of water. All the models
we are considering locate the positive charge on the hydrogen atoms and a Lennard-Jones
interaction site on the position of the oxygen. What are the differences between them? They
do differ in three significant aspects:
•
The bond geometry. By bond geometry we mean the choice of the OH bond length
and H-O-H bond angle of the model.
•
The charge distribution. All the models place one positive charge at the hydrogens
but differ in the location of the negative charge(s).
•
The target properties. By target properties we mean the properties of real water
that were used to fit the potential parameters (forcing the model to reproduce the
experimental properties).
3
In this paper we compare the performance of the following potential models: TIP3P43,
TIP4P43, TIP5P44 and TIP4P/200536. This selection presents an advantage. All these
models exhibit the same bond geometry (the d distance and H-O-H angle in all these
OH
models are just those of the molecule in the gas phase, namely 0.9572 ˚A and 104.52 degrees,
respectively). Therefore, any difference in their performance is due to the difference in their
charge geometry and/or target properties. Let us now present each of these models.
A. TIP3P
TIP3P was proposed by Jorgensen et al.43. In the TIP3P model the negative charge is
located on the oxygen atom and the positive charge on the hydrogen atoms. The parameters
of the model (i.e, the values of the LJ potential σ and ǫ and the value of the charge on the
hydrogen atom) were obtained by reproducing the vaporization enthalpy and liquid density
of water at ambient conditions. TIP3P is the model used commonly in certain force fields
of biological molecules. A model similar in spirit to TIP3P is SPC45. We shall not discuss
in this paper SPC model so extensively, since we want to keep the discussion within TIPnP
geometry (by TIPnP geometry we mean that d adopts the value 0.9572 ˚A and the H-O-H
OH
angle is 104.52 degrees). In the SPC model the OH bond length is 1 ˚A and the H-O-H is
109.47 degrees (the tetrahedral value). The parameters of SPC were obtained in the same
way as those of TIP3P. Thus TIP3P and SPC are very much similar models. Nine distances
must be computed to evaluate the energy between two TIP3P (or SPC) water models so
that their computational cost is proportional to 9.
B. TIP4P
Thekeyfeatureofthismodelisthatthesitecarryingthenegativecharge(usuallydenoted
as the M site) is not located at the oxygen atom but on the H-O-H bisector at a distance
of 0.15 ˚A. This geometry was already proposed in 1933 by Bernal and Fowler.46. TIP4P
was proposed by Jorgensen et al.43 who determined the parameters of the potential to also
reproduce the vaporization enthalpy and liquid density of liquid water at room temperature.
ItisfairtosaythatTIP4P,althoughquitepopular,isprobablylessoftenusedthanTIP3Por
SPC/E. The reader may be surprised to learn that the reasons for that are the appearance of
4
a massless site (the M center) and the apparently higher computational cost. To deal with a
masslesssitewithinaMDpackageonemusteithersolvetheorientationalequationsofmotion
for a rigid body (for instance, using quaternions) or using constraints after distributing the
forceactingontheMcenter amongtherestoftheatomsofthesystem. These optionsshould
be incorporated in the MD package and this is not always the case. Also the computational
cost of TIP4P is proportional to 16 when no trick is made (due to the presence of four
interaction sites) but it can be reduced to 10 by realizing that one must only compute the
O-O distance to compute the LJ contribution, plus 9 distances to compute the Coulombic
interaction. Thus, forsomeusers TIP4P isnota goodchoice asawater modeleither because
is slower than TIP3P in their molecular dynamics package or because the package may have
problems in dealing with a massless site. These are methodological rather than scientifical
reasons but they appear quite often. Some modern codes as GROMACS47 (just to mention
one example) solve these two technical aspects quite efficiently.
C. TIP5P
TIP5P is a relatively recent potential: it was proposed in 2000 by Mahoney and
Jorgensen.44 This model is the modern version of the models used in the seventies (ST2)
in which the negative charge was located at the position of the lone pair electrons.48 Thus
instead of having one negative charge at the center M, this model has two negative charges
at the L centres. Concerning the target properties this model reproduces the vaporization
enthalpy and density of water at ambient conditions. This is in common with TIP3P and
TIP4P. However, TIP5P incorporates a new target property: the density maximum of liq-
uid water. The existence of a maximum in density at approximately 277 K is one of the
fingerprints properties of liquid water. Obtaining density maxima in constant temperature-
constant pressure (NpT) simulations of water models is not difficult from a methodological
point of view.31,49,50 However, very long runs are required (of the order of several millions
time steps (MD) or cycles (MC)) to determine accurately the location of the density max-
imum. Thus, it is not surprising that the very first model able to reproduce the density
maximum was proposed in 2000 when the computing power allowed an accurate calculation
of the temperature and density at the maximum. Since TIP5P consists of five interaction
sites, itapparentlyrequirestheevaluationof25distances. However, usingthetrickdescribed
5
for TIP4P, it is possible to show that one just requires the evaluation of 17 distances.
D. SPC/E
In this paper we shall focus mainly on TIPnP like models. However, it is worth to
introduce SPC/E,51 considered by many as one of the best water models. SPC/E has the
same bond geometry as SPC. Concerning the charge distribution, it locates the negative
charge at the oxygen atom, as SPC and TIP3P. The target properties of SPC/E were
the density and the vaporization enthalpy at room temperature. So far everything seems
identical to SPC. The key issue is that SPC/E only reproduces the vaporization enthalpy of
real water when a polarization energy correction is included. Berendsen et al.51 pointed out
that, when using non-polarizable models, one should include a polarization correction before
comparing the vaporizationenthalpy ofthe modelto theexperimental value. This isbecause
the dipole moment is enhanced in these type of models in order to account approximately
for the neglected many body polarization forces. Thus, when the vaporization enthalpy
is calculated one compares the energy of the liquid with that of a gas with the enhanced
dipole moment. It is then necessary to correct for the effect of the difference between the
dipole moment of the isolated molecule and that of the effective dipole moment used for the
condensed phase. The correction term is given by:
E (µ−µ )2
pol gas
= (1)
N 2α
p
whereα isthepolarizabilityofthewatermolecule, µ isthedipolemomentofthemolecule
p gas
in the gas phase and µ is the dipole moment of the model. SPC/E reproduces the vaporiza-
tion enthalpy of water only if the correction given by Eq.1 is used. The introduction of the
polarization correction is the essential feature of SPC/E.
E. TIP4P/2005
The TIP4P/2005 was recently proposed by Abacal and Vega36 after evaluating the phase
diagram for the TIP4P and SPC/E models of water.34,52 It was found that TIP4P provided
a much better description of the phase diagram of water than SPC/E. Both models yielded
rather low melting points.53 For this reason, it was clear that the TIP4P could be slightly
6
modified, while still predicting a correct phase diagram but improving the description of the
melting point. The TIP4P/2005 has the same bond geometry as the TIP4P family. Be-
sides TIP4P/2005 has the same charge distribution as TIP4P (although the distance from
the M site to the oxygen atom is slightly modified). The main difference between TIP4P
and TIP4P/2005 comes from the choice of the target properties used to fit the potential
parameters. TIP4P/2005 not only uses a larger number of target properties than any of the
water models previously proposed but also they represent a wider range of thermodynamic
states. As it has become traditional, the first of the target properties is the density of water
at room temperature and pressure (but in TIP4P/2005 we have also tried to account for
the densities of the solid phases). More significantly, TIP4P/2005 incorporated as target
properties those used in the succesfull SPC/E and TIP5P models. As in SPC/E, the po-
larization correction is included in the calculation of the vaporization enthalpy. Secondly
(as in TIP5P) the maximum in density of liquid water at room pressure (TMD) was used
as a target property. Finally, completely new target properties are the melting temperature
of the hexagonal ice (Ih) and a satisfactory description of the complex region of the phase
diagram involving different ice polymorphs.
In Table I the parameters of the models SPC, SPC/E, TIP4P, TIP5P and TIP4P/2005
are shown. The dipole moment and the three values of the traceless quadrupolar tensor are
given in Table II. There, we also present the values of the quadrupole moment Q which is
T
defined as Q = (Q −Q )/2. We have chosen the H-O-H bisector as the Z axis, the X
T xx yy
axis inthedirection ofthe linejoining the hydrogen atomsso that they axis isperpendicular
to the molecular plane. It can be shown that, for models consisting of just three charges,
the value of Q is independent of the origin (when the z axis is located along the H-O-H
T
bisector).54,55,56 An interest feature appears in Table II. Although for most of water models
the dipole moment is close to 2.3 Debyes, they differ significantly in their values of the
quadrupole moment. Notice that neither the dipole nor the quadrupole moment were used
as target properties by any of the water models discussed here. Thus the charge distribution
determines the aspect of the quadrupolar tensor. One may suspect that, since quadrupolar
forces induce a strong orientational dependence, the differences in the quadrupolar tensor
between different water models will be manifested significantly in the solid phases where the
relative orientation of the molecules is more or less fixed by the structure of the solid. The
couplingbetweendipolarandquadrupolarinteractionsiswellknownsincesometimeagodue
7
to early studies about the behavior of hard spheres with dipole and quadrupole moments.57
However, the role of the quadrupole in the properties of water has been probably overlooked
although there were some clear warnings about its importance.58,59,60
III. AN EXAM FOR WATER MODELS
Guillot presented in 200233 a study of the performance of water models to describe sev-
eral properties of water. Although this is only six years ago there are several reasons to
perform this study once again. At that time TIP5P was just released. TIP4P/2005 was
proposed three years later. In these years a more precise determination of some properties
of the water models (surface tension, temperature of maximum density) has been obtained.
More importantly, the calculation in the last years of water properties that were almost
completely unknown before (as, for instance, predictions for the solid phases of water, the
determination of the melting points and the phase diagram calculations) make interesting
a new comparative of the performance of the models. It is necessary to select some prop-
erties to establish the comparison. The selected properties should be as many as possible
but representative of the different research communities with interest in water. We will
not include in the comparison properties of the gas phase (second virial coefficients, vapour
densities) because only a polarizable model can be successful in describing all the phases of
water. Non-polarizable models can not describe simultaneously the vapour phase and the
condensed phases. Thus, the models described above fail in describing vapour properties
because they ignore the existence of the molecular polarizability.61,62 Therefore, here we will
focus only in the properties of the condensed phases (liquid and solids). The ten properties
of the test will be the following ones:
•
1. Vapour-liquid equilibria (VLE) and critical point.
•
2. Surface tension.
•
3. Densities of the different ice polymorphs.
•
4. Phase diagram calculations.
•
5. Melting temperature and properties at the melting point.
8
•
6. Maximum in density of water at room pressure (TMD). Values of the thermal
expansion coefficient, α, and the isothermal compressibility, κ .
T
•
7. Structure of water and ice Ih.
•
8. Equation of state at high pressures.
•
9. Self-diffusion coefficient.
•
10. Dielectric constant.
In order to give an assessment of the quality of the predictions we will assign a score
to each of properties depending on the predictions of the model. We recognise that any
score is rather arbitrary. We do not pretend to give an absolute test but rather to give a
qualitative idea of the relative performance of these water models. So we have devised a
simple scoring scheme: for each property we shall assign 0 points to the model with the
poorest performance, 1 to the second, 2 to the third and 3 points to the model showing the
best performance.
IV. CALCULATION DETAILS
In this work we compare the performance of TIP3P, TIP4P, TIP5P and TIP4P/2005. To
perform the comparison we shall use data taken from the literature, either from our previous
works or from other authors. For the cases where no results are available we have carried
out new simulations to determine them. In all the calculations of this work the LJ part
of the potential has been truncated and standard long range corrections are used. Unless
stated otherwise all calculations of this work were obtained by truncating the LJ part of the
potential at 8.5 ˚A. The importance of an adequate treatment of the long range Coulombic
forces when dealing with water simulations has been pointed out in recent studies.63,64,65,66
In fact, in the case of water, the simple truncation of the Coulombic part of the potential
causes a number of artifacts.63,64,65,66 In our view, for water, the simple truncation of the
Coulombic part of the potential should be avoided and one should use a technique treating
adequately the long range Coulombic interactions (as for instance Ewald sums or Reaction
Field67,68,69). The Ewald sums are especially adequate for phase diagram calculations since
they can be used not only for the fluid phase but also for the solid phase. Because of this,
9
the Ewald summation technique68 has been employed in this work for the calculation of
the long range electrostatic forces. NpT and phase diagram calculations were done with
our Monte Carlo code. Transport properties were determined by using GROMACS 3.3.47
Isotropic NpT simulations were used for the liquid phase (and cubic solids) while anisotropic
Monte Carlo simulations (Parrinello-Rahman like)70,71 were used for the solid phases.
Recently, we have computed the phase diagram for the TIP4P model by means of free
energy calculations. For the solid phases, the Einstein crystal methodology proposed by
Frenkel and Ladd was used.72 Further details about these free energy calculations can be
found elsewhere.73 For the fluid phase, the free energy was computed by switching off the
chargesofthewatermodeltoarrivetoaLennardJonesmodelforwhichthefreeenergyiswell
known.74 Thefreeenergy calculationsforthefluidandsolidphases leadtothedetermination
of a single coexistence point for each coexistence line. Starting at this coexistence point,
the complete coexistence lines were obtained by using Gibbs Duhem integration.75 The
Gibbs Duhem integration (first proposed by Kofke) is just the numerical integration of the
Clapeyron equation:
dp s −s h −h
II I II I
= = (2)
dT v −v T(v −v )
II I II I
where we use lower case for thermodynamic properties per particle, and the two coexistence
phases are labelled as I and II respectively. Since the difference in enthalpy and volume
between two phases can be easily determined in computer simulations, the equation can
be integrated numerically. Therefore, a combination of free energy calculations and Gibbs
Duhem integration allowed us to determine the phase diagram of the TIP4P model.
In this work we have also computed the phase diagram for the TIP3P and TIP5P models.
Instead of using the free energy route to obtain the initial coexistence point we have used
Hamiltonian Gibbs Duhem integration76 which we briefly summarise. Let us introduce a
coupling parameter,λ, which transforms a potential model A into a potential model B (by
changing λ from zero to one):
U(λ) = λU +(1−λ) U . (3)
B A
It is then possible to write a generalized Clapeyron equation as
dT T(< u −u >II − < u −u >I )
B A N,p,T,λ B A N,p,T,λ
= (4)
dλ h −h
II I
where u is the internal energy per molecule when the interaction between particles is
B
described by U (and a similar definition for u ). If a coexistence point is known for the
B A
10
