Auxiliary-field quantum Monte Carlo calculations with multiple-projector pseudopotentials

Auxiliary-field quantum Monte Carlo calculations with multiple-projector pseudopotentials


We have implemented recently developed multiple-projector pseudopotentials into the planewave based auxiliary-field quantum Monte Carlo (pw-AFQMC) method. Multiple-projector pseudopotentials can yield smaller planewave cut-offs while maintaining or improving transferability. This reduces the computational cost of pw-AFQMC, increasing its reach to larger and more complicated systems. We discuss the use of non-local pseudopotentials in the separable Kleinman-Bylander form, and the implementation in pw-AFQMC of the multiple-projector optimized norm-conserving pseudopotential ONCVPSP of Hamann. The accuracy of the method is first demonstrated by equation-of-state calculations of the ionic insulator NaCl and more strongly correlated metal Cu. The method is then applied to calibrate the accuracy of density functional theory (DFT) predictions of the phase stability of recently discovered high temperature and pressure superconducting sulfur hydride systems. We find that DFT results are in good agreement with pw-AFQMC, due to near cancellation of electron-electron correlation effects between different structures.

71.15.-m, 02.70.Ss, 71.15.Dx, 61.50.-f

Present affiliation: ] Center for Advanced Quantum Studies and Department of Physics, Beijing Normal University, Beijing 100875, China

I Introduction

The search for new materials and their development has increasingly relied on theoretical modeling. Methods based on density functional theory (DFT) are efficient and powerful, but their predictions can break down in a number of instances. Examples range from strongly correlated materials, such as transition metal systems, to bond stretching or bond breaking in otherwise moderately correlated systems. Explicit many-body methods, which avoid the mean-field-like approximations used in standard DFT calculations, are needed in these cases. Quantum Monte Carlo (QMC) calculations have become increasingly important in this regard, because of their accuracy and favorable scaling (as a low-order polynomial of system size, similar to DFT, but with larger prefactor) compared to traditional wave function based correlated methods. Routine applications of QMC calculations in extended systems still face major challenges, however.

In diffusion QMC (DMC) Foulkes et al. (2001) and pw-AFQMC, Zhang and Krakauer (2003); Suewattana et al. (2007); Purwanto et al. (2013); Ma et al. (2015) pseudopotentials are usually used, except for some DMC calculations with the lightest elements. Pseudopotentials remove the chemically inactive core electrons, reducing the number of electrons that must be explicitly correlated and greatly reducing the computational cost. Non-local norm-conserving pseudopotentials (NCPP) are typically used in QMC. The NCPPs are usually constructed from mean-field DFT of Hartree-Fock (HF) calculations. While computationally expedient, the transferability of NCPPs is a key issue, and the neglected core-core and core-valence correlation effects may need to be considered. Even setting these many-body effects aside, transferability errors from NCPPs in QMC calculations can be significant. In DMC, moreover, the non-locality of NCPPs is handled with an additional locality approximation, whose accuracy depends on the quality of the trial wave function. Mitas et al. (1991) The overall NCPP error can be significant compared to errors from the fixed-node approximation, Hennig et al. (2010); Sorella et al. (2011); Nazarov et al. (2016) which is used to control the Fermion sign problem. In pw-AFQMC, non-local NCPPs can be used without additional approximations,Zhang and Krakauer (2003); Suewattana et al. (2007) but transferability errors can still be a problem, unless the NCPPs are made very hard, Suewattana et al. (2007); Purwanto et al. (2009, 2013); Ma et al. (2015) which requires large planewave cutoffs and increases the computational cost.

Pseudopotentials are based on the frozen-core approximation, but contain an additional layer of approximation. Frozen-core calculations are common in quantum chemistry applications, where the core orbitals are frozen at the mean-field level derived from the target system. Pseudopotentials are usually constructed for a reference atomic configuration and then used in many target systems. The accuracy (transferability) of the PP across many target systems must then be determined a posteriori. In addition to being norm-conserving, most NCPPs used in QMC calculations are of single-projector type (one per angular momentum channel), which can further contribute to transferability errors.

Recently, Hamann proposed a multiple-projector pseudopotential,Hamann (2013) based on Vanderbilt’s norm-conserving construction Vanderbilt (1990) and optimized with the Rappe-Rabe-Kaxiras-Joannopoulos pseudization scheme. Rappe et al. (1990) The resulting pseudopotential, referred to as ONCVPSP by Hamann, was shown to have accuracy comparable to all-electron (AE) and ultrasoft pseudopotentials Vanderbilt (1990) (USPP) in DFT calculations, with moderate planewave energy cutoffs. Schlipf and Gygi Schlipf and Gygi (2015) recently presented a set of automatically constructed Hamann ONCVPSPs for most of the periodic table. These were shown to be in good agreement with the all-electron results in DFT, often with cutoffs of only about 40 Ry. Schlipf and Gygi (2015); Lejaeghere et al. (2016) The ONCVPSP is of separable Kleinman-Bylander type, Kleinman and Bylander (1982) similar to NCPPs widely used in planewave DFT calculations and also in pw-AFQMC calculations. Since the treatment of one-particle Hamiltonian terms in pw-AFQMC is closely related to that in planewave DFT, the implementation of ONCVPSP into our pw-AFQMC code is straightforward.

In this paper, we show that the use of multiple-projector ONCVPSP can greatly reduce the planewave basis size in pw-AFQMC many-body calculations, while maintaining good accuracy. This results in significant reductions of computational cost, both by reducing the computing time for each step in the random walks and, at the same time, by reducing QMC statistical variance, due to the reduced number of AFQMC auxiliary fields. To test the new capability with multiple-projector ONCVPSPs, we carry out pw-AFQMC calculations of the equation-of-state in the insulator NaCl and the transition metal solid Cu. We then study the high-pressure superconducting system HS, to calibrate DFT predictions of phase stabilities. Finally we discuss the performance of the DFT- or HF-generated pseudopotentials in many-body calculations and the difference from their use in DFT calculations.

The reminder of the paper is organized as follows. Section II reviews AFQMC with a planewave basis and pseudopotentials, and discusses the implementation of multiple-projector separable pseudopotentials in pw-AFQMC. Section III presents applications of the method. Additional transferability issues and other aspects of ONCVPSP for many-body applications are discussed in Section IV. We then conclude with some general remarks in Section V.

Ii pw-AFQMC methodology

To set the context for the implementation of multiple-projector NCPPs in pw-AFQMC, we briefly review pertinent aspects of the formalism in this section. For more details about the pw-AFQMC method, see Refs. [Zhang and Krakauer, 2003; Suewattana et al., 2007; Purwanto et al., 2009].

ii.1 Hamiltonian

The electronic Hamiltonian within the Born-Oppernheimer approximation is,


where , , , and are, respectively, the kinetic energy and electron-electron, electron-ion, and classical Coulomb ion-ionYin and Cohen (1982) interactions. The pseudopotential contributions appear in the electron-ion interaction . With periodic boundary conditions and a planewave basis,


the terms in Eq. (1) can be expressed in second quantized form as


Here, () is a creation (destruction) operator, is the volume of the simulation cell, and is a reciprocal lattice vector, , is the electron spin, and is the number of electrons in the simulation cell. Both and belong to the planewave basis set whose size is controlled by the planewave kinetic energy cut-off . (When twist-averaged boundary conditions are used, is replaced by , where is within the first Brillouin zone.) The constant gives the self-interaction of an electron with its periodic images. Fraser et al. (1996) The one-body density operator is given by


where the step function ensures that , like , falls within the planewave basis set.

The local and non-local parts of the pseudopotential are defined by the planewave matrix elements and , respectively, which are discussed in more detail in Section II.3.

ii.2 Ground state projection

AFQMC uses iterative imaginary-time projection to obtain the ground state from a trial wave function (often just a single Slater determinant):


where is assumed. The projection is implemented as random walks in the space of Slater determinants. A key point in implementing this is the observation that a one-body propagator acting on a Slater determinant simply yields another Slater determinant. The AFQMC procedure is therefore to separate the propagator in Eq. (5) into one- and two-body propagators. This motivates the introduction of a small imaginary-time step :


A Trotter-Suzuki decomposition Trotter (1959); Suzuki (1976) then achieves the desired separation,


where and are the one- and two-body parts of the Hamiltonian in Eq. (1), with and denoting the remaining terms in Eq. (3) and the collection from Eqs. (3a) and (3c).

A Hubbard-Stratonovich transformation Hubbard (1959); Stratonovich (1957) allows one to express two-body propagators as a high-dimensional integral over auxiliary fields of one-body propagators:


where the are any one-body operators. Applying this to we have


where we have introduced the vector of auxiliary fields , whose dimension, , is given by the number of possible -vectors. The operators are given by linear combinations of and . Zhang and Krakauer (2003); Suewattana et al. (2007)

Our focus in this paper is on the choice of pseudopotentials, which appear only in the one-body propagator . The handling of the two-body propagator and the implementation of the AFQMC phaseless approximation are unchanged from previous applications and can be found in Refs. [Zhang and Krakauer, 2003; Suewattana et al., 2007; Purwanto et al., 2009].

The overall computing cost in QMC depends not only on the computer time to execute a single time step for each random walker, but also on the statistical variance, which controls the size of the Monte Carlo sampling required to achieve a targeted statistical uncertainty (the QMC efficiency).Ma et al. (2015) The computing cost to execute a single imaginary-time step [Eq. (6)] in pw-AFQMC is proportional to , where is the number of planewaves. [The overall scaling is , where is the number of electrons in the simulation cell.] The statistical variance depends on the number of auxiliary fields . Reducing can therefore both reduce the computing time for each step in the random walk and increase the QMC efficiency. Convergence with respect to is controlled by the pseudopotential hardness,Zhang and Krakauer (2003); Suewattana et al. (2007) so that soft accurate pseudopotentials can potentially lead to major improvements in pw-AFQMC.

ii.3 Pseudopotential

The pseudopotential appears in the electron-ion interaction, in Eq. (3c). In second-quantized form, the pseudopotential’s action is safely isolated in the planewave matrix elements of its local [] and non-local [] parts, exactly as in DFT planewave methods.Yin and Cohen (1982) Non-local potentials thus present no difficulties in AFQMC (unlike in the real-space-based DMC method Mitas et al. (1991)). The planewave matrix elements of the local part of the pseudopotential are given by


where is the position of atom in the simulation cell, and is the Fourier transform of the (spherical) local part of the atomic pseudopotential. For single-projector NCPPs, the non-local part of the atomic pseudopotential is expressed by the separable Kleinman-Bylander Kleinman and Bylander (1982) form,


where, for each partial wave (e.g., for transition elements), there is only one projector. (The pseudopotential and pseudo-orbital are both functions of radial distance only, and is the usual spherical harmonic function.) The matrix elements of the non-local part of the pseudopotential can then be expressed in a separable form as,


where , and


where is obtained from the Bessel transform of the projector . The separable form of greatly simplifies and speeds up the use of the NCPPs, just as in DFT methods.

Equation (11) can be abbreviated as , in which is the overlap between pseudo-wavefunction with constructed projector . Having only one projector for each partial wave limits the energy range over which an NCPP can reproduce the scattering properties of the all-electron potential, which reduces its transferability. Hamann generalized Eq. (11) for optimized multiple projectors (in practice, implemented for two projectors per partial wave). Hamann (2013) Written in a diagonal representation, the multiple-projector pseudopotential can be compactly expressed for atom and partial wave as Hamann (2013)


Implementing the ONCVPSP in this form requires only minor modifications in pw-AFQMC, compared to Eq. (11). The extended energy range over which the scattering properties are reproduced often allows smaller planewave with excellent transferability properties. Hamann (2013)

Iii Applications

We describe applications of pw-AFQMC with ONCVPSP in three systems, the ionic insulator NaCl, the strongly correlated metal Cu, and the recently discovered sulfur hydride high-T and high-pressure superconductors.

iii.1 Ionic insulator: NaCl

Figure 1: (Color online) NaCl DFT/LDA calculated EOS curves (fits to Murnaghan’s equation Murnaghan (1944)), comparing all-electron LAPW (green dashed line), PAW (red dot-dashed line) and ONCVPSP (blue solid line). Curves are shifted to have the same minimum energy. The experimental lattice constant is indicated by the black dotted vertical line.

NaCl is a typical ionic compound, which crystallizes in the fcc structure. Due to the large overlap of Na valence electron and semicore and states, care must be used in the choice of pseudopotentials. Relaxation of the semicore states can significantly affect valence electron and hence material properties. Neglecting these effects, for example, by treating the and electrons as core states will give underestimation of the lattice constant and overestimation of the bulk modulus of Na in DFT calculations. Pseudizing instead the and states greatly reduces the discrepancy in lattice constant, to using the local density approximation (LDA) exchange-correlation functional. Calculated equation of states (EOS) with LDA are shown in Fig. 1. ONCVPSP results are compared to those from the all-electron linearized augmented planewave (LAPW) and the projector augmented wave (PAW) methods, using ELKelk () and ABINITGonze et al. (2002), respectively. ONCVPSPs were generated with Hamann’s open source pseudopotential code. onc () Both ONCVPSP and PAW results are in excellent agreement with LAPW. The agreement can be further quantified, using the factor, which was recently introduced by Lejaeghere et al.Lejaeghere et al. (2014) for comparing two EOS curves, and . Aligning the minimum energies, the definition of is:


for a volume range . (A typical choice of is around the equilibrium volume.) The factors are  meV and  meV for ONCVPSP and PAW calculations, respectively. The Na and Cl multiple-projector ONCVPSP pseudopotentials required kinetic energy cut-offs of only  Ry, much softer than for a single-projector norm-conserving pseudopotentials, which would have required  Ry.

Figure 2: (Color online) NaCl EOS calculated from pw-AFQMC (filled red circles with statistical error bars) using the same ONCVPSPs as in Fig. 1. For comparison, the DFT EOS in Fig. 1 is reproduced (blue solid line, energy shifted for convenient display). The vertical blue dashed and black dotted lines indicate the DFT and experimental equilibrium lattice constants , respectively. The pw-AFQMC calculated is indicated by the vertical red arrow with horizontal error bar indicating the uncertainty from a fit of the statistical data to Murnaghan’s equation. Murnaghan (1944)

Figure 2 shows the calculated pw-AFQMC NaCl EOS. A four-formula cubic simulation cell was used with twist boundary condition corresponding to the special -point (, , ); one- and two-body finite-size errors were reduced, using a post-processing finite-size correction scheme. Kwee et al. (2008); Ma et al. (2011) Here and throughout the rest of this paper, Trotter errors from Eq. (7) are removed by either extrapolation to the limit or choosing sufficiently small time-step values. Our pw-AFQMC calculations used LDA-generated trial wave functions. The discrepancy of DFT/LDA with experiment is essentially eliminated by the many-body calculations. The equilibrium lattice constant and bulk modulus calculated from pw-AFQMC,  bohr and  GPa, are in excellent agreement with the experimental values,  bohr and  GPa. Schimka et al. (2011)

iii.2 Transition metal: fcc Cu

Transition metal materials have played a central role in the study of strongly correlated physics, and copper based systems have especially attracted a great deal of attention. Pickett (1989) Ab-initio many-body calculations for transition metal systems have been very limited, Foyevtsova et al. (2014) and most previous calculations have relied on DFT or related approaches. In this subsection, we present many-body pw-AFQMC results on fcc copper, a prototypical correlated metal.

For good transferability, a frozen neon-core Cu pseudopotential is required, retaining the states. Single-projector NCPPs are challenged in this regard, because the and scattering properties near the Fermi energy depend on projectors constructed at much lower energies from the semicore and states. Even the scattering properties near are difficult, due to the resonant nature of scattering. To maximize the accuracy, very hard single-projector NCPPs must be used, with large planewave  Ry. This is alleviated by the multiple-projector ONCVPSP. Projectors for and can be constructed using both the semicore and and higher-lying valence or virtual and states. Similarly, two reference energies can be used to closely reproduce the all-electron scattering. We used  Ry and radial cutoffsHamann (2013) of  bohr for , respectively. The projectors were constructed using the ONCVPSP code,onc () with the LDA exchange-correlation functional. The multiple-projector pseudopotential yields very good agreement with all-electron LAPW results at the DFT level, giving factor  meV, as shown in Fig. 3. The non-parallellity error of  mRy in the computed EOS is smaller than the targeted statistical resolution of the QMC calculations which we discuss next.

Figure 3: (Color online) Cu DFT/LDA EOS comparison of ONCVPSP and LAPW. The EOS are shifted to have the same minimum energy. The inset shows the energy difference versus lattice size.
Figure 4: (Color online) Cu EOS calculated by pw-AFQMC using the same ONCVPSP as in Fig. 3. AFQMC results are shown by filled symbols, with statistical error bars indicated. For comparison, the DFT EOS with ONCVPSP is reproduced from Fig. 3 (blue solid line, energy shifted). The vertical blue dashed and black dotted lines indicate the DFT and experimental equilibrium lattice constants , respectively. The pw-AFQMC calculated is indicated by the vertical red arrow with horizontal error bar indicating the uncertainty from a fit to Murnaghan’s equation. Murnaghan (1944)

Figure 4 shows the calculated pw-AFQMC Cu EOS. A four-atom cubic simulation cell was used. Because Cu is metallic, twist-averaging with a Monkhorst-Pack (MP) -point grid Monkhorst and Pack (1976) was applied. Small random distortions were applied to each of the -points to lift band degeneracy in the trial wave function. Additionally, post-processing one- and two-body finite-size error corrections Kwee et al. (2008); Ma et al. (2011) were applied. The residual finite-size error is not expected to affect the EOS around equilibrium significantly. This was verified with the following approximate estimate which helps to avoid many computationally costly QMC tests. Calculations with up to of primitive unit cell were carried out using the LDA+ method. The DFT+ method includes a mean-field treatment of on-site 3d electron-electron interactions on the Cu atoms. This effect is absent in standard DFT local and semilocal exchange-correlation functionals, which are based on electron gas calculations. Since the choice of is largely determined by experience and by systematic benchmarking, multiple effective values of , from to , are studied in the simulations. The same twist-averaging and post-processing finite-size techniques were applied to the LDA+ results. The equilibrium lattice constant and bulk modulus did not change up to the largest test simulation cells. The final calculated pw-AFQMC EOS in Fig. 4 yields equilibrium lattice constant and bulk modulus,  bohr and  GPa, which are in excellent agreement with experimental values  bohr and  GPa (zero-point effects removed) Schimka et al. (2011). The accuracy of the ONCVPSP compared to all-electron LAPW results at the DFT level, without partial core corrections (see Sec. IV), is thus seen to be a good predictor of its transferability at the pw-AFQMC many-body level.

iii.3 Sulfur hydride high- high-pressure superconductor: HS

In this section, we present benchmark pw-AFQMC calculations on two candidate structures for high-temperature, high-pressure superconductivity in the sulfur hydride system. Applying the multiple-projector pseudopotentials, we test DFT/GGA predictions of the structural energetics of HS and HS by comparison with many-body AFQMC results.

Since Ashcroft proposed that metallic hydrogen should exhibit superconductivity with ,Ashcroft (1968) there have been many investigations of prospective high- materials incorporating hydrogen, with a recent focus on hydrides, where reduced metallization pressures are expected. Ashcroft (2004) Recent theoretical predictionsLi et al. (2014); Duan et al. (2014) of unusually high in sulphur hydrides under high pressure were subsequently supported by experiment. Drozdov et al. (); Drozdov et al. (2015); Troyan et al. (2016); Einaga et al. (2016) Measurements of resistivity and magnetic susceptibility indicate superconducting temperatures as high as  K at pressures  GPa; Drozdov et al. (2015) this was attributed to the HS phase. A novel experiment reported Meissner effect measurements that qualitatively confirmed the finding. Troyan et al. (2016) Subsequent DFT-based calculations have led to similar conclusions regarding the central role of electron-phonon coupling in driving the superconducting transition. Bernstein et al. (2015a); Duan et al. (2015); Bernstein et al. (2015b); Papaconstantopoulos et al. (2015); Errea et al. (2015); Komelj and Krakauer (2015) These calculations support the view that the sulfur hydrides are conventional superconductors, which are well described by Bardeen-Cooper-Schrieffer (BCS) theory Bardeen et al. (1957) with strong electron-phonon coupling leading to high . This is unlike the previously known high- cuprate and iron-based superconductors, where strong electron-electron interactions are believed to play a key role, although the superconducting mechanism has not yet been established. With a  K,Drozdov et al. (2015) hydrogen sulfide is one of the highest temperature superconductors on record, although extremely high pressures are required. Their discovery has re-energized the search for new superconductors in hydrogen-based and related materials.

Little is known experimentally regarding the high-pressure stability of hydrogen sulfide compounds. There has therefore been a strong reliance on standard DFT calculations, which have examined the high pressure phase stabilities and structures of HS. Bernstein et al. (2015a); Duan et al. (2015); Errea et al. (2015). The HS structure (space group No. ) has been a leading candidate for the stoichiometry that leads to highest . Other stoichiometries like HS are predicted to have competitive but less favorable enthalpies. It is important, therefore, to test these predictions with accurate many-body calculations. Here, we focus on candidate structures for two compositions, HS and HS, and compare pw-AFQMC results of their structural energetics with DFT/GGA predictions.

Figure 5: (Color online) HS () EOS calculated with ONCVPSP, single-projector NCPP, and USPP pseudopotentials, plotted as (minimum energies aligned). OPIUM generated NCPP EOS are shown for a range of . The inset shows the the actual EOS.

The ONCVPSPs of H and S were generated with  Ry. The projectors for S used unbound scattering states.Vanderbilt (1990); Hamann (2013) Figure 5 compares calculated ONCVPSP EOS with ultrasoft (USPP) and single-projector NCPP pseudopotential calculations, using the DFT GGA/PBE xc functional. The NCPPs were generated with the OPIUMWalter () package for several values of , and USPP “Standard Solid State Pseudopotentials” (USPP-SSSP) are adopted.Lejaeghere et al. (2016) ONCVPSP is in excellent agreement with USPP over a wide volume range (). The difference is less than 0.5 mRy per formula unit. Using USPP-SSSP as the reference, is  meV for ONCVPSP over the typical choice of , around the equilibrium volume of . The volume range in Fig. 5, of , is much wider, covering a span and including the superconducting volume near at transition pressure 200 GPa. For this volume range, the is  meV for ONCVPSP. By comparison, the single-projector NCPPs have values of , , , , and  meV, for values of 50, 60, 70, 80, and 100 Ry, respectively. To achieve comparable accuracy with the ONCVPSP, the NCPP requires a  Ry, which gives a nearly three times larger planewave basis.

(eV/atom) (eV/atom)
AFQMC -0.086(7) 0.111(5)
LDA ONCVPSP -0.082 0.058
OPIUM(100Ry) -0.083 0.056
USPP -0.084 0.056
PBE ONCVPSP -0.086 0.082
OPIUM(100Ry) -0.088 0.080
USPP -0.086 0.080
PBEsol ONCVPSP -0.083 0.060
OPIUM(100Ry) -0.084 0.058
USPP -0.084 0.059
PBE0 ONCVPSP -0.082 0.077
Table 1: Calculated pw-AFQMC structural energy difference for HS and HS, using ONCVPSP, compared to DFT-based calculations for four functionals and three pseudopotentials. For each of the four crystal structures, the fully relaxed structure from ONCVPSP-PBE was used for the pw-AFQMC calculations and for the other DFT functionals and pseudopotentials.
Figure 6: (Color online) Electron density distributions in HS (top panel) and HS (bottom panel) as a function of , computed from DFT/PBE with ONCVPSP. The main plots show the difference between the two space group structures for each composition, while the insets show the actual distributions of the two structures, with one shown as negative. =0.02 is chosen as the size of histogram bin.

For the pw-AFQMC benchmarks, we selected 200 GPa structures for two compositions, guided by the DFT/PBE results of Mazin et al.Bernstein et al. (2015a) For HS, the space group structures (#26) and (# 64) were selected, both with 12-atom primitive cells. For HS, the space group structures (#66) and (#229) were selected, with 32- and 4-atom primitive cells, respectively. We first calculated DFT-based energy differences for HS and for HS. Our ONCVPSP DFT/PBE calculations are in very good agreement with the results in Ref. Bernstein et al., 2015a. DFT-based results are shown in Table 1 for combinations of three pseudopotentials (ONCVPSP, NCPP-100Ry, and USPP) and four DFT exchange-correlation functionals: LDA, PBE, PBEsol, and the hybrid PBE0 method. (Note that the sign of the energy differences does not reflect the relative structural stabilities. For example, the calculated HS DFT/PBE enthalpy is actually lowest,Bernstein et al. (2015a) making it the most stable structure at  GPa.) To facilitate comparisons, the fully relaxed 200 GPa crystal structure from ONCVPSP-PBE was used for the other functionals and pseudopotentials and for the pw-AFQMC calculations. For HS, is nearly independent of the choice of DFT functional, while for HS, varies between 0.056 to 0.082 eV/atom.

For the pw-AFQMC calculations, 24-atom simulation cells were used for HS, doubling the size of the primitive unit cell in each structure. For HS, 32-atom simulation cells were used ( for ). Twist-averaged boundary conditions with a MP grid were used. One- and two-body finite-size corrections Kwee et al. (2008); Ma et al. (2011) were then applied to the many-body results. The pw-AFQMC energy differences are also shown in Table 1. The pw-AFQMC HS energy difference, 0.111(5) eV/atom, is nearly twice that given by the LDA and PBEsol, and about 50% larger than those from PBE and PBE0, while in HS the DFT-based calculations are identical with pw-AFQMC to within its statistical uncertainty.

To understand how the better agreement in HS arises compared to HS, we investigated the electron-density distributions for each composition. The result is illustrated in Fig. 6, which plots the densities calculated from ONCVPSP DFT/PBE for the four structures on the real-space FFT grid. In both HS and HS, the distribution is largely concentrated in the high-density region . The HS composition structures, however, show larger differences, especially in the range , than for the two HS structures. This provides a possible explanation of the better agreement of the different DFT functionals for HS than for HS. Similarly, it indicates that there will be better cancellation of electron correlation effects in HS, resulting in better agreement between DFT and pw-AFQMC many-body results.

The pw-AFQMC benchmarks in Table 1 show that DFT-based predictions are semi-quantitatively correct. The DFT predictions of HS and HS enthalpy differences could be off by the order of 30 meV and 50 meV in PBE and PBEsol, respectively. However, the stabilities are dominated by independent-electron contributions to the enthalpy, which are significantly larger than these differences. This suggests that the predictions on phase stabilities and structures from recent DFT studies are likely reasonable.

Iv Discussion

The applications above show that the use of multiple-projector ONCVPSP can greatly reduce the planewave basis size in pw-AFQMC many-body calculations, while maintaining or improving accuracy compared to single-projector NCPPs. ONCVPSP uses two projectors per partial wave in our applications, which maintains fidelity to scattering properties at reduced . As discussed in Section II.2, this results in significant reductions of the computational cost, both by reducing the computing time for each step in the random walk and, at the same time, increasing the QMC efficiency because of a smaller number of auxiliary fields in Eq. (8). For example, accurate results were obtained in fcc Cu with  Ry, in contrast to an estimated value of  Ry with NCPP.

It is important to note, however, that improvement in performance in DFT calculations by the ONCVPSP over single-projector NCPP does not always correlate with improvement in QMC. There are fundamental differences in the role of DFT-generated pseudopotentials when applied in a many-body context, versus in DFT. Clearly, when core-valence (or core-core) correlation effects are non-negligible, the use of pseudopotentials generated from an independent-electron approach can incur errors in many-body calculations. This is not the case for the systems treated in this paper. For example, in NaCl, small-core pseudopotentials are taken to pseudize the , states in both DFT and AFQMC. In DFT, partial-core effects were negligible as shown by the good agreement with LAPW in Fig. 1. Similarly in AFQMC, excellent agreement is found with experiment in Fig. 2.

A case that illustrates the difference is bulk Si, where the Ne-core NCPP causes a pseudopotential error both in DFT and AFQMC. In DFT, this can be remedied using a partial-core correction, which is not available in AFQMC. One way to remove the Ne-core error in AFQMC is with the frozen-core (FC) approximation. Purwanto et al. (2013); Ma et al. (2015) A He-core pseudopotential is used to generate DFT and orbitals in the crystalline solid environment. After a unitary rotation to the Kohn-Sham basis, the and orbitals are frozen Ma et al. (2015). The corresponding FC Hamiltonian, which incorporates an effective Ne-core pseudopotential, was shown to yield excellent results Ma et al. (2015) (in which the NCPP for the He-core pseudopotential used in the DFT calculations to generate the FC orbitals and Kohn-Sham basis had an extremely high cutoff, 600 Ry). To further study the implicit treatment of core-valence interactions in the FC approximation, we repeated this procedure with a much softer He-core ONCVPSP ( 64 Ry).

At the DFT level, this ONCVPSP works as well as the 600 Ry NCPP. They both capture (treating 12 electrons/Si) the core-valence corrections and yield excellent agreement with all-electron LAPW and with partial-core-corrected Ne-core pseudopotential calculations. However, the corresponding FC AFQMC calculation is not improved, and a pseudopotential error is seen as in the Ne-core calculation. This is because the softer He-core ONCVPSP has large pseudizing radius cut-offs which affect the resulting and orbitals. In DFT, such errors can be partially recovered because the densities are properly compensated for. In QMC, when the less accurate and orbitals are frozen (treated at the Hartree-Fock level), the errors propagate into the FC many-body Hamiltonian that the QMC treats, which cannot be corrected.

A good indicator of the accuracy of ONCVPSPs in many-body calculations is thus good core-valence separation and good DFT performance of ONCVPSP (without partial-core corrections) compared to all-electron calculations. The improved ONCVPSP scattering properties and transferability then allow smaller values of , which can significantly reduce the many-body computing cost while retaining high accuracy. In intermediate cases, such as in Si, when partial-core corrections are necessary in DFT calculations, more care is required.

V Summary

We have successfully implemented the multiple-projector ONCVPSP into the many-body pw-AFQMC method. The accuracy is demonstrated by calculations of bulk properties of NaCl and the more strongly correlated fcc Cu. With this technique, we also benchmarked the structure transition energy barriers in the recently discovered high-temperature superconductor sulfur hydride systems. In these systems, modest electron-electron correlation and large cancellation effects are seen in the energies between different structures, and we find that the estimations from DFT are in reasonable agreement with the many-body AFQMC results. The implementation of multi-projectors pseudopotential allows pw-AFQMC to treat systems with smaller pseudopotential errors and at significantly lower planewave energy cut-offs, and hence to reach larger and more complicated systems.

This work is supported by DOE (DE-SC0001303), NSF (DMR-1409510), and ONR (N000141211042). An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program, using resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. We also acknowledge computing support from the computational facilities at the College of William and Mary.


  1. W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  2. S. Zhang and H. Krakauer, Phys. Rev. Lett. 90, 136401 (2003).
  3. M. Suewattana, W. Purwanto, S. Zhang, H. Krakauer, and E. J. Walter, Phys. Rev. B 75, 245123 (2007).
  4. W. Purwanto, S. Zhang, and H. Krakauer, Journal of Chemical Theory and Computation 9, 4825 (2013), pMID: 26583401.
  5. F. Ma, W. Purwanto, S. Zhang, and H. Krakauer, Phys. Rev. Lett. 114, 226401 (2015).
  6. L. Mitas, E. L. Shirley, and D. M. Ceperley, The Journal of Chemical Physics 95, 3467 (1991).
  7. R. Hennig, A. Wadehra, K. Driver, W. Parker, C. Umrigar, and J. Wilkins, Physical Review B 82, 014101 (2010).
  8. S. Sorella, M. Casula, L. Spanu, and A. Dal Corso, Physical Review B 83, 075119 (2011).
  9. R. Nazarov, L. Shulenburger, M. Morales, and R. Q. Hood, Physical Review B 93, 094111 (2016).
  10. W. Purwanto, H. Krakauer, and S. Zhang, Phys. Rev. B 80, 214116 (2009).
  11. D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  12. D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
  13. A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos, Phys. Rev. B 41, 1227 (1990).
  14. M. Schlipf and F. Gygi, Computer Physics Communications 196, 36 (2015).
  15. K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, D. Caliste, I. E. Castelli, S. J. Clark, A. Dal Corso, et al., Science 351, 1415 (2016).
  16. L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425 (1982).
  17. M. T. Yin and M. L. Cohen, Phys. Rev. B 26, 3259 (1982).
  18. L. M. Fraser, W. M. C. Foulkes, G. Rajagopal, R. J. Needs, S. D. Kenny, and A. J. Williamson, Phys. Rev. B 53, 1814 (1996).
  19. H. F. Trotter, Proc. Am. Math. Soc. 10, 545 (1959).
  20. M. Suzuki, Communications in Mathematical Physics 51, 183 (1976).
  21. J. Hubbard, Phys. Rev. Lett. 3, 77 (1959).
  22. R. D. Stratonovich, Dokl. Akad. Nauk. SSSR 115, 1907 (1957).
  23. F. D. Murnaghan, Proc. Natl. Acad. Sci. USA 30, 244 (1944).
  24. The Elk full-potential linearized augmented-plane wave, available at (accessed May 1, 2011).
  25. X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux, M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete, G. Zerah, F. Jollet, et al., Computational Materials Science 25, 478 (2002).
  26. D.R. Hamann, Open-source pseudopotential code ONCVPSP, available at
  27. K. Lejaeghere, V. V. Speybroeck, G. V. Oost, and S. Cottenier, Critical Reviews in Solid State and Materials Sciences 39, 1 (2014).
  28. H. Kwee, S. Zhang, and H. Krakauer, Phys. Rev. Lett. 100, 126404 (2008).
  29. F. Ma, S. Zhang, and H. Krakauer, Phys. Rev. B 84, 155130 (2011).
  30. L. Schimka, J. Harl, and G. Kresse, The Journal of Chemical Physics 134, 024116 (2011).
  31. W. E. Pickett, Rev. Mod. Phys. 61, 433 (1989).
  32. K. Foyevtsova, J. T. Krogel, J. Kim, P. R. C. Kent, E. Dagotto, and F. A. Reboredo, Phys. Rev. X 4, 031003 (2014).
  33. H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  34. N. W. Ashcroft, Phys. Rev. Lett. 21, 1748 (1968).
  35. N. W. Ashcroft, Phys. Rev. Lett. 92, 187002 (2004).
  36. Y. Li, J. Hao, H. Liu, Y. Li, and Y. Ma, The Journal of Chemical Physics 140, 174712 (2014).
  37. D. Duan, Y. Liu, F. Tian, D. Li, X. Huang, Z. Zhao, H. Yu, B. Liu, W. Tian, and T. Cui, Scientific reports 4, 6968 (2014).
  38. A. P. Drozdov, M. I. Eremets, and I. A. Trojan, arXiv:1412.0460.
  39. A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov, and S. I. Shylin, Nature 525, 73 (2015).
  40. I. Troyan, A. Gavriliuk, R. Rüffer, A. Chumakov, A. Mironovich, I. Lyubutin, D. Perekalin, A. P. Drozdov, and M. I. Eremets, Science (New York, N.Y.) 351, 1303 (2016).
  41. M. Einaga, M. Sakata, T. Ishikawa, K. Shimizu, M. I. Eremets, A. P. Drozdov, I. A. Troyan, N. Hirao, and Y. Ohishi, Nature Physics 12, 835 (2016).
  42. N. Bernstein, C. S. Hellberg, M. D. Johannes, I. I. Mazin, and M. J. Mehl, Phys. Rev. B 91, 060511 (2015a).
  43. D. Duan, X. Huang, F. Tian, D. Li, H. Yu, Y. Liu, Y. Ma, B. Liu, and T. Cui, Phys. Rev. B 91, 180502 (2015).
  44. N. Bernstein, C. S. Hellberg, M. D. Johannes, I. I. Mazin, and M. J. Mehl, Phys. Rev. B 91, 060511 (2015b).
  45. D. A. Papaconstantopoulos, B. M. Klein, M. J. Mehl, and W. E. Pickett, Phys. Rev. B 91, 184511 (2015).
  46. I. Errea, M. Calandra, C. J. Pickard, J. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Phys. Rev. Lett. 114, 157004 (2015).
  47. M. Komelj and H. Krakauer, Physical Review B 92, 205125 (2015).
  48. J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957).
  49. E. J. Walter, Opium: pseudopotential generation project,
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description