A comparative study using state-of-the-art electronic structure theories on solid hydrogen phases under high pressures
AARTICLE
OPEN
A comparative study using state-of-the-art electronic structuretheories on solid hydrogen phases under high pressures
Ke Liao *, Xin-Zheng Li , Ali Alavi and Andreas Grüneis * Identifying the atomic structure and properties of solid hydrogen under high pressures is a long-standing problem of high-pressurephysics with far-reaching signi fi cance in planetary and materials science. Determining the pressure-temperature phase diagram ofhydrogen is challenging for experiment and theory due to the extreme conditions and the required accuracy in the quantummechanical treatment of the constituent electrons and nuclei, respectively. Here, we demonstrate explicitly that coupled clustertheory can serve as a computationally ef fi cient theoretical tool to predict solid hydrogen phases with high accuracy. We present a fi rst principles study of solid hydrogen phases at pressures ranging from 100 to 450 GPa. The computed static lattice enthalpies arecompared to state-of-the-art diffusion Monte Carlo results and density functional theory calculations. Our coupled cluster theoryresults for the most stable phases including C2/c-24 and P2 /c-24 are in good agreement with those obtained using diffusionMonte Carlo, with the exception of Cmca-4, which is predicted to be signi fi cantly less stable. We discuss the scope of the employedmethods and how they can contribute as ef fi cient and complementary theoretical tools to solve the long-standing puzzle ofunderstanding solid hydrogen phases at high pressures. npj Computational Materials (2019) 5:110 ; https://doi.org/10.1038/s41524-019-0243-7 INTRODUCTION
Hydrogen is the lightest and most abundant element in theUniverse, yet its phase diagram at high pressures and lowtemperatures remains elusive. Due to the subtle interplay ofquantum nuclear and electronic correlation effects, – thequestion as to which state of matter is stable at high pressuresis controversial. Likely, candidates for high-pressure phasesinclude various orientationally ordered molecular crystals, – (liquid) metallic, – superconducting and super fl uid systems. These potentially exotic states of matter and their crucialimportance for astrophysical, planetary as well as materialssciences has led to intensi fi ed investigations using both experi-mental and theoretical techniques. However, currently availablecalculated as well as measured equilibrium phase boundaries varystrongly with respect to the employed methods and suffer partlyfrom uncontrolled sources of error.Experiments that seek to determine properties of hydrogenunder high pressures are hindered by various problems; forexample, the low X-ray scattering cross section of hydrogen, thesmall sample sizes and the diffusive nature of hydrogen. Recentclaims of experimentally measured metallic phases are there-fore under debate, whilst earlier experimental results havenot been able to conclusively detect metallic behaviour up to apressure of 320 –
342 GPa.Determining the Wigner – Huntington transition using theore-tical methods is extremely challenging. Despite the signi fi cantadvancements of modern ab initio theories in the past decades,the predicted metallisation pressure varies signi fi cantly in a rangeof ~150 –
450 GPa, depending on the employed method. – Most ab initio studies of solid hydrogen are based either ondensity functional theory (DFT) or quantum Monte Carlocalculations. – DFT is considered the workhorse method in computational materials science, and can be used to calculatelattice enthalpies on the level of various approximate exchangeand correlation (XC) energy functionals. Furthermore, theHellmann – Feynman theorem provides access to atomic forcesand allows for optimising structures, as well as calculatingphonons on the level of DFT. Calculated and measured infraredand Raman spectra serve as a reliable tool for a direct comparisonbetween theory and experiment. – – However, differentparameterisations of the XC functional in DFT give inconsistentpredictions, e.g., PBE predicts a too low metallisation pressurecompared with experiments, whilst other exchange functionalsproduce higher pressures than DMC.
Instead, more accurate methods including diffusion Monte Carlo(DMC) have been employed to predict more reliable pressuretemperature-phase diagrams, – which correct the underesti-mation of the metallisation pressure by DFT-PBE to a large extent.However, DMC calculations rely on the fi xed-node approximation,and most of the current studies use crystal structures optimisedusing DFT. A critical assessment of the errors introduced by theseapproximations is still missing in literature and requires computa-tionally ef fi cient and concomitantly accurate methods.In this work, we show that quantum chemical wavefunctiontheories hold the promise to serve as an ef fi cient and accuratetool for the investigation of high-pressure phases of solidhydrogen. In particular, we fi nd that coupled cluster theory achieves a good trade-off between computational cost andaccuracy when employing recently developed techniques thatallow for simulating the thermodynamic limit of periodic systemsin an ef fi cient manner. We note that these fi nite sizecorrections have paved the way for a number of ab initio studies,including predictions of molecule – surface interactions – andpressure – temperature phase diagrams of carbon and boron Max Planck Institute for Solid State Research, Heisenbergstrasse 1, 70569 Stuttgart, Germany. School of Physics, State Key Laboratory for Arti fi cial Microstructure andMesoscopic Physics, Peking University, 100871 Beijing, P.R. China. Collaborative Innovation Centre of Quantum Matter, Peking University, 100871 Beijing, P.R. China. Department of Chemistry, University of Cambridge, Lens fi eld Road, Cambridge CB2 1EW, UK. Institute for Theoretical Physics, Vienna University of Technology, WiednerHauptstrasse 8-10, 1040 Vienna, Austria. *email: [email protected]; [email protected]
Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences () :,; itride allotropes. The studies referred to above have demon-strated that coupled cluster methods achieve a similar level ofaccuracy as DMC for solid state systems that are not stronglycorrelated. Moreover, coupled cluster methods have beenbenchmarked against various more accurate methods in modelhydrogen systems, showing the high accuracy of the methods inweakly correlated situations. Furthermore, we employ full con fi g-uration interaction Quantum Monte Carlo (FCIQMC) – in thiswork for small systems to examine the validity of the coupledcluster method. RESULTS
We investigate theoretical results for the static lattice enthalpies ofsolid hydrogen phases computed on different levels of theory. Thestatic lattice enthalpy is de fi ned by H ¼ E þ PV ; (1)where P is the pressure estimated from the E (cid:2) V relation and V corresponds to the volume per atom. E refers to the total ground-state energy per atom obtained using DFT, HF or CC theory in theBorn – Oppenheimer approximation. In passing, we note that theimportance of quantum nuclear effects for transition pressures ofsolid hydrogen phases has been explored in refs – In this work,we will focus on the accuracy of the employed electronic structuretheories only, disregarding such contributions. The coupled clustersingles and doubles (CCSD) energy is de fi ned as the sum of theHartree – Fock and the corresponding electronic correlationenergy. We discuss the convergence of the electronic exchangeand correlation energy contributions with respect to theemployed k -mesh used to sample the fi rst Brillouin zone andthe one-electron basis set, as well as additional computationaldetails in the supplementary information. The pressure – volumerelation of each phase, P ð V Þ ¼ (cid:2) d E d V , is obtained in the followingmanner. The total energy retrieved as a function of the volume peratom, E ð V Þ , is fi tted with a polynomial function of V (cid:2) in anoptimal order that minimises the fi tting residual and providessmooth curves. We fi nd that a third-order polynomial fi tting isadequate for all phases, except for phase P2 /c-24 which is fi ttedusing a fourth-order polynomial. A further increase in the fi ttingorder can result in arti fi cial wiggling behaviours of the H ð P Þ curves. The derivative with respect to the volume is readilyobtained in an analytic manner using the fi tted E ð V Þ function,yielding smooth P ð V Þ curves. We present all static latticeenthalpies relative to the C2/c-24 phase unless stated otherwise.In total, we study fi ve solid hydrogen phases: Cmca-4 (Cmca-Low),Cmca-12, C2/c-24, P2 /c-24 and P6 /m-16, where we haveadopted the convention of naming the structures by theirsymmetries followed by the number of atoms in the primitivecells. Phase Cmca-4, Cmca-12 and C2/c-24 consist of layeredhydrogen molecules whose bonds lie within the plane of the layer,forming distorted hexagonal shapes. Whereas some bonds ofhydrogen molecules in phase P6 /m-16 lie perpendicularly to theplane of the layer. P2 /c-24 consists of molecules arranged on adistorted hexagonal close-packed lattice.These structures have previously been selected as potentialcandidates as the most stable high-pressure phases of hydrogen and have been studied by DMC methods. We notice that a familyof ‘ mixed ’ structures are also identi fi ed as promising candidates inref., however, for the current comparative studies among CCSD,DFT-PBE and DMC, they are not included here but could be aninteresting topic for future work.We have optimised the geometries of the structures employingthe DFT-PBE functional. The DFT calculations have beenperformed using the Vienna ab initio simulation package (VASP)employing a plane wave basis set in the framework of theprojector augmented wave method. More details about thestructures can be found in refs. We fi rst discuss results of the investigated high-pressure phaseson the level of DFT. Figure 1 depicts the DFT-PBE static latticeenthalpies relative to the C2/c-24 phase. DFT-PBE predicts the C2/c-24 phase as the most stable phase at pressures ranging from~100 –
290 GPa. In a small range of pressures ~300 GPa, the Cmca-12 phase is found to be thermodynamically stable, whereas themetallic Cmca-4 phase becomes stable at pressures exceeding~330 GPa. Experimentally, no metallic phases have been observedin this pressure range, and quantum nuclear effects do notaccount for this discrepancy either. The too low metallisationpressure can be attributed to the lacking of van der Waalsinteractions in PBE functional, resulting in underestimation ofthe stability in the molecular structures. We note that Fig. 1 alsodepicts static lattice enthalpies from ref. obtained using DFT-PBE.We attribute the minor differences between the static latticeenthalpies to small differences in the employed structures and the fi tting procedure that is employed to compute the latticeenthalpies from the total energies retrieved as a function of thevolume per atom. We stress that the computed enthalpies are verysensitive to the employed structures.In contrast to approximate XC functionals employed in DFTcalculations, quantum chemical many-electron methods allow forapproximating the electronic XC energy in a more systematicmanner, albeit at signi fi cantly larger computational cost. Thesimplest wavefunction-based method is the Hartree – Fock (HF)approximation that neglects electronic correlation effects byde fi nition, employing a single Slater determinant as Ansatz forthe electronic wavefunction. Figure 2 depicts the static latticeenthalpies computed in the HF approximation relative to C2/c-24.In contrast to DFT-PBE results, we fi nd that HF theory signi fi cantlyreduces the stability of the Cmca-4 and Cmca-12 phases, shiftingtheir transition pressures far above 400 GPa. However, the HFmethod is not a good approximation for metallic systems, despitethe fact that it is free from self-interaction errors. In particular, HFband gaps are usually signi fi cantly overestimated compared withthe experiment. Moreover, the lack of electronic correlation in theHF Ansatz leads to the neglect of van der Waals contributions thatare crucial for a correct description of relative stabilities ofmolecular crystals. We note that van der Waals contributions tothe binding energy of molecular crystals become in general largerfor smaller volumes due to the polynomial decay of the dispersioninteraction with respect to the intermolecular distance. Due to thereasons outlined above, the static lattice enthalpies calculated onthe level of HF theory are expected to exhibit signi fi cant errors Fig. 1
DFT-PBE relative enthalpies. The DFT-PBE relative enthalpiesof structures that are used for CCSD calculations in this paper(dashed lines) and that of the structures from ref. (full lines). DFTfavours the atomic phase Cmca-4 at high pressures K. Liao et al. npj Computational Materials (2019) 110npj Computational Materials (2019) 110
DFT-PBE relative enthalpies. The DFT-PBE relative enthalpiesof structures that are used for CCSD calculations in this paper(dashed lines) and that of the structures from ref. (full lines). DFTfavours the atomic phase Cmca-4 at high pressures K. Liao et al. npj Computational Materials (2019) 110npj Computational Materials (2019) 110 Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences () :,; ompared with more accurate electronic structure theories andwill serve as a reference for post-HF methods only.Here, we employ the coupled cluster singles and doubles(CCSD) method to account for electronic correlation effects usinga HF reference. Periodic CCSD theory results for the static latticeenthalpies relative to the C2/c-24 phase are shown in Fig. 3.Compared with the HF theory, CCSD stabilises the Cmca-4 phaseby ~40 meV/atom at pressures above 300 GPa. Similarly, therelative static lattice enthalpy of Cmca-12 is lowered by ~20 meV/ atom in CCSD compared with HF. For the P2 /c-24 and P6 /m-16phases, we observe an opposite effect of the CCSD correlationenergy contribution, reducing their stability relative to C2/c-24 by~30 meV/atom at pressures exceeding 250 GPa. We note thatCCSD theory reduces the differences in the relative static latticeenthalpies of the considered phases compared with the HFapproximation.The CCSD energy is the sum of the HF energy and anapproximation to the electronic correlation energy that iscomputed using an exponential Ansatz for the wavefunction.Due to the many-electron nature of the employed Ansatz, CCSDtheory is exact for two-electron systems. The coupling betweenelectron pairs is, however, approximated by truncating the many-body perturbation expansion in a computationally ef fi cientmanner and performing a resummation to in fi nite order of certaincontributions only. As a consequence, CCSD theory is expectedto yield highly accurate results for the molecular hydrogencrystals. This is con fi rmed by comparing with the correspondingDMC results from ref. for C2/c-24 and P2 /c-24 depicted in Fig. 3that agree very well with our CCSD fi ndings. Furthermore, staticlattice enthalpies obtained on the level of CCSD and DMC (onlyshown in ref. ) for Cmca-12 relative to C2/c-24 are in goodagreement as well and the transition pressure between P6 /m-16and C2/c-24 by CCSD ( (cid:3)
350 GPa) and DMC ( (cid:3) –
350 GPa onlyshown in ref. ) are in reasonable agreement. However, we notethat the DMC and CCSD results differ by ~40 meV/atom for therelative static lattice enthalpy of the Cmca-4 phase. In particular,the difference of the static lattice enthalpies of Cmca-4 and C2/c-24 at 350 GPa are ~100 meV/atom, 60 meV/atom and 20 meV/atom using HF, CCSD and DMC, respectively. DISCUSSION
We now discuss possible reasons for the discrepancy betweenDMC and CCSD results for the Cmca-4 phase. DMC calculationsemploy the fi xed-node approximation, whereas CCSD theorytruncates the particle – hole excitation operator in the exponent ofthe wavefunction Ansatz. Fixed-node DMC gives the upperbounds to the total energies of each phase. However, it is notnecessarily the case that the lower enthalpy difference betweenCmca-4 and C2/c-24 predicted by DMC is more reliable than thatby CCSD, since the fi xed-node errors in each phase do notnecessarily cancel out accurately. The fi xed-node errors in the totalDMC energy can be estimated using back fl ow transformationsand by comparing with full con fi guration interaction quantumMonte Carlo (FCIQMC) results for the uniform electron gas. It has been shown that the fi xed-node errors are ~1 mHa perelectron (27.2 meV/electron) in the high-density regime. In thecase of solid hydrogen, the authors of ref. report in theirSupplementary Material that the energy in phase C2/c-24 islowered by 1 mHa/atom (27.2 meV/atom) when employing back- fl ow transformations and ref. reports that for Cmca-4 theback fl ow transformations lower the energy by only 10 meV/atom.This indicates that back fl ow transformations can depend sig-ni fi cantly on the phases. Even though a large part of the fi xed-node errors are expected to cancel when the energy differencebetween phases is computed, the remaining errors can still be onthe scale of 10 meV/atom. On the other hand, we stress that thechange from HF to CCSD relative static lattice enthalpies is on thescale of 40 meV/atom, indicating that a better approximation tothe many-electron wavefunction than employed by CCSD theorycould be necessary to achieve the required level of accuracy. Wehave also performed calculations using higher level theories,including FCIQMC, for smaller supercells containing 24 atoms atvolumes corresponding to a DFT pressure of 400 GPa. These fi ndings indicate that post-CCSD corrections to static latticeenthalpy differences for Cmca-4 and P2 /c-24 are expected tobe ~10 meV/atom. In short, both DMC and CCSD rely on good Fig. 2
HF relative enthalpies. The HF relative enthalpies of structuresthat are used in this paper. In contrast to the DFT-PBE result, theatomic phase Cmca-4 is unfavoured at high pressures
Fig. 3
CCSD relative enthalpies. The CCSD relative enthalpies ofstructures that are used in this paper (dashed lines) and the DMCrelative enthalpies of structures from ref. (full lines). The thicknessof the full lines refer to the standard deviations of stochasticsampling of the 1st Brillouin zone while performing twist-averagingin the DMC calculations. In this work, the 1st Brillouin zone issampled using a dense regular grid such that the errors areconverged to within 1 meV/atom. See supplementary informationfor details. CCSD and DMC agree very well in the most stablemolecular phases, i.e., C2/c-24, P2 /c-24 and Cmca-12, whilst theonly discrepancy exists in the Cmca-4 phase, which is predicted byDFT-PBE to be metallic at high pressures. The phase transitionbetween P6 /m-16 and C2/c-24 predicted by CCSD happens at~350 GPa, which agrees reasonably well with the DMC transitionpressure range 250 –
350 GPa from ref. K. Liao et al. Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences npj Computational Materials (2019) 110 ancellations in errors introduced by their respective approxima-tions to produce accurate predictions, especially when phases ofdifferent physical natures are compared. In addition to theinherent errors of DMC and CCSD theory, fi nite size and basisset errors can also be signi fi cant. The latter only applies to CCSDcalculations and has been checked carefully as outlined in thesupplementary information. As regards the fi nite size error, westudy supercells containing 96 atoms and employ twist averagingas well as structure factor interpolation methods for our CCSDcalculations to achieve a level of precision that is comparable withDMC results. Despite the above considerations, we can currentlynot draw any fi rm conclusion about the reason for the discrepancybetween DMC and CCSD results for Cmca-4. However, we notethat recently developed basis set convergence accelerationtechniques will enable future studies of bigger systems usingCCSD and FCIQMC theory that can hopefully provide moreinsight.Despite the discrepancy between CCSD and DMC fi ndings forCmca-4, we point out that the good agreement for the staticlattice enthalpies of the most stable high-pressure hydrogenphases is encouraging. Achieving accurate thermodynamic limitresults for such systems on the level of CCSD theory has onlybecome possible recently due to the development of thecorresponding fi nite size corrections as outlined in refs Furthermore, we note that the computational cost of thecorresponding CCSD calculations is still moderate compared withmethods with a similar accuracy. A single CCSD ground-stateenergy calculation for a system containing 96 atoms using 400bands requires ~250 CPU hours, implying that it will becomepossible in the near future to perform structural relaxation of theemployed crystal structures rather than relying on structuresoptimised using DFT-PBE. This is necessary for truly reliablepredictions of high-pressure phases of solid hydrogen.We have presented static lattice enthalpies for high-pressurephases of solid hydrogen calculated using state-of-the-artelectronic structure methods, including coupled cluster theory.We fi nd that CCSD theory results agree well with DMC fi ndingsfrom ref.: phase C2/c-24 becomes more stable than phase P2 /c-24 at ~250 GPa; phase Cmca-4 and Cmca-12 are less stable thanphase C2/c-24 in the pressure range from 100 GPa to 400 GPa. Theonly discrepancy between CCSD and DMC is found for the Cmca-4phase, and we have discussed possible sources of error. Futurework will include the effects of the nuclei motions which arecrucial in making theoretical predictions comparable with experi-ments. Based on the presented fi ndings, the required computa-tional cost of the employed CCSD implementation and recentmethodological advancements, we conclude that prospectiveCCSD studies will make it possible to optimise structures of solidhydrogen phases at high pressures with DMC accuracy. This willenable complementary CCSD and DMC studies with a signi fi cantlyimproved level of accuracy and achieve unprecedented physicalinsight into the Wigner – Huntington transition of solid hydrogen.
METHODS
The CCSD calculations have been performed employing the coupledcluster for solids (CC4S ) code interfaced to the Vienna ab initiosimulation package (VASP). The projector augmented wavemethod, as implemented in VASP, is used for all calcula-tions. This section provides an overview of the computationalmethods and convergence techniques employed in this work. Formore details, we refer the reader to the corresponding sections inthe supplementary information under the same section titles.Geometries
The structures have been optimised using DFT-PBE and are similar to thoseemployed in ref. The forces on the atoms of the optimised structures arenot larger than 0.1 eV/ Å . With hindsight, it would have been preferable to use exactly the same structures as published in ref. However, for thepurpose of this work, the agreement between the structures suf fi ces. Forthe CCSD calculations, we employ supercells containing up to 96 atomsthat are as isotropic as possible and are obtained using the same methodas described in the supplementary Note 2 of ref. In this manner, fi nite sizeerrors can be signi fi cantly reduced. CCSD basis set convergence
For the equilibrium phase boundaries in the pressure – temperature phasediagram, only relative enthalpies are relevant. Therefore, we haveconverged the energy differences with respect to the basis set only.MP2 natural orbitals (MP2NOs) provide faster convergence thancanonical Hartree – Fock orbitals (HFOs) computed from the full planewave basis set. The convergence tests of the CCSD correlation energydifferences with respect to the number of orbitals per atom relative tophase C2/c-24 have been carried out using supercells containing 24 atomsfor all phases, except for phase P6 /m-16 which contains 16 atoms in thesupercells. The results are summarised in Table I and Fig. 1 in theSupplementary Information File. We note that the basis set incompletenesserrors are mainly due to the electronic cusp conditions, which are verylocal effects and are not dependent on the supercell size. Hartree – Fock fi nite size convergence The HF energies are converged to within 1 meV/atom using increasinglylarge supercells or dense k -meshes sampling the fi rst Brillouin zone. Therequired system sizes for all phases are summarised in Table II in theSupplementary Information. CCSD fi nite size convergence The twist-averaging technique and fi nite size corrections, based onthe interpolation of the transition structure factor, are applied on 96-atomsupercells to approximate the thermodynamic limit of the CCSDcorrelation energies. Post-CCSD error estimates
We applied some higher level theories, including DCSD,
CCSD(T) and FCIQMC, to estimate the post-CCSD error. The results are summarisedin Fig. 2, Fig. 3, Fig. 4 and Table III in the Supplementary Information File.
DATA AVAILABILITY
The data that support the fi ndings of this study are available from the correspondingauthors upon reasonable request. CODE AVAILABILITY
The VASP code is a copyrighted software and can be obtained fromits of fi cial website. The CC4S code is available from A.G. uponreasonable request and will be made open-source in the future. Received: 14 May 2019; Accepted: 3 October 2019;
REFERENCES
1. Li, X.-Z. et al. Classical and quantum ordering of protons in cold solid hydrogenunder megabar pressures.
J. Phys. Condense. Matter , 085402 (2013).2. Drummond, N. D. et al. Quantum Monte Carlo study of the phase diagram of solidmolecular hydrogen at extreme pressures. Nat. Commun. , 7794 (2015).3. Azadi, S., Foulkes, W. M. C. & Kühne, T. D. Quantum Monte Carlo study of highpressure solid molecular hydrogen. New J. Phys. , 113005 (2013).4. McMinis, J., Clay, R. C., Lee, D. & Morales, M. A. Molecular to atomic phase tran-sition in hydrogen under high pressure. Phys. Rev. Lett. , 105305 (2015).5. Morales, M. A., McMahon, J. M., Pierleoni, C. & Ceperley, D. M. Towards a pre-dictive fi rst-principles description of solid molecular hydrogen with densityfunctional theory. Phys. Rev. B , 184107 (2013).6. Morales, M. A., McMahon, J. M., Pierleoni, C. & Ceperley, D. M. Nuclear quantumeffects and nonlocal exchange-correlation functionals applied to liquid hydrogenat high pressure. Phys. Rev. Lett. , 065702 (2013).
K. Liao et al. npj Computational Materials (2019) 110npj Computational Materials (2019) 110
K. Liao et al. npj Computational Materials (2019) 110npj Computational Materials (2019) 110 Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences . Hemley, R. J. & Mao, H.-k. Phase transition in solid molecular hydrogen at ultra-high pressures.
Tech. Rep. , 857 –
860 (1988).8. Lorenzana, H. E., Silvera, I. F. & Goettel, K. A. Evidence for a structural phasetransition in solid hydrogen at megabar pressures.
Phys. Rev. Lett. , 2080 – Rev.Mod. Phys. , 671 –
692 (1994).10. Natoli, V., Martin, R. M. & Ceperley, D. Crystal structure of molecular hydrogen athigh pressure.
Phys. Rev. Lett. , 1601 – Nature , 1206 (2005).12. Pickard, C. J. & Needs, R. J. Structure of phase III of solid hydrogen.
Nat. Phys. ,473 –
476 (2007).13. Howie, R. T., Guillaume, C. L., Scheler, T., Goncharov, A. F. & Gregoryanz, E. Mixedmolecular and atomic phase of dense hydrogen.
Phys. Rev. Lett. , 125501(2012).14. Wigner, E. & Huntington, H. B. On the possibility of a metallic modi fi cation ofhydrogen. J. Chem. Phys. , 764 –
770 (1935).15. Hemley, R. J. & Mao, H.-k. Optical studies of hydrogen above 200 Gigapascals:evidence for metallization by band overlap.
Science , 1462 – fl uid molecular hydrogenat 140 GPa (1.4 Mbar). Phys. Rev. Lett. , 1860 – Nature , 632 –
635 (2000).18. Eremets, M. I. & Troyan, I. A. Conductive dense hydrogen.
Nat. Mater. , 927 – Nat. Commun. , 2064 (2013).20. Knudson, M. D. et al. Direct observation of an abrupt insulator-to-metal transitionin dense liquid deuterium. Science , 1455 – fi rst-order phase transition tometallic hydrogen. Phys. Rev. B , 155128 (2016).22. Dias, R. P. & Silvera, I. F. Erratum for the research article observation of thewigner-huntington transition to metallic hydrogen by R. P. Dias and I. F. Silvera(Science (2017) 355 (715) https://doi.org/10.1126/science.aal1579). Science ,eaao5843 (2017).23. Ashcroft, N. W. Metallic hydrogen: a high-temperature superconductor?
Phys. Rev.Lett. , 1748 – fl uid phasetransition in liquid metallic hydrogen. Nature , 666 (2004).25. Dias, R. P. & Silvera, I. F. Observation of the Wigner-Huntington transition tometallic hydrogen.
Science , 715 –
718 (2017).26. Geng, H. Y. Public debate on metallic hydrogen to boost high pressure research.
MRE , 275 –
277 (2017).27. Narayana, C., Luo, H., Orloff, J. & Ruoff, A. L. Solid hydrogen at 342 GPa: noevidence for an alkali metal.
Nature , 46 –
49 (1998).28. Loubeyre, P., Occelli, F. & LeToullec, R. Optical studies of solid hydrogen to 320GPa and evidence for black hydrogen.
Nature , 613 –
617 (2002).29. Chacham, H. & Louie, S. G. Metallization of solid hydrogen at megabar pressures:a fi rst-principles quasiparticle study. Phys. Rev. Lett. , 64 –
67 (1991).30. Tse, J. S. & Klug, D. D. Evidence from molecular dynamics simulations for non-metallic behaviour of solid hydrogen above 160 GPa.
Nature , 595 –
597 (1995).31. Azadi, S., Monserrat, B., Foulkes, W. M. C. & Needs, R. J. Dissociation of high-pressure solid molecular hydrogen: a quantum monte carlo and anharmonicvibrational study.
Phys. Rev. Lett. , 165501 (2014).32. Azadi, S., Drummond, N. D. & Foulkes, W. M. C. Nature of the metallizationtransition in solid hydrogen.
Phys. Rev. B , 35142 (2017).33. Kohanoff, J., Scandolo, S., Chiarotti, G. L. & Tosatti, E. Solid molecular hydrogen:the broken symmetry phase. Phys. Rev. Lett. , 2783 – Phys. Rev. B , 2092 – Phys. Rev. B - Condens. Matter Mater.Phys. , 184106 (2014).36. Zha, C.-S., Liu, Z. & Hemley, R. J. Synchrotron infrared measurements of densehydrogen to 360 GPa. Phys. Rev. Lett. , 146402 (2012).37. Han fl and, M., Hemley, R. J. & Mao, H.-k Novel infrared vibron absorption in solidhydrogen at megabar pressures. Phys. Rev. Lett. , 3760 – Phys. Rev. B , 134101 (2016).39. Zhang, X.-W., Wang, E.-G. & Li, X.-Z. Ab initio investigation on the experimentalobservation of metallic hydrogen. Phys. Rev. B , 134110 (2018).40. Azadi, S. & Ackland, G. J. The role of van der Waals and exchange interactions inhigh-pressure solid hydrogen. Phys. Chem. Chem. Phys. , 21829 – Phys. Scr. , 251 (1980). 42. Bartlett, R. J. & Musial, M. Coupled-cluster theory in quantum chemistry. Rev. Mod.Phys. , 291 –
352 (2007).43. Liao, K. & Grüneis, A. Communication: fi nite size correction in periodic coupledcluster theory calculations of solids. J. Chem. Phys. , 0 – Phys. Rev. X ,21043 (2018).45. Tsatsoulis, T. et al. A comparison between quantum chemistry and quantumMonte Carlo techniques for the adsorption of water on the (001) LiH surface. J.Chem. Phys. , 204108 (2017).46. Tsatsoulis, T., Sakong, S., Groß, A. & Grüneis, A. Reaction energetics of hydrogenon Si(100) surface: a periodic many-electron theory study.
J. Chem. Phys. ,244105 (2018).47. Brandenburg, J. G. et al. Physisorption of water on graphene: subchemicalaccuracy from many-body electronic structure methods.
J. Phys. Chem. Lett. ,358 –
368 (2019).48. Gruber, T. & Grüneis, A.
Ab initio calculations of carbon and boron nitride allo-tropes and their structural phase transitions using periodic coupled cluster the-ory.
Phys. Rev. B , 134108 (2018).49. Motta, M. et al. Towards the solution of the many-electron problem in realmaterials: equation of state of the hydrogen chain with state-of-the-art many-body methods. Phys. Rev. X , 031059 (2017).50. Booth, G. H., Thom, A. J. & Alavi, A. Fermion monte carlo without fi xed nodes: agame of life, death, and annihilation in Slater determinant space. J. Chem. Phys. , 054106 (2009).51. Cleland, D., Booth, G. H. & Alavi, A. Communications: survival of the fi ttest:accelerating convergence in full con fi guration-interaction quantum Monte Carlo. J. Chem. Phys. , 041103 (2010).52. Booth, G. H., Grüneis, A., Kresse, G. & Alavi, A. Towards an exact description ofelectronic wavefunctions in real solids.
Nature , 365 –
370 (2013).53. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation madesimple.
Phys. Rev. Lett. , 3865 – fi cient iterative schemes for ab initio total-energycalculations using a plane-wave basis set. Phys. Rev. B , 11169 – Phys. Rev. Lett. , 185701 (2011).56. Reynolds, P. J., Ceperley, D. M., Alder, B. J. & Lester, W. A. Fixed-node quantumMonte Carlo for molecules a)b) . J. Chem. Phys. , 5593 – fi guration interactionperspective on the homogeneous electron gas. Phys. Rev. B - Condens. MatterMater. Phys. , 81103 (2012).58. Ruggeri, M., Ríos, P. L. & Alavi, A. Correlation energies of the high-density spin-polarized electron gas to meV accuracy. Phys. Rev. B , 161105 (2018).59. Irmler, A., Hummel, F. & Grüneis, A. On the duality of ring and ladder diagramsand its importance for many-electron perturbation theories. arXiv:1903.05458[cond-mat, physics:physics] . http://arxiv.org/abs/1903.05458. (2019).60. Luo, H. & Alavi, A. Combining the transcorrelated method with full con fi gurationinteraction quantum Monte Carlo: application to the homogeneous electron gas. J. Chem. Theory Comput. , 1403 – B. Mater. Sci. , 33 –
41 (2003).62. Kresse, G. & Hafner, J. Norm-conserving and ultrasoft pseudopotentials for fi rst-row and transition elements. J. Phys.: Condens. Matter , 8245 – J. Chem. Theory Comput. , 2780 – J. Chem. Phys. , 080901(2017).65. Lin, C., Zong, F. H. & Ceperley, D. M. Twist-averaged boundary conditions incontinuum quantum Monte Carlo algorithms.
Phys. Rev. E., Statis., Nonlinear, SoftMatter Phys. , 016702 (2001).66. Kats, D. & Manby, F. R. Communication: the distinguishable cluster approxima-tion. J. Chem. Phys. , 021102 (2013).67. Kats, D. Communication: the distinguishable cluster approximation. II. The role oforbital relaxation.
J. Chem. Phys. , 061101 (2014).68. Raghavachari, K., Trucks, G. W., Pople, J. A. & Head-Gordon, M. A fi fth-orderperturbation comparison of electron correlation theories. Chem. Phys. Lett. ,479 –
483 (1989).
ACKNOWLEDGEMENTS
We thank Pablo López Ríos for providing the code to generate isotropic supercells,the original DMC static enthalpy-pressure data and some useful discussions. Thecomputational results presented have been achieved in part using the Vienna
K. Liao et al. Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences npj Computational Materials (2019) 110 cienti fi c Cluster (VSC) and the IBM iDataPlex HPC system HYDRA of the Max PlanckComputing and Data Facility (MPCDF). Support and funding from the EuropeanResearch Council (ERC) under the European Unions Horizon 2020 research andinnovation program (Grant Agreement No. 715594) is gratefully acknowledged. AUTHOR CONTRIBUTIONS
A.G. designed and led the research; K.L. performed the research; A.G. and K.L. wrotethe paper; X.Z.L. provided the crystal structures and A.A. advised and provided thetools to perform the FCIQMC calculations.
COMPETING INTERESTS
The authors declare no competing interests.
ADDITIONAL INFORMATION
Supplementary Information is available for this paper at https://doi.org/10.1038/s41524-019-0243-7.
Correspondence and requests for materials should be addressed to K.L. or A.G.
Reprints and permission information
Publisher ’ s note Springer Nature remains neutral with regard to jurisdictional claimsin published maps and institutional af fi liations. Open Access
This article is licensed under a Creative CommonsAttribution 4.0 International License, which permits use, sharing,adaptation, distribution and reproduction in any medium or format, as long as you giveappropriate credit to the original author(s) and the source, provide a link to the CreativeCommons license, and indicate if changes were made. The images or other third partymaterial in this article are included in the article ’ s Creative Commons license, unlessindicated otherwise in a credit line to the material. If material is not included in thearticle ’ s Creative Commons license and your intended use is not permitted by statutoryregulation or exceeds the permitted use, you will need to obtain permission directlyfrom the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.© The Author(s) 2019 K. Liao et al. npj Computational Materials (2019) 110npj Computational Materials (2019) 110