Computer simulation of dynamical anomalies in stretched water

Computer simulation of dynamical anomalies in stretched water

P. A. Netz^I; F. Starr^II; M. C. Barbosa^III; H. Eugene Stanley^IV

^IDepartamento de Química, ULBRA, Canoas RS, Brasil and Departamento de Química, Unilasalle, Canoas RS, Brasil

^IICenter for Theoretical and Computational Materials Science and Polymers Division, National Institute of Standards and Technology, Gaithersburg, Maryland 20899, USA

^IIIDepartamento de Física, UFRGS, Porto Alegre RS, Brasil

^IVCenter of Polymer Studies - Boston University Boston, MA 02215, USA

ABSTRACT

In this work, we describe how the anomalous diffusivity is related to the structural anomalies. For this purpose, we study how the thermodynamics and the dynamics of low-temperature water are affected by the decrease of the density.

1 Introduction

Water is one of the most important substances in nature, and its remarkably complex behavior is puzzling, specially when we consider the simplicity of the chemical structure of water's molecule [1-4]. Water behaves anomously in several ways[5]: it expands on freezing and therefore has a negative slope in the solid-liquid equilibrium line in the P-T diagram. Under atmospheric pressure, water has a density maximum at 4ºC, a minimum in the isothermal compressibility at 46ºC, and a minimum in the isobaric heat capacity at 35ºC. It has usually high melting, boiling and critical points, among several other anomalies[6, 7].

It is known that these anomalies, as well as the main properties of water, including the properties that make water the essential fluid in biological systems, are linked to the microscopic structure of water and the peculiar intermolecular interactions. Nonetheless, it remains astonishing, how does this complex behavior emerge from a simple structure. The reproduction of the anomalous behavior and the description of these complexities is a challenge not only to the computer simulations, but also to the theoreticians that are developing theories for the mechanisms that govern complex systems. And in this aspect, computer simulations have been sucessfull in reproducing a wide range of properties of water and several anomalies, starting with very simple models.

Particularly, it is already known that there are three kinds of anomalies[8]: thermodynamic, dynamical and structural. The region in the phase diagram where the first kind of anomaly occurs is entirely inside the domain of the second, which in its turn within the region of the third. Inside the region of structural anomalies, the orientational and translational order show inter-dependence[8]. A key point is, therefore, to discover how does the structure affect the mobility.

The complex behavior of water can be approached both from the thermodynamic as well as from the microscopic level. Let us begin by focusing in the thermodynamical aspects. There are three possible explanations for the anomalous increase of the thermodynamic response functions (such as isothermal compressibility, isobaric heat capacity) on cooling, namely, the stability limit conjecture [9-11], the critical point hypothesis [12-15] and the singularity-free hypothesis [16-18]. According to the stability-limit conjecture, the increase in the thermodynamic response functions on cooling can be explained supposing that the curves are close to the spinodal, the line separating the region in the phase diagram where the liquid phase is metastable from the region where it is unstable. The pressure in the liquid spinodal in water's phase diagram, according to this conjecture, decreases on cooling, but attains a minimum value (at negative pressures). The spinodal reenters the positive pressures region of the phase diagram with further cooling, being in this way called a reentrant spinodal. The critical point hypothesis explains the anomalous increase in the response functions as being due to the proximity of a second critical point, which is located probably at - 85ºC and 230 MPa[19], at the end of a line of first-order phase transition separating two liquid phases of different density. According to the singularity-free hypothesis, there is actually no divergence in the values of the response functions, and no need to postulate a proximity of spinodal or critical point. These functions remain finite and attain maximum values. The two last hypotheses are consistent with a non-reentrant spinodal.

The origin of these anomalies are in the peculiar microscopic structure of liquid water. Water can be regarded as a transient gel[16, 20], a highly associated liquid with strongly directional hydrogen bonds. These bonds induce a local order, however not strong enough, so water maintains a long-range disorder typical of liquids. The local structure and therefore also the anomalies are enhanced when water is subject to stretching, as well as strong cooling[21, 22]. In this sense, some very important insights in water's behavior can be obtained by studying supercooled water under negative pressures (stretching). Fluids under negative pressures are relevant not only from the academic point-of-view (it is not merely an academic curiosity), but also play an important role on realistic system, such as in the transport of water in plants[23]. Experimental results [24-26] as well as computer simulations [27-35] show that, starting with atmospheric conditions of pressure, the gradual increase of pressure increases the number of defects such as multiple H-bonds and interstitial water molecules, disrupting the tetrahedral local structure and leading to a weakening of the H-bonds and therefore to an increase in the mobility [32-39]. Applying very high pressures, however, leads to steric effects which lower the mobility, and as a result the diffusion constant attains a maximum at a given density r_max. Above r_max, the diffusion of water is controlled by hindrance, with the hydrogen bonds playing a secondary role. We show that a complementary behavior is found if water is submitted to a gradual decrease in pressure, from atmospheric to negative pressures, with a minimum of diffusivity at a given density r_min, as discussed in details further in this paper.

There are many different models used in the computer simulation, and the differences between them were recently reviewed[40, 41] None of them can reproduce exactly all of water properties and anomalies. Even though, the overall thermodynamic picture that emerges from different computer simulation models is roughly the same, but the dynamic properties are slightly more model-dependent. From the large diversity of potentials, the potentials with three interaction centers are computationally faster. Among them the SPC/E[42] is particularly accurate in the description of thermodynamic and dynamical properties of water and is largely used in pure water simulations as well as in biological systems. It is also known that the SPC/E water model can reproduce the maximum in diffusivity under pressure as well as the power-law behavior of dynamical properties on cooling and therefore it seems to be a good choice for analyzing the dynamics of stretched water.

In this work, we will describe how the anomalous diffusivity is related to the structural anomalies. For this purpose, we study how the thermodynamics and the dynamics of low-temperature water are affected by the decrease of the density.

2 Methods

We performed an extensive set of molecular dynamics simulations using 216 SPC/E [42] water molecules. The Newton equations of motion were integrated using SHAKE [43, 44] with time steps of 1.0 fs for T > 210 K and 2.0 fs for T = 210 K. The range of temperatures and densities covered was 210 K < T < 280 K and 0.825 g/cm³ < r < 1.30 g/cm³. Many state points in this range have negative pressure, and at low densities they are either a metastable stretched liquid, or a phase separated liquid-gas mixture. All simulations were carried out in the canonical ensemble (NVT), in a cubic simulation box using periodic boundary conditions. The rescaling of the velocities was made using the Berendsen thermostat [45], the electrostatic interactions were calculated using reaction field [46] with cut-off radius of 0.79 nm.

The diffusion coefficient D was calculated from the asymptotic slope of the mean square displacement plotted versus time:

The orientational relaxation was analyzed using the rotational autocorrelation functions [43]:

The vector e is a chosen unity vector describing the orientation of the dipole moment. Other possible choices of unity vectors such as the O-H bond direction and the vector perpendicular to the plane of the molecule give similar results[48]. The orientational relaxation time was calculated from these functions using an biexponential decay fit function. [47]

The first term correspond to a fast librational motion, and the two relaxation times take into account the possibility of a fast and a slow reorientation process. However, in the most of the cases taking only one exponential yielded a reasonable result, and thus only one relaxation time was determined.

In order to understand the effect of the structure on the dynamics we also carried out a detailed analysis of the local structure of water [49] computing the distribution of the number and angle of H-bonds (O-H ... O) and distribution of the number and angle of first neighbor molecules, analyzing the angle O ... O ... O. The number of first neighbors is calculated by computing all molecules whose distance from a given molecule is less than 0.32 nm. This distance represents the first minimum in the oxygen-oxygen radial distribution function g_OO (r) at moderate and low densities. In order to have a better comparison, the same O-O distance criterion is also applied for all state points, irrespective the density.

The distribution of the number of hydrogen bonds is computed in a similar way. A hydrogen bond between two water molecules is counted if an O-O pair within 0.32 nm and a H...O with a bound distance less than 0.25 nm are found[50]. This criterium is similar to the geometrical definition of the hydrogen bonds [51, 52] but modified in the sense that no restriction about hydrogen bond angle was applied.

3 Results

3.1 Thermodynamic Properties

Table 1 shows the thermodynamic properties obtained in our MD simulations. A crucial point in the distinction between the different thermodynamic scenarios, the stability-limit conjecture the critical point and the singularity-free hypotheses, is the shape of the spinodal. Computer simulations help us to locate the spinodal. We fit the P-r isotherms with a fifth-order polynomial and locate the minima in each isotherm, where the isothermal compressibility is zero. These minima define the stability limit beyond which no equilibrium exists and therefore it locates the spinodal P_sp(T). We found a non-reentrant spinodal[53], qualitatively similar to that found by Harrington and coworkers[31]. The same result was also found by Yamada and coworkers [54-56] using the TIP5P[57] potential. This behavior can support both the critical point hypothesis and the singularity-free scenario.

The density at the spinodal r_sp(T) was calculated as well and is located in a narrow range 0.853 < r _sp< 0.874. For our results that means that our simulated points whose density is less than 0.875 either should be ruled out or care must be taken in order to check the evidence of cavitation.

3.2 Dynamic Properties

The translational diffusion coefficient (D) and the orientational relaxation time () of SPC/E water were analyzed. Table 2 shows the dynamic properties of the simulated state points. Fig. 1 shows D along isotherms. For T < 260 K, D has a minimum at about r » 0.9 (g/cm³), which becomes more pronounced at lower temperatures. The rotational diffusion (orientational relaxation time) displays a very similar trend, as shown in Fig. 2. The behavior of along isotherms is in fact complementar to the behavior of D. The relaxation time initially increases, with decreasing density, passes through a maximum and decreases with further stretching. the maximum in can be found in the same region as the minimum in D. The product ×D (not shown) is remarkably constant irrespective the temperature or density[58, 50].

These results show that the stretching has influence both in the translational as well as in the rotational diffusion. The complementar behavior of translational and rotational diffusion should be interpreted not in terms of hydrodynamic arguments[59], but rather as a clue pointing out that both the translational diffusion as well as the orientational process involve microscopic rearrangements, i.e. both motions are correlated by a common mechanism. In order to better understand this mechanism, the detailed microscopic structure of water was also investigated.

3.3 Structural Properties

One important aspect about molecular dynamics simulations is the very detailed microscopic information that these simulations can give. For instance, it is possible to look at the neighbor structure of each water molecule at each simulation step, obtaining a statistical picture of the local structure. Several aspects can be investigated, such as the distribution of the number of neighbors and number of H-bonds and also the angle distribution of neighbor molecules (O ... O ... O) angle and angle distribution of H-bonds.

Table 3 shows the distribution of the number of neighbors and H-bonds for the simulated state points. Fig. 3 shows the H-bond and O ... O ... O angle distributions for the temperature T = 240 K, at several densities, whereas the same properties, analysed fixing the density at r = 0.90 g/cm³ and at r = 1.125 g/cm³ at several temperatures were shown in Figs. 4 and 5, respectively.

The detailed analysis of the local structure of water reveals the enhancement of the tetrahedrality at low temperatures and ice-like densities. This can be detected from the angle distributions shown in these figures. Lowering the density, the O ... O ... O angle distribution becomes more sharply peaked around the tetrahedral angle, indicating an enhancement of the ice-like structure. the peak around 60º, on the other side, is related to a fifth neighbor in the coordination shell and the presence of this fifth neighbor increases with increasing density, as shown in Fig. 3b and Table 3. At very high densities, a peak appears around 80º, related to the appearance of a 6th neighbor[50]. The O ... O ... O angle distribution, at a fixed density, seems to depend only weakly on temperature, as shown in Figs. 4b and 5b.

The H-bond angle distribution, indicated in Fig. 3a seems not to depend on density, as far as this density is not too high. Fixing the density and analyzing the effect of temperature (as in Figs. 4a and 5a) we see that the peaks become broad with increasing T, indicating only a thermal fluctuation. Only at very high densities does the H bond become distorted.

The increase of the number of molecules with higher coordination numbers at densities ranging between r_min and r_max is responsible for the increase on the diffusion coefficient with increasing density. The number of H-bonds, however, does not alter significantly, meaning that these imperfections arise from the inclusion of extra molecules in the coordination shell which either do not for H-bonds or share and H-bond with another molecule. This shared bonds are weakened and the molecule can connect to another molecule by means of a small rotation, keeping the remaining H-bonds. The net motion of water structure can be described as a slow non-oscillatory or quasi-oscillatory translational and rotational displacements[60]. This peculiar mechanism can explain the coupling between translational and rotational diffusion, without the need of invoking hydrodynamical arguments.

4 Conclusions

In conclusion, we propose that the anomalous dynamic behavior of supercooled water can be explain in terms of its structural properties: the number and the distribution of hydrogen bonds and and the number and angular distribution of the O-O-O neighbors. We calculate these quantities for a wide range of densities and temperatures. We observe that the number of molecules with four neighbors has a maximum at the density r_min where the translational diffusion has a minimum and the rotational diffusion is maximum. For r_max, the density of highest translational diffusion, the number of molecules with more than four neighbors (that is, molecules with extra neighbors) is comparable to the number of molecules with exactly four neighbors. The number of intact hydrogen bonds is almost not affected by the increase of the number of molecules with higher coordination number. This indicates the presence of imperfections in the network what is supported by the presence of distortions in the H-O-H bond angles at higher densities (see Figs. 3a, and 5a). These new neighbor molecules provide extra oxygens that will share an hydrogen with another oxygen without forming an extra hydrogen bond. As a result, the interaction becomes weaker and might break. By a small rotation the molecule can connect to another molecule and by making small rotations, it is able to diffuse. This mechanism explains the coupling between rotation and translation. This fast diffusion is, however, limited to a certain range of densities that enable the molecule to make another bond once one had broken. At r_min, for instance, there is a very low probability to find these extra molecules in the neighborhood of a given molecule, and therefore the translational and rotational diffusion are slow.

Received on 02 September, 2003

• [1] D. Eisenberg and W.J. Kauzmann. The Structure and Properties of Water Clarendon, London, 1969.
• [2] P. Ball. Life's Matrix: A Biography of Water Farrar Straus and Giroux, New York, 2000.
• [3] V. Brazhkin, S.V. Buldyrev, V.N Ryzhov, and H.E. Stanley, editors. New Kinds of Phase Transitions: Transformations in disordered systems, Dordrecht, 2002. NATO Advanced Research Workshop, Volga River, Kluwer.
• [4] O. Mishima and H.E. Stanley. Nature, 392, 164 (1998).
• [5] Martin Chaplin. Water structure and behavior, http://www.sbu.ac.uk/water.
• [6] E.W. Lang and H.-D. Lüdemann. Angewandte Chemie, International Edition in English, 21, 315 (1982).
• [7] R.C. Dougherty and L.N. Howard. Journal of Chemical Physics, 109, 7379 (1998).
• [8] J. R. Errington and P. G. Debenedetti. Nature, 409, 318 January (2001).
• [9] R.J. Speedy. Journal of Physical Chemistry, 86, 982 (1982).
• [10] R.J. Speedy. Journal of Physical Chemistry, 86, 3002 (1982).
• [11] R.J. Speedy. Journal of Physical Chemistry, 91, 3354 (1987).
• [12] P. H. Poole, F. Sciortino, U. Essmann, and H. E. Stanley. Nature, 360:324, november 1992.
• [13] P.H. Poole, F. Sciortino, U. Essmann, and H.E. Stanley. Physical Review E, 48, 3799 (1993).
• [14] F. Sciortino, P.H. Poole, U. Essmann, and H.E. Stanley. Physical Review E, 55, 727 (1997).
• [15] S. Harrington, R. Zhang, P.H. Poole, F. Sciortino, and H.E. Stanley. Physical Review Letters, 78, 2409 (1997).
• [16] H.E. Stanley and J. Teixeira. Journal of Chemical Physics, 73, 3404 (1980).
• [17] S. Sastry, P.G. Debenedetti, F. Sciortino, and H.E. Stanley. Physical Review E, 53, 6144 (1996).
• [18] L.P.N. Rebelo, P.G. Debenedetti, and S. Sastry. Journal of Chemical Physics, 109, 626 (1998).
• [19] S.B. Kiselev and J.F. Ely. Journal of Chemical Physics, 116, 5657 (2002).
• [20] A. Geiger, F. H. Stillinger, and A. Rahman. Journal of Chemical Physics, 70, 4185 (1979).
• [21] P.G. Debenedetti. Metastable Liquids. Princeton University Press, Princeton, 1996.
• [22] J.T. Fourkas, D. Kivelson, U. Mohanty, and K.A. Nelson. Supercooled Liquids: Advances and Novel Applications ACS Books, Washington, 1997.
• [23] W. T. Pockman, J. S. Sperry, and J. W. O'Leary. Nature, 378, 715 December (1995).
• [24] J. Jonas, T. DeFries, and D.J. Wilbur. Journal of Chemical Physics, 65, 582 (1976).
• [25] F. X. Prielmeier, E. W. Lang, R. J. Speedy, and H. D. Lüdemann. Physical Review Letters, 1987.
• [26] F.X. Prielmeier, E.W. Lang, R.J. Speedy, and H.-D. Ludemann. Berichte Bunsengeselschaft fur Physikalische Chemie, 92, 1111 (1988).
• [27] M. R. Reddy and M. Berkowitz. Journal of Chemical Physics, 1987.
• [28] F. Sciortino, A. Geiger, and H.E. Stanley. Nature, 354(november):218, november 1991.
• [29] F. Sciortino, A. Geiger, and H. E. Stanley. Journal of Chemical Physics, 1992.
• [30] L. A. Báez and P. Clancy. Journal of Chemical Physics, 101, 9837 (1994).
• [31] S. Harrington, P.H. Poole, F. Sciortino, and H.E. Stanley. Journal of Chemical Physics, 107, 7443 (1997).
• [32] P. Gallo, F. Sciortino, P. Tartaglia, and S.-H. Chen. Physical Review Letters, 76, 2730 (1996).
• [33] F. Sciortino, P. Gallo, P. Tartaglia, and S.-H. Chen. Physical Review E, 1996.
• [34] S.-H. Chen, P. Gallo, F. Sciortino, and P. Tartaglia. Physical Review E, 1997.
• [35] F. Sciortino, L. Fabbian, S.-H. Chen, and P. Tartaglia. Physical Review E, 1997.
• [36] T. Yamaguchi, S.-H. Chong, and F. Hirata. Journal of Chemical Physics, 119, 1021 (2003).
• [37] F. W. Starr, S. Harrington, F. Sciortino, and H. E. Stanley. Physical Review Letters, 1999.
• [38] F.W. Starr, F. Sciortino, and H.E. Stanley. Physical Review E, 60, 6757 (1999).
• [39] A. Scala, F.W. Starr, E.La Nave, F. Sciortino, and H.E. Stanley. Nature, 406, 166 (2000).
• [40] B. Guillot. Journal of Molecular Liquids, 101, 219 (2002).
• [41] A. G. Kalinichev. Reviews in Mineralogy and Geochemistry. Mineralogical Society of America, Washington DC, 2001.
• [42] H.J.C. Berendsen, J.R. Grigera, and T.P. Straatsma. Journal of Physical Chemistry, 91, 6269 (1987).
• [43] M.P. Allen and D.J. Tildesley. Computer Simulation of Liquids Clarendon Press, Oxford, 1987.
• [44] J.P. Ryckaert, G. Ciccotti, and H.J.C. Berendsen. Journal of Computational Physics, 23, 327 (1977).
• [45] H.J.C. Berendsen, J.P.M. Postma, W.F.van Gunsteren, A. DiNola, and J.R. Haak. Journal of Chemical Physics, 81, 3684 (1984).
• [46] O. Steinhauser. Molecular Physics, 45, 335 (1982).
• [47] Y. Yeh and C.-Y. Mou. Journal of Physical Chemistry B, 103, 3699 (1999).
• [48] P.A. Netz, F.W. Starr, H.E. Stanley, and M.C. Barbosa. In New kinds of phase transitions: transformation in disordered substances, volume 81 of NATO Science Series II, page 417. NATO, 2002.
• [49] T. Head-Gordon and F. H. Stillinger. Journal of Chemical Physics, 98, 3313 (1993).
• [50] P.A. Netz, F.W. Starr, M.C. Barbosa, and H. Eugene Stanley. Physica A, 314, 470 (2002).
• [51] M. Mezei and D. L. Beveridge, Journal of Chemical Physics 74, 1981 (1981).
• [52] D.C. Rappaport. Molecular Physics, 50, 1151 (1983).
• [53] P. A. Netz, F. W. Starr, H. Eugene Stanley, and M. C. Barbosa. Journal of Chemical Physics, 115, 344 (2001).
• [54] M. Yamada, S. Mossa, H.E. Stanley and F. Sciortino Phys. Rev. Lett, 88, 195701 (2002).
• [55] H. E. Stanley, M.C. Barbosa, S. Mossa, P.A. Netz, F. Sciortino, F.W. Starr and M. Yamada Physica A, 315, 281 (2002).
• [56] H. E. Stanley, M.C. Barbosa, S. Mossa, P.A. Netz, F. Sciortino, F.W. Starr and M. Yamada Water at positive and negative pressures In Liquids under Negative Pressures, volume 84 of NATO Science Series II, NATO, 2002.
• [57] M. W. Mahoney and W. L. Jorgensen Journal of Chemical Physics, 112, 8910 (2000).
• [58] P. A. Netz, F. W. Starr, M. C. Barbosa, and H. Eugene Stanley. Journal of Molecular Liquids, 101, 159 (2002).
• [59] S. Ravichandran, A. Perera, M. Moreau, and B. Bagchi. Journal of Chemical Physics, 107, 8469 (1997).
• [60] H. Tanaka and I. Ohmine. Journal of Chemical Physics, 87 6128 (1987).

Publication Dates

• Publication in this collection
  27 Apr 2004
• Date of issue
  Mar 2004

History

• Received
  02 Sept 2003