q-exponential distribution in time correlation function of water hydrogen bonds

q-exponential distribution in time correlation function of water hydrogen bonds

M. G. Campo; G. L. Ferri; G. B. Roston

Departamento de Física, Universidad Nacional de la Pampa Uruguay 151,(6300) Santa Rosa, La Pampa, Argentina

ABSTRACT

In a series of molecular dynamics simulations we analyzed structural and dynamics properties of water at different temperatures (213 K to 360 K), using the Simple Point Charge-Extended (SPC/E) water. We detected a q-exponential behavior in the history-dependent bond correlation function of hydrogen bonds. We found that q increases with T ^-1 below approximately 300 K and is correlated to the increase of the tetrahedral structure of water and the subdiffusive motion of the molecules.

Keywords: Water hydrogen bonds, q-statistics, Cage effect

1. INTRODUCTION

Water physical properties have been investigated in far greater details, because water is one of the most ubiquitous substances on Earth and most life processes are related with their properties [1-3]. In particular, structural and dynamic properties of the hydrogen bonds have been investigated using different experimental and theoretical techniques, nevertheless several issues remain unsolved. X-Ray diffraction, neutron diffraction and proton magnetic shielding tensor measurements can be related to the oxygen-oxygen radial distribution, hydrogen bonds geometry of water and other structural properties [4-6].

A number of experimental and simulation techniques were used to study the elemental dynamics process of rupturing and forming of hydrogen bonds in water, characterized by an average bond lifetime τ_HB and caused by diffusion and librational motions of water molecules on a very fast time scale [7-13]. In molecular dynamics (MD), a "history-dependent" bond correlation function P(t) is used to obtain τ_HB [14, 15]. P(t) represents the probability that an hydrogen bond formed at a time t = 0 remained continuously unbroken and will be broken at time t.

Both MD simulations and experimental data for depolarizated light scattering show an Arrhenious behavior of τ_HB, where τ_HB∝ exp(-E* /RT ) and E* is an activation energy [13-15]. E* can be interpreted as the energy required to break a hydrogen bond ( 10 KJ/mol). However, the functional behavior of P(t) has not been clearly established to date. In previous works it was found that P(t) not to verify either an exponential or a power law behavior in the liquid and supercooled region of water, regardless of the hydrogen bond definition used [14, 15].

In other works the dynamics of liquid and supercooled water was studied by analyzing the behavior of the mean-square displacement time series (M(t) )

where r(t) is the position of de oxygen atoms of water at time t, r(0) is their initial position and the brackets () denote ensemble average. At room temperature the motion of water molecules at short times is ballistic, where M(t) ∝ t^α and α 2 for pure ballistic motion. At later times the water molecules move diffusively (α 1 according to Einstein relation) and the diffusion coefficient D can be calculated. However, a plateau emerges in M(t) (0 < α < 1) at intermediate times in the supercooled region of water (T 273 K) where the movement of the water molecules is called subdiffusive [16, 17]. This phenomenon is attributed to a caging behavior of water, in which a water molecule is temporarily trapped by its neighbors and then moves in short burst due to nearby cooperative motion.

The diffusion is associated with a non-Gaussian statistics, when the Einstein relation is not satisfied [18-20]. In particular, the caging behavior of water is studied by Mazza et. al [17] calculating the non-Gaussian parameter α₂(Δt) of the displacement of water molecules.

Constantino Tsallis and collaborators introduced the q-exponential probability distribution [21]. This can be defined through their "complementary" distribution functions, also called "survival " functions:

Tsallis et al. proposed these distributions to handle statistical-mechanical systems with long-range interactions, necessitating a non-extensive generalization of the ordinary Gibbs-Shannon entropy. Following Jaynes's procedure [22] of maximizing an entropy subject to constraints on expectation values, they got the q-exponential ditributions, in wich κ enforces the constrains, and q measures the departure from extensivity, Boltzmann-Gibbs statistics being recovered as q → 1. Sometimes, a q-exponential behavior of dynamic variables has been observed in cases of particles with a subdifussive behavior [18, 23, 24]. However, this last behavior has not been considered in water dynamic variables of molecular dynamics simulations.

In this work, using a MD method, we studied the structural and dynamics behavior of the extended simple point charge (SPC/E) model of water over a wide range of temperatures (213-360 K). As a result, we will show that the bond correlation function, P(t), has a non-exponential behavior and fits to a q-exponential function. Besides, we analyzed possible correlations between this q-exponential behavior of P(t), the hydrogen bonds distributions and the displacement of water molecules.

2. THEORY AND METHOD

We performed MD simulations of the SPC/E model of water using the GROMACS package [25, 26]. This water model assumed a rigid geometry of the molecules, with an O-H distance of 0.1 nm and a H-O-H angle of 109.47º. The hydrogen atom charge is q_H = 0.4238e and the oxygen atom charge is -2q_H . The parameters of the Lennard-Jones interactions between the oxygen atoms were those of the AMBER-99φ force field [27].

We carried out MD simulations with a constant number of water molecules (1158), in which we used the Berendsen's thermostat to apply a thermal and a hydrostatic bath to the system, obtaining isobaric-isothermal ensembles at 1 atm of pressure [28]. The values of the systems' temperatures were 360, 343, 323, 313, 303, 293, 283, 273, 263, 253, 233 and 213 K.

We assigned the velocities of the molecules according to Boltzmann's distribution in the system at 360 K. The equilibration method was similar for all systems: initially, we made a 200 steps of energy minimization, continuing with an equilibration run of ~ 1 ns using the constant potential energy as stability criterion. Then, we calculated the trajectory for an additional 3 ns. The preceding higher temperature liquid configuration was used as the starting point for each successive simulation. The values of the densities were according to the bibliography [16, 29, 30]. The simulation time step was 2 fs and the trajectories were collected every 10 fs for all simulations. The sampling time is shorter than the typical time during which a hydrogen bond can be destroyed by libration movements.

We used a geometric definition of hydrogen bond [10] with a maximum distance between their oxygen atoms of 3.5 Å and 145º as the minimum angle formed by the atoms O_donor-H- O_acceptor. In order to obtain the hydrogen bonds distribution function, in each simulation, we calculated a histogram of the quantity of hydrogen bonds for each analyzed molecule and then normalized it with the number of trajectories and the number of water molecules in the system. Then, we obtained f (n) with n = 0,1,...,5, being f (n) the probability of occurrence of n hydrogen bonds.

The bond correlation function P(t) was obtained from simulations by building a histogram of the hydrogen bonds' lifetimes for each configuration. Then, we analyzed if P(t) follows a Tsallis distribution of the form

being t the hydrogen bonds lifetime and q the nonextensivity parameter. In the limit when q → 1, exp_q(t) → exp(t).

The non-Gaussian behavior of the displacement of water molecules was studied calculating the time t*, the time at which the non-Gaussian parameter α₂(t) reaches a maximum. The non-Gaussian parameter is

where r⁴(t) and r²(t) are the fourth and second moments of the displacement distribution, respectively. α₂(t) is known to be zero for a Gaussian distribution [17, 31].

3. RESULTS AND DISCUSSION

To study the possible q-exponential behavior of the bond correlation function we followed this procedure: We assigned successive values to q (1 < q 1.3, with a step of 0.01) and made a linear fitting ¹ of log_q P(t), choosing the q value with the highest correlation coefficient r². Below the 130 fs the behavior of P(t) is strongly influenced by the libration of atoms. Above 2.5 ps the statistics is poor, because of the unlikely fact that two water molecules remain bonded by a hydrogen bond at this time. Hence, we fitted all the P(t) functions in the interval from 130 fs to 2500 fs, finding that P(t) can be acceptably fitted with a q-exponential and may be written as follows

where P₀, a and b are constants. The values of q and others parameters of the fitting for different temperatures are indicated in Table I. Fig. 1 shows, as an example, the P(t) fitting of the systems at 283 and 323 K.

The q values higher than one indicate the increase of the probability that two molecules remain bonded much longer.

Although the average bond lifetimes τ_HB obtained using Eq. 1 and 6 do not coincide with those of experimental data from polarizated light scattering experiments (see Fig. 3), both curves can be fitted by an Arrhenius behavior (τ_HB∝ exp(-E*/RT )). This fact was observed in previous works [14] [15].

The hydrogen bonds distribution functions for different temperatures are shown in Fig. 4(a). In this graphic we can distinguish three regions : The first one between 213 and ~270 K in which predominates the four hydrogen bonds probability ( f (4)) that is higher than f (3) and f (2), a second region between ~270 and ~320 K, in which f (3) > f (4) > f (2) and finally, a third region above ~320 K in which f (3) > f (2) > f (4). In the lower temperature regime almost 50 % of water molecules have four hydrogen bonds, this percentage decreases when the temperature increases, reaching 20 % when T = 360 K. In Fig. 4(b) we can observe the occurrence of a reciprocal relation between f (2), f (1) and f (4) and the temperature. These results coincide only qualitatively with the ones obtained by Sutmann and Vallauri [32]. According to the water model and the hydrogen bond definition used, amongst other things, the results from different simulations may quantitatively differ.

To study the correlation between the structural transition and the change of q we analyzed the behavior of q versus f (4) (see Fig. 4 (c)). We can observe that for temperatures below 300 K there is a linear correlation between both variables. In other words, the produced structural change below approximately 320 K causes the change in the statistical behavior of P(t) and, as a consequence, a significative increase of the probability of two water molecules be bonded much longer.

Fig. 5 shows the mean square displacements for the systems at 213, 273 and 360 K. M(t) shows some evidence of the caging behavior for water at temperatures below 300 K, where we can observe the appearance of an intermediate plateau at times between 0.1 and 100 ps that indicates a subdiffusive motion of the water. At longer times, there is a recovery of diffusive motion. We can observe from Fig. 6 (a) that t* increases with the decrease of T , according to the results obtained by Mazza et al. [17]. Clearly, t*~ T ^-1at temperatures below 300 K. This fact is attributable to the increase of the subdiffusive motion of the water.

Fig. 6 (b) shows that q and t* are correlated when the temperature drops below 300 K, which is an expected result because both parameters have a similar behavior with this variable. These changes in q and t* are produced by a modification in the water structure when f (4) outnumbers f (2) at approximately 320 K, and a structure in which prevails three and four hydrogen bonds begins to be the most frequent one in the water.

Moreover, in Fig. 6 (c) we can see the correlation between t* and f (4) for values corresponding to the systems below 300 K, in which it is observed that .

Considering these results, the four hydrogen bonds probability associated to a tetrahedral structure in water is a critical variable in the phenomenon of "cage effect" that occurs at low temperatures and the dynamics of the hydrogen bonds that changes from an exponential to a non-exponential distribution.

4. CONCLUSIONS

We have shown that the temporal correlation function of hydrogen bonds P(t) has a q-exponential behavior. The nonextensivity parameter q takes values above 1 below approximately 300 K. This increase of q indicates the increase of the probability that two molecules remain bonded during a long time t. The transition of P(t) from an exponential to a q-exponential behavior occurs in parallel with a structural modification in water, where the probability of occurrence of four hydrogen bonds outnumbers the one of two hydrogen bonds. The increase of q is also correlated with the increase of the Non-Gaussian behavior of water displacement. This fact is associated which the increase of t*, the time at which the Non-Gaussian parameter α₂(Δt) reaches a maximum.

Acknowledgment

M.G. Campo, G.L. Ferri and G.B.Roston are grateful for the financial support by PICTO UNLPAM 2005 30807 and the Facultad de Ciencias Exactas de la UNLPam.

(Received on 23 December, 2008)

• [1] D. Eisenberg and W. Kauzmann, The Structure and Properties of Water (Clarendon Press, Oxford, 1969).
• [2] F. H. Stillinger, Science 209, 451 (1980).
• [3] O. Mishima and H. E. Stanley, Nature (London) 396, 329 (1998).
• [4] A. H. Narten, M. D. Danford and H. A. Levy, Faraday Discuss., 43, 97 (1967).
• [5] A.K. Soper, F. Bruni and M. A. Ricci, J. Chem. Phys. 106, 247 (1997).
• [6] K. Modig, B. G. Pfrommer and B. Halle, Phys. Rev. Lett., 90(7), 075502 (2003).
• [7] C.A. Angell and V. Rodgers, J. Chem. Phys. 80, 6245 (1984).
• [8] J.D. Cruzan, L.B. Braly, K. Liu, M.G. Brown, J.G. Loeser and R.J. Saykally, Science 271, 59 (1996).
• [9] S. Woutersen, U. Emmerichs and H. Bakker, Science 278, 658 (1997).
• [10] A. Luzar and D. Chandler, Phys. Rev. Lett. 76, 928 (1996).
• [11] A. Luzar, J. Chem. Phys. 113, 10663 (2000).
• [12] F. Mallamace, M. Broccio, C. Corsaro, A. Faraone, U. Wandrlingh, L. Liu, C. Mou, and S.H. Chen, J. Chem. Phys. 124, 124 (2006).
• [13] C.J. Montrose, J.A. Búcaro, J. Marshall-Coakley and T.A. Litovitz, J. Chem. Phys. 60, 5025 (1974).
• [14] F.W. Starr, J.K. Nielsen, and H.E. Stanley, Phys. Rev. Lett. 82, 2294 (1999).
• [15] F.W. Starr, J.K. Nielsen, and H.E. Stanley, Phys. Rev. E. 62, 579 (2000).
• [16] S. Chatterjee, P.G. Debenedetti, F.H. Stillinger and R.M. Lynden-bell, J. Chem. Phys. 128, 124511 (2008).
• [17] M.G. Mazza, N. Giovambattista, H.E. Stanley and F.W. Starr, Phys. Rev. E 76, 031203 (2007).
• [18] B. Liu and J. Goree, Phys. Rev. Lett. 100, 055003 (2008).
• [19] T. H. Solomon, E. R. Weeks, and H. L. Swinney, Phys. Rev. Lett. , 713975 (1993).
• [20] S. Ratynskaia et al., Phys. Rev. Lett. 96, 105010 (2006).
• [21] C. Tsallis, J. of Statistical Phys. 52, 479 (1988).
• [22] E. T. Jaynes, Essays on Probability, Statistics, and Statistical Physics (Reidel, London, 1983).
• [23] C. Tsallis et al., Phys. Rev. Lett. 75, 3589 (1995).
• [24] C. Tsallis, Brazilian J. of Phys., 29(1), (1999).
• [25] H.J.C.Berendsen, J.R. Grigera, T.P. Straatsma, J. Phys. Chem. 91, 6269 (1987).
• [26] H.J.C. Berendsen, D. van der Spoel and R.V. Drunen, Comp. Phys. Comm. 91, 43 (1995).
• [27] E.J. Sorin and V.S. Pande, Biophys. J. 88 (4), 2472 (2005).
• [28] H.J.C. Berendsen, J. Postma, W. van Gusteren, A. Di Nola and J. Haak, J. Chem. Phys. 81, 3684 (1984).
• [29] L.A. Báez and P. Clancy, J. Chem. Phys. 101(11), 9837 (1994).
• [30] B. Taras and A.D.J. Haymet, Molecular Simulation 30, 131 (2003).
• [31] A. Raman, Phys. Rev. 136, 405 (1964).
• [32] G. Sutmann and R. Vallauri, J. of Mol. Liq. 98-99, 213 (2002).

1

Publication Dates

• Publication in this collection
  10 Sept 2009
• Date of issue
  Aug 2009

History

• Received
  0000
• Accepted
  23 Dec 2008