Published on Web 06/12/2002
On the Debye-Waller Factor of Hexagonal Ice: A Computer Simulation Study Hideki Tanaka† and Udayan Mohanty*,‡ Contribution from the Department of Chemistry, Faculty of Science, Okayama UniVersity, 3-1-1 Tsushima-naka, Okayama 700-8530, Japan, and Eugene F. Merkert Chemistry Center, Boston College, Chestnut Hill, Massachusetts 02467 Received August 9, 2001. Revised Manuscript Received January 2, 2002
Abstract: We investigate by molecular dynamics (MD) simulations the temperature dependence of the Debye-Waller (DW) factor of hexagonal ice with 25 different proton-disordered configurations. Each initial configuration is composed of 288 water molecules with no net dipole moment. The intermolecular interaction of water is described by TIP4P potential. Each production run of the simulation is 15 ns or longer. We observe a change in slope of the DW factor around 200 K, which cannot be explained within the framework of either classical or quantum harmonic approximation. Configurations generated by MD simulations are subjected to the steepest descent energy minimization. Analysis of the local energy minimum structures reveals that water molecules above 200 K jump to other lattice sites via some local energy minimum structures which contain some water molecules sitting on the locations other than the lattice sites. As time evolves, these defect molecules move back and forth to the lattice sites yielding defect-free structures. Those motions are responsible for the unusual increase in the DW factor at high temperatures. In making a transition from an energy-minimum structure to another one, a small number of water molecules are involved in a highly cooperative fashion. The larger DW factor at higher temperature arises from jump-like motions of water molecules among these locally stable configurations which may or may not be a family of the proton-disordered ice forms satisfying the “ice rule”.
1. Introduction
In recent years, there has been renewed experimental and theoretical interest in unraveling the intricacies of describing the equilibrium and dynamical characteristics of supercooled and glassy states of matter.1-3 Some unique features of glassforming liquids include the non-Arrhenius nature of relaxation time with temperature, the apparent relation between the correlation time and the liquid entropy, and the hypothetical Kauzmann temperature (the temperature where the entropy of the supercooled liquid has an intersection with the crystalline entropy and the stretched exponential relaxation of the time dependence of measurable quantities).1-3 Recently, a phenomenon associated with a glasslike transition has been exhibited by the temperature dependence of the meansquared displacement or the Debye-Waller (DW) factor.1 A distinctive change in the slope of the DW factor at low temperatures is observed for a number of glass-forming liquids1,4 and proteins5,6 by using a variety of techniques such as neutron4 * Corresponding author. E-mail:
[email protected]. † Okayama University. ‡ Boston College. (1) Angell, C. A. Proc. Natl. Acad. Sci. U.S.A. 1995, 92, 6675. (2) Mohanty, U. AdVances in Chemical Physics; Prigogine, I., Rice, S. A., Eds.; John Wiley & Sons: New York, 1994; Vol. LXXXIX, Chapter 2, pp 89-158. (3) Hodge, I. M. J. Non-Cryst. Solids 1994, 169, 211. (4) Frick, B.; Richter, D. Phys. ReV. B 1993, 47, 14795. (5) Rasmussen, B. F.; Stock, A. M.; Ringe, D.; Petsko, G. A. Nature 1992, 357, 423. (6) Doester, W.; Cusck, S.; Petry, W. Nature 1989, 338, 754. 10.1021/ja011927h CCC: $22.00 © 2002 American Chemical Society
and X-ray scattering,5 Mossbauer spectroscopy,7 and computer simulations.8 A number of interpretations have been advanced to explain this phenomenon. The change in slope of the DW factor has been attributed to a crossover from harmonic to anharmonic dynamics,1,9,10 the notion of soft phonons,11 and the onset of inelastic processes and further elucidated by mode coupling analysis.4,12 The DW factor for crambin crystals has been measured in a range of temperatures between 100 and 240 K with considerable accuracy by high-resolution (between 0.67 and 0.89 Å) X-ray crystallography.13 These measurements indicate a change in slope of the DW factor around 200 K. The slope of the disordered side chains of crambin is indistinguishable from that of water below 200 K.13 Furthermore, the slope of the DW factor for oxygen atoms against temperature in hexagonal ice (ice Ih) is close to that of the ordered crambin atoms above 200 K, suggesting that the ordered interior of protein behaves dynamically like that of ice.13,14 The observed glass transition in proteins depends on hydration, i.e., the solvent content.4-6,8-10,12,13 (7) Achterhold, K.; Keppler, C.; Burck, U. V.; Potzel, W.; Schindelmann, P.; Knapp, E. W.; Melchers, B.; Chumakov, A. I.; Baron, A. Q. R.; Ruffer, R. et al. Eur. Biophys. J. 1996, 25, 43. (8) Kuczera, K.; Smith, J.; Karplus, M. Proc. Natl. Acad. Sci. U.S.A. 1990, 87, 1601. (9) Doster, W. In Hydration Process in Biology; Bellisent-Funel, M. C., Teixera, J.; IOS Press: Amsterdam, The Netherlands, 1998. (10) Tilton, R. F.; Dewan, J. C.; Petsko, G. A. Biochemistry 1992, 31, 2468. (11) Buchenau, U.; Zorn, R. Europhys. Lett. 1992, 18, 523. (12) Demmel, F.; Doster, W.; Petry, W.; Schulte, A. Eur. Biophys. J. 1997, 26, 327. J. AM. CHEM. SOC. 2002, 124, 8085-8089
9
8085
ARTICLES
Tanaka and Mohanty
In this paper, we investigate by molecular dynamics (MD) simulation the temperature dependence of the DW factor for hexagonal ice (ice Ih) with 25 different proton-disordered configurations. The details of the simulation and the methodology used to evaluate the DW factor are described in the Section 2. In Section 3, we present the results of the DW factor calculated over a wide range of temperature for the TIP4P15 ices within the framework of classical and quantum harmonic approximation and of the classical MD simulations that include anharmonic interactions. 2. Theoretical Methods The temperature dependence of the DW factor of ice Ih is investigated for 25 different proton-disordered configurations that were generated by a technique discussed in ref 16 and compared with experimental data.13,14 Each configuration is composed of 288 water molecules with no net dipole moment. The intermolecular interaction is mainly described by the TIP4P potential15 since it can reproduce well various properties of not only liquid water but also low-pressure ices.17-19 The apparent weakness of this potential is that it is pairwise additive although it includes higher body interactions in an effective manner. Therefore, it is important to examine DW factors by using other model potentials, which are simple but different in functional form and/or geometry. The SPC/E is one of the best potential models with only three interaction sites.20 Its negative charge is placed on the oxygen atom and the HOH bond angle is 109.5°. The CC potential is a revised version of its predecessor having same functional form in which only the pair interaction is taken into account.21 We have also examined those two potentials. A series of MD simulations containing 288 water molecules in the canonical ensemble at constant temperature, with no volume fluctuations,22 are performed for 25 proton-disordered structures of ice Ih. Temperature is controlled by the Nose method.22 The density at a given temperature is fixed to a constant value, which is determined by the preceding Monte Carlo simulation at constant temperature-pressure.23 This choice of ensemble is appropriate to calculate the DW factors since molecular displacements arise solely from the vibrational motions and do not include those corresponding to motions of the lattice sites by the volume fluctuation. The temperature is set to a range from 150 to 250 K. The density decreases with increasing temperature in the temperature range examined here. We are not interested in a range of very low temperatures where the thermal expansively of ice is negative. This characteristic of the negative thermal expansion coefficient has been attributed to quantum effects that are inherent at low temperatures.24 The time step for integration of the equations of motion is 4 × 10-16 s. Each production run is 15 ns except for an initial configuration (30 ns). The DW factor originating from the harmonic vibrations is calculated by the following statistical mechanical treatment. The system is described by a collection of harmonic oscillators with reference to the corresponding local energy minimum structure. The Hamiltonian of the ice is assumed to be harmonic as
H)
∑ (n + i
i
/2)hνi + Uq
1
(1)
(13) Teeter, M. M.; Yamano, A.; Stec, B.; Mohanty, U. Unpublished, 2001. (14) Kuhs, W. F.; Lehmann, M. S. J. Phys. Chem. 1983, 87, 4312. (15) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (16) Okabe, I.; Tanaka, H.; Nakanishi, K. Phys. ReV. E 1996, 53, 2638. (17) Ohmine, I.; Tanaka, H. Chem. ReV. 1993, 93, 2545. (18) Poole, P. H.; Essmann, U.; Sciortino, F.; Stanley, H. E. Phys. ReV. E 1993, 48, 4605. (19) Tanaka, H. J. Chem. Phys. 1996, 105, 5099. (20) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F. J. Phys. Chem. 1987, 91, 6269. (21) Carravetta, V.; Clementi, E. J. Chem. Phys., 1984, 81, 2646. 8086 J. AM. CHEM. SOC.
9
VOL. 124, NO. 27, 2002
where Uq is the potential at the minimum, h is Planck’s constant, νi is the frequency of the i-th mode, and ni is the quantum number for the i-th vibrational mode. Since the coordinates ∆r in real space are related to the collective coordinate q by
∆r ) m-1/2U†q
(2)
We use S, the square root of the mass matrix m, which may contain finite off-diagonal elements and denote the ij-element of U† and S by uij and sij, respectively, to obtain
〈∆ri2〉 )
∑
j,k,l sijsikujlukl〈ql
2
〉
(3)
The angular brackets denote thermal average. The average of the l-th normal coordinate 〈ql2〉 is defined by
〈ql2〉 ) Tr(F(H)ql2)
(4)
where F is the density matrix of the entire system. Since all the modes are completely independent, 〈qlqm〉 ) 0 for l * m, and we obtain
〈ql2〉 )
∑
nl
F(nl)〈nl|ql2|nl〉
(5a)
F(nl) ) exp(-nlβhνl)[1 - exp(-βhνl)]
(5b)
〈nl|ql2|nl〉 ) h(nl + 1/2)/4π2νl
(5c)
In the above equations, β is 1/kBT where kB is the Boltzmann constant. In classical statistical mechanics, the mean-squared displacement of the l-th vibrational mode is expressed in terms of its frequency νl: 〈ql2〉 ) kBT/(2πνl)2; substituting this relation in eq 5a leads to the desired expression for the mean-squared displacement of the individual molecules. Several comments are in order regarding our methodology. First, the configurations generated by molecular dynamics simulations are quenched to the local energy minimum structures. Second, the normalmode analysis is obtained by diagonalizing the corresponding mass weighted second derivative of the potential function but evaluated at the local energy minimum structures unless otherwise mentioned. Third, we calculate the mean-square displacements for just the oxygen atoms, since this is appropriate for comparison with X-ray diffraction data. Since eq 3 gives the mean-square displacements of the center-of-mass coordinates and Euler angles separately, these are transformed into the atomic displacements, ∆r, assuming the orientational displacements are small, which is always the case for water molecules.
3. Results and Discussion
Figure 1a shows the DW factor as a function of temperature for TIP4P ices. We observe that there is a change in slope of the DW factor around 200 K (or a little higher). The larger value for the DW factor does not mean melting of the ice. This is because at temperatures below 250 K, the mean-square displacement does not increase as the simulation proceeds, implying that the system remains a solid. The melting point of TIP4P water is known to be around 240 K, approximately 30 K lower than that of real water.23,25 However, the TIP4P ice is stable even at 250 K when started with an ice configuration. This is due to the hysteresis effect expected for a first-order transition. Ice structures are always broken into liquid water above 260 K. Thus, the break in the DW factor around 200 K is a real one and could correspond to the break observed in experimental (22) (a) Nose, S. Mol. Phys. 1984, 52, 255. (b) Nose, S. J. Chem. Phys. 1984, 81, 511. (23) Gao, G. T.; Zeng, X. C.; Tanaka, H. J. Chem. Phys. 2000, 112, 8534. (24) Tanaka, H. J. Chem. Phys. 1998, 108, 4887. (25) Karim, O. A.; Haymet, A. D. J. Chem. Phys. Lett. 1987, 138, 531.
Debye−Waller Factor of Hexagonal Ice
ARTICLES
Figure 2. The time evolution of the potential energy of local minima for 30 ns run at three different temperatures.
Figure 1. The mean-squared displacement (MSD) in Å2, i.e., the DebyeWaller factor for hexagonal ice averaged over 25 proton-disordered configurations, is plotted against temperature. The solid line in part a denotes the DW factor obtained from MD simulation with error bars, the dotted line denotes the corresponding DW factor in the classical harmonic approximation, while the dot-dash is the quantum harmonic approximation. The harmonic DW factors are magnified in part b.
data.12,13 The small difference in temperature between experimental and simulation data may be due to the smaller system and the constraint on the simulation, as well as minor deficiencies in the intermolecular interactions used. The DW factor has been calculated over a wide range of temperatures for TIP4P ices within the framework of classical and quantum harmonic approximation discussed above. Note that the use of eq 5 for a fixed density is justified since we find that small differences in density by thermal fluctuations do not give rise to noticeable differences in the DW factors. As Figure 1b shows, quantum treatment of the DW factor exhibits a fairly different characteristic near 0 K; the DW factor levels off with decreasing temperature and it remains a finite value on approaching 0 K. This is due to contributions from the zero-point vibrational energy in the low-temperature range. The classical DW factor, however, approaches the quantum behavior at temperatures higher than 150 K. Two comments deserve mentioning. First, the experimentally observed changes in the DW factor cannot be explained within the framework of the harmonic approximation. Second, the small difference in the predictions for the DW factor between the quantum and the classical treatment in the high-temperature region suggests that a classical MD simulation that includes
anharmonic interactions is likely to be adequate in examining the origin and the break in the DW factor. We have examined time evolution of one proton-disordered ice structure in more detail by performing MD simulation for 30 ns at 270, 250, 230, and 200 K. The generated configurations are recorded at an interval of 20 ps. These configurations are subjected to the steepest descent energy minimization. The local energy minimum structures so obtained are hereafter called Q-structures and correspond to the potential energy basin centers in configuration space. At temperatures such as 270 K, which is much higher than the melting point of TIP4P ice, the ice structures melt into liquid water in a few tens of picoseconds. Thus, no attention will hereafter be paid to the systems at temperatures higher than 250 K. The potential energy of the Q-structures at 250 K is plotted as a function of time in Figure 2. The potential energy gradually decreases and arrives at a limiting value for 30 ns. This is caused by rearrangement of water molecules via structures containing some defects in hydrogen bonds as discussed later. This is interesting because it may correspond to a process where proton ordering increases although the ordering is only partial. Seki and co-workers have experimentally detected freezing of the entropy decrease in ice Ih with decreasing temperature26 at around 100 K. It is expected to take more than several hundred days for the complete ordering. However, the local rearrangement at high temperature takes place in a short time period perhaps as quickly as a few nanoseconds as shown in the present study. This fast relaxation at high temperature has been achieved because it is accompanied by a fairly large energy decrease and the number of available configurations is limited for a system containing a rather small number of molecules. Even after reaching the low-energy limiting state, the system visits various Q-structures whose energy levels are a little higher than the limiting state. This limiting state is not retained for a long time and the system finds further higher energy states. In the process of lowering the potential energy, it is found that water molecules sometimes jump to other lattice sites via intervening structures in which there are a few defects in the location of molecules and the hydrogen bond network. These intervening structures have water molecules that have less than four hydrogen-bonded neighbors, and a few molecules are located at off-lattice site positions. Thus, the total number of (26) Haida, O.; Matsuo, T.; Suga, H.; Seki, S. Proc. Jpn. Acad. 1972, 48, 237. J. AM. CHEM. SOC.
9
VOL. 124, NO. 27, 2002 8087
ARTICLES
Tanaka and Mohanty
Figure 3. A snapshot of the ice configuration at 250 K. The basic cell has been divided into three sub-boxes, each of which contains 96 molecules.
hydrogen bonds is slightly less than 2N, where N is the number of water molecules in ice. Yet, these structures are mechanically stable and therefore belong to Q-structures. These structures are not a family of the proton-disordered ice forms satisfying the “ice rule” completely.27 That is, an ice form having defects in the hydrogen bond network can constitute a stable point in a limited region in configuration space. Some examples of ice structures containing a few defects are shown in Figure 3. The number of defects is very small and most of the molecules are (27) Eisenberg, D.; Kauzmann, W. The Structure and Properties of Water; Oxford University Press: London, UK, 1969. 8088 J. AM. CHEM. SOC.
9
VOL. 124, NO. 27, 2002
located at the lattice sites. Thus, the whole system can be viewed as a crystalline form with a few defects. These local energy minimum structures with some defects are not unique to the present potential, but are observed for CC and SPC/E models of water. It is also found that only a small number of water molecules are involved when a transition from a Q-structure to another one is made. The system begins to leave the initial potential basin at 250 K visiting many potential basins with no defects as well as potential basins with defects. It should be noted that the potential energy gradually decreases as time elapses. The intervening basins with defects are rather easily passed by at this temperature. After reaching the apparent limiting energy, the system sometimes visits higher energy states. The local minimum potential energy curves for five proton-disordered structures for a 15 ns run are plotted in Figure 4. A similar potential energy fluctuation to that plotted in Figure 2 and its gradual decrease with time are seen for all five proton-disordered structures. When the temperature is decreased to 230 K, molecules still move frequently from the initial potential basins to adjacent ones, which may include some defects in the hydrogen bond network but belong to one of potential basin centers. Examination of the MD trajectories at 230 K shows that the system undergoes many transitions (but less frequently than at 250 K) to other configurations either satisfying or violating the “ice rule”. The system visits nearby potential basins for a while with defects or without defects at this temperature, and then may return to the initial basin. The system has sufficient kinetic energy to surmount the activation energy barrier to lead to an adjacent basin with defects but does not have to go to the lower energy state within 30 ns. The larger DW factor at 250 K than that at 230 K is attributed not to the systematic decrease in potential energy but to the more frequent transitions. There is almost no exchange of hydrogen-bonded partners at 200 K as shown in Figure 2; it occurs a few times during the whole simulation run. Thus, molecular motions are restricted mostly to vibrational motions, which dominate the DW factor at this temperature. This does not, however, exclude the very slow proton-ordering motions below 200 K. In the solid state, the mobility of the molecules is generally zero except for very infrequent jump motions to other lattice sites. However, the jump-like motions from one lattice site to another in ice Ih lead to the larger DW factor above 200 K. If this is the case, such a motion is accompanied by a passage through the energy barrier region in configuration space. Thus, the number of imaginary vibrational modes along the MD trajectory at a given instant is an indicator of the transition between two energy basins, although some imaginary modes certainly do not lead to other energy basins.28-30 To see this, we have performed (instantaneous) normal-mode analysis of the system generated by MD simulations and identified the imaginary frequency modes. Figure 5 is a plot of the percentage of imaginary frequency modes versus temperature. The percentage of the imaginary frequency modes decreases almost linearly with a decrease in temperature. However, there is a change in the slope around 200 K (the temperature region where the break in the DW factor was observed). The change in the slope is ascribed to the onset of an appreciable number of the jumplike motions to other mechanically stable structures. (28) Gezelter, J. D.; Rabani, E.; Berne, B. J. J. Chem. Phys. 1997, 107, 4618. (29) Keyes, T.; Li, W.; Zurcher, U. J. Chem. Phys. 1998, 109, 4693. (30) Gezelter, J. D.; Rabani, E.; Berne, B. J. J. Chem. Phys. 1998, 109, 4695.
Debye−Waller Factor of Hexagonal Ice
ARTICLES
Figure 4. The time evolution of the potential energy of local minima for 15 ns run at 250 K for each of the five proton disordered structures labeled (a)-(e).
Figure 5. The percentage of imaginary frequency modes is plotted versus temperature. Normal-mode analysis is done for configurations generated by MD simulations.
4. Summary
In this paper we have carried out classical MD simulations at constant temperature with no volume fluctuations to examine the temperature dependence of the DW factor of ice Ih with 25 different proton-disordered configurations. We observe a break in slope of the DW factor around 200 K, which cannot be explained within the framework of either classical or quantum harmonic approximation. Analysis of the local energy minimum structures shows that water molecules sometimes jump to other lattice sites via intervening mechanically stable structures in which there are a few defects in the location of molecules and the hydrogen bond network. These intervening ice structures
belong to Q-structures but are not a family of the protondisordered ice forms satisfying the “ice rule” completely, that is, they contain several water molecules having less than four hydrogen-bonded neighbors, and some molecules are located in regions other than the lattice sites, yet almost all the molecules occupy their lattice sites. The break in the DW factor indicates the onset temperature above which jump-like motions between those Q-structures dominate. We do not observe frequent transitions from one configuration to another below 200 K. In this temperature range, the DW factor is small and is due mainly to the vibrational motions. Even though such a transition is very slow compared with our simulation time scale, it is experimentally observed down to 100 K. The results presented here indicate there is considerable evidence supporting that the mechanism just above 100 K is the same as that at the temperature where a change in slope of the DW factor is observed. The change in slope of the DW factor depends on the frequency of the transition (a transition that is not suppressed at low temperatures, although its occurrence is a rare event). It is confirmed from an additional MD simulation that such transitions are also observed for a larger system containing 2304 molecules (8 times larger than the systems we have examined above). We note, in passing, this break may be associated with the two-step sequence of hydrogen bond breaking, although the modeling of water is different.31 Acknowledgment. This work was supported by JSPS, NIH, and the Japan Ministry of Education. The authors are grateful to Prof. Yamamuro for stimulating discussions. JA011927H (31) Sonwalkar, N.; Yip, S.; Sunder, S. J. Chem. Phys. 1994, 101, 3216. J. AM. CHEM. SOC.
9
VOL. 124, NO. 27, 2002 8089
