Hydrogen bond topology and the ice VII/VIII and Ih/XI proton ordering phase transitions

Ice Ih, ordinary ice at atmospheric pressure, is a proton-disordered crystal that when cooled under special conditions is believed to transform to ferroelectric proton-ordered ice XI, but this transformation is still subject to controversy. Ice VII, also proton disordered throughout its region of stability, transforms to proton-ordered ice VIII upon cooling. In contrast to the ice Ih/XI transition, the VII/VIII transition and the crystal structure of ice VIII are well characterized. In order to shed some light on the ice Ih proton ordering transition, we present the results of periodic electronic density functional theory calculations and statistical simulations. We are able to describe the small energy differences among the innumerable H-bond configurations possible in a large simulation cell by using an analytic theory to extrapolate from electronic DFT calculations on small unit cells to cells large enough to approximate the thermodynamic limit. We first validate our methods by comparing our predictions to the well-characterized ice VII/VIII proton ordering transition, finding agreement with respect to both the transition temperature and structure of the low-temperature phase. For ice Ih, our results indicate that a proton-ordered phase is attainable at low temperatures, the structure of which is in agreement with the experimentally proposed ferroelectric $Cmc{2}_1$ structure. The predicted transition temperature of $98\mathrm{K}$ is in qualitative agreement with the observed transition at $72\mathrm{K}$ on KOH-doped ice samples.

Figure 1

(Color online) Four possible arrangements of H bonds within a 16-water-molecule $2\times 1\times 1$ orthorhombic unit cell of ice Ih. Here, cis and trans bonds are defined as whether protons lie on the same or opposite side of the H bond respectively, as indicated for isomer (a). The H-bond isomers are summarized mathematically by directed graphs in which directional bonds point from H-bond donor to H-bond acceptor, as illustrated for the isomer (b).

Figure 2

(Color online) (a) A $1\times 1\times 1$ hexagonal unit cell of ice Ih, space group $P{6}_3∕mmc$, containing 12 water molecules. (b) A $1\times 1\times 1$ orthorhombic unit cell of ice Ih, space group $Cmc{2}_1$, containing 8 water molecules. Bonds that appear to point straight up and down are parallel to the $c$ axis and hence called $c$-axis bonds. The remaining bonds, perpendicular to the $c$ axis, are referred to as $ab$ bonds. (c) Proposed experimental structure of proton ordered ice Ih, ice XI, as determined from diffraction and thermal depolarization experiments [13, 14, 15, 16, 17, 18]. The $ab$-layer bonds (1) are all cis with the $ab$ layers polarized parallel to the $b$ axis and alternating layers oppositely aligned. The $c$-axis bonds (2) are trans and all oriented in the same direction. cis and trans H bonds are defined, respectively, as to whether the non-hydrogen-bonded hydrogen atoms fall on the same or opposite side of the H bond.

Figure 3

(Color online) (a) An isomer of a two-water-molecule primitive unit cell of ice VII, obeying the Bernal-Fowler ice rules, is shown. The thin black lines outline the unit cell and neighboring oxygen atoms are included for clarity. The orientations of H bonds in this isomer are assigned to be the canonical arrangement of H bonds. The H bonds are labeled from 1 to 4 to provide a means to associate a bond variable, $b_r$, with each H-bond labeled $r$. (b) The H-bond configuration in (a) is summarized by a directed graph. The H bonds are taken to point from oxygen donor to oxygen acceptor as discussed in the text. (c) and (d) Additional directed graphs corresponding to other H-bond isomers of ice VII that satisfy the Bernal-Fowler ice rules and periodicity constraints. As an example, if all the bond variables $b_r$ for configuration (b) were assigned the value $+1$, then all the $b_r$’s for configuration (d) would take the value $-1$ since all H bonds are reversed.

Figure 4

(Color online) (a) The canonical orientation of H bonds, obeying the Bernal-Fowler ice rules, for a unit cell of ice VII measuring $2\times 1\times 1$ primitive cells on each side. The black lines outline the primitive unit cells, and neighboring oxygen atoms are included for clarity. H bonds are labeled so as to identify a bond variable $b_r$ with H bond $r$. The orientations of H bonds in this isomer are assigned to be the canonical arrangement of H bonds, and all bond variables are assigned the value $+1$. (b)–(d) Isomers of three symmetry-distinct H-bond configurations possible in this unit cell. The bond variables for each of the three configurations, assigned according to H-bond configuration (a), are tabulated in Table .

Figure 5

(Color online) (a) An H-bond isomer of a 16-water-molecule unit cell of ice VII measuring two primitive unit cells on each side. Bonds representative of the three second-order graph invariants used to fit the DFT energies are shown, as further described in Table . (b)The ground-state H-bond isomer of a 32-water-molecule unit cell of ice VII measuring $2\sqrt{2}\times 2\sqrt{2}\times 2$ primitive cells on each side corresponding to the experimentally determined ice VIII structure. Bond pairs representative of the second-order graph invariants, including bond pairs not possible in the smaller 16-water-molecule unit cell, used to fit the DFT energies are shown. Bond 37 connects to a water molecule in an adjoining cell.

Figure 6

(a) Graph-invariant fit to the energies of the 52 H-bond isomers of a 16-water-molecule unit cell of ice VII. (b) Calculated DFT energy of H-bond isomers of a 32-water-molecule ice VII cell plotted against energies predicted from graph-invariant parameters derived from the 16-water-molecule cell. (c) Graph-invariant fit, using second-order invariants whose generating bond pairs are farther apart than possible in the smaller 16-water-molecule unit cell, to the energies of the H-bond configurations for the 32-water-molecule unit cell. (d) Same as plot (c) except only invariants whose generating bond pairs exist in the smaller 16-water-molecule unit cell were fit to the energies. A line of slope unity is shown to indicate where points would lie for perfect agreement. Invariant coefficients for each of the three fits are listed in Table .

Figure 7

Relative DFT energy of H-bond isomers of a 16-water-molecule (●) and 32-water-molecule (엯) unit cell of ice VII plotted against fraction of trans bonds for each isomer. The lowest- and highest-energy isomers for both unit cells contain no H bonds in the trans configuration, thus indicating that features of the H-bond topologies other than cis-trans H bonds are important if physical properties are to be correctly described.

Figure 8

(a) Average energy plotted as a function of temperature from Metropolis Monte Carlo simulations for large simulation cell of ice VII/VIII. Data are presented for series of Metropolis Monte Carlo runs ascending (▵) and descending (▿) in temperature. The vertical line is located at the calculated transition temperature near $228\mathrm{K}$. (b) Entropy plotted as a function of temperature. The horizontal line is the Pauling entropy for a fully disordered ice lattice subject to the ice rules. (c) Degree of antiferroelectric ordering of the ice VII sublattices as a function of temperature. ${\mathbf{M}}_a$ and ${\mathbf{M}}_b$ are the total dipoles, in units of bond dipoles, for the two independent sublattices.

Figure 9

Relative energy of H-bond isomers calculated by periodic electronic DFT methods for 16 isomers of an eight-water-molecule orthorhombic unit cell listed in order of increasing fraction of trans H bonds. The lowest graph (dotted lines) gives the fraction of trans H bonds associated with each isomer. The energy of the H-bond isomers was calculated using the programs (●, 엯) CPMD, (∎, ◻) DMOl, and (▴, ▵) CASTEP. Solid lines: energy of H-bond isomers before geometry optimization. Dashed lines: energies after optimization of the molecular coordinates and for the CASTEP results cell dimensions as well. The six energy data sets, optimized and unoptimized, are plotted with their average taken as the zero of energy to facilitate comparison of the relative energies of the isomers. For clarity, the $Cmc{2}_1$ isomer is noted and the optimized data sets are shifted by $0.06\mathrm{kcal}∕\mathrm{mol}$.

Figure 10

Relative energy of H-bond isomers calculated by periodic electronic DFT methods for 14 isomers of a 12-water-molecule hexagonal unit cell listed in order of increasing fraction of trans H bonds. The lowest graph (dotted lines) gives the fraction of trans H bonds associated with each isomer. The energy of the H-bond isomers was calculated using the programs (●, 엯) CPMD, (∎, ◻) DMOl, and (▵, ▴) CASTEP. Solid lines: energy of H-bond isomers before geometry optimization. Dashed lines: energies after optimization of the molecular coordinates and, for the CASTEP results, cell dimensions as well. The six energy data sets, optimized and unoptimized, are plotted with their average taken as the zero of energy to facilitate comparison of the relative energies of the isomers. For clarity, the $Cmc{2}_1$ isomer is noted and the optimized data sets are shifted by $0.06\mathrm{kcal}∕\mathrm{mol}$.

Figure 11

(Color online) (a) An H-bond isomer of a 12-water-molecule primitive unit cell of ice Ih. Bonds representative of the three second-order graph invariants used to fit the DFT energies are shown, as further described in Table . All bonds used to generate second-order invariants, used to describe energy differences for H-bond fluctuations in a large simulation cell, lie perpendicular to the $c$ axis and are referred to as $ab$ bonds. (b) An H-bond isomer of a 48-water-molecule unit cell of ice Ih measuring $2\times 2\times 1$ primitive cells on each side. Both H-bond isomers shown are the lowest-energy isomer for each unit cell in agreement with the experimentally proposed ferroelectric, space group $Cmc{2}_1$, ice XI structure. Arrows indicate the direction of the relative displacement $ϵ∕2$ of the $ab$ layers which are oppositely polarized.

Figure 12

(a) Graph-invariant fit to the energies of the 14 H-bond isomers of a 12-water-molecule hexagonal (●) unit cell and the 16 H-bond isomers of an eight-water-molecule orthorhombic (▴) unit cell of ice Ih. (b) Calculated DFT energy of H-bond isomers of a 48-water-molecule hexagonal ice Ih unit cell plotted against energies predicted from graph-invariant parameters derived from the small unit cells. (c) Graph-invariant fit to the energies of the 63 “semirandomly” chosen H-bond isomers of a 48-water-molecule hexagonal unit cell of ice Ih. A line of slope unity is shown to indicate where points would lie for perfect agreement.

Figure 13

(a) Average energy plotted as a function of temperature from Metropolis Monte Carlo simulations of a large simulation cell of ice Ih. Data are presented for series of Metropolis Monte Carlo runs ascending (▵) and descending (▿) in temperature. (b) Entropy plotted as a function of temperature. The horizontal line is the Pauling entropy $k_b\ln \frac{3}{2}$ for a fully disordered ice lattice.