- Email: [email protected]
Supercooled water: A molecular dynamics simulation study with a polarizable potential
Journal of Molecular Liquids 127 (2006) 28 – 32 www.elsevier.com/locate/molliq
Supercooled water: A molecular dynamics simulation study with a polarizable potential Manuela Minozzi, Paola Gallo ⁎, Mauro Rovere Dipartimento di Fisica, Università “Roma Tre”and INFM Democritos National Simulation Centre, Via della Vasca Navale 84, 00146 Roma, Italy Available online 26 May 2006
Abstract The behaviour of supercooled water in an extended region of the phase diagram, from ambient conditions to negative pressure and low temperatures, is studied by computer simulation employing a polarizable water potential. The structural analysis of the system shows that the polarizable potential chosen well reproduces the water properties even in extreme conditions and allows us to perform a more detailed thermodynamical analysis. In order to infer features of the phase diagram relevant for the explanation of the anomalous behaviour of metastable water, the density dependence of pressure was analyzed to investigate the pathway of the spinodal line with the decreasing temperature. © 2006 Elsevier B.V. All rights reserved. PACS: 61.20.Ja; 61.20.-t; 64.60.My
1. Introduction Water is the most important liquid for life, it is the only chemical compound that occurs naturally in the solid, liquid and vapour phases, and the only naturally occurring inorganic liquid on earth [1]. Besides its fundamental role in many areas of science, it has physical properties which are often uncharacteristic of many other liquids. The peculiar behaviour of these properties is often ascribed to the structure of water, which consists of three dimensional arrangements of hydrogen-bonded molecules. Yet the real meaning and origin of the hydrogen bond interaction is still unclear. In contrast to ordinary liquids water expands on cooling below 4 C at atmospheric pressure, which corresponds to the phase diagram point at which the isobaric temperature dependence of density reaches its maximum. The density anomaly is associated with some thermodynamic anomalies, for instance it was observed [2] that the magnitude of the thermal expansion coefficient, of the isothermal compressibility and isobaric heat capacity increases markedly with decreasing temperature. For thermodynamic consistency the increase in isothermal compressibility upon isobaric cooling and the decrease in the magnitude of the negative thermal expansion ⁎ Corresponding author. E-mail address: [email protected] (P. Gallo). 0167-7322/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2006.03.006
coefficient upon isothermal compression are inseparable from the existence of a negatively sloped temperature of maximum density (TMD) [3]. The lowest temperatures at which those quantities are measurable are up to now around T = − 38 C at 1 bar, due to the strong tendency of metastable water to crystallize. If the magnitude of the thermal expansion coefficient, of the isothermal compressibility and of the isobaric heat capacity are extrapolated below this temperature they appear to diverge at Ts ∼− 45 C [4], that corresponds to the ‘‘absolute’’ limit of stability. A coherent theory of the thermodynamic and transport properties of supercooled water does not yet exist, also due to the difficulties encountered in experiments. Several theories have been developed in the past decade to explain at a thermodynamic level the anomalous behaviour of supercooled water and many numerical studies have been conducted on this subject [5]. In fact, since there are no experimental data, numerical simulations play an important role to investigate the behaviour of supercooled water. For instance, many water models were studied and compared [6] to each other in order to have a realistic ‘computer water’ and reliable data. In 1992 Poole et al. [7], performing several molecular dynamical simulations with ST2 water potential, discovered the possible existence of a second critical point in the range of T = − 50 C and P = 1 kbar. In this scenario the approach to this critical point would give rise to the experimentally observed increase in the response functions.
M. Minozzi et al. / Journal of Molecular Liquids 127 (2006) 28–32
Below this second critical point, the liquid phase would separate into two distinct phases: a low density liquid (LDL) and a high density liquid (HDL), this interpretation is consistent with the observed behaviour of amorphous water. Water, in fact, can apparently form at least two distinct forms of glass: low density amorphous ice (LDA) and high density amorphous ice (HDA). In 1985 Mishima [8] compressed LDA and observed a sharp transition at 6 ± 0.5 kbar between the two amorphous states showing the existence of a first-order transition between LDA and HDA. The two distinct supercooled liquid phases can be identified with the experimentally observed low density amorphous and high density amorphous ice. Other hypothesis were formulated on the basis of simulated data. A lattice model with nearest-neighbour and directional attractions of Rebelo et al. [9] shows a scenario where the response functions remain finite, and no underlying singularity is invoked. This scenario is called the singularity free interpretation (SFI). The SFI assumes that upon isobaric cooling the thermodynamic response functions go through a maximum but remain finite. An anomalous increase in isothermal compressibility, heat capacity, and thermal expansion is explained in the SFI by the existence in the supercooled water of a density maxima locus, which is negatively sloped in the P–T plane. According to this scenario, no phase transition or critical point occurs at low temperatures. It has been also demonstrated that the TMD lines of both ST2 [10] and TIP4P [6] water models change direction in the negative pressure region, proving that the first hypothesis formulated in literature, the stability limit conjecture (SLC), does not provide an adequate phase diagram for these water models. The SLC hypothesis is based on thermodynamic consistency arguments and it was proposed by Speedy and Angell in 1982 [11]. According with this scenario the stability limits for superheating and supercooling of water would be the observed branches of a continuous line of stability limits PS(T). This spinodal line starting from the critical point, first traces the limit of stability of superheated water, then passes through the negative pressure region of the P–T diagram, defining the limits of stability for stretched water, and finally would turn upward again to positive pressures in supercooled water. The spinodal locus estimated by Speedy and Angell turns around at the temperature of 303 K (30 °C) and pressure of − 212 MPa, and reaches positive pressures at the temperature of 228 K (− 45 °C). Debenedetti and D'Antonio [12] generalized this thermodynamic analysis showing that loss of tensile strength implies density anomalies in the vicinity of the spinodal curve.
better than the non-polarizable water model from ambient to supercritical water. It has been also shown that the BSV model reproduces the maximum density somewhat better than the nonpolarizable models and the experimental temperature dependence of the density better than other polarizable models over a broad temperature range [6–14]. For those characteristics we decided to analyze the behaviour of supercooled water using the BSV model in the unexplored region of the phase diagram: at low temperatures and negative pressures approaching the spinodal line. Th BSV model has the same geometry as the TIP4P : three sites are placed on the atoms of the molecule, the H atoms carry fractional positive charges of q = + 0.499 e, while a forth site is localized along the bisector H-O-H bond angle, at a distance of 0.15 Å from the O atom and carries the negative fractional charge of − 2q. The total energy of the system is calculated as a summation of Lennard-Jones and Coulomb interactions between point charges. An induced dipole moment located in the center-of-mass of the molecule of water is explicitly introduced to describe the effect of the environment on the electric field. The induced dipole pi = α Ei is calculated from the local electric field Ei with an iterative procedure by assuming an isotropic polarizability fixed at α = 1.44 A3. Incorporating molecular polarizability into a water model enables the molecules to react to local perturbations in the electrostatic field. All long range dipolar interactions are estimated by the reaction field approximation with suitable tapering. The system is composed of 256 molecules in a cubic box with periodic boundary conditions. The simulations have been performed in the NVT ensemble with the minimum image convention and a cut-off of the interactions at half of the box length. For every fixed value of the density ρ we decrease the temperature from T = 400 K to T = 220 K, exploring a wide region of the negative pressures phase diagram. We analyzed low density water, i.e. ρ goes from 1.0 to 0.83 gr/cm3. During the equilibration and production phases snapshots of the system were visualized to check the structural configuration of the liquid [see Fig. 1]. The time step for the integration of the molecular trajectories is δt = 1 fs.
2. Computer simulations We performed Molecular Dynamics (MD) computer simulations of water with a polarizable potential model, proposed by Brodholt–Sampoli–Vallauri (BSV) [13] and widely studied in literature [14]. According to a recent detailed analysis performed by comparing the simulation with the more recent neutron diffraction data [15], this polarizable model is able to reproduce changes in the hydrogen bond network of water
29
Fig. 1. Snapshot of the simulated system T = 220 K ρ = 0.87 gr/cm3.
30
M. Minozzi et al. / Journal of Molecular Liquids 127 (2006) 28–32
3. Structure and thermodynamic behaviour of supercooled water 3.1. Structure We have calculated the g(r) of the thermodynamical points simulated averaging over (1–3) × 103 configurations. In Fig. 2 we show the g(r) of the system for ρ = 0.9 gr/cm3. We focused our attention to the behaviour of the g(r) on decreasing temperature. At the top of Fig. 2 we report the gOO(r) at every temperature (T = 300 K to T = 210 K), the functions are translated to solve the different lines. The peak positions agree with the experimentally measured structure of water at T = 300 K [15]. At the bottom we can observe the first peak of gOH(r), 1.9 Å, due to hydrogen bonds, the second peak of gOH(r), and the first two peaks of the gHH(r), (2.4, 3.8 Å) showing the
tetrahedral arrangement. Upon decreasing temperature the peak positions remain at the same place and appear sharper, the tendency of the system to undergo a more ordered configuration is also evidenced by the behaviour of the coordination number (not shown) which converges to the limiting value of the LDA ice (n = 4). Different densities do not show substantial differences. Close to a mechanical stability limit cavitation phenomena could be present in the liquid [3]. To check the homogeneity of the system at low temperature we analysed the probability distribution of the local density. We calculated the local density by dividing the system in 64 small boxes and counting the number of particles present per each box. We then calculated the probability distribution of the local density, in fact a dishomogeneity in the liquid would cause a double peak in the distribution analogous to a phase separation.
Fig. 2. Radial distribution functions gOO(r), gOH(r), gHH(r) at different temperatures from T = 300 K to T = 210 K, for the system at r = 0.90 gr/cm3. gOO(r) are shifted on the y-axe to solve the different curves from higher (T = 300 K in the bottom) to lower (T = 210 K) temperatures. The peaks of the three atomic radial distributions become sharper as T decreases.
M. Minozzi et al. / Journal of Molecular Liquids 127 (2006) 28–32
31
Fig. 3. Probability distributions of the local density. The local density is defined as the occupation number of small sub-boxes of the simulation box (see text for details).
In Fig. 3 we report the probability distribution calculated for ρ = 0.83–0.87–0.90–1.0 gr/cm3, we show both the high and the low temperature data (T = 300 K and T = 220 K) in the same graph in order to follow the behaviour of the distributions on approaching the limit of mechanical stability of the system (see the thermodynamical results). Fig. 3 shows that approaching the lower temperatures a single maximum in the probability distribution function is still observed around the average density value, it proves that there are no holes neither small bubbles and therefore all our thermodynamical data refer to a homogeneous liquid system. A single maximum in the probability distribution function is observable around the average density value.
Fig. 4. The spinodal line (⁎) is shown with the Pρ(T) isochores for ρ = 0.90– 0.95–1.0 gr/cm3.
3.2. Thermodynamic behaviour In the phase diagram space the spinodal line is the boundary between metastable and unstable states, the limit of metastability of the system. The properties of the spinodal in a mean-field approach are illustrated by the classical van der Waals description of a liquid–gas system. We can identify the first minimum of the van der Waals loop in the P–V plane as the limit of stability of the systems. According to the meanfield approach states for which the second derivative of free energy is positive are either stable or metastable, while states for which it is negative are unstable. Furthermore for every temperature we can define the point were the isothermal compressibility turns from positive to negative values as the
Fig. 5. The statistical efficiency s is calculated by the relation shown in the figure from the plateau value to calculate the bar error (see text for details).
32
M. Minozzi et al. / Journal of Molecular Liquids 127 (2006) 28–32
limit of stability of the system, the set of these points at different temperatures corresponds to the spinodal line. If we consider the proportional relation between the inverse of the coefficient of isothermal compressibility and the slope of isotherms of P versus V, we can identify the points of the spinodal line with the minima of the isotherms in the P–V, and also in the P–ρ, plane. Then we extrapolate the liquid spinodal from our simulated data as the minima points of the isotherms of P versus ρ for T = 300–280–260–240–230– 220–210 K. In Fig. 4 we report the spinodal line and the Pρ(T) isochores for ρ = 0.90–0.95–1.0 gr/cm3. We observe that down to the lowest temperature investigated the spinodal line decreases with T and no retracing behaviour is observable, analogous to the non-polarizable water model potentials while both the singularity free and the second critical point scenario are consistent with our data. We can assume that the pressure values have a Gaussian distribution in order to have an estimate of the error bar. To evaluate the correlation time, the sequence of MD steps is broken up into nb blocks each of length tb (n = tb·nb). For each block we calculate the average and the variance, we assume that the variance will be inversely proportional to tb, for large values of tb, when the blocks are sufficiently large to appear uncorrelated. Our aim is to estimate the variance of the entire block in this situation. Therefore we define the statistical efficiency parameter s [16] as the limiting ratio of the variance observed on the blocks and the variance calculated in the hypothesis of the uncorrelated Gaussian statistics (σ (P)). We calculated s as the plateau value in Fig. 5 and then the variance σ(
Supercooled water: A molecular dynamics simulation study with a polarizable potential Manuela Minozzi, Paola Gallo ⁎, Mauro Rovere Dipartimento di Fisica, Università “Roma Tre”and INFM Democritos National Simulation Centre, Via della Vasca Navale 84, 00146 Roma, Italy Available online 26 May 2006
Abstract The behaviour of supercooled water in an extended region of the phase diagram, from ambient conditions to negative pressure and low temperatures, is studied by computer simulation employing a polarizable water potential. The structural analysis of the system shows that the polarizable potential chosen well reproduces the water properties even in extreme conditions and allows us to perform a more detailed thermodynamical analysis. In order to infer features of the phase diagram relevant for the explanation of the anomalous behaviour of metastable water, the density dependence of pressure was analyzed to investigate the pathway of the spinodal line with the decreasing temperature. © 2006 Elsevier B.V. All rights reserved. PACS: 61.20.Ja; 61.20.-t; 64.60.My
1. Introduction Water is the most important liquid for life, it is the only chemical compound that occurs naturally in the solid, liquid and vapour phases, and the only naturally occurring inorganic liquid on earth [1]. Besides its fundamental role in many areas of science, it has physical properties which are often uncharacteristic of many other liquids. The peculiar behaviour of these properties is often ascribed to the structure of water, which consists of three dimensional arrangements of hydrogen-bonded molecules. Yet the real meaning and origin of the hydrogen bond interaction is still unclear. In contrast to ordinary liquids water expands on cooling below 4 C at atmospheric pressure, which corresponds to the phase diagram point at which the isobaric temperature dependence of density reaches its maximum. The density anomaly is associated with some thermodynamic anomalies, for instance it was observed [2] that the magnitude of the thermal expansion coefficient, of the isothermal compressibility and isobaric heat capacity increases markedly with decreasing temperature. For thermodynamic consistency the increase in isothermal compressibility upon isobaric cooling and the decrease in the magnitude of the negative thermal expansion ⁎ Corresponding author. E-mail address: [email protected] (P. Gallo). 0167-7322/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2006.03.006
coefficient upon isothermal compression are inseparable from the existence of a negatively sloped temperature of maximum density (TMD) [3]. The lowest temperatures at which those quantities are measurable are up to now around T = − 38 C at 1 bar, due to the strong tendency of metastable water to crystallize. If the magnitude of the thermal expansion coefficient, of the isothermal compressibility and of the isobaric heat capacity are extrapolated below this temperature they appear to diverge at Ts ∼− 45 C [4], that corresponds to the ‘‘absolute’’ limit of stability. A coherent theory of the thermodynamic and transport properties of supercooled water does not yet exist, also due to the difficulties encountered in experiments. Several theories have been developed in the past decade to explain at a thermodynamic level the anomalous behaviour of supercooled water and many numerical studies have been conducted on this subject [5]. In fact, since there are no experimental data, numerical simulations play an important role to investigate the behaviour of supercooled water. For instance, many water models were studied and compared [6] to each other in order to have a realistic ‘computer water’ and reliable data. In 1992 Poole et al. [7], performing several molecular dynamical simulations with ST2 water potential, discovered the possible existence of a second critical point in the range of T = − 50 C and P = 1 kbar. In this scenario the approach to this critical point would give rise to the experimentally observed increase in the response functions.
M. Minozzi et al. / Journal of Molecular Liquids 127 (2006) 28–32
Below this second critical point, the liquid phase would separate into two distinct phases: a low density liquid (LDL) and a high density liquid (HDL), this interpretation is consistent with the observed behaviour of amorphous water. Water, in fact, can apparently form at least two distinct forms of glass: low density amorphous ice (LDA) and high density amorphous ice (HDA). In 1985 Mishima [8] compressed LDA and observed a sharp transition at 6 ± 0.5 kbar between the two amorphous states showing the existence of a first-order transition between LDA and HDA. The two distinct supercooled liquid phases can be identified with the experimentally observed low density amorphous and high density amorphous ice. Other hypothesis were formulated on the basis of simulated data. A lattice model with nearest-neighbour and directional attractions of Rebelo et al. [9] shows a scenario where the response functions remain finite, and no underlying singularity is invoked. This scenario is called the singularity free interpretation (SFI). The SFI assumes that upon isobaric cooling the thermodynamic response functions go through a maximum but remain finite. An anomalous increase in isothermal compressibility, heat capacity, and thermal expansion is explained in the SFI by the existence in the supercooled water of a density maxima locus, which is negatively sloped in the P–T plane. According to this scenario, no phase transition or critical point occurs at low temperatures. It has been also demonstrated that the TMD lines of both ST2 [10] and TIP4P [6] water models change direction in the negative pressure region, proving that the first hypothesis formulated in literature, the stability limit conjecture (SLC), does not provide an adequate phase diagram for these water models. The SLC hypothesis is based on thermodynamic consistency arguments and it was proposed by Speedy and Angell in 1982 [11]. According with this scenario the stability limits for superheating and supercooling of water would be the observed branches of a continuous line of stability limits PS(T). This spinodal line starting from the critical point, first traces the limit of stability of superheated water, then passes through the negative pressure region of the P–T diagram, defining the limits of stability for stretched water, and finally would turn upward again to positive pressures in supercooled water. The spinodal locus estimated by Speedy and Angell turns around at the temperature of 303 K (30 °C) and pressure of − 212 MPa, and reaches positive pressures at the temperature of 228 K (− 45 °C). Debenedetti and D'Antonio [12] generalized this thermodynamic analysis showing that loss of tensile strength implies density anomalies in the vicinity of the spinodal curve.
better than the non-polarizable water model from ambient to supercritical water. It has been also shown that the BSV model reproduces the maximum density somewhat better than the nonpolarizable models and the experimental temperature dependence of the density better than other polarizable models over a broad temperature range [6–14]. For those characteristics we decided to analyze the behaviour of supercooled water using the BSV model in the unexplored region of the phase diagram: at low temperatures and negative pressures approaching the spinodal line. Th BSV model has the same geometry as the TIP4P : three sites are placed on the atoms of the molecule, the H atoms carry fractional positive charges of q = + 0.499 e, while a forth site is localized along the bisector H-O-H bond angle, at a distance of 0.15 Å from the O atom and carries the negative fractional charge of − 2q. The total energy of the system is calculated as a summation of Lennard-Jones and Coulomb interactions between point charges. An induced dipole moment located in the center-of-mass of the molecule of water is explicitly introduced to describe the effect of the environment on the electric field. The induced dipole pi = α Ei is calculated from the local electric field Ei with an iterative procedure by assuming an isotropic polarizability fixed at α = 1.44 A3. Incorporating molecular polarizability into a water model enables the molecules to react to local perturbations in the electrostatic field. All long range dipolar interactions are estimated by the reaction field approximation with suitable tapering. The system is composed of 256 molecules in a cubic box with periodic boundary conditions. The simulations have been performed in the NVT ensemble with the minimum image convention and a cut-off of the interactions at half of the box length. For every fixed value of the density ρ we decrease the temperature from T = 400 K to T = 220 K, exploring a wide region of the negative pressures phase diagram. We analyzed low density water, i.e. ρ goes from 1.0 to 0.83 gr/cm3. During the equilibration and production phases snapshots of the system were visualized to check the structural configuration of the liquid [see Fig. 1]. The time step for the integration of the molecular trajectories is δt = 1 fs.
2. Computer simulations We performed Molecular Dynamics (MD) computer simulations of water with a polarizable potential model, proposed by Brodholt–Sampoli–Vallauri (BSV) [13] and widely studied in literature [14]. According to a recent detailed analysis performed by comparing the simulation with the more recent neutron diffraction data [15], this polarizable model is able to reproduce changes in the hydrogen bond network of water
29
Fig. 1. Snapshot of the simulated system T = 220 K ρ = 0.87 gr/cm3.
30
M. Minozzi et al. / Journal of Molecular Liquids 127 (2006) 28–32
3. Structure and thermodynamic behaviour of supercooled water 3.1. Structure We have calculated the g(r) of the thermodynamical points simulated averaging over (1–3) × 103 configurations. In Fig. 2 we show the g(r) of the system for ρ = 0.9 gr/cm3. We focused our attention to the behaviour of the g(r) on decreasing temperature. At the top of Fig. 2 we report the gOO(r) at every temperature (T = 300 K to T = 210 K), the functions are translated to solve the different lines. The peak positions agree with the experimentally measured structure of water at T = 300 K [15]. At the bottom we can observe the first peak of gOH(r), 1.9 Å, due to hydrogen bonds, the second peak of gOH(r), and the first two peaks of the gHH(r), (2.4, 3.8 Å) showing the
tetrahedral arrangement. Upon decreasing temperature the peak positions remain at the same place and appear sharper, the tendency of the system to undergo a more ordered configuration is also evidenced by the behaviour of the coordination number (not shown) which converges to the limiting value of the LDA ice (n = 4). Different densities do not show substantial differences. Close to a mechanical stability limit cavitation phenomena could be present in the liquid [3]. To check the homogeneity of the system at low temperature we analysed the probability distribution of the local density. We calculated the local density by dividing the system in 64 small boxes and counting the number of particles present per each box. We then calculated the probability distribution of the local density, in fact a dishomogeneity in the liquid would cause a double peak in the distribution analogous to a phase separation.
Fig. 2. Radial distribution functions gOO(r), gOH(r), gHH(r) at different temperatures from T = 300 K to T = 210 K, for the system at r = 0.90 gr/cm3. gOO(r) are shifted on the y-axe to solve the different curves from higher (T = 300 K in the bottom) to lower (T = 210 K) temperatures. The peaks of the three atomic radial distributions become sharper as T decreases.
M. Minozzi et al. / Journal of Molecular Liquids 127 (2006) 28–32
31
Fig. 3. Probability distributions of the local density. The local density is defined as the occupation number of small sub-boxes of the simulation box (see text for details).
In Fig. 3 we report the probability distribution calculated for ρ = 0.83–0.87–0.90–1.0 gr/cm3, we show both the high and the low temperature data (T = 300 K and T = 220 K) in the same graph in order to follow the behaviour of the distributions on approaching the limit of mechanical stability of the system (see the thermodynamical results). Fig. 3 shows that approaching the lower temperatures a single maximum in the probability distribution function is still observed around the average density value, it proves that there are no holes neither small bubbles and therefore all our thermodynamical data refer to a homogeneous liquid system. A single maximum in the probability distribution function is observable around the average density value.
Fig. 4. The spinodal line (⁎) is shown with the Pρ(T) isochores for ρ = 0.90– 0.95–1.0 gr/cm3.
3.2. Thermodynamic behaviour In the phase diagram space the spinodal line is the boundary between metastable and unstable states, the limit of metastability of the system. The properties of the spinodal in a mean-field approach are illustrated by the classical van der Waals description of a liquid–gas system. We can identify the first minimum of the van der Waals loop in the P–V plane as the limit of stability of the systems. According to the meanfield approach states for which the second derivative of free energy is positive are either stable or metastable, while states for which it is negative are unstable. Furthermore for every temperature we can define the point were the isothermal compressibility turns from positive to negative values as the
Fig. 5. The statistical efficiency s is calculated by the relation shown in the figure from the plateau value to calculate the bar error (see text for details).
32
M. Minozzi et al. / Journal of Molecular Liquids 127 (2006) 28–32
limit of stability of the system, the set of these points at different temperatures corresponds to the spinodal line. If we consider the proportional relation between the inverse of the coefficient of isothermal compressibility and the slope of isotherms of P versus V, we can identify the points of the spinodal line with the minima of the isotherms in the P–V, and also in the P–ρ, plane. Then we extrapolate the liquid spinodal from our simulated data as the minima points of the isotherms of P versus ρ for T = 300–280–260–240–230– 220–210 K. In Fig. 4 we report the spinodal line and the Pρ(T) isochores for ρ = 0.90–0.95–1.0 gr/cm3. We observe that down to the lowest temperature investigated the spinodal line decreases with T and no retracing behaviour is observable, analogous to the non-polarizable water model potentials while both the singularity free and the second critical point scenario are consistent with our data. We can assume that the pressure values have a Gaussian distribution in order to have an estimate of the error bar. To evaluate the correlation time, the sequence of MD steps is broken up into nb blocks each of length tb (n = tb·nb). For each block we calculate the average and the variance, we assume that the variance will be inversely proportional to tb, for large values of tb, when the blocks are sufficiently large to appear uncorrelated. Our aim is to estimate the variance of the entire block in this situation. Therefore we define the statistical efficiency parameter s [16] as the limiting ratio of the variance observed on the blocks and the variance calculated in the hypothesis of the uncorrelated Gaussian statistics (σ (P)). We calculated s as the plateau value in Fig. 5 and then the variance σ(
run) using the relation [17]: σ(
run) = (s / nstep)1 / 2 × σ (P). Error bars are reported in Fig. 4. 4. Conclusions We have performed several Molecular Dynamic simulations with a polarizable potential model, to investigate the role of the polar behaviour of water in a region of phase diagram experimentally unreachable. Our structural results indicate that the BSV polarizable potential describes well structural properties of liquid water even at low T and negative P. The system, in fact, is homogeneous also at lower temperatures and pressures. These results allow us to perform the thermodynamic analysis on the sample. Our thermodynamic results indicate that with decreasing temperature the spinodal line is monotonically descendent, no retracing behaviour is observable, so the BSV water behaviour
is consistent with the singularity free or with the second critical point scenario. These results are also consistent with recent thermodynamical results reported in the supercooled region for the same potential [18]. References [1] S. Franks, Water: a Matrix for Life, Royal Society of Chemistry, Cambridge UK, 2000. [2] F.X. Prielmeier, E.W. Lang, R.J. Speedy, H.D. Lu¨demann, Phys. Rev. Lett. 59 (1987) 1128; R.J. Speedy, C.A. Angell, J. Chem. Phys. 65 (1976) 851; G.S. Kell, J. Chem. Eng. Data 20 (1975) 97. [3] P.G. Debenedetti, J. Phys., Condens. Matter 15 (2003) 1670; P.G. Debenedetti, Metastable Liquids. Concepts and Principles, Princeton NJ, University Press, 1996. [4] H. Kanno, C.A. Angell, J. Chem. Phys. 70 (1979) 4008. [5] M. Yamada, S. Mossa, H.E. Stanley, F. Sciortino, Phys. Rev. Lett. 88 (2002) 195701; P.H. Poole, F. Sciortino, U. Essman, H.E. Stanley, Phys. Rev., E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 55 (1997) 727; S. Harrington, R. Zhang, P.H. Poole, F. Sciortino, H.E. Stanley, Phys. Rev. Lett. 78 (1996) 2409; P.H. Poole, F. Sciortino, U. Essman, H.E. Stanley, Phys. Rev., E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 48 (1993) 3799. [6] P. Jedlovszky, J. Richardi, J. Phys. Chem. 110 (1999) 8019; K. Kiyohara, K.F. Gubbins, A.Z. Panagiotopoulos, Mol. Phys. 94 (1998) 803; W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. [7] P.H. Poole, F. Sciortino, U. Essmann, H.E. Stanley, Nature 360 (1992) 324. [8] O. Mishima, J. Chem. Phys. 100 (1994) 5910. [9] L.P.N. Rebelo, P.G. Debenedetti, S. Sastry, J. Chem. Phys. 109 (1998) 626; S. Sastry, P.G. Debenedetti, F. Sciortino, H.E. Stanley, Phys. Rev., E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 53 (1996) 6144. [10] F.H. Stillinger, A. Rahman, J. Chem. Phys. 60 (1974) 1545. [11] H. Kanno, C.A. Angell, J. Chem. Phys. 70 (1979) 4008; R.J. Speedy, C.A. Angell, J. Chem. Phys. 65 (1976) 851. [12] M.C. D'Antonio, P.G. Debenedetti, J. Chem. Phys. 86 (1987) 2229; P.G. Debenedetti, M.C. D'Antonio, J. Chem. Phys. 84 (1986) 3339; J. Chem. Phys. 85 (1986) 4005. [13] J. Brodholt, M. Sampoli, R. Vallauri, Mol. Phys. 86 (1995) 149. [14] P. Cristofori, P. Gallo, M. Rovere, Mol. Phys. 103 (2005) 501; V. De Grandis, P. Gallo, M. Rovere, J. Chem. Phys. 118 (2003) 3646; P. Jedlovszky, R. Vallauri, J. Chem. Phys. 115 (2001) 3750; P. Jedlovszky, M. Mezei, R. Vallauri, Chem. Phys. Lett. 318 (2000) 155; P. Jedlovszky, R. Vallauri, Mol. Phys. 97 (1999) 1157; I.M. Svishchev, P.G. Kusalik, J. Wang, R.J. Boyd, J. Chem. Phys. 105 (1996) 4742. [15] A.K. Soper, J. Chem. Phys. 101 (8) (1994). [16] R. Friedberg, J.E. Cameron, J. Chem. Phys. 52 (1970) 6049. [17] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, 1987 Oxford. [18] P. Jedlovszky, R. Vallauri, J. Chem. Phys. 122 (2005) 081101; Phys. Rev., E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 67 (2003) 011201.