Stability-Ranking of Crystalline Ice Polymorphs Using Density-Functional Theory

In this work, we consider low-enthalpy polymorphs of ice, predicted previously using a modified basin-hopping algorithm for crystal-structure prediction with the TIP4P empirical potential at three pressures (0, 4 and 8 kbar). We compare and (re)-rank the reported ice polymorphs in order of energetic stability, using high-level quantum-chemical calculations, primarily in the guise of sophisticated Density-Functional Theory (DFT) approaches. In the absence of applied pressure, ice Ih is predicted to be energetically more stable than ice Ic, and TIP4P-predicted results and ranking compare well with the results obtained from DFT calculations. However, perhaps not unexpectedly, the deviation between TIP4Pand DFT-calculated results increases with applied external pressure.


Introduction
The study of water's properties, regardless of the particular condensed phase thereof, is of paramount practical and technical interest, in addition to constituting a fundamental and formidable scientific challenge [1]. One of the most intriguing properties of water is its large numbers of crystalline polymorphs (currently 18, although no doubt still set to grow in the future), which are stable over a broad range of temperature and pressure [2]. Ice Ih is one of the most frequently-observed phases of water on Earth. However, different forms of ice can be found at medium and high pressures [1,3,4]; Tammann and Bridgman discovered first high-pressure form of ice in the early 20th century [5,6]. Indeed, in this new, 21st century alone, six ice polymorphs have been found (ices XIII-XVIII) so far [7][8][9][10][11]. However, it has also been observed that most water present in the Universe is in its amorphous form [12]. Much of the detailed and technical characterization and properties-determination work done on both the crystalline and amorphous forms of water, has been, and will continue to be, experimental. However, theoretical and/or molecular-simulation studies can help to explain more insights into the rich behavior of solid (both crystalline and amorphous) water phases, especially for high-pressure ices at the electronic-structure level, where empirical models can often be found wanting (due to the inevitable electron-cloud overlap at high pressures being well outside the electronic experience of near-ambient-condition empirical-model parameterization) [13].
The theoretical pursuit, and leverage, of crystal-structure-prediction (CSP) methods can be used to explore the plausible different crystalline polymorphs, and indeed, polytypes, of ice [14]. Indeed, CSP may be a fruitful avenue for assessing the structural and energetic viability, in a prototyping or screening fashion, for potential novel, yet-to-be-discovered, ice structures. In general, CSP, done well, can provide de-novo predictions for various crystalline polymorphs of a material from a first-principles approach (either with accurate empirical force-fields or from electronic-structure approaches). CSP can be deployed, inter alia, for the development of new pharmaceutical and pigment compounds and/or for the engineering of novel molecular materials [15][16][17][18]. CSP methods often predict many possible polymorphs within a small relative-energy-difference range (<10 kJ/mol) [19]. Thus, in order to gauge and rank polymorphs' relative energetic stability, it is essential to use high-level quantum-chemical methods to capture as accurately as possible the subtlety of intermolecular interactions, and to refine the (free-) energy/enthalpy ranking from (potentially) less accurate empirical force-fields or more crude electronic-structure representations [20]. Intermolecular interactions, which are relatively weak, are typically non-directional and long-range in nature; an accurate description thereof is a challenging task for modeling molecules and materials [21]. Density-Functional Theory (DFT), combined with Grimme D3-BJ methods (e.g., the D3 method with Becke-Jonson damping), and more sophisticated many-body dispersion (MBD) techniques, are generally found to be more suitable for a more accurate description of intermolecular forces, especially in higher-pressure crystalline polymorphs [22,23]. We also have considered the density-dependent dispersion correction (dDsC) developed by Steinmann et al. [24,25] for dispersion energy, and compared the results with D3-BJ and MBD results. dDsC is used to calculate the interactions between non-overlapping densities, which cannot be described accurately by standard Density-Functional Theory.
In the present work, we consider the low-energy polymorphs of water at different pressures (0, 4 and 8 kbar), predicted previously by us using a modified basin-hopping algorithm for crystal-structure prediction; here, we re-rank the ice polymorphs using sophisticated DFT and higher-level electronic-structure calculations [14]. Some of the polymorph structures are shown in Figure 1. The structure S12/4 1 6 20 8 10 , for example, has 12 (the number following the "S") molecules per unit cell. The numbers following the forward slash ("/") give the ring count: this structure has one four-membered ring, twenty-six-membered rings, and ten eight-membered rings. In our previous work [14], basin hopping was performed, and the ice polymorphs were ranked in order of enthalpic stability at these three pressures using the 4-site transferable interaction potential (TIP4P) water model, a popular empirical model for water developed by Jorgensen et al. [26]. In this article, we compare the results of energetic-stability ranking using both TIP4P water and electronic-structure calculations.
Crystals 2020, 10, x FOR PEER REVIEW 2 of 8 CSP can be deployed, inter alia, for the development of new pharmaceutical and pigment compounds and/or for the engineering of novel molecular materials [15][16][17][18]. CSP methods often predict many possible polymorphs within a small relative-energy-difference range (<10 kJ/mol) [19]. Thus, in order to gauge and rank polymorphs' relative energetic stability, it is essential to use highlevel quantum-chemical methods to capture as accurately as possible the subtlety of intermolecular interactions, and to refine the (free-) energy/enthalpy ranking from (potentially) less accurate empirical force-fields or more crude electronic-structure representations [20]. Intermolecular interactions, which are relatively weak, are typically non-directional and long-range in nature; an accurate description thereof is a challenging task for modeling molecules and materials [21]. Density-Functional Theory (DFT), combined with Grimme D3-BJ methods (e.g., the D3 method with Becke-Jonson damping), and more sophisticated many-body dispersion (MBD) techniques, are generally found to be more suitable for a more accurate description of intermolecular forces, especially in higher-pressure crystalline polymorphs [22,23]. We also have considered the density-dependent dispersion correction (dDsC) developed by Steinmann et al. [24,25] for dispersion energy, and compared the results with D3-BJ and MBD results. dDsC is used to calculate the interactions between non-overlapping densities, which cannot be described accurately by standard Density-Functional Theory.
In the present work, we consider the low-energy polymorphs of water at different pressures (0, 4 and 8 kbar), predicted previously by us using a modified basin-hopping algorithm for crystalstructure prediction; here, we re-rank the ice polymorphs using sophisticated DFT and higher-level electronic-structure calculations [14]. Some of the polymorph structures are shown in Figure 1. The structure S12/4 1 6 20 8 10 , for example, has 12 (the number following the "S") molecules per unit cell. The numbers following the forward slash ("/") give the ring count: this structure has one four-membered ring, twenty-six-membered rings, and ten eight-membered rings. In our previous work [14], basin hopping was performed, and the ice polymorphs were ranked in order of enthalpic stability at these three pressures using the 4-site transferable interaction potential (TIP4P) water model, a popular empirical model for water developed by Jorgensen et al. [26]. In this article, we compare the results of energetic-stability ranking using both TIP4P water and electronic-structure calculations.

Computational Details
We have considered our previously-reported unit-cell geometry and lattice parameters predicted using the basin-hopping algorithm with the TIP4P water force-field [27]. We re-optimized the geometry, including the relaxation of lattice parameters to the appropriate pressures (0, 4 and 8 kbar) using DFT, as implemented in the Vienna Ab-initio Simulation Package (VASP 5.4.4) [28][29][30][31]. Here, we have considered a plane-wave energy cutoff of 800 eV, within the framework of the Generalized Gradient Approximation (GGA), with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [32]. The Monkhorst-Pack method with a gamma-centered 6 × 6 × 6 k-point grid was adopted for integration in the Brillouin zone. Non-covalent interactions largely affects energy ranking of polymorphs, and thus, we have considered the more commonly-used method of D3-BJ developed by Grimme et al. [22], in conjunction with more sophisticated MBD methods developed by Tkatchenko et al. [23], and compared the calculated results. Looking to compare this high-quality DFT with more sophisticated, accurate and expensive electronic-structure methods, we also consider, in the case of S4/6 8 (ice Ic), local-MP2 electron correlation energy using the 6-21G(d) basis set as implemented in Crystal 09 [33,34] and CRYSCOR [35,36] software packages; here, computational tractability naturally limits us to smaller unit cells and system sizes.

Results and Discussion
We compare the relative TIP4P lattice enthalpies with those of PBE-D3 and PBE-MBD (Figure 2), at different pressures (0 bar, 4000 bar and 8000 bar); i.e., the enthalpies relative to the lowest-energy structures (which were not re-ranked in going from TIP4P to either DFT approach). At 0 kbar, both PBE-D3-and PBE-MBD-relative enthalpies are well correlated with TIP4P enthalpy, with very limited re-ranking occurring for some higher-enthalpy polymorphs (i.e., local-minima points on the crystal-state potential-energy surface, as opposed to global). The slopes of the PBE-D3 (PBE-MBD) relative energy with respect to their TIP4P counterparts are 1.573 (1.573), 2.724 (2.189) and 3.871 (3.531) at 0 bar, 4000 bar and 8000 bar, respectively (cf. Figure 2 and Figure S1 in the Supplementary Materials). PBE-dDsC, which calculates the interactions between non-overlapping densities, cannot describe accurately the stability-ranking of the ice polymorphs (cf. Figure S2, Supplementary Materials).
Thus, the level of disparity, and re-ranking, increases with the application of greater external pressures. A similar deviation is also observed for the PBE-D3 and PBE-MBD cases. As mentioned previously, this would not necessarily be unexpected, given that electron-cloud overlap at high pressures is certainly not a feature of electronic-level behavior for the "comfort zone" of empirical-model parameterization at near-ambient temperatures and pressures, and typically in the liquid state also, as is the case for TIP4P [26].
The calculated relative enthalpy using the TIP4P water model and DFT (via both PBE-D3 and PBE-MBD) are shown in Tables 1-3 (see also Tables S1-S3) for different externally-applied pressures (0, 4000 and 8000 bar) along with their densities (Tables S4-S6). At low (zero) pressure, turning to Table 1, we find that ice Ic, ice Ih and ice III are the most stable structures (structures S4/6 8 , S8/6 16 and S12/5 8 7 8 8 8 , respectively) predicted using TIP4P, with lattice energy (per water molecule) of −57.104, −57.104 and −56.793 kJ/mol at zero pressure; thus, ices Ic and Ih have identical TIP4P lattice energies (−57.104 kJ/mol), as predicted using the TIP4P water model. This finding is contrary to the fact that hexagonal ice Ih is the only polytype found widely on Earth, with the exception only of a small amount of ice Ic occasionally present in the upper atmosphere [1]. The predicted energy considering local-MP2 electron correlation for S4/6 8 (ice Ic) was found to be −56.90 kJ/mol. Further, our DFT results suggest that ice Ih is energetically more stable (~0.1 kJ/mol) compared to ice Ic. We also note that the energy of ice III is comparatively higher (~3kJ/mol) via DFT approaches than ice Ih and ice Ic, which suggests that ice Ih is energetically preferred at low pressure, chiming well with experiment. and PBE-MBD) relative energy with respect to TIP4P relative energy are specified at 0 bar, 4000 bar and 8000 bar of external pressure, respectively. For visual comparison, the dash (---) line has a slope of unity, to represent perfect numerical TIP4P-DFT relative-energy agreement; the growing disparity with mounting pressure is clearly evident.
Although the PBE-D3-and PBE-MBD-predicted results are close to each other, the PBE-dDsCpredicted results are quite different. For instance, PBE-dDsC predicts, wrongly, ice Ic to be most stable structure for zero external pressure (cf. Table 1).
It is interesting to note that the H-O-H angle for ice polymorphs predicted using the TIP4P force-field are comparatively smaller (~104°). However, the H-O-H angle for ice polymorphs optimized using our DFT calculations (using both PBE-D3 and PBE-MBD) are ~107°, as observed experimentally [37]. With increasing external pressure, the average H-O-H angle decreases. Ice III is the most stable polymorph at 4000 bar and 8000 bar, and the average H-O-H angle for Ice III was ~107°. Also, as expected, the density of ice polymorphs with pressure increases (cf . Tables S4-S6).
Experimentally, ice III, which is a tetragonal-structure crystalline ice, is formed at 3000 bar (300 MPa) by lowering its temperature to 250 K [37,38]. Therefore, it is hardly surprising that we found ice III (i.e., structure S12/5 8 7 8 8 8 ) to be the most stable ice in under 4 and 8 kbar pressurization (see Tables 2 and 3). We also note that ice II and ice VI are observed at 4 kbar pressure, whereas ice VII is observed at the higher pressure of 8 kbar, which is in agreement with experimental expectations [2].    Although the PBE-D3-and PBE-MBD-predicted results are close to each other, the PBE-dDsC-predicted results are quite different. For instance, PBE-dDsC predicts, wrongly, ice Ic to be most stable structure for zero external pressure (cf. Table 1).
It is interesting to note that the H-O-H angle for ice polymorphs predicted using the TIP4P force-field are comparatively smaller (~104 • ). However, the H-O-H angle for ice polymorphs optimized using our DFT calculations (using both PBE-D3 and PBE-MBD) are~107 • , as observed experimentally [37]. With increasing external pressure, the average H-O-H angle decreases. Ice III is the most stable polymorph at 4000 bar and 8000 bar, and the average H-O-H angle for Ice III was 107 • . Also, as expected, the density of ice polymorphs with pressure increases (cf . Tables S4-S6).
Experimentally, ice III, which is a tetragonal-structure crystalline ice, is formed at 3000 bar (300 MPa) by lowering its temperature to 250 K [37,38]. Therefore, it is hardly surprising that we found ice III (i.e., structure S12/5 8 7 8 8 8 ) to be the most stable ice in under 4 and 8 kbar pressurization (see Tables 2 and 3). We also note that ice II and ice VI are observed at 4 kbar pressure, whereas ice VII is observed at the higher pressure of 8 kbar, which is in agreement with experimental expectations [2].

Conclusions
We (re)-ranked the previously-reported ice polymorphs [14] at 0, 4 and 8 kbar, predicted by us via TIP4P using a modified basin-hopping algorithm for crystal-structure prediction, in order of energetic stability using high-level quantum-chemical calculations (primarily with sophisticated DFT approaches, but also periodic local-MP2 in the case of ice Ic). Comparing against TIP4P relative energies (with no DFT re-ranking of global minima at each pressure studied), under low-(zero-) pressure conditions, both PBE-D3 and PBE-MBD predicted ice-polymorph (relative-) energy rankings consistent with TIP4P-predicted ones. The slopes of PBE-D3 (PBE-MBD) relative energy vis-à-vis their TIP4P counterparts 1.573 (1.573), 2.189 (2.189) and 3.871 (3.531) under 0, 4 and 8 kbar pressurization, respectively. We find that, although the TIP4P water model correctly predicts the ranking of ice polymorphs (especially at lower pressure, with limited higher-energy re-ranking at 4 and 8 kbar), the relative-energy deviation compared to high-level quantum-mechanical calculations increases with mounting applied external pressure.
Therefore, this study shows that high-level quantum-mechanical calculations are important for stability ranking of ice polymorphs, especially under high-pressure conditions, or if polymorphs' relative energies need to be quantified accurately. We do not consider the kinetic accessibility in the present work of the different polymorphs, or their (experimental) realization; in part, that depends upon temperature-pressure-time kinetic pathways and passage through intermediate structures between these polymorphs as a series of thermodynamic (local) minima. In addition, as [12,13] discuss more specifically for DFT for water, as well as [39] in more detail for empirical potential models and electronic-structure approaches for ice (and [40,41] for water-gas systems), increasing pressurization leads to qualitative failings of empirical-potential water models for its condensed-matter states, owing to the onset electronic-cloud overlap.
The use of advanced polarizable water models may remedy this deficit of empirical models in future for lower-pressure applications and ice-polymorph scanning [42], but accurate electronic-structure approaches are sine qua non at higher pressure if quantitative accuracy is desired [43,44].

Author Contributions:
The calculations and analysis were done by P.K.S. and C.J.B. supervised by N.J.E. The manuscript was written by P.K.S. and N.J.E. All authors have read and agreed to the published version of the manuscript.