Structural and dynamic properties of SPC/E water

 

M. G. Campo,¹*

Email: mario@exactas.unlpam.edu.ar
Universidad Nacional de La Pampa, Facultad de Ciencias Exactas y Naturales, Uruguay 151, 6300 Santa Rosa, La Pampa, Argentina.

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, using the Simple Point ChargeExtended (SPC/E) model. I performed isobaricisothermal simulations (1 bar) of 1185 water molecules using the GROMACS package. I quantified the structural properties using the oxygenoxygen radial distribution functions, order parameters, and the hydrogen bond distribution functions, whereas, to analyze the dynamic properties I studied the behavior of the historydependent bond correlation functions and the nonGaussian parameter 0.2 (i) of the mean square displacement of water molecules. When the temperature decreases, the translational (r) and orientational (Q) order parameters 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 longtime 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. Scientists 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.

The method of molecular dynamics (MD) allows to analyze the structure and dynamics of water at the microscopic level and hence to complement experimental techniques in which these properties can be interpreted only in a qualitative way (infrared absorption and Raman scattering [4], depolarized light scattering [5, 6], neutron scattering [7], femtosecond spectroscopy [811] and other techniques [1214].

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 orientational order parameter Q measures the tendency of the system to adopt a tetrahedral configuration considering the water oxygen atom as vertices of a tetrahedron, whereas the translational order parameter T quantifies the deviation of the pair correlation function from the uniform value of unity seen in an ideal gas [15, 16]. The order parameters are used to construct an order map, in which different states of a system are mapped onto a plañe tQ. The order parameters are, in general, independent, but they are linearly correlated in the región in which the water behaves anomalously [17].

The dynamics of water can by characterized by the bond lifetime, t¡jb , associated to the process of rupturing and forming of hydrogen bonds between water molecules which occurs at very short time scale [9, 18, 19, 2123]. t¡jb is obtained in MD using the historydependent bond correlation function P(t), which represents the probability that an hydrogen bond formed at time í = 0 remained continuously unbroken and breaks at time t [24, 25].

Also, the dynamics of water can be studied by analyzing the meansquare displacement time series M(t). In addition to the diffusion coefñcient calculation at long times in which M(t) oc t, in the supercooled región of temperatures and at intermedíate times M(t) oc t^a (0 < a < 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* characterizes this caging effect (see Sec. II for more details) [26, 27].

In a previous work, we found a </exponential behavior in P(t), in which q increases with T^¹ approximately 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 structural properties of water has not been clearly established to date. In this paper, I explore whether the effect that temperature has on the water dynamics 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 temperatures from 213 K to 360 K. I initialized the system at 360 K using an aleatory conñguration of water molecules, assigning velocities to the molecules according to a Boltzmann's distribution at this temperature. For stabilization, I applied Berendsen's thermal and hydrostatic baths at the same temperature and 1 bar of pressure [31]. Then, I ran an additional MD obtaining an isobaricisothermal ensemble. I obtained the other systems in a similar procedure, but using as initial conñguration that of the system of the preceding higher temperature and cooling it at the slow rate of 30 K ns¹ [17].

Table 1: Details of the simulation procedure. Duration of the stabilization period (t_est) and the MD sampling (ímd) in the different ranges of temperature

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

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 hydrogen bond can be destroyed by libration movements. I calculated the hydrogen bond distribution functions f(n) (n = 0,1, ...,5), which is the probability of occurrence of n hydrogen bonds by molecule, considering a geometric deñnition of hydrogen bond [20]. As parameters for this calculation, I used a máximum distance between oxygen atoms of 3.5 A and a minimum angle between the atoms

tl\J'acceptor Oí 14o .

The radial distribution function (RDF) is a standard tool used in experiments, theories, and simulations to characterize the structure of condensed matter. Using RDFs, I obtained the average number, N, of water molecules in the ñrst hydration layer (the hydration number)

where p is the number density. The translational order parameter, t, is deñned in Ref. [16] as

 

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 occur different relationships between /(4),/(3) and /(2).

Note the reciprocal scale for the temperatures. See the text for details.
Figure 2: Oxygenoxygen radial distribution functions for the systems at 213 K (continuous line), 293 K (dashed line), and 360 K (dotted line). Inset: The hydration number JVvs. T⁵.

distance n¹'³, and S_c corresponds to half of the simulation box size.

The orientational order parameter Q is deñned as [15]

where 6jh¡ is the angle formed by the atoms 0¿0¿ Ofc. Here, 0¿ is the reference oxygen atom, and 0¿ and Ofc are two of its four nearest neighbors. Q=í in an ideal conñguration in which the oxygen atoms would be located in the vértices 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 conñguration. Then, I ñtted this function with a Tsallis distribution of the form

being t the hydrogen bond lifetime and q the nonextensivity 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 oceurs when long lasting hydrogen bonds increase their frequeney of oceurrence.

The subdiffusive movement of water oceurs when the displacement of the molecules obeys a nonGaussian statistics. This behavior is characterized by t*, the time in which the nonGaussian parameter a.2(t) reaches a máximum [see Eq. (5)]. Then, t* is the parameter associated to the average time during which a water molecule is trapped by its environment (caging effect), and this prevenís it from reaching the diffusive state [26, 27].

III. Results and discussion

where the dimensionless variable s = rn¹'³ is the radial distance r scaled by the mean intermolecular

Three zones or ranges of temperatures can be distinguished in the graph of the hydrogen bond distributions f(n) vs. T (see Fig. 1). Zone A (T > 350 K) in which /(3) > /(2) > /(4), zone B (293 K > T > 350 K) in which /(3) > /(4) > /(2), and zone C (T < 293 K) in which /(4) > /(3) > /(2). These results indicate a predominant structure of three and two hydrogen bonds (3HB2HB) in zone A, 3HB4HB in zone B, and 4HB3HB in zone C, respectively. /(4) oc T^¹ in all ranges of temperatures, showing that the tetrahedral structure of water decreases with the increase of this variable./(3) increases with T up to 293 K, and then remains approximately constant (~ 0.4) up to 360 K. /(2) also increases with T in all range of temperatures, but only overcomes /(4) at T > 350 K.

Figure 3: Position of the ñrst minimum of the oxygenoxygen radial distribution function vs T⁴, associated to the size of the ñrst hydration layer.

Fig. 2 shows the oxygenoxygen RDFs corresponding to the systems at 213 K, 293 K and 360 K. When the temperature decreases, the minimum and máximum tend to be more deñned. This being associated with an increasing order in the system. The position of the ñrst minimum moves closer to the origin decreasing the size of the ñrst hydration layer (oc T^A see Fig. 3). Both facts can be associated to the decrease of the hydration number from N ~ 5 to N ~ 4 (see inset, Fig. 2).

The simultaneous behavior of Q and t is shown in the order map of Fig. 4, in which the location of the valúes corresponding to 293 K are indicated by an arrow. The order parameters present similar behaviors with the temperature. Upon cooling, these parameters are linearly correlated and move in the order map along a line of increasing valúes, up to reaching máximum valúes at 213 K. The slope of the line increases a little for T > 293 K, indicating that t has a response to the increase of T slightly higher than Q in this range of temperatures. The positive valúes of the slopes indicate an increasing order of the system when the temperature decrease.

Figure 4: Order map with the valúes of the order parameters corresponding to the simulated systems. Note the change in the slope of the line at T ~ 273 K

 

The /(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 describe the structure of the water occurs because the calculation of the hydrogen bond distributions includes the location of the hydrogen atoms, whereas Q only quantiñes the changes in the average angle between neighbor oxygen atoms. Although the behavior of /(4) and Q are correlated (see Fig. 5), /(4) shows a greater response to the temperature 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 approximately linear correlation between both variables also indicates a similar dependence with the temperature (oc T^_1).

Figure 6 shows the behavior of the dynamical parameters 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 signiñcantly for Q > 0.67. A similar change occurs for Q > 0.67 in the linear correlation between q and Q. Trien, the valúes Q « 0.67 and t « 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 remaining 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.

 

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

However, only the analysis of the f(n) functions reveáis the structural modiñcation 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 temperatures in which prevail a structure of four hydrogen bonds in the system.

IV. Conclusions

The molecular dynamic method allows to study the structure and dynamics of the SPC/E model of water 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 ñrst hydration layer decreases from N ~ 5 to Af ~ 4, along with a decrease in size. The increase of the tetrahedral structure of the system is also characterized by a growth of the percentage of occurrence of four hydrogen bonds and the orientational order parameter Q. However, only the analysis of the behavior of the hydrogen bond distribution allows to deduce that, when a tetrahedral structure associated to the percentage of four hydrogen bonds predominates, the behavior of the dynamical variables P(t) and tshow the occurrence of long lasting hydrogen bonds and caging effect between the molecules of the system.

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

 

Acknowledgements I am grateful for the financial support by PICTO UNLPAM 2005 30807 and Facultad de Ciencias Exactas y Naturales (UNLPam).

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

[2] F H Stillinger, Water revísíted, Science 209, 451 (1980).         [ Links ]

[3] O Mishima, H E Stanley, The relatíonshíp between liquid, supercooled and glassy mater, Nature 396, 329 (1998).         [ Links ]

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

[5] C J Montrose, J A Búcaro, J MarshallCoakley, T A Litovitz, Depolarized Rayleigh scattering and hydrogen bonding in liquid water, J. Chem. Phys. 60, 5025 (1974).         [ Links ]

[6] W Danninger, G Zundel, Intense depolarized Rayleigh scattering in Raman spectra of acids caused by large protón polarizabilities of hydrogen bonds, J. Chem. Phys. 74, 2779 (1981).         [ Links ]

[7] J Yeixeira, S H Chen, MC BellissentFunel, Molecular dynamics of liquid water probed by neutrón scattering, J. Mol. Liq., 48, 111 (1991).         [ Links ]

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

[9] S Woutersen, U Emmerichs, H J Bakker, Ferntosecond midIR pumpprobe spectroscopy of liquid water: Evidence for a twocomponent structure, Science 278, 658 (1997).         [ Links ]

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

[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).         [ Links ]

[12] A H Narten, M D Danford, H A Levy, Xray difjraction study of liquid water in the temperature range 4 °C 200 °C, Faraday Discuss. 43, 97 (1967).         [ Links ]

[13] A K Soper, F Bruni, M A Ricci, Sitesite pair correlation functions of water from 25 to 400 °C: Revised analysis of new and oíd difjraction data, J. Chem. Phys. 106, 247 (1997).

[14] K Modig, B G Pfrommer, B Halle, Temperaturedependent hydrogenbond geometry in liquid water, Phys. Rev. Lett. 90, 075502 (2003).

[15] H Tanaka, Simple physical explanation of the unusual thermodynamic behavior of liquid water, 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 Angelí, V Rodgers, Near infrared spectra and the disrupted network model of normal 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 dynamics 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 fragiletostrong 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 MarshallCoakley, T A Litovitz, Depolarizated lightscattering and hydrogen bonding in liquid water, 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).

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

[26] S Chatterjee, P G Debenedetti, F H Stillinger, R M Lyndenbell, A computational investigation of thermodynamics, structure, dynamics and solvation behavior in modified water models, 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 StokesEinstein and StokesEinsteinDebye relations in water, Phys. Rev. E 76, 031203 (2007).

[28] M G Campo, G L Ferri, G B Roston, qexponential 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 BoltzmannGibbs statistics, J. Stat. Phys. 52, 479 (1988).

[33] L A B´aez, P Clancy, Existence of a density maximum in extended simple point charge water, J. Chem. Phys. 101, 9837 (1994).