Observation of methane filled hexagonal ice stable up to 150 GPa.

Gas hydrates consist of hydrogen-bonded water frameworks enclosing guest gas molecules and have been the focus of intense research for almost 40 y, both for their fundamental role in the understanding of hydrophobic interactions and for gas storage and energy-related applications. The stable structure of methane hydrate above 2 GPa, where CH4 molecules are located within H2O or D2O channels, is referred to as methane hydrate III (MH-III). The stability limit of MH-III and the existence of a new high-pressure phase above 40 to 50 GPa, although recently conjectured, remain unsolved to date. We report evidence for a further high-pressure, room-temperature phase of the CH4-D2O hydrate, based on Raman spectroscopy in diamond anvil cell and ab initio molecular dynamics simulations including nuclear quantum effects. Our results reveal that a methane hydrate IV (MH-IV) structure, where the D2O network is isomorphic with ice Ih, forms at ∼40 GPa and remains stable up to 150 GPa at least. Our proposed MH-IV structure is fully consistent with previous unresolved X-ray diffraction patterns at 55 GPa [T. Tanaka et al., J. Chem. Phys. 139, 104701 (2013)]. The MH-III → MH-IV transition mechanism, as suggested by the simulations, is complex. The MH-IV structure, where methane molecules intercalate the tetrahedral network of hexagonal ice, represents the highest-pressure gas hydrate known up to now. Repulsive interactions between methane and water dominate at the very high pressure probed here and the tetrahedral topology outperforms other possible arrangements in terms of space filling.

Gas hydrates spontaneously form in Nature and can be produced in the laboratory by mixing ice with small hydrophobic gas molecules at modest pressures (100 bar) and temperatures close to ice melting. Under these conditions, they adopt crystalline clathrate structures, where the guest molecules are confined in nanometric polyhedral cages that are formed by the hydrogen-bonded H2O or D2O network. The amount of gas that can be stored in these structures is very significant, orders of magnitude beyond what is expected from the solubility of gas in water at normal conditions. Typical water to gas ratios are indeed close to 6 in most ice clathrates, even though their stoichiometry strongly depends on their formation conditions. Energy-relevant hydrocarbons, like methane, and greenhouse gases, like carbon dioxide, easily form clathrates, which are therefore expected to play a role for energy harvesting and in addressing the global demand for fossil energy and CO2 mitigation (1). For this reason, and for their role as model systems for the study of the hydrophobic interaction, gas clathrates have been the subject of intense research for almost 40 y (2). Rather recently (3–6) experiments have shown that clathrate structures do not collapse when subject to high pressures, but transform into denser structures with even smaller water to gas ratios. These structures, known as filled ices, consist of water networks that stem from distorted ice phases, hosting gas molecules in their voids or channels. This way, low-density ice phases, as ice Ic or ice II, can be stabilized at much higher pressure than in pure ice upon inclusion of gas species (5, 6)> and a new metastable phase of ice (ice XVII) could be obtained by emptying a high-pressure hydrogen hydrate recovered at ambient pressure (7).

Understanding the behavior of gas–water interactions in filled ices at high pressures is also crucial for accurate models of planetary interiors. For instance, Titan is expected to have volcanic eruptions of fluid methane, stored in water and ammonia ices (8). Space probes have sporadically detected the presence of methane in the atmosphere of Mars (9), and it was suggested that its origin might be linked to the destabilization of hydrates existing at depths (10). Methane is also abundant on the giant planets Uranus and Neptune where it was the product of chemical processing of the primordial solar nebula material and then included under the form of ice clathrates during the planetary accretion. In these planets methane and water coexist at pressures from gigapascals (GPa) to hundreds of GPa.

The high-pressure hydrate structure formed by methane—and also krypton, argon, and nitrogen—at pressures of a few GPa was first observed by Loveday et al. (3) and is generally referred to as methane hydrate III (MH-III). The hydrogen-bonded, proton-disordered water network of MH-III is similar to the ice Ih structure when looking along the axis, where 6-member rings of water molecules align to form hexagonal channels. In addition, 4- and 8-member rings, not existing in ice Ih, are also present along the other crystal directions (5). The water to methane ratio in MH-III is close to 2, although the precise stoichiometry is likely linked to the filling ratio of the precursor clathrate phase. Below 3 GPa, the dynamics of the enclathrated methane are not very different from those of the isolated molecule at low pressures. As pressure increases in MH-III, the O–C distances shorten and methane vibrations become more and more entangled with those of the neighboring water molecules. Above GPa, the enclathrated methane appreciably deviates from the tetrahedral symmetry and shows orientational ordering (11, 12)>.

Recent studies also suggested that under high-pressure conditions methane and water could form stoichiometric mixed molecular crystals beyond the known MH-III phase: The existence of a new high-pressure phase in methane hydrate at pressures beyond 40 to 50 GPa was first hinted by Hirai and coworkers (13–16), who observed a slow transition from MH-III into a different, yet unsolved structure during synchrotron X-ray diffraction experiments. The high-pressure structure was observed to survive up to 86 GPa, but the diffraction patterns did not allow the authors to identify it in detail.

In this paper we report results of high-resolution high-pressure Raman spectroscopy and ab initio molecular dynamics (AIMD) simulations including nuclear quantum effects (NQE) on CH4–D2O methane hydrate at pressures up to 150 GPa. Additional vibrational modes show up in the Raman spectra at GPa, which are not compatible with the MH-III phase, together with a significant change in the width of the methane stretching modes. AIMD simulations indicate the existence of a distinct high-pressure methane hydrate phase, in the following named methane hydrate IV (MH-IV), whose water skeleton is isomorphic with ice Ih. We note that 2 distinct hypothetical structures, not experimentally observed so far, were named MH-IV and MH-V in a recent work (17). Based on the Roman numbering scheme that has been adopted for the ice phases (18), and for the sake of clarity, we named the present methane hydrate phase MH-IV, which naturally follows the MH-III phase by increasing pressure. This phase is stable in the probed pressure range up to 150 GPa and accounts for the additional detected Raman modes. The MH-IV structure is also consistent with previous X-ray diffraction investigations (13–16) and represents the highest-pressure gas hydrate phase known up to now. The transition from MH-III to MH-IV is shown to proceed in 3 steps as pressure increases: 1) symmetrization of the hydrogen bonds in MH-III at 30 GPa (giving rise to a structure we name “MH-IIIs”), 2) structural change from MH-IIIs to the MH-IV phase at 40 GPa, and 3) symmetrization of the hydrogen bonds of MH-IV at 50 GPa (giving rise to “MH-IVs”).

Results and Discussion

We performed Raman spectroscopy measurements in diamond anvil cells (DACs) on a CH4–D2O methane hydrate sample at room temperature and pressures between 3 and 150 GPa, with pressure steps of about 5 GPa. Measurements of the same sample up to 45 GPa were previously performed by our group, using the same experimental apparatus but different diamonds, and are reported in ref. 11. Details of sample preparation and Raman spectra are shown in . The frequencies of the Raman-active modes, derived from a fit of the Raman spectra, are reported in Fig. 1 as a function of pressure. With increasing pressure, the 2 symmetric and antisymmetric CH4 stretching modes shifted to higher frequencies and could be clearly observed up to the maximum investigated pressure of 150 GPa. The lattice mode also shifted to higher frequencies upon compression. It weakened considerably and became undetectable above 50 GPa. Similarly, the CH4 rocking mode shifted to higher frequencies upon compression and broadened considerably above 35 GPa. It had almost no detectable intensity above 50 to 60 GPa. The 2 D2O stretching modes shifted to smaller frequencies and lost intensity swiftly between 35 and 45 GPa, as previously reported in the literature for methane hydrate (15, 16)> and in strong analogy to what is commonly observed in pure ice (19) and salt-doped ices (20). In our previous work (11), we described in detail the mode mixing occurring between the stretching mode of D2O and the rocking mode of CH4 in this pressure range. A peak at 950 appeared upon pressure increase at GPa. As shown in Fig. 1, Inset the peak has 2 components and is well fitted by 2 pseudo-Voigt functions. By its frequency range and pressure variation, the lower-frequency excitation might correspond to the O–O vibrational mode of symmetry that is a marker of the transition to the symmetric H-bonded cubic phase of pure ice (ice X) (19). It must be noted that our hydrate sample contained excess D2O—formed during the sI to sH and sH to MH-III transitions—which was in the form of ice VII or ice X in the probed pressure range and contributed substantially to the measured lattice and D2O stretching modes. However, the appearance of the higher-frequency excitation, which stems neither from pure ice nor from MH-III, suggests the existence of a high-pressure phase of methane hydrate above 100 GPa. The first indication of a high-pressure phase of methane hydrate (for GPa) was provided by X-ray diffraction experiments (14–16). They showed peaks appearing in the diffraction patterns beyond 40 GPa, which are not compatible with the MH-III structure. Moreover, MH-III presents nonnegligible angular frustration for the water and methane molecules (11), which increases upon compression and might weaken the whole structure. All of the previous results indicate that a different MH structure is formed, although experiments did not allow identifying it uniquely. To determine this phase, we performed AIMD simulations including NQE () for pressures up to 150 GPa at room temperature. We explored several candidate structures, inspired from other ice and filled-ice phases. Along the simulations of MH-III at high pressures, we observed changes of the H-bonded network during a few dynamical events. By following this route, we found a structure, the MH-IV phase, shown in Fig. 2, which becomes more stable than MH-III around 30 GPa (Fig. 3). This structure reproduces all of the features of our Raman study and previous X-ray diffraction data (Figs. 1 and 3).

Fig. 1.

Experimental Raman frequencies (circles) derived from the fitting procedure described in and theoretical frequencies (lines) upon compression of the methane stretching and rocking modes, D2O stretching modes, the lattice mode, and the MH-IV characteristic modes (IV modes). (Inset) Experimental Raman spectra showing the IV modes appearing beyond 100 GPa. Theoretical frequency values were computed through the velocity time correlation functions as extracted from AIMD trajectories, with the only exception of D2O stretching modes, which were calculated in the harmonic crystal framework.

Fig. 2.

Ball-and-stick representation of the high-pressure MH-IVs phase, seen in the 3 crystal planes. O atoms are in red, C atoms in black, and H atoms in light gray. Note in Upper panel that methane molecules have a C–H bond oriented either along (left) or (right).

Fig. 3.

(A) X-ray diffractogram at 55 GPa, showing with red arrows the peaks associated with the additional phase, along with the simulated diffractogram ( Å) for the MH-IV structure. Peak indexes can be found in . Data from ref. 15. (B) FWHM of the CH4 and stretching modes vs. pressure obtained from our Raman spectra. The FWHM changes slope around 80 GPa, which can be linked to the degeneracy lifting of the methane stretching modes upon compression. (C) Computed free enthalpy difference between MH-III and MH-IV structures at T = 0 K.

The MH-IV phase has an orthorhombic crystal cell based on an ice Ih skeleton with the same stoichiometry as MH-III, i.e., a 2:1 D2O:CH4 ratio. This indicates that the MH-III MH-IV transition can happen at constant chemical potential, without need of providing or removing atoms. The space group of MH-IV is Pmcn (n°62) where oxygen and carbon atoms occupy 4c Wyckoff sites (Table 1). MH-IV contains ordered methane molecules and differs from MH-III where both ice network and methane ordering are concerned. The water network of MH-IV is much closer to ice Ih than that of the MH-III phase, which is substantiated by a quantitative estimate of the topological distance (i.e., relying on the reduced coordinates) between the oxygen atoms of those 2 hydrate phases and ice Ih. Indeed in MH-III there are 4-, 6-, and 8-fold water rings while in MH-IV the H-bonded water network is composed of corrugated (a–b) sheets of edge-sharing 6-member rings of water molecules which are cross-linked along as in ice Ih, thus avoiding the angular frustration that 4-member water rings imply.

Table 1.

Space group and fractional coordinates of oxygen and carbon atoms in the MH-III, MH-IIIs, MH-IV, and MH-IVs phases of methane hydrate at 40 GPa, as obtained from the AIMD simulations

Phase	Group	Atom	Site	x	y	z
MH-III	Imcm	O	8i	0.250	0.080	0.830
	74	C	4e	0.250	0.680	0.000
MH-IIIs	Pmcn	O	4c	0.250	0.400	0.190
	62	O	4c	0.250	0.420	0.810
		C	4c	0.250	0.800	0.980
MH-IV	Pmcn	O	4c	0.250	0.420	0.045
MH-IVs	62	O	4c	0.250	0.420	0.455
		C	4c	0.250	0.750	0.715

For coherence with Imcm, the space group of the higher-pressure phases is given in the nonconventional representation. The uncertainty on the fractional coordinates is ±0.005. Despite the different space group between MH-III and MH-IIIs, the full oxygen positions are very similar. In MH-IVs (H-bond symmetric phase IV), the oxygen and carbon atomic positions do not change significantly from those in MH-IV. Computed lattice parameters: MH-IIIs (40 GPa), a=4:006 Ǻ , b=6:911 Ǻ , c=6:249 Ǻ ; MH-IV (40 GPa), a=4:063 Ǻ , b=6:981 Ǻ , c=6:063 Ǻ . Lattice parameter trends with pressure are reported in .

In MH-IV, the carbon atoms and a C–H bond in each methane molecule are aligned along the direction, where 2 symmetrical configurations that differ by a rotation of around alternate (Fig. 2, Lower Left). Along , methane ordering presents another alternation by a rotation of around (Fig. 2, Upper). In this arrangement, all C–H bonds of the methane point toward the hexagonal channels formed by the host O molecules, thus reducing repulsive methane–water interactions. Our optimized lattice parameters at 50 GPa for the MH-IV structure agree well with those proposed at the same pressure in ref. 14 for the hypothetical (at the time) orthorhombic structure referred to as post filled-ice–Ih (). While and increase a little after the MH-III MH-IV transition, the major change concerns , which decreases by about 3.2% due to the particular methane ordering along this axis (Table 1). In addition, Fig. 3 shows that the MH-IV structure quantitatively reproduces the location and intensity ratios of the main peaks in the X-ray diffraction data from ref. 15, which is a strong indication that the atomic positions of the C and O atoms of the simulated MH-IV structure match the experimental ones. The structural properties of MH-IV reproduce the whole information about the hypothetical phase that was invoked by Hirai and coworkers (13–16). We believe that such an excellent agreement is not accidental and that the MH-IV structure we propose corresponds to the unsolved one as obtained in ref. 15.

Fig. 3 presents the free enthalpy difference at K as a function of pressure, obtained from static calculations. Accordingly, the MH-IV structure becomes more stable than MH-III for pressures greater than 30 GPa. The transition pressure as established experimentally (14) is 40 GPa. Such a difference can be attributed to several factors: First, the chosen approximation to the exchange-correlation energy within the DFT; second, the fact that static energy calculations take into account neither thermal nor quantum effects; and last but not least, it can also be due to kinetic hindrance. Therefore, from now on we adopt 40 as the transition pressure. Upon further compression, the free enthalpy difference grows up to 400 meV per methane molecule at 150 GPa. Its trend shows that the stability of the MH-IV phase, with smaller volume, increases steadily with pressure. The transition from MH-III to MH-IV is accompanied by the formation of 6-fold water rings, at the expense of 4- and 8-fold ones. As pressure increases and hydrogen bonds contract, the angular frustration of O molecules within 4-fold rings causes a fast augmentation of the gap between the corresponding internal energies, which consistently keeps rising with pressure.

The phase transition from MH-III to MH-IV can also be detected by looking into the behavior of the vibration modes as a function of the increasing pressure. In the MH-III structure, the methane molecules undergo an increasing distortion from tetrahedral symmetry upon compression (11); at GPa, a angle diminishes while another widens by with respect to 109.. The departure from symmetry implies the loss of degeneracy of methane stretching and bending modes. After the transition to MH-IV, the angles of methane almost recover their tetrahedral value. Fig. 1 presents the computed CH4 stretching mode frequencies in MH-III and MH-IV, compared with those obtained through Raman spectroscopy. As shown in , the experimental dispersion with pressure of the CH4 stretching-mode frequencies is accounted for by the computed trends for MH-III below 30 GPa and by those calculated for the MH-IV structure above this pressure threshold. The experimental full width at half maximum (FWHM) of the 2 methane stretching modes changes slope around 80 GPa (Fig. 3). Scrutinizing the outcomes of the simulations, the methane molecules deviate from the tetrahedral symmetry in 2 respects: First, the C–H bond lengths along the axis (which point toward a neighboring methane molecule; Fig. 2) are shorter than the others; and second, the angles progressively vary from 109. as pressure increases (). Both phenomena quicken with pressure around 70 to 80 GPa and the band of the CH4 stretching modes broadens, which we relate with the significant increment of the FWHM in the Raman spectra (Fig. 3). The pressure trend of the antisymmetric stretching mode in MH-IV follows that in pure methane (21). These findings are consistent with dominant repulsive interactions and the absence of hydrogen bonds between the embedded molecules and the water network.

More strikingly, Raman spectroscopy reveals 2 additional modes in the 900- to 1,100- range, which show up at 100 GPa and increase their intensities as pressure rises (Fig. 1). Although 1 of the 2 modes is superimposed in the Raman spectra to the mode of the residual ice after transition to the symmetric phase ice X, their characteristics are specific to methane hydrate at high pressure. Our theoretical analysis, within the harmonic approximation, describes these modes as being shared among the water and methane degrees of freedom in the MH-IV structure. We denote them as IV modes, as they are characteristic of the latter structure. They are due to the strong repulsive interaction between the guest methane and the host water network in MH-IV; an analogous enhancement of mode mixing with pressure was also detected in MH-III (11), although at higher frequencies. As shown in Fig. 1, the mode dispersion upon compression extracted from our MD simulations is in agreement with Raman results. Thus, the appearance of these modes can also be considered as a signature of the MH-IV phase.

We point out that there is no signature of free methane in the measured spectra and that both the methane stretching and the IV modes are remarkably stable and recognizable with increasing pressure. These facts show that the MH-IV phase remains stable at least up to the maximum investigated pressure of 150 GPa, i.e., well beyond all previous observations in gas hydrates. Our simulations and theoretical analysis indicate that the repulsive interactions between host and guest molecules still dominate after the transition and that, up to 150 GPa, there is no evidence of hydrogen bond formation between methane and water molecules.

The transition from MH-III to MH-IV is rather complex and passes through several stages, which involve both H-bond symmetrization and structural reorganization. First, in our previous study (11), we demonstrated and characterized the methane orientational ordering occurring at GPa. The latter induces both a noticeable departure of the methane molecule from the tetrahedral symmetry and a considerable entanglement of methane and water vibrations. This structure was previously guessed and referred to as MH-GOS (guest-ordered state) (16). Then, in the range 30 to 40 GPa, the system undergoes a transition to symmetric H bonds (12). We note the 2 phases as MH-III and MH-IIIs, respectively. We used path integral molecular dynamics (PIMD) to characterize this transition, which is driven by NQE as in pure ice (22, 23)>. We also adopted as the order parameter, which measures the asymmetry of the bonds. Comparison of the probability densities, , calculated from quantum and classical simulations, shows indeed that NQE are not negligible as they downshift the transition pressure by about 15 GPa (). Then, at 40 GPa, the MH-IIIs MH-IV transition occurs. A possible transition mechanism is proposed in Fig. 4. It requires a structural reorganization of both the water skeleton and the hydrogen bond water network as well as a significant change of the methane orientation, which essentially impacts half of the unit cell (Fig. 4, Center panels). According to classical calculations within the nudged elastic band (24), some hydrogen bonds that form the 4-fold water rings in MH-IIIs break. Then, half of the molecules reorganize, through translations and rotations, with 1 C–H bond of 4 aligned along . Finally, the remaining 4-fold water rings are broken to form the MH-IV structure. Quantum corrections, computed through PIMD, change this picture and strongly decrease the barriers corresponding to H-bond breaking, mainly because of deuterium quantum spread due to tunneling and zero-point energy. Indeed, hydrogen bonds undergo massive D tunneling back and forth through the symmetric position (see distributions in Fig. 4). Therefore, deuterium ordering in MH-IV cannot be identified, as the water network is more comparable to ice Ih at the onset of the hydrogen bond symmetrization. The restructuring accompanying the MH-III to MH-IV transition increases the O–O distances in the () plane. Indeed, although the volume of the structure is reduced at the transition, only the lattice parameter decreases while the 2 others increase (Table 1), giving rise to slightly longer O–O distances, so that hydrogen bonds are almost, but not fully, symmetric. The further pressure increase therefore triggers another transition, from MH-IV to MH-IVs, where hydrogen bonds become symmetric (Fig. 4).

Fig. 4.

(A) Sketch of the main transitions from MH-III to MH-IV, according to AIMD simulations. Changes are represented in the () plane, in blue for the hydrogen-bonded water network and in orange for the methane molecules. While passing from MH-III to MH-IV, mainly 2 hydrogen bonds are impacted, changing one of the oxygen atoms they were bonded to. Along with the hydrogen bond reorganization, the oxygen atoms involved are displaced along , forming 6-fold rings as in ice Ih. The methane molecules also present several rearrangements: A displacement along and a small rotation along lead to a perfect alignment of one of their C–H bonds along ; a rotation of around gives rise to the methane orientational ordering of MH-IV. (B) Probability distributions of the order parameter (main text) of MH-III at 30 GPa, MH-IIIs at 40 GPa, MH-IV at 40 GPa, and MH-IVs at 60 GPa, as obtained by PIMD.

Finally, we tested other hypothetical structures with the MH-IV positions for the oxygen atoms and the MH-III positions for the methane, or the inverse; they always spontaneously relaxed to the MH-IV structure. Not only is the steric hindrance of the methane molecules important but also their (quasi-)tetrahedral symmetry seems to be a key parameter for this transition: Replacing them with high-radius rare gas atoms resulted in a different structure. This strengthens the emerging view of high-pressure methane hydrate as a strongly interacting inclusion compound, where a global optimization of both water and methane molecules should be fulfilled.

Conclusions

To summarize, both Raman spectroscopy experiments in diamond anvil cells and ab initio molecular dynamics simulations evidence the existence of a methane hydrate phase we named MH-IV, which is stable well beyond the probed stability limit of other hydrates (5). This phase can be pinpointed in the Raman spectra by the presence of 2 vibrational modes at frequencies between 900 and 1,100 , which are absent in the MH-III phase, as well as a marked broadening of the CH4 stretching modes. Simulations reveal that this phase consists of a skeleton isomorphic to ice Ih and orientationally ordered methane molecules forming an intercalated sublattice, which mimics 1 of the 2 D2O sublattices as found in some high-pressure ice phases. The CH4 molecules are arranged in rows along a ternary axis. Our proposed structure fully accounts for the X-ray diffractogram measured at 55 GPa in ref. 15 and the lattice parameters proposed in ref. 14, where the existence of a high-pressure hydrate phase was conjectured. Simulations depict the transition pattern from MH-III to MH-IVs according to a complex transition sequence: First, in the 30 to 40 GPa range, H bonds become symmetric in MH-III, leading to MH-IIIs. Second, the transition from MH-IIIs to MH-IV implies a reorganization of the H-bonded water network into an ice Ih-like phase, which requires 2 hydrogen bonds to move as well as a rotation of half of the methane molecules around . This structure belongs to the Pmcn space group and presents symmetric H bonds beyond 50 GPa, provided that quantum effects on the nuclear distributions are properly accounted for.

Molecular mixtures, where the single components do not, or only weakly, interact at moderate pressures, can form stable crystals when forced into a tight geometry. In the water–methane system we observe some striking pressure-induced phenomena: On the one hand, water recovers its stable low-pressure hexagonal Ih phase, promoted by the pressure-induced ordering of the guest methane molecules and by their strong repulsive interaction, which prevails over the formation of hydrogen bonds between the 2 species. On the other hand, the ice lattice could preserve the methane from dissociation, which has been observed in some (25) but not all (21) Raman experiments on pressurized methane at ambient temperature and suggested by ab initio calculations (26). We point out that the ice Ih frame matches quite well the tetrahedral symmetry and the size of the methane molecule, so that the local tetrahedral topology with the 2:1 stoichiometry simply outperforms other possible arrangements in terms of space filling. Despite the repulsive interaction, the specific volume of this methane hydrate phase remains indeed close to the sum of those of ice and methane elemental phases in the explored pressure range, which encompasses the expected pressures on the outer planets of the Solar System.

Materials and Methods

Experimental Details.

D2O–CH4 clathrate samples were prepared using the well-established method described in ref. 27. For Raman measurements, samples were loaded cold and compressed in a diamond anvil cell to have enough pressure to stabilize the sample at room temperature. All spectra were measured at room temperature, upon increasing pressure in small steps. Further details on sample preparation, Raman measurements, and analysis are given in .

Simulation Methods.

Both MH-III and MH-IV phases were simulated via ab initio molecular dynamics at room temperature. Electronic structures and interatomic forces were determined by density functional theory simulations within the generalized gradient approximation (28) provided by the Quantum Espresso package (29). Nuclear quantum effects were taken into account by quantum thermal bath (30), for dynamical properties and path integral MD, for the quantum statistics counterpart via the i-PI package (31). The reader can find further information in .