Temperature dependence of the microscopic structure and density anomaly of the SPC/E and TIP4P-Ew water models. Molecular dynamics simulation results

We have investigated temperature trends of the microscopic structure of the SPC/E and TIP4P-Ew water models in terms of the pair distribution functions, coordination numbers, the average number of hydrogen bonds, the distribution of bonding states of a single molecule as well as the angular distribution of molecules by using the constant pressure molecular dynamics simulations. The evolution of the structure is put in correspondence with the dependence of water density on high temperatures down to the region of temperatures where the system becomes supercooled. It is shown that the fraction of molecules with three and four bonds determine the maximum density for both models. Moreover, the temperature dependence of the dielectric constant is obtained and analyzed.


Introduction
It is our honor and pleasure to contribute to this special issue and to dedicate this paper to Prof.
Dr. Douglas Henderson on the occasion of his 80 th birthday. Almost four decades ago, Barker and Henderson published their seminal paper "What is "liquid"? Understanding the states of matter", that comprehensively describes basic theoretical methods to explain the equilibrium properties of simple liquids [1]. These methods have made substantial progress since that time. In particular, computer simulations techniques have become essentially powerful and common tools in the theory of liquids, mixtures and solutions. Computer simulations permitted to reach enormous progress in understanding the liquids with hydrogen bonds, more specifically, water and aqueous solutions, for needs of physical chemistry and for several inter-disciplinary areas involving complex biological molecules and membranes. In particular, Henderson with co-workers was actively interested in the properties and behavior of inhomogeneous aqueous solution, see e.g. [2][3][4].
The studies of water have long history that has been documented in an enormous number of publications. Of particular interest are the peculiar thermodynamic and dynamic properties of water manifested in a set of anomalies. In spite of formal equivalence of the relation between the microscopic structure and thermodynamic properties discussed in great detail by Barker and Henderson for simple homogeneous liquids [1], in the case of water it is inevitable to involve the notion of hydrogen bonding. Many works have dealt with the exploration of hydrogen bonds network structure, its dynamics and life-time of bonds at an early stage of computer simulations of aqueous systems, see e.g. a few original papers [5,6] as well as later reviews [7][8][9], that contain an extensive list of references on the subject. Moreover, it is known that unique properties of water specifically manifest themselves in the region of states below the freezing point corresponding to the supercooled liquid, see e.g. a recent review [10] and references therein.
It is seducing to relate at least some of water anomalies with the transformations occuring in the hydrogen bonds cooperative structure. Therefore, in this short report our focus is in the description of the evolution of the average number of H-bonds per water molecule and in the distribution of the numbers of H-bonds per molecule upon temperature changes. Moreover, we would like to relate the trends of the behavior of these properties with the density anomaly. While the properties in question describe the formation of "agglomerates" of bonded molecules at relatively short-range distances between them, we would like to get additional insights into the simultaneously occuring changes of the long-range correlations between water molecules in terms of the dielectric constant.
One of the most studied anomalies of water is its nonmonotonous temperature dependence of density, ρ(T ), see e.g. a review providing a nice historical insight [9]. Several non-polarizable water models with different number of interaction sites are capable of reproducing the presence of maximum density of water with different precision. Namely, the temperature of maximum density (TMD) at ambient pressure for most frequently used models is observed at 253 K (TIP4P), at 278 K (TIP4P/2005), at 274.15 K (TIP4P-EW), at 241 K (SPC/E) and at 182 K (TIP3P). These values for temperature are taken from the table 2 of the review by Vega and Abascal [11]. The experimental results for these and other properties used by the authors are taken from [12]. The experiment yields 277 K for the TMD of real water. On the other hand, the water density at TMD is 0.988 (TIP4P), 0.993 (TIP4P/2005), 0.994 (SPC/E), 0.98 (TIP3P), while the experimental value is 0.997. Some deviation from these values reported in different works result from technical details of each simulation.
The origin of density anomaly is commonly attributed to the evolution of hydrogen bond network with the formation of an open structure as the temperature is lowered or, on the other hand, manifests itself as the reflection of a possible liquid-liquid phase transition observed by simulations in some water models [9]. The latter and other scenarios for peculiar properties of supercooled water were documented in [10].
To summarize this brief introduction, our intention in this work is to to present our own data of molecular dynamics simulations of the temperature dependence of water density, ρ(T ), in an ample interval of temperature (from 370 K down to 170 K) at ambient pressure by using the SPC/E model. Next, we would like to compare the obtained data with the results previously reported by other authors. However, the principal focus is to follow changes of water structure in terms of different descriptors involving both the hydrogen bonds and the ρ(T ) dependence in order to establish a clear relationship between them.
Nevertheless, for the sake of comparison of our observations for the SPC/E, and to make the conclusions more sound, similar results for the TIP4P-Ew water model [13] are presented as well. Our interest in the SPC/E model is motivated by its wide use for different purposes, in particular for the description of thermodynamic properties of mixtures. Extending the previous research in mixtures, having in mind a more ample project presently being developed in this laboratory, our intention is to use the SPC/E model as a starting point in exploring water-alcohol and water-DMSO solutions by using molecular dynamics simulations.

Simulation details
All our simulations of water models were performed in the isothermal-isobaric ensemble at ambient pressure. The initial configuration was constructed with 2000 water molecules placed in a cubic array in a simulation box. Then, the system was simulated in an ample range of temperatures, from a high temperature 370 K down to 180 K. The DL-POLY Classic package was employed [14]. We used the Berendsen thermostat and barostat, the running timestep was set to 0.002 ps. Commonly, the periodic boundary conditions were used. The short range interactions were cut-off at 11 Å, whereas the long-range interactions were handled by the Ewald method with the precision 10 −4 . After equilibration, several sets of simulation runs, each for 6-8 ns, with restart from the previous configuration, were performed to obtain averages for data analysis. The overall length of simulation was dictated by the stability of internal energy, density and the dielectric constant, all being the functions of time.

Results and discussion
The temperature dependence of the density of water and the stability of several ice phases in the framework of the SPC/E model were studied in several works, see e.g. [15][16][17][18][19][20][21]. In figure 1 (a) we show the results obtained by different authors for the SPC/E water and the results of the present work. The liquid and glass branches, denominated in figure 1 (a) as Baez et al., were constructed using the points reported in table II and in table III of reference [18]. The data follow from the simulations of 360 molecules. It was assumed that the cut-off distance is at 9 Å and the effect of interactions at larger distances was neglected. The density maximum (ρ = 1.0265 g/cm 3 ) was obtained at 235 K. In the later work from that laboratory, the density values for disordered hexagonal ice (I h ) in the interval between 170 K and 220 K, were reported, see reference [20]. These data are also shown in figure 1 (a). However, it is important to mention that the results refer to the modified SPC/E model, i.e., by multiplying the interaction potential by a switching function between 6 Å and 9 Å to make the forces and the potential continuous at a cut-off distance. The modification of the potential and the application of Ewald summation technique to include the long-range interactions also resulted in the changes of the estimates for TMD, in [19,20] the improved values are ρ = 1.003 g/cm 3 and T ≈ 260 K. The algorithmic peculiarities of these series of works include a profound exploration of the effects of Ewald summation on the temperature dependence of density and other properties of the system as well as calculations of the Gibbs free energy of distinct phases which requires thermodynamic integration. The line describing the density of stable I h ice obtained by Gay et al. is reproduced from [15]. It was mentioned therein that the stability of this specific crystalline phase crucially depends on the number of molecules in the simulation cell. It is worth mentioning that the authors performed simulations in isostress, constant temperature ensemble and estimated the melting point by analyzing the time evolution of the properties of the system upon heating of the solid phase.
Quite recently, Vega et al. [22] evaluated the melting temperature of the most popular non-polarizable water models. For the SPC/E model, the melting temperature value of ice I h and the densities of liquid water and ice are shown in figure 1 (a). Our data for the liquid density of the model starting from a rather high temperature down to the region of supercooled liquid are compared with the previously published results by Bryk and Haymet [17] and with very recent results by Fennell et al. [21]. A small difference between the curve obtained in [17] and the present results in the entire interval of temperatures is presumably caused by technical details, namely by the distinct cut-off value. On the other hand, it is worth mentioning that Fennell et al. [21] slightly modified the original SPC/E model by fitting the dielectric constant value at room temperature to the experimental result. Our data locate the TMD at 245 K, i.e., in between the commonly used value 241 K, see e.g. [11,17], and the result for the optimized model by Fennell et al. 247 K. As for the value of density at TMD, our result is quite close to the data reported in [17]. In general terms, the water density values in relation to temperature grow in the interval that starts at high temperatures and is delimited below by the TMD. Next, if the temperature decreases further, the density decreases until it reaches a minimum value, T min , our data locate the minimum at T = 190 K. The minimum value of density was observed for a set of water models, e.g., for ST2 [9] and TIP4P-2005 [23]. In the case of the SPC/E model, the melting temperature T m is located in the interval between the TMD and T min , and the region between the melting point and T min actually corresponds to the supercooled states. Finally, if the temperature decreases below the T min , the density grows again as in a normal fluid, but this growth is not as rapid as in the interval of above the TMD.
In our opinion, a plausible explanation of the presence of the density minimum that involves structural changes in water was given in [9]. It refers to the behavior of a simple liquid that shrinks upon cooling due to the increasing effects of attraction between particles with a decreasing temperature. However, in a certain temperature interval, this normal behavior alters due to augmenting effects of hydrogen bonding between molecules and the resulting formation of "complexes" or "agglomerates". Cooperative behavior of these entities result in the growth of an expanded, open structure. The effects of directional hydrogen bonding overcome the effects of an attractive interaction. Consequently, the fluid density decreases in a certain interval of temperatures, evidently specific for each model. However, the directional interactions are characterized by saturation. Therefore, upon further cooling, when changes of the structural units become weaker, compared with the previous regime (as documented for example by the changes of the average number of hydrogen bonds per molecule), then the fluid density grows again similar to the temperature interval above the TMD.
However, the above discussion is incomplete and some thermodynamic arguments must be taken into account. The liquid is the stable equilibrium phase for temperatures above the melting temperature and it is unstable below a stability temperature, whereas in between, the liquid is metastable with respect to crystal ice [24]. In general, metastable liquids and specifically water, can survive even at temperatures well below the melting point until homogeneous nucleation takes place and water converts into ice [25].
Therefore, a certain part of the curve ρ(T ) presented in figure 1 (a) for the SPC/E model (below the melting point) characterized by a descending density can correspond to metastable, "long-life" states (on the simulation time scale). Similar comments concern the behavior of the TIP4P-2005 model, see figure 4 of reference [23] and the ST2 model [9,24].
Due to finite-size effects and consequent details in implementation of Ewald sums as well as due to finite sampling time, a possible spontaneous crystallization, however, is not observed in many works for a set of water models, as well as in our present simulations. Debates concerning this problem just started very recently, see e.g. [24,25]. These authors provide a comprehensive discussion of various issues relevant to the problem using their own data coming from the Monte Carlo and molecular dynamics simulations and analyzing the previously published results for different water models from the literature. Some other papers from the special issue of the Journal of Chemical Physics devoted to confined and interfacial water published in November 2014 can be useful for the reader as well.
Our final comments concerning a set of curves shown in figure 1 (a) refer to the glassy water in the framework of the SPC/E model. Actually, water, if considered below the homogeneous nucleation temperature, should be in a glassy form, see e.g., a rather comprehensive discussion of the issue in [26]. In the case of the SPC/E model, the glass transition temperature was located around 177 K [18]. At this temperature, the dependence of internal energy on temperature shows a kink (we reproduced the results of that work by performing long simulations with 2000 molecules and located the kink around T ≈ 185 K).
Thus, the augmenting density region at low temperatures refers in part to the amorphous states that cannot be distinguished in terms of the structural characteristics such as the pair distribution functions coming from the present simulations.
In figure 1 (b), we present the temperature dependence of density of the SPC/E model and compare it with our results obtained for the TIP4P-Ew model performed in the framework of the above described methodology. In addition, the results by H. Pi et al. [23] for the liquid and ice I h phases are given as well. The TIP4P-Ew model reproduces the experimental result for the TMD and an overall dependence of density at ambient pressure [27] much better compared to the SPC/E. All the models involved in this figure (SPC/E, TIP4P-Ew, TIP4P-2005) predict the existence of the density minimum and augmenting density at lower temperatures. However, it is worth mentioning that the density of ice I h substantially differs from the density of supercooled water (and of its glassy form) for the SPC/E model [see figure 1 (a)], whereas in the case of TIP4P-2005 model, the liquid phase density and that of ice are close to each other as it follows from the simulations by H. Pi et al. [23], see figure 1 (b). Similar behavior concerning the density dependences has been documented recently by J. Alejandre et al. [28] for the model from the TIP4P family developed by fitting the dielectric constant, the maximum density and the equation of state at low pressure. Unfortunately, we are unaware of the simulation results for the I h ice phase of TIP4P-Ew model at room pressure to make our discussion more general at this point. However, it seems that the differences between liquid and ice phase densities in the region of low temperatures are very sensitive to the minor changes of parameters describing different models.
The microscopic structure of a set of water models coming out from MD simulations with respect to the experimental diffraction data at room temperature and ambient pressure was discussed in great detail by Pusztai et al. [29] by using the protocol [30] combining the experimental total scattering structure factor and partial radial distribution functions from simulations. It was shown that the TIP4P-2005 model exhibits the best consistency with the experimental structure compared to other commonly used models. Nevertheless, the SPC/E, ST2 and TIP4P results are not seriously inconsistent at ambient conditions. Unfortunately, such an analysis of the structure of water has not been performed for other temperatures and pressures due to experimental difficulties.
A few examples of the radial distribution functions, g i j , where i , j stand for O and H, obtained in our NPT simulations of the SPC/E model are presented in two panels of figure 2. It can be seen that the height of the first and the subsequent maxima of both functions, g OO (r ) and g OH , increases with a decreasing temperature in the entire interval of temperatures studied. These trends are well pronounced upon cooling in the range starting from a high temperature down to around 190 K, further. At even lower temperatures, the overall shape of g i j (r ) is similar to that observed at 190 K. However, changes of the structure are much less pronounced if one moves to 180 K and 170 K. We omit these results for the sake of brevity. The functions g OO (r ) and g OH (r ) are given mostly as an illustration, because their characteristics are necessary to introduce the concept of hydrogen bonds between molecules.
In terms of the coordination number, n OO , it is appreciable that the cooling causes more pronounced shells around the chosen molecule, cf. figure 3 (a). While at high temperatures even the first coordina- tion shell is not very well pronounced, it is seen perfectly well at low temperatures, and the coordination number is around four as expected. The derivative of the coordination number is more illustrative [ figure 3 (b)]. It shows the location of the first, the second shell and even the third one, the latter being developed at very low temperatures and serving as a manifestation of a spatial order that develops in the system and involves a larger number of particles than the first neighbors.
To proceed, it is necessary to introduce the concept of H-bonds. Energetic and geometric criteria, unavoidably involving certain degree of arbitrariness, are used in the literature to interpret the structure of water and aqueous solutions. A comprehensive analysis of several geometric definitions of hydrogen bonding, leading to different values for the average number of H-bonds per molecule, has been performed in reference [31].
In this work, we use a single geometric criterion, see e.g. Zhang et al. [32] for two different water models. Three conditions determine the existence of the H-bond, namely the distance R OO between oxygen atoms belonging to two water molecules should not exceed the threshold value R C OO ; the distance R OH between the acceptor oxygen atom and the hydrogen atom connected to the donor oxygen should not exceed the threshold value R C OH . Finally, the angle O-H· · · O should be smaller than the threshold value, see e.g. references [33][34][35][36][37], concerning the analysis of the structure of liquid methanol as well as of supercritical water. If the coordinates of two molecules fulfill the above conditions, they are considered as  [32], we used the threshold constant value for the angle θ, complementary to the angle α between the vectors OH and r OH (see figure 1 of reference [31]), that is equal to 150 • . With this definition, we obtained the average number of hydrogen bonds per water molecule similar to the value by Bakó et al. [38] who used the energetic criterion for bonding.
The temperature dependence of the number of H-bonds averaged upon simulation time and normalized by the number of water molecules in the system, 〈n HB 〉 N w , following from our calculations for two water models, the SPC/E and TIP4P-Ew, is shown in figure 4 (a). This figure explains that a decreasing temperature results in a larger average number of H-bonds. It can be seen that there is a change of the slope at T ≈ 180 K for the SPC/E model where the transition into a glassy state is reported. On the other hand, the change of the slope is observed at a higher temperature, around T = 200 K, for the TIP4P-Ew model. However, this latter model is characterized by a higher number of H-bonds per molecule in the entire temperature interval studied. In both cases we see that the average number of bonds does not reach four, or, in other words, even at low temperatures, perfect tetrahedral configuration is not attained. The same conclusion can be reached by considering the "bond-angle distribution function" already studied by Bryk and Haymet [17]. It describes the probability to find certain configurations involving three oxygens belonging to different water molecules. The threshold distance is assumed to be at the first minimum of the pair distribution funcion, g OO (r ) at each temperature studied. The data were normalized by the number of molecules that one finds in the shell limited by the threshold distance. All the curves given in figure 4 (b) are characterized by two maxima. One of them at θ ≈ 53 • describes the configurations characteristic of a simple liquid [16]. On the other hand, the second maximum at θ ≈ 107 • is close to the angle θ ≈ 109.5 • describing a perfect tetrahedral configuration. At a high temperature, T = 360 K, a simple liquid type structure prevails. While the system is cooled, the first maximum decreases and becomes very 13603-7 small at T = 200 K, whereas the maximum describing the configurations close to tetrahedral substantially grows in value. It is worth mentioning that both maxima slightly shift to higher angles with a decreasing temperature. Still, even at low temperatures, the shape of the curves is not delta-like, witnessing a certain degree of disorder in the system. Moreover, we do not observe any peculiarities of this property either at TMD or at the temperature where the density is at minimum.
In order to obtain a better insight into the changes of the structure in terms of bonds with temperature, we calculate the time-averages of fractions of differently bonded molecules. This representation was already used by Bakó et al. [38] in their analysis of the simulated structures of water-methanol mixtures at room temperature. Our results concerning the SPC/E model at different temperatures are shown in figure 4 (c). We observe that at a high temperature, T = 330 K, the fractions of molecules participating in two and three bonds are the largest. Comparing this and room temperature, it can be seen that the fraction of molecules with four bonds grows at the expense of a decreasing fraction of molecules with a single bond and two bonds. The fraction of particles with three bonds remains practically unchanged.
A further decrease of temperature down to T = 240 K (note that it is very close to the TMD) leads to a substantial increase of four-bonded particles again at the expense of singly-bonded and doubly-bonded species. The fraction of triple-bonded species remains intact. If the temperature falls down to 200 K, not only the fraction of 4-bonded particles grows substantially, but the fraction of triple-bonded species exhibits a strong decrease in value. To summarize, the growth of density is accompanied by the growth of the fraction of 4-bonded particles with almost constant fraction of triple-bonded species. On the other hand, the descending density region, below the TMD, is characterized by the growth of 4-bonded species with a simultaneous decrease of fractions of triple-bonded, double-bonded and single bonded-molecules. Interestingly, the changes of the structure at even lower temperature are very small (we do not show these curves to avoid overload of the figure).
The dependence of fractions of differently bonded molecules on temperature, in the entire temperature range studied, is shown in figure 4 (d). Here, we present our results for two different models, namely for the SPC/E and TIP4P-Ew. The pattern is very similar in both cases. The TMD for each model obtained from the calculations of density on temperature is located very close to the crossing points between the curves describing the dependence of triple-bonded and 4-bonded molecules. The TMD values are marked in the figure as well. It is worth mentioning that the curves describing the fraction of triplebonded molecules is not sensitive to the temperature variable in the region between high temperatures and T ≈ 240 K. Moreover, the values for this fraction coming from the calculations of two models in question practically coincide in an ample interval of temperatures. However, the difference between these two models is well observed due to their capability of forming structures describing 4-bonded species.
In order to complete our discussion of the results given in all panels of figure 4, we would like to make the following comments. As mentioned above, the geometric definition of H-bonds permits flexibility in the sense that one can attempt to change the values and the set of parameters of H-bond definition. In this respect, the energetic definition is more restrictive. Besides the choice of a set (R OO , r OH , θ), we also tested the sets (R OO , r OH ) and (R OO , β, cf. reference [31] figure 1). These less restrictive geometric definitions lead to the results qualitatively similar to the ones discussed above. Namely, the fraction of quadruply-bonded species grows with a decreasing temperature whereas all other fractions decrease. Unfortunately, we have not observed any characteristic crossing attributable to the TMD with these definitions. Therefore, we hope to extend our study of the models in question in the framework of energetic definition of Hbonding as well.
The properties described above refer to the structure of the systems in question mostly at short interparticle separation. On the other hand, the long-range, asympthotic behavior of correlations between molecules possessing dipole moment is described by the dielectric constant, ε, see e.g. [39]. In general terms, the calculation of the dielectric constant from simulations is a difficult task [28,40,41], especially at low temperatures. Several factors can influence the result, such as e.g., the number of molecules, the type of thermostat and barostat, precision of the long-range interactions summation. Actually, very long runs are necessary to obtain reasonable estimates for this property, as it is usually calculated from the time-average of the fluctuations of the total dipole moment of the system [42] as follows: where k B is the Boltzmann constant and V is the simulation cell volume. Some of our results at different temperatures for the SPC/E model and for TIP4P-Ew are shown in figure 5 (a) and 5 (b), respectively. We would like to recall here that all the data given in figure 5 were obtained for the systems having two thousand particles. Practically all the runs lasted more than 20 ns. Still, as it can be seen, there is room for experimenting in order to get better results below the room temperature by extending the length of the runs or considering larger systems. In addition, we plan to explore the effects of changing the cutoff distance, of including the long-range corrections and the precision of Ewald summation as well as to test other thermostats. All the results for each model, however, refer to the temperatures above the respective melting point. Still, there is one delimiting factor. It manifests itself in the polarization of the samples at low temperatures. At room temperature, figure 5 (c), the projections of the total dipole moment of the system tend to zero at large times, whereas at lower temperatures, in the supercooled regime and close to the formation of glassy states, the samples inevitably fall into configuration with non-zero one or more Cartesian projections of M. Then, the calculations of the dielectric constant require different algorithms. A summarizing insight into our data for ε(T ) for two models in question are given in figure 5 (d). In general terms, the inclination of the temperature dependence of the dielectric constant is correctly described by two models, compared to the experimental results. The fitting of the dielectric constant at room temperature does not garantee the correct temperature dependence, as we can see from the results of the modified model [21].
To summarize this brief report, we have used molecular dynamic simulations in the NPT ensemble to study the density anomaly of water for the SPC/E and TIP4P-Ew models and to establish the relation between this anomaly with the microscopic structure of the system. The structure is described in terms of radial distribution functions, oxygen-oxygen coordination number and its derivative, the average number of H-bonds per molecule, angular distributions and fractions of molecules characterized by a different number of bonds. We observed that at temperature of maximum density, the fractions of molecules with three and with four bonds are almost the same, or in other words it is a reasonably good structural estimator for the TMD. In addition, we presented our preliminary results concerning the temperature dependence of the dielectric constant and compared them with experimental data. While the slope of this dependence is well described, the values are underestimated compared to the experimental results. It would be desirable in a future work to extend our structural analysis to find estimates for the cooperative bonding of molecules and relate them with the density and other anomalies of water. On the other hand, it would be interesting to establish the relationship between the structural properties describing the electric properties of the system and the dielectric constant.