The pseudopotential approach within density-functional theory: the case of atomic metallic hydrogen

Internal energies, enthalpies, phonon dispersion curves, and superconductivity of atomic metallic hydrogen are calculated. The (standard) use pseudopotentials in density-functional theory are compared with full (Coulomb)-potential all-electron linear muffin-tin orbital calculations. Quantitatively similar results are found as far as internal energies are concerned. Larger differences are found for phase-transition pressures; significant enough to affect the phase diagram. Electron--phonon spectral functions $\alpha^2 F(\omega)$ also show significant differences. Against expectation, the estimated superconducting critical temperature T$_{c}$ of the first atomic metallic phase I4$_{1}$/amd (Cs-IV) at $500$ GPa is actually higher.


I. INTRODUCTION
Hydrogen is the simplest and the most abundant element in the universe. Under pressure, it exhibits remarkable physics. First it solidifies and crystallizes, and then evolves through a series of high-density solid phases. In 1935, Wigner and Huntington predicted 1 that sufficient pressure would even dissociate hydrogen molecules, and that any Bravais lattice of such atoms would be metallic. The problem of metallic hydrogen has received considerable attention, as reviewed in Ref. 2. Herein, the structures and stabilities of atomic metallic hydrogen are considered. The background of what is known (from calculations; as motivated below) and relevant to this work will be discussed in context. Initial interest in metallic hydrogen was primarily related to astrophysical problems 3 . Subsequently (and more recently), there has been significant interest in it at relatively low temperatures. This can be attributed to the remarkable properties that are expected. This includes, for example, high-temperature superconductivity 4-6 . This will be considered herein. The possibility of a zero-temperature liquid ground-state has also been suggested 7 . In this case, hydrogen may have quantum-ordered states that represent novel types of quantum fluids 8 . Applications of the (expected) remarkable physics could revolutionize several fields. Possible scientific investigations and technological uses have been speculated on in Refs. 9 and 10.
Despite experimental advances [e.g., diamond anvil cell (DAC) 11 experiments, even coupled with direct synchrotron X-ray diffraction 12 ], it is still extremely difficult to measure the crystal structure of hydrogen under extreme conditions. Therefore, sophisticated calculations, often ab initio ones based on density-functional theory (DFT) 13 have become a powerful theoretical tool to understand high-pressure hydrogen and its physical properties.
Pseudopotentials, the focus of this work, are an essential ingredient of most of these calculations. These potentials, which are smooth and nodeless, are used to replace the 1/r Coulomb potential, in order to reach more rapidly convergent results. This same idea applies to the case of hydrogen, even though it only has one electron and lacks core ones.
For many properties, it is reasonable to assume that the pseudopotential should be almost numerically identical to the Coulomb one, as long as the cutoff radius r c is chosen to be (extremely) small. Under high pressure, the distance between nearest-neighbor protons in atomic metallic hydrogen is approximately twofold of the Wigner-Seitz radius r s [V = (4π/3)r 3 s a 3 0 , where V is the volume per electron and a 0 the Bohr radius]. According to the evolution of shortest (interatomic) H-H distance under pressure, r s changes from 3.12 to 1.23 when the pressure increase from 1 atm to 500 GPa 14 . The concern comes to be that if the pseudopotential with cutoff radius is suitable to ensure minimal core overlap.
The validity of the pseudopotential approximation in the above contexts has been discussed by McMahon and Ceperley 15 . The internal energies of two structures, with Hermann-Mauguin space-group notation I4 1 /amd (c/a > 1) (the family of structures to which this belongs will be considered further below) and R3m, with different cutoff radii (0.5 and 0.125 a.u.) of normconserving Troullier-Martins pseudopotentials 16 , were compared. Their study indicated that this approximation has a very small impact on these calculations (subject to the above constraint). In another study 17 , the energy difference between face-(Fm3m; f cc) and bodycentered cubic (bcc) phases were compared, by using a projector-augmented wave (PAW) pseudopotential 18,19 and an all-electron method. This work showed that the error introduced for these calculations is insignificant. Note though that the structures considered in these studies have very high symmetry. Another important consideration is whether using a pseudopotential will influence the calculation of properties, such as the superconducting critical temperature T c . This was made long ago by Gupta and Sinha 20 , suggesting that the estimate of T c may be considerably reduced by screening effects. This is based on the idea 21 that, in the vicinity of the proton, the electron wavefunction is rigidly displaced together with the proton, and hence is not involved in the electronphonon interaction. That is, the screening of the bare Coulomb potential should result in a decrease of cou-arXiv:1912.01802v1 [cond-mat.mtrl-sci] 4 Dec 2019 pling constant λ. This will be discussed in more detail further below.
There are still several open comments and questions concerning the use of the pseudopotential method; some specific ones are as follows: Compared with the f cc and bcc phases, which both belong to cubic system of crystal structures, lower-symmetry ones may be more representative and convincing. How are the internal energies of these affected? Are transition pressures (being a function of both energy and its change to first order) affected? And, is the superconductivity-physics affected?
The purpose of this work is to make a thorough analysis of the error made using pseudopotentials, using modern calculation techniques. Calculations of internal energies, (the first) phase-transition pressure, and superconducting properties of atomic metallic hydrogen under high pressures are performed. Structures that come from different crystal systems (cubic, rhombohedral, tetragonal, and orthorhombic) are considered. These quantities will be compared as calculated within the pseudopotential method to the all-electron full-potential linearised augmented-plane wave (LAPW) 22,23 one.

II. COMPUTATIONAL METHODOLOGY
Both the pseudopotential and all-electron calculations were performed from first principles. These were based on DFT 13 . Exchange-correlation effects were described using the generalized gradient approximation (GGA), according to the Perdew-Burke-Ernzerhof (PBE) 24 form. (Other) settings were chosen similarly between the two methods, for as direct comparisons as (reasonably) possible later; these are described in the following.
The pseudopotential calculations were performed using Quantum ESPRESSO (QE) 25 . PAW pseudopotentials 18,19 with a cutoff radius of 0.75 a.u. was used to describe the region near the nucleus of hydrogen. Convergence tests (energy to within 1 meV/proton) for this pseudopotential required 57.5 and 345.5 Ry for the plane-wave basis-set cutoffs (kinetic energy) for the wavefunction and charge density, respectively.
All-electron calculations were performed selfconsistently using the full-potential LAPW as implemented in the Elk code 26 . A plane-wave cutoff of |G + K| max = 9/R M T min (R M T min is the average of the muffin-tin radii in the unit cell) was used for the expansion of the wavefunction in the interstitial region. The muffin-tin radii for H is 0.9 a.u. (comparable to that in PAW). The cutoff for charge density, which is the maximum length of |G| for expanding the interstitial density was considered to 2 |G + K| max + ε where ε = 10 −6 .
It is important to briefly recognize the difference between the all-electron (LAPW) and PAW methods. Both consider a plane-wave basis set, but augmented in the region near the nucleus to more accurately (while re-taining or increasing efficiency) describe the atomic-like wavefunction. For the PAW method, however, inside the augmentation region, the (pseudo) wavefunction will be much smoother than the all-electron one. That is, the physics in this region, for this method, are similar to what happens in the pseudopotential approximation.
Convergence (to the same criterion as above) with respect to the number of k points needed to sample (integrate over) the irreducible Brillouin zone were tested individually between QE and Elk. Values obtained for the considered structures were as follows: I4 1 /amd (26 3 both), Cmcm (26 3 and 20 3 for QE and Elk, respectively), I43d (26 3 and 28 3 ), and Fm3m (32 3 ). Smearing was used to improve convergence (of the integrations). In QE, the scheme of Methfessel-Paxton 27 was used, with a value of 0.02 Ry; in Elk, that of Fermi-Dirac 28 , with a suggested value 29 of 0.001 Ha.
For the phonon calculations, the GGA functional is implemented with the finite-displacement method (supercell method), but not with density-functional perturbation theory (DFPT) 30 in the current version of Elk (6.3.2). To make the comparison under the same conditions, phonon dispersions were calculated using the former approach with 4×4×4 supercell in both QE and Elk, combined with the phonopy code 31 . Note that such a grid is sufficient for a quantitative determination of the phonons in this system 5,6 . For phonon dispersions, paths between high-symmetry points (covering all special points and lines necessarily and sufficiently) in the Brillouin zone were determined automatically, using the SeeK-path tool 32 .
For the superconductivity calculations, again in order to use the GGA functional, electron-phonon coupling calculations were carried out using DFPT in QE, and the supercell method in Elk. The two methods should give (numerically) the same results, as long as the sampling in reciprocal space (former method) is consistent with the supercell size (latter method); that is, the difference is one of computational efficiency 33 . Considering this, a 4×4×4 q-point grid and supercell were used for all calculations. This should be sufficient to make a quantitative comparison between the two methods, even if only calculate approximate values of the superconducting parameters themselves 5,6 .
T c is estimated by numerically solving the two (complete) nonlinear Eliashberg equations. Detailed derivation of the isotropic Eliashberg gap equations have been presented by Allen and Mitrovic 34 . The following corresponding numerical method has been explained in Refs. 35 and 36. These are for the superconducting order parameter ∆ n ≡∆(i ω n ) along the imaginary frequency axis (i = √ −1), the maximum value of which corresponds to the wavefunction of the superconducting condensate, and wavefunction renormalization factor Z n ≡Z (i ω n ), and where β=1/k B T where k B denotes the Boltzmann constant and T the temperature, µ * is the Coulomb pseudopotential, θ is the Heaviside function, ω c is the phonon cut-off frequency, ω c = 3ω max where ω max is the maximum phonon frequency, ω n =(π/β)(2n-1) is the n th fermion Matsubara frequency with n = 0, ±1, ±2, . . ., the pairing kernel for electron-phonon interaction has the form λ(ω n − ω m ) = 2 where ω is the phonon frequency, and α 2 F (ω) is the Eliashberg spectral function where F (ω) is the density of states of lattice vibrations (the phonon spectrum), and α 2 describes the coupling of phonons to electrons on the Fermi surface. Ashcroft demonstrated 37 , via an ab initio calculation, that µ * = 0.089 in metallic hydrogen, which is similar to the (rather) standard value for a high-density system of µ * ≈ 0.1. The former value is used herein. These two equations are solved iterative self-consistently at a certain temperature T. T c is then defined as the temperature at which the Matsubara gap ∆ n become zero. Herein, 2201 Matsubara frequencies (M = 1100) have been used.
The most stable structures of atomic metallic hydrogen from 500 to 3000 GPa, as predicted by calculations, were considered. These include the I4 1 /amd (Cs-IV) 15 , Cmcm 38 , and I43d 38 . Lower-symmetry, related structures, essentially the same up to a distortion(s), (such as F ddd 17 and C222 1 39 for the first two structures, respectively) were not considered; Fm3m was also considered, for reference. The considered pressures cover the range from approximately the expected molecular-toatomic phase transition 40,41 to just above the first predicted atomic phase transition I4 1 /amd → Cmcm 38 .

III. RESULTS AND DISCUSSION
The structures (themselves) of high-pressure hydrogen are extremely difficult to determine by experiment. Based on first-principles calculations 15 , a body-centered tetragonal (BCT) is considered to be the most promising candidate, for the first atomic phase(s). Representations of structures from this family are shown in Fig. 1. This family of structures can be characterized in terms of their c/a ratio, and they are often done so using an "elemental" naming scheme. These are c/a < 1 (β-Sn type), ≈ √ 2 (diamond), and > 1 (Cs-IV).

A. Internal Energies
Internal energies as a function of c/a ratio were calculated at six (constant) volumes. Note that zero-point energies were not directly included in these (or below) calculations. This ratio was varied from 0.05 to 10, for (a) c/a < 1 (β-Sn type) each volume. Volumes were determined by geometry optimizations with QE over the considered pressure range (see above) in steps of 500 GPa. These (volumes) were then fixed, and used for both QE and Elk. The results are shown in Fig. 2. For both sets of calculations, there are four energy minima: a shallow one at c/a 1, the deepest one at c/a > 1 (Cs-IV type), and two deep ones at c/a < 1 (β-Sn) and c/a 1. Notice that c/a ≈ √ 2 (diamond) is always unstable. From the difference plot [ Fig. 2(c)], a few meV/proton difference (the PAW pseudopotential energies are, in general, higher) occur on both sides of c/a ≈ 3.5. While this difference does not change the relative stabilities of the (BCT) structures [see Fig. 2(a)], it is still significant, considering the magnitude of energies.
Consider also the changes as a function of volume. The global energy minimum is always for Cs-IV. As the volume decreases, c/a increases. For the pseudopotential calculations, this ranges from 2.53 to 3.03. For the allelectron ones, from 2.6 to 3.05. These ranges are in very good agreement. For both sets of calculations, the energies of β-Sn and diamond decrease with increasing volume.
The above results show that, as far as (relative) energies, structures, and both qualitative and quantitative changes with volume are concerned, the replacement of Coulomb potential by a pseudopotential appears to be reasonable. This is consistent with previous results 17 focused on structures with very high symmetries. In addition (in a way) to verifying the approach, the results here extend (together, generalize) these for structures with low(er) symmetries.

B. Phase Diagram
In order to quantify the aforementioned considerations with volume, the pressure-volume (pV ) phase diagram was constructed. This is a more sensitive measure [than internal energies (above)], as the free energy (enthalpy H, in this case) depends on both the energy and its firstorder changes via the (hydrostatic) pressure, where U is the internal energy. Note that pressures were calculated according to Eq. (3); by derivatives of the equation of state (EoS) with respect to volume (instead of directly calculating the trace of external stress tensor). Specifically, once the volume dependence is known, the energy as a function of volume can be constructed, then this data is fitted with the 3 rd -order Birch-Murnaghan EoS 42 , and derivatives are calculated. Results for the first (predicted) phases of atomic hydrogen are shown in Fig. 3. Now using Hermann-Mauguin space-group notation (as common), these are I4 1 /amd (Cs-IV), Cmcm, and I43d. Enthalpies relative to F m3m are shown. Note that values were calculated every 250 GPa. The pseudopotential results are in both qualitative and quantitative agreement with earlier work 38 . The allelectron ones show some important differences, however.
Consider first the trends in relative enthalpy differences. These are consistent with earlier work. In particular, I4 1 /amd becomes very unstable with increasing pressure, relative to a set of structures with much flatter enthalpy changes.
Consider now the phase transition pressures. That of the (first) I4 1 /amd → Cmcm transition is 2300 GPa, which is in agreement with the approximate value of > 2100 GPa calculated in Ref. 38 (the latter based on a lessdense pressure grid). For the all-electron calculations, the transition occurs at 2410 GPa. Compared to the above results (for the two pseudopotential calculations -herein and in Ref. 38), this difference (increase) is relatively small. But this trend appears consistent with the next (potential) phase transition, discussed below.
Consider now the latter structures. It appears that a phase transition Cmcm → I43d will occur. (Indeed, but with consideration of zero-point energy. This is predicted 38 above 3.5 TPa.) Considering this next phase transition, a significant difference can be seen. Consistent with the first transition, it appears that this one will also be pushed to even higher pressures. In this case, however, it is enough such that this transition may not occur. Consider the difference in enthalpy between these two structures, ∆H = H I43d − H Cmcm . The (maximum) value with the pseudopotential approximation is 4.5 meV/proton at 1700 GPa; and this decreases to 2.2 meV/proton by 3000 GPa. This is (even) qualitatively much different in the all-electron calculations, where ∆H increases from 7.6 to 7.8 meV/proton at these pressures. That is, a phase transition, in this case, seems unlikely.
Considering the results together, all-electron calculations seem to (at least, in this region of the phase diagram considered) push transition pressures higher. Relative stabilities may also change. These results may be significant enough to change the phase diagram.

C. Phonon Dispersion
An important consideration for phase stabilities (by zero-point energy), properties (e.g., superconductivity), etc. is lattice vibrations. Throughout reciprocal space, these are illustrated most clearly by phonon dispersions.
The case of I4 1 /amd (Cs-IV) at 500 GPa is considered, as an example. These results are shown in Fig. 4. Phonon dispersion curves along various symmetry directions were calculated. Comparison shows that the dispersion relations calculated by the methods are similar. The most significant difference is near the Γ point, where the frequencies of the optical modes as calculated with the pseudopotential approximation are much flatter. A possible explanation for this is that with the pseudopotential approximation, the electrons near the proton are not (as) bound with its motion (unlike the all-electron methodsee below); this means that the change in the electronic charge density become noticeable, and hence its phonon density of state is large and phonon dispersion flat.

D. Superconductivity
Superconductivity of atomic metallic hydrogen is considered, in this section.
That of I4 1 /amd at 500 GPa is again used, as an example. Figure 5 shows a detailed comparison of the electron-phonon spectral function α 2 F(ω) and coupling parameter λ. There are significant differences in the quantities as calculated by the two methods, both qualitatively and quantitatively. The pseudopotential calculations display significant (and "peaked") electron-phonon interaction at both (relatively) low and high frequencies, but much less at intermediate ones. This can be compared to the broad spectral function, centered at intermediate frequencies, as calculated by the all-electron method.
Only at high frequencies is this result is consistent with the work of Gupta and Sinha 20 . This interaction (Color online) Full dependence of the maximum value of the order parameter ∆m=1 on temperature for I41/amd at 500 GPa. Calculations are shown by both the all-electron method and with the PAW pseudopotential approximation.
is mainly due to that near the proton (in metallic hydrogen). (Consider the change in the bare Coulomb interaction with r; this scales as 1/r 2 , and hence is largest for small r.) However, the electrons in this vicinity are not at all free-electron-like; their motion is bound with that of the proton. These electrons will therefore not (significantly) participate in the electron-phonon interaction.
Unlike the earlier expectation 20 of a decrease in λ, it is actually found to increase by the all-electron calculation. This can be attributed to the increase contribution to intermediate frequencies to the α 2 F(ω).
The dependence of the maximum value of the order parameter ∆ m=1 on temperature is shown in Fig. 6. The superconducting transition temperature is defined as that at which this parameter vanishes, ∆ m=1 (T c , µ * ) = 0. The obtained T c is 352 K in the all-electron calculation, compared to 339 K with the pseudopotential approximation.

IV. CONCLUSIONS
In conclusion, the reliability of (the standard use of) pseudopotentials to simulate atomic metallic hydrogen was studied. This was done for calculations of internal energy, enthalpy, phonon dispersion spectrum and superconductivity, by comparing pseudopotential to allelectron calculations. In the case of calculating internal energy, as has been considered to some extent, the accuracy that can be obtained by PAW pseudopotentials is sufficient. Differences occur for enthalpy and phonon dispersion relations, however. These may significant enough to affect the phase diagram, by both pushing transition pressures higher and changing relative stabilities. Significant differences also occur for the calculation of (at least, some) properties. For superconductivity, for example, the magnitude of the electron-phonon spectral function at both (relatively) low and high frequencies is considerably smaller as calculated by the all-electron method than with the pseudopotential approximation, while that at intermediate frequencies is increased. Together, these changes actually increase the value of λ, which causes the calculated superconducting critical temperature to be higher. These results are important for understanding metallic hydrogen; and will be so for future calculations of this system.