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

Abstract

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.

ARTICLE

Received 9 Oct 2014 | Accepted 10 Jun 2015 | Published 28 Jul 2015

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 identied 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) nd 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.

DOI: 10.1038/ncomms8794 OPEN

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

N.D. Drummond1, Bartomeu Monserrat2, Jonathan H. Lloyd-Williams2, P. Lpez Ros2, Chris J. Pickard3

& R.J. Needs2

1 Department of Physics, Lancaster University, Lancaster LA1 4YB, UK. 2 Theory of Condensed Matter Group, Cavendish Laboratory, University of Cambridge,J. J. Thomson Avenue, Cambridge CB3 0HE, UK. 3 Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK. Correspondence and requests for materials should be addressed to N.D. (email: mailto:[email protected]

Web End [email protected] ).

NATURE COMMUNICATIONS | 6:7794 | DOI: 10.1038/ncomms8794 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 1

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8794

Hydrogen (H) is the simplest and most abundant of all elements and yet it displays amazing richness in its phase behaviour1,2: it is observed to form a quantum crystalline

state and orientationally ordered molecular phases, and it has been predicted to exhibit a liquid-metal phase at high pressures and low temperatures35, metallic superuid and superconducting superuid states6,7 and high-temperature superconductivity810. Several crystalline phases of solid molecular H have been observed in diamond anvil cell experiments carried out at pressures up to over 300 GPa (refs 1119). The low-pressure phase I, which is a hexagonal close-packed structure formed of freely rotating molecules, transforms to a broken-symmetry phase II, in which the molecular rotations are restricted, at low temperatures1,2. The transition pressure decreases strongly with isotopic mass1,2023 and also depends on the total spin of the molecules1,22. As the pressure is increased at low temperatures, there is a further transition from phase II to a phase III at about 160 GPa, with the transition pressure for deuterium (D) exceeding that for H by about 12 GPa (ref. 23). Experimental studies have also demonstrated the existence of a phase IV at temperatures above a few 100 K and pressures above 220 GPa (refs 1214,17,16). Some constraints on the structures of the observed phases have been obtained from X-ray diffraction experiments24,25, but the low X-ray scattering cross section of H and the small sample sizes available limit the possible resolution. Infrared (IR) and particularly Raman spectroscopic measurements have yielded valuable information about the vibrational modes of H at high pressures1124, but the available experimental data are insufcient to determine the structures of phases II, III and IV.

Candidate structures for phases II, III and IV have been suggested by structure searches based on density functional theory (DFT)2632, although it should be emphasized that none of these structures has been identied as being unambiguously correct. The candidate structures for phase II consist of packings of molecules with bond lengths almost identical to the zero-pressure value26,27. We have modelled phase II using a molecular structure of P21/c symmetry with 24 atoms in the primitive unit

cell, which we refer to as P21/c-24; see Fig. 1a. (We adopt the convention of labelling structures by their symmetry followed by the number of atoms per primitive cell.) P21/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 pattern26; see Fig. 1b. This very stable structure can account for the high IR activity of phase III26. We also consider a molecular Cmca-12 structure26, which is similar to C2/c-24, but slightly denser; see Fig. 1c. We model phase IV by a Pc-48 structure28,29, shown in Fig. 1d, which consists of alternate layers of strongly bonded molecules and weakly bonded graphene-like sheets. This type of structure was predicted by Pickard and Needs26. Pc-48 can account for the occurrence of stiff and soft vibronic modes in phase IV, and its stabilization by temperature. Finally, we consider the Cmca-4 structure33, which has weaker molecular bonds than C2/c-24 and Cmca-12, and is shown in Fig. 1e. The main goals of our present work are to obtain accurate theoretical results for the relative stabilities of the P21/c-24, C2/c-24, Cmca-12, Pc-48 and Cmca-4 structures of H at pressures of 100400 GPa and temperatures of 0500 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 localized 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

a b c

d e

Figure 1 | Atomic structures of the ve H phases considered in this work. (a) P21/c-24, (b) C2/c-24, (c) Cmca-12, (d) Pc-48 and (e) Cmca-4.

The blue dumbbells indicate short bonds between atoms (o0.8 ). The white dumbbells indicate long bonds between atoms (o 0.9 ). The red lines indicate close contacts between atoms (o1.2 ) in the layered structures. P21/c-24 consists of molecules arranged on a distorted hexagonal close-packed lattice. C2/c-24, Cmca-12 and Cmca-4 consist of layers of molecules whose bonds lie within the planes of the layers, forming distorted hexagonal patterns, and we show topdown views of single layers. Pc-48 consists of alternate layers of isolated strongly bonded molecules and weakly bonded graphene-like sheets, and we show a topdown view of four layers. The structures are shown at a common DFTPBE pressure (250 GPa).

2 NATURE COMMUNICATIONS | 6:7794 | DOI: 10.1038/ncomms8794 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8794 ARTICLE

a

b

c

25

25

Relative enthalpy (meV per atom)

25

Relative enthalpy (meV per atom)

Relative enthalpy (meV per atom)

C 2/c-24 Cmca-12

Cmca- 4 P 21/c-24

Pc-48

C 2/c-24 Cmca-12

Cmca- 4 P 21/c-24

Pc-48

C 2/c-24 Cmca-12 Cmca- 4 P 21/c-24

Pc- 48

20

20

20

15

15

15

10

10

10

5

5

5

0

0

0

5

5

5

10 100 150 200 250 300 350 400Pressure (GPa)

10 100 150 200 250 300 350 400Pressure (GPa)

10 100 150 200 250 300 350 400Pressure (GPa)

Figure 2 | DFT and DMC static-lattice enthalpypressure relations for the different H structures relative to C2/c-24. (a) DFTPBE, (b) DFTBLYP and (c) DMC. The relative DFTenthalpies are converged to better than 0.1 meV per atom. The widths of the DMC lines indicate the estimated uncertainties in the enthalpies due to single-particle nite-size errors, which are greater than the uncertainties due to random sampling in the Monte Carlo algorithm, as explained in Supplementary Note 2.

240 100 150 200 250 300 350 400Pressure (GPa)

a

b

340

Harmonic ZP enthalpy (meV per atom)

20

Anharmonic correction (meV per atom)

320

C 2/c-24 Cmca-12 Cmca- 4

P 21/c-24 Pc-48

C2/c-24 Cmca-12 Cmca- 4 P 21/c-24

Pc-48

15

300

10

280

5

260

0 100 150 200 250 300 350 400 Pressure (GPa)

c

10

d

10

8

8

Relative harmonic enthalpy

(meV per atom)

Relative anharmonic correction

(meV per atom)

6

6

4

4

2

2

0

0

2

2

C 2/c-24 Cmca-12 Cmca- 4

P 21/c-24 Pc- 48

4

4

6

C 2/c-24 Cmca-12 Cmca- 4

P 21/c-24 Pc- 48

6

8

8

10 100 150 200 250 300 350 400Pressure (GPa)

10 100 150 200 250 300 350 400Pressure (GPa)

Figure 3 | DFTPBE vibrational contributions to the enthalpies of the H structures. (a) Harmonic zero-point (ZP) contributions to enthalpies,(b) anharmonic ZP corrections to enthalpies, (c) harmonic ZP enthalpies relative to C2/c-24, and (d) anharmonic ZP corrections relative to C2/c-24. P21/c-24 is destabilized by both harmonic vibrations and anharmonic corrections, relative to C2/c-24. Cmca-12, Cmca-4, and Pc-48 are all stabilised by harmonic vibrations but destabilized by anharmonic corrections, relative to C2/c-24.

above 300 GPa, in contradiction with experiment26,28,29,3436. We have instead used the diffusion quantum Monte Carlo (DMC) method37 to calculate static-lattice energyvolume relations for the different H phases. DMC is generally regarded as the most accurate rst-principles method available for carrying out such studies3840. Furthermore, the low mass of the H atom means that a full treatment of quantum nuclear vibrational motion, including anharmonic effects35,40, is crucial for an accurate description of the energetics. We have therefore used a DFT-based vibrational self-consistent eld approach41 to calculate

anharmonic vibrational energies. We nd 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.

ResultsRelative enthalpies. Figure 2 shows the static-lattice enthalpies of the structures relative to C2/c-24. In Fig. 2a,b we report DFT enthalpies calculated using the Perdew-Burke-Ernzerhof (PBE)42

NATURE COMMUNICATIONS | 6:7794 | DOI: 10.1038/ncomms8794 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 3

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8794

a

10

b

10

c

10

d

10

C 2/c-24 P 21/c-24

C 2/c-24 P 21/c-24

C 2/c-24 Pc-48

C 2/c-24 Pc- 48

Relative Gibbs free energy

(meV per atom)

Relative Gibbs free energy

(meV per atom)

Relative Gibbs free energy

(meV per atom)

Relative Gibbs free energy

(meV per atom)

5

5

5

5

0

0

0

0

5

5

5

5

0 K

150 K

300 K

400 K

10 150 200 250 Pressure (GPa)

10 150 200 250 Pressure (GPa)

10 250 300 350 Pressure (GPa)

10 250 300 350 Pressure (GPa)

Figure 4 | Relative Gibbs free energies of the different H structures. (a) 0 K, (b) 150 K, (c) 300 K and (d) 400 K. The Gibbs free energies were calculated using static-lattice DMC calculations together with DFTPBE harmonic and anharmonic vibrational calculations. The transition from P21/c-24

to C2/c-24 occurs at around 23510 GPa between 0 and 150 K. Pc-48 is stabilised by temperature with respect to C2/c-24. The complete set of relative enthalpies is shown in Supplementary Fig. 2.

0 150 200 250 300 350Pressure (GPa)

and BeckeLeeYangParr (BLYP) density functionals43,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 DFTPBE and DFTBLYP relative enthalpies arises chiey from the energetics and not from the slightly different structures obtained from geometry optimization calculations performed at xed external pressures using the two different functionals: see Supplementary Note 1 and the accompanying Supplementary Fig. 1. In Fig. 2c, we report DMC enthalpies, which were obtained by tting polynomials to the extrapolated innite-system-size DMC energies as a function of volume, and differentiating the polynomials to obtain pressures. The structures used for the DMC calculations were obtained from DFTPBE geometry optimization 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 signicant consequences for the static-lattice relative enthalpies of the candidate structures. Compared with both DFTPBE and DFTBLYP, Cmca-4 and Cmca-12 are destabilized with respect to C2/c-24, whereas P21/c-24 is stabilized with respect to it, but in each case the

DFTBLYP results are closer to the DMC enthalpies, as also found in ref. 36. For Cmca-4 and P21/c-24 the difference between the DMC and the DFTBLYP results is greater than the difference between the DFTBLYP and DFTPBE results, while for Cmca-12 these differences are of similar size. Although DFTBLYP 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 contributions to the enthalpies of the H phases increase sublinearly with pressure, as shown in Fig. 3a, while the anharmonic corrections tend to decrease with pressure; see Fig. 3b. The harmonic zero-point enthalpies are roughly 30 times larger than the anharmonic corrections. However, the differences between the harmonic zero-point energies of the ve phases considered at xed pressure are similar in magnitude to the differences between the anharmonic corrections, both being about 10 meV per atom, as shown in Fig. 3c,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.

Structural phase transitions. Figure 4 shows the two structural phase transitions that we have determined in this work, and our theoretical temperaturepressure phase diagram for solid molecular H is shown in Fig. 5. At 0 K, we nd a transition from P21/c-24 to C2/c-24 at around 23510 GPa. The corresponding transition pressure for D is 13 GPa higher (note that the difference between H and D is purely due to the DFT vibrational free energy and hence the difference in transition pressures between H and D is relatively precise). Our transition pressure is around 75 GPa greater than those observed experimentally for the transition between phases II and III, but the 13 GPa difference between the transition pressures for H and D agrees well with the experimentally measured value23. We note that the theoretical transition pressures between H and D would only differ by around 6 GPa without the inclusion of anharmonic effects.

As shown in Fig. 4, we also nd a temperature-driven transition from C2/c-24 to Pc-48 at pressures above 250 GPa

500

Pc- 48 Models phase IV

C 2/c-24

Models phase III

Phase I

P 21/c-24 Models phase II

400

Temperature (K)

300

200

100

Figure 5 | Theoretical temperaturepressure phase diagram for H. The solid black lines show the phase transitions calculated in this work, that is, the set of points at which the relative Gibbs free energy of two phases is zero. The dotted lines show the set of points at which the relative Gibbs free energy is one error bar from zero, and hence indicate the uncertainty in the phase boundaries. At pressures in excess of 350375 GPa the Gibbs free energies of the C2/c-24 and Pc-48 structures are within error bars of each other. The grey region indicates the temperaturepressure conditions under which phase I is found to exist in experiments.

4 NATURE COMMUNICATIONS | 6:7794 | DOI: 10.1038/ncomms8794 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.

NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8794 ARTICLE

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/c24 to Pc48 transition temperature as 32020 K.

DiscussionOur theoretical H phase diagram is in reasonable quantitative agreement with experiment, indicating that the P21/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 signicant disagreement of about 75 GPa between the experimental and theoretical phase IIIII transition pressure at 0 K. There are several possible reasons for our substantially larger phase IIIII transition pressure. First, 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 infrared spectra of the observed phase III. If a signicantly 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 IIIII transition pressure could be that we have neglected a signicant contribution to the energy of P21/c-24 (our model for phase II).

In particular, our calculations do not account for nuclear exchange effects, which are known to have a signicant effect on the phase I-II transition pressure22. However, 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 IIIII transition pressure in both H and D. Our analysis of different nite-size corrections in Supplementary Note 2 (see also the accompanying Supplementary Figs 3 and 4) indicates that nite-size effects in our relative enthalpies are well-controlled, but it is always possible that nite-size effects may be larger than anticipated. Finally, the xed-node approximation is an uncontrolled source of error in our DMC calculations and, although xed-node errors should largely cancel when relative energies are calculated, it cannot be ruled out that xed-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,29 because its Raman spectrum agrees well with the experimental one and because its weakly bonded layers lead to soft vibrational modes that thermally stabilize it. However, DFT static-lattice calculations together with the harmonic approximation for nuclear motion (used in refs 28,29) predict the Cmca-4 structure to be energetically favoured at all temperatures in the relevant pressure range (note that Cmca-4 is stabilised signicantly by harmonic zero-point 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 destabilizes Cmca-4, and we nd 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 nd that Cmca-12 is unstable at the pressures and temperatures studied in this work. We have found an important discrepancy between our calculated phase IIIII 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 method37,45 is capable of delivering much higher accuracy than DFT, and the scaling of the computational cost with system size enables the simulation 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 DFTPBE at a given external pressure. In DMC, the ground-state component of a trial wave function is projected out by simulating the Schrdinger equation in imaginary time, subject to the constraint that the nodal surface of the wave function is xed to be that of the trial wave function37,45. We used SlaterJastrow wave functions as implemented in the CASINO code46. Full technical details of our calculations can be found in Supplementary Note 2. The single-particle orbitals were obtained from the CASTEP code47 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 exible Jastrow factor48 whose parameters were optimized using variational Monte Carlo49. Variational Monte Carlo and DMC simulations were performed using 96 and 768 atoms, and the results were extrapolated to innite system size. Using the resources of the Oak Ridge Leadership Computing Facility, we achieved statistical error bars of o0.3 meV per atom in all our DMC calculations.

Anharmonic vibrational calculations. We have calculated harmonic vibrational free energies by using the nite-displacement method to construct the matrix of force constants and diagonalising the corresponding dynamical matrices over a ne 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 self-consistent eld method40,41,50, sampling the low-energy part of the DFTPBE BornOppenheimer energy surface along harmonic normal modes to large amplitudes. The resulting anharmonic Schrdinger equation for the nuclear motion was solved by expanding the wave function in a basis of simple harmonic oscillator eigenstates. Thermal occupation of excited states allowed us to calculate free energies at arbitrary temperatures. The vibrational free energy differences between the structures were converged to better than 1 meV per atom. Our approach does not describe possible melting.

Data availability. All relevant data present in this publication can be accessed at: http://www.repository.cam.ac.uk/handle/1810/248864

Web End =http://www.repository.cam.ac.uk/handle/1810/248864 .

References

1. Mao, H.-K. & Hemley, R. J. Ultrahigh-pressure transitions in solid hydrogen. Rev. Mod. Phys. 66, 671692 (1994).

2. McMahon, J. M., Morales, M. A., Pierleoni, C. & Ceperley, D. M. The properties of hydrogen and helium under extreme conditions. Rev. Mod. Phys. 84, 16071653 (2012).

3. Deemyad, S. & Silvera, I. F. Melting line of hydrogen at high pressures. Phys. Rev. Lett. 100, 155701 (2008).

4. Bonev, S. A., Schwegler, E., Ogitsu, T. & Galli, G. A quantum uid of metallic hydrogen suggested by rst-principles calculations. Nature 431, 669672 (2004).

5. Chen, J. et al. Quantum simulation of low-temperature metallic liquid hydrogen. Nat. Commun. 4, 2064 (2013).

6. Babaev, E., Sudb, A. & Ashcroft, N. W. A superconductor to superuid phase transition in liquid metallic hydrogen. Nature 431, 666668 (2004).

7. Babaev, E., Sudb, A. & Ashcroft, N. W. Observability of a projected new state of matter: a metallic superuid. Phys. Rev. Lett. 95, 105301 (2005).

NATURE COMMUNICATIONS | 6:7794 | DOI: 10.1038/ncomms8794 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications 5

& 2015 Macmillan Publishers Limited. All rights reserved.

ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/ncomms8794

8. Ashcroft, N. W. Metallic hydrogen: a high-temperature superconductor? Phys. Rev. Lett. 21, 1748 (1968).

9. Cudazzo, P. et al. Ab initio description of high-temperature superconductivity in dense molecular hydrogen. Phys. Rev. Lett. 100, 257001 (2008).

10. McMahon, J. M. & Ceperley, D. M. High-temperature superconductivity in atomic metallic hydrogen. Phys. Rev. B 84, 144515 (2011).

11. Loubeyre, P., Occelli, F. & LeToullec, R. Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen. Nature 416, 613617 (2002).12. Eremets, M. I. & Troyan, I. A. Conductive dense hydrogen. Nat. Mater. 10, 927931 (2011).

13. Howie, R. T., Guillaume, C. L., Scheler, T., Goncharov, A. F. & Gregoryanz, E. Mixed molecular and atomic phase of dense hydrogen. Phys. Rev. Lett. 108, 125501 (2012).

14. Howie, R. T., Scheler, T., Guillaume, C. L. & Gregoryanz, E. Proton tunneling in phase IV of hydrogen and deuterium. Phys. Rev. B 86, 214104 (2012).

15. Zha, C.-S., Liu, Z. & Hemley, R. J. Synchrotron infrared measurements of dense hydrogen to 360 GPa. Phys. Rev. Lett. 108, 146402 (2012).

16. Loubeyre, P., Occelli, F. & Dumas, P. Hydrogen phase IV revisited via synchrotron infrared measurements in H2 and D2 up to 290 GPa at 296K. Phys.

Rev. B 87, 134101 (2013).17. Zha, C.-S., Liu, Z., Ahart, M., Boehler, R. & Hemley, R. J. High-pressure measurements of hydrogen phase IV using synchrotron infrared spectroscopy. Phys. Rev. Lett. 110, 217402 (2013).

18. Goncharov, A. F., Howie, R. T. & Gregoryanz, E. Hydrogen at extreme pressures (review article). Low Temp. Phys. 39, 402408 (2013).

19. Zha, C.-S., Cohen, R. E., Mao, H.-K. & Hemley, R. J. Raman measurements of phase transitions in dense solid hydrogen and deuterium to 325 GPa. Proc. Natl Acad. Sci. USA 111, 47924797 (2014).

20. Cui, L., Chen, N. H., Jeon, S. J. & Silvera, I. F. Megabar pressure triple point in solid deuterium. Phys. Rev. Lett. 72, 3048 (1994).

21. Moshary, F., Chen, N. H. & Silvera, I. F. Remarkable high pressure phase line of orientational order in solid hydrogen deuteride. Phys. Rev. Lett. 71, 3814 (1993).

22. Mazin, I. I., Hemley, R. J., Goncharov, A. F., Hanand, M. & Mao, H.-K. Quantum and classical orientational ordering in solid hydrogen. Phys. Rev. Lett. 78, 1066 (1997).

23. Goncharov, A. F., Hemley, R. J. & Mao, H.-K. Vibron frequencies of solid H2 and D2 to 200 GPa and implications for the P-T phase diagram. J. Chem. Phys. 134, 174501 (2011).

24. Goncharenko, I. & Loubeyre, P. Neutron and x-ray diffraction study of the broken symmetry phase transition in solid deuterium. Nature 435, 12061209 (2005).

25. Akahama, Y. et al. Evidence from x-ray diffraction of orientational ordering in phase III of solid hydrogen at pressures up to 183 GPa. Phys. Rev. B 82, 060101 (2010).

26. Pickard, C. J. & Needs, R. J. Structure of phase III of solid hydrogen. Nat. Phys.

3, 473476 (2007).

27. Pickard, C. J. & Needs, R. J. Structures at high pressure from random searching. Phys. Status Solidi B 246, 536540 (2009).

28. Pickard, C. J., Martinez-Canales, M. & Needs, R. J. Density functional theory study of phase IV of solid hydrogen. Phys. Rev. B 85, 214114 (2012).

29. Pickard, C. J., Martinez-Canales, M. & Needs, R. J. Erratum: Density functional theory study of phase IV of solid hydrogen. Phys. Rev. B 86, 059902 (2012).

30. McMahon, J. M. & Ceperley, D. M. Ground-state structures of atomic metallic hydrogen. Phys. Rev. Lett. 106, 165302 (2011).

31. Liu, H., Zhu, L., Cui, W. & Ma, Y. Room-temperature structures of solid hydrogen at high pressures. J. Chem. Phys. 137, 074501 (2012).

32. Liu, H., Wang, H. & Ma, Y. Quasi-molecular and atomic phases of dense solid hydrogen. J. Phys. Chem. C 116, 92219226 (2012).

33. Johnson, K. A. & Ashcroft, N. W. Structure and bandgap closure in dense hydrogen. Nature 403, 632635 (2000).

34. Azadi, S. & Foulkes, W. M. C. Fate of density functional theory in the study of high-pressure solid hydrogen. Phys. Rev. B 88, 014115 (2013).

35. Morales, M. A., McMahon, J. M., Pierleoni, C. & Ceperley, D. M. Towards a predictive rst-principles description of solid molecular hydrogen with density functional theory. Phys. Rev. B 87, 184107 (2013).

36. Clay, R. C. et al. Benchmarking exchange-correlation functionals for hydrogen at high pressures using quantum Monte Carlo. Phys. Rev. B 89, 184106 (2014).

37. Foulkes, W. M. C., Mitas, L., Needs, R. J. & Rajagopal, G. Quantum Monte Carlo simulations of solids. Rev. Mod. Phys. 73, 3383 (2001).

38. Natoli, V., Martin, R. M. & Ceperley, D. Crystal structure of molecular hydrogen at high pressure. Phys. Rev. Lett. 74, 16011604 (1995).

39. Azadi, S., Foulkes, W. M. C. & Khne, T. D. Quantum Monte Carlo study of high pressure solid molecular hydrogen. N. J. Phys. 15, 113005 (2013).

40. Azadi, S., Monserrat, B., Foulkes, W. M. C. & Needs, R. J. Dissociation of high-pressure solid molecular hydrogen: a quantum Monte Carlo and anharmonic vibrational study. Phys. Rev. Lett. 112, 165501 (2014).

41. Monserrat, B., Drummond, N. D. & Needs, R. J. Anharmonic vibrational properties in periodic systems: energy, electron-phonon coupling, and stress. Phys. Rev. B 87, 144302 (2013).

42. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 38653868 (1996).

43. Lee, C., Yang, W. & Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 37, 785789 (1988).

44. Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 38, 30983100 (1988).

45. Ceperley, D. M. & Alder, B. J. Ground state of the electron gas by a stochastic method. Phys. Rev. Lett. 45, 566 (1980).

46. Needs, R. J., Towler, M. D., Drummond, N. D. & Lpez Ros, P. Continuum variational and diffusion quantum Monte Carlo calculations. J. Phys. Condens. Matter 22, 023201 (2010).

47. Clark, S. J. et al. First principles methods using CASTEP. Z. Kristallogr. 220,

567570 (2005), 5/6/2005.

48. Drummond, N. D., Towler, M. D. & Needs, R. J. Jastrow correlation factor for atoms, molecules, and solids. Phys. Rev. B 70, 235119 (2004).

49. Umrigar, C. J., Toulouse, J., Filippi, C., Sorella, S. & Hennig, R. G. Alleviation of the fermion-sign problem by optimization of many-body wave functions. Phys. Rev. Lett. 98, 110201 (2007).

50. Monserrat, B., Drummond, N. D., Pickard, C. J. & Needs, R. J. Electron-phonon coupling and the metallization of solid helium at terapascal pressures. Phys. Rev. Lett. 112, 055504 (2014).

Acknowledgements

We thank Dominik Jochym for help with the implementation of the BLYP density functional. Financial support was provided by the Engineering and Physical Sciences Research Council (EPSRC), UK. B.M. thanks Robinson College, Cambridge and the Cambridge Philosophical Society for a Henslow Research Fellowship. This research used resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Ofce of Science of the US. Department of Energy under Contract No. DE-AC05-00OR22725. Additional calculations were performed on the Cambridge High Performance Computing Service facility Darwin and the N8 high-performance computing facility provided and funded by the N8 consortium and EPSRC (Grant No. EP/K000225/1).

Author contributions

R.J.N. conceived and managed the project, C.J.P. implemented the BLYP functional in the castep code for our work, P.L.R. optimized the trial wave functions for our DMC calculations, J.H. L.-W. carried out the DFT calculations and calculated the DFT and DMC equations of state, B.M. carried out the anharmonic vibrational calculations, and N.D.D. performed the DMC calculations. All authors contributed to the manuscript.

Additional information

Supplementary information accompanies this paper at http://www.nature.com/naturecommunications

Web End =http://www.nature.com/ http://www.nature.com/naturecommunications

Web End =naturecommunications

Competing nancial interests: Chris Pickard is an author of the CASTEP code, which is marketed commercially by Accelrys, and receives royalties on the sales. The remaining authors declare no competing nancial interests.

Reprints and permission information is available online at http://npg.nature.com/reprintsandpermissions

Web End =http://npg.nature.com/ http://npg.nature.com/reprintsandpermissions

Web End =reprintsandpermissions/

How to cite this article: Drummond, N. D. et al. Quantum Monte Carlo study of the phase diagram of solid molecular hydrogen at extreme pressures. Nat. Commun. 6:7794 doi: 10.1038/ncomms8794 (2015).

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the articles Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/3.0/

Web End =http://creativecommons.org/licenses/by/4.0/

6 NATURE COMMUNICATIONS | 6:7794 | DOI: 10.1038/ncomms8794 | http://www.nature.com/naturecommunications

Web End =www.nature.com/naturecommunications

& 2015 Macmillan Publishers Limited. All rights reserved.