Quantum Monte Carlo study of the phase diagram of solid molecular hydrogen at extreme pressures

Establishing the phase diagram of hydrogen is a major challenge for experimental and theoretical physics. Experiment alone cannot establish the atomic structure of solid hydrogen at high pressure, because hydrogen scatters X-rays only weakly. Instead, our understanding of the atomic structure is largely based on density functional theory (DFT). By comparing Raman spectra for low-energy structures found in DFT searches with experimental spectra, candidate atomic structures have been identified for each experimentally observed phase. Unfortunately, DFT predicts a metallic structure to be energetically favoured at a broad range of pressures up to 400 GPa, where it is known experimentally that hydrogen is non-metallic. Here we show that more advanced theoretical methods (diffusion quantum Monte Carlo calculations) find the metallic structure to be uncompetitive, and predict a phase diagram in reasonable agreement with experiment. This greatly strengthens the claim that the candidate atomic structures accurately model the experimentally observed phases.

Candidate structures for phases II, III, and IV have been suggested by structure searches based on density functional theory (DFT) [26][27][28][29][30][31][32] , although it should be emphasised that none of these structures has been identified as being unambiguously correct.The candidate structures for phase II consist of packings of molecules with bond lengths almost identical to the zero-pressure value 26,27 .
We have modelled phase II using a molecular structure of P2 1 /c symmetry with 24 atoms in the primitive unit cell, which we refer to as P2 1 /c-24; see Fig. 1(a).(We adopt the convention of labelling structures by their symmetry followed by the number of atoms per primitive cell.)P2 1 /c-24 is the most stable structure found to date in static-lattice DFT within the pressure range appropriate for phase II, and its vibrational characteristics are also compatible with those of phase II.We model phase III using a C2/c-24 structure consisting of layers of molecules whose bonds lie within the planes of the layers, forming a distorted hexagonal pattern 26 ; see Fig. 1(b).This very stable structure can account for the high IR activity of phase III 26 .We also consider a molecular Cmca-12 structure 26 , which is similar to C2/c-24, but slightly denser; see Fig. 1(c).We model phase IV by a Pc-48 structure 28,29 , shown in Fig. 1(d), which consists of alternate layers of strongly bonded molecules and weakly bonded graphene-like sheets.This type of structure was predicted by Pickard and Needs 26 .Pc-48 can account for the occurrence of stiff and soft vibronic modes in phase IV, and its stabilisation by temperature.Finally, we consider the Cmca-4 structure 33 , which has weaker molecular bonds than C2/c-24 and Cmca-12, and is shown in Fig. 1(e).The main goals of our present work are to obtain accurate theoretical results for the relative stabilities of the P2 1 /c-24, C2/c-24, Cmca-12, Pc-48, and Cmca-4 structures of H at pressures of 100-400 GPa and temperatures of 0-500 K, and to use these data to construct a temperature-pressure phase diagram of H.We have not considered phase I in our calculations, which is stable at low pressures, because an accurate description of this phase would require a full quantum treatment of the proton spin.
Instead we focus our attention on the phase behaviour at higher pressures, where the candidate structures are such that the nuclei are highly localised and hence the motion of the protons is likely to be well-described by collective bosonic vibrational modes.
Useful theoretical descriptions of solid H require very accurate calculations with an energy resolution of a few meV per atom.Various studies have shown that DFT currently cannot provide such accuracy for H structures, as evidenced by the disagreement of results obtained with different exchange-correlation functionals and the fact that DFT predicts H to be metallic at pressures above ∼ 300 GPa, in contradiction with experiment 26,28,29,[34][35][36] .We have instead used the diffusion quantum Monte Carlo (DMC) method 37 to calculate static-lattice energy-volume relations for the different H phases. DMC is generally regarded as the most accurate first-principles method available for carrying out such studies [38][39][40] .Furthermore, the low mass of the H atom means that a full treatment of quantum nuclear vibrational motion, including anharmonic effects 35,40 , is crucial for an accurate description of the energetics.We have therefore used a DFT-based vibrational selfconsistent field approach 41 to calculate anharmonic vibrational energies.We find that the use of DMC (and to a lesser extent the treatment of phonon anharmonicity) renders the metallic Cmca-4 structure that is favoured in DFT energetically uncompetitive, leaving us with a phase diagram in reasonable quantitative agreement with experiment.

Results
Relative enthalpies Figure 2 shows the static-lattice enthalpies of the structures relative to C2/c-24.In Figs.2(a) and 2(b) we report DFT enthalpies calculated using the Perdew-Burke-Ernzerhof (PBE) 42 and Becke-Lee-Yang-Parr (BLYP) density functionals 43,44 .The relative DFT enthalpies are converged to better than 0.1 meV per atom with respect to k-point sampling and plane wave cutoff energy.The difference between the DFT-PBE and DFT-BLYP relative enthalpies arises chiefly from the energetics and not from the slightly different structures obtained from geometry optimisation calculations performed at fixed external pressures using the two different functionals: see Supplementary Note 1 and the accompanying Supplementary Fig. 1.In Fig. 2(c) we report DMC enthalpies, which were obtained by fitting polynomials to the extrapolated infinite-systemsize DMC energies as a function of volume, and differentiating the polynomials to obtain pressures.
The structures used for the DMC calculations were obtained from DFT-PBE geometry optimisation calculations.We truncate the DMC enthalpy curves at the highest and lowest pressures at which we have performed calculations.
The use of the DMC method has significant consequences for the static-lattice relative enthalpies of the candidate structures.Compared with both DFT-PBE and DFT-BLYP, Cmca-4 and Cmca-12 are destabilised with respect to C2/c-24, whereas P2 1 /c-24 is stabilised with respect to it, but in each case the DFT-BLYP results are closer to the DMC enthalpies, as also found in Ref. 36.
For Cmca-4 and P2 1 /c-24 the difference between the DMC and the DFT-BLYP results is greater than the difference between the DFT-BLYP and DFT-PBE results, while for Cmca-12 these differ- ences are of similar size.Although DFT-BLYP happens to be relatively accurate in the pressure range of interest, it is clear that DFT is unable to provide a consistent, quantitative description of the relative enthalpies of the phases of H.

Vibrational results
The harmonic zero-point (ZP) contributions to the enthalpies of the H phases increase sublinearly with pressure, as shown in Fig. 3(a), while the anharmonic corrections tend to decrease with pressure; see Fig. 3(b).The harmonic ZP enthalpies are roughly thirty times larger than the anharmonic corrections.However, the differences between the harmonic ZP energies of the five phases considered at fixed pressure are similar in magnitude to the differences between the anharmonic corrections, both being about 10 meV per atom, as shown in Figs.3(c) and 3(d).This demonstrates that the variations in the anharmonic vibrational corrections are as important as those of the harmonic contributions to the enthalpies in determining the relative stabilities of phases in this system.GPa difference between the transition pressures for H and D agrees well with the experimentally measured value 23 .We note that the theoretical transition pressures between H and D would only differ by around 6 GPa without the inclusion of anharmonic effects.

Structural phase transitions
As shown in Fig. 4, we also find a temperature-driven transition from C2/c-24 to Pc-48 at pressures above 250 GPa and temperatures above 300 K, in good agreement with the experimentally observed transition between phases III and IV.In Fig. 4 we show the relative free energies of C2/c-24 and Pc-48 at 300 and 400 K.At the lower temperature, C2/c-24 is marginally more stable, but at 400 K, Pc-48 has clearly become the more stable structure.The variation in the transition temperature with pressure is smaller than the uncertainty in that quantity, and so we report the C2/c-24-Pc-48 transition temperature as 320 ± 20 K.

Discussion
Our theoretical H phase diagram is in reasonable quantitative agreement with experiment, indicating that the P2 1 /c-24, C2/c-24, and Pc-48 structures provide satisfactory models for phases II, III, and IV.These model structures reproduce the experimental Raman and IR spectra quite well.However, there is a significant disagreement of about 75 GPa between the experimental and theoretical phase II-III transition pressure at 0 K.There are several possible reasons for our substantially larger phase II-III transition pressure.Firstly, the actual structure of phase III may be more stable than the C2/c-24 model structure.However, C2/c-24 is the most stable non-metallic structure found in DFT searches over a wide range of pressures, and it is compatible with the Raman and IR spectra of the observed phase III.If a significantly more stable structure than C2/c-24 were to be found for phase III, the excellent description of the transition from phase III to IV with increasing temperature obtained with our calculated data would be spoilt.Another possible explanation for the discrepancy with experiment regarding the phase II-III transition pressure could be that we have neglected a significant contribution to the energy of P2 1 /c-24 (our model for phase II).In particular, our calculations do not account for nuclear exchange effects, which are known to have a significant effect on the phase I-II transition pressure. 22However, reliable estimates of the size of nuclear exchange effects in solid H at high pressure are not currently available.Furthermore, nuclear exchange effects are expected to be much smaller in D than H, because deuterons are bosons, whereas protons are fermions, and each deuteron has twice the mass of a proton.This suggests that nuclear exchange effects cannot be entirely responsible for the discrepancy in the phase II-III transition pressure in both H and D. Our analysis of different finite-size corrections in Supplementary Note 2 (see also the accompanying Supplementary Figs. 3 and 4) indicates that finite-size effects in our relative enthalpies are well-controlled, but it is always possible that finite-size effects may be larger than anticipated.Finally, the fixed-node approximation is an uncontrolled source of error in our DMC calculations and, although fixed-node errors should largely cancel when relative energies are calculated, it cannot be ruled out that fixed-node errors may be larger in one phase than another.
The results we obtain by combining our DMC static-lattice energies and harmonic and anharmonic vibrational energies resolve a discrepancy between DFT and experiment for the transition between phases III and IV.The Pc-48 structure was proposed as a candidate for phase IV in Refs.28 and 29 because its Raman spectrum agrees well with the experimental one and because its weakly bonded layers lead to soft vibrational modes that thermally stabilise it.However, DFT static-lattice calculations together with the harmonic approximation for nuclear motion (used in Refs.28 and 29) predict the Cmca-4 structure to be energetically favoured at all temperatures in the relevant pressure range.(Note that Cmca-4 is stabilised significantly by harmonic ZP energy; at the static-lattice DFT level it is not competitive, as shown in Fig. 2.) The metallic nature of Cmca-4 contradicts experiment, in which insulating structures containing strong molecular bonds are found up to pressures in excess of 300 GPa.The phase diagram predicted by DFT is shown in Supplementary Fig. 5.The use of static-lattice DMC energies and anharmonic vibrational energies destabilises Cmca-4, and we find that it is thermodynamically unstable over the entire pressure and temperature range considered here.Our results establish that DFT does not provide even a qualitatively correct description of the phase behaviour of hydrogen.We also find that Cmca-12 is unstable at the pressures and temperatures studied in this work.We have found an important discrepancy between our calculated phase II-III transition pressure and experiment, which is currently unresolved, although we have described possible physical reasons for the disagreement.Our calculations demonstrate that anharmonic vibrational effects are crucial for determining the relative stabilities of the phases.

Methods
Quantum Monte Carlo calculations The DMC method 37,45 is capable of delivering much higher accuracy than DFT, and the scaling of the computational cost with system size enables the simula-tion of the hundreds of atoms required for accurate calculations.We have used the DMC method to calculate static-lattice energies using H structures relaxed within DFT-PBE at a given external pressure.In DMC, the ground-state component of a trial wave function is projected out by simulating the Schrödinger equation in imaginary time, subject to the constraint that the nodal surface of the wave function is fixed to be that of the trial wave function 37,45 .We used Slater-Jastrow wave functions as implemented in the CASINO code 46 .Full technical details of our calculations can be found in Supplementary Note 2. The single-particle orbitals were obtained from the CASTEP code 47 using the PBE exchange-correlation functional.The nuclei were represented by bare Coulomb potentials and appropriate cusp corrections were applied to the orbitals.We used a flexible Jastrow factor 48 whose parameters were optimised using variational Monte Carlo (VMC) 49 .VMC and DMC simulations were performed using 96 and 768 atoms, and the results were extrapolated to infinite system size.Using the resources of the Oak Ridge Leadership Computing Facility, we achieved statistical error bars of less than 0.3 meV per atom in all our DMC calculations.

Anharmonic vibrational calculations
We have calculated harmonic vibrational free energies by using the finite-displacement method to construct the matrix of force constants and diagonalising the corresponding dynamical matrices over a fine vibrational Brillouin-zone grid, as described in Supplementary Note 3 (with accompanying data presented in Supplementary Figs. 6 and 7).
We determined anharmonic corrections to the harmonic free energies using a vibrational selfconsistent field method 40,41,50 , sampling the low-energy part of the DFT-PBE Born-Oppenheimer (BO) energy surface along harmonic normal modes to large amplitudes.The resulting anharmonic

Figure 1 :
Figure 1: Atomic structures of the five H phases considered in this work.(a) P2 1 /c-24, (b)

Figure 2 :
Figure 2: DFT and DMC static-lattice enthalpy-pressure relations for the different H struc-

Figure 4 Figure 3 :
Figure 3: DFT-PBE vibrational contributions to the enthalpies of the H structures. (a) Har-