19634
J. Phys. Chem. B 2004, 108, 19634-19639
Structure of Ice-VII and Ice-VIII: A Quantum Mechanical Study† Jer-Lai Kuo* and Michael L. Klein Department of Chemistry and Center for Molecular Modeling, UniVersity of PennsylVania, Philadelphia, PennsylVania 19104 ReceiVed: April 22, 2004; In Final Form: June 30, 2004
Ab initio quantum mechanical methods have been used to investigate the structures of Ice-VII and Ice-VIII. In the case of proton-ordered Ice-VIII, accurate structural data enable us to benchmark the methodology. The disorder in H-bond orientations of Ice-VII causes difficulties in understanding its structure from either an experimental or a theoretical viewpoint. In this work, we explore the consequence of randomness of H-bond arrangements in Ice-VII by examining all ice-rule-allowed configurations in a unit cell of 16 water molecules. Noticeable relative displacement () between the two sublattices, previously thought to exist only in Ice-VIII, was found in many of these configurations. We find a correlation between the energy of a particular configuration of Ice-VII and , an energy gap of 1.5 mhartree that separates a unique ground state from all the others. Our results have proven that ab initio calculations can be used to resolve ambiguities in interpreting neutron diffraction data of Ice-VII.
I. Introduction The rich and complicated behavior of water has long fascinated scientists. For pressures above 2GPa, the phase diagram is dominated by two phases of ice, namely, Ice-VII and Ice-VIII.1-3 The structures of these phases differ only in H-bond orientations and slight distortion in lattice constants. Ice-VIII is known as a proton-ordered phase with a tetragonal unit cell, while Ice-VII has a cubic unit cell and H-bond orientations that are thought to be random subject to the ice rules.4-6 In the P-T phase diagram, the phase boundary between IceVII and Ice-VIII contains many intriguing transitions of very different nature.1,2 Under a wide pressure range (between 2 and 10 GPa), a thermally induced transition between high-temperature disordered Ice-VII and low-temperature ordered Ice-VIII is characterized by a pressure-insensitive transition temperature (∼260 K). Upon further increase in pressure, the phase boundary falls rapidly at about 60 GPa. This pressure-induced transition (from low-pressure Ice-VIII to high-pressure Ice-VII) was predicted by Schweizer and Stillinger,7-9 and they attributed its cause to the quantum tunneling of protons. Furthermore, recent experiments have discovered the existence of a lowtemperature metastable form of Ice-VII (LT-IceVII)10 inside the stable region of Ice-VIII and both thermal and pressure induced transitions between LT-IceVII and Ice-VIII have been suggested.11 A better understanding of these fascinating phenomena in the high-pressure phases of ice hinges upon a detailed characterization of the structures, which is the focus of this work. Despite the chronological order that the existence of Ice-VIII was confirmed after that of Ice-VII, the latter remains less well understood. In the 1980s, neutron diffraction experiments gave solid evidence to validate the proton-ordered structures of IceVIII.12,13 The same experimental setup, however, encountered difficulties in data analysis caused by the disorder in H-bond arrangements in Ice-VII and yielded either surprisingly short oxygen-hydrogen distances (ROH ∼ 0.89 Å in ref 12 and †
Part of the special issue “Frank H. Stillinger Festschrift”.
Figure 1. Structure of Ice-VII. Oxygen atoms are represented by larger spheres, and hydrogen atoms are represented by smaller yellow spheres. For easy reference, oxygen atoms in different sublattices are colored differently. Hydrogen atoms are placed in the middle between two oxygen atoms to reflect the 50% occupancy. The cubic cell is drawn to emphasize the diamondlike form of a sublattice. It should be noted that the origin of the unit cell shown here is different from that used in experimental structure refinements (both single-site model and multisite model) in which the origin resides at the inversion center. The green arrows near the central oxygen indicate the symmetry-allowed directions of distortion described in the multisite model. The blue arrows labeled by /2 give some possible directions of the relative displacement of the two sublattices.
∼0.943 Å in ref 14) or improbable H-O ˆ -H angles (∠HOH: 98.2-114° in ref 14). In those studies, the “single-site” model used in structural refinement assumes the oxygen atoms are located on the perfect diamondlike lattice (often denoted as 1/ ,1/ ,1/ along the [111] axis) and the hydrogen atoms take a 4 4 4 special position (x, x, x also along [111] axis, with x ∼ 0.42). Details of this model are discussed in the caption of Figure 1. Recently, a “multisite” model has been proposed to reexamine the neutron diffraction data,15 in which oxygen atom positions are allowed to deviate from their single-site locations to six other sites allowed by symmetry (see Figure 1). This model attributes the deviation of oxygen from its mean position solely to the distortion occurring within a sublattice, and the amount of deviation (δ) is found to be ∼0.1 Å. The dashed lines in Figure 1 demonstrate how the multisite model succeeds in resolving the difficulties encounted in the single-site model. But, the resulting structure from the multisite model indicates
10.1021/jp0482363 CCC: $27.50 © 2004 American Chemical Society Published on Web 09/10/2004
Structure of Ice-VII and Ice-VIII that Ice-VII contains a mixture of two kinds of H-bonds with different lengths (∼2.75 and ∼2.93 Å). Another alternative to the multisite model is to consider the relative displacement between the two sublattices. This kind of displacement (shown by the arrows labeled as /2 in Figure 1) has been known to exist in Ice-VIII along its c axis. As a result of the symmetry of Ice-VII (space group Pn3m), the overall average of has to be 0, but locally a nonzero value of is possible. Throughout this manuscript, we will use δ to stand for the averaged amount of distortion in oxygen positions described in the multisite model and use to represent the amount of relative displacement between two sublattices. Two more elaborate models were proposed by Kuhs et al. to account for the disorder in oxygen positions in analyzing their powder neutron diffraction data.12 Unfortunately, using experimental means, it has proven difficult to separate the contribution from δ, , and other sources.12,14,16 Hence, one of our aims in this work is to use ab initio calculations to shed some light on these ambiguities. Theoretical modeling of Ice-VII is also hampered by the random orientation of the H-bonds. In this work, we explore all possible disordered proton arrangements in Ice-VII allowed in a unit cell of 16 water molecules. As a result of the high symmetry of Ice-VII, only 52 symmetrically distinct structures are needed for ab initio calculations compared to a total of 8100 ice configurations permitted in the above-mentioned unit cell. The remainder of this paper is structured as follows. In section II we will briefly discuss the enumeration of all ice-rule-allowed configurations and the generation of their initial geometries. We will provide a detailed description on our geometry optimization procedure and use the better understood structural data of IceVIII as a benchmark for our ab initio calculations. In section III, we will present the energetics and structural data of all 52 ice configurations, followed by a detailed examination on the contributions of and δ to the deviation of oxygen atoms from their perfect lattice. This paper concludes with a discussion of our results and its relevance to existing models used in interpreting neutron diffraction data. II. Methodologies A. Enumeration of Ice-Rule-Allowed Configurations. IceVII consists of two inter-penetrating Ice-Ic sublattices and the H-bond networks in these two sublattices are independent. The smallest cubic unit cell of Ice-VII (shown in Figure 1) consists of only two water molecules with a lattice constant of 3.337 Å at 2.4 GPa and 263 K. This smallest unit cell (denoted as 1 × 1 × 1) is so small that only one ice-rules-allowed H-bond arrangement exists. Hence, in this work, we focus on a unit cell that is 8 times (extended in all three directions, hence, denoted as 2 × 2 × 2) larger than the 1 × 1 × 1 unit cell. The 2 × 2 × 2 unit cell along with some ice-rules-allowed configurations are shown in Figure 2, and we will use 2 × 2 × 2 unit cell through this manuscript unless specified otherwise. The fact that the H-bond networks in the two inter-penetrating Ice-Ic sublattices are independent means that the number of possible proton configurations in Ice-VII will be the square of the number of possible proton configurations in Ice-Ic. The latter has been studied by Lekner who found 90 possible proton arrangements of Ice-Ic in a unit cell of eight water molecules.17 Our straightforward enumeration of all ice-rules-allowed structures in the 2 × 2 × 2 unit cell of Ice-VII gives 8100 possible proton configurations.
J. Phys. Chem. B, Vol. 108, No. 51, 2004 19635
Figure 2. Four representative ice-rules-allowed H-bond arrangements in the 2 × 2 × 2 unit cell. The ground-state structure of Ice-VII is shown in part a, and it has an H-bond arrangement (anti-ferroelectric) identical to the structure of Ice-VIII. The relative displacement of the two sublattice is not visible from this figure, but its direction is indicated by the arrow. Part b shows the second most stable structure, part c is the structure exhibiting the maximal deviation in ∠OOO, and part d shows a ferroelectric structure with the highest energy among all 52 ice configurations. Their relative energies are shown in Figure 3 by corresponding letters, and their bond lengths and angles are detailed in Table 2.
Examining all 8100 different configurations via ab initio methods is not only prohibitively expensive but also unnecessary. As a result of the high symmetry of the lattice, most of the 8100 configurations are symmetrically related to one another. For all “chemically identical” (that is, symmetrically related) structures, only one representative is needed for detailed examination. The elimination of symmetrically related configurations can be quickly done by employing the graph invariance algorithm that has been detailed in previous works.18,19 The number of “chemically distinct” configurations in our 2 × 2 × 2 unit cell is 52, and in the following we will examine all of them by ab initio methods. B. Geometry Optimization. The initial oxygen-oxygen distances (ROO) are taken to be equal to x3/4a2×2×2, where a2×2×2 is the lattice constant of the 2 × 2 × 2 unit cell. No prior distortion between the two sublattices is assumed. The initial positions of the hydrogen atoms are located along the straight line between two oxygen atoms, and the initial oxygenhydrogen bond length (ROH) is assumed to be one-third of ROO. In the geometry optimization all degrees of freedom are allowed to relax, except the lattice constants, which are fixed to the experimental values (a2×2×2 ∼ 6.674 Å at 2.4 GPa and 263 K). All ab initio calculations are carried out using the CarParrinello methods20,21 based on the Becke-Lee-Yang-Parr (BLYP) functional, using norm-conserving pseudo-potentials (Troullier and Matrins form)22 with a plane-wave cutoff of 70 Ry, and considering only the Γ point of the Brillouin zone. These are the standard conditions empolyed in numerous investigations of aqueous systems.23-28 C. Structure of Ice-VIII: A Benchmark Calculation. Even though the above-mentioned procedure and parameters have been well tested, we present in this section the results from our geometry optimization on Ice-VIII in which detailed structural data from neutron diffraction12,13 and other ab initio methods29-31 are available for benchmarking purposes.
19636 J. Phys. Chem. B, Vol. 108, No. 51, 2004
Kuo and Klein
TABLE 1: Comparisons of the Bond Distances (Angstroms) and Angles (Degrees) in Ice-VIIIa RHB OO Rnon OO ROH ∠HOH ∠OHO
expt
calcd (BLYP)
calcd (HF)
calcd (PW)
2.879(1) 2.743(9) 0.9685(71) 105.61(1.06) 178.25(70)
2.8797 2.7724 0.9877 105.74 177.26
2.88 2.74 0.9522 107.2 179.6
2.875 2.756
a
For the ROO distance, “HB” stands for H-bonded and “non” stands for non-H-bonded. The experimental data shown in the first column is from neutron diffraction in Table 3 of ref 12. Note that the three ab initio calculations were carried out with different methods or parameters (HF ) ref 30, PW ) ref 31, and BLYP ) this work).
Figure 3. Potential energy (solid circles) of all 52 symmetrically distinct configurations of Ice-VII are shown in the left figure along with the average energy (dashed line). The unit of the energies is mhartree, and the letters a-d indicate the energies of the corresponding structures in Figure 2. The right figure shows the heat capacity in the reduced and dimensionless form. The average energy and heat capacity clearly show that a “phase transition” occurs at ∼150 K.
The lattice constants of Ice-VIII were fixed at the experimental values (a ) 4.656 Å and c ) 6.775 Å at 2.4 GPa and 10 K).12 Comparison of selected bond distances (in angstroms) and angles (in degrees) are summarized in Table 1, from which it is obvious that all three ab initio methods give a reasonable description of the structure of Ice-VIII. We estimate by taking the difference between the centers of mass of the two oxygen sublattices from the optimized geometry. Our calculation shows a displacement of ∼ 0.19 Å, which is in excellent agreement with the commonly accepted value of ∼ 0.2 Å.16 To complete the comparison, we also calculate the H-bonded O‚‚‚O‚‚‚O angles (denoted as ∠OOO). We adapt the convention of Kuhs et al.12 and categorize all ∠OOO into three classes, namely, O‚‚‚H-O-H‚‚‚O, O‚‚‚H-O‚‚‚H-O, and O-H‚‚‚O‚‚‚H-O. In the following, we shall use ∠OOOA, ∠OOOB, and ∠OOOC to represent these three kinds of ∠OOO, respectively. The experimental (ab initio) values [∠OOOA ∼ 107.92° (107.88°), ∠OOOB ∼ 110.25° (110.25°), and ∠OOOC ∼ 107.92° (107.88°)] again demonstrate the excellent performance of the methodology we employ herein. A larger plane-wave cutoff (up to 100 Ry) and an increase in k-point mesh via the Monkhorst-Pack scheme (up to 32 special k points) were carried out to address the accuracy of the methodology given above. It was found that the geometry of Ice-VIII is insensitive to the change in plane-wave cutoff and k-point sampling. Hence, we conclude that the procedure used herein is sufficiently accurate for the present purposes. III. Results and Discussions A. Ground State and Energy “Gap”. The energies of all 52 symmetrically distinct ice configurations of Ice-VII are shown in Figure 3, where the zero of the energy scale is set to the energy of the ground-state structure (see Figure 2a). The energy difference between the most and least stable configurations is about 5 mhartree, and from Figure 3, it is clear that there is a “gap” (of roughly 1.5 mhartree for an unit cell of 16
water molecules) that separates the ground-state structure from the other 51 ice configurations. To demonstrate the physical relevance of this energy gap, we calculate “thermal properties” estimated via the following configurational partition function, M
QN )
fie-βE ∑ i)1
i
(1)
where Ei is the potential energy of the ith ice configuration, fi is the number of symmetry-related configurations that are represented by the ith configuration, and M is the total number of symmetrically distinct representatives. The temperature dependence of the average energy (〈E〉) and the heat capacity (Cv), derived from the configurational partition function, are shown in Figure 3. It is important to note that results shown in Figure 3 are only rough estimates, because the “ensemble” used here consists of only the 8100 ice configurations permitted in a 2 × 2 × 2 unit cell. The true partition function (and, thus, the thermal properties) would, of course, depend on the entire energy landscape, which includes ice configurations with random H-bond networks permitted only in larger unit cells. This direction can be pursued by using graph invariants methods19 and will be discussed in a future publication. Nevertheless, as indicated in the heat capacity and average energy there is clearly a “phase transition” at about 150 K. This temperature is in reasonable agreement with the phase transition temperature of the thermally induced transition between LTIceVII and Ice-VIII.11 The ground-state structure deserves some remarks. By examining its H-bond arrangement, we found that this groundstate structure has an H-bond network identical to the structure of Ice-VIII. Like the structure of Ice-VIII, the ground-state structure of Ice-VII also reveals a significant amount of displacement of the two sublattices ( ) 0.14 Å). Opposite to Ice-VIII, which has different averaged values for the three kinds of ∠OOO, the two sublattices of this ground state structure maintain a nearly perfect tetrahedral arrangement. In this section, we will just mention that the nonzero and extremely small value of δ are by no means the unique signature of the groundstate structure of Ice-VII and leave the further discussions on the structures of all 52 ice configurations to section 3.2. From our studies on the energetics, we also notice an interesting correlation between and the energy (shown in Figure 4). In many theoretical calculations of ice (most of them based on either empirical force fields or simple electrostatic models), various “parameters” (such as the percentage of trans/ cis conformations between a pair of water molecules32-36 and the degree of ferro-electricity)17,37 were proposed to link with the energetics. However, none of these parameters have been proven to be effective and none of these correlations have been tested in other phases of ice. In Figure 4, the correlations between the energetics and the above-mentioned parameters are shown. It is clear that only exhibits a strong correlation with the energy. Plain electrostatic models, which in general favor anti-ferroelectric structures, cannot fully explain the trends seen in Figure 4b because structures with no net dipole moments can have very different energetics. For a water dimer, the concept that a trans conformation is energetically more favorable than a cis structure is commonly accepted, but in the solid state its validity is not very convincing. For example, in Ice-Ih this notion has been debated for many years in the literature.32-36 From Figure 4c, we demonstrate that there is no clear dependence between ntrans and energy and, in fact, the most and
Structure of Ice-VII and Ice-VIII
J. Phys. Chem. B, Vol. 108, No. 51, 2004 19637 TABLE 3: Average Bond Distances and Angles of Ice-VII from ab Initio Calculationsa X RHB OO Rnon OO ROH
∠HOH ∠OHO ∠OOOA ∠OOOB ∠OOOC
calcd (BLYP)
〈σX〉i
〈σX〉all
2.8907 2.8086 0.9885 105.38 177.73 107.77 109.45 111.14 0.0514
0.018 0.081 0.005 0.26 2.17 1.92 1.37 2.13
0.010 0.014 0.004 0.286 1.057 0.921 0.936 0.969 0.037
a Notations and units are the same as those in Table 1. The second column (〈σX〉i) indicates the maximal standard deviations in an ice configuration, and the third column (〈σX〉all) gives the overall standard deviation averaged over all 52 ice configurations. When comparing the results in Table 1 and Table 2, note that different lattice constants were used, which lead to different averaged values of ROO.
Figure 4. Correlation between energy (in mhartree) and some structural parameters of Ice-VII. is the relative displacement of the two sublattices (in angstroms), |µ| is the magnitude of the bond dipole (calculated on the basis of geometrical criteria with arbitrary units), and ntrans is the number of trans conformations.
TABLE 2: Bond Lengths and Angles of Four Configurations (Shown in Figure 2) of Ice-VII Calculated by ab Initio Methodsa calcd (BLYP) RHB OO Rnon OO ROH
∠HOH ∠OHO ∠OOO
a
b
c
d
2.890(0) 2.892(81) 0.9881(0) 106.11(0) 177.44(0) 109.47(0) 0.140
2.890(8) 2.892(72) 0.9883(4) 105.71(26) 177.49(89) 109.47(83) 0.116
2.891(12) 2.891(54) 0.9885(5) 105.32(18) 178.27(81) 109.45(1.64) 0.003
2.890(0) 2.890(0) 0.9887(0) 105.30(0) 176.623(0) 109.46(0) 0.000
a Standard deviations are given in parentheses. Notations and units are the same as those in Table 1.
the least stable structures both contain 100% trans bonds. These observations indicate that the energy of an ice configuration is a delicate balance among H-bond, dipole-dipole, and nonbonded repulsive interactions. Further studies on different phases of ice will likely provide useful hints to a better understanding of these interactions. B. Structure of Ice-VII. In this section, we will first use four selected configurations (shown in Figure 2) as examples to demonstrate the difference in the structures of individual ice configurations. In the second half of this section, we will discuss the averaged bond lengths and angles obtained from all 52 ice configurations we have studied and their relevance to previous interpretations of neutron diffraction data. 1. Four Selected Configurations. The H-bond networks of four selected configurations are shown in Figure 2, and their individual bond lengths and angles are summarized in Table 2. A quick comparison of Table 1 and Table 2 suggests that the structures Ice-VIII and Ice-VII are very similar. Significant difference in ∠HOH and ∠OHO (that is, the O-H‚‚‚O angle) among the four configurations of Ice-VII (see Table 2) are in agreement with the large experimental error bars in Ice-VIII (see Table 1). As we mentioned before, configuration a is the ground-state structure of Ice-VII, and it has very small variations in ∠OOO and the largest . We would like to emphasize that there is no simple correlation (if any) between and δ; the latter is associated with a larger deviation in RHB OO and ∠OOO. In fact, configuration d, the least stable structure, also poses almost
perfect diamondlike sublattices and no noticeable variation in ∠OOO, but it has no relative displacement (that is, ) 0). Unlike configurations a and d, configurations b and c show rich diversities in the bond lengths and angles as one can see from the larger standard deviations. Configuration c is chosen because it has the largest variation in ∠OOO among all 52 allowed configurations, and configuration b is selected to demonstrate that significant amounts of and δ can coexist. It is worth noting that, in both configurations b and c, large variations in ∠OOO are in agreement with the larger variations in RHB OO, but unexpected short and long H-bonds, suggested in the multisite model, were not observed. This set of four ice configurations illustrates the effect of H-bond arrangements on the structure of Ice-VII and serves as a good example to remind us of the necessity to consider many different ice-rules-allowed configurations when modeling proton disordered phases of ice. In most ab intio calculations of ice, only one or two configurations are chosen for geometry optimization and highly symmetrical structures are typically used to save computational cost. As demonstrated by the four configurations we presented here, misleading conclusions could be drawn if only a few configurations were examined. 2. AVeraged Structure. Structural data on Ice-VII are compiled in Table 3. The numbers in the first column are appropriately weighted averages over all distances and angles of the 52 configurations. The standard deviations in the second column give the maximum fluctuation existing in an ice configuration, and these numbers should not be mistaken as error bars in our calculations. Global distributions of the corresponding quantities among the set of 52 ice configurations are shown in Figure 5. By comparing the ab initio results in Table 1 and Table 3, it is evident that the structures of Ice-VII and Ice-VIII are very similar. This finding confirms a traditional wisdom that proton ordered and disordered phases often share a similar underlying structure1-3 (other examples are Ice-Ih and Ice-XI38 and Ice-III and Ice-IX39). This conventional view, however, has been recently challenged by Nelmes et al., and in their multisite model,15 the structure of Ice-VII is different and more complicated than that of Ice-VIII. They concluded that Ice-VII should contain a mixture of two kinds of H-bonds with different lengths (∼2.75 and ∼2.93 Å). Our calculation of RHB OO shows a much smaller variation in all configurations we studied (with the greatest standard deviation being less than 0.012 Å). The distribution of all RHB OO in Figure 5 shows no sign of any bimodal distribution. To further support our finding and rule out the possibility that the absence of a mixture of short and long H-bonds in our calculations is caused by the finite size of
19638 J. Phys. Chem. B, Vol. 108, No. 51, 2004
Kuo and Klein show strong asymmetry. The fact that ∠OOOA tends to be smaller than the tetrahedral angle and ∠OOOC tends to be larger than the tetrahedral angle implies a preferred direction of deviation (against the direction of the dipole moment in a water molecule). IV. Conclusions
Figure 5. Distribution of various bond distances (angstroms) and angles (degrees) obtained from analyzing all 52 configurations in Ice-VII. The averages and standard deviations of each quantity are summarized in Table 3.
the 2 × 2 × 2 unit cell, we re-optimize the four configurations (shown in Figure 2) in a 4 × 4 × 4 unit cell (that is, a unit cell with 128 water molecules). It is found that the unit cell size has no noticeable effect on all the bond distances and angles. Hence, we conclude that the existence of two kinds of RHB OO in Ice-VII is unlikely. Differences in averaged ∠HOH angles exist among the configurations we examined, but we do not find any improbable ∠HOH angles (from Figure 5 we can see that all angles are between 104.8 and 106.2°). This observation is a direct confirmation of the notion that the improbable ∠HOH angle reported in various experiments is due to the deficiency in the single-site model used for structural refinement. In the following, we will examine ∠OOO because of its ability to distinguish the different contributions to the deviation of oxygen positions from δ and . Assuming that both oxygen sublattices remain as perfect diamondlike structures, all ∠OOO would have perfect tetrahedral angles (∼109.47°) even with the presence of a nonzero value of . On the other hand, the multisite model with δ ∼ 0.1 Å would imply that ∠OOO should fluctuate between 106.29 and 112.77° (the angles were obtained under the assumption that 〈RHB OO〉 ) 2.89 Å). In Figure 5, we can see that the majority of ∠OOO angles are located between 108 and 111° (rather than 106.29 and 112.77°), and this distribution indicates a smaller amount of fluctuation in ∠OOO than what the multisite model implies. Nevertheless, a noticeable amount of deviation of ∠OOO from perfect tetrahedral geometry exists and the overall standard deviation in ∠OOO is ∼1.3°. A closer examination of the ∠OOO distributions in Figure 5 reveals more detailed structural information. Although the overall distribution of ∠OOOall has a symmetric form with respect to the perfect tetrahedral angle, individual distributions
We have examined the structures of both Ice-VII and IceVIII by ab initio density functional methods. The comparison of various bond lengths and angles of Ice-VIII with experimental values confirms the utility of density functional methods based on the procedures described in section 2. The same protocol has been applied in numerous previous investigations of aqueous systems, and here we have again demonstrated the transferability of this methodology. The poor understanding of the structure of Ice-VII caused by the randomness in the H-bond orientation is investigated by considering all ice-rule-allowed configurations for a given unit cell. In all, we have studied a total of 52 configurations, which represent a set of 8100 structures permitted in a periodically replicated 2 × 2 × 2 unit cell of 16 water molecules. Significant differences in structural data (bond distances and angles) and energetics were noted among these configurations. This finding emphasizes that it is critical to pay close attention to H-bond orientations to avoid misleading conclusions especially in modeling proton-disordered phases of ice. This is a rather important issue because among all 13 known phases of crystalline ice, 8 of them are proton-disordered and most of the amorphous phases of ice also have disordered H-bond networks.1-3 Our ab initio calculations show that great similarities exist between the structures of Ice-VII and Ice-VIII. Anomalous structural results of Ice-VII (such as surprisingly short ROH, improbably H-O ˆ -H angles, and a mixture of two kind of RHB OO bond lengths) arising from previous interpretations of neutron diffraction data are not found in our calculations. We give direct confirmation that these anomalies were caused by the inadequacy of the simple single-site model or a revised multisite model to include the possible causes of deviation of oxygen from its mean position. The former model assumes that both sublattices remain in a perfect diamondlike structure, and the latter model attributes the deviation solely to distortion within one sublattice (δ). In this manuscript, we also consider the possibility that the relative displacement of the two sublattices () can contribute to the deviation in oxygen position. Both possible causes (δ and ) were critically examined by ab initio data. Although it is difficult to distinguish their contributions by experimental means (not to mention the additional contribution from thermal and zero-point motion, which is not present in our calculation), ab initio calculations offer a fresh view to resolve this ambiguity. We found the existence of both and δ in many of the configurations we examined, and statistically both factors will contribute. The picture emerging from our calculations is that the complexity in the structure of Ice-VII is caused by the disorder in H-bond orientations that yields many ice-rules-allowed structures, each of which possess different bond distances and angles. It is, thus, important to account for distortion in oxygen positions in interpreting neutron diffraction data. Furthermore, our results have also proven that ab initio calculations are useful complements to experimental methods, especially for cases in which high-quality experimental results are difficult to obtain.
Structure of Ice-VII and Ice-VIII Acknowledgment. The work was supported by the National Science Foundation under NSF CHE02-05146, and the calculations were performed at the Pittsburgh Supercomputer Center and NCSA provided by NPACI. References and Notes (1) Hobbs, P. V. Ice Physics; Clarendon Press: Oxford, 1974. (2) Petrenko, V. F.; Whitworth, R. W. Physics of Ice; Oxford University Press: New York, 1999. (3) Chaplin, M. Water Structure and Behavior. http://www.lsbu.ac.uk/ water (accessed 2004). (4) Bernal, J. D.; Fowler, R. H. J. Chem. Phys. 1933, 1, 515. (5) Pauling, L. J. Am. Chem. Soc. 1935, 57, 2680. (6) Rahman, A.; Stillinger, F. H. J. Chem. Phys. 1972, 57, 4009. (7) Stillinger, F. H.; Schweizer, K. S. J. Phys. Chem. 1983, 87, 4281. (8) Schweizer, K. S.; Stillinger, F. H. Phys. ReV. B 1984, 29, 350. (9) Schweizer, K. S.; Stillinger, F. H. J. Chem. Phys. 1984, 80, 1230. (10) Klotz, S.; Besson, J. M.; Hamel, G.; Nelmes, R. J.; Loveday, J. S.; Marshall, W. G. Nature 1999, 398, 681. (11) Song, M.; Yamawaki, H.; Fujihisa, H.; Sakashita, M.; Aoki, K. Phys. ReV. B 2003, 68, 24108. (12) Kuhs, W. F.; Finney, J. L.; Vettier, C.; Bliss, D. V. J. Chem. Phys. 1984, 81, 3612. (13) Jorgensen, J. D.; Beyerlein, R. A.; Watanabe, N.; Worlton, T. G. J. Chem. Phys. 1984, 81, 3211. (14) Jorgensen, J. D.; Worlton, T. G. J. Chem. Phys. 1985, 83, 329. (15) Nelmes, R. J.; Loveday, J. S.; Marshall, W. G.; Hamel, G.; Besson, J. M.; Klotz, S. Phys. ReV. Lett. 1998, 81, 2719. (16) Nelmes, R. J.; Loveday, J. S.; Wilson, R. M.; Besson, J. M.; Pruzan, P.; Klotz, S.; Hamel, G.; Hull, S. Phys. ReV. Lett. 1993, 71, 1192. (17) Lekner, J. Physica B 1997, 240, 263. (18) Kuo, J.-L.; Coe, J. V.; Singer, S. J.; Band, Y. B.; Ojama¨e, L. J. Chem. Phys. 2001, 114, 2527.
J. Phys. Chem. B, Vol. 108, No. 51, 2004 19639 (19) Kuo, J.-L.; Singer, S. J. Phys. ReV. E 2003, 67, 16114. (20) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55, 2471. (21) Marx, D.; Hutter, J. Modern Methods and Algorithms of Quantum Chemistry; John von Neumann Institute for Computing: 2000, p 301 (also available at http://www.theochem.ruhr-uni-biochem.de/go/cprev.html). (22) Troullier, N.; Martins, J. L. Phys. ReV. 1991, 43, 1993. (23) Silvestrelli, P. L.; Parrinello, M. Phys. ReV. Lett. 1999, 82, 3308. (24) Silvestrelli, P. L.; Parrinello, M. J. Chem. Phys. 1999, 111, 3572. (25) Mantz, Y. A.; Geiger, F. M.; Molina, L. T.; Molina, M. J.; Trout, B. L. J. Chem. Phys. 2000, 113, 10733. (26) Chen, B.; Ivanov, I.; Klein, M. L.; Parrinello, M. Phys. ReV. Lett. 2003, 91, 215503. (27) Benoit, M.; Bernasconi, M.; Focher, P.; Parrinello, M. Phys. ReV. Lett. 1996, 76, 2934. (28) Benoit, M.; Romero, A. H.; Marx, D. Phys. ReV. Lett. 2002, 89, 145501. (29) Knuts, S.; Ojama¨e, L.; Hermansson, K. J. Chem. Phys. 1993, 99, 2917. (30) Ojama¨e, L.; Hermansson, K.; Dovesi, R.; Roetti, C.; Suanders, V. R. J. Chem. Phys. 1994, 100, 2128. (31) Tse, J. S.; Klug, D. D. Phys. ReV. Lett. 1998, 81, 2466. (32) Bjerrum, N. Science 1952, 115, 386. (33) Pitzer, K. S.; Polissar, J. J. Am. Chem. Soc. 1956, 60, 1140. (34) Tse, J. S.; Klug, D. D. Phys. Lett. A 1995, 198, 464. (35) Morrison, I.; Li, J.-C.; Jenkins, S.; Xantheas, S. S.; Payne, M. C. J. Phys. Chem. B 1997, 101, 6146. (36) Buch, V.; Sandler, P.; Sadlej, J. J. Phys. Chem. B 1997, 102, 8641. (37) Lekner, J. Physica B 1998, 252, 149. (38) Jackson, S. M.; Nield, V. M.; Whitworth, R. W.; Oguro, M.; Wilson, C. C. J. Phys. Chem. B 1997, 101, 6142. (39) London, J. D.; Kuhs, W. F.; Finney, J. L. J. Chem. Phys. 1993, 98, 4878.
