J. Phys. Chem. 1981, 85, 1108-1119
1108
FEATURE ARTICLE A Random Network Model for Water Stuart A. Rice" The Department of Chemistry and The James Franck Institute, The University of Chlcago, Chicago, Illinois 60637
and Mark 0. Scealst The Department of Chemistry and Institute of Optics, The University of Rochester, Rochester, New York 14627 (Received: October 30, 1980; In Final Form: January 13, 1981)
A random network model for liquid water is discussed. This model, which uses information derived from studies of crystalline ices, amorphous solid water, and liquid water, is based on two assumptions, one of which refers to structure and the other to the dynamics of the molecular motion. These assumptions are (i) there is a continuous,albeit distorted, network of hydrogen bonds in the liquid, and (ii) it is meaningful to describe the molecular motion by using two distinct time scales. It is shown how adoption of these assumptions leads to the definition of a random network potential, and to rather good descriptionsof the 00 distribution function, the dielectric constant, and the bulk thermodynamic properties of liquid water. Despite its successes, the random network model fails to describe adequately some of the properties of liquid water; these discrepancies are discussed from the point of view of modifications needed to improve the model.
I. Introduction This article describes recent efforts to develop a simple, accurate, description of the structural, spectroscopic, and thermodynamic properties of liquid water. It is our assertion that such a description can be based on a random network model (RNM), in which there is essentially continuous hydrogen bonding throughout the configurations occupied in the liquid state. By a random network we mean a spatial structure which has the following features on a time scale which is short with respect to diffusion processes (- lo-'' s) yet long with respect to the characs): teristic periods of intermolecular vibrations ( (a) The distribution of nearest-neighbor separations is very narrow and is centered at a value close to that in a crystalline phase of the substance. (b) The bonding from molecule to molecule, which defines the connectivity of the network, is irregular, so that there are a variety of even- and odd-membered rings of molecules randomly interwoven throughout the network. At least some of the rings do not appear in crystalline phases of the substance. (c) The distribution of values of the angle between nearest-neighbor bonds has a nonzero, nontrivial, width. The notion that hydrogen bonding in water is continuous, although distorted, is not new. The essential characteristics of this model were described by Bernal and Fowler1 in 1933, while the use of the random network concept to represent the structure of a disordered system was introduced by Zachariasen2 in 1932. These ideas were combined by P ~ p l ein , ~1951, in what is usually called the bent hydrogen bond model. At that time the quantitative considerations presented were limited to a fitting of the model to the observed oxygen-oxygen pair correlation function and a calculation of the static dielectric constant N
'Alfred P. Sloan Fellow. Currently at the Department of Physical Chemistry, University of Sydney, Sidney, Australia. 0022-3654/81/2085-1108$01.25/0
of the liquid. Although the new developments to be described share with the Pople model the notion that a continuous distorted hydrogen bond network exists, our analysis of the consequences of that notion, and of the properties of the liquid, is fundamentally different from his. The most important change involves exploitation of another old idea, namely, that there exists a separation of time scales for various classes of molecular motion in liquid water.4 In particular, we assert that it is useful to regard libration and hindered translation as occurring in a quasistatic network of hydrogen-bonded molecules, and diffusion and related relaxation processes as occurring much more slowly. We also introduce a random network potential (RNP), the properties of which are inferred from spectroscopic and other data. The RNP incorporates the structural notions inherent to a distorted locally tetrahedral hydrogen bond geometry, yet depends only on the average oxygen-oxygen separation and the deviation of the hydrogen bond from linearity. Finally, our version of the RNM incorporates those aspects of the quantum dynamics of the system which are important even in the room temperature liquid. These are, primarily, the dependence of the molecular zero-point energy on the hydrogen bonding, the alteration of intramolecular interactions induced on condensation to the liquid or solid phases, and the proper weighting of thermal motions with frequencies from tens to many hundreds of cm-'. What can reasonably be expected of a theory of liquid water? A fundamental theory would start with only a given water-water interaction and predict all of the properties of the liquid in a specified thermodynamic state. Despite the outstanding contributions to our under(1)J. D. Bernal and R. H. Fowler, J . Chem. Phys., 1, 515 (1933). (2) W. H. Zachariasen, J. Am. Chem. Soc., 54, 3841 (1932). (3) J. A. Pople, Proc. R. SOC.London, Ser. A, 205, 163 (1951). (4) D. Eisenberg and W. Kauzmann, "The Structure and Properties of Water", Oxford University Press, Oxford, 1969.
0 1981 American Chemical Society
Feature Article
standing derived from computer simulations of water,5 a fundamental theory of the type referred to has not yet been developed. The random network model of water is an intermediate stage theory. It starts with the two assumptions that the dominant structural feature of the liquid is a locally tetrahedral continuous hydrogen-bonded network and that molecular motion in that network occurs, characteristically, on two very different time scales. The liquid structure is described, primarily, in terms of two distributions, namely, that of the separation of molecular centers and that of the deviation of the hydrogen bonding from linearity. The effective water-water interaction, the high-frequencymolecular motions, etc., are then expressed as functions of the rms widths of these distributions, and the effective molecular parameters used to calculate the properties of the liquid. In a model theory such as this one there are two classes of approximation. As an example of the first class we note that, obviously, the two assumptions concerning, respectively, structure and time scales for molecular motion, oversimplify the physical situation. Assumptions in this class are essential to the model and define its range of validity. Second, the use of these assumptions in practical calculations sometimes requires further simplifications not inherent to the random network model, for example, in the representation of the vibrational density of states of the liquid. At this stage of development of the theory of liquid water we believe it is reasonable to ask for modest accuracy in the prediction of thermodynamic properties, to seek organizing concepts which unify the description of all the condensed phases of water, and to be provided with clues to the interpretation of heretofore puzzling observations. As will be shown, the random network model of water fulfills these goals.
11. Low-Density Amorphous Solid Water and Ice Ih The interaction between two water molecules is strong (- lOkBT a t room temperature) and highly directional; it is responsible for a complex many-body potential energy surface which supports no fewer than 10 polymorphs of solid water, some of which contain unit cells with as many as 28 molecules. As background for our description of liquid water we consider first some properties of amorphous solid water and the ordinary form of ice which is stable near 1 atm and for temperatures below 273 K, namely, ice Ih. It is now well-known that very slow condensation of water vapor onto a variety of substrates at sufficiently low temperature leads to formation of an amorphous solids6 In the vast majority of cases the material so formed is the same irrespective of substrate composition and temperature, provided the latter is low enough. We call this material low-density HzO(as) for reasons that will be apparent soon. In one case, namely, deposition on an oriented single crystal of copper at 10 K, a different, high-density, form of H20(as)is made. Very little is known of the properties of high-density H,O(as). X-ray' and neutron diffraction8 studies of low-density H,O(as) establish the existence of tetrahedral coordination with nearest-neighbor separation of 2.76 A, very nearly the same as in ice Ih (2.751 A) (see Figure 1). The distribution (5)F. H. Stillinger, Adu. Chem. Phys., 31, 1 (1975). (6) E. F. Burton and W. F. Oliver, h o c . R. SOC.London, Ser. A , ,153, 166 (1935);P. Bondot, C. R. Acad. Sci. Paris,Ser. B , 265,316 (1967):268, 933 (1969); D. Olander and S. A. Rice, Proc. Natl. Acad. Sci. U.S.A., 69, 98 (1972). (7) C. G. Venkatesh, S. A. Rice, and A. H. Narten, Science, 186, 927 (1974); A. H. Narten, C. G. Venkatesh, and S. A. Rice, J. Chem. Phvs.. " , 64, 1106 (1976). (8)J. Wenzel, C. Linderstrom-Land, and S. A. Rice, Science, 187,428 (1975); J. Wenzel, Ph.D. Dissertation, University of Chicago, 1976.
The Journal of Physjcal Chemistry, Vol. 85,No. 9, 198 1 1109
-
-
2
1
0
cs,
0
2
6
4
R
(a)
8
IO
Flgure 1. The oxygen-oxygen pair correlation functions of polycrystalline ice Ih, H,O(as), and liquid water.
of nearest neighbors is quite narrow and is not the source of disorder in the solid. That disorder is generated by a distribution of 000 angles about 109.5', the value characteristic of an ice Ih lattice. A distribution of 000 angles with width of the order of 8' is inferred from the experimental data, and this is sufficient to destroy correlation between the orientations of molecules separated by more than 7 A. The density of this form of H,O(as), both as measured directlyg and as inferred from fitting the X-ray diffraction data,' is 0.94 g ~ r n -to ~ ,be compared with the density of ice Ih,lo namely, 0.93 g ~ m - ~ . The differences between low-density HzO(as)and ice Ih are analogous to the differences between amorphous Si and Ge and their parent crysta1s.l' In each case the parent crystal has the diamond structure (for ice Ih this refers to the location of the 0 atoms), the amorphous phase retains local tetrahedral coordination, the differences in nearestneighbor separation and in density between the crystalline and amorphous solids are of the order of 1%,and the disorder in the amorphous solid is generated by a distribution of angles about the tetrahedral value, the width of that distribution being of the order of 7-8". A variety of model studies12suggest that amorphous Si and Ge have a fully bonded random network structure; we adopt the (9) J. A. Ghormley, Science, 171, 62 (1971). (IO) P. V. Hobbs, "Ice Physics", Oxford University Press, London, 1974. (11) See W. Paul, R. Temkin, and G. N. Connell, Adu. Phys., 22,529, 581, 643 (1973). (12) D. E. Polk, J. Non-Cryst. Solids,5, 365 (1971); D. E. Polk and D. S. Boudreaux, Phys. Reu. Lett., 31, 92 (1973); M. G. Duffy, D. S. Boudreaux, and D. E. Polk, J. Non-Cryst. Solids, 15, 435 (1974); D. Henderson, ibid., 16, 317 (1974).
1110
The Journal of Physical Chemistry, Vol. 85,
AMORPHOUS SOLID
No. 9,
Rice and Sceats
1981
H20
DEPOSITION AT IO’K
i, POLYCRYSTALLINE
H20
DEPOSITION AT 150’K
c V” (cm-’) Flgure 2. The OH-stretching regions of the Raman spectra of polycrystalline ice Ih, H,O(as), and liquid water.
same model for low density H20(as). The higher density form of H20(as)also has local tetrahedral coordination, a nearest-neighbor separation of 2.76 A with very little dispersion, and disorder generated by a distribution of 000 angles.’ In addition, there is a next-nearest-neighbor peak at 3.3 A, and a third peak at 4.5 A, whereas in ice Ih and low-density H20(as)the second neighbor peak is at 4.5 A. A consequence of the difference in structure is a higher density; it is inferred to be p = 1.1 g from a fitting to the X-ray diffraction data. No other data are available for this form of HzO(as), so speculation about its nature is poorly founded. Given the apparent structural parallelism between the pairs of solids ice Ih and low-density H20(as), and crystalline and amorphous Si, Rice13 suggested that high-density H20(as) is derived from ice I1 or ice I11 by introduction of a distribution of 000 angles and a distribution of irregular rings analogous to that in HzO(as). Ices I1 and I11 have more complicated crystal structures than ice Ih. Each has 12 molecules per unit cell, a second nearest neighbor at about the right distance, and a density not far from that observed.1° Moreover, in each of these higher ices there is a large range of 000 angles, a feature which must also exist in the high-density HzO(as);it is the bending of the 000 angles which generates the higher density and the different kinds of second near neighbors (compared to ice Ih) in the single completely hydrogen-bondednetworks of ices I1 and III.14 In view of the distortions that lead to (13)S.A. Rice, Top. Curr. Chem., 60, 109 (1975). (14)B. Kamb, Acta Crystallogr., 17,1437(1964);B.Kamb and S. K. Datta, Nature (London),187,140(1960);B.Kamb, W.C. Hamilton, S. J. La Placa, and A. Prakash, J. Chem. Phys., 65,1934 (1971).
1
0
20
40
60
80
100
120
140
160
T (K) Flgure 3. The temperature dependence of the frequency of OH stretching. Data for all eight posslbilitles (OH and OD stretches In polycrystalline ice Ih and H,O(as), both fully coupled and mechanlcally decoupled) fall on one curve. The lower panel compares the observed shift in frequency wlth that expected because of thermal expansion. The squares, representing the calculated values, are shown without error bars. I f the uncertainty in the coefficient of expansion is displayed as a set of error bars, the two sets of points overlap.
loss of orientational and positional ordering in the amorphous solid, fine distinctions between the possible parents of a random network structure will be meaningful only if they lead to predictions of observable differences associated with that parentage; for the moment such distinctions have not been tested. For the remainder of this article we will discuss only the properties of low-density H20(as),hence will drop the term “low density”. The vibrational spectra of the OH-stretching regions of ice Ih and H20(as)have been extensively studied.I5 As shown in Figure 2, even a cursory examination of the spectra reveals striking similarities, suggesting that the origins of particular features are likely to be the same in the two condensed phases. The stretching spectrum of strongly coupled OH oscillators spans several hundred cm-l. In contrast, the OH-stretching spectrum of a mechanically decoupled OH oscillator, in DzO ice Ih or DzO(as), is a single narrow peak. (The decoupling is a consequence of the OH/OD frequency mismatch.) It is interesting that the shifts with temperature of the decoupled OH-stretching frequency and of the lowest peak in the coupled OH-stretching spectrum are the same, but the width of the decoupled oscillator transition is independent of temperature, whereas that of the peak in the coupled oscillator spectrum mentioned is not. Also, the temperature dependences of the widths of the lowest peaks in the Raman spectra of the coupled OH stretches in ice Ih and H20(as) are the same, though the limiting values at zero temperature are different. A t the phenomenological level the shift and broadening of the low-frequency peak in the OH-stretching spectrum (15)T.C. Sivakumar, S. A. Rice, and M. G. Sceats, J. Chem. Phys., 69,3468 (1978).
Feature Article
The Journal of Physical Chemist& Vol. 85, No. 9, 198 1 1111
TABLE I: Calculated and Observed Overtones of HDO Transitions in Ice Ih overtones, cm-' transition -
, /T
I
I
(000)- ( 0 2 0 )
( 0 4 0j
(101) (021) (002)
0
obsd 2875 3845
5600 5720, 5 7 6 0 6100 6330
calcd 2870 3852 4701 4690 5258 5621 5678 6108 6290
of ice Ih can be accounted for as follows. It is well established that the OH-stretching frequency in a hydrogen-bonded system depends on the average donor-acceptor separation.16 For water Falkl' has found that dijoH/dR = 1843 cm-l/A when R = 2.76 A. When this correlation is used the temperature dependence of the peak shift is accurately reproduced by allowing for the indirect influence of thermal expansion (Figure 3). Moreover, the difference in low-frequency peak positions between the ice Ih and H,O(as) spectra agrees with the shift predicted from the measured difference in nearest-neighbor separations. A second element of the phenomenological analysis is the recognition that proton disorder, and even protondeuteron mass disorder, does not destroy the low-frequency lattice vibrations of ice Ih.16 Given that the density of vibrational states for OH stretching is high only for frequencies equal to or greater than the low-frequencypeak of the Raman spectrum,19the temperature dependence of the width of that peak is accounted for by phonon scattering, the frequency of the dominant lattice phonon being 200 cm-' (Figure 4). This is consistent with other observations, since there is a peak in the density of lattice vibrational states of ice Ih at 200 cm-'. Finally, the temperature-independent width of the decoupled oscillator transition, and the zero temperature width of the coupled oscillator transition, are associated with a distribution of local environments generated by proton disorder.% That is, it is suggested that the differing force fields resulting from proton disorder lead to a small dispersion of 00 separations; the spectroscopic data, coupled with the stretching frequency-separation correlation mentioned, lead to the inference that the fwhm of that distribution of separations is 0.014 A. By virtue of the similarities in behavior of the spectral features of the OH stretch in ice Ih and HzO(as),the same phenomenological analysis applies to both solids. The central problem in the microscopic interpretation of the OH-stretching spectrum is the determination of the extent to which the free molecule spectrum is altered by interaction between OH oscillators on different molecules.
A straightforward description of the molecular dynamics of OH stretching in ice Ih and H20(as)can be based on a spectroscopic model, that is one that fits a force field to the vibrational spectrum. In the simplest of these models the OH motion is taken to be harmonic, using an effective force constant that subsumes any effect of anharmonicity and the frequency shift (relative to a free oscillator) due to hydrogen bonding; any interaction between OH stretching and HOH bending motions, or OH stretching and 0-0 lattice motions, is neglected. Surprisingly, this simple model fits the observed OH-stretching spectrum reasonably well, but not perfectly.l92l An improved model, which includes anharmonic terms in the potential energy and the interaction between HOH bending overtones and the OH-stretchingmotions, fits the observed spectrum very well (see Figure 5).22 In addition to the expected importance of the intermolecular coupling, the following strong inferences follow from the analysis of the OH-stretching region of the spectrum: (i) Hydrogen bonding in the condensed phase alters some of the intramolecular interactions. For example, the intramolecular OH oscillator stretch-stretch interaction is changed on passage from the free molecule to the solid. In the free molecule the intramolecular coupling constant is k+ = -0.101 mdyn/A, but in the solid it is 0.06 mdyn/A, of opposite sign. The reality of this sign reversal is s u p ported by the Bertie and Bates analysis of the spectra of proton ordered ices I1 and IX,23where it is also found to occur, and by the existence of a correlation between the value of the stretching force constant k, and the stretchstretch interaction force constant krf in symmetrically bonded which correlation predicts a value for kr,,in substantial agreement with that cited as characteristic of ice Ih. (ii) It is a good approximation to assume that the anharmonic coefficients are unchanged in going from the free molecule to the condensed phase. This implies that the change in anharmonic contribution to the vibrational energy induced by hydrogen bonding is traceable to the large shift in the diagonal OH-stretchingforce constant, k,. The dependence of k, on the 00 separation is known from an internally consistent analysis of the OH-stretching spectra of a wide variety of systems.17 A convenient vehicle for testing this assumption concerning the anharmonicity of the potential surface is the analysis of the overtone spectrum of ice Ih.25 Sceats and Rice have made such an
(16)M. Falk and 0. Knop, Water: Compr. Treat., 2, 55 (1973). (17) M.S.Falk, "Proceedingsof the Electrochemical Society Symposium on Chemistry and Physics of Aqueous Solutions", Toronto, 1975, p 79. (18)J. C. Bellows and P. N. Prasad, J. Chem. Phys., 64,3674(1976). (19)R. McGraw, W.G. Madden, M. S. Bergren, S. A. Rice, and M. G. Sceats, J. Chem. Phys., 69,3483 (1978). (20)J. E.Bertie and E. Whalley, J. Chem. Phys., 69, 3497 (1978).
(21)W. Madden, M. Bergren, R. McGraw, S. A. Rice, and M. G. Sceats, J. Chem. Phys., 69,3497 (1978). (22) M. Bergren and S. A. Rice, to be submitted for publication. (23)J. E. Bertie and F. E. Bates, J. Chem. Phys., 67, 1511 (1976). (24)K. Kuchitsu and Y. Moreno, Bull. Chem. SOC.Jpn., 38, 814 (1965);J. Schiffer, M.Intenzo, P. Hayward, and C. Calabrese, J. Chem. Phys., 64,3014 (1976);R. A. Fifer and J. Schiffer, ibid., 52,2664 (1970); S. C.Mohr, W. D. Wilk, and G. M. Barrow, J. Am. Chem. SOC.,87,3048 (1965).
Flgure 4. Observed (0)and calculated (-) temperature dependence of the width of the 3100-cm-' peak in the Raman spectrum of polycrystalline ice Ih. The theoretical curve, which is fitted to the datum at 150 K,corresponds to an average effective lattice frequency of 200 cm-l.
Rice and Sceats I
s-
I
!
(Dru
m
O 0 '
t
l
I
I
I
l
l
n
A
& A 3000.
l
+
1
1
EXPERIMENT
=
e3200. 7 -
FREQUENCY
h
-
1
3400.
1
1
I
-
,,I
3600.
ICM-')
"LU
zd
3000.
3200.
FREQUENCT
3400.
3600.
ICM-'l
Figure 5. Comparison of the observed and predicted OH-stretching regions of the Raman and infrared spectra of polycrystalline ice Ih. A and C: H,O. Band D: D,O.
analysis;26the predicted and observed overtone spectra are in good agreement (see Table I). (iii) Just as the intramolecular stretch-stretch interaction constant, k r f , is altered by hydrogen bonding, so also is the coupling constant which determines Fermi resonance, kraa. It appears that k,,, = 0 in ice Ih, whereas it has the value 0.15 mdyn/A in the free molecule. Inferences (i), (ii), and (iii) can be checked by analysis of the observation, by Ritzhaupt and D e ~ l i n , ~of' the (25)C. Haas and D. F. Hornig, J. Chem. Phys., 32,1763(1960);W.A. P.Luck and W. Ditter, 2. Naturforsch. B , 24,482 (1969);J. D. Worley
and I. M. Klotz, J. Chem. Phys., 45,2868(1966);D. Kroh and A. Ron, Chem. Phys. Lett., 36, 527 (1975). (26)M.G. Sceats and S. A. Rice, J. Chem. Phys., 71, 973 (1979). (27)G.Ritzhaupt and J. P. Devlin, J. Chem. Phys., 67, 4779 (1977).
spectrum of intact D20 isolated in H20 ice Ih and H20(as). The splitting in this spectrum between the symmetric and antisymmetric stretches of D20 appears, superficially, rather like that in the free molecule, for which k,, < 0, so it is tempting to draw that conclusion. However, the observed spectrum is very much broader than expected from the observed widths of the isolated OD oscillator transitions in H20 ice Ih and H20(as).15 Sceats, Stavola, and Rice2$have shown that a more consistent interpretation, which reproduces both the splittings and widths in the spectra of D20 in both H 2 0 ice Ih and H20(as),follows from recognizing that there is an inhomogeneous broad(28)M.G.Sceata, M. Stavola, and s. A. Rice, J. Chem. Phys., 71,983 (1979).
The Journal of Physlcal Chemlstty, Vol. 85, No. 9, 1981 1113
Feature Article
I
D20 in H2O (as)
--
!300
2400
k raa
100 cm-'
1-94,+ cm-
2500
F R EQ u E N CY
( c rn-' )
Figure 6. Comparison of the observed and predicted absorption spectra of intact DO , in polycrystalline H,O ice Ih and in H,O(as).
ening of the bending mode transition in both solid media. This inhomogeneous distribution leads, in turn, to a distribution of splittings and intensities in the ensemble of intact H 2 0 or DzO molecules that constitute the sample studied. The analysis used by Sceats, Stavola, and Rice adopts inferences (i) and (ii) and, from the fitting of the observed spectrum (see Figure 6 ) , deduces (iii). Given the above description of the OH-stretching spectrum of ice Ih, we must expect that the OH-stretching spectrum of H20(as)will differ from that of ice Ih because of the spread in OH-stretching force constants generated by the distribution of local environments; Madden et invert the stretching spectrum of the isolated OH oscillator in D20(as)to obtain the distribution of stretching force constants needed as input for the model calculations. Otherwise, as in the case of ice Ih, the features of the OH-stretching region of the spectrum are dominated by the effects of intermolecular coupling of the oscillators. Calculations have, to date, only been carried out by using the cruder quasiharmonic model mentioned first above. The quality of the agreement between the predicted and observed spectra is the same as for ice Ih at the same level of approximation. In essence, these calculations support Whalley's contentionz9that the features of the vibrational spectrum of ice Ih, and by extension of HzO(as), arise from strong intermolecular coupling. Simple interpretations of the spectrum in terms of the molecular modes of the free molecule are not very useful since the consequences of the intermolecular coupling are quite complex. The fact that the density of liquid water is not very different from the densities of ice Ih and H20(as)suggests there cannot be a radical change in the importance of intermolecular coupling as one goes from, say, ice Ih to the liquid, and the gross similarities between the vibrational spectra, dielectric constant, and other properties of these condensed media supports this notion (see Figures 1and 2). Although these studies of H,O(as) and ice Ih are incomplete, and the models used to interpret the experimental data are approximate to varying degrees, considerable information of use in the construction of a model of liquid water can be culled. It is clear that the distribution of 000 angles plays an important role in deter(29) E. Whalley, Dev. Appl. Spectrosc., 6 , 277 (1968).
mining both the structural and dynamical properties of condensed phases of water, and that some key properties of the water-water interaction (e.g., the anharmonicity of the OH-stretching motion) are the same in the free molecule and in the condensed phases, while others are very different (e.g., the intramolecular coupling constant ke). These deductions serve as guideposts in the development of the random network model, and the accompanying random network potential, for water. The construction of the random network potential is particularly important, since it provides the means for converting the structural and dynamical implications of the random network model into predictions of the properties of liquid water. 111. A Random Network Potential Sceats and Rice30 have developed an effective potential to represent the interaction between hydrogen-bonded water molecules. By virtue of the methodology used that potential can only be valid in a small domain of configuration space, correspondingto oxygen-oxygen separations from 2.7 to 3.0 A, but it is this domain which is most important in the determination of the properties of the condensed phases of water. The Sceats-Rice random network potential is based on five ideas: (i) Because the region of configuration space sampled in a condensed phase is small, the effective water-water interaction can be generated as a function of displacement from a reference configuration via a power series expansion; the coefficients in the expansion can be evaluated with the aid of spectroscopic data. (ii) It is asserted that only two coordinates are needed to define an accurate effective potential, namely, the 00 separation and the hydrogen bond angle, a simplification which is supported by a variety of subsidiary model calculations. (iii) It is assumed that the diagonal anharmonicity in the potential surface is independent of the strength of the hydrogen bonding (note that this is the intermolecular potential energy surface), so that the cubic term in the expansion of the potential energy can be extracted from the Gruneisen constant of the solid.1° The three assumptions we have described define the variation in the potential relative to a hydrogen-bonded reference configuration, which is taken to be the tetrahedrally coordinated RNM described for H20(as). To complete the effective potential account must be taken of the interaction between molecules which are not hydrogen bonded, and also of nonadditivity of the pair interactions. (iv) It is assumed that the nonbounded interactions between water molecules are accurately described by a combination of the Kamb potential (developed originally for the interaction of interpenetrating but disjoint hydrogen-bonded sublattices in ice VI131) and the electrostatic interaction of a dipolar molecule and its surroundings. (v) The nonadditive interactions so evident in the quantum mechanical calculations of the energies of water polymers32are assumed to be similar in the liquid and in ice Ih, since in the RNM the average local environment of each water molecule has tetrahedral coordination geometry. An accurate estimate of this nonadditive component of the energy is obtained from the difference between the experimental value of the lattice energy of ice (30) M. G. Sceats and S. A. Rice, J. Chem. Phys., 72, 3236 (1980). (31) B.Kamb, J. Chem. Phys., 43, 3917 (1965). (32) H. J. C. Berendsen and G. A. vander Velde in CECAM Report of Workshop on Molecular Dynamics and Monte Carlo Calculations on Water, 19 June-11 August, 1972, p 77; J. E. Del Bene and J. A. Pople, J. Chem. Phys, 52, 4858 (1970); 58, 3605 (1973).
1114
The Journal of Physical Chemistry, Vol. 85, No. 9, 1981
TABLE 11: Calculated and Observed Librational and Hindered Translational Frequencies in Several Ices at Low Temperature
e 2 lI2,
@L), cm-'
I
ice I ice11 ice I11 iceV
(Fs),cm"
deg
a, A
obsd
calcd
obsd
calcd
0.0 12.0 11.6 13.0
2.571 2.80 2.78 2.80
810 694 743 704
810 664 673 642
224 151 172 169
225 185 187 178
Ih33and that computed from the first four interaction terms mentioned above; it is found to have the form ENA= -18.41 110.46(AR/2.751) kJ mol-'
Rice and Sceats
TABLE 111: Distribution of Hydrogen Bonds no. of bonds1 molecule
configurations, % time averaged instantaneousa
0.00 0.77 6.33 25.62 67.28
0 1 2 3 4
-5
0 2 11 36 50 2 3.4
3.6
~ H B
a Reference 1 3 .
+
with R in angstroms. As applied to the liquid, these contributions to the energy lead to the following effective potential:
I
E,dAR, e,? = E s ( W + EB(es2) + EEL+ ENB with E,(hR) = -52.37 - 22.68(AR/2.751) + 919.6(AR/2.751)2 - 1195.5(AR/2.751)3 k J mol-' EB(es2) =
16.18[1 - exp(-96,2)] kJ mol-'
ENB+ EEL= -6.09
- 0.0043(T -
114) kJ mol-]
The random network reference configuration is characterized by two parameters. These are, first, the distribution of hydrogen bond angles in the quasistatic network, measured by the mean square width Os2 = 211.6 + 1.33(T - 273) deg2 which is determined30from the measured librational frequency in the liquid, the functional dependence of the librational frequency on predicted by theory, and the measured value of for H , O ( ~ S )and, , ~ second, the mean molecular separation R = 2.835 + 0.00052(T - 273) A determined from diffraction studies. The effective potential is different for the crystalline ices by virtue of variations in the value of Em. It should be noted that the true bending potential is asymmetric since it matters whether the direction of bending is such as to lead to another H or a lone pair being inserted between the oxygens. The Sceats-Rice RNP neglects this asymmetry. One consequence of this simplification is a gross overestimate of the effective anharmonicity of the librational motion, to the extent that for large amplitude displacements the RNP does not yield accurate librational frequencies. Nevertheless, the functional dependence of the librational frequency on 8," is correctly predicted. This is an important point, because the RNM of liquid water postulates the existence of a broad distribution of 000 angles, and it is necessary for the RNP to weight accurately the contributions to the energy from water molecules in configurations in which the 000 angles deviate markedly from 109.5O. The higher ices are a well-defined set of materials in which large deviations of the 000 angles from 109.5' occur. As shown in Table 11, the predicted and measured librational3* and translationaP5 frequencies of several of the higher ices are in reasonably good agreement. (33)E. Whalley, in "Physics and Chemistry of Ice", E. Whalley, S. J. Jones, and L. W. Gold, Ed., University of Toronto Press, Toronto, 1973, p 73. (34)J. E. Bertie and E. Whalley, J.Chem. Phys.,40,1637,1646(1964). (35)M.J. Taylor and E. Whalley, J. Chem. Phys., 40, 1660 (1964).
2.60
2.80
3.00 R
3.20
cir
Figure 7. A comparison of the 00 separation distribution function for hydrogen-bonded pairs in time-averaged configurations (TAC) with the predicted static distribution of nearest neighbors from the Sceats-Rice random network model (RNM) at 283 K.
IV. Comparison of Molecular Dynamics Simulation Data with RNM Predictions The static and dynamical assumptions of the SceatsRice RNM have an implicit connection. An instantaneous configuration in a molecular dynamics simulation, or in an equilibrated Monte Carlo simulation, captures the distributions of positions and orientations which include all the thermal excitations. In the RNM it is assumed, via the separation of time scales, that the configurations averaged over the librational and hindered translational motions will have a much narrower distribution than those including the thermal excitations of these motions. Moreover, the static random network is considered to be tetrahedrally coordinated, and moments of the 000 angle distribution can be predicted from the distribution of hydrogen bond angles. Belch, Rice, and S c e a t P have tested several of the characteristics of the RNM using time-averagedmolecular dynamics data from a simulation of water in which the molecules are represented as rigid bodies which interact via the ST-2 p ~ t e n t i a l .For ~ purely technical reasons, the time averaging of the simulation data could be carried out only over one and one-half librational periods, which is really too short to fully damp out the postulated oscillatory motion about the static, distorted, random network. Nevertheless, the following features of the RNM are confirmed: Analysis of the hydrogen bond connectivity, using a geometric criterion to define the hydrogen bond, shows the existence of a continuous network of hydrogen bonds without any small clusters of free water molecules. The average number of hydrogen bonds per molecule is 3.6 for this potential, density, and temperature. Table I11 displays the reduction in the fractions of 1, 2, and 3 hydrogenbonded molecules in the time-averaged configurations relative to the values found in a recent Monte Carlo sim~lation.~'Although different intermolecular potentials (36)A. Belch, S.A. Rice, and M. G. Sceats, submitted to Chem. Phys. Lett. (37)W. L.Jorgensen, Chem. Phys. Lett., 70,326 (1980).
Feature Article
The Journal of Physical Chemistry, Vol. 85, No. 9, 198 1
303 r
I
+ L
-
L
22
40
/
L
I
~
I
80 ICC 120 140 160 COO Apq e ir Degrees 60
I
I80
Flgure 8. A comparison of the 000 angle distribution functions in time-averageconfigurations (TAC) with the prediction of the SceatsRice random network model (RNM) at 283 K.
were used in the two simulations, we attribute the differences in hydrogen bonding to the change in perception generated by averaging over the librational motions of the molecules. Indeed, it is likely that averaging for a longer time will generate a set of configurations with more than 3.6 hydrogen bonds per molecule. The 00 separation distribution function for the timeaveraged molecular dynamics simulations is shown in Figure 7, along with that deduced by Sceats, Stavola, and Rice (see next section). As implied by the separation of time scales in the RNM, the static 00 distribution function is much narrower than that found from diffraction studies of water because the latter includes the broadening due to thermal motion. Given that the simulation results correspond to the ST2 potential, which is at best only a good approximation to the real water-water potential, and that Sceats, Stavola, and Rice use an empirical relationship to obtain the distribution plotted, quantitative agreement between the two distributions is not to be expected. In fact, they are surprisingly similar, and the widths and asymmetries are nearly the same. The 000 angle distribution deduced from analysis of the time-averaged molecular dynamics simulations is shown in Figure 8, along with the 000 angle distribution predicted from the RNM. The computer-generated distribution peaks at 107', is very nearly Gaussian with a full-width at half-height of 48', and is weakly skewed toward 180'. The RNM distribution peaks at 105.5', has a full-width at half-maximum of 42.3', and correctly predicts the skewing toward 180'.
V. Liquid Water We now illustrate how the inferences concerning molecular dynamics and local ordering, drawn from the studies of H20(as) and embodied in the RNM, can be used to describe liquid water. It must be emphasized that the ideas put forward in the following text describe only a zeroth order model of liquid water, hopefully one which includes most of its important characteristics and which will be robust under the refinements necessary to account for the details of particular properties of the liquid. A. T h e Radial Distribution Function of Liquid Water. In the RNM, liquid water has a structure with local tetrahedral coordination of the molecules, and broad distributions of 000 angles and 00 separations. In the version used by Pople3 the full widths of the observed peaks in the 00 pair correlation function were fitted to widths of peaks of overlapping distributions of first, second, and third neighbor molecules. There is, in this Pople model, no separation of the time scales for different classes of molecular motion. In contrast, Sceats, Stavola, and Rice38
1115
exploit the separability of slow and fast molecular motions. The structure is viewed by them as consisting of a distorted static hydrogen-bonded network with thermal excitation of librational and hindered translational motions. There is, in addition, an interaction between nonhydrogenbonded molecules, omitted in the Pople calculations,which Sceats, Stavola, and Rice include in their analysis of the RNM predictions. To determine the properties of the static continuous random network referred to above, Sceats, Stavola, and Rice invert the empirically established relationship between OH-stretching frequency and 00 separation, using as input the isotropic component of the Raman spectrum of the OH vibration^.^^ The nearest-neighbor distribution so determined is much narrower than is the first peak in the 00 correlation function (see Figure 7), in agreement with the idea that hydrogen-bonded water molecules oscillate about quasistatic positions in the random network for at least several periods of vibration, and that relaxation of the local structure that defines the quasistatic centers of oscillation is slow relative to vibration. Nonbonded interactions between water molecules influence the static distribution of neighbors. To include this effect Sceats, Stavola, and Rice use the already mentioned Kamb p ~ t e n t i a l .The ~ ~ average effect of the interaction between nonbonded molecules is displacement along the 00 directions and alteration of the hindered translation mode force constant, kRR, since that force constant is a function of the average 00 separation. Given that on the time scale of interest the dynamics of the continuous random network can be represented as a superposition of oscillatory motions, and assuming that these motions are quasiharmonic, calculation of the effects of thermal excitation is straightforward. The lattice dynamics of ice I in the 60-200-~m-~ region can be reasonably well modeled by the use of only two force constants, kRR for 00 stretching and k , for 000 bending. The modeling of the same spectral region for liquid water is more complicated by virtue of the much greater importance of anharmonic contributions to the potential energy (these determine the temperature dependence of the quasiharmonic frequencies), but not fundamentally different. It is a simple matter, then, to determine the thermal distribution of nearest-neighbor oxygens from the convolution of the amplitude distribution and the static distribution inferred from the isotropic component of the Raman spectrum. Proceeding outward from the central molecule, the amplitudes of motion of the second and third neighbor molecules are primarily determined by the 000 angle bending. The thermal distributions of the separations of these neighbors are calculated from the respective convolutions of the amplitude distributions with the static network distributions. The latter are generated from the properties of the nearest-neighbor static distribution and modified by the effect of interactions between molecules not directly hydrogen bonded. The predicted 00 pair correlation functions are compared with experimental data in Figure 9. The agreement between the RNM and experiment is clearly rather good, but not perfect. B. T h e Dielectric Constant of Liquid Water. The assertion that the librations and hindered translations of hydrogen-bonded molecules are rapid on the time scale of diffusion and network relaxation has implications for the interpretation of the static dielectric constant of liquid (38)M.G.Sceats, M. Stavola, and S. A. Rice, J. Chem. Phys., 70,3927 (1979). (39)J. R. Scherer, M. K. Go, and S. King, J.Phys. Chem., 78, 1304 (1974).
Rice and Sceats
TABLE V: Lattice Energies of Several Ices at 0 K Liquid Water
E , kJ mol-’
2 .o
ice I ice I1 ice I11 ice V
goo I .o
I
Liquid Water’ 50°C
I
Liquid Wder
c
2.0
20
40
60
Rill Figure 9, Comparison of the random network model prediction (-) and the observed (..--) oxygen-oxygen pair correlation function of liquid water at several temperatures.
TABLE IV: Calculated and Observed Values of the Dielectric Constant of Water €0
t, “C
obsd
calcd
10 90 200 300
83.80 58.32 34.50 19.57
86.46 59.41 35.34 21.27
water. For example, the high-frequencyoscillations about equilibrium positions do not contribute to the static dielectric constant eo; only displacementsproportional to the field strength contribute to eg. We expect, then, that only those dipolar correlations which are represented by the static distribution function of the random network contribute to e@ Sceats, Stavola, and Rice use this assumption and the Kirkwood-Onsager-Dupius40’41formulation to calculate e@ For ice Ih a direct multipole expansion and subsequent evaluation of the electric field strength at a molecule leads to the effective molecular dipole moment p = 2.60 D.42But p is dependent on molecular configuration, and so varies with the angular disorder in the static random network. Sceats, Stavola, and Rice suggest that p = 2.60 cos2 ( c ) l / z D . Although the necessity for this assumption concerning the functional form of the dependence of p on disorder limits the value of the computation of eo as a test of the RNM, the results are still of considerable interest. As shown in Table IV, there is excellent agreement between model predictions and measured values (40) J. G.Kirkwood, J. Chern. Phys., 7,911 (1939). (41) L. Onsager and M. Dupuis in “Electrolytes”, B. Proce, Ed., Pergamon Press, London, 1962, p 27. (42) C. A. Couleon and D. Eisenberg, Proc. R. SOC.London, Ser. A , 291, 445,454 (1966).
obsd
calcd
-47.3 -46.2 -46.3 -46.0
-45.05 -44.43 -44.36 -43.68
over a wide temperature range. If the dipole moment for ice Ih calculated by Eisenberg and C o u l s ~ n(p~=~ 2.60 D) were scaled downward even better agreement would be obtained. These results are also consistent with the dielectric constant of ice 111, which has a distribution of hydrogen bond angles with an rms width of 16.4’. The predicted value of eo is 126, the measured value 117,@both at 253 K. C. The Internal Energy. The separation of time scales in the RNM suggests that it is convenient also to separate the internal energy, entropy, etc. into vibrational and configurational components. Sceats and Rice do this, and then approximate the density of states of low-frequency vibrational modes of the RNM by three librational, two hindered translational, and one 000 flexing mode per molecule,43 An assessment of the accuracy of the approximate representation of the density of vibrational states can be obtained from a comparison of the computed and observed moments (P); it is found to be satisfactory.3o It must be remembered that these frequencies, and their parent RNP, are dependent on the distribution of hydrogen bond angles (to be specific on the mean-squared width of that distribution). Given these frequencies it is straightforward to calculate the vibrational component of the energy at temperature T. A second contribution to the vibrational energy arises from the molecular distortion induced by hydrogen bonding, which changes the zero-point energies of the intramolecular OH-stretching motions. Finally, the configurational contribution to the energy is calculated from the random network potential, the properties of which were described in section 111. The calculated@and measured energiesu of several ices are compared in Table V; the differences amount to less than 5%. The measured values of the energies of the ices, and the calculated values, increase monotonically with increasing hydrogen bond bending, but at a rate that is much less than that described by the bending contribution itself. This observation was first made by Kamb.“ In agreement with his comment, our analysis shows that the destabilization induced by increased hydrogen bond bending is partly overcome by a change in the contribution of nonbonded interactions to the energy. However, a still larger stabilization arises from the decrease in the zero-point energy with hydrogen bond bending, and the dominant contribution to this arises from the decrease in librational frequencies. This balancing implies that the molecular configurations of ices I-IV have, in part, similar energies by virtue of quantum effects. It is then reasonable to expect to find that quantum effects play a significant role in determining the properties of liquid water. A comparison of the calculated43and measured energies of liquid water is shown in Figure 10; the agreement is very good. Examination of the calculations shows that destabilization of the liquid relative to ice Ih at 273 K is dominated by the change in hydrogen bond bending, with a (43) M. G. Sceats and S. A. Rice, J. Chern. Phys., 72, 3248 (1980). (44) B. Kamb in ‘‘Structural Chemistry and Molecular Biology”, A. Rich and N. Davidson, Ed., Freeman, San Francisco, 1968, p 507.
Feature Article -25
I
The Journal of Physical Chemisity, Vol. 85, No. 9, 7981
1117
TABLE VI: Entropy of Liquid Water S, J K-' mol-'
water (273K ) water (373K)
obsd
calcd
63.46 86.97
63.72 86.99
The configurational entropy of the random network can be divided into three contributions. The proton disorder S, = R In (3/2), is taken to be the same for all configurations likely to be found in the liquid or H20(as). If we assume that the 00 rms zero-point amplitudes of the random network and ice Ih are the same, the quasistatic distribution of 00 -separations - - introduces a contribution SR = 2R In [(8R2 6RZp2)/8R~p2]1/2, since there are two hydrogen bonds per water molecule. The volume of space made available by distortion of the hydrogen bond angle is difficult to calculate because of the bonding constraints in the network. Sceats and Rice46argue that it is proportional to (@)1/2 because the three 0's defining the angle 4 always define a plane. Therefore they take the contribution to the entropy from distortion of the 000 angles in the network - to be S, = R In + @)/ 8&p2]1/2,with (8q5Zp2)1/2 the rms zero-point bending amplitude in ice Ih. The value of (?@)1/2 can be obtained by inverting the isotropic component of the Raman spectrum, as explained in section VA, and the value of (@?1/2 from the Sceats-Rice correlation between mean square hydrogen bond angle and temperature. The needed zero-point amplitudes are calculated from the same force constants and frequencies used to model the density of states of lowfrequency vibrations (see section VC). Now a perfect amorphous solid does not possess a configurational entropy because there is not the molecular mobility necessary to explore the possible configurations. Nevertheless, it is usual to associate a configurational entropy with the amorphous solid where it is understood that this is the entropy of the supercooled liquid at a temperature infinitesimally above that for which molecular mobility is frozen out. The Sceats-Rice estimate for S, for H20(as)is 6.05 J K-' mol-l, in very reasonable agreement with the upper bound of Bell and Dean,47namely, 5.76 J K-l mol-l, established for vitreous silica by consideration of the possible configurations of an assembly of tetrahedra. Table VI displays calculated and measured entropies for liquid water; the agreement is very good. E . The Free Energy. At first sight it appears that to calculate the RNM free energy at some temperature all that need be done is to combine the already calculated energy and entropy at the same temperature. However, it is commonly found that approximate theories of the liquid state are thermodynamically inconsistent, e.g., the pressures calculated by two different but thermodynamically equivalent methods are often not the same. In the RNM the heat capacities C,@) = T ( d S / d r ) and Cp(H) (dH/dT), differ by about 20% on average a d h a v e slightly different temperature dependences although, of course, they should be the same. Therefore, a simple combination of the RNM energy and entropy to generate the free energy is subject to magnification of the errors inherent in the individual functions. Sceats and Rice&have addressed this problem and shown that the RNM can be made self-con-
+
t ("C) Comparison of the random network model prediction and the observed energy of liquid water. Flgure 10.
smaller contribution from hydrogen bond stretching. Dissected one step further, it is the smaller contribution of nonbonded interactions to the stabilization, and the larger angular distortions of the hydrogen-bonded network, that account for the high energy of the liquid relative to ice Ih. At this point it is necessary to remark that an important characteristic of the RNM, namely, the temperature dependence of the angular disorder, is known only along the liquid-vapor coexistence curve, and not at constant volume. Accordingly, there is some ambiguity concerning the nature of the computed energy; it is not precisely an energy, nor an enthalpy, although probably closer to the latter. Because the thermal expansion of water is small, and the product of pressure and volume is small compared to the energy, this uncertainty does not have serious consequences. Nevertheless, its existence does highlight a deficiency of the random network model, namely, there does not appear to be any completely satisfactory way of specifying its density. We shall return to this point later. For the present we use energy and enthalpy interchangeably, accepting the inherent ambiguity. D . The Entropy. To calculate the entropy of liquid water according to the RNM we again consider separately the vibrational and configurational contributions. Given the approximate density of states described in section VC, and the dependence of the librational, hindered translational, and flexing modes on the width of the distribution of hydrogen bond angles, computation of the vibrational contribution to the entropy is straightforward. The configurational entropy is determined by the number of ways of constructing a random network with a given energy, number of molecules, and volume. In the RNM, as we have formulated it, the principal variables defining the system's properties are the average 00 separation, R, the mean square deviation of the distribution of 000 angles, $, and the mean square deviation of the distribution of 00 separations, 6R2. In the ice Ih lattice from which the RNM is derived, the only contribution to 3 and 6R2 come from the inhomogeneities induced by proton disorder, so both are small. Then, in the absence of proton disorder, the configurational entropy of the hypothetical perfect ice Ih lattice is defined to be zero, and the configurational entropy of the random network defined by the increase in the volume of phase space occupied relative to that occupied by the hypothetical perfect ice Ih configuration.
[(G2
(45) L. Pauling, J. Am. Chem. SOC.,67, 96 (1935). (46) M. G. Sceata and S. A. Rice, J. Chem. Phys., 72, 3260 (1980). (47)R.J. Bell and P. Dean,Phys. Chem. Glasses, 9,125 (1968);10,
164 (1969). (48)M. G. Sceata and S. A. Rice, J. Chem. Phys., 72, 6183 (1980).
1118
Rice and Sceats
The Journal of Physlcal Chemistry, Vol. 85, No. 9, 1981
TABLE VIII: Average Ring Distribution per Time-Averaged Configuration
, I80 - [
ring size
no. of rings/ configuration
4 5 6 7
16 47 66 83
,
-10
io
30
50
70
90
Temperature (‘Cc)
Figure 11. The average 00 stretching frequency (F,) as a function of temperature for the self-consistent calculation (- - .) and in the original zeroth order RNM (-). The temperature dependence of the measured peak of the Raman spectrum is also shown (- - -).
.
TABLE VI1
ice water
a
T, K 233 269
calcd
obsda
0.387 0.270
0.395 0.253
283 323 363
0.158 0.119 0.097
0.178 0.103 0.060
Reference 50.
sistent by requiring that Cis) Cpm)via adjustment of the librational and hindered translational frequencies in the model density of states. Only minor adjustments to the frequencies are required to achieve self-consistency (see Figure ll),so that the agreement with the various properties of water discussed in sections VA-VD is maintained. They find that the Gibbs free energy generated in the self-consistent procedure is a minimum with respect to variation of the width of the hydrogen bond angle distribution. To test the accuracy of the RNM free energy, Sceats and Rice calculate the ratio of vapor pressures of HzO and DzO, pH/pD. This ratio deviates from unity because the potential and kinetic energy operators do not commute. To order h2 is is found that In (PHIpD) depends on the mean squared force acting on a molecule,49and so can be used as a sensitive test of the accuracy of the RNP and of the suggestion that the important coordinates in the potential energy are the average 00 separation and the hydrogen bond angle. The results obtained are shown in Table VII. The agreement between theory and experiment50is seen to be excellent for ice I at 233 and 269 K, and for liquid water at 283 and 323 K; the agreement is slightly worse at 363 K. In absolute terms, the error at 363 K amounts to 0.11 kJ/mol in the difference of the chemical potentials of HzO and D20, which could easily arise from the inadequacies of the approximations introduced as conveniences for ease of calculation.
VI. Concluding Remarks We have presented evidence in support of the assertion that the random network model provides a simple, modestly accurate, description of liquid water. It is now important to draw attention to the deficiencies of that model. Despite the several successes of the RNM, it is an incomplete description of liquid water. The RNM is constructed from a central configuration which has the connectivity (49) J. Bigeleisen, J. Chem. Phys., 34, 1485 (1961). (50) Data taken from G . Jansco and W.A. van Hook, Chem. Rev., 74, 689 (1974).
of a continuously branching loop-free tree with the local topology of a distorted ice I structure. The possibility of ring formation in the structure is ignored, as are possible contributions from regions of configuration space which correlate with the symmetries of the higher ices. Are the regions of configuration space ignored by the RNM important for the description of water? There is good reason to believe that cooperative motions neglected in the RNM, such as those which generate fluctuations with high local density and grossly distorted hydrogen bonding (possibly multiply distributed about several mean bonding angles), do occur in liquid water. Sceats and Rice48have shown that, at constant temperature, the Gibbs free energy of the RNM has only a weak dependence on and the known thermodynamic and structural properties of the higher ices show that very different molecular arrangementa, with very different densities, can have nearly the same energy. Some information on cooperative fluctuations can also be obtained from computer simulation data. We take the formation of rings of more or less than six water molecules to be one (incomplete) measure of the extent of the cooperative fluctuations in the liquid, since these represent departures from the ice I like central configuration of the RNM. Belch, Rice, and S ~ e a t have s ~ ~analyzed time-averaged molecular dynamics simulation data to determine the distribution of sizes of simple rings of hydrogen-bonded water molecules. Their results are shown in Table VIII. It is worthwhile pointing out that in ice I1 each of the two inequivalent oxygens is involved in 7 six-membered rings and 42 eight-membered rings, while in ice I11 the two inequivalent oxygens are each involved in five-, six-, seven-, eight-, and nine-membered rings. The existence of these bonding topologies suggests that the system potential energy surface has many minima, and that transitions of the system from one to another minimum involve cooperative motion of a large (-12) number of water molecules. The insufficiency of the RNM with respect to description of the ring statistics of water, or possible cooperative fluctuations, should not be overemphasized. As the temperature increases becomes large enough that gross distortions of the central configuration overlap some of the regions of configuration space corresponding to distorted higher ice structures, so those fluctuations not involving cooperative motions of molecules should be adequately sampled. Accordingly, we propose that the RNM be viewed as a description of those properties of liquid water which are independent of the existence of large fluctuations induced by cooperative motions of the molecules. The RNM should, therefore, accurately describe water in the region where the so called “anomalous” contributions to its properties can be neglected. As an example of this expectation we show in Figure 12 the RNM prediction for the heat capacity of liquid water, and the experimental values. It is clear that the RNM prediction is much better outside the temperature region where cooperative fluctuations lead to a A-like anomaly; the latter occurs in supercooled water at 228 K. It can be shown that the difference between the observed and RNM heat capacities can be modeled by a cooperative contribution of the form A(T - Tx)-’f, although the data are not of sufficient pre-
e,
Feature Article
The Journal of Physical Chemistty, Vol. 85, No. 9, 198 1 I
I
I
I
ik
iooc
20t Ii ,
160
240
I
I
I
4
G~(RNM)+
,
, 320
I
1 400
Temperature ( K ) Flgure 12. A comparison of the RNM heat capacity and experimental data. Curve A is obtained from the data of C. A. Angeli, J. Shuppert, and J. C.Tucker, J. Phys. Chem., 77,3092 (1973). The constructlon of the heat capacity curve including the effect of cooperative fluctuations by addffion of the term 1000/( T - 228)"' to CJRNM) is intended to be suggestive; no significance should be attached to the values of the coefficient or the exponent.
cision to permit an unambiguous determination of the critical exponent y'; one possible set of values, namely, A = 1000, y' = 5/4, leads to the contribution sketched in Figure 12. The fact that this functional form is found to fit the difference in heat capacities supports the notions that the contribution to the thermodynamic properties arising from cooperative fluctuations are (a) omitted from the RNM and (b) can be grafted onto the RNM by superposition. It is, furthermore, strongly suggested that, in addition to being the principle internal coordinate of the RNM, has some of the characteristics of a collective coordinate. It is worthwhile reemphasizing that the RNM provides a framework for the analysis of the properties of liquid water which permits unusual insights into the character of the molecular dynamics. As a last example of this, consider the heat capacity of liquid water, which can be related to the fluctuations in the enthalpy, C, 3)/ R P . Given the already used subdivision of the energy into configurational and vibrational contributions, the mean square fluctuation in energy is the sum of the mean
(z
1119
square fluctuations of the configurational and vibrational energies and of their correlated fluctuations. By virtue of the dependence of the vibrational frequencies on the distortion of the hydrogen-bonded network, configurational and vibrational contributions to the energy are not independent. An exact calculation of these contributions to the mean square fluctuation in energy is not at present possible, but Sceats and Rice43have shown how simple arguments can be used to obtain qualitative information about them. They find that the fluctuations in the configurational energy are dominated by the fluctuations in the bending contribution to the energy, and that the correlated cross fluctuations of the vibrational (essentially librational) and configurational energies give rise to a large negative contribution to C,. The existence of this negative contribution is particularly interesting in that it implies that the heat capacity of liquid water is kept low because a fluctuation in configurational energy is balanced by an opposite change in the vibrational energy resulting from a large change in frequency, hence also of zero-point energy. The fact that this cross correlation is not zero also implies that quantum mechanical effects cannot be neglected in the liquid. Further improvement of the RNM must focus attention on the following obvious questions: (i) How is related to the collective coordinates that can be used to describe the fluctuations that determine the anomalous properties of liquid water? (ii) How can one calculate the density of a random network system? (iii) How are the collective coordinates mentioned in (i) related to the density of the random network? (iv) How are the properties of the random network influenced by changes in pressure? Of course, at a more fundamental level, it is necessary to establish how the RNM can be obtained by systematic approximation from the partition function of liquid water, and how that approximation can be improved.
Acknowledgment. This research has been supported by the USPHS and the NSF. We have also benefitted from the use of facilities provided by the NSF for materials research at The University of Chicago.
