Structural and dynamic properties of SPC/E water

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 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 parameters, 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 alpha_2(t) of the mean square displacement of water molecules. When the temperature decreases, the translational (tau) 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 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][2][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 (infra-red absorption and Raman scattering [4], depolarized light scattering [5,6], neutron scattering [7], femtosecond spectroscopy [8][9][10][11] and other techniques [12][13][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 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 τ 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 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][22][23]. τ HB is obtained in MD using the history-dependent bond correlation function P (t), which represents the probability that an hydrogen bond formed at time t = 0 remained continuously unbroken and breaks at time t [24,25].
Also, the dynamics of water can be studied by analyzing the mean-square displacement time series 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 intermediate 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 * characterizes this caging effect (see Sec. II for more details) [26,27].
In a previous work, we found a q-exponential behavior in P (t), in which q increases with T −1 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 configuration 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 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 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 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 O donor -H-O acceptor of 145 • .
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 first hydration layer (the hydration number) where ρ is the number density.
The translational order parameter, τ , is defined in Ref. [16] as f (2) f (3) f (4) f ( where the dimensionless variable s ≡ rn 1/3 is the radial distance r scaled by the mean intermolecular distance n 1/3 , and S c corresponds to half of the simulation box size. The orientational order parameter Q is defined as [15] where θ jik is the angle formed by the atoms Here, O i is the reference oxygen atom, and O j and O k 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 configuration. Then, I fitted 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 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 parameter α 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 environment (caging effect), and this prevents it from reaching the diffusive state [26,27].

III. Results and discussion
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 >    2 shows the oxygen-oxygen RDFs corresponding 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 associated 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 behaviors 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 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 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 quantifies the changes in the average angle 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 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 (∝ 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 significantly for Q ≥ 0.67. A similar  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.
change occurs for Q ≥ 0.67 in the linear correlation 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 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. 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 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 first hydration layer decreases from N ∼ 5 to N ∼ 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 t * show the occurrence of long lasting hydrogen bonds and caging effect between the molecules of the system.