Papers in Physics, vol. 2, art. 020001 (2010)

Received: 13 July 2009, Revised: 29 December 2009, Accepted: 7 February 2010
Edited by: S. A. Cannas
Reviewed by: P. Netz, Inst. de Qúımica, Univ. Federal do Rio Grande do Sul, Brazil
Licence: Creative Commons Attribution 3.0
DOI: 10.4279/PIP.020001

www.papersinphysics.org

ISSN 1852-4249

Structural and dynamic properties of SPC/E water

M. G. Campo,1∗

I have investigated the structural and dynamic properties of water by performing a series
of molecular dynamic simulations in the range of temperatures from 213 K to 360 K, us-
ing the Simple Point Charge-Extended (SPC/E) model. I performed isobaric-isothermal
simulations (1 bar) of 1185 water molecules using the GROMACS package. I quantified
the structural properties using the oxygen-oxygen radial distribution functions, order pa-
rameters, and the hydrogen bond distribution functions, whereas, to analyze the dynamic
properties I studied the behavior of the history-dependent bond correlation functions and
the non-Gaussian parameter α2(t) of the mean square displacement of water molecules.
When the temperature decreases, the translational (τ ) and orientational (Q) order param-
eters are linearly correlated, and both increase indicating an increasing structural order
in the systems. The probability of occurrence of four hydrogen bonds and Q both have
a reciprocal dependence with T , though the analysis of the hydrogen bond distributions
permits to describe the changes in the dynamics and structure of water more reliably.
Thus, an increase on the caging effect and the occurrence of long-time hydrogen bonds
occur below ∼ 293 K, in the range of temperatures in which predominates a four hydrogen
bond structure in the system.

I. Introduction

Water is the subject of numerous studies due to its
biological significance and its universal presence [1–
3]. The thermodynamic behavior of water presents
important differences compared with those of the
other substances, and many of the characteristics of
such behavior are often attributed to the existence
of hydrogen bonds between water molecules. Scien-
tists have found that the water structure produced
by the hydrogen bonds is peculiar as compared to
that of other liquids. Then, the advances in the
knowledge of hydrogen bond behavior are crucial
to understanding water properties.

∗E-mail: mario@exactas.unlpam.edu.ar

1 Universidad Nacional de La Pampa, Facultad de Ciencias
Exactas y Naturales, Uruguay 151, 6300 Santa Rosa, La
Pampa, Argentina.

The method of molecular dynamics (MD) allows
to analyze the structure and dynamics of water at
the microscopic level and hence to complement ex-
perimental techniques in which these properties can
be interpreted only in a qualitative way (infra-red
absorption and Raman scattering [4], depolarized
light scattering [5, 6], neutron scattering [7], fem-
tosecond spectroscopy [8–11] and other techniques
[12–14].

Among the usual methods to study the short
range order in MD simulations of water are the
calculus of radial distribution functions, hydrogen
bond distributions and order parameters. The ori-
entational order parameter Q measures the ten-
dency of the system to adopt a tetrahedral con-
figuration considering the water oxygen atom as
vertices of a tetrahedron, whereas the translational
order parameter τ quantifies the deviation of the
pair correlation function from the uniform value of

020001-1



Papers in Physics, vol. 2, art. 020001 (2010) / M. G. Campo

unity seen in an ideal gas [15, 16]. The order pa-
rameters are used to construct an order map, in
which different states of a system are mapped onto
a plane τ -Q. The order parameters are, in general,
independent, but they are linearly correlated in the
region in which the water behaves anomalously [17].

The dynamics of water can by characterized by
the bond lifetime, τHB, associated to the process of
rupturing and forming of hydrogen bonds between
water molecules which occurs at very short time
scale [9, 18, 19, 21–23]. τHB is obtained in MD
using the history-dependent bond correlation func-
tion P (t), which represents the probability that an
hydrogen bond formed at time t = 0 remained con-
tinuously unbroken and breaks at time t [24, 25].

Also, the dynamics of water can be studied by
analyzing the mean-square displacement time se-
ries M (t). In addition to the diffusion coefficient
calculation at long times in which M (t) ∝ t, in the
supercooled region of temperatures and at interme-
diate times M (t) ∝ tα (0 < α < 1). This behavior
of M (t) is associated to the subdiffusive movement
of the water molecules, caused by the caging effect
in which a water molecule is temporarily trapped
by its neighbors and then moves in short bursts due
to nearby cooperative motion. A time t∗ character-
izes this caging effect (see Sec. II for more details)
[26, 27].

In a previous work, we found a q-exponential be-
havior in P (t), in which q increases with T −1 ap-
proximately below 300 K. q(T ) is also correlated
with the probability of occurrence of four hydrogen
bonds, and the subdiffusive motion of the water
molecules [28].

The relationship between dynamics and struc-
tural properties of water has not been clearly es-
tablished to date. In this paper, I explore whether
the effect that temperature has on the water dy-
namics reflects a more general connection between
the structure and the dynamics of this substance.

II. Theory and method

I have performed molecular dynamic simulations of
SPC/E water model using the GROMACS package
[29, 30], simulating fourteen similar systems of 1185
molecules at 1 bar of pressure in a range of temper-
atures from 213 K to 360 K. I initialized the system
at 360 K using an aleatory configuration of water

Table 1: Details of the simulation procedure. Du-
ration of the stabilization period (test) and the MD
sampling (tM D) in the different ranges of tempera-
tures

Temp. range (K) test (ns) tM D (ns)
213 - 243 20.0 10.0
253 - 273 16.0 10.0
283 - 360 16.0 8.0

molecules, assigning velocities to the molecules ac-
cording to a Boltzmann’s distribution at this tem-
perature. For stabilization, I applied Berendsen’s
thermal and hydrostatic baths at the same tem-
perature and 1 bar of pressure [31]. Then, I ran
an additional MD obtaining an isobaric-isothermal
ensemble. I obtained the other systems in a similar
procedure, but using as initial configuration that
of the system of the preceding higher temperature
and cooling it at the slow rate of 30 K ns−1 [17].
Stabilization and sampling periods for the systems
at different temperatures are indicated in Table 1.
Simulation and sampling time steps were 2 fs and
10 fs, respectively. The sampling time step was
shorter than the typical time during which a hydro-
gen bond can be destroyed by libration movements.

I calculated the hydrogen bond distribution func-
tions f (n) (n = 0, 1, ..., 5), which is the probabil-
ity of occurrence of n hydrogen bonds by molecule,
considering a geometric definition of hydrogen bond
[20]. As parameters for this calculation, I used a
maximum distance between oxygen atoms of 3.5 Å
and a minimum angle between the atoms Odonor–
H–Oacceptor of 145◦.

The radial distribution function (RDF) is a stan-
dard tool used in experiments, theories, and sim-
ulations to characterize the structure of condensed
matter. Using RDFs, I obtained the average num-
ber, N , of water molecules in the first hydration
layer (the hydration number)

N = 4πρ
∫ rmin

0

g(r)r2dr (1)

where ρ is the number density.
The translational order parameter, τ , is defined

in Ref. [16] as

020001-2



Papers in Physics, vol. 2, art. 020001 (2010) / M. G. Campo

Figure 1: Hydrogen bonds distribution functions
f (n) (n = 0, ..., 5) versus T . The zones A, B and C
correspond to ranges of temperatures in which oc-
cur different relationships between f (4), f (3) and
f (2). Note the reciprocal scale for the tempera-
tures. See the text for details.

Figure 2: Oxygen-oxygen radial distribution func-
tions for the systems at 213 K (continuous line),
293 K (dashed line), and 360 K (dotted line). In-
set: The hydration number N vs. T 5.

τ ≡
∫ Sc

0

|g(s) − 1|ds (2)

where the dimensionless variable s ≡ rn1/3 is the
radial distance r scaled by the mean intermolecular

distance n1/3, and Sc corresponds to half of the
simulation box size.

The orientational order parameter Q is defined
as [15]

Q =

〈
1 −

3
8

N∑
i=1

4∑
j=1

4∑
l=j+1

[
cosθjik +

1
3

]2〉
(3)

where θjik is the angle formed by the atoms Oj –Oi–
Ok. Here, Oi is the reference oxygen atom, and Oj
and Ok are two of its four nearest neighbors. Q=1
in an ideal configuration in which the oxygen atoms
would be located in the vertices of a tetrahedron.

I obtained the bond correlation function P (t)
from the simulations by building a histogram of
the hydrogen bonds lifetimes for each configura-
tion. Then, I fitted this function with a Tsallis
distribution of the form

expq(t) = [1 + (1 − q) t]
1/(1−q) (4)

being t the hydrogen bond lifetime and q the nonex-
tensivity parameter [28, 32]. If q = 1, Eq. (4)
reduces to an exponential, whereas if q > 1, P (t)
decays more slowly than an exponential. This last
behavior occurs when long lasting hydrogen bonds
increase their frequency of occurrence.

The subdiffusive movement of water occurs when
the displacement of the molecules obeys a non-
Gaussian statistics. This behavior is characterized
by t∗, the time in which the non-Gaussian param-
eter α2(t) reaches a maximum [see Eq. (5)]. Then,
t∗ is the parameter associated to the average time
during which a water molecule is trapped by its en-
vironment (caging effect), and this prevents it from
reaching the diffusive state [26, 27].

α2(t) =
3〈r4(t)〉
5〈r2(t)〉

− 1 (5)

III. Results and discussion

Three zones or ranges of temperatures can be dis-
tinguished in the graph of the hydrogen bond dis-
tributions f (n) vs. T (see Fig. 1). Zone A (T >
350 K) in which f (3) > f (2) > f (4), zone B (293
K > T > 350 K) in which f (3) > f (4) > f (2), and
zone C (T < 293 K) in which f (4) > f (3) > f (2).
These results indicate a predominant structure of
three and two hydrogen bonds (3HB-2HB) in zone

020001-3



Papers in Physics, vol. 2, art. 020001 (2010) / M. G. Campo

A, 3HB-4HB in zone B, and 4HB-3HB in zone C,
respectively. f (4) ∝ T −1 in all ranges of tempera-
tures, showing that the tetrahedral structure of wa-
ter decreases with the increase of this variable.f (3)
increases with T up to 293 K, and then remains
approximately constant (∼ 0.4) up to 360 K. f (2)
also increases with T in all range of temperatures,
but only overcomes f (4) at T > 350 K.

Figure 3: Position of the first minimum of the
oxygen-oxygen radial distribution function vs T 4,
associated to the size of the first hydration layer.

Fig. 2 shows the oxygen-oxygen RDFs corre-
sponding to the systems at 213 K, 293 K and 360
K. When the temperature decreases, the minimum
and maximum tend to be more defined. This being
associated with an increasing order in the system.
The position of the first minimum moves closer to
the origin decreasing the size of the first hydration
layer (∝ T 4 see Fig. 3). Both facts can be associ-
ated to the decrease of the hydration number from
N ∼ 5 to N ∼ 4 (see inset, Fig. 2).

The simultaneous behavior of Q and τ is shown
in the order map of Fig. 4, in which the location of
the values corresponding to 293 K are indicated by
an arrow. The order parameters present similar be-
haviors with the temperature. Upon cooling, these
parameters are linearly correlated and move in the
order map along a line of increasing values, up to
reaching maximum values at 213 K. The slope of
the line increases a little for T > 293 K, indicating
that τ has a response to the increase of T slightly
higher than Q in this range of temperatures. The

Figure 4: Order map with the values of the or-
der parameters corresponding to the simulated sys-
tems. Note the change in the slope of the line at
T ∼ 273 K

positive values of the slopes indicate an increasing
order of the system when the temperature decrease.

The f (n) functions allow to obtain a more
detailed picture of the structural orientational
changes at shorter ranges between water molecules
than the orientational order parameter. While a
small change at 293 K occurs in the order map, the
structures of two, three and four hydrogen bonds
are alternated in importance when the temperature
changes. The ability of f (n) to more reliably de-
scribe the structure of the water occurs because the
calculation of the hydrogen bond distributions in-
cludes the location of the hydrogen atoms, whereas
Q only quantifies the changes in the average an-
gle between neighbor oxygen atoms. Although the
behavior of f (4) and Q are correlated (see Fig.
5), f (4) shows a greater response to the temper-
ature than Q, indicating that the main change in
the tetrahedral structure with the decrease of the
temperature occurs mainly in the orientation of
the bonds between water molecules. The approx-
imately linear correlation between both variables
also indicates a similar dependence with the tem-
perature (∝ T −1).

Figure 6 shows the behavior of the dynamical pa-
rameters t∗ and q with Q. The characteristic time
t∗ has an exponential response to Q ≥ 0.58 (T ≤
360 K), but the slope of the semilog plot of t∗ vs.
Q increases significantly for Q ≥ 0.67. A similar

020001-4



Papers in Physics, vol. 2, art. 020001 (2010) / M. G. Campo

Figure 5: Q vs f (4). The change in f (4) is higher
than that of Q in the range of temperatures stud-
ied. See the text for details.

change occurs for Q ≥ 0.67 in the linear correla-
tion between q and Q. Then, the values Q ≈ 0.67
and τ ≈ 1.1 of the order map can be associated
to changes in the dynamics of the system. The
transition of q ≈ 1 to q > 1 indicates the increase
of the probability of two water molecules remain-
ing bonded by a hydrogen bond during an unusual
long time, whereas the increase of t∗ is associated to
the increase of the time during which the molecules
remain in a subdiffusive regime.

However, only the analysis of the f (n) functions
reveals the structural modification that explains
the structural and dynamic changes in the system.
The changes in the increase of the order map, t∗(Q)
and q(Q) occur below 293 K, in the range of tem-
peratures in which prevail a structure of four hy-
drogen bonds in the system.

IV. Conclusions

The molecular dynamic method allows to study the
structure and dynamics of the SPC/E model of wa-
ter in the range of temperatures from 213 K to 360
K.

Lowering the temperature of the system from
360 K to 213 K, the number of water molecules
in the first hydration layer decreases from N ∼ 5
to N ∼ 4, along with a decrease in size. The in-
crease of the tetrahedral structure of the system is
also characterized by a growth of the percentage

Figure 6: (a) Semilog plot of t∗ vs. Q. (b) q vs. Q.
See the text for details.

of occurrence of four hydrogen bonds and the ori-
entational order parameter Q. However, only the
analysis of the behavior of the hydrogen bond dis-
tribution allows to deduce that, when a tetrahedral
structure associated to the percentage of four hy-
drogen bonds predominates, the behavior of the dy-
namical variables P (t) and t∗ show the occurrence
of long lasting hydrogen bonds and caging effect
between the molecules of the system.

Acknowledgements - I am grateful for the finan-
cial support by PICTO UNLPAM 2005 30807 and
Facultad de Ciencias Exactas y Naturales (UNL-
Pam).

[1] D Eisenberg, W Kauzmann, The structure and
properties of water, Clarendon Press, Oxford
(1969).

020001-5



Papers in Physics, vol. 2, art. 020001 (2010) / M. G. Campo

[2] F H Stillinger, Water revisited, Science 209,
451 (1980).

[3] O Mishima, H E Stanley, The relationship be-
tween liquid, supercooled and glassy mater, Na-
ture 396, 329 (1998).

[4] E W Castner, Y J Chang, Y C Chu, G E Wal-
rafen, The intermolecular dynamics of liquid
water, J. Chem. Phys. 102, 653 (1995).

[5] C J Montrose, J A Bucaro, J Marshall-
Coakley, T A Litovitz, Depolarized Rayleigh
scattering and hydrogen bonding in liquid wa-
ter, J. Chem. Phys. 60, 5025 (1974).

[6] W Danninger, G Zundel, Intense depolarized
Rayleigh scattering in Raman spectra of acids
caused by large proton polarizabilities of hydro-
gen bonds, J. Chem. Phys. 74, 2779 (1981).

[7] J Yeixeira, S H Chen, M-C Bellissent-Funel,
Molecular dynamics of liquid water probed
by neutron scattering, J. Mol. Liq., 48, 111
(1991).

[8] R Laenen, C Rauscher, A Laubereau, Dy-
namics of local substructures in water observed
by ultrafast infrared hole burning, Phys. Rev.
Lett. 80, 2622 (1998).

[9] S Woutersen, U Emmerichs, H J Bakker, Fem-
tosecond mid-IR pump-probe spectroscopy of
liquid water: Evidence for a two-component
structure, Science 278, 658 (1997).

[10] H K Nienhuys, S Woutersen, R A van Santen,
H J Bakker, Mechanism for vibrational relax-
ation in water investigated by femtosecond in-
frared spectroscopy, J. Chem. Phys. 111, 1494,
(1999).

[11] G M Gale, G Gallot, F Hache, N Lascoux, S
Bratos, J C Leickman, Femtosecond dynamics
of hydrogen bonds in liquid water a real time
study, Phys. Rev. Lett. 82, 1068, (1999).

[12] A H Narten, M D Danford, H A Levy, X-ray
diffraction study of liquid water in the temper-
ature range 4 �– 200 �, Faraday Discuss. 43,
97 (1967).

[13] A K Soper, F Bruni, M A Ricci, Site-site pair
correlation functions of water from 25 to 400
�: Revised analysis of new and old diffraction
data, J. Chem. Phys. 106, 247 (1997).

[14] K Modig, B G Pfrommer, B Halle,
Temperature-dependent hydrogen-bond ge-
ometry in liquid water, Phys. Rev. Lett. 90,
075502 (2003).

[15] H Tanaka, Simple physical explanation of the
unusual thermodynamic behavior of liquid wa-
ter, Phys. Rev. Lett. 80, 5750 (1998).

[16] J R Errington, P G Debenedetti, Relationship
between structural order and the anomalies of
liquid water, Nature 409, 318 (2001).

[17] N Giovambattista, P G Debenedetti, F
Sciortino, H E Stanley, Structural order in
glassy water, Phys. Rev. E 71, 061505 (2005).

[18] C A Angell, V Rodgers, Near infrared spec-
tra and the disrupted network model of nor-
mal and supercooled water, J. Chem. Phys. 80,
6245 (1984).

[19] J D Cruzan, L B Braly, K Liu, M G Brown, J
G Loeser, R J Saykally, Quantifying hydrogen
bond coopertivity in water: VRT spectroscopy
of the water tetramer, Science 271, 59 (1996).

[20] D. C. Rapaport, Hydrogen bonds in water,
Mol. Phys. 50, 1151 (1983).

[21] A Luzar, Resolving the hydrogen bond dynam-
ics conundrum, J. Chem. Phys. 113, 10663
(2000).

[22] F Mallamace, M Broccio, C Corsaro, A
Faraone, U Wandrlingh, L Liu, C Mou,
S H Chen, The fragile-to-strong dynamics
crossover transition in confined water: nuclear
magnetic resonance results, J. Chem. Phys.
124, 124 (2006).

[23] C J Montrose, J A Búcaro, J Marshall-
Coakley, T A Litovitz, Depolarizated light-
scattering and hydrogen bonding in liquid wa-
ter, J. Chem. Phys. 60, 5025 (1974).

[24] F W Starr, J K Nielsen, H E Stanley, Fast
and slow dynamics of hydrogen bonds in liquid
water, Phys. Rev. Lett. 82, 2294 (1999).

020001-6



Papers in Physics, vol. 2, art. 020001 (2010) / M. G. Campo

[25] F W Starr, J K Nielsen, H E Stanley,
Hydrogen-bond dynamics for the extended sim-
ple point charge model of water, Phys. Rev. E.
62, 579 (2000).

[26] S Chatterjee, P G Debenedetti, F H Stillinger,
R M Lynden-bell, A computational investiga-
tion of thermodynamics, structure, dynamics
and solvation behavior in modified water mod-
els, J. Chem. Phys. 128, 124511 (2008).

[27] M G Mazza, N Giovambattista, H E Stanley,
F W Starr, Connection of translational and
rotational dynamical heterogeneities with the
breakdown of the Stokes-Einstein and Stokes-
Einstein-Debye relations in water, Phys. Rev.
E 76, 031203 (2007).

[28] M G Campo, G L Ferri, G B Roston, q-
exponential distribution in time correlation
function of water hydrogen bonds, Braz. J.
Phys. 39, 439 (2009).

[29] H J C Berendsen, D van der Spoel, R
V Drunen, GROMACS: a message passing
parallel molecular dynamics implementation,
Comp. Phys. Comm. 91, 43 (1995).

[30] H J C Berendsen, J R Grigera, T P Straatsma,
The missing term in effective pair potentials,
J. Phys. Chem. 91, 6269 (1987).

[31] H J C Berendsen, J Postma, W van Gunsteren,
A Di Nola, J Haak, Molecular dynamics with
coupling to an external bath, J. Chem. Phys.
81, 3684 (1984).

[32] C Tsallis, Possible generalization of
Boltzmann-Gibbs statistics, J. Stat. Phys.
52, 479 (1988).

[33] L A Báez, P Clancy, Existence of a density
maximum in extended simple point charge wa-
ter, J. Chem. Phys. 101, 9837 (1994).

020001-7