title: Equation of state of atomic solid hydrogen by stochastic many-body wave function methods
authors: Sam Azadi George H. Booth Thomas D. Kühne
text: 
  
   
  
  
  
  Department of Physics, King's College London, Strand, WC2R 2LS London, United Kingdom
  sam.azadi@kcl.ac.uk 
  
  
  Department of Physics, King's College London, Strand, WC2R 2LS London, United Kingdom
  
  
  Department of Chemistry, Paderborn Center for Parallel Computing, Paderborn University, 33098 Paderborn, Germany
  
  
  
   We report a numerical study of the equation of state of crystalline body-centered-cubic (BCC) hydrogen, tackled with a variety of complementary many-body wave function methods. These include continuum stochastic techniques of fixed-node diffusion and variational quantum Monte Carlo, and the Hilbert space stochastic method of full configuration-interaction quantum Monte Carlo. In addition, periodic coupled-cluster methods were also employed. Each of these methods is underpinned with different strengths and approximations, but their combination in order to perform reliable extrapolation to complete basis set and supercell size limits gives confidence in the final results. The methods were found to be in good agreement for equilibrium cell volumes for the system in the BCC phase, with a lattice parameter of 3.307 Bohr.
   
  
  Equation of state of atomic solid hydrogen by stochastic many-body wave function methods
      Thomas D. Kühne
      December 30, 2023
  ========================================================================================
  
  
  
  
  
  § INTRODUCTION
  
  
  A stochastic description of quantum mechanics has significant advantages in the understanding of quantum systems, especially when a large number of degrees of freedom are involved. The main advantage of this approach relies on the exploitation of well-established mathematical bounds derived from probability theory and stochastic processes to control the convergence of these properties. In this picture, quantum particles move along stochastic trajectories, and expectation values can be formulated as ensemble averages over the space of these trajectories. In practice, we exploit the similarity between the Schrödinger equation in imaginary time, which is a linear, parabolic partial differential equation, and the diffusion equation. The efficiency of Monte Carlo approaches relies on the use of random numbers to sample the 3N-dimensional phase-space space of configurations, where N is the number of variables <cit.>. 
  
  
  Quantum Monte Carlo (QMC) approaches to large, ab initio systems with realistic many-body Hamiltonians have provided some of the most accurate and reliable descriptions of both Fermionic and Bosonic quantum matter <cit.>. The most established QMC techniques include variational Monte Carlo (VMC) <cit.>, and diffusion Monte Carlo (DMC) <cit.>, have been successfully applied to a variety of realistic quantum materials, including chemical systems <cit.>, and solids <cit.>. The key advantages of these QMC methods is their efficient large-scale parallelization, low scaling with system size, and predictive power. At present, DMC is the most commonly used QMC technique for high-precision wave function-based predictions of material properties <cit.>. DMC works with an ensemble of (almost) independent 'walkers' or 'configurations', whose stochastically realised dynamics ensures that they explore the Fermionic configuration space of the system. Beginning with a trial wave function, samples of the ground state wave function of the system are projected onto by applying the imaginary time operator exp(-τH). The dominant open problem preventing the exact numerical calculation of many-electron systems by DMC is the Fermionic sign problem, which arises from the antisymmetric nature of many-body wave function with respect to electron exchange. A widespread practical solution for this problem is the so-called fixed-node (FN) approximation <cit.>. In this, a many-electron nodal surface is defined as the one coinciding with that given by a trial many-electron wave function. This trial state is a function of 3N variables and the trial nodal surface is the (3N-1)-dimensional hypersurface on which the function is zero, and across which the wave function amplitude changes sign. The approximation ensures that the hypersurface of the sampled wave function in DMC coincides with that of the trial wave function. This FN algorithm gives the lowest-energy many-electron state with the given nodal surface <cit.>, which ensures that it is a variational approximation. In principle, this approximation would be exact if the applied nodal boundaries coexisted with the exact nodal surface of the many-Fermionic wave function. In practice, however, the errors in FN DMC energies are usually about 5% of the correlation energy for commonly used trial states. 
  
  
  An alternative approach, to tackle the Fermionic sign problem, is to allow the exact nodal surface to be an emergent property of the underlying algorithm. Such a simulation would therefore not require initial information of a trial nodal surface. It was demonstrated that such property can be obtained in full configuration-interaction quantum Monte Carlo (FCIQMC) <cit.>. This technique can be applied to Fermionic systems after projection into a discrete basis set familiar to conventional quantum chemistry approaches. It was demonstrated that this method can converge to capture the complete correlation energy, or full configuration-interaction (FCI) solution, for the given basis set. There are a number of similarities and differences between the DMC and FCIQMC approaches to stochastic realisation of quantum systems. Both techniques run a long-time integration of the imaginary time Schrödinger equation. However, while DMC works in continuous real space, FCIQMC spans the Hilbert space of Slater determinants. In DMC walkers follow the diffusion equation, while in FCIQMC the propagation step is based on a fictitious population dynamics of creation and death processes. In DMC and FCIQMC, the wave function is rendered by walkers representing a specific configuration in their respective configurations, which enables the methods to stochastically sample the wave function without storing the exponential number of amplitudes in the space. The key step of the FCIQMC algorithm, which allows the nodal structure of the wave function to emerge, is walker annihilation. Since each walker has a defined sign (or phase for complex-valued wave functions), if two configurations with opposite signs simultaneously occupy the same determinant, both walkers are deleted from the simulation. The walker annihilation mechanism has also been explored in DMC and Green's function Monte Carlo <cit.>.
  
  
  Solving realistic many-body Hamiltonians is the main challenge in condensed matter physics and quantum chemistry. Traditional quantum chemical methods, including coupled-cluster (CC) theory <cit.> and configuration-interaction <cit.>, have been developed for solving the Schrödinger equation, primarily for chemical systems described by a given one-particle basis set. These approaches truncate the wave function to a specific ansatz, which can be relaxed to define a systematic hierarchy of approximations to exactness. Due to the computational complexity of these quantum chemical methods, applying this systematic hierarchy of methods for extended systems and solids is in its early stages of research<cit.>. The application of FCIQMC to realistic solids also demonstrated a promising route for providing reference ground state many-electron energies to benchmark quantum-chemical techniques, including the CC ansatz. On the other hand, many alternate approaches for dealing with the high dimensionality of real extended systems have been developed, including local truncation, single-particle Green's function methods, novel Monte Carlo algorithms, and embedding techniques, all of which also benefit from comparison to higher accuracy approaches rather than experiment <cit.>. 
  
  In this work, we compute the equation of state (EOS) of atomic solid hydrogen in the body centered cubic (BCC) lattice using CC, FCIQMC and DMC techniques. Studying the BCC phase of solid hydrogen is critical to understand the origin of metallic magnetism for this system. Many theoretical and numerical investigations have concentrated on this atomic phase of solid hydrogen<cit.>. Despite the fact that at low densities the BCC atomic phase is not the most stable structure of solid hydrogen<cit.>, investigations on this simple but realistic system can provide qualitative insights into metal-insulator transitions, and also supply a reference for commonly used mean-field approximations. This is why we consider the ground state energy of the paramagnetic BCC phase of solid hydrogen within the density regime 1.3<r_S<2.4, where the paramagnetic phase is more stable than antiferromagnetic and ferromagnetic states. We compare our FCIQMC results with CC singles-doubles (CCSD) calculations, which are obtained using the same basis set and system size. We then extend our investigation to include DMC to allow for larger system sizes and without a substantial basis set dependence. Since DMC results are not constrained by basis size, the DMC and FCIQMC values can not be directly compared due to their different model assumptions. However, by careful extrapolation of their intrinsic model assumptions, these complementary descriptions can be combined to provide accurate estimates of the correlation energy in the thermodynamic and complete basis set (CBS) limits.
  
  
  
  § COMPUTATIONAL METHODS
  
  We briefly describe the FCIQMC and DMC approaches. Both methods have been presented in more detail in previous works <cit.>, but less often alongside each other.
  
  
  
   §.§ FCIQMC Method
  
  
  In FCIQMC, we first choose a basis of 2M one-particle spin-orbitals ϕ_n, from which the space of all possible N-electron determinants can be constructed via 
  
  
      D_ i   =    D_i_1, ..., i_N = 1/√(N!) | ϕ_i_1⋯ϕ_i_n | 
         =   1/√(N!)ϕ_i_1(r_1)    ⋯   ϕ_i_1(r_N) 
      ϕ_i_2(r_1)    ⋯   ϕ_i_2(r_N) 
      ⋯   ⋯   ⋯
      ϕ_i_N(r_1)    ⋯   ϕ_i_N(r_N) 
      .
  
  
  With this set of determinants as an N-body basis, the wave function can be expanded as
  
  
      | Ψ⟩  = ∑ c_ i | D_ i⟩ ,
  
  
  with the optimal coefficients c_ i being the ones that variationally minimise
  
  
      E(c_ i) = ⟨Ψ|H|Ψ⟩/⟨Ψ|Ψ⟩ .
  
  
  These coefficients are found from the solution of the matrix eigenvalue problem 
  
  
      ∑_j H_ ij c_ j = E_0^FCI c_ i ,
  
  
  where H_ ij = ⟨ D_ i|H|D_ j⟩. This FCI approach captures all possible correlation energy within the basis of orbitals, and is therefore in principle systematically improvable as the basis is enlarged, albeit at exponential cost. Due to the basis incompleteness error, the calculated total energies are therefore often higher than FN-DMC energies, but cancellation of errors is in general more reliable in appropriately optimized quantum chemical basis sets. However, the total number of determinants increases exponentially with system and basis size, hence the appeal of a stochastic approach which can deal with such high-dimensional spaces. 
  
  The FCIQMC algorithm converges to the FCI (lowest energy) eigenvector of the FCI matrix eigenproblem of Eq. <ref>, via the solution of the imaginary-time Schrödinger equation. The master equation governing the stochastic walker dynamics is given by
  
  
      dc_ i(t)/dt   =    -∑_ j [H_ij -(E_HF +S) δ_ij] c_ j(t) 
         =    -∑_j (K_ij - Sδ_ij) c_ j(t) .
  
  
  The amplitude of each determinant is then coarse-grained with a variable resolution, with the vast majority of amplitudes then represented at any single FCIQMC iteration by a zero amplitude. By representing the wave function at any single iteration just by the walkers (necessarily with non-zero weights) rather than the amplitudes, there is a significant compression of the wave function storage requirement to describe this snapshot of the wave function, which is formally decoupled from the size of the underlying Hilbert space. The success of the algorithm therefore relies on the ability to perform this compression, whilst maintaining a faithful realisation of the state, ultimately exploiting its inherent sparsity afforded by the choice of representation. The (signed) population of walkers then sample the configuration space q_ i, updating each iteration depending on a choice of time step Δ t, through a series of `spawning' steps onto connected configurations and `death' steps, which generally reduce the local population of walkers.
  
  These spawning and death steps stochastically update the amplitude on a determinant via
  
      q_ i^(n+1) = q_ i^(n) - Δ t ∑_ j (K_ ij -Sδ_ i,j) q_ j^(n) ,
  
  where q_ i is the random variable denoting the instantaneous walker population on configuration |D_ i⟩.
  
  These walker dynamics are not dissimilar to a first-order approximation to those of DMC, where the propagator is G̃ = I - (H - SI) Δ t instead of G = exp(-(H-SI)Δ t), but without importance sampling and the FN approximation.
  As long as Δ t ≤ 2/(E_max - S), FCIQMC will formally yield the exact ground sate without any time-step error.  The energy can be extracted as 
  
  
      E(t)    =   ⟨exp(-tĤ) D_0|Ĥ|D_0⟩/⟨exp(-tĤ)D_0|D_0⟩
         =    E_HF + ∑_ j≠ 0c_ j(t) ⟨ D_ j|Ĥ|D_0⟩/c_0(t)
         =    E_HF + ⟨∑_ j≠ 0 q_ j(t) H_ j0⟩/⟨ q_0(t) ⟩ ,
  
  
  where |D_0 ⟩ represents a trial state <cit.>.
  
  In keeping with all general Fermionic QMC techniques, there is also a sign problem that afflicts FCIQMC <cit.>. In order to ensure that the annihilation is sufficient in order to allow the exact nodal structure of the FCI wave function to emerge, the number of walkers must be large enough, and while this number is generally much smaller than the dimensionality of the entire Hilbert space, it still grows exponentially with system size. To reduce the number of walkers required for high accuracy, the FCIQMC algorithm can be improved with the `initiator' approximation <cit.>. Initiator FCIQMC is a systematically improvable approximation <cit.>, whereby the determinants are divided into two classes, labelled `initiator' and `non-initiator'. Initiator determinants are allowed to create new walkers on unoccupied determinants. However non-initiator determinants are only allowed to create new walkers on already occupied determinants. The label of initiator or non-initiator is chosen depending on the current population of walkers on any determinant, with the threshold n_ add used to determine this label. This dynamic adaptation helps the walker population to stabilize a sign structure at far lower walker numbers, at the cost of introducing a small systematic error. This error can be systematically reduced as the number of walkers is increased, which enlarges the set of initiators. Other adaptations in recent years have also improved the accuracy and scope of the method, including the computation of excited states<cit.>, unbiased molecular properties <cit.> and its use as a solver within an active space framework<cit.>. 
  
  
  
   §.§ Diffusion Monte Carlo
  
  
  The DMC method is a stochastic technique for many-electron systems with a much longer history than FCIQMC. Any solution of the time-dependent Schrödinger equation can be expanded in the eigenfunctions of the Hamiltonian Φ_n( R), i.e.
  
  
      Ψ ( R, τ) = Σ_n c_n Φ_n ( R) exp(-E_nτ),
  
  
  where E_n is the eigenvalue corresponding to Φ_n and τ = it. One can find that Ψ ( R, τ→∞) ≃ c_0 Φ_0 ( R) exp(-E_0 τ), which is proportional to the ground state wave function. In principle, the Schrödinger equation can be solved by propagating an arbitrary wave function in imaginary time for long enough. This propagation can be seen in the differential form of the imaginary-time Schrödinger equation
  
  
      -1/2∇^2 Ψ ( R, τ) + [V( R) -E_T] Ψ( R, τ) 
         =    -∂/∂τΨ( R, τ),
  
  
  where ∇^2 = Σ_i ∇_i^2 acts over all coordinates within the vector R, whereas E_T is a constant energy offset. This equation is equivalent to a diffusion equation, in which Ψ ( R) represents the density of particles at point R. The particles diffuse with a diffusion coefficient D = 1/2, and are absorbed with rate V( R) - E_T. Assuming Ψ ( R, τ) is a probability density, we distribute an initial set of walkers with probability density given by Ψ ( R, τ). The walkers then diffuse and can be removed or created accordingly. This can be simulated via a stochastic process, whereby in the limit τ→∞, the walkers would be distributed according to the ground sate wave function. 
  
  For a Fermionic system, a wave function Ψ( R) must have both positive and negative regions to be antisymmetric with respect to particle exchange. Hence, it can not be used as a probability density. This problem can be overcome by using a guiding function f( R, τ) = Ψ_T( R) Ψ( R, τ).  Provided Ψ_T( R) and Ψ( R, τ) have the same nodal surface, f( R, τ) has the same sign over all configuration space, and can be interpreted as a probability distribution function. By multiplying both sides of the imaginary-time Schrödinger equation by Ψ_T, we obtain
  
  
      -∂ f/∂τ = -1/2∇^2 f +∇· [f v] + [E_L - E_T] f,
  
  
  where v( R) = ∇Ψ_T ( R)/Ψ_T( R) and E_L( R) = ĤΨ_T( R)/Ψ_T( R). If we consider f as a probability distribution, this equation is the Fokker-Planck equation describing the diffusion of non-interacting classical particles, with an imposed drift velocity v( R, τ) and absorption coefficient [E_L( R) - E_T]. We can therefore distribute a set of particles according to an initial distribution f( R, 0) = |Ψ_T( R)|^2, and let them evolve according to the Fokker-Planck equation. In the limit of τ→∞, the walkers will be distributed according to the minimal energy wave function with the same nodal surface as the trial wave function<cit.>. The quality of the nodal surface of the trial wave function therefore determines the error of DMC in an uncontrolled way. In practice, the trial wave function is optimised by VMC before being used in DMC. Generally, lower VMC energy imply a better nodal surface <cit.>. Further details of the implementation of the DMC algorithm are discussed in Ref. casino. 
  
  
  
   §.§ Simulation setup
  
  
  FCIQMC calculations were performed using the NECI<cit.> package. The periodic Hartree-Fock (HF) calculations, from which the single particle orbitals were extracted, and CCSD simulations were carried out using PySCF <cit.> with norm-conserving pseudopotentials <cit.>. Gaussian basis sets of SZV, DZVP, and TZVP quality were used for all basis set-based calculations<cit.>, with density fitting employed in the computation of two-electron matrix elements. For the FCIQMC simulations, up to 1.5×10^8 walkers were used for the most demanding calculations, at lower densities and larger simulation cells <cit.>.
  
  Our VMC and DMC calculations were performed using the CASINO QMC package <cit.> and a trial wave function of the Slater-Jastrow form was employed. The one-electron orbitals defining the Slater determinant, were extracted from density functional theory (DFT) calculations using the Quantum Espresso code <cit.>, with a DFT plane-wave cutoff of 5,000 eV. The norm-conserving DFT pseudopotential with the Perdew-Zunger parameterization<cit.> of the local density approximation was used for our DFT, VMC and DMC calculations. The Jastrow term J( R), which captures most of the dynamical correlation between electrons is a positive, symmetric, explicit function of interparticle distances and consisted of polynomial one-body electron-nucleus (en), two-body electron-electron (ee), and three-body electron-electron-nucleus (een) terms. Whereas the parameters of the Jastrow were optimized by variance minimization at the VMC level <cit.>, the Slater determinant was taken from DFT, but not reoptimized. The main approximation in the DMC results is the FN approximation, which can be improved by including backflow transformation in the trial wave function <cit.>. The BCC unit cell, which was used to build a supercell for all calculations, includes two hydrogen atoms located at the corner and centre of the cell that were fixed for all the investigated densities and only the lattice parameter was changed for each r_S. To minimize the time-step error, a small time step of τ = 0.005 a.u. was used in all the DMC simulations.
  
  
  
  § RESULTS AND DISCUSSION
  
  
  
  
   §.§ FCIQMC and CCSD
  
  
  We first consider the convergence of the FCIQMC correlation energy in a restricted simulation cell. This is given in Fig. <ref>, where representative convergence of the BCC crystalline hydrogen for two densities is given within a DZVP basis and a simulation cell of 16 hydrogen atoms and 80 orbitals. The total number of walkers in the calculation for each unit cell volume was grown in stages, to check the convergence of the energy estimator. The number of walkers required to achieve convergence with respect to the initiator error varied between the cell sizes, from 40 million walkers at compressed geometries, to 150 million walkers for more expanded geometries, where stronger static correlation effects are expected to be prevalent. Remaining systematic errors are expected to be sub-mHa per atom, resulting in confidence in its ability to benchmark other approaches at these restricted cell and basis sizes.
  
  
  
  
  We also consider the convergence of the energy with respect to the basis set size, considering both the DZVP (with 5 orbitals per atom) and the smaller SZV basis (with only one basis functions per atom). Table <ref> gives the corresponding FCIQMC energies of the system at different densities within a statistical uncertainty of around 0.1 mH/atom from a blocking analysis, with likely remaining systematic error also sub-mH/atom, giving confidence in the ability of the FCIQMC to be used as a finite-basis benchmark. We compare these results to those of CCSD, obtained using the pyscf simulation package. Table <ref> gives the CCSD energies of BCC atomic hydrogen with two atoms per primitive unit cell, which were obtained using SZV, DZVP, and larger TZVP basis sets and a 2×2×2 k-point mesh, at different densities. 
  The CCSD total energy at the CBS limit are also given by a three-point extrapolation of the HF energy (E_HF^CBS) via a exponential form, as well as the CCSD correlation energy (Δ E_CCSD^CBS), via the standard two-point (DZVP and TZVP) inverse-cubic form to the complete basis set (CBS) limit<cit.>.  
  The values of the CCSD total energy in the CBS limit given by E_tot^CBS = E_HF^CBS + Δ E_CCSD^CBS are reported in Table <ref> as a function of r_S. Comparing the FCIQMC-DZVP (Table <ref>) with the corresponding CCSD energies (Table <ref>) in the same basis set can quantify the systematic error in the latter. The close agreement of these results shows the ability of CCSD to recover the dominant correlated physics in this system, in particular for small r_S values, and underlines the potential of a composite scheme which will be used later, where the lower-cost CCSD can be used to correct the FCIQMC results to account for the effect of remaining basis set incompleteness (or alternatively, for the FCIQMC to correct for static correlation errors in the CCSD).
  
  
  
  
  
  
  
  
  
  
  From these results, we can now obtain the EOS of BCC atomic hydrogen for these methods, both for finite basis sets, and also for extrapolated basis results, within 2×2×2 simulation cells. This is shown in Fig. <ref> for the FCIQMC and CCSD methods, respectively. Agreement between these two methods is excellent at higher densities, where the single-reference nature of the CCSD ansatz is expected to perform well. 
  Around the equilibrium cell volumes, at r_S=1.6282 and r_S=1.8608 a.u, the difference between FCIQMC-SZV and CCSD-SZV energies are -0.97, and -2.4 mHa/atom, respectively.  However, for expanded cell volumes, at the r_S=2.3261 a.u., the correlation energy captured by FCIQMC increases significantly over the CCSD results for both basis sizes, with the FCIQMC-SZV energy being 9.8 mHa/atom lower than the CCSD-SZV value. This is anticipated, due to the increasing levels of stronger correlation effects present in the system, as evidenced by the increasingly multiconfigurational nature of the FCIQMC calculations.
  
  
  Using the larger DZVP basis set introduces a major energy gain within the FCIQMC and CCSD calculations. At the equilibrium cell volume, the difference between FCIQMC-DZVP and FCIQMC-SZV is 23.3 mHa/atom. Once again, the agreement between FCIQMC and CCSD at compressed cells is excellent, with the agreement getting worse as the cell expands. However, using the larger basis set, the discrepancy between the CCSD and FCIQMC values for more expanded geometries is much smaller, pointing to an overestimation of the relative importance of strong correlation effects in small basis sets, while larger basis sets are able to more effectively screen these strong correlations. This improved description of the screening available in the larger basis results in a qualitatively different shape to the EOS, reducing the compressibility of the system, and predicting an equilibrium volume of r_S=1.6282 a.u., corresponding to a BCC lattice parameter of 3.307 Bohr. 
  
  We were also able to conduct CCSD calculations using the more complete TZVP basis set. Employing this increased basis did not qualitatively change the EOS, with the energies only deviating from those using the DZVP basis by just over 1 mHa/atom across all considered densities. Furthermore, extrapolating the CCSD energies to CBS limit marginally lowers the total energies compared to the CCSD-TZVP results, with this incompleteness being more notable at larger cell volumes. However, the DZVP basis is still found to be accurate for all cell sizes. 
  
  
  
  
  
  
  
   §.§ Diffusion Monte Carlo
  
  
  The results of our VMC and DMC calculations at a cell volume of r_s=1.3956 a.u. with 16 atoms are shown in Table. <ref>. These show that the dominant error in the DMC values derives from the FN approximation, since the accuracy of the DMC is almost independent of the flexibility afforded by the Jastrow component of the trial wave function, which does not affect the nodal structure. Other systematic errors in DMC, such as time-step errors or non-local pseudopotential errors can be affected by the quality of the Jastrow component, with these results demonstrating that these errors are small. This suggests that while the importance of the two-body Jastrow is significant for the VMC energies, it may be more efficient to consider a trial wave function for DMC that has not included the 2-body Jastrow terms, which can result is a substantial speedup in the calculation. 
  
  The 1-body Jastrow terms are also found to be negligible in this system, with the 3-body terms also found to be less important to the DMC energy than the inclusion of backflow correlations.
  It is found that the backflow reduces the DMC energy by 0.8 mHa/atom, which is more than twice the energy gain originating from the inclusion of the 2-body Jastrow. The importance of backflow in high-density matter has been seen before, where in the homogeneous three-dimensional electron gas with r_S < 5 the effects of backflow were found to dominate over those introduced by three-body correlations<cit.>. Indeed, for this pure hydrogen system, the inclusion of all Jastrow terms and backflow transformation into the wave function results in only a ∼1.36 mHa/atom improvement in the DMC total energy compared to just using a single Slater determinant trial state. This energy gain is smaller than what is generally known as chemical accuracy, and so is unlikely to be worthwhile in general for pure hydrogen systems.
  
  
  
  
  The comparison of DMC, FCIQMC-DZVP, and CCSD-CBS energies in this reduced 16 atom simulation cell are shown in Fig. <ref>. We find a qualitative difference between the Fock space approaches of FCIQMC and CCSD compared to the real-space DMC approach. The former methods exhibit a sharper minimum around r_S=1.6 a.u., while the DMC EOS has a shallow minimum in the density window of 1.6 < r_S < 1.86. This discrepancy is not due to basis set or strong correlation effects, but rather the differing rates of convergence of these methods to the thermodynamic limit. It is clear that any comparison and agreement between these methods will have to ensure that the effects of these small supercell sizes are accounted for, which is considered in the following section.
  
  
  
  
   §.§ Finite size convergence
  
  
  Using a finite simulation cell to mimic the properties of a periodic system introduces finite-size (FS) errors, which is one of the main challenges in the application of many-body techniques to extended systems. These FS errors have a number of different origins, including the description of the kinetic energy, the periodic Ewald interaction, and exchange energy <cit.>. The convergence of these different terms to the thermodynamic limit can vary between methods and their representation, as indicated in the previous section. Therefore, careful control to mitigate the impact of finite simulation cells and ensure that all desired quantities are converged with respect to these errors is essential for reliable and comparable results. Here, we employ the standard FS extrapolation technique introduced by Ceperley and co-workers for reducing the FS errors<cit.>. Specifically, we employ a 1/N extrapolation form of
  
  
      E_MB, ∞∼ E_MB, N + α (E_0, ∞ - E_0, N) + β/N,
  
  
  where α and β are fitting parameters, and E_0 is the energy of system obtained via a single-particle, or mean-field approach. Therein, E_MB, N represents the many-body energy of the finite system of N interacting electrons. Using this form, we extrapolate the CCSD-SZV and DMC energies to the infinite system size limit. For the CCSD-SZV calculations, the E_0, ∞ value is approximated to be the HF energy of the system for that density, obtained using an 8×8×8 k-point mesh. For the DMC extrapolation, the local density approximation was used with a 24×24×24 k-point mesh to obtain E_0, ∞.
  
  The top panel of Fig. <ref> shows CCSD-SZV energies, which are calculated at four supercell sizes and different densities, up to a maximum supercell size of 4×4×4. 
  The bottom panel shows the extrapolation of the DMC energies at each density (computed with the Slater-Jastrow trial wave function including one- and two-body terms), up to the largest supercell size with 432 atoms. The final CCSD-SZV and DMC energies extrapolated to the infinite system size limit, estimated from Eq. <ref>, are given in Table <ref> along with the standard deviation in the fit. 
  
  We find that the finite size correction of the DMC results is always positive, regardless of density. Hence, the DMC energies of finite systems are lower than the DMC energy at the thermodynamic limit. In contrast, the CCSD FS correction lowers the energy of the system for the density range of r_S < 2.0, while is positive for lower densities. 
  However, the comparison between these approaches is fundamentally limited by the small basis size (SZV) of the CCSD energies used for the extrapolation, which was required in order to reach the large supercell sizes. While this basis is qualitative different to the CBS limit, we aim to use these results as a finite size correction to the more complete basis set results.
  
  
  
  
  
  
  To correct the finite-size extrapolated CCSD results to mitigate for this small basis, we can assume that the basis set error between SZV and DZVP is independent of supercell size, with these errors therefore being additive. This allows us to compute a final basis set correction, which can be applied to the extrapolated CCSD results at the SZV basis level. Furthermore, we can similarly compute the energetic correction between CCSD and FCIQMC in the DZVP basis sets, and also consider this as a correction for the correlated physics beyond the CCSD ansatz in the thermodynamic limit. Leveraging these strengths of these different methods, we are able to compute a final EOS for this composite approach.
  
  
  
  
  
  
   §.§ BCC atomic solid hydrogen equation of state
  
  
  Figure <ref> illustrates the final equation of state for BCC atomic solid hydrogen, as obtained from our FCIQMC, CCSD, and DMC calculations. The FCIQMC and CCSD results, calculated within a DZVP basis set, are shown for both a finite supercell of 16 hydrogen atoms, and for the extrapolation to the infinite system size limit. Moreover, the CCSD-CBS energies, obtained using a finite simulation cell of 16 hydrogen atoms, are also illustrated. The shown DMC energies have also been extrapolated to the thermodynamic limit. Encouragingly, all three many-body wave function-based methods (FCIQMC, CCSD and DMC), which operate under very different assumptions and approximations, give the same equilibrium density for this phase, which is found to be close to r_S=1.65 a.u., with the finite-size correction for the DMC energies qualitatively reducing the equilibrium cell size. This agreement arises despite the fact that FCIQMC and CCSD are built from a Fock-space orbital representation, while DMC operates in the continuum, with the nodal surface being built on the Kohn-Sham DFT determinant. Yet, all three methods agree that the BCC atomic solid hydrogen with two atoms per cell has a minimum ground state energy at lattice parameter 3.307 Bohr. 
  
  However, while the agreement with respect to the equilibrium geometry is good, it is clear that discrepancies between the approaches remain, with the 2×2×2 supercell CCSD and FCIQMC results being in closer agreement with the thermodynamic limit DMC results than their extrapolated counterparts. This discrepancy is likely to be due to the extrapolation to the thermodynamic limit of CCSD, which was calculated only for the small SZV basis sets, which have previously been shown to be rather different to the larger DZVP basis sets. Further investigation into this effect and its mitigation is likely to be a continuing research direction.
  
  
  
  § CONCLUSION
  
  We report a numerical study of the EOS for BCC solid hydrogen within the density range of 1.3<r_S<2.4 a.u., using FCIQMC, VMC, and DMC as many-body stochastic wave-function approaches, as well as the deterministic quantum chemical periodic CCSD method. We find that although sizable static correlation effects were captured at the FCIQMC level of theory for expanded cell sizes, at compressed geometries, however, these are less important and the results of CCSD and DMC agree rather well with those of FCIQMC. These complementary approaches can therefore be combined in order to mitigate their respective weaknesses, to provide a tractable route for the analysis of the electronic structure of extended systems. With careful consideration of basis set incompleteness and finite size errors, all methods demonstrated good agreement in the prediction of the equilibrium lattice parameter of 3.307 Bohr. However, discrepancies remain in the broader shape of the EOS, likely arising from the restricted basis sizes required for the thermodynamic limit extrapolation. Nevertheless, agreement between these methods, and demonstration of a viable composite scheme in these extended settings is likely a reliable approach to study further phases of these high-pressure hydrogenic materials.
  
  
  
  § ACKNOWLEDGEMENT
  
  S.A. gratefully acknowledges the generous allocation of computing time by the Paderborn Center for Parallel Computing (PC^2) on OCULUS and the FPGA-based supercomputer NOCTUA. G.H.B. has received funding from the Royal Society via a University Research Fellowship, as well as funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 759063. T.D.K. would like to the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (grant agreement no. 716142).
  
  The data that support the findings of this study are available from the corresponding author upon reasonable request.
  
  
  9
  Landau D.P. Landau, and K. Binder,
  "A Guide to Monte Carlo simulations in statistical Physics"
  Cambridge University Press (2014)
  Kalos  M. H. Kalos,
   Phys. Rev. 128, 1791 (1962) ; J. Com. Phys. 1, 257 (1966)
  Ceperley95 D. M. Ceperley, 
  Rev. Mod. Phys. 67, 279 (1995)
  Becca F. Becca, and S. Sorella,
  "Quantum Monte Carlo approaches for correlated systems",
  Cambridge University Press (2017)
  Matthew W.M.C. Foulkes, L. Mitas, R.J. Needs, and G. Rajagopal
  Rev. Mod. Phys. 73, 33 (2001)
  VMC W. L. McMillan, 
  Phys. Rev. 138, A442 (1965)
  Umrigar007 C. J. Umrigar, J. Toulouse, C. Filippi, et al.,
  Phys. Rev. Lett. 98, 110201 (2007)
  
  Anderson J.B. Anderson,
  J. Chem. Phys. 63 1499, (1975)
  Ceperley80 D. M. Ceperley and B. J. Alder, 
  Phys. Rev. Lett. 45, 566 (1980)
  Reynold P.J. Reynold, D.M. Ceperley, B.J. Alder, and W.A. Lester,
  J. Chem. Phys. 77, 5593 (1982)
  Grossman J.C. Grossman, and L. Mitas,
  Phys. Rev. Lett. 74, 1323 (1995)
  Ozone S. Azadi, R. Singh, and T. D. Kühne, 
  Int. J. Quantum Chem. 115, 1673 (2015)
  Benzene S. Azadi and T. D. Kühne, 
  Phys. Rev. B 97, 205428 (2018) 
  Kolorenc J. Kolorenč, and L. Mitas,
  Report on Progress in Physics 74, 026502 (2011)
  Dubecky M. Dubecký, L. Mitas, and P. Jurečka, 
  Chem. Rev. 116, 5188 (2016)
  SignSWF F. Calcavecchia, F. Pederiva, M. H. Kalos, and T. D. Kühne,
  Phys. Rev. E 90, 053304 (2014)
  Disulfide S. Azadi, and T. D. Kühne, 
  J. Chem. Phys. 146, 084503 (2017)
  Luke L. Shulenburger, and T. R. Mattsson,
  Phys. Rev. B 88, 245117 (2013) 
  Hydrogen S. Azadi, W. M. C. Foulkes, and T. D. Kühne,
  New J. Phys. 15, 113005 (2013)
  Wagner Lucas K. Wagner, 
  Phys. Rev. B 92, 161116(R) (2015)
  Antonietti S. K. Sahoo, et al., 
  Sci. Rep. 10, 5832 (2020) 
  AndersonFN J. B. Anderson, 
  J. Chem. Phys. 63, 1499 (1975) ;J. Chem. Phys. 65, 1421 (1976); Int. J. Quantum Chem 15, 109 (1979)
  Ceperley91 D. M. Ceperley,
  J. Stat. Phys. 63, 1237 (1991)
  GB009 G.B. Booth, A.J.W. Thom, and A. Alavi,
   J. Chem. Phys. 131, 054106 (2009)
  Booth10 G.H. Booth, and A. Alavi,
   J. Chem. Phys. 132, 174104 (2010)
  Cleland D. Cleland, G.H. Booth, and A. Alavi,
   J. Chem. Phys. 132, 041103 (2010)
  Booth11 G.H. Booth, D. Cleland, A.J.W. Thom, and A. Alavi,
   J. Chem. Phys. 135, 084104 (2011)
  Booth13 G.H. Booth, A. Grüneis, G. Kresse, and A. Alavi,
  Nature 493, 365-370  (2013)
  Arnow D. M. Arnow, M. H. Kalos, M. A. Lee, and K. E. Schmidt,
  J. Chem. Phys. 77, 5562 (1982)
  Ceperley84 D. M. Ceperley and B. J. Alder, 
  J. Chem. Phys. 81, 5844 (1984)
  Purvis G.D. Purvis III, R.J. Bartlett,
  J. Chem. Phys. 76, 1910 (1982)
  Bartlett81 R.J. Bartlett,
  Annual Rev. Phys. Chem. 32 359 (1981)
  Bartlett07 R.J. Bartlett, and M. Musial,
  Rev. Mod. Phys. 79, 291 (2007)
  Stanton J.F. Stanton, and R.J. Bartlett,
  J. Chem. Phys. 98, 7029 (1993)
  Cizek J. Čížek, 
  J. Chem. Phys. 45, 4256 (1966)
  Jeziorski B. Jeziorski, and H. J. Monkhorst,
  Phys. Rev. A 24, 1668 (1981)
  Sherrill1999 C. D. Sherrill, and H. F. Schaefer,
  "Advances in Quantum Chemistry"  34, edited by P.-O. Löwdin,
  Academic Press, New York, (1999)
  Muller C. Müller, B. Paulus, 
  Phys. Chem. Chem. Phys. 14, 7605 (2012)
  Marsman M. Marsman, A. Grüneis, J. Paier, and G. Kresse, 
  J. Chem. Phys. 130, 184103 (2009)
  Gruber T. Gruber, K. Liao, T. Tsatsoulis, F. Hummel, and A. Grüneis, 
  Phys. Rev. X 8, 021043 (2018)
  Pulkin20 A. Pulkin, and G. Kin-Lic Chan,
  Phys. Rev. B 101, 241113(R) (2020)
  Gao20 Y. Gao, Q. Sun, J. M. Yu, M. Motta, J. McClain, A. F. White, A. J. Minnich, and G. Kin-Lic Chan,
  Phys. Rev. B 101, 165138 (2020)
  McClain J. McClain, Q. Sun, G. Kin-Lic Chan, and T. C. Berkelbach,
  J. Chem. Theory Comput. 13, 1209  (2017)
  Booth16 G. H. Booth, T. Tsatsoulis, G. Kin-Lic Chan, and A. Grüneis
  J. Chem. Phys. 145, 084111 (2016)
  Setten M. J. van Setten, F. Caruso, S. Sharifzadeh, X. Ren, M. Scheffler, F. Liu, J. Lischner, L. Lin, J. R. Deslippe, S. G. Louie, C. Yang, F. Weigend, J. B. Neaton, F. Evers, and P. Rinke,
  J. Chem. Theory Comput. 11, 5665 (2015)
  H2SWF F. Calcavecchia, and T. D. Kühne,
  Eur. Phys. Lett. 110, 20011 (2015) 
  Lejaeghere K. Lejaeghere et al.,
  Science 351, aad3000 (2016)
  Motta M. Motta et al.,
  Phys. Rev. X 7, 031059 (2017)
  Williams K. T. Williams et al., 
  Phys. Rev. X 10, 011041 (2020) 
  Ceperley87 D. M. Ceperley and B. J. Alder, 
  Phys. Rev. B 36, 2092 (1987)
  Ashcroft89 N. W. Ashcroft, 
  Nature (London) 340, 345 (1989)
  Mao89  H. K. Mao and R. J. Hemley, 
  Science 244, 1462 (1989) 
  Barbee T. W. Barbee, III, A. Garcia, and M. L. Cohen, Phys. Rev.
  Lett. 62, 1150 (1989); Nature (London) 340, 369 (1989)
  Wang90 X. W. Wang, J. Zhu, S. G. Louie, and S. Fahy, 
  Phys. Rev. Lett 65, 2414 (1990)
  Natoli V. Natoli, Richard M. Martin, and D. M. Ceperley
  Phys. Rev. Lett. 70, 1952 (1993) 
  McMahonJ. M. McMahon, M. A. Morales, C. Pierleoni, and D. M. Ceperley,
  Rev. Mod. Phys. 84, 1607 (2012)
  PRL14 S. Azadi, B. Monserrat, W.M.C. Foulkes, and R.J. Needs,
  Phys. Rev. Lett. 112, 165501 (2014)
  HydrogenBCC F. Calcavecchia, and T. D. Kühne, 
  Zeitschrift für Naturforschung A 73, 845 (2018)
  KeLiao K. Liao, X. Li, A. Alavi, and A. Grüneis,
  npj Comput Mater 5, 110 (2019)
  PRB19 S. Azadi, and T.D. Kühne,
  Phys. Rev. B 100 155103 (2019)
  Petruzielo F. R. Petruzielo, A. A. Holmes, Hitesh J. Changlani, M. P. Nightingale, and C. J. Umrigar
  Phys. Rev. Lett. 109, 230201 (2012)
  Blunt15 N. S. Blunt, S. D. Smart, J. A. F. Kersten, J. S. Spencer, G. H. Booth, and A. Alavi,
  J. Chem. Phys. 142, 184107 (2015)
  JSpencer J. S. Spencer, N. S. Blunt, and W.M.C. Foulkes, 
  J. Chem. Phys. 136, 054110 (2012)
  Blunt2015 N. S. Blunt, S. D. Smart, G. H. Booth, and A. Alavi,
  J. Chem. Phys. 143, 134117 (2015)
  Blunt2017 N. S. Blunt, G. H. Booth, and A. Alavi,
  J. Chem. Phys. 146, 244105 (2017)
  Overy14 C. Overy, G. H. Booth, N. S. Blunt, J. J. Shepherd, D. Cleland, and A. Alavi,
  J. Chem. Phys. 141, 244117 (2014)
  Thomas15 R. E. Thomas, D. Opalka, C. Overy, P. J. Knowles, A. Alavi, and G. H. Booth,
  J. Chem. Phys. 143, 054108 (2015)
  Samanta18 P. Samanta, N. S. Blunt, and G. H. Booth,
  J. Chem. Theory Comput. 14, 3532 (2018) 
  Thomas2015 R. E. Thomas, Q. Sun, A. Alavi, and G. H. Booth,
  J. Chem. Theory Comput. 11, 5316 (2015)
  Anderson2020 R. J. Anderson, T. Shiozaki, and G. H. Booth,
  J. Chem. Phys. 152, 054101 (2020)
  Ceperley D.M. Ceperley, and L. Mitas
  "New methods in computational quantum mechanics", 1-38 (1995)
  Umrigar C.J. Umrigar, M.P. Nightingale, K.J. Runge, 
  J. Chem. Phys. 99, 2865 (1993)
  casino R. J. Needs, M. D. Towler, N. D. Drummond, P. López Ríos, and J. R. Trail
  J. Chem. Phys. 152, 154106 (2020). 
  Turbo K. Nakano, C. Attaccalite, M. Barborini, L. Capriotti, M. Casula, E. Coccia, M. Dagrada, C. Genovese, Y. Luo, G. Mazzola, A. Zen, and S. Sorella1,
  J. Chem. Phys. 152, 204121 (2020)
  NECI K. Guther et al., 
  J. Chem. Phys. 153, 034107 (2020)
  pyscf1 Q. Sun et al
  J. Chem. Phys. 153, 024109 (2020)
  pyscf2 Q. Sun et al
  WIREs,  https://doi.org/10.1002/wcms.1340 (2017)
  gthPP S. Goedecker, M. Teter, and J. Hutter,
  Phys. Rev. B 54, 1703 (1996)
  Vondele J. VandeVondele, and J. Hutter,
  J. Chem. Phys. 127 114105 (2007)
  Pitarke R. Gaudoin and J. M. Pitarke, 
  Phys. Rev. Lett. 99, 126406 (2007)
  QE P. Giannozzi et al.,
  J. Phys.: Condens. Matter 21, 395502 (2009). 
  lda J. P. Perdew, and A. Zunger, 
  Phys. Rev. B 23, 5048 (1981).
  PRB10 S. Azadi, C. Cavazzoni, and S. Sorella, 
  Phys. Rev. B 82, 125112 (2010).
  Umrigar88 C. J. Umrigar, K. G. Wilson, and J. W. Wilkins, 
  Phys. Rev. Lett. 60, 1719 (1988).
  Neil05 N. D. Drummond, and R. J. Needs, 
  Phys. Rev. B 72, 085124 (2005). 
  BF P. López Ríos, N. D. Drummond, M.D. Towler, and R. J. Needs, 
  Phys. Rev. B 74, 066701 (2006)
  Feller D. Feller, J. Chem. Phys. 96, 6104 (1992)
  HelgakerT. Helgaker, W. Klopper, H. Koch, and J. Noga, 
  J. Chem. Phys. 106, 9639 (1997)
  Tajti A. Tajti, P.G. Szalay, A.G. Császár, M. Kállay, J. Gauss, E.F. Valeev, B.A. Flowers, J. Vázquez, and J. Stanton,
  J. Chem. Phys. 121, 11599 (2004)
  GTH S. Goedecker, M. Teter, and J. Hutter, 
  Phys. Rev. B 54, 1703 (1996)
  Kwon Y. Kwon, D. M. Ceperley, and R. M. Martin
  Phys. Rev. B 58, 6800, (1998)
  Neil08 N.D.  Drummond,  R.J.  Needs,  A. Sorouri, and  W.M.  C.Foulkes,
  Phys.Rev.B 78, 125106 (2008)
  FS15 S. Azadi, and W.M.C. Foulkes,
  J. Chem. Phys. 143 102807 (2015)
  Holzman16 M. Holzmann, R. C. Clay, M. A. Morales, N. M. Tubman, D. M. Ceperley, and C. Pierleoni,
  Phys. Rev. B 94, 035126 (2016)
  FS19 S. Azadi, and W.M.C. Foulkes, 
  Phys. Rev. B 100, 245142 (2019)
  Ceperley-FS D. M. Ceperley,  and B. J. Alder, 
  Phys. Rev. B 36, 2092 (1987)
  
  