Quantum Monte Carlo study of high pressure solid molecular hydrogen

Diffusion quantum Monte Carlo [30] (DMC) is a stochastic projector method for solving the imaginary-time many-body Schrodinger equation

$$\begin{equation} -\partial_{t} \Phi(\mathbf{R},t)=(\skew3\widehat{H}-E_{T}) \Phi(\mathbf{R},t), \end{equation} \tag{ 1 }$$

where t is a real variable measuring the progress in imaginary time, R = (r₁,r₂,...,r_N) is a 3N-dimensional vector specifying the coordinates of all N electrons in the system, and E_T is an energy offset. Defining the Green function

$$\begin{equation} G(\mathbf{R} \longleftarrow \mathbf{R}^\prime, \tau)=\langle \mathbf{R} \vert \exp[-\tau(\skew3\widehat{H}-E_{T})]\vert \mathbf{R}^\prime \rangle, \end{equation} \tag{ 2 }$$

one can rewrite the imaginary-time many-body Schrödinger equation in the integral form

$$\begin{equation} \Phi(\mathbf{R},t+\tau)=\int G(\mathbf{R} \longleftarrow \mathbf{R}^\prime, \tau) \Phi(\mathbf{R}^\prime ,t) \mathrm{d}\mathbf{R}^\prime . \end{equation} \tag{ 3 }$$

Utilizing the spectral expansion,

$$\begin{equation} \exp(-\tau \skew3\widehat{H})=\sum_{i} \vert \Psi_{i} \rangle \exp(-\tau E_{i}) \langle \Psi_{i} \vert , \end{equation} \tag{ 4 }$$

the Green function can be expressed as

$$\begin{equation} G(\mathbf{R} \longleftarrow \mathbf{R}^\prime, \tau)=\sum_{i} \Psi_i (\mathbf{R}) \exp[-\tau (E_i - E_T)] \Psi_{i}^{*} (\mathbf{R}^\prime) , \end{equation} \tag{ 5 }$$

where Ψ_i and E_i are eigenfunctions and eigenvalues of $\skew 3\widehat {H}$, respectively, and the sum is over the complete set of eigenfunctions. It is straightforward to show that, as τ → ∞, the operator $\exp [-\tau (\skew 3\widehat {H} - E_T)]$ projects out the lowest eigenstate |Ψ₀〉 that has non-zero overlap with the starting state |Φ(R,t = 0)〉. This fundamental property of the projector $\exp [-\tau (\skew 3\widehat {H}- E_T)]$ is the basis of the DMC method and similar projector-based approaches.

In a practical DMC simulation, equation (1) is expressed in real space

$$\begin{equation} \partial_t \Phi(\mathbf{R},t) = \frac{\hbar^2}{2m}\mathbf{\nabla}^2_{\mathbf{R}} \Phi(\mathbf{R},t) - (V(\mathbf{R})-E_T) \Phi(\mathbf{R},t), \end{equation} \tag{ 6 }$$

where $\mathbf {\nabla }^2_{\mathbf {R}} = \sum _{i=1}^{N}\mathbf {\nabla }_{\mathbf {r}_i}^2$ is the 3N-dimensional Laplacian. Equation (6) may be interpreted as describing the time evolution of the probability density of a population of abstract 'walkers'—points R in the 3N-dimensional phase space—diffusing under the influence of the Laplacian term whilst simultaneously multiplying or dying out at a rate dependent on the local sign and magnitude of V (R) − E_T . This intuitive picture can be derived mathematically by using the Trotter decomposition [31] to obtain an approximation to G(R←R',τ) valid in the small τ limit. The short-time propagator consists of a Gaussian diffusion term multiplied by a reweighing term that describes the multiplication/death process. A DMC simulation uses this prescription to evolve a population of a few hundred walkers explicitly, time step by time step. After many time steps the walker population settles down to sample the ground-state wave function and ground-state expectation values can be accumulated.

Due to the antisymmetry property of fermions, the many-body electronic wave function cannot be positive everywhere and has to take both positive and negative values. The weak point of the naive DMC algorithm described above is that it can handle only non-negative wave functions. Fixed-node DMC is the most widely used method for dealing with the fermionic antisymmetry. The wave function Φ(R,t) is multiplied by a fermionic trial wave function Ψ_T (R) and equation (6) is replaced by an equivalent equation for the imaginary-time evolution of the product f(R,t) = Ψ_T (R)Φ(R,t). The initial wave function Φ(R,t = 0) is set equal to Ψ_T (R), ensuring that f(R,t = 0) is non-negative. The evolution equation for f is then solved subject to the constraint that f(R,t) is zero wherever Ψ_T (R) is zero, which ensures that f never becomes negative. It may be shown [30] that this procedure yields the energy expectation value of the variationally optimal wave function with the same nodal surface (surface of zeros) as the chosen many-electron trial function. The introduction of the trial wave function (the form of which is discussed below) also vastly reduces the statistical noise in the results. Indeed, if Ψ_T is equal to the exact ground-state Ψ₀, the fixed-node method yields the exact ground-state energy with zero statistical noise. If Ψ_T differs from Ψ₀ but has the same nodal surface, the exact ground-state energy is obtained with a non-zero statistical error bar. If the nodal surface of Ψ_T differs from that of Ψ₀, the expected value of the fixed-node ground-state energy is a variational upper bound on the true ground-state energy.

The fixed-node DMC method has a long history and has been thoroughly benchmarked [30, 32, 33]. For small systems, where even more accurate results can be obtained using quantum chemical methods, total energies calculated using fixed-node DMC are consistently more accurate than energies calculated using DFT and comparable in quality to energies obtained using the coupled cluster singles and doubles with perturbative triples method often referred to as the 'gold standard' of quantum chemistry [34–37]. Moreover, unlike DFT, fixed-node DMC provides a good description of weak interactions [38–44] such as those acting between the H₂ molecules in solid molecular hydrogen. DMC also describes the dependence of the total energy on coordination number (very important when comparing the energies of different crystal structures) and the form of the potential surface during bond-breaking processes more accurately than DFT. For large systems such as the ones studied in this paper, fixed-node DMC is the most accurate known ground-state total energy method and the only possible benchmarks are against experiment. Again, however, DMC consistently outperforms DFT. For example, the measured cohesive energy of crystalline Si is 4.62 eV atom⁻¹; fixed-node DMC simulations using a trial function similar in form to the trial functions used in this paper yield 4.63 eV atom⁻¹ [45]; and DFT in the local spin-density approximation yields 5.28 eV atom⁻¹. In this particular (unusually favourable) case, the DMC result is 65 times more accurate than the DFT result. The ability of DMC to describe changes in coordination is illustrated by the example of the Si self-interstitial defect [45]. The fixed-node DMC value of the defect formation energy, which is in good agreement with experiment, is about 1 eV larger than the PW91 GGA value and 1.5 eV larger than the LDA value. The main drawback of DMC, which explains why it is not more widely used, is the computational cost: a DMC simulation typically costs about 10⁴ times more than a DFT simulation of the same system.

Our fixed-node DMC simulations used the CASINO QMC code [32] and a trial function of Slater–Jastrow (SJ) form,

$$\begin{equation} \Psi_{\rm SJ}({\bf R})=\exp[J({\bf R})]\det[\psi_{n}({\bf r}_i^{\uparrow})]\det[\psi_{n}({\bf r}_j^{\downarrow})], \end{equation} \tag{ 7 }$$

where R is a 3N-dimensional vector that defines the positions of all N electrons, r^↑_i is the position of the ith spin-up electron, r^↓_j is the position of the jth spin-down electron, $\exp [J({\bf{ R}})]$ is the Jastrow factor and $\det [\psi _{n}({\bf{ r}}_i^{\uparrow })]$ and $\det [\psi _{n}({\bf{ r}}_j^{\downarrow })]$ are Slater determinants of spin-up and spin-down one-electron orbitals. These orbitals were obtained from DFT calculations using the plane-wave-based Quantum Espresso code [46]. A norm-conserving pseudopotential constructed within DFT using the PBE XC functional was employed [47]. Previous work has shown that the pseudopotential approximation has a very small impact on the results of calculations of high-pressure solid hydrogen [25]. We chose a very large basis-set cut-off of 300 Ryd to guarantee converge to the complete basis set limit [48]. The plane-wave orbitals were transformed into a blip polynomial basis [49]. Our Jastrow factor consists of polynomial one-body electron–nucleus (en) and two-body electron–electron (ee) terms, the parameters of which were optimized by variance minimization at the variational Monte Carlo (VMC) level [50, 51].

Geometry and cell optimizations were carried out using DFT with a dense 16 × 16 × 16 k-point mesh. The BFGS quasi-Newton algorithm, as implemented in the Quantum Espresso package, was used for both cell and geometry optimization, with convergence thresholds on the total energy and forces of 0.01 and 0.1 mRy Bohr⁻¹, respectively, to guarantee convergence of the total energy to less than 1 meV proton⁻¹ and the pressure to better than 0.1 GPa. An investigation of the proton–proton distance of accurately relaxed structures in the pressure range of 100–450 GPa reveals two main features. Firstly, as the pressure increases, the intra-molecular H₂ bond length (nearest-neighbour distance) increases slightly, rising from 0.735 to 0.778 Å. This could be a consequence of a weakening of the H₂ bonds as the coordination of the H atoms increases. Second, unlike the intra-molecular H₂ bond length, the inter-molecular distance decreases smoothly as the pressure rises. Details of the behaviour of the inter-atomic distance with pressure and of the experimental evidence are presented in [26].

Band gaps obtained from DFT simulations are not formally equivalent to the quasi-particle band gaps measured experimentally and depend strongly on the choice of XC functional. Local and semi-local functionals, in particular, yield DFT band gaps much smaller than experimental gaps. DMC does not suffer from the band-gap underestimation problem and can be used to obtain limited information about excited states. However, other approaches such as the GW method used in this work are generally better for excitation spectra. Many-body perturbation theory within the GW approximation has been widely tested and shown to yield highly accurate bandgaps of semiconductors and weakly correlated insulators. The results of GW calculations of the band gaps of sp-bonded materials generally match experiment to within a few tenths of an eV [52].

The quasi-particle spectrum of a solid can be determined from the Dyson equation [53],

$$\begin{equation} (T+V_{\mathrm{n}}+V_{\mathrm{H}})\psi_{\mathrm{n}\mathbf{k}}(\mathbf{r})+\int \mathrm{d}\mathbf{r}^\prime \, \Sigma (\mathbf{r}, \mathbf{r}^\prime ;E_{\mathrm{n}\mathbf{k}}) \psi_{\mathrm{n}\mathbf{k}}(\mathbf{r}^\prime) = E_{\mathrm{n}\mathbf{k}} \psi_{\mathrm{n}\mathbf{k}}(\mathbf{r}), \end{equation} \tag{ 8 }$$

where E_nk are the quasi-particle energies, T is the kinetic-energy operator, V_n is the potential due to the nuclei, V_H is the Hartree potential and Σ is the electron self-energy operator. In the GW approach, the self-energy operator is approximated by the first term of an expansion in the dynamically screened Coulomb interaction [53]

$$\begin{equation} \Sigma_{GW}=(\mathrm{i}/2\pi)\int \mathrm{d}\omega \, \exp(-\mathrm{i}\delta \omega) G(\mathbf{r}, \mathbf{r}^\prime, E-\omega) W(\mathbf{r}, \mathbf{r}^\prime ,\omega) , \end{equation} \tag{ 9 }$$

where W = ⁻¹V_C is the screened Coulomb interaction, V_C is the bare Coulomb interaction and is the dielectric matrix. The electron Green function G is given in the quasi-particle approximation by

$$\begin{equation} G(\mathbf{r}, \mathbf{r}^\prime ,E-\omega)=\sum_i [(E-\omega)-E_i-\mathrm{i}\delta_i]^{-1} \psi^{*}_i(\mathbf{r}) \psi_i(\mathbf{r}^\prime) , \end{equation} \tag{ 10 }$$

where ψ_i is the DFT wave function obtained by solving the Kohn–Sham equations, δ_i = 0⁻ for a quasi-electron and δ_i = 0⁺ for a quasi-hole. The static dielectric matrix is obtained within the random-phase approximation, _RPA = 1 − V_Cχ₀, where

$$\begin{equation}\fl \chi_{\mathbf{G}\mathbf{G}^\prime}^0(\mathbf{q})=\frac{4}{\Omega}\sum_{c,\upsilon , \mathbf{k}} { \frac{ \langle \upsilon ,\mathbf{k} \vert \exp(-\mathrm{i}(\mathbf{q}+\mathbf{G})\cdot\mathbf{r}) \vert c, \mathbf{k}+\mathbf{q} \rangle \langle c, \mathbf{k}+\mathbf{q} \vert \exp(\mathrm{i}(\mathbf{q}+\mathbf{G}^\prime)\cdot\mathbf{r}^\prime) \vert \upsilon , \mathbf{k} \rangle }{\varepsilon_{\upsilon ,\mathbf{k}}-\varepsilon_{c,\mathbf{k}+\mathbf{q}}}} \end{equation} \tag{ 11 }$$

is the independent-particle polarizability, |υ,k〉 and ε_υ,k are the valence-band DFT eigenvectors and eigenvalues and |c,k〉 and ε_c,k are the conduction-band DFT eigenvectors and eigenvalues. The static dielectric matrix (ω = 0) is extended to finite frequencies using a generalized plasmon-pole model [53].

In this work, all GW quasi-particle band gaps were calculated using the Yambo code [54]. Non-self-consistent GW calculations were performed starting from DFT wave functions, using a plasmon pole approximation for the dielectric function [55, 56]. Following convergence checks, the numbers of unoccupied bands and k-points were set to 200 and 512, respectively. We also checked convergence with respect to the number of plane waves, which was set to 1800 for the P6₃/m structure and 2700 for the Cmca-12 and C2/c structures. It is important to point out that our GW and DMC results are independent of one another. DMC is used to calculate the enthalpy as a function of pressure, whilst GW is used to calculate the quasi-particle band gap as a function of pressure. Both methods use the same relaxed DFT crystal structures, but they are otherwise unrelated.

Our previous DFT work [20] showed that hybrid XC functionals including a small fraction of exact exchange favour insulating phases relative to metallic phases. Very recently [28], we have also shown that DFT results for high-pressure solid molecular hydrogen, at both the static (ignoring the effects of proton ZPE) and dynamic (including the effects of proton ZPE) levels, depend strongly on the choice of XC functional. Since DFT orbitals are used to construct our DMC trial wave functions and thus to define the fixed nodal surface, it is important to investigate how the choice of density functional affects the DMC results. For this purpose, we compared the VMC and DMC energies of the insulating P6₃/m phase (which has the largest band gap of the four structures considered) using orbitals calculated using the semi-local PBE GGA functional [29] and the PBE0 hybrid functional including 25% exact exchange [57]. The VMC and DMC calculations were performed at the Γ point of the Brillouin zone of a simulation cell containing 128 hydrogen atoms and 2 × 2 × 2 replicas of the P6₃/m unit cell. The VMC and DMC energies at pressures 100, 200 and 300 GPa, are given in table 1.

The difference between the PBE and PBE0 energies is about 6 meV atom⁻¹ at the VMC level and 1 meV atom⁻¹ at the DMC level. Similar calculations for the Cmca-12 structure at two different pressures gave similar results. This demonstrates that, in the case of solid molecular hydrogen, DMC energies calculated using PBE and PBE0 orbitals differ negligibly, even though the change of functional has a large effect on the DFT phase diagram [20, 28]. The very weak dependence of the DMC energy on the choice of single-particle orbitals suggests, furthermore, that the fixed-node errors are small. This contrasts with the case of 3d transition metal compounds [58], where changing the DFT XC functional does affect the DMC results.

Computational cost prevented us from using backflow trial wave functions [59–62] throughout this work, but we did test the accuracy of the trial nodal surface by carrying out a few backflow DMC simulations for the metallic Cmca phase in a 64 electron simulation cell. Based on an old idea due to Feynman [63], the backflow transformation allows the one-electron orbitals in the determinant to depend on the positions of the other electrons in the simulation. In practice this dependence is parameterized and the parameters optimized variationally, thus improving the trial nodal surface. The introduction of backflow decreases the DMC total energy by approximately 10 meV per proton. Much of this reduction is expected to cancel when comparing the energies of different phases, so the effect of backflow on energy differences is unlikely to be larger than a few meV per proton.

In DFT, the limit of infinite system size can be approached by improving the k-point sampling. Since the Coulomb energy depends only on the charge density, which is periodic in an infinite system, it can be calculated exactly given the charge density of one unit cell. In a QMC simulation, however, we sample random snapshots of the positions of the N point-like electrons within the simulation cell. The Coulomb energy of a snapshot is defined by replicating it periodically and solving Poisson's equation for the infinite periodic lattice. In a real crystal, the electrons in neighbouring cells would be in different positions, so this procedure introduces spurious correlations. Hence, although it is still possible to use the many-body equivalent of k-point sampling (often referred to as twist averaging [64]), increasing the number of k-points no longer produces exact results for an infinite system. The resulting FS effects are a significant source of error [30]. The following QMC results used C2/c, Cmca-12, P6₃/m and Cmca simulation cells containing 96, 96, 128 and 144 atoms, respectively, equivalent to 2 × 2 × 2, 2 × 2 × 2, 2 × 2 × 2 and 3 × 2 × 3 repetitions of the corresponding unit cells.

QMC FS errors are conventionally separated into one-body and many-body contributions. One-body (independent-particle) errors arise from the non-interacting kinetic, potential and Hartree energies and include shell-filling effects. These errors can be removed by twist averaging [64]. Many-body errors arise from the effects of exchange and correlation on the Coulomb and kinetic energies and are not removed by twist averaging. Various techniques may be used to reduce or cancel these errors, but none is entirely successful and care is required. The oldest approach is extrapolation, which remains useful. The use of the 'modified periodic Coulomb interaction' [65–67] reduces the Coulomb errors but not the kinetic energy errors and must therefore be combined with other techniques. The Kwee–Zhang–Krakauer (KZK) approach [68] applies corrections obtained from DFT calculations carried out using a modified XC functional explicitly designed to mimic the DMC many-body errors. Most recently, structure-factor-based approaches that reduce both the kinetic and potential energy many-body errors have been developed [69, 70].

To remove single-particle errors and eliminate shell effects in the kinetic energies of metallic systems, we use twist-averaged boundary conditions (TABC) [64]. A twist k_s is imposed by insisting that the many-electron wave function Ψ_ks obeys the boundary condition

$$\begin{equation} \Psi_{\mathbf{k}_s}(\mathbf{r}_1, \ldots, \mathbf{r}_i + \mathbf{L},\ldots, \mathbf{r}_N) = \exp(\mathrm{i}\mathbf{k}_s\cdot\mathbf{L}) \Psi_{\mathbf{k}_s}(\mathbf{r}_1, \ldots, \mathbf{r}_i,\ldots, \mathbf{r}_N) \end{equation} \tag{ 12 }$$

for all electrons i. Expectation values of observables are obtained by averaging over twist vectors k_s uniformly distributed over the simulation-cell Brillouin zone:

$$\begin{equation} \langle \widehat{O}\rangle = \frac{1}{N_{\rm twists}} \sum_{\mathbf{k}_s} \langle \Psi_{\mathbf{k}_s} \vert \widehat{O}\vert \Psi_{\mathbf{k}_s}\rangle. \end{equation} \tag{ 13 }$$

The twists can be chosen from a uniform Monkhorst–Pack grid, preferably offset from Γ, or can be chosen randomly, as in this work. The number of twists should be as large as computational resources allow.

We tested three different methods for dealing with the many-body FS errors. The first was the KZK approach [68], which uses a system-size-dependent LDA fitted to the results of DMC simulations of finite cubic simulation cells of uniform electron gas. DFT energies calculated using the KZK functional incorporate DMC FS errors within an approximation analogous to the LDA. To estimate the FS error in the DMC total energy of a given simulation cell, the DFT total energy of exactly the same simulation cell is calculated using the KZK functional. The difference between this value and the DFT energy of an infinite simulation cell calculated using the standard LDA provides an estimate of the DMC FS error. The KZK functional was not originally combined with twist averaging, but the combination is easy to implement (see below) and very successful.

The second and third methods tested made use of the structure-factor-based corrections proposed by Chiesa et al [69] and Drummond et al [70]. The second approach used the standard Ewald form of the periodic Coulomb interaction and Chiesa corrections for the kinetic and potential energies. The third approach combined the use of the modified periodic Coulomb interaction [65–67] with Chiesa corrections for the kinetic energy only.

To check the accuracy and convergence of the DMC energy as a function of the number of protons N in the simulation cell, we performed calculations for different N (and thus also different L). Table 2 compares the results obtained using all three finite-size correction methods for the metallic Cmca phase. The convergence of the KZK-corrected DMC energy with increasing N is excellent and there is no difficulty in reaching the accuracy of a few meV per proton necessary for calculations of the hydrogen phase diagram at high pressure and low temperature. Our DMC-TABC energies were averages over 12 random twists in the simulation-cell Brillouin zone, corresponding to 12 randomly translated Monkhorst–Pack k-point meshes in the primitive Brillouin zone. In the case of the metallic Cmca structure, we also carried out DMC-TABC calculations using 24 twists; the change in the total energy was about 1.4 meV proton⁻¹ at P = 400 GPa, implying that DMC energies obtained by combining TABC and KZK FS-corrections converge very rapidly as the number of twists is increased. We emphasize the importance of our approach to the treatment of FS errors as an application of TABC to real systems [71].

The KZK scheme removes the FS errors in DMC calculations of cubic homogeneous systems exactly by construction. Our systems are neither cubic nor homogeneous, but the KZK approach still works remarkably well. Its use of DFT could be considered a drawback, but it is important to bear in mind that any valid FS-correction method must yield the same total energy in the limit of a large enough simulation cell. The important question is not whether the unattainable limiting value is correct but how rapidly it is reached. For solid H we find that the system-size dependence of the KZK-corrected DMC total energy is noticeably weaker than the system-size dependence of energies corrected using other methods. The KZK-corrected energy also settles down to a system-size-independent constant more quickly than energies corrected using other methods.

The correct combination of TABC and KZK corrections has not previously been described in the literature so we provide a brief explanation here. We write the QMC ground-state energy of the infinite simulation cell as

$$\begin{equation} E^{\rm QMC}(\infty) = E^{\rm QMC}_{\rm TAV}(L) + \Delta E^{\rm FS}_{\rm TAV}(L), \end{equation} \tag{ 14 }$$

where

$$\begin{equation} E^{\rm QMC}_{\rm TAV}(L) = \frac {1}{N_{\rm twist}} \sum_{\mathbf{k}_s} E^{\rm QMC}(L,\mathbf{k}_s) \end{equation} \tag{ 15 }$$

is the twist-averaged DMC energy for a simulation cell of size L and ΔE^FS_TAV(L) is the required FS correction. This correction is approximated using the formula

$$\begin{equation} \Delta E^{\rm FS}_{\rm TAV}(L) \approx E^{\rm DFT,LDA}(\infty) - E^{\rm DFT,FS}_{\rm TAV}(L), \end{equation} \tag{ 16 }$$

where E^DFT,LDA(∞) is the DFT energy computed within the LDA using a fully converged k-point mesh (which in this work means 24 × 24 × 24) and

$$\begin{equation} E^{\rm DFT,FS}_{\rm TAV}(L) = \frac {1}{N_{\rm twist}} \sum_{\mathbf{k}_s} E^{\rm DFT,FS}(L,\mathbf{k}_s). \end{equation} \tag{ 17 }$$

The FS DFT energy, E^DFT,FS(L), is obtained from DFT-KZK calculations employing the same simulation cell and twists as our DMC-TABC results. The FS correction, ΔE^FS_TAV(L), therefore accounts both for the many-body FS errors and for any one-body FS errors not removed by the limited twist averaging employed in the DMC-TABC simulations.

The quasi-particle band gap is given by

$$\begin{equation} E_{\rm gap} = \epsilon_{N+1,N+1} - \epsilon_{N,N} = \Delta_{\mathrm{xc}}+\epsilon_{N,N+1}-\epsilon_{N,N} , \end{equation} \tag{ 18 }$$

where _N,M is the energy of the Mth Kohn–Sham orbital in an N-electron system. The band-gap discontinuity, Δ_xc = _N+1,N+1 − _N,N+1, is by definition the energy difference between the highest occupied Kohn–Sham orbital in an (N + 1)-electron system and the lowest unoccupied Kohn–Sham orbital in an N-electron system. The last two terms in equation (18) give the band gap of the N-electron Kohn–Sham system. The non-zero value of Δ_xc results from the existence of a discontinuous jump [72] in the XC potential V_xc on addition of an electron: Δ_xc = V_xc,N+1( r) − V_xc,N( r). The lack of explicit position dependence on the left-hand side of this equation indicates that Δ_xc is a position-independent constant.

Although some of the band-gap underestimation observed in Kohn–Sham DFT calculations is due to the inadequacies of local and semi-local XC functionals such as LDA and PBE, a significant part is due to Δ_xc. Hybrid functionals such as PBE0 [57], which is formed from PBE by replacing 25% of the semi-local exchange energy by a fully non-local Hartree–Fock-like expression, yield better band gaps for most systems, including solid molecular hydrogen [16], but not necessarily for the right reasons. The GW method described in section 2.2 calculates true quasi-particle energies and yields consistently accurate band gaps for weakly-correlated sp-bonded solids.

In the following, we present the pressure dependence of the band gap of the Cmca-12 phase as calculated using the PBE, PBE0 and GW approximations. To obtain a very simple estimate of the metallization pressure of phase III, we linearly extrapolate the band gap as a function of pressure and look for the pressure at which the gap closes. Figure 1 illustrates band-gap closure in the Cmca-12 phase. The PBE, PBE0 and GW results show that the band gap of Cmca-12 vanishes at 245, 309 and 373 GPa, respectively. The PBE metallization pressure is in agreement with the value estimated by Pickard and Needs [24]. From band-structure calculations we see that the Cmca-12 phase has an indirect band gap; the linear decrease of the band gap with pressure is not unexpected in indirect gap structures and is also reported experimentally [9].

Zoom In Zoom Out Reset image size

As usual, band gaps calculated using the semi-local PBE functional are much too small. A significant contributor to this band-gap underestimation is the inability of any analytic functional of the density to represent the discontinuous change in XC potential that occurs as a single electron is added to a system with a filled band [72]. By including a fraction of the exact non-local orbital dependent exchange energy, the PBE0 hybrid functional circumvents this problem to some extent. However, even the PBE0 band gap is considerably smaller than the GW band gap, which is expected to be more accurate.

Our GW results show that the band gaps of the P6₃/m, C2/c and Cmca-12 structures vanish at the static pressures of 446, 398 and 373 GPa, respectively. Semilocal DFT, hybrid-DFT [20] and GW calculations all show that P6₃/m has the largest band gap and, consequently, the largest metallization pressure of the three insulating structures considered. The GW indirect band gap and its pressure dependence agree very well with recently reported experimental results [73].

The agreement with experiment is gratifying but could be misleading, since the large ZP fluctuations of the proton positions in solid H may reduce the band gap significantly. Indeed, Morales et al [74] recently investigated the influence of ZP motion on the DFT band gap of solid hydrogen and observed a dramatic reduction. Their paper did not provide details of the method used, but it was probably as described in [75]. Unfortunately, as Morales et al point out, their results disagree with experiment, for reasons that remain unclear. A similar approach could in principle be applied to the GW band gap of solid H, but the computational cost would be very high.

In the following, we present enthalpy-pressure diagrams obtained from DMC simulations at the static level of theory, which ignores proton ZPE. Calculating the pressure accurately in DMC is difficult. However, we note that the enthalpy H(P) is related to E(V ) via a Legendre transformation,

$$\begin{equation} H(P) = \mbox{Min}_{V} ( E(V) + PV ) , \end{equation} \tag{ 19 }$$

where the equilibrium volume V_eq(P) minimizes the right-hand side at fixed P. Hence, if we choose a pressure P and evaluate the quantity H(P,V ) = E(V ) + PV at an approximate volume V =V_eq(P) + ΔV , we obtain

$$\begin{equation} H(P,V)= H(P) + \frac{1}{2} \left . \frac{\mathrm{d}^2E}{\mathrm{d}V^2} \right |_{V_{\rm eq}(P)} (\Delta V)^2 + \cdots. \end{equation} \tag{ 20 }$$

The leading contribution to the difference between H(P,V ), which is easy to evaluate, and H(P), which is not, is B(ΔV )²/2nV² per proton, where B = −V dP/dV =V d²E/dV² is the bulk modulus of high-pressure hydrogen [76, 77] and n is the number of protons per unit volume. In practice, we set V to the DFT volume, V_DFT(P). Since DFT volumes are normally quite accurate, we estimate that ΔV/V ≃3%. By plotting P as a function of V , one finds that −dP/dV ≃4 GPa a⁻³₀ at 300 GPa. The volume of a unit cell containing 16 protons at this pressure is approximately 159 a³₀. The error in the enthalpy is thus approximately 3 meV per proton.

Figure  2 illustrates H_stat as a function of pressure. The static lattice enthalpies of three different insulating crystal structures are shown relative to the enthalpy of the metallic Cmca structure. Within the static approximation, the P6₃/m structure is the most stable up to a pressure of approximately 250 GPa, but is less stable than both C2/c and Cmca-12 by 350 GPa. The error bars make it difficult to establish exactly where the transition from P6₃/m to C2/c or Cmca-12 takes place or to say which of C2/c or Cmca-12 is favoured. The transition from C2/c or Cmca-12 to the metallic Cmca phase happens around 375 GPa.

Zoom In Zoom Out Reset image size

Thus, within the static approximation, DMC predicts that the hcp P6₃/m structure is more stable than the other structures investigated over the entire pressure range at which phases I and II exist experimentally, and that it remains more stable well above the measured transition to phase III. As already discussed above, recent x-ray results [14] at pressures up to 183 GPa suggest that hcp symmetry is a strong candidate not only for phase II but also for phase III.

We claim that our static DMC phase diagram is much more accurate than the static DFT phase diagram. One reason is the importance of vdW interactions between hydrogen molecules, which significantly influence the phase diagram, especially at low densities. QMC calculations have recently shown that the vdW interaction has a substantial effect on the transition pressures between the crystalline phases of ice at high pressure [40] and one would expect them to matter even more in solid molecular hydrogen. An important consequence is that transition pressures obtained from standard DFT functionals, which neglect vdW forces, are likely to be unreliable. Because of this recognition, some very recent work [74, 78] has used vdW density functionals to study solid and liquid H. There is no need to emphasize that DMC provides a highly accurate description of exchange and correlation in relatively weakly correlated materials and that it includes an excellent quantitative treatment of vdW forces.

Combining our static DMC phase diagram and GW results suggests: (a) that the well known rule of thumb 'the lower the energy the wider the gap' [79] is valid here (with the word 'energy' replaced by 'enthalpy'); and (b) that metallization occurs via a transition from the zero band-gap Cmca-12 structure to the metallic Cmca structure at around 375 GPa. The agreement between the pressure at which the GW band gap of the Cmca-12 structure closes and the pressure at which DMC predicts a phase transition from Cmca-12 to the metallic Cmca structure is remarkable but fortuitous, since including the effects of phonon ZPE changes the picture.

The above scenario is based on the static approximation. In reality, even in the ground state, H₂ molecules are quantum oscillators with a natural frequency of  ∼ 10¹⁴ Hz (for the H–H bond stretching) and a ZPE of several hundred meV per atom. Following standard practice in the field, our treatment of ZPE uses the harmonic approximation.

Estimating the ZPE within the harmonic approximation is feasible and was until recently believed to produce accurate results for molecular solid hydrogen [24]. A few months ago, however, Morales et al [74] published a paper reporting the results of DFT-based path-integral molecular dynamics simulations of solid hydrogen beyond the harmonic approximation. The quasi-harmonic and fully anharmonic total energies agreed well at high temperature but began to differ below about 300 K. At 200 K, the error in the quasi-harmonic energy had grown to approximately 8 meV per proton for the Cmca-12 phase at  ∼ 265 GPa. Path-integral molecular dynamics simulations become more difficult as the temperature decreases, so temperatures lower than 200 K were inaccessible.

The quasi-harmonic energy was seen to be only weakly dependent on temperature below 300 K because the vibrational modes were frozen out, but the anharmonic energy continued decreasing almost linearly down to 200 K. The specific heat of k_B/2 per atom was consistent with the assumption that every H₂ molecule was acting as a free quantum rotor at a temperature well above the rotational energy level spacing. To estimate the anharmonic ZPE it was necessary to extrapolate the path-integral molecular dynamics results from T = 200 to 0 K. This extrapolation was carried out on the assumption that the hydrogen molecules continued to act like free quantum rotors. The difference between the quasi-harmonic ZPE and the extrapolated anharmonic ZPE was more than 20 meV per proton.

This surprisingly large difference is probably an overestimate. It is clear from the data that the molecules behave like free rotors at 200 K. However, since the Cmca-12 structure is orientationally ordered, the molecules must in fact be rotating in a shallow potential well. This implies the existence of a temperature between 200 and 0 K below which the effects of the potential well become important and the rotational modes begin to freeze out. This freezing temperature is lower than predicted by the quasi-harmonic approximation but probably higher than the freezing temperature of a free rotor as determined by its moment of inertia alone. All we can say for sure is that the difference between the harmonic ZPE and the anharmonic ZPE is between 8 and 20 meV per proton. A more relevant question is how the anharmonicity affects the relative stabilities of different crystal structures. According to Morales et al [74], the anharmonic contributions to energy and enthalpy differences are of order 5–10 meV per proton. This may also be an overestimate.

Another issue is the limited accuracy of DFT. We showed recently [28] that the quasi-harmonic ZPE depends on the choice of XC functional, increasing by up to 20 meV per proton as the functional is changed from LDA to PBE to Becke–Lee–Yang–Parr. The dependence of energy differences on functional is considerably smaller, but it remains unclear whether DFT is accurate enough to provide reliable results. The neglect of the fermionic nature of the protons may also be problematic.

It is clear that further work is required to assess the effects of the quasi-harmonic approximation, the choice of XC functional, and the neglect of the quantum statistics of the protons on calculations of the ZPE of solid molecular hydrogen. It may prove necessary to go beyond DFT and carry out DMC simulations with both electrons and ions treated quantum mechanically. This is possible in principle, although the computational cost would be enormous. For this work we use quasi-harmonic ZPEs calculated within DFT. Reassuringly, although errors at the 10 meV level change calculated transition pressures significantly, they are unlikely to change the qualitative structure of the phase diagram. Ignoring the ZPE altogether introduces much larger quantitative errors, but a comparison of figures 2 and 3 shows that this too leaves the qualitative structure of the DFT phase diagram more or less intact.

Zoom In Zoom Out Reset image size

The quasi-harmonic ZPE was obtained from phonon frequencies calculated using density functional perturbation theory as implemented in Quantum Espresso [46]. The ZPE per proton was estimated as $E_{\mathrm { ZPE}}(V) =3\hbar \overline {\omega }/2$, where $\overline {\omega } = \sum _{\bf{ q}}\sum _{i=1}^{N_{\mathrm { mode}}} \omega _{i}({\bf{ q}})/(N_{\bf{ q}} N_{\mathrm { mode}})$. Here N_mode and N _q are the number of vibrational modes in the unit cell and the number of phonon wave vectors q, respectively, ω_i( q) is the angular frequency of the phonon in band i with wave vector q and the summation over q includes all k-points on a 2 × 2 × 2 grid in the Brillouin zone. McMahon [25] demonstrated that a q-point grid of this density is sufficient to converge ZPE differences between structures to within a few percent.

In agreement with previous DFT results [24] we find that the ZPE is large, structure dependent, and density dependent, ranging from 257 meV per proton at low density to 351 meV per proton at high density. The P6₃/m structure has the largest and the Cmca structure the smallest ZPE over the whole range of pressures considered. The C2/c and Cmca-12 structures have almost the same ZPE, which lies between the P6₃/m and Cmca values.

By defining H_dyn = E_dyn + PV , where E_dyn = E_stat + ZPE, we obtain the dynamic phase diagram of solid molecular hydrogen illustrated in figure 3. The dynamic lattice enthalpies of three different insulating crystal structures are shown relative to the enthalpy of the metallic Cmca structure. Bearing in mind the statistical error bars, we see that P6₃/m is the most stable phase up to a pressure of approximately 220 GPa. This is in good agreement with experimental results [14]. As discussed above, because of the tiny enthalpy difference between the C2/c and Cmca-12 phases, it is difficult to determine their pressure stability ranges accurately. However, we are able to conclude that we expect to observe the following structural transitions with increasing pressure: P6₃/m →C2/c → Cmca-12 →Cmca.

Comparing the static and dynamic results, we see that the proton ZPE has two main effects. Firstly, although the order in which the phases appear with increasing pressure is unaffected, their enthalpies and regions of stability change. If we were to take an optimistic view and ignore the error bars, we would say that P6₃/m, C2/c, Cmca-12 and Cmca phases are now stable in the pressure ranges <220, 220–360, 360–430 and >430 GPa, respectively. As in the static case, however, the statistical and systematic errors make it difficult to discriminate between the C2/c and Cmca-12 structures over a wide range of pressures, this time from 300 to approximately 450 GPa. Secondly, the inclusion of ZPE increases the pressure of the transition to the metallic Cmca phase by 14%, which is close to the 10% increase previously reported [16].

Note that the static GW band gap still closes before the transition to the Cmca structure takes place. This leaves a region of 20–30 GPa where the Cmca-12 phase persists but is no longer insulating in the static approximation.

A comparison of the densities of different structures at the same pressure shows that, in phases II and III: density(Cmca) > density(Cmca-12) > density(C2/c) > density(P6₃/m). As the pressure is increased, the PV term in the enthalpy rises faster for the low-density insulating phases than for the metallic phase, reducing the relative stability of the insulating phases and leading eventually to the metal–insulator transition. The difference between the densities of the insulating phases and the metallic Cmca phase reduces at high pressure, especially in the case of the Cmca-12 structure.

The differences between the DMC and DFT phase diagrams are primarily due to the internal energy term in the enthalpy. The ZPE contributions to the pressure are quite small and have no significant effect on the PV term. For pressures up to 220 GPa, the DMC energy of the low density P6₃/m phase lies further below the energies of the other phases, especially the high-density metallic Cmca phase, than does the DFT energy. Our results, at both the static and dynamic levels, predict that the P6₃/m phase is the most stable of the structures considered for a large range of pressures up to 220 GPa in the dynamic case or 280 GPa in the static case. This prediction is consistent with the latest low temperature x-ray powder-diffraction results [14], which indicate that the hydrogen molecules in phases II and III are still in a hexagonal structure and that phase III remains hexagonal up to the maximum experimental pressure of 183 GPa. The nature of the difference between phases II and III remains to be established, but may be related to a classical orientational ordering of the H₂ molecules [80].