Prediction of a Phase Transition to a Hydrogen Bond Ordered Form of

Oct 19, 2005 - Department of Chemistry, Ohio State University, 100 West 18th Avenue, ... ice XI, near 72 K for H2O and 76 K for D2O.1-5 The oxygen ato...
0 downloads 0 Views 159KB Size
Download PDF
21040

J. Phys. Chem. B 2005, 109, 21040-21046

Prediction of a Phase Transition to a Hydrogen Bond Ordered Form of Ice VI Chris Knight and Sherwin J. Singer* Department of Chemistry, Ohio State UniVersity, 100 West 18th AVenue, Columbus, Ohio 43210 ReceiVed: July 22, 2005; In Final Form: August 31, 2005

Ice VI is a hydrogen bond disordered crystal over its known region of stability. In this work, we predict that ice VI will transform into a hydrogen bond ordered phase near 108 K, and have identified the likely lowtemperature phase as ferroelectric (space group Cc) with an antiferroelectric structure (space group P212121) close by in energy. Electronic density functional theory calculations provide input to our calculations, which are extended to cells large enough for statistical simulations by using graph invariants. A significant decrease in the configurational entropy is predicted as hydrogen bonds exhibit partial order above the transition, provided that the hydrogen bonds can equilibrate on an experimental time scale. Conversely, partial disorder is predicted at temperatures below the transition. Although some evidence for ordering of ice VI has been observed in experiments, a low-temperature proton ordered phase has not been identified experimentally.

In the phase diagram of water, the liquid phase is surrounded by a series of hydrogen bond disordered ice phases. Ice Ih (ordinary hexagonal ice) is the hydrogen bond disordered solid phase adjacent to the liquid at low pressures. Then, in order of appearance with increasing pressure, ice III, ice V, ice VI, and then ice VII are adjacent to the liquid phase. Several of these phases give way to a hydrogen bond ordered version of themselves at lower temperatures. Ice Ih, when suitably doped with hydroxide, transforms to what is thought to be a ferroelectric phase, ice XI, near 72 K for H2O and 76 K for D2O.1-5 The oxygen atom positions in ice XI are very close to those of ice Ih,6-9 and the transition principally involves selection of one particular H-bond arrangement in the low-temperature phase. Ice III, when cooled at about 1 K per minute or faster, transforms to a metastable H-bond ordered version known as ice IX.10-14 In the 2.1-12 GPa range, the hydrogen bond ordered form of ice VII, known as ice VIII, appears near 0 °C.15,16 At low temperature, the hydrogen bond ordered form of ice V or ice VI has not been clearly identified, although some experiments, as described below, give some preliminary indications of low temperature, fully ordered phases. The purpose of this work is to predict the structure of a low-temperature, H-bond ordered form of ice VI, and the temperature at which the transformation should be found. In accordance with past practice, this new phase should be granted a roman numeral of its own, but for now, following Kamb,17,18 we use primes to distinguish the hydrogen bond ordered phase and call this phase VI′.19 Our predictions are based on a method we have developed that enables us to use data from electronic structure calculations from small unit cells to parametrize a Hamiltonian for a cell large enough for statistical simulations.20-22 The energy difference between various H-bond isomers of ice is known to be rather small.23 However, we have shown that periodic density functional theory successfully predicts the low-temperature structure and location of the phase transition for ice Ih/XI and ice VII/VIII.22 More recent calculations also indicate that the ice III/IX phase transition is predicted well by these methods.24 There is partial experimental evidence for the transformation of ice V and ice VI to hydrogen bond ordered forms, but

complete diffraction data that would identify the hydrogen bond ordered phases of ice V and VI are not available. In 1965, Kamb noticed an X-ray reflection at 77 K that was incompatible with the P42/nmc space group of ice VI and could have signaled the formation of a hydrogen bond ordered version of ice VI.17 Later, he reported neutron diffraction data taken at 100 K on a sample previously equilibrated at high pressure and 77 K, which indicated antiferroelectric ordering.18 In 1976, Johari and Whalley predicted an ordering transition at 47 K to a ferroelectrically ordered state in ice VI based on observed CurieWeiss behavior of the low-frequency dielectric constant.25 Later they concluded that the high-frequency permittivity of ice VI at 0.9 GPa indicated a very slow phase transition occurs in the temperature range 123-128 K, but these experiments did not reveal the structure of the low-temperature phase. Kuhs et al.16 obtained neutron diffraction data on ice VI under temperature and pressure conditions where ice VI is stable, unlike earlier diffraction experiments where the diffraction experiments were performed at ambient pressure on samples recovered from highpressure cells. They took data at 225, 125, and 8 K, but their data were not sufficient to fully determine the structure. They found no evidence of the transformation observed by Johari and Whalley,26 although the transition might be too slow compared to their experimental time scales. The situation with regard to ice V is also unsettled. A brief report of neutron scattering evidence27 and calorimetric data28 indicates that ice V transforms to an ordered structure in the range of 100-130 K. However, Lobban, Finney, and Kuhs failed to observe a transition in later work.29 Johari and Whalley report that discrepancies in the dielectric constant hint at the partial antiferroelectric ordering in ice V at low temperatures.30 The low-temperature behavior of ice V will be treated in a subsequent work.24 In this work we report periodic density functional theory calculations for many hydrogen bond isomers of ice VI, and the prediction of the low-temperature phase transition that results from those calculations. Our calculations predict that hydrogen bond disordered ice VI should give way to H-bond ordered ice VI′ near 108 K. At this point, ice VI′ should be more stable than ice VI, and at least be a metastable, if not the globally stable, phase. Factors that might affect global stability are

10.1021/jp0540609 CCC: $30.25 © 2005 American Chemical Society Published on Web 10/19/2005

Hydrogen Bond Ordered Form of Ice VI

J. Phys. Chem. B, Vol. 109, No. 44, 2005 21041 partition function for ice can be written as a sum over M symmetry-distinct H-bond topologies M

Q)

Figure 1. Examples of different H-bond topologies possible for a 10water unit cell of ice VI. After the cell is periodically replicated, each water molecule is hydrogen bonded to four others, and the oxygen atoms are in nearly the same positions. The black bonds are covalent, and the white bonds are hydrogen bonds.

discussed in the final section. In section 1, we briefly review the graph invariant theory by which electronic structure calculations on small systems can be extrapolated to predict statistical behavior in ice. The periodic density functional theory31,32 results are presented in section 2, and our predictions for the thermodynamic limit are developed there. We discuss the implications of our results in section 3. 1. Graph Invariant Theory for Ice VI The hydrogen bond order/disorder phase transitions of ice share certain features. First, the oxygen positions are only slightly changed at the phase transition. For example, the a and c lattice parameters of ice VIII at 263 K differ by only -1.0% and +2.0% from their values in the perfect cubic structure of ice VII.29 Upon transformation of ice Ih to ice XI, a and b change by -0.75% and +0.84%, and c changes by -0.36%.8 While these changes are highly significant in certain respects, i.e., crystal strain upon transformation, they have a very slight effect on our calculated energies. This is not surprising because, to the extent that the experimental lattice parameters used in our calculations are the actual minimum energy configuration of our electronic structure model, the energy is variational with respect to the cell parameters. Therefore, we neglect the very small lattice parameter change upon hydrogen bond ordering. Second, isotope effects seem to have only minor effect on the transition temperature. For example, the ice Ih/XI transition in D2O occurs at 76 K, shifted upward by only 4 K from the transition temperature in H2O.2-4 Mass does not affect phase transition temperatures in classical statistical mechanics and is a signature of quantum effects. Evidently, quantum effects do not significantly affect the relative free energies of the hydrogen bond ordered and disordered phases, even for a relatively low temperature transition. Accordingly, we use classical statistical mechanics in our calculations. If we stay away from extremely high pressure where the hydrogen bond becomes symmetrical, the potential energy surface for ice Ih exhibits a number of deep minima, each corresponding to a different hydrogen bond topology. Examples of a few local minima corresponding to different H-bond topologies are shown in Figure 1. Working within the framework of classical statistical mechanics, the partition function can be written as a sum of contributions from each of M symmetry-distinct local minima of the potential energy surface.33-41 We have previously discussed20,21 how the

fie-β(E +A ∑ i)1 i

vib,i)

(1)

where fi is the number of symmetry-related configurations which are represented by one symmetry-distinct configuration, Ei is the potential energy at the ith potential minimum, and Avib,i is the vibrational free energy associated with movement near that minimum. At sufficiently low temperature the classical procedure could be modified to incorporate some quantum effects, for example, by calculating Avib,i quantum mechanically. However, the use of classical statistical mechanics seems warranted for H-bond order/disorder transitions in ice. The graph invariants we have introduced20-22 provide a feasible way to circumvent the need to calculate all the Ei and Avib,i values for the billions of hydrogen bond topologies found in a “simulation cell”, a unit cell large enough to approximate the thermodynamic limit. Each Ei and Avib,i in eq 1 corresponds to a local minimum of the potential surface defined by an H-bond topology, like the four examples shown in Figure 1. Graph invariants link quantities such as Ei and Avib,i to H-bond topology, and allow us to construct statistical models with a kind of bootstrap procedure.20,21 Each hydrogen bond topology is specified by assigning a variable br to the rth bond, which can take the values +1 or -1 according to whether the H-bond points in the same or opposite to an arbitrarily chosen reference direction for that bond. Each H-bond topology coincides with a symmetry distinct graph in which the directionality of H-bonds is indicated by directed bonds. If scalar physical quantities such as Ei and Avib,i depend on the H-bond topology as summarized by a directed graph, they must depend on combinations of the bond variable that are invariant to symmetry operationsshence the term graph inVariant. (Whether H-bond topology is indeed sufficient to describe properties of ice will be confirmed by analysis of the data in section 2.) The construction of graph invariants is accomplished by application of the group theoretical projection operator for the totally symmetric representation on bonds br, or products of bonds brbs, brbsbt, .... Projecting onto bonds gives an invariant linear combination of bonds, which we call a first-order invariant. For systems with sufficient symmetry, there are no first-order invariants.42 Next, projection onto bond pairs brbs produces quadratic polynomials in the bond variables which we call second-order invariants. When the first-order invariants vanish, second-order invariants provide the leading dependence of scalar physical properties on H-bond topology. To date, in our examinations of water clusters20,43 and ice21,22 we have found that sufficient accuracy for our needs was obtained by truncating the expansion with second-order invariants. In other words, third and higher order polynomials in the bond variables (third, fourth, ... order invariants) were not required. In our work on ice to date on the ice VII/VIII, Ih-XI, and III/IX transitions,22,24 we have obtained good agreement with available experiments assuming that Avib,i, the vibrational free energy in each of the H-bond isomers, is approximately the same: Avib,i ≈ A h vib. Then the problem reduces to estimating the energy Ei of each of the symmetry distinct isomers. We determine the Ei values by expressing them as a sum of secondorder graph invariants.

E(b1, b2, ...) ) E0 +

RrsIrs(b1, b2, ...) ∑ rs

(2)

21042 J. Phys. Chem. B, Vol. 109, No. 44, 2005

Figure 2. (a) A fragment of ice VI, as viewed perpendicular to the c-axis, depicting the 10-water primitive unit cell of ice VI, and some additional waters that help clarify the bonding pattern. Labels for the bonds used to generate the graph invariants, as described in Table 1, are indicated in the figure. The dashed lines represent H-bonds emanating from waters not explicitly shown in the figure. (b) A larger unit cell, measuring x2 × x2 × 2 primitive cells on each side, containing 40 water molecules as viewed down the c-axis. The lowest energy isomer, that of the ice VI′ phase, is shown. It consists of two interpenetrating, but not interconnected lattices. Bonds which appear horizontal in the figure all point to the left, while bonds that are vertical all point downward.

In the above equation, E0 is an overall constant for the energy, and the sum is over second-order invariant functions of the bond variables, Irs(b1, b2,...). We identify the second-order invariants with a double subscript “rs” to indicate a bond pair brbs that generates the invariant by action of the projection operator for the totally symmetry representation. Equation 2 assumes that the first-order invariants vanish either because of symmetry or the constraints of the ice rules for most of the common phases of ice. (Notable exceptions are ice III and V, for which the term ∑r RrIr(b1, b2,...) would appear in eq 2.) The coefficients Rrs are determined by fitting to the results of periodic density functional theory calculations. If the approximation Avib,i ≈ A h vib would break down, then an alternative would be to fit the Avib,i to invariants with an expression similar to the energy expansion in eq 2. The details of constructing the Irs(b1, b2, ...) given the crystalline space group are provided in another publication.21 The smallest unit cell we examined was the primitive unit cell containing 10 water molecules (Figure 2a). The projection operator for the totally symmetric representation, when applied to some of the bonds in Figure 2a, G ˆ (br), did not always yield zero, indicating that, algebraically, some first-order invariants exist. However, all the first-order invariants vanished when evaluated for bond configurations allowed by periodicity and the ice rules. Applying the projection operator to bond pairs, G ˆ (brbs), we find 22 linearly independent second-order invariants. Some of them, as is obvious for the three generated from projection on the square of one bond, G ˆ (brbr), evaluate to constants. Their effect is included by the constant E0 in eq 2. The remaining second order invariants could be sorted into groups based on geometrical features of the generating bond pairs, such as whether both of the generating bond pair lay in the same sublattice, and the distance between the two bonds of the generating pair. We found that only invariants generated by bonds in close proximity had significant weight in the energy fit, and were the only invariants necessary to achieve acceptable accuracy in fitting energy as a function of H-bond topology, as in eq 2. The invariants which were significant are given in Table 1, where the generating bond pairs are defined in reference to Figure 2a. One main result from our previous work21 is that invariants for a small unit cell are also invariants for a larger

Knight and Singer unit cell. A larger unit cell will, in general, also contain additional invariants which are generated by projection on bond pairs brbs which are separated further than possible in the smaller cell. At some point, we expect that addition of further invariants involving bond pairs far from each other will make no improvement in the energy eq 2. At that point, the invariant expansion is converged and electronic structure calculations on larger cells are no longer necessary. In practice, we have found for ice VII, VIII, Ih, XI,22 III, and IX24 that calculations on relatively small unit cells are sufficient. Furthermore, even within a small unit cell we only require bond pairs that either share a common vertex or for whom br and bs contain vertexes that are nearest neighbors. This is confirmed below for ice VI. Therefore, the expansion for which eq 2 is the leading term appears to converge quite rapidly. Once the coefficients Rrs have been determined from electronic structure calculations, the energy expression in eq 2 serves as the Hamiltonian for what is effectively a spin-lattice model. While it is possible to derive analytic approximations to solve this model, we have simply used relatively inexpensive Metropolis Monte Carlo simulations to obtain predictions without introducing further approximations. Extrapolation of electronic structure calculations to large cells with graph invariants is essential. For example, Kuo and Klein44 predicted that the ice VII/VIII transition should occur at 150 K based on energetics of a cell consisting of 16 water molecules. With use of similar energetics as input to the graph invariant theory, a transition temperature of 228 K was calculated, in better agreement with experimental reports of 274 K15,45 and 263 K.16 2. Results Ice VI consists of two interpenetrating, but not interconnected lattices with a tetragonal unit cell. The space group is P42/nmc. The lattice constants used in the following calculations were a ) 6.181 Å and c ) 5.698 Å as determined by diffraction experiments16 in the region of stability of ice VI, at 1.1 GPa and 225 K. Kuhs et al. report that the a- and c-axis lattice constants change by -0.15 Å (-2.4%) and 0.10 Å (+1.8%), respectively, upon decreasing temperature from 225 to 8 K at 1.1 GPa.16 The smallest unit cell we examined was the primitive unit cell containing 10 water molecules. The 45 symmetrydistinct H-bond isomers possible in this unit cell were enumerated, using previously described methods.20,44,46 Geometry optimizations were performed on all 45 H-bond isomers, using the Car-Parrinello31,47 method with periodic boundary conditions. We employed the Becke-Lee-Yang-Parr (BLYP) gradient correction48,49 to the local density approximation and Martins-Troullier norm-conserving pseudopotentials50 in the electronic density functional theory (DFT) calculations. The electronic wave function was expanded in plane waves up to a cutoff of 70 Ry. Sampling of the Brillouin zone was either restricted to the Γ point or extended by way of a 2 × 2 × 2 k-point grid, using the Monkhorst-Pack scheme.51 The ability of the expansion of eq 2, using invariants 1-6 of Table 1, to fit the energies of the isomers of the 10-water unit cell is evaluated in Figure 3a. The fitted value is plotted as a function of the actual DFT energy. If a perfect fit was achieved, all the points would lie on the diagonal. As can be seen from the figure, even though the energy difference between isomers is quite small, the typical error fitting Γ-point energies of the 45 symmetry-distinct isomers is a small fraction of the energy range. Furthermore, we recalculated some of the energies using 2 × 2 × 2 k-point sampling and, among those isomers, the fit was even better.

Hydrogen Bond Ordered Form of Ice VI

J. Phys. Chem. B, Vol. 109, No. 44, 2005 21043

TABLE 1: Geometrical Features and Contribution to the Description of the Energy of H-Bond Isomers of the Second Order Graph Invariantsa invariant

generating bond pair

distance (Å)

same sublattice?

10-water (Γ only)

10-water (2 × 2 × 2)

40-water (Γ only)

1 2 3 4 5 6 7

5,37 2,17 10,21 22,38 9,20 9,21 10,12

3.56 3.29 3.29 2.77 0.00 3.29 5.70

false false false true true false true

0.048 475 0 0.012 861 7 0.050 216 6 -0.010 345 5 -0.008 829 4 -0.032 844 6

0.046 790 6 0.022 426 7 0.055 241 3 -0.010 527 0 -0.045 172 5 -0.039 177 3

0.046 984 0 0.016 148 0 0.045 166 4 -0.014 333 2 -0.036 714 5 -0.038 350 6 -0.004 970 8

a Invariant 7 does not appear as a linearly independent invariant in the 10-water, primitive 1 × 1 × 1 cell. The indices of the generating bond pair refer to Figure 2. The distance associated with each bond pair is the distance between the closest vertices from each bond in an ideal structure before geometry optimization. The last three columns give the fitting coefficients for each of the invariants as used in eq 2 for the energy in units of kcal mol-1 water-1.

Figure 3. (a) Test of the ability of the graph invariant expansion, eq 2, to fit the energies of the 45 H-bond isomers of a 10-water unit cell of ice VI. The filled symbols are the energies calculated with only the Γ point, and the open symbols are a fit to the energies of 15 isomers calculated with 2 × 2 × 2 k-point sampling. (b) Test of the ability of graph invariant parameters derived from the 10-water cell to predict energies of a larger 40-water cell. The filled symbols are a comparison for parameters derived from Γ point calculations for the 10-water cell. The open symbols are a similar prediction with extended 2 × 2 × 2 k-point sampling. In both plots, a line of slope unity is shown to indicate where points would lie for perfect agreement.

As mentioned in section 1, any invariant present in a small unit cell is also an invariant for a larger unit cell whose cell vectors are multiples of those of the smaller cell. Figure 3b shows that by the time we move beyond the 10-water cell, we are very close to the point where the additional invariants that arise from larger cells are not important in describing the energy of the H-bond isomers. We performed DFT for 54 isomers of a 40-water cell measuring x2 × x2 × 2 primitive cells on each side. The 40-water cell energies predicted by using invariant functions from the smaller 10-water cell with the coefficients determined from the smaller cell calculation are compared with the calculated energies for the 40-water cell in Figure 3b. Using the coefficients from the Γ-point calculation for the small cell, we find that the prediction from the small cell has the correct trend, but the prediction underestimates the range of energies in the larger cell. Actually, this effect does not indicate lack of convergence of the invariant expansion, but rather the inadequacy of taking only the Γ-point for the small cell. Coefficients determined by 2 × 2 × 2 k-point sampling of the small cell, the open symbols in Figure 3, yield an excellent prediction of the energies in the 40-water cell. The larger cell, even though calculated at the Γ-point, effectively includes more k-point sampling than a Γ-point calculation for the small cell. Therefore, the small-cell coefficients obtained with extended k-point sampling do a better job predicting the energetics of the large cell. Ideally, one needs to converge the energies with respect to k-point sampling for each cell size. Among the new invariants possible for the larger 40-water cell that are not present in the 10-water cell, all the first order

Figure 4. The smallest repeating unit for each of the independent lattices that generate (a) the lowest energy and (b) second lowest energy isomer, as determined from DFT calculations on a 40-water unit cell of ice VI. The configurations are viewed perpendicular to the c-axis. The H-bonds for each lattice in part b are antiferroelectrically oriented, making the overall structure antiferroelectric. This arrangement of H-bonds has tetragonal symmetry and is assigned the space group P212121 as determined with the FINDSYM52 program. The H-bonds parrallel to the a- and b-axes point in a counterclockwise fashion as one looks down the c-axis. The bonds that must be reversed to interconvert the two structures are circled on the left. Spatial coordinates for these structures are available upon request from the authors.

invariants evaluated to zero for configurations allowed by periodicity and the ice rules, just as they did for the smaller cell. We included one of the new second order invariants, invariant 7 in Table 1, for the following reason. In the large unit cell, the two lowest energy H-bond isomers, shown in Figure 4, are nearly degenerate. Invariant functions 1-6 in Table 1 have exactly the same value for these two isomers. We therefore included a new second-order invariant from the 40water cell that broke the degeneracy between these isomers. As can be seen from the last column in Table 1, the coefficient of this invariant is quite small. The lowest energy structure of the 10-water unit cell was a ferroelectric structure (Figures 2 and 4a) whose space group is Cc, as confirmed by using the FINDSYM program.52 We attempted to identify alternative low-energy structures for the larger 40-water cell using Monte Carlo sampling based on invariants obtained from the 10-water cell, much like the Monte Carlo simulations described below. Of the lowest energy structures identified in this way, the ferroelectric Cc structure was still the ground state. The energy of an antiferroelectric structure of P212121 symmetry (right side of Figure 4) was only 4 K per water molecules higher in energy, as determined from DFT calculations on the 40-water unit cell. The two structures actually have many H-bonds oriented in the same way. Those bonds that must be reversed to go from the ferrroelectric to

21044 J. Phys. Chem. B, Vol. 109, No. 44, 2005

Knight and Singer

Figure 5. (a) Average energy plotted as a function of temperature from Metropolis Monte Carlo simulations for a simulation cell of ice VI containing 1250 waters. (b) Ferroelectric ordering of the ice VI lattices as a function of temperature. Data are presented for a series of Metropolis Monte Carlo runs ascending (4) and descending (3) in temperature.

antiferroelectric structures are circled in Figure 4. While the electronic structure methodology we use has proved to be remarkably accurate for the ice VII/VIII, Ih-XI, and III/IX transitions,22,24 the very small energy difference between the two structures in Figure 4 implies that both structures should be considered as candidates for the ground state of ice VI′. We explored the effect of reversing the energetic ordering of the two structures shown in Figure 4 by artifically changing the sign of the cofficient of invariant 7 in Table 1, and report those results below. Using the invariant parameters obtained from the 40-water unit cell and eq 2 for the energy, we have an effective Hamiltonian describing H-bond fluctuations. In essence, we have turned the H-bond order/disorder problem into a spin-lattice model. Rather than introduce further approximations, we choose an essentially exact numerical solution using the standard Metropolis Monte Carlo algorithm (e.g., ref 53) for which the computational cost is relatively modest. The only nonstandard feature is the algorithm for generating trial moves. We must choose trial rearrangements of H-bonds that do not violate the ice rules. We use the algorithm invented by Rahman and Stillinger54 to find closed loops of H-bonds that are reversed to generate configurations in ice that obey the ice rules. Trial moves that preserve the ice rules are generated by chains of H-bonds that either begin and end on the same oxygen atom or else, in a simulation with periodic boundary conditions, end on a periodic replica of the starting point. Rahman and Stillinger only accepted loops where the starting and end points lay on the same oxygen atom. In an approximation where the total dipole is a sum of bond dipoles it is easily seen that closed loops do not change the total dipole moment of the system. However, as Rahman and Stillinger noted, loops that begin in one periodic simulation cell and terminate in another cell will change the dipole moment. We allow both types of moves because we want to sample all H-bond configurations, including both ferroelectric

and antiferroelectric configurations, and allow exact statistical simulations to identify the equilibrium properties of the system. Metropolis Monte Carlo simulations were performed on a simulation cell measuring five primitive unit cells on each side containing 1250 waters. The seven invariant coefficients obtained from the 40-water unit cell were used to evaluate eq 2 for the simulation cell. A series of simulations were performed for both decreasing and increasing temperature. The lowest energy isomer, as determined from electronic DFT calculations, was the starting configuration for the initial low-temperature simulation. The final H-bond configuration in a preceding simulation was always the initial configuration in the next simulation. A highly disordered configuration obtained from simulation at an extremely large temperature (∼107 K) was used to initialize the sequence of simulations descending in temperature. The simulations predict a first-order phase transition near 108 K to the ferroelectric ground-state identified from calculations on the 10- and 40-water unit cells. As shown in Figure 5, negligible hysteresis, 1-2 degrees, is observed from the series of simulations with decreasing and increasing temperature. The predicted proton ordering transition temperature is located in a region of the phase diagram where the ice II-VI and VI-VIII extrapolated phase boundaries have not intersected. Partial H-bond ordering above the transition and disordering below the transition is observed. Shown in Figure 5b is 〈cos θab〉, the average cosine between unit vectors in the direction of the total dipole for each sublattice of the ice VI structure.

〈cos θab〉 )

〈

Ma·Mb |Ma||Mb|

〉

(3)

In the above equation, Ma and Mb are the total dipole moments of the two sublattices of the ice VI lattice calculated by using a bond dipole approximation. Ferroelectric ordering of the H-bonds takes place over a wide temperature range gradually starting from above 300 K.

Hydrogen Bond Ordered Form of Ice VI

J. Phys. Chem. B, Vol. 109, No. 44, 2005 21045

Figure 6. Entropy plotted as a function of temperature calculated from Metropolis Monte Carlo simulations on a 1250-water simulation cell. The horizontal line is the Pauling entropy for a fully disordered ice lattice. With decreasing temperature, 29.8% is lost before the transition, 45.5% at the transition, and 24.7% as the ferroelectric ground state is formed.

Even though the low-energy ferroelectric and antiferroelectric configurations shown in Figure 4 are quite close in energy, the simulations rapidly convert to the ferroelectric configuration at all temperatures below the transition temperature. We found no evidence of a barrier to this conversion. Given that the lowest energy ferroelectric and next-lowest energy antiferroelectric structures are quite close in energy, we sought to determine the effect of reversing the energetic ordering of these two structures. We could accomplish this conveniently by changing the sign of the small coefficient of invariant 7 in Table 1, the one added to distinguish the energy of these two structures. After this change, the energies of the other configurations are barely shifted. We find that reversing the energetic order of the two low-energy structures causes ice VI′ to be antiferroelectric, but the transition occurs within 2 K of the transition to a ferroelectric phase found with the original parameters. Using either the original value of the coefficient for invariant 7 or the signreversed value, we never found competition between ferroelectric and antiferroelectric order that would produce more evidence of partial disorder below the transition. That the system might encounter two low-temperature phases upon cooling, one antiferroelectric and the other ferroelectric, is conceivable. It is also possible that formation of one phase as a metastable phase might be kinetically favored. Neither of these possibilities occurred in our simulations. Entropy as a function of temperature is shown in Figure 6. The entropy of the low-temperature phase was calculated via thermodynamic integration from 0 K and that of the hightemperature phase from infinite temperature. The entropy at infinite temperature is taken to be the Pauling entropy for a fully disordered arrangement of H-bonds in ice. With decreasing temperature, 29.8% of entropy is lost before the transition while 24.7% is lost after the transition. The calculated entropy at the transition is 45.5% of the ideal value for a fully disordered ice phase. 3. Discussion In this work, we have presented the results of electronic DFT calculations on H-bond isomers for two unit cells of ice VI. The lowest energy H-bond isomer in each case was ferroelectric, space group Cc. From calculations on a larger unit cell, an antiferroelectric H-bond isomer was identified lying 4 K per molecule higher in energy. This configuration was assigned the space group P212121. From application of graph invariants and electronic DFT calculations on small unit cells, the energy differences in H-bond fluctuations for a large simulation cell, containing 1250 waters, could be calculated. A first-order phase transition to the proton ordered ferroelectric ground state (Figure 4a) was observed near 108 K. We note that an antiferroelectric structure (Figure 4b)

lies very close to the ferroelectric structure in energy. If more accurate electronic structure calculations would indicate that the antiferroelectric structure is actually lower in energy, our explorations with the antiferroelectric structure as the ground state reported in section 2 suggest that we would find a transition at almost the same temperature, but to an antiferroelectrically ordered phase. We have shown that hydrogen bond disordered ice VI gives way to a proton-ordered phase near 108 K. While we have shown that the ordered phase is more stable than ice VI, establishing this phase as the global minimum in this portion of the phase diagram is another, and more difficult, matter. To establish whether ice VI′ is globally stable we would have to determine whether there exists another phase with a different underlying oxygen lattice that supersedes ice VI′ as the stable phase in this region. In fact, this scenario does occur in the case of ice III, which can order into metastable ice IX which is superseded by the more stable ice II. There is also a possibility that ice VI might give way to either ice II or ice VIII at low temperatures. In the phase diagram of ice, the range of pressure where ice VI is stable seems to shrink as temperature is lowered, as the ice II and ice VIII regions expand. If the ice II and ice VIII regions would expand and meet at low temperature, providing a lower bound to the ice VI phase, then our ice VI′ phase would be metastable with respect to either ice II or ice VIII. At our calculated transition temperature of 108 K, ice II and ice VIII are typically drawn as separate,23 but the phase diagram is not well characterized in this region. In summary, we cannot exclude another phase intervening and rendering ice VI′ only metastable. However, even in such a case (like ice III/IX) experimental evidence for the metastable phase can be obtained. In our calculations, significant proton ordering was observed above the transition over a wide temperature range. This should be observable in calorimetric experiments provided that hydrogen bond arrangements can equilibrate on an experimental time scale. (The H-bonds in phases such as ices VII and VIII seem to reach equilibrium on an experimental time scale, while those in ices Ih and XI do not.) Even below the transition, the degree of residual proton disorder is much larger than we found in ice XI or ice VIII. We were somewhat surprised that, despite the near degeneracy of ferroelectric and antiferroelectric structures, fluctuations to antiferroelectric configurations do not persist to lower temperatures. Simulations initialized with the H-bond configuration of the second lowest energy antiferroelectric isomer rapidly transformed to the ferroelectric ground state, indicating that factors beside energetics (i.e., entropic factors) seem to also favor the ferroelectric state. In Figure 5b, we see that the dipole moment of the sublattices has nearly complete ferroelectric order at all temperatures below the transition.

21046 J. Phys. Chem. B, Vol. 109, No. 44, 2005 In conclusion, we have predicted that ice VI will transform into a hydrogen bond ordered phase near 108 K, and have identified the likely low-temperature phase as ferroelectric with an antiferroelectric structure close by in energy. Hopefully our results will motivate further characterization of low-temperature ice. Acknowledgment. We gratefully acknowledge NSF Grant No. CHE-0109243 for support of this work, and the Ohio Supercomputer Center for resources needed to perform the calculations reported here. References and Notes (1) Kawada, S. J. Phys. Soc. Jpn. 1972, 32 (5), 1442. (2) Tajima, Y.; Matsuo, T.; Suga, H. Nature 1982, 299, 810. (3) Tajima, Y.; Matsuo, T.; Suga, H. J. Phys. Chem. Solids 1984, 45 (11-12), 1135. (4) Matsuo, T.; Tajima, Y.; Suga, H. J. Phys. Chem. Solids 1986, 47 (2), 165. (5) Matsuo, T.; Suga, H. J. Phys., Colloq. 1987, C1, 477. (6) Leadbetter, A. J.; Ward, R. C.; Clark, J. W.; Tucker, P. A.; Matsuo, T.; Suga, H. J. Chem. Phys. 1985, 82 (1), 424. (7) Howe, R.; Whitworth, R. W. J. Chem. Phys. 1989, 90 (8), 4450. (8) Line, C. M. B.; Whitworth, R. W. J. Chem. Phys. 1996, 104 (24), 10008. (9) Jackson, S. M.; Nield, V. M.; Whitworth, R. W.; Oguro, M.; Wilson, C. C. J. Phys. Chem. 1997, B101 (32), 6142. (10) Kamb, B.; Prakash, A. Acta Crystallogr. 1968, B24, 1317. (11) Whalley, E.; Heath, J. B. R.; Davidson, D. W. J. Chem. Phys. 1968, 48 (5), 2362. (12) LaPlaca, S. J.; Hamilton, W. C.; Kamb, B.; Prakash, A. J. Chem. Phys. 1973, 58 (2), 567. (13) Nishibata, K.; Whalley, E. J. Chem. Phys. 1974, 60 (8), 3189. (14) Londono, J. D.; Kuhs, W. F.; Finney, J. L. J. Chem. Phys. 1993, 98 (6), 4878. (15) Johari, G. P.; Lavergne, A.; Whalley, E. J. Chem. Phys. 1974, 61 (10), 4292. (16) Kuhs, W. F.; Finney, J. L.; Vettier, C.; Bliss, D. V. J. Chem. Phys. 1984, 81 (8), 3612. (17) Kamb, B. Science 1965, 150, 205. (18) Kamb, B. In Physics and chemistry of ice; Symposium on the Physics and Chemistry of Ice; Whalley, E., Jones, S. J., Gold, L. W., Eds.; Royal Society of Canada: Ottawa, Canada, 1973; p 28. (19) Actually, Kamb uses VI′ to indicate a partially ordered phase, and VI′′ to denote the fully ordered form of ice VI. (20) Kuo, J.-L.; Coe, J. V.; Singer, S. J.; Band, Y. B.; Ojama¨e, L. J. Chem. Phys. 2001, 114 (6), 2527.

Knight and Singer (21) Kuo, J.-L.; Singer, S. J. Phys. ReV. 2003, E67 (1), 016114. (22) Singer, S. J.; Kuo, J.-L.; Hirsch, T. K.; Knight, C.; Ojama¨e, L. P.; Klein, M. L. Phys. ReV. Lett. 2005, 94 (13), 135701. (23) Petrenko, V. F.; Whitworth, R. W. Physics of Ice; Oxford University Press: New York, 1999. (24) Knight, C.; Singer, S. J., miscellaneous work. (25) Johari, G. P.; Whalley, E. J. Chem. Phys. 1976, 64 (11), 4484. (26) Johari, G. P.; Whalley, E. J. Chem. Phys. 1979, 70 (5), 2094. (27) Kamb, B.; Placa, S. J. L. Trans. Am. Geophys. Union 1974, 56 (V39), 1202. (28) Handa, Y. P.; Klug, D. D.; Whalley, E. J. Phys. (Paris) 1987, 48, C1. (29) Lobban, C.; Finney, J. L.; Kuhs, W. F. J. Chem. Phys. 2000, 112 (16), 7169. (30) Johari, G. P.; Whalley, E. J. Chem. Phys. 2001, 115 (7), 3274. (31) Car, R.; Parrinello, M. Phys. ReV. Lett. 1985, 55 (22), 2471. (32) CPMD, Copyright IBM Corp. 1990-2001, Copyright MPI fu¨r Festko¨rperforschung Stuttgart 1997-2004. (33) McGinty, D. J. J. Chem. Phys. 1971, 55 (2), 580. (34) Burton, J. J. J. Chem. Phys. 1972, 56 (6), 3133. (35) Hoare, M. R. AdV. Chem. Phys. 1979, 40, 49. (36) Stillinger, F. H.; Weber, T. A. Phys. ReV. 1982, A25 (2), 978. (37) Stillinger, F. H.; Weber, T. A. Phys. ReV. 1983, 28 (4), 2408. (38) Mezey, P. G. Studies in Physical and Theoretical Chemistry; Potential Energy Hypersurfaces, Vol. 53; Elsevier: New York, 1987. (39) Franke, G.; Hilf, E. R.; Borrmann, P. J. Chem. Phys. 1993, 98 (4), 3496. (40) Wales, D. J. Mol. Phys. 1993, 78 (1), 151. (41) Doye, J. P. K.; Wales, D. J. J. Chem. Phys. 1995, 102 (24), 9659. (42) The conditions necessary and sufficient for vanishing of all firstorder invariants were presented in reference KuoCoeSingerBandOjamae2001. (43) Belair, S. D.; Francisco, J. S.; Singer, S. J. Phys. ReV. 2005, A71 (1), 013204. (44) Kuo, J.-L.; Klein, M. L. J. Phys. Chem. 2004, B108 (51), 19634. (45) Johari, G. P.; Lavergne, A.; Whalley, E. J. Chem. Phys. 1980, 73 (8), 4150. (46) McDonald, S.; Ojama¨e, L.; Singer, S. J. J. Phys. Chem. 1998, A102 (17), 2824. (47) Car, R.; Parrinello, M. In Simple Molecular Systems Very High Density; Polian, A., Loubeyre, P., Boccara, N., Eds.; Plenum: New York, 1989; Vol. 186, NATO ASI Ser. B, p 455. (48) Becke, A. D. Phys. ReV. 1988, A38 (6), 3098. (49) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. 1988, B37, 785. (50) Troullier, N.; Martins, J. L. Phys. ReV. 1991, 43 (3), 1993. (51) Monkhorst, H. J.; Pack, J. D. Phys. ReV. B 1976, 13 (12), 5188. (52) Stokes, H. T.; Hatch, D. M. 2004, FINDSYM, stokes.byu.edu/ isotropy.html. (53) Binder, K., Ed. Monte Carlo Methods in Statistical Physics, 2nd ed.; Topics in Current Physics, Vol. 7; Springer: New York, 1986. (54) Rahman, A.; Stillinger, F. H. J. Chem. Phys. 1972, 57 (9), 4009.