A new theoretical model for hexagonal ice, Ih(d), from first principles investigations

1 Introduction

Water and its frozen form in the solid-state, ice, have been widely investigated for many years due to their great importance for science, technology, and life in general. Yet there are still many open questions, and new insights continue to be obtained [1], [2], [3]. Many phases of ice have been observed experimentally; most of these are thermodynamically stable under certain temperature and pressure conditions, but some are only metastable [4], [5], [6], [7], [8], [9], [10], [11], [12]. In addition to these experimental observations, there is a growing set of theoretical ice models ranging from high symmetry structures to amorphous ones [13], [14], [15], [16], [17], [18], [19], [20], [21]. Hexagonal ice is the most predominant phase of ice on Earth, and although a large body of research deals with this structure, still new and unusual features are being discovered. Hexagonal ice has been observed to exist in three different modifications exhibiting different space groups: Ih(a) (space group P6₃cm, no. 185) [22], Ih(b) (space group P6₃mc, no. 186) [23] and Ih(c) (space group P6₃/mmc, no. 194) [6], [7], [8], [10], [24]. Only the first modification Ih(a) has been fully determined structurally including atomic positions, while the other two hexagonal modifications are described using partial occupancy of hydrogen atom sites.

2 Results and discussion

The Ih(c) modification of ice is the most common and most investigated structure [6], [7], [8], [10], [24]; however, its properties are difficult to calculate due to the partial occupancy of atomic sites. In order to overcome this problem, many groups have tried to construct large supercells, employ combinations of various low-symmetry structure models or even perform calculations in space group P1 (i.e. with no symmetry at all) [13], [14], [15], [16], [17], [18], [19], [20], [21], [25], [26]. In this study, we endeavor to obtain a theoretical model that would fit the high-symmetry space group P6₃/mmc, for which we have used the experimental data from the Inorganic Crystal Structure Database (ICSD) [27], [28] as a starting point.

Historically, the crystallographic characterization of the ice Ih(c) structure (space group P6₃/mmc, no. 194) was carried out by two complementary methods, XRD (X-Ray Diffraction) and neutron diffraction. In 1929, Barnes and Bragg [5] determined the arrangement of the oxygen atoms via the XRD method: they form an open structure, and they are expected to be surrounded by four hydrogen atoms in tetrahedral coordination. However, the precise location of the hydrogen atoms is not straightforward to obtain, and typically the models involve partial occupancy to some degree.

Thus, as a starting point for our calculations, we have used an experiment-based model constructed from the neutron diffraction data by Fortes et al. [8], where the proposed structure involves partial occupancy of 50% for hydrogen atoms as shown in Fig. 1a. Using crystallographic analysis and keeping the symmetry-equivalent positions, we have generated a theoretical model, called Ih(d), where we essentially replace pairs of partially occupied H-atom sites by single fully present H-atoms. This structure model was then optimized on ab initio level, without any restrictions of symmetry, using four different functionals: DFT-LDA, GGA-PBE, and hybrid B3LYP and PBE0 in order to establish the robustness of the structure with respect to the choice of the functional. All four functionals yielded the same result; the new structure turned out to be stable under arbitrary small perturbations away from the optimized atom positions without symmetry restrictions. Even when we started from a structure model close to the high-pressure structure ice II, we reached the Ih(d) structure as a closest local minimum. The resulting structure shown in Fig. 1b exhibits a four-fold (tetrahedral) coordination of the oxygen atoms by hydrogen atoms and possesses the symmetries of the space group P6₃/mmc. The cell parameters obtained with the four functionals are shown in Table 1. Clearly, this new model constitutes a kinetically stable, thermodynamically metastable, alternative structure for hexagonal ice in the space group symmetry P6₃/mmc.

Fig. 1:

(a) Model for ice Ih(c) by Fortes et al. [8]; all hydrogen atom sites (marked in white) are partially occupied at 50%; oxygen atoms are marked in red. (b) Visualization of the theoretically proposed structure of ice Ih(d) (space group P6₃/mmc, no. 194) from this work; here, all hydrogen atom sites (marked in white) are fully occupied at 100%.

Next, we have compared our new model with the experimental observations for ice Ih(c). In Fig. 2 we show simulated neutron diffraction patterns for the experiment-based model for ice Ih(c) by Fortes et al. [8] (Exp_ICSD) and for the theoretical model Ih(d) calculated using the LDA functional (Theo_LDA_model). These two models exhibit the same space group (P6₃/mmc, no. 194) but the positions of the hydrogen atoms are different, which is reflected in the differences in the diffraction patterns. We can conclude that the chosen method of neutron diffraction can be used for identifying and distinguishing between these two different models for hexagonal ice in space group P6₃/mmc. In contrast, due to the fact that hydrogen is the lightest element and is difficult to detect via XRD, ice Ih(c) and Ih(d) would not be distinguishable by X-ray diffraction. As mentioned earlier, there is no well-defined accurate way to compute the energy of model systems that exihibit partial occupancy on ab initio level without introducing various (possibly rather drastic) approximations and simplifications. This also applies to the model by Fertes et al., which is based on partial occupancy of the H-atom sites; as a consequence, we cannot compare its energy with the one of our model Ih(d) in a straightforward fashion.

Fig. 2:

Simulated powder neutron diffraction pattern of the model for ice Ih(c) by Fortes et al. [8] (Exp_ICSD), and of the new theoretical model Ih(d) calculated using the LDA functional (Theo_LDA_model) samples, depicted in black and red, respectively.

The full crystallographic data for the two models of hexagonal ice in space group (P6₃/mmc, no. 194), Ih(c) and Ih(d), are presented in Table 2.

A comparison between the new theoretical model of the Ih(d) structure and the experiment-based model of the Ih(c) phase of hexagonal ice shows that the new theoretical model has slightly smaller unit cell parameters than the ones of Ih(c), regardless of the computational approach. However, the major difference between the two models is the partial occupancy of hydrogen in the experiment-based model at the 4f and 12a Wyckoff positions, while our new model shows full occupancy of hydrogen atoms at the 2c and 6g Wyckoff positions. The rather large volume difference between the (experiment based) Ih(c) model and the Ih(d) model makes it unlikely that the Ih(d) model can serve as a structure solution model for ice Ih(c) without modifications. But we expect that this new theoretical model can aid in future calculations of hexagonal ice and its properties, as well as provide an alternative kinetically stable model for structure solution in studies of phase transitions of ice. In particular, the high density of the Ih(d) structure suggests that it would be a good candidate for an intermediary structure on the transition route between hexagonal ice under standard pressure and ice II at elevated pressures.

3 Computational details

The ab initio calculations in this work were performed using the Crystal17 program package [29]. We optimized the H_3-1p1G_gatti1994 basis set for the hydrogen atoms [30], [31], and the O_8-411-towler1994 basis set for the oxygen atom [32], [33]. These basis sets are used in the Linear Combination of Atomic Orbitals (LCAO) method in Crystal17, and they are built from Contracted Gaussian Type Orbitals (CGTO). We added a p orbital to the basis set of the hydrogen atom in order to take polarization effects into account. In this study, the geometry optimizations were performed on the level of Density Functional Theory (DFT), employing four different functionals: LDA (Local Density Approximation), GGA-PBE (Generalized Gradient Approximation following Perdew, Burke and Ernzerhof), the hybrid functional B3LYP (Becke, three-parameter, Lee-Yang-Parr), and the PBE0 functional (mixes the Perdew-Burke-Ernzerhof and Hartree-Fock exchange energy). Predicted and experimentally observed structures were analyzed using the Kplot [34] program and visualized using the Vesta [35] program.