Skip to content
Publicly Available Published by De Gruyter November 30, 2019

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

  • Dušica Jovanović , Dejan Zagorac , J. Christian Schön EMAIL logo , Branislav Milovanović and Jelena Zagorac EMAIL logo


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.

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 P63cm, no. 185) [22], Ih(b) (space group P63mc, no. 186) [23] and Ih(c) (space group P63/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 P63/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 P63/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 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.

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 P63/mmc, no. 194) from this work; here, all hydrogen atom sites (marked in white) are fully occupied at 100%.
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 P63/mmc, no. 194) from this work; here, all hydrogen atom sites (marked in white) are fully occupied at 100%.

Table 1:

Optimized cell parameters of the proposed theoretical model of ice Ih(d), calculated using four different functionals, DFT-LDA, GGA-PBE, and hybrid B3LYP and PBE0, respectively.

Computational method
Unit cell parametersa=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 Å3V=88.04 Å3V=88.68 Å3V=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. [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 (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.

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.
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 (P63/mmc, no. 194), Ih(c) and Ih(d), are presented in Table 2.

Table 2:

Unit cell parameters and atom positions of the partial occupancy model and of the theoretically predicted model with full occupancy for the ice Ih(c) and Ih(d) modification, respectively.

Partial occupancy model adjusted to experiment Ih(c)Proposed full occupancy model Ih(d)
Space groupP63/mmc (no. 194)P63/mmc (no. 194)
Unit cell parametersa=4.497 Å

c=7.322 Å

V=128.27 Å3
a=3.959 Å

c=6.484 Å

V=88.04 Å3
Atomic positionsO1 (1/3, 2/3, 0.06140)O1 (2/3, 1/3, 0.93708)
Wyckoff position; occupancy4f; 1.04f; 1.0
Atomic positionsH1 (1/3, 2/3, 0.19640)H1 (2/3, 1/3, 3/4)
Wyckoff position; occupancy4f; 0.52c; 1.0
Atomic positionsH2 (0.45100, 2x, 0.01790)H2 (1/2, 1/2, 0)
Wyckoff position; occupancy12a; 0.56g; 1.0
  1. 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 [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.

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.


[1] N. H. Fletcher, Rep. Prog. Phys.1971, 34, 913–994.10.1088/0034-4885/34/3/303Search in Google Scholar

[2] E. A. Zheligovskaya, G. G. Malenkov, Russ. Chem. Rev.2006, 75, 57–76.10.1070/RC2006v075n01ABEH001184Search in Google Scholar

[3] 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

[4] P.W. Bridgman, Proc. Am. Acad. Arts Sci.1912, 47, 441–558.10.2307/20022754Search in Google Scholar

[5] W. H. Barnes, W. H. Bragg, Proc. Royal Soc. London. Series A Math1929, 125, 670–693.10.1098/rspa.1929.0195Search in Google Scholar

[6] A. Goto, T. Hondoh, S. Mae, J. Chem. Phys.1990, 93, 1412–1417.10.1063/1.459150Search in Google Scholar

[7] L. G. Dowell, A. P. Rinfret, Nature1960, 188, 1144–1148.10.1038/1881144a0Search in Google Scholar

[8] A. D. Fortes, I. G. Wood, D. Grigoriev, M. Alfredsson, S. Kipfstuhl, K. S. Knight, R. I. Smith, J. Chem. Phys.2004, 120, 11376–11379.10.1063/1.1765099Search in Google Scholar PubMed

[9] B. Kamb, A. Prakash, Acta Crystallogr.1968, B24, 1317–1327.10.1107/S0567740868004231Search in Google Scholar

[10] K. Rottger, A. Endriss, J. Ihringer, S. Doyle, W. F. Kuhs, Acta Crystallogr.1994, B50, 644–648.10.1107/S0108768194004933Search in Google Scholar

[11] Y. Wang, H. Zhang, X. Yang, S. Jiang, A. F. Goncharov, J. Chem. Phys.2018, 148, 044508.10.1063/1.5017507Search in Google Scholar PubMed

[12] T. Kubo, W. B. Durham, L. A. Stern, S. H. Kirby, Science2006, 311, 1267–1269.10.1126/science.1121296Search in Google Scholar PubMed

[13] C. Lee, D. Vanderbilt, K. Laasonen, R. Car, M. Parrinello, Phys. Rev. B1993, 47, 4863–4872.10.1103/PhysRevB.47.4863Search in Google Scholar PubMed

[14] A. D. Fortes, I. G. Wood, J. P. Brodholt, L. Vočadlo, J. Chem. Phys. 2003, 119, 4567–4572.10.1063/1.1593630Search in Google Scholar

[15] G. A. Tribello, B. Slater, C. G. Salzmann, J. Am. Chem. Soc.2006, 128, 12594–12595.10.1021/ja0630902Search in Google Scholar PubMed

[16] Z. Raza, D. Alfè, C. G. Salzmann, J. Klimeš, A. Michaelides, B. Slater, Phys. Chem. Chem. Phys.2011, 13, 19788–19795.10.1039/c1cp22506eSearch in Google Scholar PubMed

[17] E. A. Engel, B. Monserrat, R. J. Needs, Phys. Rev. X2015, 5, 021033.10.1103/PhysRevX.5.021033Search in Google Scholar

[18] S. Casassa, M. Calatayud, K. Doll, C. Minot, C. Pisani, Chem. Phys. Lett.2005, 409, 110–117.10.1016/j.cplett.2005.04.068Search in Google Scholar

[19] B. Santra, J. Klimeš, A. Tkatchenko, D. Alfè, B. Slater, A. Michaelides, R. Car, M. Scheffler, J. Chem. Phys.2013, 139, 154702.10.1063/1.4824481Search in Google Scholar PubMed

[20] S. Jenkins, I. Morrison, J. Phys.: Condens. Matter2001, 13, 9207–9229.10.1088/0953-8984/13/41/312Search in Google Scholar

[21] A. Hermann, P. Schwerdtfeger, Phys. Rev. Lett.2008, 101, 183005.10.1103/PhysRevLett.101.183005Search in Google Scholar PubMed

[22] J. D. Bernal, R. H. Fowler, J. Chem. Phys.1933, 1, 515–548.10.1063/1.1749327Search in Google Scholar

[23] R. E. Rundle, J. Chem. Phys.1953, 21, 1311–1311.10.1063/1.1699206Search in Google Scholar

[24] A. D. Fortes, Powder Diffr.2015, 30, 149–157.10.1017/S0885715615000123Search in Google Scholar

[25] T. K. Hirsch, L. Ojamäe, J. Phys. Chem. B2004, 108, 15856–15864.10.1021/jp048434uSearch in Google Scholar

[26] J. J. Shephard, B. Slater, P. Harvey, M. Hart, C. L. Bull, S. T. Bramwell, C. G. Salzmann, Nat. Phys.2018, 14, 569–572.10.1038/s41567-018-0094-zSearch in Google Scholar

[27] G. Bergerhoff, I. D. Brown, Crystallographic Databases, International Union of Crystallography, Chester (U.K.) 1987.Search in Google Scholar

[28] D. Zagorac, H. Muller, S. Ruehl, J. Zagorac, S. Rehme, J. Appl. Crystallogr.2019, 52, 918–925.10.1107/S160057671900997XSearch in Google Scholar PubMed PubMed Central

[29] 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

[30] C. Gatti, V. R. Saunders, C. Roetti, J. Chem. Phys.1994, 101, 10686–10696.10.1063/1.467882Search in Google Scholar

[31] M. Corno, C. Busco, B. Civalleri, P. Ugliengo, Phys. Chem. Chem. Phys.2006, 8, 2464–2472.10.1039/b602419jSearch in Google Scholar PubMed

[32] M. D. Towler, N. L. Allan, N. M. Harrison, V. R. Saunders, W. C. Mackrodt, E. Aprà, Phys. Rev. B1994, 50, 5041–5054.10.1103/PhysRevB.50.5041Search in Google Scholar PubMed

[33] D. Zagorac, J. C. Schön, J. Zagorac, M. Jansen, RSC Adv.2015, 5, 25929–25935.10.1039/C4RA16574HSearch in Google Scholar

[34] R. Hundt, Kplot, A Program for Plotting and Analyzing Crystal Structures, Technicum Scientific Publishing, Stuttgart (Germany) 2016.Search in Google Scholar

[35] K. Momma, F. Izumi, J. Appl. Crystallogr.2011, 44, 1272–1276.10.1107/S0021889811038970Search in Google Scholar

Received: 2019-10-15
Accepted: 2019-10-25
Published Online: 2019-11-30
Published in Print: 2020-02-25

©2020 Walter de Gruyter GmbH, Berlin/Boston

Downloaded on 1.12.2023 from
Scroll to top button