Studies of Hydrogen Bond Vibrations of Hydrogen-Disordered Ice Ic

The hydrogen-disordered structure of ice, Ic, makes it difficult to analyze the vibrational normal modes in the far-infrared region (i.e., the molecular translation band). To clarify the origin of the energy-splitting of hydrogen bond vibrations in this area, a 64-molecule supercell was constructed and calculated using first-principles density functional theory. The results were in good agreement with inelastic neutron scattering experiments and our previous study of a hydrogen-ordered ice Ic model. Assisted by analytic equations, we concluded that the origin of the two hydrogen bond peaks in real ice Ic is consistent with that of hydrogen-ordered ice Ic: the peaks originate from two kinds of normal mode vibration. We categorize the four peaks in the far-infrared region recorded from inelastic neutron scattering experiments as the acoustic peak, the superposition peak, the two-hydrogen bond peak and the four-hydrogen bond peak. We conclude that the existence of two intrinsic hydrogen bond vibration modes represents a general rule among the ice family, except ice X.


Introduction
Water/ice is among the Earth's most abundant materials. Despite its simplicity, at least 18 ice phases have been identified experimentally at specific pressures and temperatures [1][2][3][4][5][6][7][8][9][10][11]. In 2019, a new ice phase, 'ice XVIII', was found to exist at pressures greater than 100 GPa and temperatures above 2000 K [11]. Most ice phases exist in pairs, such as ice Ih and ice XI, which have the same oxygen positions but differ in terms of hydrogen order or disorder. The vibrational spectrum of ice may be detected by infrared (IR) absorption, Raman scattering and inelastic neutron scattering (INS). These spectra can also be simulated using the first-principles density functional theory (DFT) to reveal the lattice structure and the vibrational energies of atoms and molecules. However, the complex hydrogen-disordered structure obstructs phonon analysis. First, we must construct a large supercell to mimic the real lattice structure. However, because the number of simulated normal modes is very large, it is difficult to discern the insight vibration rules. In this study, we demonstrate how to determine the intrinsic vibrational types of the hydrogen bonds (H-bonds) of ice Ic using a comparative method with a hydrogen-ordered ice Ic structure.
König was the first to identify and deduce the structure of ice Ic [12]. The density of ice Ic is nearly the same as ice Ih [13], and its oxygen atoms are arranged in a diamond structure. In other word, the stacking mode of Ic is A, B, C while that of Ih is A, B, A. Ice Ic has a face-centred cubic structure. It can be formed in various ways, such as by the homogeneous crystallisation of pure water droplets at approximately 235 K [14,15], by the decompression of high-density ice at 77 K [16] or from amorphous ice [13,17], ice II, ice III or ice V under high pressure [1]. Under normal conditions, cubic ice Ic is metastable and will gradually convert to hexagonal ice Ih. Whalley first discovered the presence of ice Ic in nature in the refractive halo of sunlight [18]. However, the experimental synthesis of large, high-purity ice Ic crystals is difficult [16,19,20]. A large number of methods have been developed to prepare high-purity ice Ic in recent decades [21][22][23], and disordered accumulation of hydrogen ice with purity in excess of 70% can be prepared [24]. In 2020, Ulivi, del Rosso [25] and Komatsu et al. [26] proposed two methods for the preparation of ice Ic with 100% purity without stacking defects.
Although the ice Ic phase is hydrogen-disordered, Geiger et al. concluded that the ferroelectric hydrogen-ordered structure of ice Ic is the most stable geometric structure based on the results of spectral experiments and first-principles simulations [27]. They proposed that their experimental spectra provided evidence for a partially hydrogenordered structure in ice Ic, but no unique assignment of the hydrogen order was possible. We simulated the vibrational spectrum of the ferroelectric hydrogen-ordered structure of ice Ic [28]. We later used that hydrogen-ordered model to identify two kinds of H-bond vibrational modes and confirmed that this is a general rule for ice structures. In this study, we optimised the hydrogen-disordered phase of ice Ic with a large supercell. Our comparison with the former hydrogen-ordered ice Ic allowed us to verify the two main H-bond vibrational types in the lattice and analyse the physical mechanism with the aid of DFT simulations.

Methods
Despite the hydrogen atoms, the primitive cell can be experimentally measured. Due to its hydrogen-disordered structure under static equilibrium conditions, ice Ic has no periodic repeated unit. Based on the consistency between simulation results and experiments in our previous work [29][30][31][32][33], we constructed a 64-molecule supercell with GenIce [34], in which the hydrogen atoms were rearranged randomly in accordance with the Bernal-Fowler rules [35] to mimic the real structure. The phonon density of states (PDOS) was calculated using the CASTEP module [36] of Materials Studio 6.0 based on the first-principles DFT method. For such a supercell, it was treated as a periodic primitive cell and can be calculated with basic set of plane waves. Meanwhile, the inner hydrogen atoms were disordered so that the calculated normal modes could represent the massive non-degenerate states. Since the electron density has a large gradient change, the exchange-correlation functional of generalized gradient approximation (GGA) was appropriate. Consistent with our previous work [28,31], and to reproduce the relevant features in experimental spectra, we chose the GGA RPBE [37] functional. The calculation of phonons requires high-accuracy geometry optimization. The energy threshold and self-consistent field tolerance were set to 1 × 10 −9 eV per atom for geometry optimization. The energy cut-off was 830 eV, based on our experience. Due to the large supercell, the first Bouillon Zone was small, so that the k-point can be set as 1 × 1 × 1. Norm-conserving pseudopotentials were used to calculate the PDOS with the linear response method. We used our self-written code [31] to classify the main vibration vectors of oxygen, discussed in detail below. Figure 1 shows the good agreement between the PDOS spectra of the hydrogenordered ice Ic and the supercell ice Ic. All four bands have identical widths in the hydrogenordered ice Ic and supercell ice Ic spectra. As the primitive cell of hydrogen-ordered ice Ic has two water molecules, there are only 2 × 3 × 3 − 3 = 15 optical normal modes; thus, the PDOS curve shows sharp pronounced peaks. In contrast with this, the 64-molecule supercell has 64 × 3 × 3 − 3 = 573 normal modes. The vibrational modes are all nondegenerate, and many adjacent modes overlap to form a smooth curve. We focus on the translation band, which comprises the optical and acoustic normal modes. Two sharp H-bond vibrational peaks are obvious in this region and show good agreement with the INS experimental data, as shown in Figure 1.

Results and Discussion
band, which comprises the optical and acoustic normal modes. Two sharp H-bond vibrational peaks are obvious in this region and show good agreement with the INS experimental data, as shown in Figure 1. The simulated PDOS is generated by integrating the phonon signals from the entire Brillouin zone. The INS data also reflect the phonons throughout the Brillouin zone, so the accuracy of PDOS can be verified by the INS spectrum. The PDOS of the supercell ice Ic shows two peaks at 215 and 293 cm −1 , which correspond to those at 225 and 298 cm −1 in the INS spectrum. We have shown that the two characteristic peaks in many ice phases in the translation region are formed from two kinds of H-bond vibrational modes. The consistency between the supercell and hydrogen-ordered ice Ic shows that the main features of H-bond vibrations can be deduced from the hydrogen-ordered structure of ice Ic. In the study of hydrogen-ordered ice Ic, we found that the two peaks originate from three normal modes in total [28]. The molecular configurations of any ice phase are subject to the Bernal-Fowler rules, in which one water molecule links four neighbours via H-bonds to form a tetrahedral geometry [35]. For the strong peak, a molecule vibrates along the H-O-H angle bisectors with all four neighbours. Thus, all four H-bonds participate in the vibration, which can thus be named the 'four-bond' mode. For the weak peak, a molecule vibrates along two H-bonds while the other two remain static, so we refer to this kind as the 'two-bond' mode. In ice Ic, there are two degenerate modes of the two-bond kind due to symmetry in the basal plane. Hence, there are three normal modes in total.
According to the simulation results of hydrogen-ordered ice Ic, we analyse the dynamic process as follows.
The water molecule is treated as a point mass, and the oxygen is regarded as the approximate mass center. Details for this model can be seen from Reference [31]. In that work, we consider the O-O-O bending energies to explain the H-bond energy splitting. However, without the H-O-H bending, the angle change between two H-bonded O-O- The simulated PDOS is generated by integrating the phonon signals from the entire Brillouin zone. The INS data also reflect the phonons throughout the Brillouin zone, so the accuracy of PDOS can be verified by the INS spectrum. The PDOS of the supercell ice Ic shows two peaks at 215 and 293 cm −1 , which correspond to those at 225 and 298 cm −1 in the INS spectrum. We have shown that the two characteristic peaks in many ice phases in the translation region are formed from two kinds of H-bond vibrational modes. The consistency between the supercell and hydrogen-ordered ice Ic shows that the main features of H-bond vibrations can be deduced from the hydrogen-ordered structure of ice Ic. In the study of hydrogen-ordered ice Ic, we found that the two peaks originate from three normal modes in total [28]. The molecular configurations of any ice phase are subject to the Bernal-Fowler rules, in which one water molecule links four neighbours via H-bonds to form a tetrahedral geometry [35]. For the strong peak, a molecule vibrates along the H-O-H angle bisectors with all four neighbours. Thus, all four H-bonds participate in the vibration, which can thus be named the 'four-bond' mode. For the weak peak, a molecule vibrates along two H-bonds while the other two remain static, so we refer to this kind as the 'two-bond' mode. In ice Ic, there are two degenerate modes of the two-bond kind due to symmetry in the basal plane. Hence, there are three normal modes in total.
According to the simulation results of hydrogen-ordered ice Ic, we analyse the dynamic process as follows.
The water molecule is treated as a point mass, and the oxygen is regarded as the approximate mass center. Details for this model can be seen from Reference [31].  where k is the H-bond force constant, d is the O-O distance at equilibrium, ( d ij − d) represents the relative displacement between two adjacent molecules, i indicates the two molecules in one cell, and j indicates four neighbours. The factor 1/4 is used to subtract the double-counted potential energy. Regardless of the hydrogen atoms, the two-point masses in one primitive cell have six degrees of freedom. We use index α = 1, 2, . . . , 6 for the generalized coordinates. Expanding Φ as a function of quadratic terms of molecule displacements u l i where l, l indicates different cells. Assuming that the included angle of the H-bond and the vibrational direction are cosθ, the vibration energy is proportional to cos 2 θ. Thus, due to rotational symmetry, we obtain A locally tetrahedral structure, such as diamond, has three energy-degenerate optical normal modes. In fact, the tetrahedral symmetry of ice is broken by hydrogen atoms, so the vibrational normal modes show energy level splitting. To consider the two types of vibrational modes, we set a coefficient c and obtain The Lagrangian of each cell is and the dynamic matrix is Based on the simulations, the hydrogen-ordered model comprises two kinds of molecular translational mode. According to the results of hydrogen-ordered ice Ic, the two frequencies are 320.76 and 229.96 cm −1 .
These results indicate that the contributions of H-bond force constants of the fourbond mode are double those of the two-bond mode and that the vibrational energy ratio between the two modes is approximately √ 2. The results verify that the H-bond energy level splitting is due to the existence of two kinds of H-bond vibrational mode in the tetrahedral structure of ice. It is then natural to deduce that the two characteristic peaks in the supercell ice Ic result from these two kinds of modes. We use a simple model of intermolecular connecting spring to analyze the two kinds of vibrational modes in Reference [28]. We obtained the simple ratio √ 2. In this work, we obtained the vibrational frequencies from first-principles DFT calculations and verified that the classification of four-bond and two-bond modes is rational.
The stronger vibrational band ranging from 160 to 320 cm −1 comes from H-bondrelated monomolecular vibrations, so we classified the normal modes (gamma point only) into two categories and fitted the signals into two curves, as shown in Figure 2, which shows that the sum curve is consistent with the PDOS. The principle of classification is based on the vibrational direction of the maximum amplitude of oxygen in one mode. The source code of this self-written code was published in Reference [31]. As shown in Figure 2, the four-bond modes have higher energies at 293 cm −1 , whilst the greatest distribution of the two-bond modes lies at 215 cm −1 . However, it is notable that the two kinds of modes overlap below 280 cm −1 . As all of the modes are non-degenerate, many coupled modes are located between the peaks in the four-bond and two-bond modes, so our simple self-written code may be unable to distinguish them clearly. In addition, the use of oxygen as the point-mass center may produce some errors. The main conclusions from Figure 2 are that the strong peak at 293 cm −1 in the PDOS comes from four-bond mode vibrations and that the weak peak at 215 cm −1 comes from two-bond mode vibrations. Thus, based on the hydrogen-ordered model, we verify that the real ice Ic phase follows the same rule as that derived from the hydrogen-ordered model. Considering our series of investigations into ice phases, we conclude that this is a general rule among the ice family, except ice X [28][29][30][31][32][33][39][40][41][42][43][44][45][46]. The vibrational modes discussed above comprise the vibration of a single molecule against its neighbors. In contrast, the modes in the lower-energy regions are cluster vibrations, in which many molecules bonded together present a collective vibration. In this case, the number of participating H-bonds decreases while the mass of the cluster increases. For a simple harmonic motion, the vibrational frequency = ⁄ ; therefore, the cluster vibrations present lower energies than the monomolecular vibrations, and many of them may overlap with acoustic phonons. This is the origin of peaks below 160 cm −1 in the PDOS. Figure 3 portrays several typical examples of the dynamic processes of The vibrational modes discussed above comprise the vibration of a single molecule against its neighbors. In contrast, the modes in the lower-energy regions are cluster vibrations, in which many molecules bonded together present a collective vibration. In Crystals 2021, 11, 668 6 of 9 this case, the number of participating H-bonds decreases while the mass of the cluster increases. For a simple harmonic motion, the vibrational frequency υ = √ k/m; therefore, the cluster vibrations present lower energies than the monomolecular vibrations, and many of them may overlap with acoustic phonons. This is the origin of peaks below 160 cm −1 in the PDOS. Figure 3 portrays several typical examples of the dynamic processes of these normal modes. Figure 3a,b shows two cluster vibrations. Figure 3c In the lowest-energy translation region, a distinct peak was observed at 57 cm −1 via the INS experiments. This peak represents phonons that belong to acoustic branches and corresponds to the peak at 95 cm −1 in the simulated PDOS curve. Another INS peak at 153 cm −1 should correspond to the PDOS peak at 168 cm −1 . We attribute this peak to the superposition of acoustic phonons with optical phonons, as seen in the dispersion curves from our previous study of hydrogen-ordered ice Ic [28]. Thus, from low to high energy, the four pronounced peaks in Figure 1 in the translation band of ice Ic are attributed to the acoustic peak, the superposition peak, the two-bond peak and the four-bond peak.

Conclusions
We previously confirmed the origin of energy splitting in the H-bond vibrations of ice Ic via simulation of a hydrogen-ordered ice Ic model [28]. However, the real ice Ic crystal is hydrogen-disordered. To simulate the real structure, we constructed a 64-molecule supercell. The calculated PDOS shows good agreement with the INS spectrum and the hydrogen-ordered model of ice Ic, implying that the origin of the two H-bond peaks in real ice Ic are consistent with that inferred from simulations of hydrogen-ordered ice Ic. We confirmed the existence of two kinds of intrinsic H-bond vibration modes in hydrogen-ordered ice Ic due to the local tetrahedral structure. The strong mode was named the 'four-bond' mode because the molecule vibrates against four neighbours, so that all In the lowest-energy translation region, a distinct peak was observed at 57 cm −1 via the INS experiments. This peak represents phonons that belong to acoustic branches and corresponds to the peak at 95 cm −1 in the simulated PDOS curve. Another INS peak at 153 cm −1 should correspond to the PDOS peak at 168 cm −1 . We attribute this peak to the superposition of acoustic phonons with optical phonons, as seen in the dispersion curves from our previous study of hydrogen-ordered ice Ic [28]. Thus, from low to high energy, the four pronounced peaks in Figure 1 in the translation band of ice Ic are attributed to the acoustic peak, the superposition peak, the two-bond peak and the four-bond peak.

Conclusions
We previously confirmed the origin of energy splitting in the H-bond vibrations of ice Ic via simulation of a hydrogen-ordered ice Ic model [28]. However, the real ice Ic crystal is hydrogen-disordered. To simulate the real structure, we constructed a 64molecule supercell. The calculated PDOS shows good agreement with the INS spectrum and the hydrogen-ordered model of ice Ic, implying that the origin of the two H-bond peaks in real ice Ic are consistent with that inferred from simulations of hydrogen-ordered ice Ic. We confirmed the existence of two kinds of intrinsic H-bond vibration modes in hydrogen-ordered ice Ic due to the local tetrahedral structure. The strong mode was named the 'four-bond' mode because the molecule vibrates against four neighbours, so that all four H-bonds work together. The weak mode was named the 'two-bond' mode because the molecule vibrates along only two H-bonds, while the other two remain static. Based on the two kinds of H-bond vibrations, we classified the H-bond modes into two corresponding groups and obtained consistent results. Thus, the most intense H-bond vibrations correspond to these two kinds of modes, giving rise to two distinct peaks.
Due to the simplicity of the primitive cell of hydrogen-ordered Ic, there are only one strong H-bond vibration mode and two degenerate weak H-bond vibration modes. For the 64-molecule hydrogen-disordered model, there are 3 × 64 − 3 = 189 normal modes in the translational band. As the normal modes are all non-degenerate due to the hydrogendisordered structure, many coupled modes can be found away from the two H-bond peaks. Thus, a complex set of vibration modes coexist in a hydrogen-disordered ice Ic lattice. In the region below 160 cm −1 of real Ic, there are many cluster vibration modes, which have much lower vibrational energies. However, the PDOS curve shows a very similar pattern to the hydrogen-ordered model. This means that the two kinds of H-bond vibration modes represent the major feature of molecular vibrations in the far-IR region.
We classify the four distinct peaks in the INS spectrum of ice Ic as the acoustic peak, the superposition peak, the two-bond peak and the four-bond peak, as in our previous investigations into ice Ih [31]. We conclude that the existence of two intrinsic H-bond vibration modes represents a general rule among the ice family, except ice X.