Due to their great importance in science, technology, and the life sciences, water and ice have been extensively investigated over many years. In particular, hexagonal ice Ih has been of great interest since it is the most common form of ice, and several modifications, Ih(a), Ih(b) and Ih(c) are known, whose structural details are still under discussion. In this study, we present an alternative theoretical model, called Ih(d), for the hexagonal ice modification in space group P63/mmc (no. 194), based on first-principles calculations that have been performed using DFT-LDA, GGA-PBE, and hybrid B3LYP and PBE0 functionals.
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 , , . 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 , , , , , , , , . In addition to these experimental observations, there is a growing set of theoretical ice models ranging from high symmetry structures to amorphous ones , , , , , , , , . 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 P63cm, no. 185) , Ih(b) (space group P63mc, no. 186)  and Ih(c) (space group P63/mmc, no. 194) , , , , . 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 , , , , ; 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) , , , , , , , , , , . In this study, we endeavor to obtain a theoretical model that would fit the high-symmetry space group P63/mmc, for which we have used the experimental data from the Inorganic Crystal Structure Database (ICSD) ,  as a starting point.
Historically, the crystallographic characterization of the ice Ih(c) structure (space group P63/mmc, no. 194) was carried out by two complementary methods, XRD (X-Ray Diffraction) and neutron diffraction. In 1929, Barnes and Bragg  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. , 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 P63/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 P63/mmc.
|Unit cell parameters||a=3.989 Å||a=3.959 Å||a=3.967 Å||a=3.942 Å|
|c=6.542 Å||c=6.484 Å||c=6.504 Å||c=6.461 Å|
|V=90.19 Å3||V=88.04 Å3||V=88.68 Å3||V=86.98 Å3|
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.  (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 (P63/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 P63/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.
The full crystallographic data for the two models of hexagonal ice in space group (P63/mmc, no. 194), Ih(c) and Ih(d), are presented in Table 2.
|Partial occupancy model adjusted to experiment Ih(c)||Proposed full occupancy model Ih(d)|
|Space group||P63/mmc (no. 194)||P63/mmc (no. 194)|
|Unit cell parameters||a=4.497 Å|
|Atomic positions||O1 (1/3, 2/3, 0.06140)||O1 (2/3, 1/3, 0.93708)|
|Wyckoff position; occupancy||4f; 1.0||4f; 1.0|
|Atomic positions||H1 (1/3, 2/3, 0.19640)||H1 (2/3, 1/3, 3/4)|
|Wyckoff position; occupancy||4f; 0.5||2c; 1.0|
|Atomic positions||H2 (0.45100, 2x, 0.01790)||H2 (1/2, 1/2, 0)|
|Wyckoff position; occupancy||12a; 0.5||6g; 1.0|
O1: oxygen atom 1; H1: hydrogen atom 1; H2: hydrogen atom 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 . We optimized the H_3-1p1G_gatti1994 basis set for the hydrogen atoms , , and the O_8-411-towler1994 basis set for the oxygen atom , . 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  program and visualized using the Vesta  program.
Dedicated to:Professor Arndt Simon on the occasion of his 80th birthday.
The authors thank Professor R. Dovesi, Professor K. Doll and Crystal Solutions for software support. The authors are also grateful to B. Matovic, M. Prekajski and J. Lukovic for valuable discussions. This work has been supported by the grant No. 45012 from the Ministry of Education, Science and Technological Development of the Republic of Serbia.
 T. Bartels-Rausch, V. Bergeron, J. H. E. Cartwright, R. Escribano, J. L. Finney, H. Grothe, P. J. Gutiérrez, J. Haapala, W. F. Kuhs, J. B. C. Pettersson, S. D. Price, C. I. Sainz-Díaz, D. J. Stokes, G. Strazzulla, E. S. Thomson, H. Trinks, N. Uras-Aytemiz, Rev. Modern Phys.2012, 84, 885–944.10.1103/RevModPhys.84.885Search in Google Scholar
 G. Bergerhoff, I. D. Brown, Crystallographic Databases, International Union of Crystallography, Chester (U.K.) 1987.Search in Google Scholar
 R. Dovesi, A. Erba, R. Orlando, C. M. Zicovich-Wilson, B. Civalleri, L. Maschio, M. Rérat, S. Casassa, J. Baima, S. Salustro, B. Kirtman, Wiley Interdiscipl. Rev.: Comp. Mol. Sci.2018, 8, e1360.Search in Google Scholar
 R. Hundt, Kplot, A Program for Plotting and Analyzing Crystal Structures, Technicum Scientific Publishing, Stuttgart (Germany) 2016.Search in Google Scholar
©2020 Walter de Gruyter GmbH, Berlin/Boston