Computing Investigations of Molecular and Atomic Vibrations of Ice IX

The normal modes of ice IX were investigated using the CASTEP code package, which is based on density functional theory. We found that the translational modes could be divided into three categories: four-bond vibrations, which possessed the highest energy; two-bond vibrations, which possessed the medium energy; and cluster vibrations with the lowest energy. The former two categories represent monomers vibrating against neighbors and present as two distinct peaks in many ice phases recorded in inelastic neutron-scattering experiments. It is typically difficult to assign the molecular vibration peaks in the far infrared region. The method we developed to analyze the normal modes, especially in the translation band of ice IX, provided physical insights into the vibrational spectrum of ice.


INTRODUCTION
Water is one of the most important materials in the human body and the most valuable resource on Earth. 1 The solid state of water, ice, has been experimentally confirmed to exist in more than 17 crystalline and amorphous phases, depending on pressure and temperature. 2−11 These phases can be divided into two main categories according to whether their structure is hydrogen long-range disordered or ordered. At present, there are six known disordered/ordered pair phases of ice: Ih/XI, III/IX, V/XIII, VI/XV, VII/VIII, and XII/XIV. 12 This work focused on interpreting the vibrational spectra of ice IX.
Ice III has a hydrogen-disordered structure and exists at 24− 44°C and 2.4−3.4 kbar. Whalley et al. found that, at 65−108°C , ice III can gradually transform into a new hydrogenordered phase that is designated ice IX. 8 Many infrared (IR) absorption, 13 Raman scattering, 14 and neutron scattering 15 experiments have been conducted to characterize ice III and ice IX. Garg and other groups of researchers found by Raman analysis that pressurizing ice Ih at 150 K transformed it into ice IX. 16−19 We performed a series of studies 20−27 of the vibrational spectra of ice by developing a method to analyze its lattice vibrations. We constructed a dual-track approach to verify the validity for the computational phonon spectrum: collate the simulated spectrum with inelastic neutron scattering (INS) experiments and assign the photon absorption/scattering peaks according to the calculated normal vibration frequencies. 28 The two distinct peaks in the translational band were found to correspond to two molecular vibration modes. In this study, we continued to investigate ice IX and confirmed this finding. In addition, we analyzed the optic phonons in the far IR region and provided physical insights into molecular vibrations in the translational band of ice IX.

RESULTS AND DISCUSSION
The computed spectra of Raman scattering, IR absorption, and phonon density of states (PDOS) are shown in Figure 1, which are divided into four parts according to the four separate vibrational regions. Due to the wide range of the intensity in different regions, we adjusted their proportions for comparison. As the interactions between photons and phonons are restricted near the Brillouin zone (BZ) center, the peaks of Raman and IR spectra were close to the frequencies of normal modes. The spectrum of INS was proportional to the calculated PDOS curve because INS could collect all the phonon signals throughout the first BZ. 20 Comparisons between PDOS and INS and between the normal modes and Raman/IR peaks are shown in Table 1. Because of the hydrogen-disordered nature of ice III, their vibrational spectrum does not show sharp features as ice IX. We supplemented the discussions in the corresponding context.
A primitive cell contained 12 molecules. Thus, there were 12 × 3 × 3 − 3 = 105 optic normal modes, and 33 of them were in the translation band at 66−307 cm −1 . From Figure 1, we can see that the PDOS curve had three main bands at approximately 100, 200, and 300 cm −1 . A comparison of the obtained INS measurements with those reported by Li shows that the set exactly fitted the experimental curve. 29 A detailed comparison of peaks is given in Table 1.
It is typically difficult to assign the peaks of a vibrational spectrum in the fingerprint region. To clarify the physical features of these spectra, we needed to review our previous study of the simplest phase of the ice family, the hydrogenordered structure of ice Ic. We found two kinds of intrinsic translational modes in this lattice. 21 In the strong mode, the molecule vibrates against with four connecting H-bonds, which is denoted the four-bond mode. In the weak mode, two perpendicular degenerate modes involve only two oscillating H bonds, and these are denoted the two-bond mode. By treating the two H-bonded molecules as a spring, a simple harmonic oscillator model yielded a strength ratio of 2 . Later, we also observed this phenomenon in ices XIV, XVI, XVII, VII, and VIII 22−25 and concluded that the two distinct peaks of ice in the translation band are from two types of vibrational modes. In this work, we not only confirmed this rule but also determined the origin of the third lower peak.
As shown in Figure 2, the normal mode at 307 cm −1 was a typical four-bond mode, in which all the molecules in one primitive cell vibrated along its HOH angular bisector. This was the highest frequency found for the translational modes, and there were nine such strong modes from 265 to 307 cm −1 , as listed in Table 1.
The mode at 245 cm −1 was an example of two-bond modes that possessed a medium energy of 119−245 cm −1 . These two kinds of modes corresponded to vibrations of monomer and thus were independent of movements. The dynamic vibration of the mode at 116 cm −1 (see Supporting Information) was very similar to that of the strong mode in which all molecules vibrated along their bisector. However, this was a collective vibration in which three molecules combined to form a cluster. The inner H bonds did not stretch, so the vibrational frequency was much lower than that of the four-bond mode. Based on the analysis above, we plotted the distributions of the above three categories to mimic the PDOS curve, as shown in Figure 3. Compared with the insertion of PDOS, one can see a minor discrepancy due to phonon dispersion throughout the BZ as explained in our previous work. 28 There were 36 normal modes in the intermolecular libration band at 580−1003 cm −1 , and the bandwidth was 423 cm −1 , which was in good agreement with the bandwidth at 440 cm −1 described in Li's report. 18 There were three types of vibration modes: rocking or rotation of the whole molecule around an axis perpendicular to the molecular plane; twisting or rotation of the molecule around an axis coincident with the HOH angle bisector; and wagging or rotation of the molecule around an axis in the molecular plane, perpendicular to the bisector of the HOH angle. In an optic mode, the mass center of a primitive cell remains static, while the molecules within it may vibrate differently. Two typical modes are illustrated in Figure 2.
All the molecules displayed rocking vibrations at 580 cm −1 and twisting at 1003 cm −1 . As all the molecules vibrated in wagging mode at 771 cm −1 , we cannot display them all in Figure 2. This normal mode analysis clearly demonstrated that the order of strength was rocking, wagging, and finally twisting.
Three sharp peaks were observed in the libration region of the INS spectrum by Li et al. at 88.9 (711 cm −1 ), 116 (928 cm −1 ), and 90 meV (720 cm −1 ). 18 These peaks were observed in the PDOS spectrum at 702, 720, and 916 cm −1 , respectively. However, the Raman/IR data of ice IX were missing in the libration band. For ice III, four IR peaks were detected in this region by Bertie et al. at 600, 660, 734, and 925 cm −1 . 30 These peaks can match well with normal modes at 593, 662, 737, and 934 cm −1 , respectively.
In the bending band, there were 12 normal modes at 1639− 1716 cm −1 . All the monomer vibrations were the same. However, the collective vibrations could be classified as inphase or out-of-phase vibration modes. The mode at 1639 cm −1 with the lowest energy in this band had an out-of-phase mode, while the strongest mode at 1716 cm −1 had an in-phase vibration. On the contrary, an increasing energy trend from inphase bending to out-of-phase bending was observed in ice XIV. 22 Li et al. found a peak at 210 meV (1680 cm −1 ) in the INS spectrum, 18 which matched well with the PDOS curve. Although the simulated Raman and IR spectra both contained some peaks in this band, their intensities were too weak to permit proper detection. We also found an IR medium peak of ice III at 1690 cm −1 , 30 which can be observed in normal modes at 1702 cm −1 .
There were 24 normal vibration modes in the stretching band at 3177−3391 cm −1 . In this region, the O−H covalent bonds of each molecule showed symmetric stretching or asymmetric stretching. Two examples are illustrated in Figure 2 to show that the energy increasing trend was from symmetric stretching to asymmetric stretching, consistent with trends in other ice phases. We also found many modes in this band with one O−H bond exhibiting "isolated stretching", while the other remained static, as was reported by Whale et al. 31 This is very common in high-pressure ices in our previous studies of ice XIV, XVI, XVII, and VII. We attributed this phenomenon

ACS Omega
Article (NM). 12 According to their work, two main Raman peaks were calculated at 3154 and 3243 cm −1 . Furthermore, our results were at 3177 and 3285 cm −1 . Compared with the experimental Raman spectrum, we found that the weak peak calculated by MB-MD matched well with the experimental data. In contrast, our work of the second peak at 3285 cm −1 was better. Our experiences indicated that the lattice dynamics simulation based on harmonic approximation is good enough to present the vibrational spectrum of ices. The INS spectrum had a peak at 418 meV (3344 cm −1 ), 18

CONCLUSIONS
In summary, we investigated the vibrational spectrum of ice IX using a DFT-based method. We analyzed the normal modes in this model, especially in the translational band. We divided the translational modes into three categories: four-bond vibrations, which possessed the highest energy; two-bond vibrations, which possessed the medium energy; and cluster vibrations, which had the lowest energy. This may be a general rule among ice phases. In ice XV, VII, and VIII, the cluster modes were relative vibrations between two sublattices. The four-bond and twobond modes constituted the two distinct peaks seen in many ices, while the cluster modes often had a hump in lower bands or merged with acoustic phonons.
The normal modes of ice IX displayed significant regularities. In the libration band, the energy-increasing trend was from rocking to wagging to twisting; in the bending band, this trend was from out-of-phase bending to in-phase bending; and in the stretching band, it was from symmetric stretching to asymmetric stretching. Thus, ice IX was confirmed to have a typical phase that could represent the nature of ice in terms of molecular and atomic vibrations.

ACS Omega
Article 4. COMPUTATIONAL METHOD Using CASTEP, 33 a package of first-principles density functional theory (DFT), we performed geometrical optimization and phonon calculation of ice IX with the method of linear response. The linear response was used to calculate the PDOS, which is one of the most popular methods to provide an analytical way of computing the harmonic approximation results. 34 As the electron density has a large distribution, we chose generalized gradient approximation (GGA) for the ice crystal. However, predictions made by different exchange− correlation (XC) functionals can have large disagreement and may not be accurate enough as mentioned by Gillan et al. 35 To test the (XC) functional, we tried RPBE, PBE, PBE + D, PW91, PW91 + D, WC, and PBESOL for geometric optimization of ice IX. The results showed that the densities are 1.01, 1.15, 1.28, 1.18, 1.58, 1.24, and 1.28 g/cm 3 , respectively. Although the densities from PBE and PW91 results are close to the experiment (1.17 g/cm 3 ), the consequent simulations of the strongest translation modes are at 370 and 409 cm −1 , quite bigger than the experimental data of 317 cm −1 . Since the best matchable frequency was at 307 cm −1 of RPBE, we thus chose XC RPBE 36 for this work. The norm-conserving pseudopotential was used to calculate the PDOS. The energy convergence and SCF tolerance were set at 1 × 10 −9 eV/atom. The energy cutoff was 750 eV, and the k-point sampling was 2 × 2 × 2. The hydrostatic pressure was set to 0.3 GPa.

Notes
The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS
The numerical calculations were done on the supercomputing system in the Supercomputing Center, Shandong University, Weihai.