Ice Ih anomalies: Thermal contraction, anomalous volume isotope effect, and pressure-induced amorphization

Hexagonal, proton-disordered ice Ih,^1,2 the predominant solid phase of water on Earth, displays several anomalous thermodynamic properties. It thermally contracts at low temperatures (<70 K), unlike most solids, which expand upon heating. Substitution of hydrogen by deuterium in ice Ih results in a volume expansion by 0.09% at 10 K,^3,4 which is also the opposite behavior from most other solids. Furthermore, this anomalous volume isotope effect (VIE) grows with temperature and persists even in the liquid phase, despite the fact that isotope-dependence of volume is a quantum effect and is thus expected to vanish in the high-temperature (classical) limit.

Both temperature- and isotope-dependence of volume are quantum as well as anharmonic effects. Anharmonicity usually causes the frequencies of lattice vibrations to decrease upon expansion, since longer atom-atom distances tend to experience a softer part of the anharmonic potential. It, in turn, makes more vibrational states thermally accessible to the solid for an entropy gain, energetically favoring expansion. The fact that the opposite is the case with ice Ih means that many of its thermally accessible lattice vibrations have reverse anharmonicity, where frequency increases upon expansion^5,6 (i.e., their mode Grüneisen parameters are negative). Specifically, the modes displaying such characteristics are the hydrogen-bond bending modes (transverse acoustic phonons)^7,8 that squeeze hollow hydrogen-bond cages of ice Ih into denser structures.

With pressure loading, one expects the frequencies of these modes to decrease further and ultimately reach zero, causing a structural instability in compressed ice Ih. This is believed to be the process underlying the “mechanical melting”⁹ (to be distinguished from thermodynamic melting) of ice Ih into a noncrystalline, glass-like, denser phase known as high-density amorphous (HDA) ice^5,9–11 at 1.0 GPa and 77 K. Both the negative thermal expansion and this pressure-induced amorphization therefore share a common lattice-dynamical origin: the reverse anharmonicity of quantized hydrogen-bond bending modes, which strongly modulate the volume of hollow cages in ice Ih. Corroborating this is the observation that other tetrahedrally coordinated hollow crystals such as silica¹² also exhibit thermal contraction. On the other hand, which modes are responsible for the anomalous VIE is a far more complex issue.

These anomalous behaviors of ice Ih have striking parallels with those of liquid water, underscoring the broader significance of the issue beyond just ice. Liquid water also thermally contracts at lower temperatures (<4 °C) before it expands. Deuterated liquid water has a greater molar volume than normal liquid water. Compression of ice Ih between 250 and 273 K causes melting to the denser liquid phase. Indeed, the Ih-HDA transition pressure and temperature can be roughly extrapolated from the melting curve of ice Ih, which has a negative slope reflecting the volume contraction of water upon melting. These comparisons suggest that the anomalous behaviors of liquid water, including the negative slope of the melting curve, may have a similar (if not the same) origin as the ice Ih anomalies,¹³ the understanding of which sheds light on the former. On this note, it is worth mentioning that the existence of two distinct phases of liquid water, the low- and high-density liquids, have been proposed¹⁴ as a continuation of the low- and high-density amorphous (LDA and HDA) phases,^9,15 respectively, though this remains to be a highly controversial issue.^16–18 

The ice anomalies are macroscopic manifestations of the quantum effects of nuclei moving in anharmonic potentials. These potentials are, in turn, created by intra- or inter-molecular covalent, hydrogen-bond, and dispersion interactions in the proton-disordered or amorphous solid environment. Therefore, quantitative or even just the correct descriptions of these effects are expected to require high-accuracy quantum many-body methods for both electrons and nuclei.^19–24 Nevertheless, Ramírez et al.²⁵ were able to reproduce the temperature-dependence of the volume of ice Ih accurately for a wide range of temperatures using the empirical force field, q-TIP4P/F. The nuclear quantum effects were taken into account with either path-integral molecular dynamics (PIMD) or the quasiharmonic approximation (QHA), which were shown to differ from each other to an insignificant extent. Strässle et al.⁵ experimentally demonstrated the pressure-induced softening of transverse acoustic phonons in ice Ih and, combining the experimental data with a simple lattice-dynamical model, showed that it was responsible for the negative thermal expansion and pressure-induced amorphization at 2.5 GPa.

Pamuk et al.⁸ used density-functional theory (DFT) with and without dispersion correction in the QHA and determined the primary source of the anomalous VIE as reverse anharmonicity of the intramolecular O–H/O–D stretching modes, though they did not offer an explanation of its growth with temperature. Remarkably, on the basis of the calculations, Pamuk et al.⁸ predicted the normality of the VIE of ice Ih upon ¹⁶O/¹⁸O substitution, which was borne out by experiment. On the other hand, when Murray and Galli²⁶ used another dispersion-corrected DFT, the anomalous VIE of H/D substitution was not reproduced.

These prior studies suggest that the negative thermal expansion and pressure-induced amorphization of ice Ih are relatively well understood and robustly reproducible computationally. However, the anomalous VIE is a much subtler effect and its origin is still not fully understood. Its increase with temperature has not been explained, either.

In this article, therefore, we present embedded-fragment ab initio second-order many-body perturbation (MP2) calculations of the structures of ice Ih as a function of temperature and pressure, including the quantum effects of anharmonic lattice vibrations in the QHA. We show that they too reliably reproduce the thermal contraction of ice Ih at low temperatures and thermal expansion at higher temperatures, though the latter is underestimated. Our MP2 calculations confirm that the contraction is caused by reverse anharmonicity (negative mode Grüneisen parameters) of the hydrogen-bond bending modes, whereas the subsequent expansion is due to the hydrogen-bond stretching modes.

Our MP2 calculations show that the net VIE at 0 K is the result of an extremely close competition between volume-expanding and contracting effects involving all fundamental lattice vibrations (not just the O–H/O–D stretching modes) that exist in the zero-point state. Consequently, even MP2 gives contradicting predictions for the VIE, depending on minor calculation details (specifically, the choice of the embedding field). However, the thermal growth of the anomalous VIE is reliably reproduced and primarily ascribed to the isotope-dependent volume-expanding effect of the librational modes.

The calculations also detect near-abrupt softening of the hydrogen-bond bending (acoustic) modes at 2.35–3 GPa, accompanied by difficulty in converging the geometry optimization. This signals a mechanical instability in ice Ih under pressure, leading to the pressure-induced amorphization, which has been observed experimentally.⁹ This is further verified by Born–Oppenheimer molecular dynamics (BOMD) using on-the-fly MP2 forces; at 2.5 GPa and 75 K, the structure of ice degrades from crystalline Ih to amorphous with an initial 20% reduction in volume in a short time.

The internal energy, equilibrium structure, and phonon dispersion and density of states (DOS) of an infinite crystal of ice Ih were determined with an implementation of the embedded-fragment scheme^27–29 known as the binary-interaction method (BIM).^29–32 It was combined with MP2 with the aug-cc-pVDZ or aug-cc-pVTZ basis set. The energy was calculated also at the coupled-cluster singles and doubles (CCSD) level with the aug-cc-pVDZ or aug-cc-pVTZ basis set. Here, the energy per unit cell of a crystal was obtained as a sum of MP2/CCSD energies of water dimers and monomers embedded in the self-consistently determined electrostatic field of the crystal. See the supplementary material of Ref. 33 for the formalisms and algorithmic details of BIM. Two types of embedding field were considered: atomic partial charges³⁴ calculated by the electrostatic potential (ESP) module³⁵ of NWChem³⁶ and molecular dipole moments,³⁷ both at the Hartree–Fock (HF) level with the corresponding basis set. These embedding fields are hereafter referred to as “atomic” and “dipolar.” First and second derivatives of the energy with respect to atomic coordinates and lattice constants were also obtained efficiently as sums of the corresponding derivatives of the fragment energies, which were then utilized to determine the equilibrium crystal structures and phonon dispersion curves and DOS. The electronic-structure backend for the BIM calculations was a version of NWChem³⁶ modified for calculating embedding-field derivative contributions to analytical gradients.

Orthorhombic unit cells of ice Ih containing 32 or 64 water molecules were generated by a Monte Carlo method.³⁸ They had quasi-random proton configurations obeying the ice rules with net zero unit-cell dipole moment. Various lattice-sum truncation radii of the BIM calculations in the notation of Ref. 33 were as follows: R_QM = 8 Å, R_EF = 10 Å, and R_LR = 200 Å. At a typical density, each monomer interacted quantum-mechanically with its 66 nearest-neighbor monomers, going beyond the third coordination shell, at which the MP2 correlation interaction decayed to a negligible level.¹⁹ Each of these monomers was embedded in the field of about 128 nearest neighbors. Each dimer was embedded in the union of the fields of the two monomers. Classical electrostatic interactions beyond the embedding field with an additional 1 140 000 monomer-monomer pairs were also included in the internal energy.

The crystal structure was optimized at 0 GPa and 0 K with tolerances for the atomic and lattice gradients being 6 × 10⁻⁴ and 3 × 10⁻⁴ a.u., respectively. The unconstrained optimization introduces an imperceptibly small positional disorder of the oxygen atoms, which is a physical outcome of the proton disorder neglected in most structural refinements.³⁹ Furthermore, as discussed below, a more visible loss of symmetry and difficulty converging the geometry optimization occurs at higher pressures, owing to the pressure-induced amorphization.

The interaction force constants were computed in the 3 × 3 × 3 unit cells. The normal modes and harmonic frequencies of both H₂O and D₂O ices were calculated at 9 × 9 × 9 evenly spaced k points in the reciprocal unit cell. A regularization of acoustic phonon branches at the Γ point was invoked.³¹ 

The thermal expansion and VIE of H₂O and D₂O ice Ih were determined from their vibrationally corrected volumes (V₀) as a function of pressure (P) and temperature (T). The vibrationally corrected volume V₀ is the one that minimizes the Gibbs energy evaluated in the QHA.^7,8,40–42 In atomic units, the Gibbs energy per unit cell is written as

with the Helmholtz energy being

where E(V) is the internal (i.e., electronic) energy per unit cell as a function of volume (V), K is the number of k points in the reciprocal unit cell, ω_nk(V) is the harmonic frequency of the nth (fundamental) phonon branch with wave vector k, and β = (k_BT)⁻¹. The residual entropy due to proton disorder is negligible for our problems and was therefore neglected.

The Helmholtz energy F(V, T) was evaluated at 0 to 250 K with a 1-K interval at 29 evenly spaced volumes about the equilibrium structure at 0 GPa and 0 K. In the MP2/aug-cc-pVDZ calculations with the atomic-embedding field, these volumes ranged from 27.3 to 32.0 ų per molecule, corresponding to electronic (virial) pressures from 1.25 to −1.4 GPa, respectively. The force constants were evaluated and the harmonic frequencies were calculated at each volume, wherein no approximate linearization of frequencies⁸ was invoked. The volumes were obtained by first isotropically scaling the lattice constants and atomic coordinates and then relaxing them under a fixed pressure for each volume. This relaxation step caused slight changes in the unit-cell shape, except when an abrupt and dramatic change occurred near the limit of mechanical stability leading to amorphization (see below). The structures near the stability limit were not included in the QHA calculations. The vibrationally corrected volume (V₀) at finite pressure (including both virial and phonon pressures) and temperature was then obtained by solving

using a third-order polynomial fit of F(V, T). V₀(P, T) thus obtained directly provides both thermal expansion and VIE.

The QHA relies on the analytical expression of the partition function of a harmonic insulating crystal, and the effect of anharmonicity, responsible for thermal expansion and VIE, is implicit in the parametric volume-dependence of the harmonic frequencies, ω_nk(V). It is qualitatively superior to Grüneisen’s theory,⁴³ which expresses the isochoric thermal expansion coefficient (α_V) in terms of the total Grüneisen parameter (γ), isothermal bulk modulus (B₀), and isochoric heat capacity (C_V) as

Here, γ is related to the mode Grüneisen parameters by Eq. (2) of Li et al.⁴³ The latter are defined by

which is positive if the frequency increases upon compression and negative if the frequency decreases upon compression, corresponding to reverse anharmonicity.

Grüneisen’s theory neglects phonon pressure and also assumes constancy of γ_nk with volume. It therefore cannot account for VIE and may give poorer results for thermal expansion.^42,43 Nonetheless, the mode Grüneisen parameters defined above remain useful for understanding the origin of these effects semi-quantitatively. Owing to the congestion of phonon dispersion curves for a large unit cell, a tiny fraction of γ_nk are not calculated correctly by a finite-difference method because of reordering of phonon dispersion branches between two volumes. Conclusions drawn from the analyses using them are unaffected. The QHA calculations are not affected at all by reordering of phonon dispersion branches.

At a pressure beyond the limit of mechanical stability, geometry optimization of ice Ih was found to become increasingly difficult, signaling the pressure-induced amorphization. To confirm this, a short (320 fs), constant-pressure (2.5 GPa), constant-temperature (75 K) BOMD simulation was carried out starting from the 64-molecule-unit-cell crystalline structure optimized at 2.35 GPa, using on-the-fly MP2/aug-cc-pVDZ forces with the Nosé–Hoover thermostat^44,45 and Berendsen pressure coupling algorithm.⁴⁶ A time step of 1 fs was used. Additional algorithmic details of MP2 BOMD can be found in Ref. 33.

First, we assess the performance of MP2 for the structural prediction of ice Ih.^19–22,28 Table I compares the structural parameters obtained with MP2 using different basis sets and embedding fields against the experimental data.^2,4 Differences between 32- and 64-molecule unit cells are minor, indicating that the 32-molecule unit cell is large enough to represent an average proton-disordered structure. From here on, all the aug-cc-pVDZ calculations are based on the 64-molecule unit cell, while the aug-cc-pVTZ results are obtained with the 32-molecule unit cell.

In all cases, MP2 underestimates the average oxygen-oxygen distance 〈R_OO〉, lattice constant a, and the volume V_e, underscoring the well-known tendency of the theory to overbind molecular crystals.^{20,23,31,32,43,47} The dipolar embedding suffers less from the overbinding, displaying a systematic improvement of the calculated quantities with a basis-set extension.

Table II compares the calculated and observed⁴⁸ binding energies of ice Ih. Consistent with the structural data, MP2 overbinds, but not significantly more than various approximations of DFT.^8,26 Both the equilibrium and vibrationally corrected binding energies of CCSD/aug-cc-pVTZ with the dipolar-embedding field are in accurate agreement with the observed (within 4%), while the same method using the atomic-embedding field overshoots the latter significantly. We thus conclude that the atomic-embedding field tends to overestimate the water-water interaction strength in ice.

The vibrationally corrected volumes of H₂O and D₂O ice Ih are plotted as a function of temperature in Fig. 1. The equilibrium and vibrationally corrected volumes, VIE, and bulk moduli are compiled in Table III. The figure shows that the MP2/aug-cc-pVDZ calculation with atomic-embedding field correctly reproduces the thermal contraction at low temperatures and thermal expansion at higher temperatures occurring in both H₂O and D₂O ices, though the expansion at higher temperatures is underestimated considerably.

To analyze their origin, we first categorize in Fig. 2 the phonons into five groups with different volume-dependence: the hydrogen-bond bending (HB), hydrogen-bond stretching (HS), librational (L), intramolecular HOH/DOD bending (B), and intramolecular O–H/O–D stretching (S) modes. The hydrogen-bond stretching (HS) modes increase their frequencies upon compression and thus have positive mode Grüneisen parameters. This is because the water molecules have greater difficulty pushing one another away in a more crowded environment. In contrast, the hydrogen-bond bending (HB) modes (transverse acoustic modes) act to alleviate the congestion in a hollow hydrogen-bonded cage of ice Ih, decreasing the frequencies upon compression and thus having negative mode Grüneisen parameters. Since the strengths of an O–H/O–D bond and of its hydrogen bond have inverse correlation, the O–H/O–D stretching (S) frequencies decrease with pressure and have negative mode Grüneisen parameters. The librational (L) modes increase their frequencies upon compression, while the bending-mode frequencies seem nearly stationary with pressure.

The volume-contracting effect of the hydrogen-bond bending (HB) modes has been implicated^5,7,51 in the negative thermal expansion. Our analysis confirms this. Figure 3 plots the phonon-group contributions to C_V(T) γ(T), a factor in the thermal expansion coefficient according to Grüneisen’s theory [Eq. (4)]. Although our calculations are based on the more accurate QHA than Grüneisen’s theory, the total C_V(T) γ(T) function (black curve) semi-quantitatively reproduces the temperature-dependence of volume in Fig. 1. We find that the thermal contraction is entirely due to the hydrogen-bond bending (HB) modes, whereas the subsequent thermal expansion is mostly caused by the hydrogen-bond stretching (HS) modes and to a much lesser extent by the librational (L) modes.

The hydrogen-bond bending (HB) and the hydrogen-bond stretching (HS) modes make opposite temperature-dependent contributions to the volume, while their DOS adjoin each other. These two groups are, however, clearly distinguishable by the mode-stretching parameter⁷ defined as

where the average is taken over all modes labeled by k with frequencies in the window of ω to ω + Δω, $uik$ is the displacement of the oxygen atom in the ith molecule and the kth normal mode, and n_i,j is the unit vector from the jth to ith molecule. The i and j summation runs over all hydrogen-bonded pairs in the unit cell.

Figure 4 shows that the mode-stretching parameter sharply rises from near zero to about 0.5 at 125 cm⁻¹, which demarcates the lower-frequency hydrogen-bond bending (HB) modes and the higher-frequency hydrogen-bond stretching (HS) modes. This demarcation accurately coincides with the sign change in the mode Grüneisen parameter. In the hydrogen-bond bending (HB) region, its values range from −0.5 to −3 with an average of −1.2, whereas in the hydrogen-bond stretching (HS) region, its values fall between +1 and +2.

The thermal expansion of ice Ih at higher temperatures is underestimated by a factor of 3 (see Fig. 1), causing also the temperature of maximum density to be shifted to a somewhat higher value than the observed (ca. 70 K). It is unlikely that the underestimation is caused by higher-order anharmonicity neglected in the QHA, since Ramírez et al.⁵² showed that the QHA results for the temperature-dependence of volume of ice Ih agree accurately with the results of PIMD, which is exact in the limit of an infinite number of beads. Furthermore, the QHA results seem to slightly overestimate (instead of underestimate) the thermal expansion as compared with the PIMD results. Interestingly, the dispersion-corrected DFT calculation in the QHA by Pamuk et al.⁸ also displayed a similar underestimation of thermal expansion.

A more probable source of this discrepancy is the overestimation of the bulk modulus (B₀), as evidenced in Table III. MP2 systematically overbinds ice Ih (see Sec. III A) and molecular crystals in general, and thus underestimates their volumes.^{20,23,31,32,43,47} The error in the volumes is amplified in B₀ because the latter involves the derivative of pressure with respect to volume.⁴³ Since this quantity enters Grüneisen’s expression [Eq. (4)] of the thermal expansion coefficient (α_V) as a denominator (implicitly in the case of the QHA), its overestimation results in the underestimation of α_V. In other words, MP2 tends to predict too stiff a molecular crystal, which therefore cannot expand thermally as much as it should. In our previous study⁴³ of solid CO₂-I, the underestimation of α_V by a factor of nearly two by MP2 was ascribed to the overestimation of the bulk modulus by ca. 100%. This error was subsequently erased by Heit et al.,⁴² who used more accurate volumes calculated by CCSD with noniterative triples [CCSD(T)] in the complete-basis-set (CBS) limit in conjunction with the QHA.

Figure 5 plots the MP2 equation of state of H₂O ice Ih, which illustrates this point. The equilibrium bulk modulus (B_e) is proportional to the reciprocal of the slope of the red curve at 0 GPa, which plots the equilibrium (non-vibrationally corrected) volume. The reciprocal slope is greater than the corresponding reciprocal slope of the experimental curve (black curve). Clearly, this is caused by the overbinding of ice Ih by MP2/aug-cc-pVDZ. The inclusion of anharmonic vibrational corrections at 150 K in the QHA (blue curve) softens the crystal and brings the value of B_e from 19.4 GPa down to B_H2O = 14.8 GPa, which is much closer to the experimental value of 9.8 GPa at 145 K,⁵⁰ but is still too large. This is largely responsible for the underestimation of the thermal expansion around that temperature. In other words, this discrepancy is primarily traced to the smallness of the basis set.

Next, we seek to explain the anomalous VIE and its temperature dependence. Figure 1 and Table III show that the MP2/aug-cc-pVDZ calculation with the atomic-embedding field reproduces the anomalous VIE at 0 K and also its subsequent increase with temperature, both nearly quantitatively. Experimentally,^3,4 the VIE at 0 K is −0.09% (the negative value means anomalous behavior), which remains nearly the same at −0.08% at 100 K, but increases in magnitude to −0.21% at 220 K. Computationally, we find the corresponding values to be −0.14% (0 K), −0.16% (100 K), and −0.25% (220 K), which are in more accurate agreement than other previous calculations.^8,26

However, as will be discussed in the  Appendix, the calculated results change qualitatively once such a minor detail as the choice of the embedding field is altered; the calculations with the dipolar-embedding field do not reproduce the anomalous (negative) VIE. This means that the agreement in the sign of the VIE (if not its temperature-dependence) from the atomic-embedding field is likely accidental. Computationally reproducing it for the right reason remains to be an exceedingly difficult task even for ab initio MP2. This difficulty can be understood by analyzing the origin of the VIE.

At 0 K, the Gibbs energy reduces to

At V = V_e, the enthalpy is at a minimum, and so we can write

where c₁ and c₂ > 0 are constants. In this vicinity, ω_nk(V) can be approximated⁸ as

where γ_nk is the mode Grüneisen parameter. The vibrationally corrected volume V₀ that minimizes the Gibbs energy then satisfies⁸ 

Table IV lists phonon-group contributions to the anharmonic vibrational correction to the volume of H₂O and D₂O ices as estimated by 〈ωγ〉. Its total VIE does not agree with the value from the QHA, but it reproduces its sign and thus the anomalous VIE. This table paints an exceedingly complex picture of this anomaly. No single phonon group determines even the sign of the VIE. This is because this effect exists even at 0 K, at which it is caused by anharmonicity of zero-point vibrational levels involving all fundamental modes. This is in striking contrast to thermal expansion, in which only low-frequency thermally accessible lattice vibrations play a dominant role and they tend to share the same volume-expanding or contracting effect. Consequently, none of the phonon groups is negligible in VIE. Even the bending (B) modes, whose contribution to 〈ωγ〉 is the smallest, has a VIE of a similar order of magnitude as the total VIE. It is true that ice Ih exhibits this anomaly owing to the large volume-contracting effect of the O–H/O–D stretching (S) modes, as pointed out by Pamuk et al.,⁸ but whether the VIE is negative or merely small positive is determined by close competition of volume-expanding and contracting effects of all phonon groups. It is doubtful if any method short of a nearly exact one, such as CCSD(T) in the CBS limit, can reproduce it for the right reason.

Our MP2 calculation, however, reproduces the thermal increase of the anomalous VIE. This can be unequivocally attributed to the volume-expanding effect of the librational (L) modes, which are the lowest-energy phonons whose frequencies are strongly isotope-dependent (see Fig. 3). Unlike the VIE, this property is not the result of cancellation of competing effects and is thus more reliably reproduced by theory.

In 1984, Mishima et al.⁹ extrapolated the negative-slope melting curve of ice Ih to low temperatures, predicting and then observing the pressure-induced transition of ice Ih to HDA at 1.0 GPa and 77 K. Combining inelastic neutron scattering from ice Ih under pressure and a lattice-dynamical model calculation, Strässle et al.^5,6 argued that the Γ–M transverse acoustic phonon dispersion branch softens with pressure, ultimately becoming flat zero at ca. 2.5 GPa, whereupon a mechanical instability exists, causing the pressure-induced amorphization. These phonons are believed to be the same as what we call the hydrogen-bond bending (HB) modes with the volume-contracting effect. Hence, both the negative thermal expansion and pressure-induced amorphization originate from this particular phonon group.

The phase diagram of H₂O is drawn in Fig. 6. Simulations with empirical potentials⁵³ and inelastic neutron scattering measurements⁶ have established two amorphization regimes: mechanical melting at low temperatures and thermal melting at high temperatures. The thermal melting of ice Ih in the range of 160–250 K results in the metastable formation of supercooled water (with subsequent nucleation of other crystalline phases) at transition pressures that extrapolate smoothly from the melting curve of ice Ih.¹⁰ Below 160 K, the slope of the transition pressure curve changes abruptly, signaling the entry into the regime of pressure-induced amorphization by mechanical melting. The latter (as opposed to thermal melting) occurs at low temperatures because thermal energy is insufficient for ice to reach structures with lower Gibbs energies. In other words, with increasing pressure, the Gibbs energy of HDA becomes lower than that of ice Ih, where the latter is thus metastable and is unable to cross an energy barrier for thermal melting. With further pressure loading, flattening of transverse acoustic phonons occurs and Born’s stability condition is violated,⁵³ and ice Ih thus undergoes mechanical melting to HDA.

Our MP2 calculations detect a sudden reduction in volume, a loss of symmetry⁵⁴ and long-range order in structure, and softening of the acoustic phonons, all in the pressure range of 2.35–3 GPa, strongly suggesting the occurrence of pressure-induced amorphization during the computer simulations.

Figure 7 plots the pressure-dependence of the structural parameters of H₂O ice Ih calculated by MP2/aug-cc-pVDZ with the atomic-embedding field. With increasing pressure, the geometry optimization starting from an initial crystalline structure became more difficult at >2.35 GPa. By 3 GPa, it did not converge within an allotted number of iterative cycles. The nonconverged geometry at 3 GPa shows a sudden drop in volume by 15% relative to the ambient-pressure value and also a loss of symmetry, as evidenced by an elevated root-mean-square-deviation of 0.04 in the oxygen fractional coordinates from the experiment values (it is nearly constant at 0.01 below 2.5 GPa). The volume reduction of 15% is in good agreement with the corresponding value of 15% at 2.5 GPa predicted by the lattice-dynamical model based on the inelastic neutron scattering data.⁶ Our calculated mechanical transition pressure (2.35–3 GPa) is also roughly consistent with the observed value of 1.0 GPa by Mishima et al.,⁹ 1.4 GPa predicted by the q-TIP4P/F simulation,^52,55 or 2.5 GPa quoted above.⁶ 

Figure 8 draws the crystalline and disordered structures at 0 and 3 GPa, respectively. While the latter is not converged, it shows that molecules permeate through the otherwise hollow spaces of the hydrogen-bonded crystalline cage and the solid becomes denser. It is an initial structural change towards complete amorphization.

Superimposed in Fig. 6 is the calculated thermodynamic phase boundary (blue line) between Ih and HDA. It was calculated as the temperature and pressure at which the Gibbs energies (minus the vibrational contributions) of the two phases, i.e., E(V) + PV − TS_conf., agree with each other at the MP2/aug-cc-pVDZ level with the atomic-embedding field. Here, S_conf. = 0 (configuration entropy) was assumed for ice Ih and S_conf. = 2 J K⁻¹ mol⁻¹ for HDA,⁵⁶ and E(V) and V for HDA corresponded to the snapshot structure in Fig. 8. The calculated thermodynamic transition pressure (blue curve) is about 1.6 GPa at 0 K with little temperature-dependence, above which ice Ih is metastable. This transition pressure is a crude estimate, as we neglected vibrational contributions and used a single nonconverged HDA structure, whose energy and volume were not ensemble-averaged. Nevertheless, it is consistent with the extrapolated melting curve of ice Ih (dotted-dashed curve) at 0.6 GPa with little temperature-dependence. The calculated mechanical transition pressure is no higher than 2.35 GPa at 0 K, above which ice Ih is no longer metastable but becomes unstable and transforms into HDA. This is also consistent with the experimental transition pressures (dots) in the mechanical melting regime (<160 K).

To further confirm that the failure of geometry optimization has a physical origin and is not a computational artifact, we performed a BOMD simulation using on-the-fly MP2/aug-cc-pVDZ forces. Figure 9 attests to a rapid reduction in volume by nearly 20% once the initial crystalline structure is allowed to relax according to classical equation of motion. The lattice constant ratio (c/a) also fluctuates, unequivocally demonstrating that ice Ih is unstable at 2.5 GPa. With the Berendsen pressure coupling algorithm⁴⁶ employed here, the relaxation time of unit cell volume does not have much physical significance, but the dynamics and ensemble averages are meaningful.

Concomitant with the detection of the structural instability, our MP2 calculations also observe pressure-induced red-shifting of DOS due to acoustic phonon branches at ca. 2 GPa, as shown in Figure 10. A clearer piece of evidence is given in Fig. 11, where rapid softening of the frequency of a representative acoustic mode near the edge of the Brillouin zone is recorded. The softening becomes complete at 2.35 GPa, coinciding with the onset of the geometry optimization difficulty. This pressure is an upper bound for the mechanical instability predicted by MP2/aug-cc-pVDZ. The lag in the transition pressure (2.35–3 GPa) suggested by the geometry optimization procedure is likely due to an insensitivity of the structure to a small number of zero or imaginary frequencies. It is, however, interesting to note that, experimentally also,⁹ it takes an additional 0.5 GPa of pressure loading for the amorphization to complete after it starts at 1.0 GPa.

We have applied the embedded-fragment ab initio MP2 method to anomalous thermodynamic properties of proton-disordered ice Ih. In combination with the QHA, we have computed the Gibbs energies of ice Ih across a range of pressures and temperatures. The calculations have taken into account the quantum effect of nuclei moving in the anharmonic potentials created by the covalent, hydrogen-bond, and dispersion interactions within and between the water molecules, both from first principles without any empiricism or nonsystematic approximations to the electronic or vibrational Hamiltonian. This has enabled a prediction of the vibrationally corrected volume of ice Ih as a function of temperature and isotope, reproducing quantitatively the thermal contraction at low temperatures and qualitatively the thermal expansion at higher temperatures. An analysis using the mode Grüneisen parameters has confirmed that the negative thermal expansion is almost exclusively ascribed to the hydrogen-bond bending (HB) modes (transverse acoustic phonons), while the subsequent expansion is mostly due to the hydrogen-bond stretching (HS) modes and to a lesser extent to the librational (L) modes. The underestimation of the rate of thermal expansion at higher temperatures is traced largely to a general overestimation of the bulk modulus by MP2/aug-cc-pVDZ, although the same difficulty has been encountered in dispersion-corrected DFT.⁸ 

The anomalous VIE is exceedingly difficult to predict correctly, owing to its small magnitude determined by subtle cancellation among many opposing phonon-group contributions. The atomic-embedding field succeeds in reproducing its correct sign, but the dipolar-embedding field does not, highlighting the sensitivity of the result to minor calculation details. This effect was ascribed largely to the negative mode Grüneisen parameters of the O–H/O–D stretching (S) modes, which indeed have the largest contribution at 0 K,⁸ but our detailed analysis has shown that there is another phonon group (libration) with a contribution of nearly equal magnitude, and even the smallest phonon-group (bending) contribution is on a similar order of magnitude as the final sum. It is, therefore, rather doubtful if any method short of a nearly exact one such as CCSD(T) in the CBS limit has the ability to reproduce this effect for the right reason. It would also be interesting to see if any of the more recently developed many-body potentials of water^57–64 obtained by fitting to ab initio data can reproduce this extremely subtle effect. The temperature-dependence of the VIE, on the other hand, has been reproduced by our MP2 calculations and is ascribed to the librational (L) modes. This conclusion is more robust and largely independent of calculation details, as this property does not rely on any cancellations.

MP2 has experienced difficulty optimizing the geometry of ice Ih at high pressures (>2.35 GPa) and a loss of symmetry and a volume collapse was observed at 3 GPa. These are highly suggestive of the mechanical instability of ice Ih, leading to the pressure-induced amorphization to HDA detected experimentally. The MP2 BOMD simulation at 2.5 GPa, starting with an initial crystalline geometry, has also shown a loss of symmetry and long-range order and a 20% volume reduction in a short time. This computational detection of amorphization has been further corroborated by the red-shifting of acoustic phonons accelerating around 2.0 GPa, consistent with the experimental observation of the same by Strässle et al.⁵ Our crude estimate of the thermodynamic Ih-HDA phase boundary is 1.6 GPa at 0 K with small temperature-dependence, which is reasonably close to, if larger than, the experimentally inferred value of ≈0.6 GPa. The onset of the mechanical instability is higher both in experiment (1.0–2.5 GPa) and our calculation (<2.35 GPa).

This material is based on work supported by the National Science Foundation under Award No. CHE-1361586. It is also part of the Blue Waters sustained-petascale computing project, which is supported by the National Science Foundation (Award Nos. OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint effort of the University of Illinois at Urbana-Champaign and its National Center for Supercomputing Applications. S.H. and S.Y.W. are also supported by the CREST program of the Japan Science and Technology Agency.

Here, we examine how robust the qualitative or quantitative predictions about the ice Ih anomalies are, upon changing embedding fields and basis sets.

Figures 12 and 13 show the temperature-dependence of the ice Ih volume calculated by MP2 with aug-cc-pVDZ and aug-cc-pVTZ, respectively, and the dipolar-embedding field. In both calculations as well as ones shown in Fig. 1, the thermal contraction at low temperatures followed by thermal expansion at higher temperatures is reproduced (though the latter is underestimated), indicating that this qualitative behavior is neither subtle nor difficult to describe computationally. A similar analysis (not shown) based on mode Grüneisen parameters given in Sec. III D confirms that the thermal contraction is caused almost exclusively by the hydrogen-bond bending (HB) modes (transverse acoustic phonons), while the thermal expansion is caused by the hydrogen-bond stretching (HS) modes. The underestimation of the thermal expansion is to a great extent caused by the overestimation of bulk modulus. Here, we do not observe a systematic improvement upon a basis-set extension, suggesting that approximations in the electronic structure theory, fragmentation scheme, and/or embedding field are as significant as the basis-set incompleteness.

As already mentioned in Sec. III D, MP2 with the dipolar-embedding field does not reproduce the anomalous VIE; the predicted sign of the VIE is wrong. This is traced to the fact that the sign is determined by an extremely close competition between the volume-expanding and contracting effects of all fundamental modes even at 0 K (and of their overtones and combinations at higher temperatures). Even ab initio MP2 with the aug-cc-pVTZ basis set does not seem to have sufficient fidelity to reproduce this effect qualitatively correctly, although the basis-set extension seems to bring the calculated VIE closer to the experiment. Note that dispersion-corrected DFT calculations by Pamuk et al.⁸ and by Murray and Galli²⁶ gave varied results for this quantity as well. On the other hand, the greater rate of thermal expansion in D₂O ice Ih, causing the temperature growth of the VIE (provided its anomalous value at 0 K is correctly reproduced), is consistently reproduced by all three calculations, as this is a relatively straightforward consequence of the isotope-dependence of the volume-expanding librational (L) modes.

One may at this point wonder if the atomic-embedding field is simply more accurate than the dipolar-embedding field and, therefore, the former is predicting the anomalous VIE correctly for the right reason (not accidentally). Unfortunately, this does not seem to be the case. First, as we have seen in Sec. III A, MP2 with the atomic-embedding field overbinds ice Ih to a greater extent than MP2 with the dipolar-embedding field, the latter yielding more accurate structures and binding energies. Figure 14 testifies the tendency of the atomic-embedding field to overbind the water clusters to an increasing degree with size. This is, in turn, traced to the overestimation of the dipole moment in the atomic-embedding field, as illustrated in Fig. 15. The atomic partial charges that make up the atomic-embedding field are determined so as to reproduce a short-range electrostatic field accurately³⁴ and, therefore, they do not necessarily have accurate dipole moments or the correct long-range electrostatic field.

Since MP2 with the dipolar-embedding field predicts weaker binding in ice Ih, its acoustic modes have lower, near-zero frequencies, which are, therefore, subjected to a greater degree of numerical errors such as those arising from the truncation of long-range interactions and the finite precision of optimized geometries. The volume-temperature curves in Figs. 12 and 13 should be understood to contain considerable numerical uncertainty and are used here only to demonstrate the sensitivity of the calculated VIE on the embedding fields and basis sets. The results of the atomic-embedding field discussed in the main text should be converged and reliable.