J. Phys. Chem. 1983, 87, 4190-4197
4190
estimate f$/f: = 0.55 and use the Ross and Andersson value of 0.4122 W m-l K-’ for the clathrate conductivity at T = 100 K, the conductivity at 250 K is estimated to be 0.507 W m-I K-l as compared with the experimental value of 0.504 W m-I K-l. A definitely higher value, e.g., 0.6, for fi/f: gives 0.55 W m-l K-l. Thus it is meaningless to look for better or worse agreement without a more elaborate model. But it is evident that the simple model contains all the essential physics. It is important to note that the T dependence of X(T) is controlled by at least three competing factors, viz., (i) an almost linearly increasing T dependence in the heat capacity, (ii) an almost linear decreasing T dependence in the density, (iii) a complicated T dependence in the factor f ( T ) which corrects for rotational and other “nonconductive” contributions to the heat capacity, which become important for high T. The approximate form exp(-To/ T)observed by Ross and Andersson is merely an average fit to these several concurrent effects and has, in our opinion, no direct physical significance. V. Concluding Discussion We have shown that the thermal conductivity of icelike systems can be predicted, in most cases, to good accuracy via the mean-free path formula if unit cell number, density, and the mean bond length were known. It is also clear that the partitioning of the heat capacity between rotational and vibrational modes has to be taken into account if temperature-dependent effects are to be understood. Since the thermal conductivity of ice Ih is fairly well understood from a microscopic sense,2 and since the other ices scale according to eq 7, it may be said that the microscopic
considerations of ref 2 are applicable to the other ices, and icls. The agreement also c o n f i i s the theory of Roufosse and Klemens. The extreme case of the clathrates has shown that in sufficiently large crystals the effective phonon mean-free path drops to a limiting value (typical of the lattice constant) which is independent of n. This explains the previously intriguing fact that icl I and ecl I1 have nearly the same thermal conductivities. It does not call for an assumed special “glassy“behavior in the clathrates, or a special guest-host interaction effect. It is clear that thermal data for the icelike NHIF phasesz0 could be examined within the formalism used here. Further, it may be possible to form clathrate-like structures with NH4F. Their thermal conductivities could provide added insight into the behavior of h ( T ) as a function of n, p , etc. Acknowledgment. The author acknowledges many discussions with John Cook and Marek Laubitz. Discussions with D. W. Davidson, E. Whalley, John Tse, and Derek Leiast, of the Chemistry Division, and the assistance of their research groups are gratefully acknowledged. Finally, the author thanks the referee for many valuable comments. Registry No. Water, 7732-18-5. (20)R. G.Ross and 0. Sandberg, J. Phys. C, 11,667(1978). (21)B. Kamb in ‘The Physics and Chemistry of Ice”, E. Whalley, S. J. Jones, and L. W. Gold Ed., The Royal Society of Canada, Ottawa, 1973, p 28. (22)R. K. McMullan and G. A. Jeffrey, J. Chem. Phys., 42, 2725 (1965).
A Molecular Dynamics Study of Water Clathrates P. L. M. Plummer” and T. S. Chen Department of Physics and the Oreduate Center for Cloud Physlcs Research, Unhwslty of Mlssouri, Rolla, Mlssouri 6540 1 (Received: September 27, 1982; In Final Form: December 13, 1982)
Molecular dynamics calculations were carried out for a series of temperatures between 67 and 315 K for a 20 water molecule microcluster. These studies provided geometricalproperties such as pair distribution functions, dynamic properties such as power spectra, and thermodynamic properties such as the caloric equation of state. The molecules were originally arranged in a pentagonal dodecahedral clathrate confiiation. The cluster showed a change in structure and other properties around 200 K but did not dissociate even above 300 K.
I. Introduction In the attempt to gain insight into the growth and decay of small aggregates of water, molecular dynamics simulation of cluster systems has been undertaken. Understanding of the molecular aspects of growth, decay, and stability (or lifetime) of these aggregates is necessary for modeling the nucleation process. Of particular interest to this study are the stability of a given hydrogen bonded configuration as a function of time and temperature, the evolution of an isolated cluster in terms of its radial distribution function and power spectrum, barriers to structural interconversion, and energy distribution among the molecules. That water is exceedingly complex has become increasingly apparent. Its permanent dipole moment, its 0022-3654/83/2087-4190$01.50/0
internal degrees of freedom, the fact that it is an equally good donor and acceptor of protons in hydrogen bonds, and the flexibility of these hydrogen bonds, all contribute to the difficulty of designing a satisfactory potential to describe its interactions. However, during the past decade, several potentials’ have been proposed which yield reasonable results for many of water’s properties.2 These (1)(a) F.H.Stillinger and A. Rahman, J. Chem. Phys., 68,666(1978); (b) 0.Matsuoka, E. Clemeneti, and M. Yoshimine, ibid., 64,1351(1976); (c) J. C. Owicki, L. L. Shipman, and H. A. Scheraga, J.Phys. Chem., 79, 1794 (1975); (d) R. 0.Watts, Chem. Phys., 26, 367 (1977); (e) F. H. Stillinger and C. W. David, J. Chem. Phys., 69,1473 (1978);(0 W.L. Jorgensen, J.Am. Chem. Soc., 101,2011(1979);( 9 ) P. Barnes, “Progress in Liquid Physics”, C. A. Croxton, Ed., Wiley, New York, 1978,Chapter 4; P. Barnes, J. L. Finney, J. D. Nicholas, and J. E. Quinn, Nature (London), 282,459 (1979).
0 1983 American Chemical Society
A Molecular Dynamics Study of Water Clathrates
The Journal of Physical Chemistty, Vol. 87, No. 21, 1983 4191
TABLE I: Comparison of Calculated (CFM)a and Experimental Water Properties
CFM
expt
~
1A
Figure 1. Initial configuration of the 20-molecule pentagonal dodecahedron,
potentials, together with Monte Carlo and molecular dynamics simulation^,^ have greatly refined our knowledge of water and water systems. During this same period of time, there has also been an increasing interest in the properties of water clathrate ices. Clathrate structures have served as models for nucleation phenomena4 and for ion cluster stability5 in molecular beam systems. The clathrate ices appear to act as traps for large amounts of natural gas in the Arctic regions and beneath the sea floor.6 Unfortunately, existing physical data on clathrate properties, either experimental or theoretical, are very meager. The current study was initially undertaken to examine the assumed stability of the water clathrate in the gas phase. The system chosen to model the water clathrates is the 20-molecule pentagonal dodecahedron shown in Figure l. This microcluster was used previously in the development of a molecular model for prenucleation water cluster^.^ Castleman and others5p7 have reported unusual stability for water clusters in molecular beam experiments at sizes corresponding to single clathrate cages. Since our studies examine the dynamics and average properties of an isolated water microcluster, the results may not be directly applicable to bulk clathrate hydrates, but should provide some insights into bulk properties, just as "spherical" face-centered cubic microclusters provide information about both surface and energy properties of rare gas solids.8 In calling our system a microcluster, we use the terminology due to Hoareg who defines a microcluster as an aggregate of molecules so small that at any given time an appreciable percent of its molecules are in the surface. The term, microcluster, also implies a lifetime for the aggregate which is considerably longer than the characteristic inverse (2) A recent paper, M. D. Morse and S. A. Rice, J.Chem. Phys., 76, 650 (1982), discusses many of these potentials as applied to the polymorphs of ice. (3) See ref la, c, d, and g as examples. (4) P. L. M. Plummer and B. N. Hale, J. Chem. Phys., 56,4329 (1972). (5) P. M. Holland, and A. W. Castleman, J. Chem. Phys., 72, 5984 (1980); S. S. Lin, Rev. Sci. Instrum., 44,516 (1973); J. Q. Searcy and J. B. Finn, J. Chem. Phys., 61, 5282 (1974). (6) D. W. Davidson, 'Water, A Comprehensive Treatise'", F. Franks, Ed., Plenum, New York, 1973, Chapter 2. (7) J. L. Kasaner, Jr., and D. E. Hagen, J . Chem. Phys., 64, 1861 11976). . (8)(a) W. D. Kristensen, E. J. Jensen, and R. M. Cotterill, J. Chem. Phys., 60, 416 (1974); (b) D. J. McGinty, ibid., 58, 4733 (1973); (c) C. L. Briant and J. J. Burton, ibid., 63, 2045 (1975). (9) M. R. Hoare, Adu. Chem. Phys., 40, 49 (1979).
equilibrium properties (water monomer) b o n d length, A b o n d angle, deg dipole m o m e n t , D dissociation energy, kcal/mol intramolecular frequencies, c m - ' asymmetric stretch symmetric stretch bend cluster binding energy, kcal/mol dimer clathrate
0.9584 104.45 1.86 134.522
0.9584 104.45 1.86 219.34
3805.70 4266.01 1369.33
3755.79 3656.65 1594.59
-5.838 -176.660
-5.446
'
Central force model; Stillinger and Rahman.'" L. A. Curtiss, D. J. Frurip, and M. Blander, J. Chem. Phys., 71,
2703 (1979).
frequencies for its internal motions. Thus, such an aggregate can be identified with a bounded region of phase space within which time averages over its natural motions are meaningful. In section I1 we describe the potential model used and the details of the computational procedures. Section I11 gives the results of the simulations. The last sections discuss and comment on the observations and summarize our findings. 11. Method
In molecular dynamics the classical equations of motion for the molecules are numerically solved to produce the position and velocity of the particles as functions of discrete time steps. The central force pair potential model of Stillinger and Rahman'O is used to describe the interaction between each pair of atoms. At normal temperatures and pressures, the water monomer produced by this potential is stable and has properties'l in reasonable agreement with experiment (Table I). A central difference method due to Verlet12 is used to integrate the resulting equations of motion. Since we are simulating a microcanonical system, the total energy of the system serves as a valuable aid to assess the choice of time step size and algorithm stability. We found that a time step of 0.00025 ps provided energy stability of 0.001 kcal/mol per molecule. This time step together with the Verlet algorithm provided time reversability of approximately 50000At.13 We emphasize that these simulations are for an isolated cluster. There are no periodic boundary conditions, and the cluster is not confined. Molecular evaporation and cluster disintegration are possible and not restricted. The microcluster consisting of twenty water molecules was felt to be of sufficient size to have properties approaching those of a microscopic1* surface and was approximately the largest cluster affordable for long simulations. Our choice of initial configuration was determined (10)F. H. Stillinger and A. Rahman, J. Chem. Phys., 68,666 (1978). (11) In addition, our studies on water pentamers show the CF model predicts the same trends in stability of open vs. closed structures as the polarization model of Barnes**and ab initio quantum mechanical calculations. (12) L. Verlet, Phys. Reu., 159, 98 (1967). (13) T. S. Chen and P. L. M. Plummer in "AtmosphericAerosols and Nuclei", A. F. Roddy and T. C. O'Conner, Ed., University of Galway Press, 1981, pp 75-80; T. S. Chen, Ph.D. Dissertation, University of Missouri-Rolla, Dec 1980. (14) 0. Sinanoglu, J. Chem. Phys., 75, 463 (1981).
4192
The Journal of Physical Chemistry, Vol. 87, No. 21, 1983
e
Plummer and Chen
I
INPUT C O N F i G U R h T I f 3 4
I
O P T I M I Z E STRUCTURE
--
I I
-138.
0
E
-7-
\
7;
m -139.2
I
E W I LL I BRAT I ON, DATA FOR AVERAGE PROPERTIES
-
-4
x w I
140.1
i
,/’J
1141.
I DISSOCIATION OB FRAGMENTATIm
I
, ,+ I
50
150
250
350
T (K)
Flgure 2. Flow chart for the heating-equilibration procedure.
by the interest in clathrate stability mentioned earlier. A previous molecular dynamics study13showed the hydrogen arrangement depicted in Figure 1 to be somewhat more stable against disruptions due to thermal fluctuations than a more random placement of hydrogens. Initially, the nearest-neighbor 0-0distance is chosen as 2.76 A, characteristic of ice. The intramolecular 0-H distance is 0.9584 A. The H-0-H angle is 104.45 if only one hydrogen of the molecule is involved in hydrogen bonding and 108.0 if both are donated. For simplicity the initial velocities of the atoms are set to zero. The cluster then evolves for 1000At. After this initial equilibration, physical properties of the cluster are obtained by time averaging. If the time is sufficiently long, according to the ergodic hypothesis,15 this is equivalent to the ensemble average. We used 32768At, or approximately 8 ps, for computing average properties at each temperature examined. The temperaturesb is calculated from the time average of the kinetic energy
In eq 1 N is the total number of atoms; mi and ui are the mass and speed of the ith particle, respectively. Six degrees of freedom are subtracted since the cluster as a whole has neither translational nor rotational motion. The brackets denote the time average for the entire simulation. The temperature of the cluster is increased by adjusting the velocity of the atoms. A new temperature, or cluster energy, is achieved by alternately heating and equilibrating the cluster. The heating rate employed depended on the initial temperature, but ranged from 10 to 40 K/ps.16 After the desired energy is reached, the system is equilibrated and then the 8-ps simulation run made. A schematic for this procedure is shown as Figure 2. The power spectrum17 is also calculated at each of the temperatures from the Fourier transformation of the velocities:
(15) K. Huang, ‘Statistical Mechanics”,Wiley, New York, 1963, p 203. (16) If the heating is too rapid or there is insufficient equilibration prior to adding more energy, excessive distortions of the structure occur.
(17) P. 0. Esbjorn, E. J. Jensen, W. D.Kristensen, J. W. Martin, and L. B. Pedersen, J. Comput. Phys., 12, 289 (1973).
Figure 3. Energy (per molecule in kcal/mol) vs. temperature for the 20-molecule cluster. Bars indicate temperature flucuations.
Here Cij(uk)is the discrete velocity spectrum of the j t h component of the ith atom. The power spectrum provides information about the vibrational motion and we can estimate the diffusion constantl8 of the water molecules as (3)
where T is the total simulation time. Three pair distribution functions, g(OO), g(OH), and g(HH), are calculated for each of the temperatures. In addition to the radial distribution functions and the power spectrum, averages for the vibrational and rotational energies and, for each temperature, the total internal potential energy and kinetic energy for the cluster are determined. In section 111, we present the quantitative results of these temperature studies, and in section IV we discuss the more qualitative information obtained from direct viewing of the temporal sequences of “snapshots” of the system as aids to interpreting these results. 111. Results
The heating-equilibration process was carried out at eleven temperatures between 67 and 315 K. Equilibrationlg at each temperature was substantially complete within 500 time steps after the last heating cycle. The rapid equilibration indicates the ease with which energy flows in the cluster. This rapid redistribution of energy has also been found in a preliminary study13of the collision of a monomer with the 20-molecule cluster. In that study, the incoming molecule quickly transferred its impact energy to the cluster normal modes and was absorbed. The plot of energy vs. temperature is shown in Figure 3. The slope changes in the vicinity of 200 K, having a value of 11 cal/mol per molecule per Kelvin below 200 K and a value of 26 cal/mol per molecule per Kelvin above. The mean interaction energy per molecule in the cluster changed from -8.83 kcal/mol at 5 K to -6.78 kcal/mol at 315 K. For bulk liquid water, the CF potential model gives (18) C.L.Briant and J. J. Burton, J. Chem. Phys., 63, 3327 (1975). (19) Evidence for completeness of equilibration is obtained by comparing average properties calculated for the first 4 ps with the last 4 pa, and with those determined from the entire simulation period.
The Journal of Physical Chemistry, Vol. 87, No. 27, 1983 4193
A Molecular Dynamlcs Study of Water Clathrates
0.U6
-0.3
3 iw 0.23 v
0.11
0.
407.
0.
611.
814.
1018
(CM-’) 0.39
.
:
:
:
:
:
:
:
,
:
:
:
:
:
;
:
;
:
1.01,:
:
b 1
0.80
0.31
:
h
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
:
T-3 15K
I
:
I t
0.119
0.39
-0.29
3
v
%20
0.
(107.
611.
814.
1018
(CM-’)
Figure 4. Power spectrum of the 20-molecule cluster for the 0-900-cm-’ range at (a) 67 K, (b) 159 K, (c) 194 K, (d) 276 K, and (e) 315 K.
-9.48 kcal/mol at 29.5 OC.l0 Following Briant and Burton18 we calculated the surface energy of the cluster as ENs(r) = ENc(r)
- NEB + 6kT
(4)
where N is the number of water molecules, ENc is the energy of the cluster, EBis the energy per molecule in the bulk, and 6kT adjusts the kinetic energy of the cluster for the missing translational and rotational degrees of freedom. A t 298 K we get 58 kcal/mol for the surface energy of the cluster. This is remarkably close to the 62 kcal/mol calculated by Briant and Burton1* who used a differential potential model. The power spectrum for the cluster at each of several temperatures is plotted in Figure 4. The spectra show three groups of peaks in the 0-900-cm-1 range. Isotopic substitution of the hydrogen atoms by deuterium in the
calculation of the spectrum shows that the lowest peaks can be attributed to the motions of the water molecule as a whole or to cooperative motions of several water molecules. The peaks at the higher frequencies of this range are primarily due to librational motions of the molecules about their ”lattice sites”. A more detailed assignment of the lowest two groups may be obtained by comparison with a previous analysism of clathrate normal modes using point molecules and a two-force constant, harmonic potential. Comparison with these studies permits us to assign the group of peaks around 51 cm-’ to normal modes of the pentagonal rings, the peak at about 160 cm-’ to the symmetric stretch of the (20) P. L. M. Plummer, ‘Physics and Chemistry of Ice”, E. Whalley, S. J. Jones, and L. W. Gold, Ed., Royal Society, Canada, 1973, p 103.
4194
Plummer and Chen
The Journal of Physical Chemistry, Vol. 87, No. 21, 1983
i
:
:
:
:
:.
:
:
:
:
:
:
:
:
:
:
:
:
a
1
12.0
32.0
9.00
2u.o
0 0 06.00
3.00
0.
0.0 6.00
3.00
0.
12.0
R ( A 1 9*00
3.00
0.0
6.00
12.0
R ( A )9*O0 T49UK
12.0
1
32.0
9.00
0 6.00 0)
3.00
t
0.0 0.0
1 t
12.0
3.00
6.00
R(A :
:
:
:
:
:
:
:
:
Fo0 :
:
12.0
, , ::
:
:
T-276K
:
.
'I
0.0
3.00
6.00
12.0
R ( A ,9*O0
,
32.0
1j
T-276K
/
1
0 09'00
1
0.0
3.00
6.00
R ( A 1 '*0°
12.0
0.0
3.00
9.00
6.00
12.0
R(A)
Flgue 5. Oxygen-oxygen pair dlstribution function for the 20-mdecule cluster at (a) 67 K, (b) 194 K, and (c) 276 K.
Flgure 6. Oxygen-hydrogen pair distribution function for the 20molecule cluster at (a) 67 K, (b) 194 K, and (c) 276 K.
entire clathrate (pure hydrogen bond stretching), and the group at 200 cm-l to more complex torsional-type group motions. Diffusional-like motion, rather than oscillatory motion, is revealed by modes at zero frequency. With eq 3, diffusional constants are obtained which range from 1.9 X cm2/s at 67 K (no diffusional motion) to 2.5 X cm2/s at 250 K. The molecular dynamics result for bulk liquid at 29.5 "C is 1.12 X cm2/s.l0 Ohr result is similar to that for the simulation of argon in that a value characteristic of the bulk liquid is achieved at a much lower temperature in the molecular dynamics. For each temperature the 0-0, the 0-H, and the H-H distribution functions were calculated. In Figures 5 and 6 we plot the 0-0and 0-H distributions for several temperatures. These two distribution functions best illustrate the change in structure and hydrogen bonding patterns
with temperature. The first peak of g(0H) is due to the intramolecular O H pair. The second peak of g(0H) and the first peak of g(O0) are attributed to hydrogen bonds. Integration over these peaks provides estimates for the average number of bonds present. Both distributions show 30 bonds at 67 K, pass through a maximum at 194 K, and remain at high values of about 26 bonds at temperatures above 300 K. In Table I1 we list the number of hydrogen bonds for the cluster as determined by either a bond energy or a bond distance criterion. Because of the rotational motion of the water molecule, the bond energy can change continuously from negative to positive values at a fixed bond distance. Thus, a distance criterion alone is insufficient to define the formation of a hydrogen bond. On the other hand, there is considerable uncertainty in choosing a threshold energy for the hydrogen bond. The table il-
The Journal of Physical Chemistry, Vol. 87, No. 21, 1983 4195
A Molecular Dynamics Study of Water Clathrates
TABLE 11: Number of Hydrogen Bonds” in the Cluster as a Function of Temperature for a Range of Energy Criteria temp,K
67
105 1 5 6 1 9 4 226 276 293 315
V H B~ 1.0 VHB - 2.0 VHB- 3.0 VHB- 4.0
30 30 30 30
30 30 30 30
39 31 30 23
37 30 26 23
35 28 26 23
33 27 24 19
29 24 22 17
26 23 19 12
g( OH)
30
30
31
32
30
28
27
26
-~
” T h e n u m b e r of b o n d s as determined b y an energy criterion is representative of a n instantaneous cluster structure, whereas t h e n u m b e r of b o n d s predicted f r o m t h e radial distribution f u n c t i o n represents a t i m e average. kcal/mol.
lustrates the range in the number of hydrogen bonds for the cluster at several different temperatures as determined from the distribution function, and different distance and energy criteria.21 The average number of hydrogen bonds provides only a partial picture of the cluster’s structure. The lifetime of these bonds is equally important. For example, at 67 K no bond is broken during the 8-ps simulation, while at 276 K no bond had a lifetime which exceeded 6 ps. This point will be discussed more in the following sections. To begin to assess the role of a “guest” molecule in the clathrate structure, a series of parallel calculations were made in which a twentyfirst water molecule is introduced at the center of the cluster at 67 K. The system is then equilibrated and time averages are computed. Heating, equilibration, and simulation runs on the twentyonemolecule system were carried out to provide data at several additional temperatures. In contrast to other studies6 which concluded the guest molecule has only a minor role in clathrate hydrate properties, these preliminary studies show the water guest molecule disrupting the bonding pattern of the single-cage microcluster and producing structures which differ from the unoccupied clathrate at the same temperature. Additional simulations for this system, including a variety of guest molecules, will be the subject of a future report.
IV. Discussion The data presented in the previous section, together with qualitative observations (of successive plots) permit a rather detailed picture to be formed of the microcluster dynamics as a function of temperature. By using the molecular dynamics method, we have gained information about energy flow processes in a microcluster as the structure evolved, whether cooperative molecular motions can either preserve or disrupt the cluster, and how internal structure of the individual units of the microcluster affect the observed properties. Both the rotational asymmetry and the ability of the individual units (water molecules) to distort, in response to a changing environment, could decrease the sharpness of any change of phase. We also hoped that the study of the isolated microcluster, its geometry, its power spectrum, and its reaction to energy injection would aid in the understanding and analysis of bulk ice, clathrate hydrates, and amorphous ice, and provide information to help identify the origin of many of the features of the bulk water systems. It is of particular interest to understand the causes and significance of the changes in slope of the E vs. T curve. Cotterill et al.& have found from two-dimensional studieszz (21) A pair of molecules is counted as hydrogen bonded if the interaction energy exceeds VHBand R ( O 0 ) is less than 3.2 .&. (22) R. M. J. Cotterill, W. D. Kristensen, J. M. Martin, L. B. Pedersen, and E. J. Jensen, Comput. Phys. Commun., 5, 28 (1973).
that many indicators of change of phase are no longer reliable in the presence of a free surface. They found that E vs. T curve, did remain the caloric equation of useful for describing some aspects of the phase transition. However, the meaning of phase transition in a microcluster is not entirely clear. Thus, we are seeking additional insights from examination of bonding patterns, strength of bonds, pair distribution functions, or other time-averaged properties of the system to clarify the change in the caloric equation of state. A t the lowest temperature, 67 K, the molecules have insufficient energy to change the bonding pattern and remain “locked” into the lattice sites of the pentagonal dodecahedron. Thus, the hydrogen-bonding pattern did not change during the 8 ps. The H-0-H angle for the central water molecules, the double hydrogen donors, had decreased to l o l o , and the 0-H bond lengths differed between donor and nonbonded hydrogens by about 0.15 A. At this temperature, the number of hydrogen bonds determined by either the radial distribution functions, energy criteria for a wide range of energies, or distance criteria, are the same. (See Table 11.) Below 100 K, the cluster is executing only normal-mode motions. The existence of the three peaks centered around 51 cm-’ in the power spectrum (Figure 4a) has been associated with the normal-mode motion of the pentagonal rings. These same three peaks are also present in the spectrum for the bulk clathrate hydrate23when calculated by using the Clementilb potential. At the next higher temperature, 105 K, the hydrogen bond configuration remains unchanged while the amplitudes of the lowest frequency modes are increased more than the higher. When the temperature is increased to 159 K, some hydrogen bond disruption is observed, and the peaks in the g(O0) and g(0H) associated with hydrogen bonding begin to broaden. The mechanism for the bond disruption is an increase in energy in the librational motion. Bond disruption by stretching only requires more kinetic energy. Additional distortion of the hydrogen bonds is shown by the analysis of Table 11. The number of hydrogen bonds estimated by the distribution function increases to 31.4 for g(O0) and 31.2 for g(0H). Further examination of the detailed motion associated with bond disruption or distortion due to librational motion shows an intermediate transition state with a “bent” hydrogen bond. If molecule A is the donor molecule in a linear hydrogen bond with molecule B, then a simultaneous rotation of both A and B can result in either a bifurcated or cyclic hydrogen bond. Further rotation2* of A and/or B can change B into the donor molecule (Figure 7). As the temperature is increased to 194 K, the second and third peaks of g(O0) begin to merge. The separation between the two lowest sets of frequencies in the power spectrum became less pronounced. The center of the librational band shifted down approximately 50 cm-l, and a smaller shift occurred in the lowest two bands. The number of hydrogen bonds, as determined from the distribution functions, is at a maximum of 32.9 for g(O0) and 31.7 for g(0H) at this temperature. There is also a change in slope in the E vs. T curve (Figure 3). For the next three temperatures studied, 226,249, and 262 K, the number of hydrogen bonds decreases mono(23) J. s. Tse, M. L. Klein, and I. MacDonald, Symposium on the Physics and Chemistry of Ice, Rolla, Missouri, 1982. Abstracts, p 120, and private communication with Tse. (24) These cooperative motions can result in a change in bonding pattern at a lower energy than that required by breaking and re-forming bonds.
4196
The Journal of Physical Chemlstry, Vol. 87, No. 21, 1983
1
I
Flgure 7. The top sequence iluatrates a cyclic bond with molecule A changing from a donor to an acceptor molecule. The bottom sequence illustrates a “proton switching” reaction which proceeds through a bifurcated configuration.
tonically. Disruption of the hydrogen bond occurs more frequently. This is indicated in the increased width of the hydrogen bond peak in the g(0H) distribution, and its shift to longer distances. More peaks appear in the power spectra in the O-lOO-cm-’ band, and the amplitude of this band increased to three times the other bands. The cluster expanded somewhat with the molecules, on the average farther from the center of mass. Often two molecules simply move away from each other with a third molecule moving into the vacant spot, giving rise to a diffusional type of motion. As the cluster temperature is raised to 276 K, the maximum of the libration band of the power spectrum shifts from the higher to lower end of the band. This shift is attributed to quasi-rotational motion of some of the molecules which occasionally have only one or two bonds. At this temperature, we also see an occasional molecule acquire a translational kinetic energy greater than its interaction energy with the rest of the cluster and start to escape. However, within a few tens of time steps, reorientation of the molecules increases the interaction energy, and the escaping molecule is reabsorbed. We have seen evidence of similar cooperative behavior in collision studies.13 We did not see a complete evaporation event even at the highest temperature, 315 K. From the examination of “snapshots” of the time steps surrounding a %ear“ evaporation, we see that, in addition to reorientation, the nearest molecules also distort slightly, thereby increasing the interaction. Thus, the internal structure of the monomers appears to increase the stability of the cluster and permits it to persist to higher energies. This is not a large effect, as “freezing” the internal modes decreased the bond lifetime by 20-30% above 200 K. At 315 K, the number of hydrogen bonds determined is 27.4 from g ( O 0 ) and 26.8 from g(0H). The cluster is somewhat more compact than at 250 K, and for hydrogen bonds characterized by an attractive interaction and a 0-H distance of less than 2.5 A,the number of molecules participating in four hydrogen bonds has increased. However, the av-
Plummer and Chen
erage energy per bond is less than 3 kcal/mol. There are also increasing numbers of nonlinear hydrogen bonds, many involving two hydrogens in either cyclic or bifurcated configurations. The maximum lifetime of any specific hydrogen bond at this temperature is about 6 ps. Many are broken and re-formed very quickly with no change in the overall bonding pattern. However, above 276 K, broken bonds often lead to an altered structure. Montrose et al.25have estimated lifetimes of about 1 ps at similar temperatures, a value consistent with our results for a more stringent bond definition of VHB > 3.1 kcal/mol.
V. Comments and Conclusions Examination of the radial distribution functions, the power spectra, and the energy vs. temperature behavior indicates a change in the microcluster properties around 200 K. Below this temperature, both the instantaneous and the average structure have the clathrate bonding pattern. Above 230 K, the cluster bonding pattern varies widely, with instantaneous structures which can differ substantially from each other and from the average structure. The change in slope in the E w. T curve at -200 K can be considered as a change in the microcluster from the open, semirigid clathrate structure to a more compact, less-ordered system. That this change is not sharp is consistent with the finite size of the system and a change from an ordered to a disordered configuration. The increase in heat capacity is also consistent with more disrupted or distorted hydrogen bonds and a decrease in long-range order. The changing bond pattern is also reflected in the loss of definition in the second and third peaks in the 0-0distribution which has an increasingly liquidlike appearance at 220 K and above. A structural change in this temperature range may be a general characteristic of the clathrate bonding, as Tse et al.23found their bulk clathrate structure “melted“ or became unstable at -220 K. If the changes in cluster properties above 200 K are attributed to a type of transition point, the fact that this transition occurs at a much lower temperature than the bulk reflects the more open bonding pattern of the clathrate and the greater free surface for the microcluster. Both provide for a greater amplitude of motion at a given temperature for the molecules in the microcluster. The frequencies in the power spectrum shift in a similar way with temperature as in microclusters of Lenard-Jones atoms. The Lenard-Jones clusters also become liquidlike at temperatures substantially below bulk values. However, we did not observe a sharp or abrupt change of phase, as was the case for the clusters of atoms.8a Neither did we see an abrupt increase in the diffusional constant as was observed for argon ~1usters.S~ The small size of the clusters, together with an initial configuration which consisted only of surface atoms, effectively precludes a “two-phase” system such as observed by Weber and StillingerZ6for a 250-molecule cluster. Our observations of the microcluster gave no evidence in support of either the iceberg or mixture model for liquid water.n We found the hydrogen bonds to be quite flexible with a strong tendency to maintain the maximum number of near neighbors, as opposed to maintaining a fewer number of optimal hydrogen bonds. This conclusion is (25) C. J. Montroae, J. A. Bucaro, J. Marshall-Coakley,and T. Litovitz, J . Chem. Phys., 60, 5025 (1974). (26) T. Weber hnd F. H. Stillinger, J . Phys. Chem., this issue. (27) D. Eisenberg and W. Kausmann, ‘The Structure and Properties of Water”, Oxford University Press, London, Chapter 4, 1974.
A Molecular Dynamics Study of Water Clathrates
further supported by the nearest-neighbor peak in the 0-0 distribution function. The first peak remains sharp and clearly separated from the rest, indicating that the immediate environment of the water molecule is not drastically altered. An increase in energy and, hence, in libration motion decreases the hydrogen bond strength or number at any single instant in time, but despite the “flickering”on and off of the hydrogen bonds, the number of nearest neighbors remains essentially constant. Examination of the bonding patterns found the cluster to consist of networks of cross-linked chains of hydrogen bonds. We seldom observe distinct polyhedra forming within the microcluster. Monomers and dimers with an occasional larger polymer would appear only to be reabsorbed within a few hundred time steps or less. Also, at the higher temperatures distorted bonds or those which involve two hydrogen bonds between a pair of molecules were found with increasing frequency. Since, within the time of the simulation at temperatures above 250 K, every hydrogen bond in the cluster has broken even if it subsequently re-forms, the average properties are being computed for more than one structure or bonding pattern. Shorter time averages provide information about instantaneous or I structure while the &ps averages provide information about an average structure. At the higher temperatures, large regions of phase space are sampled in relatively short times. This leads to large temperature fluctuations as the number of hydrogen bonds change. Qualitatively, the structure changes from one of low potential energy to one of high potential energy with a corresponding change in the kinetic energy and, hence, in the cluster temperature. Analysis of the X-ray spectra of both crystalline and amorphous iceB suggests that the latter contains molecules which are “properly” bound in the sense of having optimal hydrogen bonds, as well as quasi-free molecules which have either broken or highly distorted hydrogen bonds. Thus, the time-average structure as revealed by the molecular dynamics simulations of microclusters may aid in identifying the origin of these features in the bulk systems. There is increasing evidence29that many of the observations on bulk systems reflect the time average of the local environment. An example is the change in position and shape of the peaks in the power spectra of the microcluster with temperature which closely mimic those observed in Raman spectra29for bulk water.30 Also, the activation energy required for molecular rotation and hydrogen bond switching is approximately the same as assumed by Montrose et al.25for liquid water and that predicted by (28) E. Gilberg, M. J. Hanus, and B. Foltz, J. Chem. Phys., 76, 5093 (1982). (29) R. Basil, J. Wiefe-Akenten, and J. L. Taaffe, J. Chem. Phys., 76, 2221 (1982), and private communication with Basil. (30) Giguere and G u i l l ~ have t ~ ~ shown that the addition of acid affects spectra in the same way as increasing temperature and have offered this as evidence for a decrease in the population of strong hydrogen bonds in favor of weaker ones (bent?) between water molecules”. (31) P. A. Giguere and J. G. Guillot, J . Phys. Chem., 86, 3231 (1982).
The Journal of Physical Chemistry, Vol. 87, No. 21, 1983 4197
neutron defraction studies on ice.32
VI. Summary For the 20 water molecule microcluster, we have used molecular dynamics to examine geometric quantities such as pair distribution functions; dynamic quantities such as power spectra, diffusion coefficient, structural changes, and hydrogen bond lifetimes; thermodynamic quantities such as potential and kinetic energies and caloric equation of state; and the pseudodynamic quantities of surface energy and transition point. Consideration of the results of these calculations provide the following observations: (1)The clathrate-type bonding pattern is stable to about 200 K. (2) The energy vs. temperature curve shows a change in slope in the same temperature range as the hydrogenbonding pattern of the clathrate begins to exhibit considerable disruption and the open structure collapses. (3) The power spectrum has three bands in the 0-900cm-’ region. These are centered at -50,200, and 450 cm-’ at 67 K and shift to lower frequencies with a loss of definition and an increase in width as the temperature increases. (4) The average number of hydrogen bonds as determined by the OH distribution function remained high (-25) above 300 K.33 The large number of bonds was maintained at the expense of the energy per bond and the lifetime of a specific bond. (5) The average local environment of the individual water molecules as reflected by the number of nearest neighbors persisted even though the number of bonds per molecule fluctuated widely at the higher temperatures. (6) The ability of the individual units (monomers) to distort appears to increase the cooperative effects of hydrogen bonding and, hence, stability of the cluster to dissociation. (7) The microcluster did not separate into smaller clusters nor were polygons observed which persisted beyond 0.05 ps at the temperatures studied. (8) The structure of the microcluster at the higher temperature is best described as a hydrogen-bonded network with many of the bonds highly distorted at any instant in time. Acknowledgment. The authors thank 0. R. Plummer and S. H. Suck-Salk for helpful comments, and the University of Missouri-Rolla for providing much of the computer time needed for this study. This work was supported in part by NASA Grant NAS8-31150 and the Atmospheric Science Section, National Science Foundation, under NSF ATM77-72614 and NSF ATMSO-19752. Registry No. Water, 7732-18-5. (32) M. R. Chowdhury, J. C. Dore, and J. T. Wenzel, J.Noncrystall. Solids, in press. (33) This cluster temperature corresponds to an added energy of 100 kcal/mol equivalent to 20 hydrogen bonds to the 67 K cluster.
