Computational Investigation of Order, Structure, and Dynamics in

J. Phys. Chem. B 2005, 109, 6527-6534

6527

Computational Investigation of Order, Structure, and Dynamics in Modified Water Models† R. M. Lynden-Bell*,‡ and Pablo G. Debenedetti§ UniVersity Chemical Laboratory, Lensfield Road, Cambridge CB2 1EW, U.K., and Department of Chemical Engineering, Princeton UniVersity, Princeton, New Jersey 08544 ReceiVed: September 14, 2004; In Final Form: October 29, 2004

Model liquids have been constructed to study the role of local structure in the anomalous properties of liquid water. The intermolecular potentials were modified by increasing the weight of the Lennard-Jones term relative to the electrostatic term in the SPC/E model for water. The resulting family of liquids varies from SPC/E water to a Lennard-Jones-like liquid. Properties were measured as a function of density and temperature. The local structure was described by two order parameters, one measuring the tetrahedral order and the other measuring the translational order. The translational order parameter was found to be large for both tetrahedral and Lennard-Jones liquids, but to go through a minimum as the potentials were modified, demonstrating that the two types of structure are incompatible. Just as in water several properties (e.g., the translational diffusion coefficient, entropy) exhibit anomalous density dependence as a result of the breakdown of local tetrahedrality, we observed nonmonotonic behavior of the translational diffusion constant and reorientational relaxation rate as the fluids were transformed from tetrahedral to Lennard-Jones-like. This is also an indication of the incompatibility between Lennard-Jones and water-like structure.

1. Introduction Liquid water has many unusual properties when compared with simple organic liquids and these are associated with the distinctive local structure around each molecule. In water there are 4-5 nearest neighbors and the four nearest occupy tetrahedral positions determined by the hydrogen-bonding between electropositive H atoms on one molecule and electronegative oxygen atoms on another. In liquid argon there is no local orientational order and the number of molecules in the first solvation shell is much larger (12-13). The unusual properties of liquid water, such as the density maximum as a function of temperature, the increase in specific heat (Cp) as the temperature is reduced, and an increasing diffusion constant as the density is increased,1,2 can be attributed to the breakdown of tetrahedral order as either the temperature or density is increased. Classical simulation methods allow one to relate macroscopic to microscopic phenomena and to investigate changes in both as a function of density and temperature. One can also investigate the effects of changing the potential parameters to try to understand the microscopic origin of different phenomena. In this paper we describe an investigation of the changes in properties of liquids in which the relative contributions of the Lennard-Jones and the electrostatic contributions to the potential are varied systematically. An important observation from this work is that as the weight of the Lennard-Jones contribution to the intermolecular potential is increased with respect to the electrostatic terms, the local order in the liquid evolves from tetrahedral to spherically symmetric. Because these types of order are topologically incompatible, the evolution from one to the other gives rise to an intermediate regime of minimum order. These statements will be given precise, quantitative meaning

after we introduce the metrics for structural order. The reason for the minimum in translational order is the breakdown of structure outside the first shell. This work complements earlier work done by one of us3 in which the anomalies of water were mapped in the temperaturedensity plane, and work done by the other author4 on the local structural properties of geometrically modified water models. We choose as the basic system for our study the SPC/E model of water,5 and vary systematically the relative strength of its Lennard-Jones and electrostatic terms. Water-like anomalies in the SPC/E model are shifted to lower temperatures compared to experiment (e.g., the temperature of maximum density at atmospheric pressure is +4 °C for water and -31 °C for the SPC/E model6). Accordingly, the SPC/E model is said to be “understructured” relative to real water. In contrast, the ST27 model is “overstructured”, and anomalies persist to higher temperatures than in real water (e.g., the temperature of maximum density at atmospheric pressure is shifted to 49 °C6). In this sense, these two models bracket the behavior of real water, and we choose as a starting point for this investigation the understructured limit. 2. Simulation Details 2.1. Model. The intermolecular energy Φ in the SPC/E model5 is written as a sum of pairwise intermolecular interactions. This can be expressed as

Φ)

[VLJ(rij) + VES(rij,Ωij)] ∑ pairs

(1)

where rij is the separation between oxygen atoms on molecules i and j and Ωij describes the relative orientation of the molecules. The Lennard-Jones interaction is

†

Part of the special issue “David Chandler Festschrift”. * Corresponding author E-mail: [email protected]. ‡ University Chemical Laboratory. § Princeton University.

[(σr ) - (σr ) ]

VLJ(r) ) 4

10.1021/jp0458553 CCC: $30.25 © 2005 American Chemical Society Published on Web 12/30/2004

12

6

(2)

6528 J. Phys. Chem. B, Vol. 109, No. 14, 2005

Lynden-Bell and Debenedetti

and the electrostatic interaction acts between all pairs of charges zm, zn on oxygen and hydrogen atoms belonging to different molecules

VES(r) )

zmzn 4π0r

(3)

where 0 is the vacuum permittivity. The charges on oxygen and hydrogen atoms are -0.8476e and +0.4238e, respectively, the values of the Lennard-Jones parameters are ) 0.6502 kJ mol-1 and σ ) 3.169 Å, the OH bond lengths are 1 Å, and the HOH bond angle is 109.47°. In this model the first term of the Lennard-Jones interaction provides short-range repulsion, preventing the molecules from getting too close to each other. The second term, which lacks orientation dependence, describes the intermolecular attraction due to the dispersion forces. The hydrogen bonds are described by the electrostatic terms that favor local tetrahedral structure. In our family of models, the relative weights of the LennardJones and electrostatic components are changed by scaling the Lennard-Jones term by a factor γ, so that the potential energy Φ becomes

Φ ) γVLJ + VES

(4)

A related transformation can be obtained by rescaling the charges by a factor λ with λ2 ) γ-1, in which case the potential energy would be

V(λ) ) λ2Φ ) VLJ + λ2VES

(5)

The configurational properties of the system are determined by the ratio of the potential energy to kT, and can be written in alternative ways:

Φ/kT ) (γVLJ + VES)/kT ) (VLJ + λ2VES)/λ2kT

(6)

showing that the static properties of the simulations can be interpreted as describing a system in which the Lennard-Jones contribution is scaled up by a factor of γ ()λ-2) or one where the electrostatic interactions and the temperature are scaled down by λ2 (λ < 1). The reason for choosing the scaling by γ rather than by λ at constant temperature is illustrated in Figure 1, which shows the variation of Φ(r) and V ) λ2Φ(r) for two water molecules approaching to form a hydrogen bond. The curves show how the potentials change as the scaling is varied. It can be seen that the well depth of Φ in the upper curves changes slowly, whereas the well depth of λ2Φ in the lower curves decreases rapidly. The result is that if λ is decreased at constant temperature, the liquid changes from a tetrahedral liquid to a supercritical Lennard-Jones liquid, whereas if γ is increased at constant temperature to values of 4 or 5 the liquid changes to a normal dense Lennard-Jones liquid. It is this latter transformation that we require. The main effect of the Lennard-Jones term in SPC/E water is to provide enough repulsion to prevent the hydrogen and oxygen sites approaching each other too closely. As the weighting, γ, of the Lennard-Jones term increases, the optimum OH distance increases and the hydrogen bonds are weakened. The interactions become more isotropic and the liquid changes continuously from being tetrahedral (or water-like) to spherical (or Lennard-Jones-like). The state points and models can be described by the three parameters, γ, number density (F ) N/V), and temperature T. Simulations of liquids at various densities and values of the

Figure 1. Potential energy curves of a water dimer: above, V ) Φ; below, V ) λ2Φ. In each graph the curves from left to right are for λ ) 1/xγ ) 1, 0.9, 0.8, 0.7, 0.6 and 0.15, respectively. Note the rapid variation in well depth in the lower part of the figure.

scaling parameter γ were carried out at temperatures of T ) 272 K and T ) 220 K. These temperatures were chosen to include examples where SPC/E water shows anomalies in diffusion3,8,9 and, at the lower temperature, a negative thermal expansion. 2.2. Details of Runs. A number of molecular dynamics runs were performed using a program adapted from the DLPOLY code.10 Each run, specified by the three parameters γ, T, and F/F0, was first equilibrated and then average quantities collected over 0.5 ns. Here F0 denotes the reference density of water, 1 g cm-3 ) 55 mol L-1 ) 0.033 Å-3. Long range corrections were included for both the Lennard-Jones and the electrostatic interactions, the latter being calculated using the Ewald method. In all simulations the molecular dynamics cell contained 256 molecules. Constant volume FCC periodic boundary conditions were used with distances between a molecule and its nearest image ranging from 23.22 Å (F/F0 ) 0.85) to 20.42 Å (F/F0 ) 1.25). Densities are reported in terms of F/F0. The time step was 2 fs, and a Nose´ thermostat with a relaxation time of 0.5 ps was applied to maintain the temperature. Average values of energies, pressure, orientational order parameter, and radial distribution functions were accumulated during the run. In addition, a trajectory file was written every 0.5 ps for later processing to determine diffusion constants, reorientational rates, translational order parameters, and other properties. A few

Order, Structure, and Dynamics in Modified Water Models

Figure 2. Variation of the tetrahedral order parameter q (points connected by solid lines) and the translational order parameter t (points connected by dashed lines) as a function of γ at T ) 272 K for two densities.

shorter runs in which the trajectory file was written more frequently were made for determining linear and angular velocity time correlation functions. Runs were carried out for various values of the scaling parameter γ ranging from 1 (SPC/E water) to 2.78 (λ ) 0.6) and covering the (γ, F) plane for the two temperatures, 272 K and 220 K. 2.3. Order Parameters. The local tetrahedral order at each oxygen site i is defined by a sum over the angles subtended by the four nearest neighbors at the oxygen site3,11 6

qi ) 1 -

[

3

∑ j 0, also have (∂s/∂F)T > 0. This is the result of the thermodynamic identity

F

(∂F∂s) ) V(∂F∂p) (∂T∂F) T

T

(10)

p

As (∂F/∂P)T > 0 for mechanical stability, the signs of (∂F/∂T)P and (∂s/∂F)T are the same. The value of the partial derivative was determined from the simulation data from the relationship

1 ∂U/N ) [( - pF ] (∂S/N ∂F ) T ∂F ) -2

T

T

(11)

where U/N is the potential energy per molecule. Values of its derivative with respect to density were determined by fitting a polynomial to results from state points with different densities. 3. Results 3.1. Structural Order. Figure 2 shows the tetrahedral and translational order parameters as a function of the scaling parameter, γ, for two different densities at T ) 272 K. As expected the tetrahedral order parameter, q, decreases as γ increases. The translational order parameter, t, initially decreases as the tetrahedral order is destroyed, but then increases as the translational order associated with the Lennard-Jones liquid increases. Note that in the low-gamma regime (γ ≈ 1, tetrahedral behavior) both order metrics decrease with increasing density, whereas for γ > ≈1.4 the translational order increases and the tetrahedral order decreases on compression. This change in translational order can be understood from the radial distribution functions shown in Figure 3. The lowest curve shows the oxygen-oxygen g(r) for SPC/E water at 272 K and F/F0 ) 1. The structure is typical of a tetrahedral liquid with the second shell at about 1.7 times the distance of the first shell. The top curve corresponds to a liquid at the same temperature and density but with γ ) 2.78 and is typical of a Lennard-Jones liquid with the second shell at twice the distance of the first shell. The intervening curves show that, for



Figure 5. Plot of the γ-density plane showing regions of tetrahedral and spherical order at two temperatures. The solid lines refer to T ) 272 K and the dashed lines to T ) 220 K. The loci of maximum q are also shown.

Figure 4. Variation of q (solid lines) and t (dashed lines) as a function of density for various values of γ.

intermediate values of γ, the second peak in the position characteristic of a tetrahedral liquid disappears and is gradually replaced by a second peak in the position characteristic of a Lennard-Jones liquid. When γ ) 1.23 and 1.38, there is very little structure beyond the first peak, and consequently, the t order parameter is low. The minimum value of t at this density is near γ ) 1.23 although the structure beyond the first peak is slightly less at γ ) 1.38. These changes can be interpreted as showing the breakdown of the tetrahedral order characteristic of SPC/E water and its replacement by the characteristic order of the Lennard-Jones liquid. As these two types of order are incompatible, t goes through a minimum. The variation of the order parameters as a function of density at constant γ is shown in Figure 4. As found in Errington and Debenedetti’s work on SPC/E water,3 q has a maximum and t a minimum as a function of density at constant temperature. The maximum value of q occurs at low densities and negative pressures. At lower densities the liquid can no longer sustain its optimum tetrahedral structure and soon reaches the point where it becomes unstable with respect to cavitation. The minimum in t as a function of density has the same explanation as the minimum in t seen as a function of γ: both are a consequence of the breakdown of liquid structure as the liquid transforms from tetrahedral to Lennard-Jones in character. In (γ, F) planes at constant temperature, the loci of minima in t that result from either varying γ at constant density and

Figure 6. Distributions of numbers of nearby hydrogens (squares) and nearby oxygen atoms (circles) for a tetrahedral liquid (γ ) 1; points connected by solid lines) and a Lennard-Jones liquid (γ ) 2.78; points connected by dashed lines). The state points both have T ) 272 K and F/F0 ) 0.85. Note the change in the distribution of neighboring oxygen atoms.

temperature or from varying density at constant γ and temperature coincide. The line of of these minima divide the (γ, F) plane into regions where the liquid is tetrahedral and regions where it is a Lennard-Jones liquid. Examples are shown for two temperatures in Figure 5. Note that the region where the liquid displays predominantly tetrahedral behavior shrinks in size as the temperature is increased. The loci of maximum q resulting from isothermal compression at constant γ are also shown in Figure 5. As the liquid becomes less tetrahedral the number of nearby oxygens increases. Figure 6 shows histograms of the distribution functions of numbers of hydrogens within 2.5 Å of an oxygen site and the number of oxygens within 3.5 Å for liquids at the same density and temperature but with γ values characteristic of tetrahedral (γ ) 1) and Lennard-Jones (γ ) 2.78) behavior, respectively. The tetrahedral liquid has q ) 0.67 whereas the Lennard-Jones liquid has q ) 0.51. At the state point with γ )


J. Phys. Chem. B, Vol. 109, No. 14, 2005 6531

Figure 8. Variation of the diffusion coefficient D with γ at two densities at T ) 272 K. Note that the maximum in D occurs at higher values of γ than the minimum translation order.

Figure 7. Plots of F0(∂(s/k)/∂F)T as a function of F/F0 for a number of values of γ and T. Note that, although this quantity usually shows a maximum, it only reaches positive values for low values of T and γ.

1 the most probable number of nearby hydrogen atoms is 2 and that of oxygen atoms is 4, as expected for tetrahedral coordination. At the state point with higher γ and lower q the most probable number of oxygen atoms increases to 6 and the width of the distribution becomes larger. The distribution of numbers of neighboring hydrogen atoms changes less. 3.2. Density Anomalies and Entropy. Figure 7 shows plots of (∂s/∂F)T as a function of density for a number of values of γ and two temperatures. Although (∂s/∂F)T only reaches positive values at the lowest temperature, there is a maximum in this function for low values of γ even at 272 K. At 220 K there is a region in the (γ, F) plane where the overall order measured by the entropy decreases with density. In normal liquids the entropy decreases with density at constant temperature as the free volume for each molecule is reduced. The anomalous behavior in the liquids studied here is due to the disruption of tetrahedral order on compression. At high enough temperatures or high enough values of γ the tetrahedrality is so weakened as to prevent (∂s/∂F)T from becoming positive despite the initial increase in this quantity on compression at γ ) 1 and γ ) 1.23. As stated earlier, where (∂s/∂F)T > 0, the density behaves anomalously, increasing with temperature at constant pressure. Thus the region of density anomalies in these liquids is where the curves in Figure 7 are on the positive side of the y axis.The density maxima are where the curves recross the y axis. 3.3. Dynamical Processes. Figure 8 shows the variation of the diffusion constant, D, as a function of γ for two densities.

Figure 9. Variation of the diffusion constant D with density at T ) 220 K and a number of values of γ. The anomalous region in which D increases with density gets smaller and moves to lower densities as γ increases.

As γ increases, initially the molecules become more mobile as a result of the disruption of tetrahedral order, and the diffusion rate increases. Eventually, the liquids become more structured and the molecules less mobile, with D reaching a maximum value. This maximum, however, occurs well after the point at which t goes through a minimum, showing that the LennardJones type of structure needs to be well-established before translational motion is inhibited. At sufficiently low temperatures9 the SPC/E model shows both minima and maxima in the diffusion constant as a function of density. Normally diffusion decreases with density, so the region between the diffusion minimum and maximum can be identified as the region of diffusion anomaly. Figure 9 shows that this anomalous region moves to lower densities as γ is increased at T ) 220 K. At T ) 272 K the region is smaller and has disappeared at γ ) 1.56. This behavior as γ is varied at constant temperature is very similar to the behavior found by Errington and Debenedetti3 as temperature is increased at γ ) 1.



Figure 10. Variation of the relaxation rate τ2-1 as a function of density at T ) 220 K. Note the anomalous increase of this rate on compression for most of the density range shown.

Figure 11. Velocity correlation functions for molecules initially with high (full line) and low (dashed line) tetrahedral order. Note the difference in the region of the first minimum.

The reorientation rate also increases as γ increases and the liquid becomes less tetrahedral. Figure 10 shows the variation of the reorientation rate τ2(z)-1 of the dipole axis (z) as a function of density at T ) 220 K for three values of γ. τ2 is defined in section 2.4. Although generally the rate increases with γ, eventually it goes through a maximum as a function of density at constant γ as the liquid takes on the Lennard-Jones type of structure. The correlation of increased mobility with local order may be probed by investigating the mobility of individual molecules as a function of their site orientational order parameter qi. Figure 11 shows the velocity correlation functions for a sample of molecules with low tetrahedral order (qi < 0.6) and a sample with high tetrahedral order (qi > 0.75) for SPC/E water at 220 K and density F/F0 ) 0.85. The main difference is the depth of the first minimum of the velocity autocorrelation function, which almost disappears for the molecules with low qi but is quite deep for those molecules with high tetrahedral order. This implies that the dynamics of molecules with higher local

Figure 12. Anomalies in the (γ, F) plane. These plots show regions where different types of anomalies are found. Solid symbols denote maxima and open symbols minima in diffusion (squares) and reorientation rates (triangles). Open and filled circles denote the lower and upper density boundaries of the region within which the liquid exhibits a negative coefficient of thermal expansion (density anomaly). Lines are a guide to the eye.

tetrahedral order tend to be more solidlike at short times, and those with less local tetrahedral order are more mobile. As the mean square displacement at time t is proportional to the integral of this correlation function, molecules with high local tetrahedral order move less far on average than do those with low tetrahedral order, at short times. 3.4. Anomalies. We have identified three types of anomalies in this work, namely, regions where the density increases with temperature, regions where the diffusion constant increases with density, and regions where the reorientation rate increases with density. Figure 12 shows the anomalous regions on the (γ, F) plane for T ) 220 K and T ) 272 K. At the lower temperature one can see regions where the density behaves anomalously, where the diffusion behaves anomalously and where the reorientation behaves anomalously. It can be seen that the regions of reorientation, diffusion, and density anomalies are progressively smaller and enclose one another. The line of minimum t is close to the upper bound of the diffusion anomaly. At the higher temperature, the anomalous regions have all shrunk and the region of density anomaly has disappeared.


Figure 13. Space of possible state points for the modified water models. The axes cross at γ ) 1, F ) 0.8, and T ) 220 K. Labeled regions: W, tetrahedrally ordered liquids; L, liquids with LennardJones order; D, fluids with little order. The boundary between regions W and L is defined by a surface where (∂t/∂F)γ,T and (∂t/∂γ)F,T are both equal to zero.

4. Discussion and Conclusions In this work we have constructed a family of artificial liquids that vary from tetrahedral to spherical (Lennard-Jones). This was done in a controlled way by varying the parameter γ, the factor by which the Lennard-Jones contribution to the pair potential is scaled up relative to the electrostatic terms. The state points of individual simulations can be plotted in a threedimensional space defined by γ, F and T. The degree of order was measured by two parameters, q and t, the first of which measures the local tetrahedrality of any site, and the second measures the translational order, which is high in both tetrahedral and Lennard-Jones-like liquids. This latter order parameter goes through a minimum value as one moves from a state point with high tetrahedral order (low γ, F, and T) to one with high Lennard-Jones order (high γ, F, and T). The surface of such minima can be used to separate regions where liquids exhibit tetrahedrally and spherically ordered structure and behavior. Figure 13 shows a sketch of part of the γ, F, T space of state points illustrating these ideas. The back face is a γ-T plane corresponding to a density of about 0.8F0, which is near the limit of metastability of these liquids under tension. The lefthand side of the figure is a F-T plane with γ ) 1, which is standard SPC/E water. Three regions are distinguished by the letters W, L, and D. In region W near the origin (γ ) 1, F ) 0.8, and T ) 220 K) the liquid shows tetrahedral properties (water-like). In region L the liquid is translationally ordered in the Lennard-Jones structure. As the temperature is raised, the liquids become less structured and more akin to supercritical fluids. This region is denoted by D. The surface of minimum translational order (t) is shown separating the W region from the L region. A diffuse boundary is also shown separating the L region from the D region. The region of anomalies is close to the surface separating the W and L regions, but mostly on the tetrahedral side. 4.1. Order. We have used three measures of order, q, the tetrahedral order, t, the translational order and S, the entropy. These measure different aspects of order and a comparison of their behavior demonstrates the peculiarities of water. The entropy provides a very broad measure of the degree of disorder in a liquid. In general for liquids one expects it to increase as the temperature is increased and the molecules explore more of the possible configurational space and to decrease as the density is increased and the accessible configurational space decreases. The tetrahedral order parameter, q, decreases with density, γ and temperature as the space in Figure 13 is traversed from

J. Phys. Chem. B, Vol. 109, No. 14, 2005 6533 points in the vicinity of the label W through the boundary between regions W and L. It also decreases as the density is reduced toward the limit of metastability of the liquid under negative pressures. The translational order parameter t has contours that are approximately parallel to the surface of minimum t. The values decrease toward this surface and increase beyond it. The behavior of these three measures of order shows that at different state points the liquid has different types of order, orientational and translational order in the W region and only translational order in the L region. As the space is traversed from the W region to the L region, there is a breakdown of the orientational order and a decrease in translational order. It is this breakdown of structure that gives rise to anomalies. In water, this breakdown gives rise to an anomalous density dependence of the entropy and of the translational and reorientational diffusion rates; in the family of liquids studied here, the diffusion rates exhibit a nonmonotonic dependence on γ due to structural breakdown. Though it is not surprising that the orientational structure breaks down as density or γ are increased, it is unexpected that the translational structure breaks down. The radial distribution functions shown in Figure 3 show that the tetrahedral and Lennard-Jones liquids have incompatible structures in the region of the second shell and, as one moves from region W to region L, the structure outside the second shell disappears and t goes through minimum values. It is this incompatibility of the shell structures allied to the change in the local orientational structure that gives rise to the unusual properties of these liquids. 4.2. Dynamics. For a normal liquid dynamical processes slow as the density is increased and speed up as the temperature is increased. However, the breakdown of local structure in the liquids studied here gives anomalous behavior.The decrease of tetrahedral structure with γ and density allows the molecules to reorient more freely, so that the reorientational rates generally increase with density at constant temperature. However, a point is reached at high densities where the reorientational rate starts to decrease again as a function of density, although the order parameter q continues to decrease. The translational diffusion constant D initially increases as γ increases, but eventually reaches a maximum, indicating the slowing of diffusion on progressive Lennard-Jones ordering. The diffusivity shows both minima and maximum with respect to density, the maximum being near the boundary between the regions W and L. The correlation of motion at short times with the instantaneous values of the local order parameter qi shown in Figure 11 suggests that the movement of molecules is closely related to the instantaneous local order, so that molecules move most when the surroundings have less structure but become trapped as the structure increases. However, the time scale for such structural changes is short. 4.3. Conclusions. In this work we investigated the effects of varying the relative weights of the Lennard-Jones and Coulombic contributions to the energy of the SPC/E model on its waterlike structural, dynamic, and thermodynamic anomalies. The SPC/E model is “understructured” with respect to real water, and hence its anomalies occur at lower temperatures. Our analysis suggests that a model with γ < 1 (a parameter range not explored in this work) would display anomalies at higher temperatures. In fact ab initio calculations with density functionals suggest that the true molecular dipole moment of liquid water is larger than in most classical water models.13,14 As the weight of the Lennard-Jones component is increased, the liquid evolves from tetrahedral to spherically symmetric.

6534 J. Phys. Chem. B, Vol. 109, No. 14, 2005 We found that important structural elements characteristic of these two types of order are incompatible. For example, the second solvation shell of water is located much closer to the central molecule than for the Lennard-Jones fluid. In fact, the radial location of the second peak in the oxygen-oxygen radial distribution function (tetrahedral liquid) coincides roughly with the position of the first minimum of the spherically symmetric Lennard-Jones liquid. This incompatibility results in an intermediate regime (sampled upon varying γ at fixed F or vice versa) of minimum translational order, where the fluid is less structured than at either the tetrahedral or spherically symmetric extremes. Barker et al.15 have compared the structure factors and void distributions in SPC/E water at densities of 0.94, 1, and 1.17 g cm-3 and show that breakdown of tetrahedral structure on compression can be associated with changes in the first (or pre-) peak of the void-void and the oxygen-oxygen S(q). Spherically symmetric models with oscillations, changes of curvature, or steps in their pair potential can also exhibit water-like anomalies (see, for example, ref 16 and references therein). It would be interesting to investigate the behavior of such models upon varying their interaction potential systematically between the Lennard-Jones and “core-softened”17 limits. It is not clear how far the anomalous properties seen in tetrahedral fluids depend on the fact that the local coordination is tetrahedral. Our results suggest that, for a liquid to show negative thermal expansion and anomalous density dependence in its transport properties, it must have a local structure that differs from and is inconsistent with simple Lennard-Jones type spherical symmetry.18 This may be possible with other types of local orientational structure. Acknowledgment. We thank J. Errington for supplying data. R.M.L.B. thanks the Department of Chemical Engineering at

Lynden-Bell and Debenedetti Princeton for hospitality and the Leverhulme Trust for support. P.G.D. gratefully acknowledges the support of the Department of Energy, Division of Chemical Sciences, Geosciences and Biosciences, Office of Basic Energy Science (grant DE-FG0287ER13714). References and Notes (1) Prielmeier, F. X.; Zang, E. W.; Speedy, R. J.; Lu¨demann, H.-D. Phys. ReV. Lett. 1987 59, 1128. (2) Angell, C. A.; Finch, E. D.; Woolf, L. A.; Bach, P. J. Chem. Phys. 1976, 65, 3063. (3) Errington, J. R.; Debenedetti, P. G. Nature 2001, 409, 318. (4) Bergman, D. L.; Lynden-Bell, R. M. Mol. Phys. 2001, 99, 1011. (5) Berendsen, H. J. C.; Grigera, R. J.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (6) Harrington, S.; Poole, P. H.; Sciortino, F.; Stanley, H. E. J. Chem. Phys. 1997, 107, 7443. (7) Stillinger, F. H.; Rahman. A. J. Chem. Phys. 1974, 60, 1545. (8) Scala, A.; Starr, F. W.; La Nave, E.; Sciortino, F.; Stanley, H. E. Nature 2000, 406, 166. (9) Ruocco, G.; Sampoli, M.; Torcini, A.; Vallauri, R. J. Chem. Phys. 1993, 99, 8095. (10) Smith, W.; Forester, T. R. 1996, The DL_POLY manual, Daresbury Laboratory. http://www.cse.clrc.ac.uk/msi/software/DL_POLY. (11) Chau, P.-L.; Hardwick, A. J. Mol. Phys. 1998, 93, 511. (12) Truskett, T. M.; Torquato, S.; Debenedetti, P. G. Phys ReV E 2000, 62, 993. (13) Silvestrelli, P.; Parrinello, M. Phys. ReV. Lett. 1999, 82, 3308. (14) Delle Site, L.; Alavi, A.; Lynden-Bell, R. M. Mol. Phys. 1999, 96, 1683. (15) Barker D. R.; Wilson, M.; Madden, P. A.; Medvedev, M. N.; Geiger, A. Phys. ReV. E 2000, 62, 1427. (16) G. Malescio, G. Franzese, G. Pellicane, A. Skibinsky, S. V. Buldyrev; Stanley, H. E. J. Phys.: Condens. Matter 2002, 14, 2193. (17) Debenedetti, P. G.; Raghavan, V. S.; Borick, S. S. J. Phys. Chem. 1991, 95, 4540. (18) Tanaka, H. J. Chem. Phys. 2000, 112, 799.

Computational Investigation of Order, Structure, and Dynamics in

Recommend Documents