Toward Large-Scale AFQMC Calculations: Large Time Step Auxiliary-Field Quantum Monte Carlo

We report modifications of the ph-AFQMC algorithm that allow the use of large time steps and reliable time step extrapolation. Our modified algorithm eliminates size-consistency errors present in the standard algorithm when large time steps are employed. We investigate various methods to approximate the exponential of the one-body operator within the AFQMC framework, distinctly demonstrating the superiority of Krylov methods over the conventional Taylor expansion. We assess various propagators within AFQMC and demonstrate that the Split-2 propagator is the optimal method, exhibiting the smallest time-step errors. For the HEAT set molecules, the time-step extrapolated energies deviate on average by only 0.19 kcal/mol from the accurate small time-step energies. For small water clusters, we obtain accurate complete basis-set binding energies using time-step extrapolation with a mean absolute error of 0.07 kcal/mol compared to CCSD(T). Using large time-step ph-AFQMC for the N2 dimer, we show that accurate bond lengths can be obtained while reducing CPU time by an order of magnitude.


■ INTRODUCTION
The accurate description of electron−electron correlation in large systems remains a fundamental challenge in solid-state physics and quantum chemistry.When beyond DFT accuracy is needed, one commonly resorts to coupled-cluster singles, doubles, and perturbative triples (CCSD(T)) 1,2 or fixed-node diffusion Monte Carlo (FN-DMC). 3,4−25 Extensive tests on small systems show that these methods yield nearly exact FCI solutions.As a result, they serve as valuable benchmark tools or high-level solvers.However, their application is currently limited to modest system sizes due to high computational costs and other inherent limitations.Besides the mentioned methods, phaseless auxiliary-field quantum Monte Carlo (ph-AFQMC) 26 is one of the most promising alternatives to both CCSD(T) and FN-DMC because of its favorable scaling with system size and because it is readily applicable to molecular and extended systems.
C ̌i ́zěk and Paldus 1,2 introduced the coupled-cluster theory 50 years ago, and since then, it has found extensive use in correlation-consistent calculations.The gold standard of quantum chemistry�coupled-cluster singles, doubles, and perturbative triples (CCSD(T)) 27−29 �has been carefully benchmarked on a wide range of small molecules, achieving nearly chemical accuracy. 29,30−35 However, CCSD(T) is often not accurate enough, especially in the presence of strong static correlation and spin contamination. 29For instance, the total energy of the CN radical in a double-ζ basis set converges to within 0.3 kcal/mol only at the CCSDTQP level of theory. 36ore importantly, the CCSD(T) theory already exhibits adverse computational scaling N ( ( )) 7 and higher coupledcluster expansions, such as CCSDTQP, are accessible only for the smallest systems (20 electrons in 100 orbitals).Similar to stochastic and selected versions of CI, analogous coupled-cluster variants have been proposed.−42 More extensive benchmark data for these methods is necessary to judge the performance and accuracy relative to CCSD(T).
Fixed-node diffusion Monte Carlo 3,4,43 (FN-DMC) is a commonly used method for studying solids and molecules containing up to 1000 electrons.It uses a random walk in real space to sample the Hamiltonian expectation values and scales as N N ( ) ( ) 3 4 with system size.DMC requires accurate trial wave functions to guide the walker density to the important positions in real space.Typically, Jastrow factors 44−46 on top of the Slater determinant impose cusp conditions.While Jastrow factors yield faster convergence to the complete basis set (CBS) limit, FN-DMC becomes comparable to other quantum chemistry methods only in the CBS limit.−49 Petruzielo et al. 49 demonstrated that a variational optimization of orbitals within the Slater-Jastrow wave function reduces the error to 2.1 kcal/mol.Furthermore, adding multiple Slater determinants from small complete active spaces reduces the mean absolute error to 1.2 kcal/mol for the G2 set.To mitigate the Fermionic sign problem, DMC relies on the fixednode approximation.This constrains the nodal structure of the exact many-body ground-state wave function to match the nodes of the trial wave function.However, the trial wave function is the primary source of errors in DMC. 3,50More accurate trial wave functions such as backflow wave functions, 51 geminal wave functions, 52 pfaffians, 53 or recently introduced deep learning models, 54,55 provide a route for further improvements of the method.
In this work, we will focus on the phaseless auxiliary-field quantum Monte Carlo (ph-AFQMC).Although originally developed as the path integral MC algorithm, 56−59 it experienced a renaissance with the introduction of openended random walks and the phaseless approximation. 26Since then, numerous applications for molecules 60−69 and solids 70−74 have confirmed that ph-AFQMC is a promising alternative to CCSD(T) and FN-DMC.−79 ph-AFQMC exhibits low polynomial scaling N N ( ( ) ( )) 3 4 with system size that can even be reduced to cubic scaling by utilizing tensor hypercontractions, 80,81 stochastic resolution of identity, 82 or employing a plane wave basis. 70The phaseless approximation controls the Fermionic sign problem in ph-AFQMC, however, it also introduces systematic errors that are challenging to control. 26,83Mahajan et al. developed a fast method for the evaluation of the local energies for multideterminantal trial wave functions and demonstrated that the phaseless approximation can be systematically improved using better trial wave functions. 84,85Recently, there have been attempts to directly enhance the phaseless approximation. 69,86,87n this work, we focus on reducing the computational prefactor of ph-AFQMC by utilizing large time steps.Larger time steps facilitate faster equilibration, diminish the autocorrelation between subsequent local energy evaluations, and thereby significantly reduce the computational cost to achieve fixed statistical errors (see Sec. Theoretical Background).Therefore, every projector QMC method employs the largest time step not compromising the accuracy.For typical systems, ph-AFQMC utilizes time steps around 5 × 10 −3 Eh −1 without significant time-step errors. 69,83,88−91 In FCIQMC, the time step is limited by the maximal excitation energy. 8ecifically, we modify the AFQMC algorithm to enable sizeconsistent simulations with large time steps.We demonstrate that the size-consistency errors in the ph-AFQMC method do not originate from the cosine projection, as previously suggested, 92 but are due to hybrid energy capping.To eliminate the remaining time-step errors, we investigate various approaches to compute the action of the one-body propagator on a Slater determinant.The optimal method minimizes the required number of H ̂|Ψ⟩ operations for a set of representative systems.Further, we adapt various propagators within the AFQMC framework and analyze their corresponding leading error terms.We show that the time-step errors agree with the theoretical predictions and that a reliable and robust time-step extrapolation to a zero time-step limit is possible.Finally, we show the accuracy of time-step extrapolated AFQMC energies for a variety of molecular systems and Gaussian basis sets.
The remainder of the paper is organized as follows: in the next section we briefly summarize the ph-AFQMC method, and we explore the advantages associated with the use of large time steps.We investigate different propagators within the ph-AFQMC framework in Sec.AFQMC Propagators, followed by a comparison of methods to compute the exponential of the onebody operator in Sec.Exponential of the One-Body Operator within AFQMC.Then, practical recipes for ph-AFQMC simulations with large time steps are presented in Sec.Modified AFQMC Algorithm.We validate our method through various applications in Sec.Applications and Discussion.Finally, we summarize our findings and outline the perspectives of AFQMC in the simulation of real materials in Conclusion.
■ THEORETICAL BACKGROUND AFQMC Review.In this section, we will briefly introduce the ph-AFQMC formalism.For a more detailed overview of the theory, we recommend the reviews in refs 26, 83, 93.Auxiliary-field quantum Monte Carlo projects an initial state |Ψ I ⟩ onto the exact ground-state wave function |Φ⟩ using the imaginary time Schrodinger equation where âp † and âq are Fermionic creation and annihilation operators, respectively and the h pq matrix elements represent the one-body Hamiltonian H ̂1. The two-body terms V pqrs are usually the electron repulsion integrals ⟨pq|rs⟩ and must be positive definite to decompose them to using a spectral decomposition, 26 an iterative Cholesky decomposition, 83,94−96 density fitting techniques, 97−99 or plane waves. 70,74,100The indices p, q, r, and s go over N basis functions, while the index g goes over the size of the decomposition N g .In the Gaussian basis, we usually have N g

Journal of Chemical Theory and Computation
≈ 10N.In this work, we use a Cholesky decomposition and refer to the L g,pq tensors as Cholesky vectors.Specifically, we opt for the iterative Cholesky procedure because a single threshold parameter controls the accuracy of the decomposition and the algorithm stops when reaching the desired accuracy, i.e., does not compute superfluous L g,pq tensors.
To treat the two-body Hamiltonian in eq 1, AFQMC employes the Hubbard-Stratonovich transformation 101,102 e x p e x d ( ) ( ) to map the interacting system onto a noninteracting one coupled to a fluctuating random field.p(x) denotes an N g -dimensional standard Gaussian probability density with the random numbers x g , while L ̂g is the one-particle operator formed from Cholesky vectors Equation 4 is the core of the method and the integration is performed using Monte Carlo sampling.
where the mean-field shift L̅ g and the force bias f g w minimize the variance of the sampling The with the hybrid energy i k j j j j j y and the importance sampling factor The choice of the propagator U ̂in eq 10 is thoroughly explored in Sec.AFQMC Propagators.Equations 10 and 11 define the free-projection AFQMC (fp-AFQMC).Although formally exact, fp-AFQMC suffers from the Fermionic phase problem, where walkers can acquire any phase e iθ with θ ∈ [0, 2π).It is equivalent to the sign problem in diffusion Monte Carlo except that for DMC walkers only e iθ = ± 1 are possible.As a consequence, the Monte Carlo estimator of the quantity ⟨Ψ T |Φ⟩ averages to zero.To circumvent this problem, Zhang et al. 26 proposed the phaseless approximation (ph-AFQMC), i.e. a modification of the eq 11 where Δθ is the phase change in the overlap that one walker acquires after one time step The phaseless approximation keeps the algorithm stable, but it also introduces a systematic error of order ( ) known as the phaseless error.It is widely assumed that the phaseless constraint contributes to significant size-consistency errors at larger time steps in ph-AFQMC. 92We investigated modifications to the phaseless approximation that significantly reduce the overcorrelation issues frequently encountered in ph-AFQMC, but they did not demonstrate systematic improvement compared to the standard phaseless approximation. 69However, the phaseless errors can be systematically reduced by improving the trial wave functions 84,85,104 or by release constraint techniques. 86,87e perform T Monte Carlo steps with the time step τ and T energy measurements where The time average E̅ and the variance σ E 2 are defined as where E̅ is normally distributed around a mean E with the variance The initial T eq samples define the equilibration phase and are excluded from the average to avoid excited state mixing.The variance in eq 20 is scaled by a correlation length κ, representing the number of steps after which two samples E t and E t+κ are statistically independent.

Journal of Chemical Theory and Computation
Computational Advantage of the Large Time-Step AFQMC.The first advantage of ph-AFQMC with large time steps is the reduction of the number of steps to equilibrate the systems.We need to propagate each walker for T eq steps so that the initial state |Ψ I ⟩ is projected to the ground state.The number of required equilibration steps T eq decreases like 1/τ with the size of the time step.
The second advantage of large time-steps is the reduction of the autocorrelation between subsequent samples.After taking a sample of the energy, one needs to propagate for κ steps to avoid autocorrelation between the samples.Ideally, we would increase the time step until κ = 1 so that we can sample the energy after every time step.
In practice, the benefits of large time steps are reduced by the fact that the energy does not have to be measured at each time step, if the samples are correlated.However, other tasks that constitute the random walk (creation of the effective Hamiltonian, AFQMC propagation, and the force bias calculation) are still needed at each time step. 69

■ AFQMC PROPAGATORS
The exact many-body propagator is represented by where H̅ 1 and H̅ 2 are the one-and two-body terms defined in eq 2. In the AFQMC method, it is formally necessary to employ a split-operator method to separate e −τH ̂2 from e −τ(H ̂1 + H ̂2) , for example using first order split technique (Split-1) For e −τH ̂2 the Hubbard-Stratonovich transformation is then used: with the effective interaction introduced in eq 8.After some algebra, we derive the leading error term of this transformation , where {A, B} = AB + BA and [A, B] = AB − BA.The exact propagator is reproduced in the limit N w → ∞ and in case that L̅ g operators commute.Note that in eq 24, we only show the systematic errors remaining after integrating over the random fields and not the statistical errors from the Monte Carlo sampling.
After applying the Hubbard-Stratonovich transformation, one can choose to either evaluate the split propagator directly or refactorize it back to derive a full effective Hamiltonian and then apply any approximation to e w .In the subsequent analysis, we examine two split propagators where the matrix exponentials are evaluated by applying methods outlined in the next section.Additionally, we examine two more propagators where the refactorized e w is approximated directly.Split-1 Propagator.The simplest propagator splitting of the one-body Hamiltonian and the effective interaction is This has the advantage of reducing the number of necessary operations, however, the split-1 propagator leads to a sizable time-step error with the leading term 28)   where E ̂HS is the error arising from the Hubbard-Stratonovich transformation (eq 24).Since, the one-body propagator e −τH ̂1 remains constant during the simulation, we precomputed it exactly using matrix diagonalization.Conversely, e w int is generated anew for each walker and time step.Sec.Exponential of the One-Body Operator within AFQMC discusses efficient techniques for the evaluation of this term.
Split-2 Propagator.Split propagators are more accurate with second-order decomposition The advantage of this approach is that it does not introduce additional errors beyond the Hubbard-Stratonovich transformation to quadratic order Similar to the Split-1 propagator, the term e −τH ̂1/2 is precomputed exactly in our implementation.
Taylor Propagator.The simplest propagator within the AFQMC formalism is obtained by the Taylor expansion of the refactorized time propagator The leading error term in this expansion is , , provided that the Taylor expansion is truncated at an order k ≥ 4.
Crank-Nicolson Propagator.The Crank-Nicolson propagator replaces the matrix exponential by a linear problem The corresponding leading error term is with Summary.In summary, the error terms specified above suggest that the Split-2 propagator is the optimal propagator within the AFQMC framework with the leading error term dictated by the Hubbard-Stratonovich transformation.We will show below that this is also confirmed by numerical examples,

Journal of Chemical Theory and Computation
albeit Taylor expansion will also yield excellent results almost on par with the Split-2 propagator.Conversely, the split-1 and the Crank-Nicolson propagator show significantly worse results.

■ EXPONENTIAL OF THE ONE-BODY OPERATOR WITHIN AFQMC
To propagate using the Split-2 propagator, we need to compute the action of the matrix exponential e w int and then apply it to the orbitals in |Ψ w ⟩.One could use exact matrix diagonalization 105,106 to compute e w int ; however, this is computationally expensive especially because w int is non-Hermitian (see eq 8).
While of limited practical utility, spectral decomposition still serves as a valuable benchmark tool.We will use it to evaluate the accuracy of iterative methods.
An alternative is to use iterative methods to directly compute e Taylor Expansion.The Taylor expansion offers a simple and robust scheme to approximate matrix exponentials k is the order of the approximation and the number of A|Ψ⟩ operations needed to approximate the exponential.The scaling of the algorithm is dominated by A|Ψ⟩ operations and scales as kN 2 N e , where k is the expansion order, N denotes the number of basis functions, and N e is the number of electrons in the system.Taylor expansion with k = 6 is usually used in the AFQMC community. 85,107hebyshev Expansion.Ezer et al. 108 demonstrated that, although the error terms of different polynomial expansions are of the same order, the Chebyshev polynomials T n (x) minimize the prefactor of the leading error term compared to other polynomial expansions.They are conveniently defined using the following recurrence relation The Chebyshev expansion is then formally given by For the matrix exponential, the expansion coefficients where J n (z) are the Bessel functions of the first kind.The k-th order Chebyshev expansion requires k A|Ψ⟩ operations and scales as kN 2 N e , same as the Taylor expansion.Krylov Subspace Projection.The Krylov subspace A ( , ) k for the matrix A and an orbital vector ψ is In our case, ψ is the N-dimensional vector of the orbital expansion coefficients in the Slater determinant |Ψ⟩ = |ψ 1 , ψ 2 , ..., ψ Nd e ⟩.For the non-Hermitian A, one can use the Arnoldi algorithm 105,106 to produce an orthonormal basis

., b k ] and the operator A k in the subspace
A k is a k × k upper Hessenberg matrix. 105B k is a N × k matrix with the first column b 1 = ψ/∥ψ∥ and Using this property, one can approximate the matrix exponential as follows: where e 1 is the first unit vector in k and B k e 1 = b 1 .The Krylov approximation to the matrix exponential of the order k is then simply given by Like the polynomial expansions, the A|Ψ⟩ applications dominate the computational cost of the Krylov method with kN 2 N e operations.In addition, the Arnoldi algorithm scales as kNN e 2 and one needs to compute the matrix exponential in subspace e −τA k .The latter is typically not expensive because k is very small compared to N; one computes it using matrix diagonalization, a higher-order polynomial expansion, or a Padeá pproximation.
Block-Krylov Subspace Projection.Contrary to the previous algorithm that expands each orbital in the Slater determinant individually, the block-Krylov projection simultaneously computes the action of the exponential operator on all orbitals.As a consequence, the Krylov subspace increases by N e elements in each Arnoldi iteration, and k-th order Krylov subspace consists of kN e vectors.The Krylov subspace is where Ψ denotes the N × N e dimensional matrix of the expansion coefficients.The orthonormal basis B k = [B̅ 1 , B̅ 2 , ..., B̅ k ] is now a N × kN e matrix, where Ψ = B̅ 1 R, and can be computed using a QR decomposition.Note that the initial orthogonalization of orbitals is necessary due to the imaginarytime propagation employed in AFQMC.The projected operator A k is computed analogously to eq 41 and is an upper block Hessenberg matrix of the size kN e × kN e .The approximation of the exponential is similar to eq 44 where E Nd e is the kN e × N e part of the kN e -dimensional identity matrix.For a given order k, the blocked Krylov method requires the same number of A|Ψ⟩ operations as the nonblocked version.The advantage of the blocked algorithm is that it may converge in fewer iterations because the space spanned by other orbitals is used to expand the exponential as well.Furthermore, matrix−matrix multiplications often achieve more floating-point operations per second than matrix-vector multiplications.The

Journal of Chemical Theory and Computation
disadvantage is that the dimension of the Krylov subspace scales with the number of electrons N e so that the Arnoldi iteration requires additional k 2 NN e 2 operations and the exponentiation in subspace requires roughly k 3 N e 3 operations.The total computational cost to perform k K -th order block Krylov expansion increases by the factor i k j j j y { z z z ( ) compared to the k T -th order Taylor expansion.This can add computational overhead for small basis sets, but for larger basis sets the computational overhead is usually negligible, especially when compared to other more expensive tasks in AFQMC. 69erformance Comparison.We compare the different matrix-exponential methods for a short ph-AFQMC simulation using 10 time steps and 2 400 walkers and calculate the error in the ph-AFQMC total energy induced by the errors in calculating the matrix exponential.Then we determine the minimal number of A ̂|Ψ⟩ operations k min required to achieve a specified accuracy ε = 10 −5 Eh.The exact diagonalization described at the beginning of the section yields the reference energies.We average k min over 22 different systems for time steps from τ = 0.01 Eh −1 to τ = 0.3 Eh −1 and for two different basis sets (cc-pVDZ and cc-pVTZ).In ref 109, we provide the list of considered compounds, ranging from small molecules to benzene-like cyclic compounds.
Figure 1 compares the different exponentiation methods and demonstrates that Krylov methods need less A ̂|Ψ⟩ operations than the Taylor and the Chebyshev expansion.Furthermore, k min increases with the time step and this is more pronounced for the Taylor and the Chebyshev expansion than for the Krylov methods.k min is also more sensitive to the system for the Taylor and the Chebyshev expansion leading to larger error bars for these methods.Hence, one needs to verify the order of the expansion for every system in particular for larger timesteps.Finally, the larger cc-pVTZ basis set needs slightly larger k min than the smaller cc-pVDZ basis set.

■ MODIFIED AFQMC ALGORITHM
In this section, we present modifications of the commonly used AFQMC algorithm 69,83,88 that allow large time steps in AFQMC and reliable time step extrapolation.
Local-Energy Capping.Inspired by Zen et al., 90 we implement the local energy capping as Here, 0 is the best estimate of the ground-state energy and the energy window is E N N e . This window increases with the number of electrons as opposed to E 2 = typically used in AFQMC. 71We also add a term proportional to to avoid excessively strong capping for large time steps.The energy capping introduces a bias, potentially affecting the accuracy of the time step extrapolation.To strike a balance, we opt to cap the local energy to control its variance only�a parameter far more sensitive to rare events than the mean value.
No Hybrid-Energy Capping.The energy capping is usually applied for the hybrid energy (eq 12), too.We turn off hybrid energy capping because it strongly affects the time-step errors.In particular, we observe that hybrid-energy capping can introduce strong size-consistency errors.We will discuss this in more depth later in this section.
Reweighting Factor Capping.To compensate for the lack of hybrid energy capping, we restrict the reweighting factor in eq 0 else This is as effective as the hybrid-energy capping and does not cause size-consistency errors.
Force-Bias Capping.Force-bias capping is essential for stable AFQMC simulations.We recommend-capping the force bias with to minimize the fluctuations and time-step errors.We find that removing the force-bias components with large absolute values is more reliable for large time steps than the conventional choice of capping them at an absolute value of 1.
No Energy Mixing.Many implementations modify eq 11 to include a mixture E H mix of the current and previous hybrid energy 85,93 However, we find that energy mixing deteriorates the timestep extrapolation, so we do not use it in this work.
Accurate Matrix Exponentials.Reliable time-step extrapolation requires minimizing remaining errors caused by approximating the exponential of the effective interaction int .We recommend using either the Krylov projection method with k = 5 or the Block-Krylov method with k = 4 for this purpose (see Sec. IV).
Size-Consistency Errors.To quantify the effect of the listed modifications, we test the size consistency of the algorithms similar to Zen et al. 90 and Lee et al. 92 The former group considers a CH 4 −H 2 O dimer at a separated C−O distance of 11.44 A using the cc-pVDZ basis set.We measure the sizeconsistency error using Lee et al. 92 investigates a chain of N 2 molecules at large separation.They considered the AFQMC correlation energy per N 2 molecule E corr /k as a measure of size consistency.
We expect that E s is less than 0.01 mEh for size-consistent methods.For the CH 4 −H 2 O dimer (Figure 2, top), the standard algorithm and our modifications yield small size-consistency errors E s at small time steps τ = 0.05 Eh −1 .However, the standard algorithm exhibits sizable errors at larger time steps (≈ 0.4 mEh at τ = 0.20 Eh −1 ), whereas the modified algorithm is sizeconsistent for all time steps.These size-consistency errors directly impact the binding energy where the bonded C−O distance is 0.63 A. Figure 2(bottom) shows similar-size deviations of the energy when the dimer is bound.Hence, the standard algorithm requires much smaller time-steps to achieve the same level of accuracy as the algorithm presented in this work.For the N 2 chain, we revisit the work of Lee et al. 92 Using the traditional algorithm, we can reproduce the size-consistency errors for large time steps 0.05 Eh −1 (see Figure 3).Previously, it has been suggested 92 that this error originates from the cosine projection in the phaseless approximation.However, deactivating the hybrid energy capping suppresses most size-consistency errors.The errors are not larger in magnitude than the sizeconsistency errors at the small time step τ = 0.005 Eh −1 .

Time-Step Extrapolation of ph-AFQMC Energies.
Generally, the time-step errors in the total energy follow a quadratic relation for all propagators and trial wave functions.The quadratic term directly emerges from the analysis of the propagators (Sec.III), whereas the linear term is attributed to the cosine projection (eq 14) and perturbations introduced in the Monte Carlo algorithm to maintain the stability of the random walk (Sec.V.A−D).The top plot in Figure 4 illustrates the time-step errors for the absolute ph-AFQMC energies of the 2C S water dimer with the cc-pVDZ basis set.Notably, the Split-2 propagator exhibits the smallest time-step errors predominantly characterized by the quadratic term.The Taylor propagator shows similar behavior with slightly larger time-step errors.Conversely, both the Crank-Nicolson and the Split-1 propagator show larger time-step errors primarily governed by the linear term.The bottom plot of Figure 4 displays the time-step extrapolation for relative energies, specifically the binding energy of the 2C S water dimer.For relative energies, time-step errors cancel by 1−2 orders of magnitude compared to the absolute energies.Again, the Split-2 and Taylor propagators show the smallest time-step errors with vanishingly small linear terms.For these propagators and this system, the time-step errors are so small (<0.1 mEh) that extrapolation is not necessary and the binding energy could be determined by a single calculation.Conversely, for the Split-1 and Crank-  We compare our data with the data from Lee et al. 92 Reference results at τ = 0.005 Eh −1 (in gray) are size-consistent, while at the larger time step τ = 0.05 Eh −1 , sizeconsistency is maintained only when hybrid energy capping is not employed (in blue).In contrast, our results using hybrid energy capping show significant size-consistency errors (in yellow) and are consistent with previous results (in magenta).
Nicolson propagators, the linear term remains significant, leading to larger time-step errors.Based on these observations, we recommend a general rule for extrapolation of the binding energies: use for propagators where quadratic errors dominate (Split-2 and Taylor), and for those propagators with the predominant linear term (Split-1 and Crank-Nicolson).As a rule of thumb, we observe that timestep extrapolation decreases time-step errors by an order of magnitude.

■ APPLICATIONS AND DISCUSSION
In the following three subsections, we (i) compare accurate ph-AFQMC energies obtained with small time steps 69 to extrapolating the energies from large time-steps, (ii) extrapolate the binding energies of small water clusters in the heavyaugmented double-ζ basis set and in the CBS limit, and (iii) calculate the equilibrium bond length of the N 2 dimer using large time-step ph-AFQMC calculations.
HEAT Set Total Energies.We calculate the total energies for the molecules in the HEAT set 110 using the Split-2 propagator and time-step extrapolation.We include Cholesky vectors up to a threshold of 10 −6 Eh and sample with the Hartree−Fock trial determinant.We use Dunning's cc-pVDZ basis set and the frozen-core approximation (1s shell for B−F atoms).To verify the reliability of the time-step extrapolation, we utilized 4 equidistant time steps ranging between 0.05 Eh −1 and 0.20 Eh −1 .We perform each ph-AFQMC simulation using 2 400 walkers, with 2 000 equilibration steps, and 10 000 time steps during the sampling phase.To have accurate reference values for comparison, we perform small time-step ph-AFQMC simulations with τ = 0.002 Eh −1 , 4 800 walkers, 100 000 equilibration steps, and 400 000 sampling steps.This setup guarantees roughly identical statistical errors for the small timestep and extrapolated energies, but the CPU time to obtain extrapolated energies is reduced by a factor of 20 using large time steps.All calculations are performed using the QMCFort code and a single CPU node (dual-socket Intel Skylake Platinum 8174) per system.Figure 5 compares the large time-step and extrapolated ph-AFQMC energies with the reference ones at τ = 0.002 Eh −1 .The ph-AFQMC energies at τ = 0.10 Eh −1 exhibit sizable time-step errors with a mean error of 1.89 mEh.The time-step errors lead to a consistent under-correlation.Reducing the time step to τ = 0.05 Eh −1 reduces the mean error to 0.53 mEh (already a small fraction of the desired chemical accuracy of 1.59 mEh).The extrapolated energies are accurate with a mean absolute error of 0.31 mEh.The mean signed error of 0.06 mEh shows that the extrapolated energies are not systematically biased compared to the reference small time-step energies.
Water Clusters.As a second example, we apply time-step extrapolation to determine the binding energies of the water clusters.We utilize the Split-2 propagator at three distinct time steps (τ = 0.10, τ = 0.15, and τ = 0.20 Eh −1 ) to extrapolate the binding energies using eq 54.To make reference calculations with a small time step possible, we employ the heavy-augmented double-ζ basis set.Furthermore, we obtain the CBS limit using the aug-cc-pVNZ (N = D, T, Q) basis sets and a 4−5 inverse polynomial extrapolation scheme 111 due to its behavior compared to the standard inverse cubic extrapolation scheme.Here, E CBS is the energy in the complete basis set, N is the largest angular momentum in the aug-cc-pVNZ basis set (N = 2−4 for N = D, T, Q), and E N is the corresponding energy.Each ph-AFQMC calculation is performed using 2 400 walkers and 50 000 time steps, except for the water molecule, where we used 9 600 walkers and 200 000 time steps to minimize binding energy errors.The computational setup is  described in more detail in our previous work. 69As an illustrative example, the largest calculation was performed on the 5CYC cluster using the aug-cc-pVQZ basis set, comprising 855 orbitals, 40 electrons, and 5685 Cholesky vectors.Executed on 4 CPU nodes (dual-socket AMD Epyc 7713), this calculation took 8 days to complete (average performance of 40 GFLOPS/core).Table 1 presents time-step extrapolated ph-AFQMC binding energies for water clusters and compares them with MP2, CCSD(T), and small time-step (τ = 0.01 Eh −1 ) ph-AFQMC values reported in our previous work. 69The time-step extrapolated binding energies are within the error bars in agreement with the small time-step ph-AFQMC values and also agree with MP2 and CCSD(T) values within chemical accuracy.Compared to the small time-step ph-AFQMC calculations, the total number of samples is reduced by a factor of 2, while concurrently reducing the statistical errors, too.We used a simple expression to estimate the statistical error σ x for the extrapolated energies where N τ is the number of time steps used for the extrapolation and σ i denotes the estimated 1σ error for each of the time steps.
Table 2 lists the MP2, CCSD(T), 111 and our time-step extrapolated ph-AFQMC binding energies in the CBS limit.
Small time-step ph-AFQMC calculations are not possible in the CBS limit due to the large computational cost.The mean absolute error between our results and MP2 amounts only to 0.15 kcal/mol.The ph-AFQMC results deviate even less from CCSD(T) with a mean absolute error of 0.07 kcal/mol.
We estimate the statistical errors of the extrapolated binding energies using eq 57.Moreover, for the statistical errors after the basis set extrapolation, we used the corresponding errors observed at the aug-cc-pVQZ basis set.note that the CCSD(T) values are approximated by adding the difference between CCSD(T) and MP2 at the aug-cc-pVDZ basis set to the CBS MP2 binding energies.
As an example, we show the time-step extrapolation of the binding energy for the 3UUD water cluster at the aug-cc-pVQZ basis set in Figure 6(top).Figure 6(bottom) illustrates the basis set extrapolation for the same cluster.
Equilibrium Geometry of the N 2 Molecule.Unlike the total energies, properties such as energy differences, bond lengths, and lattice constants exhibit significantly lower sensitivity to time-step errors.Using the N 2 molecule, we demonstrate that accurate bond lengths are successfully obtained with single-point, large time-step ph-AFQMC calculations.Using the aug-cc-pVQZ set and a RHF trial wave function, we calculate ph-AFQMC energies for five equidistant bond lengths ranging from 0.9 to 1.3 A. We compare the results with a small time step of 0.004 Eh −1 to a larger one of 0.1 Eh −1 .Fitting total energies to the Morse potential energy curve determines the equilibrium bond length R 0 , the dissociation energy D, and the parameter a controlling the width of the potential.
Figure 7 shows both ph-AFQMC and the reference CCSD(T) potential energy curves.The absolute CCSD(T) and small time step ph-AFQMC energies agree very well, however, the large time step ph-AFQMC energies are shifted by 10 mEh−14 mEh (inset in Figure 7).The slight increase of the ph-AFQMC time-step errors with increasing bond length may result from the deteriorating quality of the trial RHF wave function.Nevertheless, the observed change in time-step errors is not large enough to compromise the accuracy of the predicted bond length; the predicted bond lengths listed in Table 3   Frozen-core approximation is employed.MP2, CCSD(T), and timestep extrapolated ph-AFQMC values in kcal/mol are reported.Figure 6.(Top) Time-step extrapolation of the binding energy for the 3UUD water cluster at the aug-cc-pVQZ basis set.Eq 54 is used to estimate the zero time-step limit.(Bottom) Basis set extrapolation of the binding energy for the 3UUD water cluster using aug-cc-pV(D,T,Q) basis sets and 4−5 inverse polynomial extrapolation technique.the computational time for each calculation by a factor of τ L /τ S = 25.

■ CONCLUSION
We implemented modifications of the ph-AFQMC algorithm in the QMCFort code to reduce the size-consistency errors at large time steps.The modifications are mainly related to the control of rare events and the weight update procedure.Addressing sizeconsistency errors is crucial as they have a consequential impact on the accuracy of binding energies.Using the CH 4 −H 2 O dimer (Figure 2) and chains of infinitely separated N 2 molecules (Figure 3), we showed that the modified approach markedly diminishes size-consistency errors across all time steps.Revisiting results of Lee et al., 92 we could identify hybrid energy capping as the primary source of size-consistency errors.
Employing large time steps required better control of the errors in the matrix exponentiation of the one-body operators.Utilizing a diverse set of 22 different systems using Dunning's cc-pVDZ and cc-pVTZ basis sets, we demonstrated that the common sixth-order Taylor expansion may prove insufficient.Larger time steps require on average up to 10−12 | operations for adequate accuracy.In contrast, Krylov methods� particularly a blocked variant �are more robust, with negligible time-step errors while employing merely 4 | operations even for large time steps.The block-Krylov method is also less sensitive to system and basis set choice than the Taylor expansion.Moreover, it introduces negligible computational overhead, especially for larger systems, making it a promising alternative.
While all AFQMC propagators yield consistent results for small time steps, we have shown that the Split-2 propagator is the optimal propagator within the AFQMC formalism, with the leading error term dictated by the Hubbard-Stratonovich transformation.We demonstrated that reliable time-step extrapolation is possible for all propagators via a second-order polynomial fit.Notably, the Split-1 and Crank-Nicolson propagators exhibit significant linear time-step errors compared to the Taylor and Split-2 propagators that are dominated by quadratic errors.For the HEAT set molecules, using a Split-2 propagator reduced the computational cost by an order of magnitude while retaining almost the same accuracy (mean absolute error of mEh) as the small time-step reference.For water clusters, the time-step extrapolated binding energies at the heavy-augmented double-ζ basis set agree well with those obtained by small time-step calculations.Similarly, the complete basis-set time-step extrapolated binding energies are in excellent agreement with the reference CCSD(T) values (root-meansquare error of 0.08 kcal/mol).
Finally, we showed that certain observables, such as the equilibration bond length, are much less sensitive to time-step errors than total energies.Using the N 2 dimer as an example, we demonstrated that the ph-AFQMC calculation with τ = 0.1 Eh −1 yields the same bond length as a calculation with τ = 0.002 Eh −1 , while reducing the computational cost by a factor of 25.
In conclusion, we highlighted the ability to employ larger time steps in ph-AFQMC, resulting an order of magnitude speedup compared to standard ph-AFQMC calculations.We plan to validate the efficacy of large time-step ph-AFQMC on a wider range of systems where accurate small time-step computations are very expensive.

Notes
The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS
This research was funded in whole by the Austrian Science Fund (FWF) within the project 10.55776/P33440.All calculations were performed on the VSC4/VSC5 (Vienna scientific cluster).
that are the most expensive part of the procedure.Consequently, the optimal method for a fixed time step τ is the one that minimizes the number of w w int | operations.To simplify the notation, we will use A w int = throughout the remainder of this section.

Figure 1 .
Figure 1.Number of A|Ψ⟩ operations, i.e., the order of the method k min required to achieve an accuracy of ε = 10 −5 Eh in the ph-AFQMC energy for different matrix-exponentiation methods.The local energy is averaged over 10 time steps and 2 400 walkers.k min is averaged over 22 systems.Note that the error bars correspond to the standard deviation of the average over the 22 systems and not the statistical fluctuations of the AFQMC simulation.Krylov methods (nonblocked, brown triangles; blocked, magenta diamonds) outperform the Taylor (blue circles) and the Chebyshev (yellow squares) expansion typically utilized in ph-AFQMC codes.The cc-pVDZ basis set (top) needs fewer A|Ψ⟩ operations than the cc-pVTZ basis set (bottom).

Figure 2 .
Figure 2. (Top) The size-consistency error E s and (Bottom) the binding energy E b for a CH 4 −H 2 O dimer using a cc-pVDZ basis set.We compare the standard ph-AFQMC algorithm (yellow squares), our modified algorithm (blue circles), and CCSD(T) results (magenta line).

Figure 3 .
Figure3.Correlation energy per N 2 molecule as a function of the number of infinitely separated N 2 molecules.We compare our data with the data from Lee et al.92 Reference results at τ = 0.005 Eh −1 (in gray) are size-consistent, while at the larger time step τ = 0.05 Eh −1 , sizeconsistency is maintained only when hybrid energy capping is not employed (in blue).In contrast, our results using hybrid energy capping show significant size-consistency errors (in yellow) and are consistent with previous results (in magenta).

Figure 4 .
Figure 4. Time-step errors in the total energy (top) and the binding energy (bottom) of the 2C S cluster for different AFQMC propagators.Time steps <0.05 Eh −1 are used for the linear extrapolation of the binding energy with the Crank-Nicolson (brown triangles) and the Split-1 (magenta diamonds) propagators, while all time steps are used for the quadratic extrapolation with the Taylor (yellow squares) and the Split-2 (blue circles) propagator.All calculations use the frozen-core approximation and the cc-pVDZ basis set.Figure 5. Comparison of time-step extrapolated and large time-step energies (τ = 0.05 Eh −1 and τ = 0.10 Eh −1 ) against reference values for 26 molecules in the HEAT set.Reference values are obtained using small time-step ph-AFQMC with τ = 0.002 Eh −1 .All computations use the cc-pVDZ basis set and the frozen-core approximation.

Figure 5 .
Figure 4. Time-step errors in the total energy (top) and the binding energy (bottom) of the 2C S cluster for different AFQMC propagators.Time steps <0.05 Eh −1 are used for the linear extrapolation of the binding energy with the Crank-Nicolson (brown triangles) and the Split-1 (magenta diamonds) propagators, while all time steps are used for the quadratic extrapolation with the Taylor (yellow squares) and the Split-2 (blue circles) propagator.All calculations use the frozen-core approximation and the cc-pVDZ basis set.Figure 5. Comparison of time-step extrapolated and large time-step energies (τ = 0.05 Eh −1 and τ = 0.10 Eh −1 ) against reference values for 26 molecules in the HEAT set.Reference values are obtained using small time-step ph-AFQMC with τ = 0.002 Eh −1 .All computations use the cc-pVDZ basis set and the frozen-core approximation.

Figure 7 .
Figure 7. Potential energy curves E obtained for CCSD(T), ph-AFQMC with τ = 0.004 Eh −1 , and ph-AFQMC at τ = 0.1 Eh −1 using aug-cc-pVQZ basis set.Five equidistant bond lengths in the range from 0.9 to 1.3 A are used to fit the Morse potential.The equilibrium bond lengths agree closely with the experimental value.Furthermore, the inset illustrates the behavior of ph-AFQMC time-step errors ΔE with increasing bond length.

Table 1 .
are nearly identical and align well with the experimental value within statistical accuracy.Using large-time step ph-AFQMC reduces Binding Energies of the Four Most Stable Water Clusters Calculated Using a Heavy-Augmented cc-pVDZ Basis Set and All-Electron Wavefunctions a

Table 3 .
Equilibrium Bond Lengths of the N 2 Molecule for Different Levels of Theory: CCSD(T), ph-AFQMC with Small Time Step, and the ph-AFQMC at Large Time Step