Toward a systematic improvement of the fixed-node approximation in diffusion Monte Carlo for solids—A case study in diamond

While Diffusion Monte Carlo (DMC) is in principle an exact stochastic method for ab initio electronic structure calculations, in practice, the fermionic sign problem necessitates the use of the fixed-node approximation and trial wavefunctions with approximate nodes (or zeros). This approximation introduces a variational error in the energy that potentially can be tested and systematically improved. Here, we present a computational method that produces trial wavefunctions with systematically improvable nodes for DMC calculations of periodic solids. These trial wavefunctions are efficiently generated with the configuration interaction using a perturbative selection made iteratively (CIPSI) method. A simple protocol in which both exact and approximate results for finite supercells are used to extrapolate to the thermodynamic limit is introduced. This approach is illustrated in the case of the carbon diamond using Slater–Jastrow trial wavefunctions including up to one million Slater determinants. Fixed-node DMC energies obtained with such large expansions are much improved, and the fixed-node error is found to decrease monotonically and smoothly as a function of the number of determinants in the trial wavefunction, a property opening the way to a better control of this error. The cohesive energy extrapolated to the thermodynamic limit is in close agreement with the estimated experimental value. Interestingly, this is also the case at the single-determinant level, thus, indicating a very good error cancellation in carbon diamond between the bulk and atomic total fixed-node energies when using single-determinant nodes. © 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0021036., s

agreement with experiment. 9,10 In sharp contrast, for strongly correlated systems (e.g., transition metal oxides) for which electrons can no longer be treated as weakly interacting quasiparticles, the approaches mentioned above may dramatically fail, even qualitatively. 11,12 An alternative path to describe correlation effects in solids is to resort to wave-function-based correlated methods, as done in a number of recent works {e.g., second order Møller-Plesset perturbation theory (MP2), 13 coupled cluster single and double (CCSD) with perturbative triple [CCSD(T)], and the equation of motion-CCSD (EOM-CCSD) 14,15 }. However, as single-reference approaches, they are of limited use for strongly correlated materials. Designing efficient and general computational approaches with controlled and testable approximations is, therefore, of great interest.
Quantum Monte Carlo (QMC) methods have been developed to treat weakly through strongly correlated systems while invoking few approximations. These include variational Monte Carlo (VMC) and diffusion Monte Carlo (DMC) in the continuous space, 16,17 and full configuration interaction QMC (FCIQMC) 18,19 and auxiliary field QMC (AFQMC) 20,21 in the determinant space. These methods utilize stochastic integration to help treat the complexity of the many-body electronic structure problem and are applicable to both solids and isolated molecules. Importantly, it is becoming possible to compare these and other many-body methods with each other on equivalent problems and to both validate their results and focus improvements. 22,23 In this study, we will focus on DMC, which is one of the most widely used QMC approaches for calculating the electronic structure of solids. 16,24 DMC enjoys a relatively low computational scaling with system size and is well adapted to massive parallelism. [25][26][27][28][29] For bosons, DMC can deliver the exact energy to an arbitrarily small error. However, for fermionic problems, because of the uncontrolled fluctuations of the wavefunction sign, a DMC implementation that does not address the sign problem is not stable in time. 30,31 To address the fermion sign problem and correctly obtain the required fermionic symmetry, practical DMC calculations must rely on the fixed-node approximation: the fermionic constraints are fulfilled by imposing the zeros (nodes) of the wavefunction to be those of an approximate antisymmetric (fermionic) trial function. However, because these nodes are approximate, the DMC energy has a residual variational error, referred to as the fixednode error. Devising an efficient strategy to improve the trial wavefunction nodes in a systematic way, thus better controlling the only source of error of the method, brings us closer to achieving a genuine general from-first-principles theory, which is one of the major challenges of DMC.
One advantage of real space QMC methods is their ability to flexibly use different forms of trial wavefunctions. Since only the values and derivatives of the wavefunction are required, there is no need, in principle, to choose forms that are suited for conventional numerical integration. While it is now standard to utilize sophisticated multideterminant trial wavefunctions in molecules, most solid state applications of DMC use a single Slater determinant (SD) trial wavefunction. A Jastrow factor is used to incorporate dynamical correlations and improve the wavefunction overall, but it does not change the nodes. Although SD-DMC calculations have proven valuable in addressing a wide range of problems in chemistry and materials science, for example, obtaining near chemical accuracy in van der Waals binding, 32,33 this accuracy is far from general. A potentially general strategy to address this problem, which has been successfully applied to atoms and molecules, is to employ multideterminant (MD) trial functions. [34][35][36][37] However, this requires an efficient method to generate the multideterminant expansion.
In this work, we present a scheme to perform multideterminant DMC calculations in solids using configuration interaction (CI) based trial functions. We employ a selected CI (SCI) approach in which only the most important Slater determinants from the full CI (FCI) space are identified and selected. Although SCI does not suppress the exponential wall of conventional CI calculations, it may be deferred sufficiently to allow practical calculations. In the case of atomic and molecular systems, several successful applications of SCI trial functions in DMC calculations have been reported. 34,[37][38][39][40][41][42][43][44] Here, we use the configuration interaction using a perturbative selection made iteratively (CIPSI) algorithm, 45 as implemented in the QUANTUM PACKAGE code 46 to generate the trial functions for DMC calculations on periodic systems. This approach is illustrated by computing the cohesive energy of carbon diamond in the thermodynamic limit. We use the QMCPACK quantum Monte Carlo package 28,29 for the DMC calculations. Although diamond is not a particularly challenging system with respect to the role of electron correlation, it is a valuable test system since its cohesive energy is accurately known and pseudopotential accuracy is not a major concern. This work may, thus, be considered as a first step toward the study of more difficult solids where greater correlation effects and the generation of CIPSI expansions may be more challenging.
The paper is organized as follows: In Sec. II, the theoretical aspects of the approach are presented. A brief summary of the CIPSI method is given, and its extension to the case of periodic systems is described in depth. The computational details are provided in Sec. III. Results for diamond, both at the single-and multideterminant levels, are reported and discussed in Sec. IV. Finally, some concluding remarks are given in Sec. V, including areas where further improvements are desired. Unless otherwise stated, atomic units are used throughout this paper.

A. The CIPSI algorithm
In the present manuscript, the CIPSI algorithm developed for finite-size atomic or molecular systems is generalized to the case of periodic solids. CIPSI is one of the variants of the broad class of methods known as SCI. Selecting determinants in the CI space is a natural idea, and many SCI variants have been developed under various acronyms and implementations over the last five decades. 34,37,[39][40][41]45, SCI methods are ordinary CI approaches except that determinants are not chosen based on predetermined occupation or excitation criteria but are instead selected step by step based on their estimated contribution to the FCI energy and/or wavefunction. It is a particularly successful approach since it is well recognized that, in a predefined subspace of determinants, for example, single and double excitations with respect to the Hartree-Fock (HF) reference, only a small fraction of them make a non-negligible contribution to the wavefunction. 59,81 Therefore, an on-the-fly selection of determinants has been proposed in the late 1960s by Bender and Davidson, 47 as well as by Whitten and Hackmeyer. 48 The main advantage of SCI methods is that no a priori The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp assumption is made on the type of determinants needed to describe the electronic correlation effects. Therefore, at the potentially greater price of a blind and automated calculation, a SCI calculation is less biased by the user's understanding of the problem's complexity (for example, when choosing the active space (AS) orbitals for a particular problem). The CIPSI approach was developed in 1973 by Huron, Malrieu, and Rancurel. 45 In recent years, two of us have revived this approach 34 and developed a very efficient and massively parallel version of the algorithm that has been implemented in QUANTUM PACKAGE. 46 Briefly, at each iteration n, CIPSI selects some external determinants |Dα⟩ not present in the current zeroth-order wavefunction, where |DI⟩ is an internal determinant. Starting, for example, with the HF determinant at n = 0, the external determinants are selected using a criterion based on the estimated gain in correlation energy evaluated using the second-order perturbation theory that would result from the inclusion of |Dα⟩. Denoting the second-order correction as where is the variational energy of the wavefunction at the nth iteration, a number of external determinants associated with the greatest e var . In practice, the size of the variational wavefunction is roughly doubled at each iteration. Next, the second-order Epstein-Nesbet energy correction to the variational energy (denoted as E and the total CIPSI energy is given by The algorithm is then iterated until some convergence criterion (for example, |E PT2 | ≤ ϵ) is met. For simplicity, in the following, the superscript denoting the largest iteration number attained will be dropped from various expressions. When the number of determinants in the variational space gets large, computing the second-order correction E PT2 by adding up all contributions eα to the sum becomes computationally unfeasible. To perform this formidable task, a simple and efficient hybrid stochastic-deterministic algorithm has been developed. 62 In short, the leading contributions to the sum in Eq. (4) are exactly computed (deterministic part), while the very large number of small residual contributions are sampled via a Monte Carlo algorithm (stochastic part).

B. CIPSI for periodic systems
We now present the extension of CIPSI to periodic solids. As in any CI calculation, we must define (i) the system (number of electrons, charge, and positions of the nuclei), (ii) the Hamiltonian, and (iii) the one-electron basis set. The calculation of the one-and twoelectron integrals required for the computation of the Hamiltonian matrix elements then needs to be discussed.

Supercells
In practice, an infinite solid is modeled by a finite supercell obtained by replicating a given primitive cell ni times in each Cartesian direction (i.e., i = x, y, z). Accordingly, we will label a supercell as nx × ny × nz. In the present study, we will restrict ourselves to the simplest case of a cubic primitive cell of side L with arbitrary nx, ny, and nz values. The supercell is then a rectangular cuboid of volume Ω = N × L 3 , where N = nxnynz is the total number of primitive cells replicated to build the supercell.
The supercell being defined, we are led back to an ordinary CI calculation consisting of the set of electrons and nuclei present in the supercell and subject to an external periodic electric potential. We denote RI and ZI the position and charge of the of Ith nuclei of the primitive cell (I = 1, . . ., Nn). We emphasize that, in contrast to effective one-body theories such as HF or DFT, many-body electronic structure calculations explicitly taking into account the electron-electron interaction, such as the ones performed here, cannot be restricted to the primitive cell if accurate properties are to be obtained. 16 In one-body theories, the thermodynamic limit can be reached by indefinitely improving the k-point sampling of the first Brillouin zone of the primitive cell.
In the presence of electron-electron interaction, the translation of individual electrons is no longer a symmetry of the Hamiltonian (k is no longer a good quantum number) and the use of supercells is mandatory (see, e.g., Refs. 16 and 82 for a discussion of this aspect).

Hamiltonian
The supercell electronic Hamiltonian HN =T +Vne +Vee +Vnn (6) has a standard form, except that the electron-electron, electronnucleus, and nucleus-nucleus Coulombic potentials,Vee,Vne, and Vnn, respectively, are now periodized to model the interaction of the where ri is the position of the ith electron, and the prime symbol indicates that the self-interaction term i = j or I = J has to be excluded when T = 0 = (0, 0, 0).
As is well known, periodic Coulomb potentials are mathematically ill-defined due to the long-range character of the Coulomb interaction. The periodic infinite sums in Eq. (7) not only converge very slowly, but they are also conditionally convergent, meaning that the result depends on the order of summation. A specific and careful mathematical treatment has to be introduced in order to provide a meaningful answer. The standard solution, and the one we employ here for the three Coulomb potentials in Eq. (7), is to resort to the Ewald summation technique.
Applying the Ewald summation technique, for example, to the nucleus-nucleus termVnn in Eq. (7) enables one to compute this term as a sum of two contributions: a short-range contribution in the real space and a short-range in the reciprocal space. Both contributions are expressed as rapidly converging infinite sums, thus, leading to a very fast and efficient calculation of the potential. Explicitly, it is given byV where erfc(x) is the complementary error function, G refers to the set of quantized reciprocal vectors of the supercell defined by the condition e iG⋅T = 1, and σ is a small positive parameter chosen to facilitate convergence. Additional details can be found, for example, in Refs. 83 and 84.

Basis functions
In this work, the one-electron basis functions are chosen to be crystalline Gaussian-based atomic orbitals, i.e., the periodized (or translationally symmetry-adapted) version of the (localized) Gaussian atomic orbitalsχμ(r) from the supercell. In Eq. (9), the crystal momentum vector k is chosen within the first Brillouin zone of the primitive cell, and it is sampled from a Monkhorst-Pack grid, 85 an evenly spaced rectangular grid in the reciprocal space (see Sec. II C). The atomic index μ is referred to as the band index after periodization.
The molecular orbitals of the system are then defined as where N bas is the number of basis functions, and the molecular orbital coefficients Cμp(k) are now momentum-dependent due to the translational-symmetry adaptation of the basis functions. Once again, we emphasize that, because of the explicit treatment of the electron-electron interaction, k is no longer a good quantum number and the orbitals defined above do not have the correct translational symmetry of the problem. However, choosing such a representation is particularly convenient in practice since it allows us to take full advantage of the techniques and codes developed for the effective one-particle theories, such as HF and DFT.

One-and two-electron integrals
The Hamiltonian and the one-electron basis set having been defined, the next step in a CI calculation is to evaluate the Hamiltonian matrix elements between Slater determinants. This requires the evaluation of one-and two-electron integrals. It is easy to see from the expressions of the Hamiltonian [see Eqs. (7) and (8)] and basis functions [see Eq. (9)] that this can be accomplished by the calculation of the integrals over Fourier transforms of the product of Gaussian functions for which fast and reliable algorithms have been developed. 14,86,87 The only new aspect with respect to standard CI implementations is the use of complex-valued orbitals, integrals, and Hamiltonian matrix elements. Let us denote the two-electron integrals as Since the one-particle basis functions are invariant with respect to the primitive lattice translation vectors t, the two-electron integrals must conserve crystal momentum, i.e., kp + kq = kr + ks + g, where g is a reciprocal lattice vector of the primitive cell. When real-valued orbitals are used, a given two-electron integral is symmetric with respect to permutations of the orbital indices p and r, or q and s, as well as with the permutation of the electron labels 1 and 2, thus, resulting in a 8-fold symmetry of the set of twoelectron integrals. Consequently, the N 4 bas four-index integrals can be mapped to a set of approximately N 4 bas /8 unique integral values. 88 However, when one considers complex-valued orbitals, the integral ⟨ pkpqkq|rkrsks⟩ is no longer invariant with the permutations of p and r, or q and s, so the number of unique integral values scales as N 4 bas /2 (exchanging the electron indices 1 and 2 still leaves the value of the integral unchanged). When storing these complex-valued integrals, it is useful to recognize that exchanging p with r, or q with s changes the integral value in a non-trivial way. However, exchanging both of these pairs simultaneously yields the complex conjugate of the original value, i.e., ⟨pkpqkq|rkrsks⟩ * = ⟨rkrsks|pkpqkq⟩. This allows one to save an additional factor of two in the storage requirements; one then only needs to store N 4 bas /4 complex-valued two-electron integrals.

Hamiltonian matrix elements
In order to implement FCI or any of its lower-cost alternatives, one must be able to evaluate matrix elements of the Hamiltonian in the space of Slater determinants. If the determinants are all built from the same set of orthonormal spin-orbitals (which is the case here), then one can evaluate matrix elements via the usual Slater-Condon rules (see, e.g., Refs. 88 and 89 for an efficient implementation in a determinant-driven context).

C. Boundary conditions and twist averaging
In order to reduce finite-size effects and to accelerate the convergence to the thermodynamic limit (i.e., N → ∞), it is advantageous to exploit the freedom in the type of periodic boundary conditions of the supercell. By judiciously choosing the electron momenta k of the basis functions [see Eq. (9)], the boundary conditions can be easily implemented. Translating one electron of the supercell by a superlattice vector (say, A = La) generates a phase factor, e iLk⋅a , for each of the orbitals of the CI determinants. Accordingly, a global phase factor common to all determinants is obtained whenever all these individual phase factors are made equal, that is, where θ is some arbitrary angle (or twist) between −π and π. These conditions are fulfilled when momenta are chosen uniformly spaced in the first Brillouin zone of the lattice (that is, by using a Monkhorst-Pack grid 85 ) and shifted by a common vector K, with i l = 0, . . ., t l (l = a, b, c) being the twist index. In this case, we have |D⟩ → e iθ |D⟩ with θ = K ⋅ a. Boundary conditions can be varied with θ from −π to +π, θ = 0, and θ = ±π corresponding to periodic and anti-periodic boundary conditions, respectively. For a given system size, a property O can be computed by averaging out its values for different K values (or twists θ). In the limit of an infinite set of sampling values, we have where Ψ θ is the exact wavefunction of the system for the corresponding boundary condition. In practice, the integral is computed as a finite sum over shifted Monkhorst-Pack grids, as expressed by Eqs. (13a)-(13c). The main effect of twist-averaging is to suppress the major part of the one-body shell effects in the filling of singleparticle states. Note that each value of θ requires an independent calculation, the total computational cost is then proportional to the number of twists.

III. COMPUTATIONAL DETAILS
The DMC method is employed to compute ground-state energies within the fixed-node approximation. The application of DMC to solids being now well documented, we refer the interested reader to the existing literature for the theoretical background (see, for example, Refs. 16, 28, and 82). Calculations were performed using QMCPACK. 28,29 The integrals over periodized Gaussian-type orbitals (GTOs) are carried out using the PySCF program, 90 which has been interfaced with QUANTUM PACKAGE. 46 Plane wave trial wavefunctions were generated with the QUANTUM ESPRESSO package. 91 To model carbon diamond, we used a primitive cell with two carbon atoms separated by 1.545 Å (∼3.568 Å for the lattice constant). In order to compute the cohesive energy in the thermodynamic limit and to correct energies for finite-size effects, we consider supercells made of 16 (2 × 2 × 2), 54 (3 × 3 × 3), and 128 (4 × 4 × 4) atoms with converged twist-averaged (TA) boundary conditions with a total of 216, 64, and 27 twists, respectively. 92 All calculations are carried out using the Burkatzki-Filippi-Dolg (BFD) 93,94 effective core potentials and the associated 2s2p1d or 3s3p2d1f contracted Gaussian-type orbital (GTO) basis sets (BFD-vDZ or BFD-vTZ, respectively). For the plane wave trial wavefunction, consistently with Ref. 95, a 200Ry cutoff was applied to the kinetic energy. Both SD-DMC and MD-DMC calculations used a time step of 0.001 a.u. and Casula's T-moves for pseudopotential evaluation. 96 The trial wavefunctions for DMC consist of either a singledeterminant or multideterminant expansion determined using CIPSI multiplied by a Jastrow factor with one-, two-, and three-body terms. The parameters for the one-and two-body terms were represented by B-splines and the three-body term by a polynomial. These parameters were optimized in VMC through a variant of the linear method of Umrigar et al. 97 Note that the coefficients of the CIPSI expansion were not optimized in VMC, and therefore, in this study, the nodal surface is wholly determined by the CIPSI expansion.

A. Single-determinant fixed-node DMC
To establish a reference for the multideterminant studies, we first investigate the dependence of single-determinant DMC energies on the initial orbitals and their basis sets. In Fig. 1, we show the dependence of the fixed-node DMC energy for single-determinant trial wavefunctions built with local-density approximation (LDA), 98 Perdew-Burke-Ernzerhof (PBE), 99 Becke, 3-parameter, Lee-Yang-Parr (B3LYP), [100][101][102][103] and HF orbitals calculated with the BFD-vDZ and BFD-vTZ basis sets. These calculations were performed for the simplest case of a single primitive cell consisting of two carbon atoms (1 × 1 × 1 supercell), with the SD-DMC energy being computed at the Γ point. As an additional reference, the SD-DMC energy obtained using LDA orbitals expanded in a large plane-wave (PW) basis set with a high energy cutoff of 175 hartree and the same pseudopotential is also reported. Using a larger energy cutoff of 250 hartree leads to a difference in DMC energies below 1 mhartree, and therefore, our initial choice of 175 hartree is close to the complete basis set limit. Note that the fixed node SD-DMC of Chemical Physics energies using nodal surfaces of the Kohn-Sham determinants obtained from different DFT approximations differ by less than the variation in the pure DFT energies, 6.8(6) mhartree vs 13.05 mhartree, respectively. The SD-DMC energies obtained with LDA orbitals expanded in the BFD-vTZ and PW basis sets coincide within 1.0(6) mhartree, thus, indicating the suitability of the BFD-vTZ basis set. On the other hand, an appreciably higher SD-DMC energy is obtained when using the BFD-vDZ basis set. Consequently, in the following, all calculations are performed with the BFD-vTZ basis set.
In Table I, we report twist-averaged SD-DMC total energies (per primitive cell) as a function of the supercell size N, as well as the corresponding cohesive energies. We also report the SD-DMC energies computed by Shin et al. in Ref. 95, which were obtained using a TABLE I. Twist-averaged SD-DMC total energies (in hartree/cell) and cohesive energies (in hartree) of diamond for supercells of increasing size N = n 3 (n = 2, 3, and 4). The (extrapolated) thermodynamic values and the cohesive energies (with and without ZPE correction) are also reported. PW basis set and a single Slater determinant wavefunction of PBE orbitals and using the same twist values. 104 Results are compared to our SD-DMC data using B3LYP orbitals and the BFD-vTZ basis set. The SD-DMC total energies from the two calculations are very close, with a difference never greater than 3 mhartree. Electronic cohesive energies calculated by subtracting from the crystal energy the SD-DMC energy of the isolated atoms calculated using the same approaches [ −5.4226(2) hartree/atom for PBE/PW and −5.421 38(2) hartree/atom for B3LYP/GTO] lead to a difference of only 3.0 ± 0.4 mhartree once extrapolated to the thermodynamic limit. Adding the zero-point vibrational energy (ZPE) correction estimated to be 6.0 mhartree/atom 105 brings the SD-DMC(B3LYP/GTO) cohesive energy to 1.3 ± 0.4 mhartree of the estimated experimental value 106 of 0.2699 hartree (see Table I). The result indicates that there is very good error cancellation between the bulk and atomic total energies at the single-determinant level in carbon diamond.

B. Multideterminant trial wavefunctions
In Fig. 2, we show the convergence of the CIPSI variational energy Evar and its second-order corrected value Evar + E PT2 as a function of the number of determinants in the reference space (see Sec. II A) for the 1 × 1 × 1 diamond primitive cell at the Γ point and using the BFD-vTZ basis set. Independent calculations using both HF and B3LYP orbitals are shown. Despite the large size of the Hilbert space (about ∼10 11 determinants), energy convergence is achieved with a reasonable accuracy using less than ∼10 7 determinants. In addition, the FCI limit is independent of the orbital set employed in the SCI calculation, as it should be. In the variational calculations, convergence to millihartree accuracy is achieved at about 3 × 10 6 determinants using whether HF or B3LYP orbitals.
When the perturbative corrections to the energy are included, millihartree accuracy is reached with only ∼2 × 10 3 and ∼4 × 10 3 determinants for HF and B3LYP orbitals, respectively.
Although CIPSI rapidly converges to the FCI limit for the 1 × 1 × 1 supercell (Γ point calculations), the exponential character of the FCI approach is still present, and for larger supercells, the number of determinants required to reach the FCI limit to the required accuracy rapidly becomes computationally intractable. To alleviate this issue and to enable treating larger supercells, we limit the SCI calculations to a subset of orbitals belonging to a pre-defined active space, thus, performing, in practice, a CIPSI calculation for a set of Ne electrons distributed among No orbitals, denoted as CIPSI (Ne, No), instead of a FCI one. Of course, there are many ways of choosing such an active space (AS). A natural choice is to build the AS using the natural orbitals obtained in a preliminary CIPSI calculation. This is what has been systematically done for molecular systems in a number of works published by some of the authors. [39][40][41]43,62,67,68,[78][79][80] Here, in view of the large number of orbitals needed for solids, an alternative choice could be to employ instead the natural orbitals of a preliminary MP2 calculation. Another possibility is to use Kohn-Sham orbitals instead of Hartree-Fock ones. Here, we have chosen to postpone the investigation of such important aspects to a future detailed study, and we have considered the simplest approach where only virtual orbitals with one-electron energies below a given energy threshold are included in the active space. Although this choice is rather crude (and certainly not optimal), it is important to emphasize that the choice of the active space for constructing the multideterminant trial wavefunction is expected to be less critical in DMC than when just limiting the calculation to pure CI. High-energy orbitals are, indeed, expected to have a small impact on the wavefunction and, thus, on the location of the nodes. Figure 3 illustrates the use CIPSI(Ne, No). The error in DMC total energy (in hartree/atom) computed at the Γ point as a function of the number of orbitals belonging to the active space is presented. The error is calculated with respect to the converged DMC value obtained with 58 orbitals. For the 1 × 1 × 1 primitive cell, results are reported for both B3LYP and HF orbitals and for a number of orbitals in CIPSI (8,No) up to the maximum of 58. For the 2 × 1 × 1 supercell, only the B3LYP results are presented, and the maximum number of orbitals in CIPSI (16,No) corresponds to only half of the total number of orbitals (116). For the 1 × 1 × 1 supercell, a wellconverged SD-DMC energy is achieved when using only 30 of the 58 B3LYP orbitals, but all orbitals must be retained to obtain a wellconverged result when using HF orbitals. For the larger 2 × 1 × 1 cell for which the CIPSI (16,No) is built with B3LYP orbitals, convergence to within 1.0 ± 0.3 mhartree/atom is reached between 30 and 40 orbitals out of a total of 116. Due to the better convergence of the calculations with truncated orbital spaces when using B3LYP than HF orbitals, B3LYP orbitals are used in all subsequent calculations.  Table II. Here, ϵ is a threshold associated with the number of determinants retained in the expansion, and ϵ = 0 corresponds to the full CIPSI wavefunction (which is much smaller than the FCI space due to a limit on the total number of determinants retained, see Ref. 107). In view of the computational cost of using a very large number of determinants in the trial wavefunction in DMC calculations, it is important to limit the determinants used to the smallest number possible. As seen in Table II, the fixed-node energies are essentially converged within statistical errors at ϵ = 10 −5 and, thus, considering the full CIPSI wavefunction (ϵ = 0) in DMC is not necessary here. Note that the number of determinants needed to reach a given threshold ϵ is smaller for the FCI space than the CIPSI (8,30) space. This results from the greater flexibility available with the FCI space. The convergence curves of the MD-DMC energy both for the FCI and CIPSI (8,30) spaces are similar, and at near-convergence, the two DMC energies differ by only about 2 ± 1 mhartree. The DMC intrinsic variance in the energy shows a systematic improvement with the number of determinants in the expansion when using the FCI space. As seen from Fig. 4, the DMC energy decreases monotonically and smoothly when increasing the number of determinants in the trial wavefunction. As discussed in Ref. 41, there is no guarantee that increasing the number of Slater determinants in the trial wavefunction lowers the DMC energy because the selection procedure does not explicitly optimize the nodal surface and DMC energy. However, in all applications performed so far-atoms, 34,42,107 molecules, [38][39][40][41]61 and now solids-a systematic decrease of the fixed-node DMC energy is observed whenever the SCI trial wavefunction is improved variationally, upon increasing the number of determinants, the size of the basis set, or the size of the active space.

C. Multideterminant fixed-node DMC
The DMC total energy (in hartree/cell) of diamond as a function of 1/N, where N is the size of the system, i.e., the total number of primitive cells replicated to create the supercell (N = 8, 18, 27, and 64 for 2 × 2 × 2, 3 × 3 × 2, 3 × 3 × 3, and 4 × 4 × 4, respectively), computed with various methods and approximations is reported in Fig. 5. The two upper curves (solid lines) report SD-DMC (red) and MD-DMC (black) energies computed at the Γ point. The two lower curves (dashed lines) are SD-DMC (red) and MD-DMC (black) energies obtained by twist-averaging, as described in Sec. II C. To minimize the statistical fluctuations, we re-optimized the Jastrow factor at each twist. The number of determinants in the multideterminant trial wavefunction varies from about 600 000 to 1 200 000.
SD-DMC calculations (at the Γ point and twist-averaged) have been performed for the 2 × 2 × 2, 3 × 3 × 3, and 4 × 4 × 4 supercells. MD-DMC calculations are much more computationally demanding than their SD-DMC analogs. Consequently, MD-DMC calculations have been limited to the 2 × 2 × 2 and 3 × 3 × 3 supercells, with the energy of the 4 × 4 × 4 supercell being estimated as explained below. In all cases, we only computed the symmetry inequivalent twists in the Brillouin zone, using symmetries to reduce the number of expensive DMC calculations. There are 16, 8, and 4 inequivalent  To further reduce the computational effort, an approximate version of the twist-averaged MD-DMC energies was employed for the 3 × 3 × 3 and 4 × 4 × 4 superlattices. In this procedure, we estimate the energy at a given twist (other than Γ) via where wi denotes the ith twist (w 0 being the Γ point) and E B3LYP HF is the energy computed with the HF energy functional and B3LYP orbitals. This approximation relies on the B3LYP band structure being sufficiently accurate but can be further improved by considering how these data are extrapolated to the thermodynamic limit, as we will detail below. Figure 6 reports the SD-and MD-DMC energies computed both exactly and using Eq. (15) for the 16 independent twists of the 2 × 2 × 2 supercell. To facilitate the visualization of the actual differences between data, the 16 DMC energies are artificially connected with straight lines. The two upper (red) curves and the two lower (black) curves correspond to SD-DMC and MD-DMC calculations, respectively. As one can see, there is an excellent agreement between the exact (solid lines) and approximate (dashed lines) treatment of twist averaging in the case of SD-DMC, introducing an error of only 1.7(2) mHa in the twist averaged energy. For MD-DMC, the agreement is not as good but remains satisfactory for the present purpose, introducing an error of 8(4) mHa to the twist averaged energy.
We can estimate the twist-averaged DMC energy of the 4 × 4 × 4 supercells by combining the extrapolated value of the DMC energy at the Γ point and a twist-averaged contribution from Eq. (15). The extrapolated value reported in Fig. 5 has been obtained as the value resulting from a quadratic fit of the three values corresponding to the Γ point energies of the 2 × 2 × 2, 3 × 3 × 2, and 3 × 3 × 3 supercells. Note that, as is common in solid state calculations, a more precise fit function based on the theoretical basis could be employed (see, for example, Refs. 82 and 108). However, here, in view of the remaining errors that are certainly much larger than the difference in error between various extrapolation fit functions, we have used the simplest possible fit function that reproduces the approximate parabolic shape of the curves. Fitting our data only within the domain of size where a quasi-linear regime is reached would clearly be preferable but the high computational cost of the MD-DMC data restricts the range of accessible sizes.
Finally, in the spirit of extracting the maximum of information from our limited set of data, we propose to exploit the fact that both Γ point and twist-averaged energies must converge to the same value in the thermodynamic limit (i.e., N → ∞). This exact property is used as a constraint when fitting simultaneously both sets of energies with quadratic expressions, i.e., The five parameters in Eqs. (16a) and (16b), i.e., E N→∞ DMC and the four ci's, are obtained by minimizing the χ 2 -type function, where δE DMC 's are the corresponding statistical errors, and the sum runs over Ni = 8, 18, 27, and 64 for the Γ point energies (centered grid) and Ni = 8, 27, and 64 for the twist-averaged (TA) energies.
The quantity E N→∞ DMC represents the best estimate of the DMC energy in the thermodynamic limit.
Following this extrapolation procedure, the SD-DMC total energy in the thermodynamic limit is found to be E N→∞ SD-DMC = −11.3986(2) hartree. Combined with the atomic SD-DMC total energy of −5.4214(2) hartree, it yields a cohesive energy (including the ZPE contribution) of 0.2719(3) hartree, an estimate close to the value of 0.2712(4) hartree obtained by a simple quadratic fit of the DMC energies computed at the Γ point (see Table I). In the multideterminant case, we find E N→∞ MD-DMC = −11.4246(4) hartree, and using the atomic MD-DMC energy computed at −5.4335(1) hartree, we obtain a ZPE-corrected cohesive energy of 0.2729 (1) hartree. This value compares quite well with the estimated experimental cohesive energy of 0.2699 hartree extracted from Ref. 106. While this value is also similar to that obtained with singledeterminant wavefunctions, there has been a significant lowering in the supercell and atomic energies due to use of multideterminant trial wavefunctions. In particular, the energy of the atom obtained from the MD-DMC calculations is essentially exact (for the chosen pseudopotential).
It should be emphasized that comparing calculated and experimental cohesive energies is subject to the error made in estimating the ZPE contribution. The value of 6.0 mhartree/atom employed in this work was obtained by Schimka et al. using a zero-point anharmonic expansion based on DFT energy calculations. 105 19 In our calculations, because the cohesive energy of the solid is slightly too large and the pseudoatom is solved essentially exactly, the remaining difference with the experimental data must lie with a combination of the uncertainties in the extrapolation, time step error in the DMC calculations, or pseudopotential construction and evaluation. The ZPE error is also sizable compared to the residual difference from experiment in all of the literature predictions.
For the 2 × 2 × 2 supercell containing 16 carbon atoms, it is also instructive to compare our results to those of Ríos et al. 110 using the nodes of an optimized backflow (BF) trial wavefunction. The pseudopotentials used are different, making a direct comparison of the energies not possible, although they are both around −11.4 hartree per primitive cell (see Table X

V. CONCLUSIONS
We have demonstrated the feasibility of fixed-node DMC calculations for periodic solids using large multideterminant trial wavefunctions built from SCI expansions. Using as an example carbon diamond, this procedure is shown to be able to improve systematically the nodal surface of the trial wavefunction. In particular, we have found that the fixed-node DMC energy decreases monotonically and smoothly as a function of the number of determinants in the trial wavefunction.
Performing DMC calculations using large CI expansions together with large supercells is not feasible without further approximations due to overall cost. Here, this issue was addressed by introducing an approximate, yet well-defined, protocol combining both exact and approximate results for finite supercells and a controlled fitting procedure to reach the thermodynamic limit. Although approximate, our estimate of the diamond cohesive energy is, to the best of our knowledge, the first example of a fully ab initio MD-DMC calculation of a periodic solid. In our protocol, only the Jastrow parameters are optimized in the VMC step, with the linear coefficients of the CIPSI expansions being kept fixed. The main motivation for not optimizing the coefficients is to exploit the property of SCI methods to provide in a simple, well-defined (linear optimization), and systematic way of generating sequences of wavefunctions of increasing quality, leading to a systematic reduction of the fixed-node error (see Figs. 3 and 4). Defining such sequences is important, for example, when extrapolating fixed-node energies at different supercell sizes or computing an energy difference such as the cohesive energy. However, the computational cost of using (very) large CI expansions in DMC calculations is high, and strategies for generating more compact expansions are certainly desirable. Based on multideterminant DMC studies for isolated molecules (see, for example, Refs. 44, 111, and 112), a practical solution consists of optimizing the CI coefficients in the presence of the Jastrow factor and keeping only the most important determinants. However, such optimizations are challenging because the parameters in the Jastrow factors enter non-linearly and the objective function that one must minimize is evaluated stochastically. From a more general perspective, it would also be desirable to investigate the comparative performance with other types of trial wavefunctions including those with backflow as introduced by Ríos et al. 110 (a first numerical comparison is presented in this work) or of geminal forms. 113 Finally, we note that, to exploit the full potential of the present approach, more challenging materials must be investigated. In addition, it would be instructive to compute physical properties other than the cohesive energy that would potentially display a different dependence on the nodal error. We hope to be able to address these points in future studies.

SUPPLEMENTARY MATERIAL
See the supplementary material for the raw data associated with each figure of the manuscript. This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-publicaccess-plan).

DATA AVAILABILITY
The data that support the findings of this study are available within the article and its supplementary material.