Computational Analysis of Hydrogen Bond Vibrations of Ice III in the Far-Infrared Band

Abstract

The hydrogen-disordered structure of ice III makes it difficult to analyze its vibrational spectrum theoretically. To clarify the contribution of hydrogen bonds (HBs), we constructed a 24-molecule supercell to mimic the real structure and performed first-principles density functional theory calculations. The calculated curve of phonon density of states showed good correspondence with the experimental data. Based on the theory of two kinds of HB vibrational modes, we analyzed the distributions of two-bond modes and four-bond modes. The energy splitting of these modes results in a flat vibrational band, which is a common phenomenon in high-pressure ice phases. These findings verified the general rule that there are two types of HB vibrations in ice, thereby furthering our understanding of HB interactions in water ice and their broad role in nature.

1. Introduction

As a condensed state of water, ice is the most prevalent molecular crystal on Earth. By 2020, the year in which Jiang et al. [1] experimentally confirmed the existence of two-dimensional ice, over 20 amorphous and crystal ices had been revealed [2] in various environments. Most of these crystal phases can be classified into two categories according to whether their hydrogen atoms are long-range disordered or ordered in their lattices. There are six established pairs of disorder–order mated isomers: Ih/XI, III/IX, V/XIII, VI/XV, VII/VIII, and XII/XIV [3]. Each of these counterparts has identical oxygen sublattices but different hydrogen arrangements [4].

Ice III, one of the hydrogen-disordered (H-disordered) phases of ice, exists at 240–260 K and 0.2–0.35 GPa [5,6]. Because the oxygen locations are analogous to those of the silicon atoms in the SiO₂ phase keatite, ice III is also called keatite ice [7]. The stable region of ice III covers a relatively small range of the phase diagram [5]. Commonly, metastable ice III is prepared in an ice II phase environment by raising the compression rate in the density-driven phase transition from hexagonal ice at 170 K [8]. Tammann was the first to identify the existence of ice III in 1900, with his work later confirmed by Bridgman [9,10]. In 1968, Whalley et al. found a gradual disorder–order transition from ice III to an antiferroelectric phase, denoted as ice IX [11]. In 1973, based on measurements of the heat and entropy of the ice III–IX transformation, Whalley and Nishibata [12] inferred that ice III transforms into partially disordered ice IX within a wide range of temperatures. Furthermore, the space group P4₁2₁2 allows the hydrogen atoms to become partially ordered or disordered in ices III and IX, as proposed by Knight et al. [13]. In 2021, Salzmann reported a directly reversible phase transition between ices III and IX, which can be achieved by using ammonium fluoride as an anti-ice-II agent [14].

Since two distinct peaks in the far-infrared (IR) region of ice were first probed in the 1990s by inelastic neutron scattering (INS) experiments [15,16,17], the origin of these two vibrational peaks has remained a controversial issue [18,19,20]. Using first-principles density functional theory (DFT), we developed an approach to determine the lattice dynamic process of the crystal vibrational spectrum. Two kinds of intrinsic HB vibrational modes, denoted two-bond and four-bond modes, were found for the ideal ice Ic [21]. A series of investigations of the ice phases has revealed that this is a general rule for the ice family and may successfully explain their far-IR band [18,19,20]. In the present study, we theoretically analyzed the vibrational modes in the translation band to identify the characteristic HB vibrations of ice III in comparison with those of the H-ordered ice IX.

2. Computational Method

For the H-disordered ice III, there is no periodically repeating unit. Using the source code GenIce [22], we constructed a $1\times 2\times 1$ 24-molecule supercell of ice III, with hydrogen atoms randomly rearranged on the basis of Bernal–Fowler rules [23], to mimic the periodic primitive cell. From these stochastic results, we selected the cell with the lowest polarizability to use in our ideal model. This approach can improve the calculation efficiency, and it proved effective in our previous study on ice XII [24].

First-principles DFT calculations were conducted with the Cambridge Serial Total Energy Package (CASTEP) [25]. We used the generalized gradient approximation to handle the complicated spatial distribution of electron density. The revised Perdew–Burke–Ernzerhof (RPBE) exchange–correlation functional was selected because it gives rise to the most satisfactory accuracy for the vibrational spectrum calculation of ice [18]. The energy threshold and self-consistent field tolerance were both set at $1\times {10}^{-9}$ eV/atom for the geometry optimization. Norm-conserving pseudopotentials were introduced to calculate the phonon density of states (PDOS) with the linear response method. The energy cutoff was 830 eV, and the K-point mesh was $2\times 1\times 2$. The environmental pressure was set at 0.3 GPa to simulate real experimental conditions.

3. Results and Discussion

In general, thermal neutrons exchange energy with phonons throughout the first Brillouin zone in an INS experiment, and the calculated PDOS curve reflects the entire PDOS. Thus, the PDOS curve can be directly verified by the INS spectrum, as shown in Figure 1.

In the far-IR region, the INS spectrum of ice III shows an absorption band from 14.2 meV to 40.1 meV (114 cm⁻¹ to 327 cm⁻¹), which matches the band from 118 cm⁻¹ to 317 cm⁻¹ of the PDOS. Compared with the Raman scattering results for ice III [27,28], there is a peak at 188 cm⁻¹, which matches the peak at 198 cm⁻¹ in PDOS.

For a 12-molecule primitive cell of ice IX, there are $12\times 3\times 3-3=105$ optical normal modes. By comparison, a 24-molecule supercell of ice III contains $24\times 3\times 3-3=213$ optical normal modes. Because ice III is H-disordered, all the modes are nondegenerate, and many close-lying modes merge to produce a smoother curve in the spectrum than that of ice IX. The PDOS of ices III and IX are consistent with their respective INS spectra. In the translation band, it is clear that the two typical HB peaks are at approximately 200 cm⁻¹ and 300 cm⁻¹ for ice IX [20], whereas there is a flat vibrational band for ice III in the same region. To analyze the HBs of ice III, we first introduce the theory of two kinds of intrinsic HB vibrational modes in ice.

Subject to Bernal–Fowler rules, every water molecule in an ice lattice is linked to four water molecules via HBs [23]. The simplest primitive cell is H-ordered ice Ic. Our theoretical analysis of HB vibrational modes previously found that there are only two kinds of normal modes, which are responsible for the two typical peaks in the translational band [21]. For the strong vibrational mode, the molecule vibrates along the H–O–H angle bisectors with four HBs stretching against the four neighbors—this is the four-bond mode. For the weak vibrational mode, the molecules oscillate along two HBs while the other two remain static—this is the two-bond mode. The vibrational frequency of the four-bond mode is higher than that of the two-bond mode.

Hereafter, we give the analytical solution according to the simulated results. Based on an ideal tetrahedral structure of H-ordered ice Ic, the lattice dynamics analysis is conducted as follows. We simplify a water molecule as a point mass with the position of oxygen taken as the mass center. Thus, there are two atoms per primitive cell [21]. Under the harmonic approximation, the vibrational energy is a sum of eight O–O oscillators expressed as

$$\Psi =\frac{1}{4}k{\displaystyle \sum }_{i=1}^2{\displaystyle \sum }_{j=1}^4{\left(\left|d_{ij}\right|-d_0\right)}^2,$$

(1)

where k is the HB force constant, d₀ is the O–O equilibrium distance, i indicates that there are two point masses in the primitive cell, and j indicates that there are four neighbors. The factor 1/4 is introduced to subtract the double-counted potential energy. Next, we replace the displacement term $\left(\left|d_{ij}\right|-d_0\right)$ with $\mathrm{u}\left(\begin{array}{c}m\\ n\end{array}\right)$. There are six degrees of freedom per primitive cell, which are expressed by the generalized coordinate indices α, β = 1, 2,…, 6. Performing a Taylor series expansion on $\Psi$ produces

$$\Psi =0+\frac{1}{2}{\displaystyle \sum }_{\mathsf{\alpha },\mathsf{\beta }=1}^6{\Phi }_{\mathsf{\alpha }\mathsf{\beta }}{u}_{\mathsf{\alpha }}\left(\begin{array}{c}m\\ n\end{array}\right)u_{\mathsf{\beta }}\left(\begin{array}{c}{m}^{\prime }\\ n^{\prime }\end{array}\right)+ \cdots ,$$

(2)

$${\Phi }_{\mathsf{\alpha }\mathsf{\beta }}=\left(\frac{{\partial }^2\Psi }{\partial u_{\mathsf{\alpha }}\left(\begin{array}{c}m\\ n\end{array}\right)\partial u_{\mathsf{\beta }}\left(\begin{array}{c}{m}^{\prime }\\ n^{\prime }\end{array}\right)}\right)$$

(3)

$$=\left(\begin{array}{cccccc}\frac{4}{3}k& 0& 0& -\frac{4}{3}k& 0& 0\\ 0& \frac{4}{3}k& 0& 0& -\frac{4}{3}k& 0\\ 0& 0& \frac{4}{3}k& 0& 0& -\frac{4}{3}k\\ -\frac{4}{3}k& 0& 0& \frac{4}{3}k& 0& 0\\ 0& -\frac{4}{3}k& 0& 0& \frac{4}{3}k& 0\\ 0& 0& -\frac{4}{3}k& 0& 0& \frac{4}{3}k\end{array}\right).$$

(4)

Considering the included angle cosθ between the HB and vibrational orientations, the constant 4/3 is generated from the rotational symmetry:

$$co{s}^2{\theta }_1+co{s}^2{\theta }_2+co{s}^2{\theta }_3+co{s}^2{\theta }_4=\frac{4}{3}.$$

(5)

Without hydrogen atoms, we can obtain three degenerate eigenvalues of frequencies, which represent three optical vibrational modes. The structure of the diamond is suitable for this type of model. However, the DFT calculation only presents two frequencies. This is due to the breaking of the tetrahedral symmetry by hydrogen atoms, which results in two force constants. Ideal ice Ic contains three intermolecular optical vibrational modes—one is strong and the other two are weak and degenerate [21]. Based on the computational results, the two vibrational frequencies are at 320.76 cm⁻¹ and 229.96 cm⁻¹. Consider the relationship between the frequency and force constant,

$$\begin{array}{c}{\mathsf{\omega }}_1=\sqrt{\frac{1}{m}\left(\frac{8k_s}{3}\right)}=2\pi \times 320.76\times 3\times {10}^{10}{s}^{-1}\\ {\mathsf{\omega }}_2={\mathsf{\omega }}_3=\sqrt{\frac{1}{m}\left(\frac{8k_w}{3}\right)}=2\pi \times 229.96\times 3\times {10}^{10}{s}^{-1}.\hfill \end{array}$$

(6)

The ratio between the force constants of the two vibrational modes is, therefore

$$\frac{k_s}{k_w}={\left(\frac{{\mathsf{\omega }}_1}{{\mathsf{\omega }}_2}\right)}^2\approx 2.$$

(7)

This result shows that the force constant of the four-bond mode is double that of the two-bond mode, with a vibrational energy ratio of $\sqrt{2}$. The two intrinsic HB vibrational modes originate from the local tetrahedral geometry of the ice lattice. This finding has proved to be a general rule for the ice family [18]. For ice phases under atmospheric pressure, such as ices Ih, XI, and Ic, the two characteristic HB peaks can be easily explained by the existence of two kinds of HB vibrational modes. However, for high-pressure ices, the PDOS for their HBs show complicated curves that need to be analyzed individually.

To determine the two kinds of HB vibrational modes in ice III, we employed a self-compiled program to classify the normal modes in the translation band [19]. The criterion for classification is based on the included angle between the vibrational direction and the H–O–H angle bisector. As shown in Figure 2, the HB vibrational modes are categorized into two groups: two-bond modes and four-bond modes. If we sum the two curves, the result largely matches the PDOS. The two curves both show energy splitting at above 200 cm⁻¹. Thus, all the phonons show a superposition band without two characteristic peaks as in ice Ih. Figure 3 presents three example modes, of which the mode at 217 cm⁻¹ is a two-bond vibration, and the other two are four-bond vibrations. To observe the dynamic processes, please see the three animations in the Supplementary Materials.

Notably, the mixed HB vibrational modes lead to a broad overlapping zone below 280 cm⁻¹. Because the HB strength is determined by the bond length and included angle [29], we analyze the intramolecular HB bond lengths and the included angles between HBs around one central molecule in the tetrahedron, as presented by the distributions in Figure 4. The included angles are divided into two groups at $105° \pm 5°$ and 130° (109°28′ in an ideal tetrahedron), and the bond lengths range from 1.82 Å to 1.89 Å (1.79 Å to 1.87 Å in ice Ic [21]). This demonstrates how high pressure deforms the tetrahedral structure of ice III. Similar to our previous discussions on ice V [30], the tetrahedral symmetry in ice III is remarkably distorted by this high pressure. Unlike the six-membered rings that exist under atmospheric pressure, five-, seven-, and eight-membered rings occur in ice III [31]. Moreover, the five-membered rings combine to form bicyclo-heptamers [32]. The broad distributions of bond lengths and deformed angles result in various HB strengths. This is the origin of the energy splitting in both the two-bond and four-bond modes. By comparison, ice IX, the H-ordered counterpart of ice III, has two relatively well-defined HB peaks, as shown in Figure 1. This is because the simple primitive cell of ice IX reduces the phonon complexity, making the spectrum simpler to interpret.

4. Conclusions

To clarify the vibrational characteristics of H-disordered ice III in the far-IR region, we calculated its normal modes and PDOS curve. The good agreement between the PDOS and INS experimental data enabled us to further investigate the HB normal modes. Although its H-ordered counterpart, ice IX, presents two HB peaks, the experimental and simulated results for ice III presented multiple peaks in the HB region. Based on the theory of two HB vibrational modes, we found that the deformation of local tetrahedra by high pressure may be responsible for the energy splitting for both two-bond and four-bond modes. The overlapping of nondegenerate phonons produced a flat vibrational band.

We first found the two kinds of intrinsic HB vibrational modes in H-ordered ice Ic, which has only two molecules per primitive cell [21]. Then, we confirmed that the two distinct HB peaks in ice Ih originated from these two kinds of HB modes [19]. By now, we have investigated almost all the ice phases, thereby determining that this is a general rule for ice because of its local tetrahedral structure [18]. Based on these new findings, we have obtained a deeper understanding of the HB interactions in water ice and their broad role in nature.