Chemical Transformations Approaching Chemical Accuracy via Correlated Sampling in Auxiliary-Field Quantum Monte Carlo
James Shee, Shiwei Zhang, David R. Reichman, Richard A. Friesner
CChemical Transformations Approaching Chemical Accuracy via Correlated Samplingin Auxiliary-Field Quantum Monte Carlo
James Shee, ∗ David R. Reichman, and Richard A. Friesner
Department of Chemistry, Columbia University, 3000 Broadway, New York, NY, 10027
Shiwei Zhang
Department of Physics, College of William and Mary, Williamsburg, Virginia 23187-8795
The exact and phaseless variants of Auxiliary-Field Quantum Monte Carlo (AFQMC) have beenshown to be capable of producing accurate ground-state energies for a wide variety of systemsincluding those which exhibit substantial electron correlation effects. The computational cost ofperforming these calculations has to date been relatively high, impeding many important applica-tions of these approaches. Here we present a correlated sampling methodology for AFQMC whichrelies on error cancellation to dramatically accelerate the calculation of energy differences of rel-evance to chemical transformations. In particular, we show that our correlated sampling-basedAFQMC approach is capable of calculating redox properties, deprotonation free-energies, and hy-drogen abstraction energies in an efficient manner without sacrificing accuracy. We validate thecomputational protocol by calculating the ionization potentials and electron affinities of the atomscontained in the G2 Test Set, and then proceed to utilize a composite method, which treats fixed-geometry processes with correlated sampling-based AFQMC and relaxation energies via MP2, tocompute the ionization potential, deprotonation free-energy, and the O-H bond disocciation energyof methanol, all to within chemical accuracy. We show that the efficiency of correlated samplingrelative to uncorrelated calculations increases with system and basis set size, and that correlatedsampling greatly reduces the required number of random walkers to achieve a target statistical error.This translates to CPU-time speed-up factors of 55, 25, and 24 for the the ionization potential of theK atom, the deprotonation of methanol, and hydrogen abstraction from the O-H bond of methanol,respectively. We conclude with a discussion of further efficiency improvements that may open thedoor to the accurate description of chemical processes in complex systems.
I. INTRODUCTION
The pursuit of chemical accuracy (as defined by errorsnot exceeding 1 kcal/mol) in the ab initio computation ofthe energetic properties of generic many-electron systemsis a long-standing goal that has yet to be reached.[1–5] Traditional methods such as full configuration in-teraction (FCI)[6] and single-reference coupled clusterwith single, double, and perturbative triple excitations(CCSD(T))[7, 8] scale exponentially and with the seventhpower of the system size, respectively, and the latter canbreak down for systems which exhibit sufficiently strongelectron correlation.[9, 10] Multi-reference methods suchas CASSCF,[11–13] often supplemented with second or-der perturbation theory,[14, 15] and FCI-QMC[16–19]have been shown to yield benchmark-quality results evenfor strongly correlated electronic systems, an import classof problems that includes many transition metal (TM)-containing systems.[20, 21] However, while a judiciouschoice of the active space can make such calculationstractable for small systems, exponential scaling prohibitsthe use of these methods in studying most realistic sys-tems of interest in biology, materials science, and chemi-cal catalysis. As a result, with the exception of isolatedstudies using more accurate, specialized methods,[22–24]these systems can only be investigated feasibly with less ∗ [email protected] accurate but more economical approaches such as Den-sity Functional Theory (DFT).[25] Given the plethoraof investigations using this tool to generate hypothesesconcerning reaction mechanisms, equilibrium structures,etc., the importance of having a systematically improv-able ab initio method which is both accurate and feasiblebecomes readily apparent.Among the class of Quantum Monte Carlo (QMC)methods used in electronic structure theory,[26–28]Auxiliary-Field QMC, and in particular its “phaseless”variant (ph-AFQMC)[29, 30] which controls the signproblem at the cost of introducing a bias which in princi-ple can be systematically reduced, has produced state-of-the-art, benchmark-quality results for systems suchas first- and second- row d elements,[31] the chromiumdimer,[9] cobalt adatoms on graphene,[32] and vari-ous transition metal oxides,[33] while exhibiting low-polynomial scaling ( M for Gaussian basis sets). How-ever, the widespread use of ph-AFQMC in quantumchemistry has not yet taken place. One major reasonis undoubtedly its relatively high computational cost, asthe favorable scaling is masked by a large prefactor.In this work we present an approach to greatly re-ducing this prefactor which involves the use of corre-lated sampling for a particular class of important pro-cesses. The general idea is that for sufficiently simi-lar systems, energy differences are expected to convergemore rapidly, i.e. with smaller error bars, than total en-ergies when the errors or statistical fluctuations in the a r X i v : . [ phy s i c s . c h e m - ph ] M a y calculations are biased in the same direction. Indeed, er-ror cancellation is largely responsible for the success ofmany approximate methods such as DFT in the com-putation of energy differences. Correlated sampling haspreviously been adapted to reduce the statistical errorsin QMC approaches via the use of the same set of con-figurations sampled for both the primary and secondarysystems.[26] This technique is often referred to in theliterature as “Differential” QMC, and the details of itsimplementation can vary depending on the type of QMCbeing used. The potential energy curves of H and BHhave been calculated using correlated sampling with Vari-ational QMC,[34] and similarly that of the H cationwith Differential Green’s Function QMC.[35] The lat-ter method has been used to compute the dipole mo-ment of LiH, [36] and to calculate infinitesimal energydifferences from which forces and various polarizabilitieshave been obtained.[37] This idea has also been extendedto Diffusion Monte Carlo, which has been used to com-pute forces and potential energy surfaces for the first rowdiatomics.[38] Correlated sampling has also been used tocalculate energy differences between ground and excitedstates of the same Hamiltonian, as illustrated by a Varia-tional QMC study of particle-hole excitations in the two-dimensional electron gas.[39] In addition, the concept hasbeen extended to enable concerted propagation of a sys-tem with different time steps, in order to extrapolate theTrotter error in Differential Diffusion QMC.[40]Correlated sampling is, in fact, also well-suited tomodel the energetics of myriad chemical reactions, since only energy differences, as opposed to total energies,are relevant. In this paper we present a novel corre-lated sampling-based AFQMC approach, and show thatit is capable of computing ground-state energy differ-ences corresponding to redox, deprotonation, and hy-drogen abstraction reactions to a given statistical er-ror in a fraction of the time previously required, with-out any loss of accuracy. Redox reactions (often in-volving TMs) abound, for example, in metabolic andphotosynthetic processes,[41–43] battery chemistry,[44]and catalysis ( e.g. CO and O reduction).[45, 46]A reliable ab initio method to calculate deprotonationfree energies would provide an improvement upon ex-isting computational approaches to determining pKa’sand protonation states,[47–49] which would have signif-icant ramifications for drug discovery,[50, 51] materialsscience,[52–54] and the structural determination of bi-ological complexes.[55, 56] In the context of chemicalcatalysis, proton removal is known to be the rate-limitingstep in many important reactions, e.g. oxygen reduc-tion on the surface of TiO .[57] Hydrogen abstractionreactions are ubiquitous and play a major role in com-bustion and the oxidation of hydrocarbons,[58, 59] di-amond growth via chemical vapor deposition,[60] bio-chemical processes involving e.g. the antioxidant vita-min E [61] and various metalloenzymes,[62, 63] indus-trial processes,[64] and organic synthesis.[65] AFQMCwould be a useful benchmark for previous ab initio studies[63, 66–73] in predicting thermodynamic proper-ties of this difficult class of chemical reactions. Thus theclass of applications we target is large and important. Wehighlight the fact that the cost of our correlated samplingapproach, relative to the uncorrelated method, decreaseswith increasing system and basis set size, opening thedoor to the treatment of realistic large chemical systemswith correlated sampling-based AFQMC in the near fu-ture.This work is organized as follows: Section II.A will re-view the exact and phaseless variants of AFQMC, whileSection II.B will present our correlated sampling ap-proach. We justify this approach for modeling molecularsystems in Section II.C, and Section II.D will disclosefurther computational details. Section III.A will presentcalculations of the ionization potentials (IPs) and elec-tron affinities (EAs) of the 1st row atoms in the G2 IonTest Set,[74] while Section III.B will consider adiabaticmolecular properties, taking the IP, deprotonation en-ergy, and O-H bond dissociation energy of methanol asexamples. The efficacy of correlated sampling as a func-tion of both basis set size and number of random walkerswill be explored in Sections III.C. The latter subsectionwill provide an assessment of the reduction in CPU-timeafforded by the use of correlated sampling. In Section IVwe conclude, emphasizing opportunities for further gainsin computational efficiency that may be possible and fu-ture targets of investigation. II. METHODSA. Overview of AFQMC - Exact Method and thePhaseless Constraint
The excited states of a many-body state | Φ (cid:105) canbe projected out via imaginary time propagation, i.e. lim N →∞ ( e ∆ τ ( E − ˆ H ) ) N | Φ (cid:105) → | Φ (cid:105) , with (cid:104) Φ | Φ (cid:105) (cid:54) = 0 andthe general electronic Hamiltonianˆ H = M (cid:88) ij T ij (cid:88) σ c † iσ c jσ + 12 M (cid:88) ijkl V ijkl (cid:88) σ,τ c † iσ c † jτ c lτ c kσ , (1)where M is the size of the orthonormal one-particle ba-sis, and c † iσ and c iσ are the second-quantized fermioniccreation and annihilation operators with particle andspin labels. The two-body matrix elements, V ijkl , canbe expressed in terms of Cholesky vectors as V ijkl = (cid:80) α L αik L αjl .[75] Defining the one-body operator ˆ v α ≡ i (cid:80) ik L αik (cid:80) σ c † iσ c kσ , and subtracting the expectationvalue with respect to the trial wavefunction (cid:104) ˆ v α (cid:105) fromˆ v α , the Hamiltonian can be written as the sum of allone-body operators, ˆ H , plus the following two-body op-erator ˆ H = − (cid:88) α (ˆ v α − (cid:104) ˆ v α (cid:105) ) . (2)Use of the Trotter-Suzuki decomposition[76, 77] gives e − ∆ τ ˆ H = e − ∆ τ ˆ H / e − ∆ τ ˆ H e − ∆ τ ˆ H / + O (∆ τ ) . (3)The exponential terms involving ˆ H may be decomposedusing a Hubbard-Stratonovich (HS) transformation[78,79] e ∆ τ (ˆ v α −(cid:104) ˆ v α (cid:105) ) = (cid:90) ∞−∞ dx α (cid:18) e − x α √ π (cid:19) e √ ∆ τx α (ˆ v α −(cid:104) ˆ v α (cid:105) ) , (4)which expresses the exponential of a two-body operatoras the exponential of a one-body operator integrated overauxiliary-fields (AFs). This transformation allows forpractical propagation in terms of the Thouless theorem,which states that the application of an exponential of aone-body operator on a Slater determinant produces an-other Slater determinant,[80] which can be implementedvia a simple matrix multiplication.[81, 82]The propagator (3) now takes the form of a multi-dimensional integral e − ∆ τ ˆ H = (cid:90) d x P ( x ) ˆ B ( x ) , (5)where x = ( x , x , . . . , x α ), P ( x ) is a normaldistribution with unit variance, and ˆ B ( x ) = e − ∆ τ ˆ H e √ ∆ τ x · ( ˆv −(cid:104) ˆv (cid:105) ) e − ∆ τ ˆ H .The integral in (5) may be approximated using a MonteCarlo scheme, with walkers whose propagation in thespace of Slater determinants is guided by the compleximportance function (cid:104) φ T | φ (cid:105) which is proportional to thewalker weights. The representation of the total wavefunc-tion is thus a weighted sum over walker determinants | Φ (cid:105) = (cid:88) k w k | φ k (cid:105)(cid:104) φ T | φ k (cid:105) , (6)yielding essentially a multi-reference description. The en-ergy is calculated at intervals using the mixed-estimator (cid:104) φ T | ˆ H | Φ (cid:105)(cid:104) φ T | Φ (cid:105) = (cid:80) k w k E L ( φ k ) (cid:80) k w k , (7)where the “local energy” is given by E L ( φ k ) = (cid:104) φ T | ˆ H | φ k (cid:105)(cid:104) φ T | φ k (cid:105) ,in analogy to Diffusion QMC.[83]The method as described above is called the “FreeProjection” (FP) approach.[84] While it is formally ex-act and can yield excellent results for small system sizes,this method suffers from the “phase problem,” which isa generalization of the Fermionic “sign problem” to thecomplex plane.[29, 85–87] For the standard Coulomb in-teraction the ˆ v α operators are purely imaginary, and eachapplication of e √ ∆ τ x · ( ˆv −(cid:104) ˆv (cid:105) ) can be thought of as a ro-tation of the Slater determinant | φ (cid:105) , causing an evolu-tion of the overlap (cid:104) φ T | φ (cid:105) in the complex plane. Overthe course of the random walk, a determinant accumu-lates a phase e iθ , and the infinitely many possible values of θ (cid:15) [0 , π ) result in the possibility of infinitely manyindistinguishable determinants. Furthermore, over thecourse of the propagation, walkers will populate the ori-gin where (cid:104) φ T | φ (cid:105) = 0, and subsequent propagation yieldsonly noise from signal cancellation effects and divergencesin the weights and local energies.The ph-AFQMC employs a multifaceted strategy tocontrol this problem. The weights are initialized to apositive real constant, and after each propagation stepare projected back onto the real axis, i.e. the ro-tated weights are multiplied by max { , cos (∆ θ ) } , wherewe have defined the phase of the overlap ratio ∆ θ ≡ Im { ln (cid:104) φ T | φ ( τ +1) (cid:105)(cid:104) φ T | φ ( τ ) (cid:105) } .[88] For this phase projection to work,the AFs are shifted by a force bias (FB) ¯x ,[30, 89] theoptimal choice of which is obtained by minimizing thefluctuations of the weights with respect to the AFs attheir average value: ∂∂x α (cid:20) (cid:104) φ T | e √ ∆ τ ( x α − ¯ x α )(ˆ v α −(cid:104) ˆ v α (cid:105) ) | φ k (cid:105)(cid:104) φ T | φ k (cid:105) e − ¯ x α / x α ¯ x (cid:21) x α =0 = 0 . (8)This is, in essence, a stationary-phase approximation.Expanding the expression inside the brackets to O ( √ ∆ τ )and taking the derivative gives¯ x α = −√ ∆ τ (cid:2) ¯ v α − (cid:104) ˆ v α (cid:105) (cid:3) , (9)where ¯ v α ≡ (cid:104) φ T | ˆ v α | φ (cid:105)(cid:104) φ T | φ (cid:105) . The introduction of the FB does not add any additional approximations, as the integra-tion variable in (4) is merely shifted by a constant, yetit is crucial for two reasons. First, it diverges when the“nodal surface” (as defined in the complex plane of over-laps) is approached, pushing the walker away from theorigin. Second, since ˆ v α is complex, Im [ ¯x ] reduces theamount of physical information discarded in the phaseprojection. Similarly, the subtraction of (cid:104) ˆ v (cid:105) from ˆ v in thepropagator also greatly reduces the severity of this pro-jection, as the smaller diagonal matrix elements of theresulting propagator cause milder rotations of the phasesof the orbitals.[90]The choice of FB allows the weight factor which mul-tiplies the previous weight after a propagation step tobe written as W ( φ ) = e − ∆ τE L ( φ ) . The second approx-imation in ph-AFQMC, which is much milder than thefirst, takes the real part of the local energy in the weightabove: E L ( φ ) ≡ Re { (cid:104) φ T | ˆ H | φ (cid:105)(cid:104) φ T | φ (cid:105) } . (10)The severity of the phaseless constraint in ph-AFQMCcan be reduced with the use of more accurate trialwavefunctions, especially those with the correct sym-metry properties.[84] This can be seen from (10), forwhen | φ T (cid:105) = | Φ (cid:105) , the local energy which determinesthe weights and energy measurements is a real constant(equal to the exact ground state energy). B. Correlated Sampling Methods for AFQMC
Given that the statistical error in an AFQMC calcula-tion arises solely from the MC evaluation of the integralover AFs, a natural way to implement correlated sam-pling in the calculation of an energy difference is to pairthe walkers of the two systems, and use the same set ofAFs to propagate each pair of walkers. To be precise, thecorrelation is established on the level of each term in (2),for each value of α . This, in essence, matches Choleskyvectors of the same iteration, and becomes more effectivethe more similar the interactions are in the two systems.Sampling in this way, however, requires one to re-linquish optimal importance sampling, which is usuallyimplemented via a population control (PC) scheme inwhich walkers with large (small) weights are replicated(stochastically purged), since performing independentPC for each system separately would quickly destroythe walker-pair correlation. Alternatively, one can im-plement PC with respect to the weights of the primarysystem for both systems . This, however, will only be ef-fective if the two systems are essentially identical. In theabsence of a PC scheme, the noise from the accumula-tion and persistence of divergent walkers inevitably growswith propagation time. We find that the immediate re-duction in statistical error following the correlation ofthe AFs, augmented as needed by what we call the “pre-liminary equilibration scheme” below, allows convergedaverages to be obtained at short and intermediate pro-jection times. In light of the fact that the stability of therandom walks at long times, as afforded by a branchingscheme, appears to be only marginally relevant when ourcorrelated sampling approach is used, we simply choose not to implement PC when the AFs are correlated. Inthe event that a walker’s weight becomes zero or negative(in the latter case the phase projection sets the weightto zero), the walker is no longer propagated and the ran-dom number stream is updated if necessary such that thecorrelation of the other walker-pairs is unaffected.The data shown in blue in Fig. 1 illustrates featuresassociated with the typical correlated sampling protocolused in this work. In this example, the propagation inimaginary time (measured in Ha − ) is performed usingcorrelated AFs and repeated 11 times using different ran-dom number seeds. The mean energy difference and asso-ciated standard error at each τ point is computed amongthe repeats (see left plot). We then choose the imaginary-time at which the energy difference is seen to stabilize (inthis case at τ ∼ cu-mulative average at each τ , which represents the runningaverage of the energy measurements taken after the endof the equilibration period up to the given value of τ . Toobtain the final result, we compute the mean and stan-dard error of the cumulative averages among the repeats(see right plot), and choose the value corresponding tothe τ at which the standard error reaches a minimum (ifthere are multiple equivalent minima the one occurringearliest is chosen) or falls below a target error level. We note that our choice of 11 repeats is arbitrary; a largernumber could be used to reduce the standard error asnecessary. Clearly in this example correlating the AFsfor τ ≥ τ = 80 ph-AFQMCruns of the neutral and cationic species. This benchmark,indicated by the solid black line, was found to have a neg-ligible standard error of 0.2 mHa after employing a re-blocking analysis which corrects for auto-correlation.[91]As it is often the case that a very small population sizeis sufficient to achieve a desired statistical error via thecorrelated sampling approach, we choose not to employPC in the uncorrelated comparisons since this would re-sult in a bias that typically goes as 1/ N wlk .[92] A secondreason is that over the relatively short imaginary-timescales relevant for the correlated sampling method, PCis expected to have little effect on the uncorrelated com-parisons as the walker weights usually do not stray farfrom unity. This is confirmed by the similarity of theerror curves plotted in the insets of Fig. 1, correspond-ing to uncorrelated runs with (dotted green) and without(red) PC. The former is obtained by using a large enoughpopulation size such that the bias from the PC algorithmis negligible (360 walkers per repeat), and rescaling theresulting standard error by (cid:113) .If one or both of the comparative systems requires along equilibration time, which can be the case for, e.g. ,initial populations which are severely spin-contaminatedor for strongly correlated systems in which the “guiding”trial function poorly describes the true ground state, thenwalker pair correlation can be lost prior to convergence,and the associated noise growth can make measurementsimpossible. The simplest way to overcome this problemis to use a better trial function, e.g. a multi-determinantCASSCF wavefunction instead of the single Hartree-Fock(HF) determinant. Alternatively, we have devised whatwe will refer to as the “preliminary equilibration scheme”(PES). First, one of the two systems is equilibrated us-ing PC for the required interval, then the walkers of thesecondary system are initialized with the resulting deter-minants of the equilibrated primary system, with weightsscaled by (cid:104) φ secondaryT | φ (cid:105)(cid:104) φ primaryT | φ (cid:105) (and any resulting phases pro-jected) to reflect the appropriate importance sampling.Finally, both systems are propagated using correlatedsampling with their respective Hamiltonians without PCfor a short period of time after which energy measure-ments are collected. A schematic of this procedure isshown in Fig. 2. This protocol is repeated with differentrandom number seeds to obtain satisfactory statistics onthe energy difference between the ground states. We notethat a similar scheme was published many years ago byTraynor et al.[35]Inspection of the ph-AFQMC propagatorˆ B ( x − ¯x ) = e − ∆ τ ˆ H / e √ ∆ τ ( x − ¯x )( ˆv −(cid:104) ˆv (cid:105) ) e − ∆ τ ˆ H / , (11)demonstrates that stricter approaches to correlated sam- (a) Averaged IPs (circles) among the repeats at each τ along the imaginary-time propagation.(b) Mean values (circles) of the cumulative averages takenfor τ > τ = 0 .
01, 12walkers per repeat, and a HF reference state. The error bars givethe standard errors (the standard deviation times √ N r , where N r is the number of repeats) of the mean values among the repeatsat each τ point in (a), and of the cumulative averages in (b).These standard errors are plotted in the insets, along with thescaled standard error resulting from an uncorrelated run in whichPC was used with 360 walkers per repeat (dotted green). pling exist. In one such alternative to correlating only theAFs, each pair of walkers is propagated using both thesame AFs and FBs. However taking the simple average ofthe FBs of the primary and secondary systems was foundto cause an early onset of noise even for marginally differ-ent systems. This issue likely arises from divergences inthe local energies and weights which can occur more fre-quently due to the use of sub-optimal FBs, a fact whichfollows from the discussion centered around Eq. (9). Sep-arately, we choose not to correlate the (cid:104) ˆv (cid:105) since any gainin sampling efficiency would come at the cost of an in-creased bias from the phaseless constraint.[90]Finally, we note that our method of correlated sam-pling is not limited to the phaseless version of AFQMC, FIG. 2: Schematic of the PES version of correlated sampling.Preliminary equilibration of the primary system is shown in red.This phase can be as long as necessary due to the use of PC. Theresulting walkers are used to initialize the secondary system, afterwhich the two systems are propagated with their ownHamiltonians using correlated sampling. Equilibration of thesecondary system is relatively rapid, allowing measurements to betaken before the noise growth becomes prohibitive. and can be used in exactly the same manner for FP cal-culations. In fact, using the same AFs in the latter caseresults in a relatively stricter form of correlated sampling,since the FB is not present in the post-HS propagatorused in FP.
C. Utilizing Optimal Correlation in MolecularApplications
Having justified our choice to correlate only the AFs,we claim that maximal correlation between a pair ofwalkers is achieved when the primary and secondary sys-tems use the same set of basis functions. To see this,recall the definition of the two-electron matrix elementsin the Hamiltonian (1): V ijkl = (cid:90) d r d r φ i ( r ) φ j ( r ) 1 r φ k ( r ) φ l ( r ) . (12)When { φ } primary = { φ } secondary , the V ijkl and, in turn,the one-body operators ˆv in the propagators (11) forboth systems will be identical. This ideal condition isobviously satisfied when the two systems are an atomor molecule with the same geometry but e.g. differentcharges, since atomic basis functions such as Gaussian-type orbitals are usually centered on the positions of thenuclei. In fact, for the calculation of vertical IPs andEAs, exactly the same T ij and V ijkl elements are usedin the imaginary-time propagation of both systems. Foradiabatic redox processes, however, the ground-state ge-ometries of the neutral and ionized species are, in gen-eral, different. In most cases the differences in the V ijkl are slight, and correlated sampling is still found to beeffective, albeit with reduced efficiency. For example,Fig. 3 compares the correlated and uncorrelated stan-dard errors corresponding to the vertical and adiabaticIPs of methanol. While the errors from the uncorre-lated runs are roughly similar in magnitude, that fromthe correlated vertical case is significantly smaller andmore constant compared to that from the correlated adi-abatic case. This suggests the use of a two-step processto compute adiabatic energy differences, in which thefixed-geometry transition is calculated with correlatedsampling-based AFQMC, and the geometry relaxationenergy is obtained from a lower level of theory. Furtherdetails will be presented in Section III.B. (a) Vertical IP(b) Adiabatic IPFIG. 3: Standard errors resulting from uncorrelated (red) andcorrelated (blue) propagation of the repeats for methanol in thecc-pVTZ basis. ∆ τ = 0 .
01 and 24 walkers per repeat were used.
In the same spirit, we have devised a protocol whichenables the calculation of the change in energy corre-sponding to the removal of a proton while the rest ofthe geometry remains fixed, in which the optimal condi-tion mentioned above is realized with the use of so-called“ghost” basis functions. In this scheme, the deprotonatedspecies uses the same set of basis functions as the acid, i.e. the position of the removed proton still serves as acenter for hydrogen basis functions but not as a centerof nuclear charge. As a result, the number of basis func-tions remains the same for the primary and secondarysystems as required by our correlated sampling method.Moreover the V ijkl and thus the ˆv are, by construction,identical. We note that even though the T ij now dif- fer due to the altered electron-nucleus attraction terms,the statistical error in the calculated energy difference is not exacerbated since it results only from the MC eval-uation of the integral over AFs. Fig. 4 illustrates theefficacy of this procedure in a calculation of the deproto-nation energy for methanol, showing a drastic reductionin the statistical errors when the AFs are correlated. The (a) Averaged deprotonation energy (circles) among therepeats at each τ .(b) Mean values (circles) of the cumulative averages takenat τ > OH in the cc-pVTZbasis with “ghost” basis functions. We use ∆ τ = 0 .
01, 12 walkersper repeat, and a HF reference state. The error bars, plotted inthe insets, in (a) represent the standard error among the repeatsat each τ ; those in (b) give the standard error of the cumulativeaverages. “ghost” basis function strategy can also be directly ap-plied effectively to hydrogen abstraction reactions. Fig.5 demonstrates the error reduction afforded in a calcula-tion of the MeOH → MeO · reaction energy.Thus we have presented a correlated sampling ap-proach which has the potential to reduce the statisticalerror for redox, deprotonation, and hydrogen abstractionreactions. The extension of this protocol to obtain adia-batic energy differences will be described in Section III.B. (a) Averaged H abstraction energy (circles) among therepeats at each τ .(b) Mean values (circles) of the cumulative averages takenat τ > · from theO-H bond of CH OH, with the rest of the geometry fixed. Thecc-pVTZ basis is used, with ∆ τ = 0 .
01, 12 walkers per repeat, anda CASSCF reference state. The error bars, plotted in the insets,in (a) represent the standard error among the repeats at each τ ;those in (b) give the standard error of the cumulative averages. D. Computational Details
One-electron and overlap integrals, Cholesky vectors,restricted open-shell and unrestricted HF trial wave-functions were all obtained from a modified version ofNWChem.[75, 93] The maximum residual error in theCholesky decomposition was chosen to be 1 × − Ha.This cutoff has a negligible effect on the QMC energies( i.e. orders of magnitude smaller than the statistical er-ror bars), and is utilized to decrease the number of AFsrequired, which in turn decreases both the error from therandom noise and the computational cost of the calcula-tions.For calculations using a single-determinant trial wave-function with unrestricted reference, spin-contaminationcan occur due to the presence of higher multiplicity spin configurations in the trial function, which can prolongequilibration times and reduce accuracy. We use a spin-projection technique[94] in AFQMC to minimize such ef-fects. In this scheme, the walkers are initialized with arestricted open-shell HF determinant, which is an eigen-function of ˆ S , such that propagation with e √ ∆ τ x · ( ˆv −(cid:104) ˆv (cid:105) ) preserves spin symmetry. The unrestricted HF determi-nant, while not an eigenfunction of the total spin opera-tor, is known in most cases to provide a better descrip-tion of the ground-state energy, and is therefore used toimplement the phaseless constraint and to estimate theenergy.Multi-determinant trial states were obtained fromCASSCF calculations performed with PySCF.[95] Theresulting expansions of the wavefunctions in CI spaceare truncated such that determinants associated with co-efficients below a specified threshold are discarded. Inwhat follows, CASSCF( X , Y ) will denote an active spacewith X electrons and Y orbitals; for cases in whichan energy difference is calculated we use the notation X ≡ N neutral /N ion and Y ≡ M neutral /M ion to denotethe numbers of electrons and orbitals in the neutral andcharged systems, respectively.In all AFQMC calculations the spin -up and -down sec-tors of the walker determinants are separately orthonor-malized via a modified Gram-Schmidt procedure[96] ev-ery 0.05 Ha − to preserve the anti-symmetry of thewalker determinants. When required, we use a PC al-gorithm in which the number of walkers is fixed through-out the entire calculation.[92] PC (when used) and energymeasurements are performed every 0.1 Ha − .Use of the “hybrid” formulation[97] of ph-AFQMC al-lows for the calculation of the local energy, the bottle-neck in the algorithm’s scaling, at intervals rather thanat every step since E L is not required to compute theweight factors, which are now calculated explicitly as (cid:104) φ T | φ ( τ +1) (cid:105)(cid:104) φ T | φ ( τ ) (cid:105) e x · ¯x − ¯x / . We chose to implement the hybridmethod primarily because it is much faster than the localenergy method, and offers more flexibility when devisingcorrelated sampling strategies.The comparison of QMC energies with experimentalresults in general requires extrapolations to mitigate er-rors associated with the use of a finite time step and basisset. The Trotter error, due to the decomposition in (3)and three bounding conditions,[97] can be estimated froma linear extrapolation of the energy differences from threeindependent simulations at ∆ τ = 0 . , . , and 0 . − . For the IPs and EAs of the G2 atomic test set,the energies resulting from the 0.01 Ha − time step werefound in all cases to be well within 1 mHa of the ∆ τ -extrapolated value. Hence, we use only the ∆ τ = 0 . x = 3 (where x is the cardinal number of the basis set); UHF calcula-tions were run at x = 2-5 and MP2 calculations at x = 3,4 using GAMESS[98] or, for the reactions which utilize“ghost” basis functions, NWChem. IPs were computedusing the cc-pVxZ basis,[99] while EAs and the deproto-nation energies to be compared with experimental datawere computed using the aug-cc-pVxZ basis sets.[100]QMC statistical error bars were propagated throughthe data analysis procedure. In assessing the statisticalerror on the deviation of QMC from experiment, we as-sume that there are negligible uncertainties associatedwith 1) the experimental measurements, 2) the exponen-tial (HF) and linear (MP2) fits in the CBS extrapolationprocedure where we found that the former error is an or-der of magnitude or more smaller than the QMC errorin representative cases, and 3) the scaling factor for thezero-point energies. With regard to 1), experimental un-certainties can be found in the NIST database,[101] andare ignored in order to isolate the statistical error associ-ated with the QMC measurements, since the quantitiespresented in this work are differences between calculatedenergies and experimental measurements.Our Fortran 90 code is linked with OpenBLAS[102,103] and Expokit.[104] Random numbers were gener-ated with the 48 bit Linear Congruential Generator withPrime Addend, as implemented in SPRNG5.[105] III. RESULTSA. Atomic IPs and EAs
In this section, we use our correlated sampling-basedAFQMC approach to calculate the IPs and EAs of the1st row atoms included in the G2 Ion Test set, which haveexperimental uncertainties of less than 0.05 eV.[74] Thedeviations from experiment within calculated AFQMCresults are presented in Tables I and II alongside the de-viations resulting from G2 theory,[106] and DFT with theB3LYP exchange-correlation functional[107, 108] in the6-311+G(3df,2p) basis.[74]For the calculations which used the unrestricted HFdeterminant as the trial function, PES was used for theIPs of B, C, N, and O with equilibration times of 35, 20,15, and 15 Ha − , respectively, and similarly for the EAsof C and F with 15 and 10 Ha − , respectively. For the re-maining species, equilibration was facile and thus the AFswere correlated from the beginning of the imaginary-timepropagation. Notable deviations from the experimentalvalues are found for the IPs of B and N, and the EA of Fwhen the UHF state is used as the trial function (nearlyidentical errors have been previously reported within ph-AFQMC[90]). These discrepancies are resolved in thecorrelated sampling-based FP results. Since the energiesfrom FP are not biased by the phaseless constraint andthus insensitive to the trial function used, (any) smallresidual errors can be attributed to the MP2-assistedCBS extrapolation scheme.The long equilibration times and inaccuracies encoun- tered in the above cases are manifestations of the factthat single-determinant trial functions obtained frommean-field calculations are generally not well-suited todescribe the open-shell systems involved in redox reac-tions. For instance in the IP of B, comparing FP and ph-AFQMC calculations for both the neutral and cationicspecies exposes the fact that the error in the IP stemsfrom inaccuracies in the computation of the total energyof the neutral B atom, which has a single unpaired elec-tron in the triply-degenerate p orbital manifold. The useof trial wavefunctions with proper symmetry propertieshas previously been shown to lead to improved accuracywithin the phaseless approximation,[84] so we now con-sider trial functions which are eigenfunctions of ˆ S foruse with the phaseless constraint. While employing a re-stricted open-shell HF trial affords no appreciable gainin accuracy, the use of a truncated CASSCF trial witha very modest active space size and a small number ofdeterminants is sufficient to eliminate most of the errorfor the IPs of both B and N.The accurate description of the EA of F proves to bemore demanding. It is possible that the phase prob-lem is particularly problematic for the case of F- dueto the small atomic radius which may lead to very strongelectronic correlations. Such sizable correlations man-ifest themselves in the algorithm as a large imaginarycomponent of the propagator, which in turn leads tothe elimination of crucial physical information during theprojection to the real axis if the trial function does notadequately provide the gauge information on the Slaterdeterminants of the ground state. In such cases, exci-tations into a large number of virtual orbitals will makesignificant contributions to the correlation energy, andtherefore a large active space is required to generate theCASSCF trial function so that the resulting phase pro-jections are sufficiently benign.In general, the implementation and efficacy of the cor-related sampling scheme presented in this work remainunaltered in the case of a multi-determinant trial func-tion since walker determinants are propagated in exactlythe same manner. We do find a reduction in the requiredequilibration times for all multi-determinant ph-AFQMCcalculations that we have performed, relative to single-determinant calculations of the same systems, renderingthe use of PES unnecessary. In addition, the use of multi-determinant trial functions results in drastically smallererror bars, even when more than 4x fewer walkers areused. Increasing the quality of the CASSCF trial func-tion is a promising way to systematically reduce the errorfrom the phaseless constraint, and we find that in all casesthe use of a reasonable number of determinants producesredox energies within the maximum experimental errorof 0.05 eV for the G2 Test Set. TABLE I: Experimental IPs and the deviations of various calculated results(theory - experiment) for atoms in the G2 Test Set in eV. QMC statistical errors inthe two right-most digits are shown in parenthesis. QMC calculations usingsingle-determinant trial functions (phaseless and FP) use 5040 walkers per repeat,while those using multi-determinant trial functions use 1056 walkers per repeat.
Atom Expt. ∆ph-HF/QMC ∆FP ∆ph-CAS/QMC ∆G2 ∆B3LYPB 8.2980 -0.156(10) * -0.0012(50) -0.0162(26) a * b * c * -0.0360(24) d e f -0.05 -0.21 * PES a QMC trial from CASSCF(5/4,8) with minimum CI coefficient of 0.034 (20/14determinants) b QMC trial from CASSCF(4/3,8) with minimum CI coefficient of 0.034 (14/7determinants) c QMC trial from CASSCF(5/4,8) with minimum CI coefficient of 0.01 (29/32determinants) d QMC trial from CASSCF(6/5,13) with minimum CI coefficient of 0.01 (67/62determinants) e QMC trial from CASSCF(7/6,8) with minimum CI coefficient of 0.033 (7/2 de-terminants) f QMC trial from CASSCF(8/7,16) with minimum CI coefficient of 0.0085(101/104 determinants)TABLE II: Experimental EAs and the deviations of various calculated results(theory - experiment) for atoms in the G2 Test Set in eV. QMC statistical errorsin the two right-most digits are shown in parenthesis. QMC calculations usingsingle-determinant trial functions (phaseless and FP) use 5040 walkers perrepeat, while those using multi-determinant trial functions use 1056 walkers perrepeat.
Atom Expt. ∆ph-HF/QMC ∆FP ∆ph-CAS/QMC ∆G2 ∆B3LYPB 0.2797 -0.0422(31) -0.0090(30) a * b c * -0.033(20) 0.0474(49) d -0.08 -0.13 * PES a QMC trial from CASSCF(3/4,8) with minimum CI coefficient of 0.01 (34/46determinants) b QMC trial from CASSCF(4/5,8) with minimum CI coefficient of 0.01 (30/53determinants) c QMC trial from CASSCF(6/7,8) with minimum CI coefficient of 0.01 (39/77determinants) d QMC trial from CASSCF(7/8,16) with minimum CI coefficient of 0.0075(145/243 determinants)
B. The Case of Methanol
Motivated by the arguments set forth in Section II.C,we now describe the details of a composite methodfor calculating the adiabatic
IPs, deprotonation free-energies, and hydrogen-dissociation energies of molecu-lar systems, and illustrate the accuracy of our approachon the case of methanol. We utilize a stepwise processconsisting of: 1) the fixed-geometry process calculatedwith correlated sampling-based ph-AFQMC, and 2) a geometry-relaxation step calculated within MP2. Op-timal correlation can be achieved in 1) since the sameset of basis functions is used, while in 2) we expect largeerror cancellation resulting from the fact that the initialand final ground-state geometries are typically very sim-ilar in redox, deprotonation, and hydrogen abstractionreactions. The calculations are performed in a triple-zeta basis with additional diffuse basis functions for allspecies involved in the deprotonation reaction. The re-sulting energy differences from these two steps are addedtogether, and the endpoints are extrapolated to the CBS0limit.For all reaction types, zero-point energies are calcu-lated at the level of HF/6-31G* and scaled by a factorof 0.899 to account for anharmonicity and the knownshortcomings of HF theory, following the G2 protocol.Deprotonation free-energy results incorporate the valueof -6.28 kcal/mol as the free-energy of a proton at 298K,[109] and we use the exact ground-state energy of H · (-0.5 Ha) to calculate the bond dissociation enthalpy ofthe O-H bond.Table III illustrates the accuracy of our correlated sam-pling protocol with respect to experimental results for allthree reaction types. The quality of our IP and deproto-nation free-energy results surpass that of the more costlyG2 method, while the O-H bond dissociation result hasan error of comparable magnitude. We note that thesingle-determinant trial wavefunction is sufficient to pro-duce a near-exact deprotonation free-energy, which we at-tribute to the fact that in this reaction type both the pro-tonated and deprotonated species are closed-shell. More-over, our correlated sampling-based ph-AFQMC calcula-tions, using an extremely modest number of walkers, pre-dicts the energy differences corresponding to all of thesereaction types to within chemical accuracy, which is notthe case for the G2 method. C. Basis Set Size, Number of Random Walkers,and CPU-time Reduction
1. Basis Set Size
The statistical noise of individual AFQMC runs is ex-pected to increase with the number of AFs, α , yet theprecise scaling is subtle. Each AF contributes additionalnoise. On the other hand, the magnitude of the con-tribution is moderated by a degree which depends onthe form of the interaction and the decomposition pro-cedure leading to Eq. (5). For example, the error baris seen to change little beyond a modest cutoff in plane-wave AFQMC.[112] In principle α grows as M , but our(rather conservative) truncation of the Cholesky decom-position via the cutoff mentioned in Section II.D resultsin α ∼ M . As a result, we expect the statistical errorin an uncorrelated calculation to increase with M beforesaturating toward the CBS limit. This could be prob-lematic given that large basis sets containing functionsassociated with high angular momenta are frequently re-quired to accurately describe the correlation energy,[99]and also given the fact that most interesting applicationsinvolve large systems.Here we show that the magnitude of the reduction instatistical error enabled by the use of correlated samplinggrows with M , such that the error is nearly independentof M . Fig. 6 shows the standard errors resulting fromcalculations of the IP of the K in the 6-31G*, 6-31+G*,and 6-311+G* basis sets (which consist of 23, 35, and45 basis functions, respectively); Fig. 7 illustrates the same effect for fixed-geometry deprotonation reactionsof water, methanol, and ethanol (58, 116, and 174 ba-sis functions). For both reaction types, while the errorsfrom the uncorrelated calculations increase significantlywith M , we find that those resulting from correlated sam-pling remained roughly constant. We anticipate that thisfinding will hold in general, and will be of crucial impor-tance in future applications of correlated sampling-basedAFQMC to larger molecules. FIG. 6: Comparison of the standard errors of the cumulativeaverages resulting from correlated and uncorrelated sampling incomputing the IP of the K atom in three different basis sets with5040 walkers per repeat. The inset zooms in on the lower region,for clarity.FIG. 7: Comparison of the errors resulting from correlated anduncorrelated sampling in computing the fixed-geometrydeprotonation energies of H O (M=58), CH OH (M=116), andC H OH (M=174) with 192 walkers per repeat.
2. Reduction in the Number of Random Walkers andCPU-Time
In an uncorrelated QMC calculation for a fixed lengthof propagation time, the resulting standard error can be1
TABLE III: Adiabatic Reaction Energies for Methanol.Experimental values and the deviations of the correlatedsampling-based ph-AFQMC and G2 results (theory - experiment) ineV. QMC statistical errors in the two right-most digits are shown inparenthesis. All QMC calculations use 192 walkers per repeat.
Reaction Expt. ∆ph-QMC ∆G2Ionization Potential 10.84 0.034(27) † -0.11Deprotonation Free-Energy 16.2695 -0.005(21) †† a O-H Bond Dissociation Energy 4.5359 0.039(14) † b † QMC trial from CASSCF(10/11, 10/11) with minimum CI coef-ficient of 0.02 (27/32 determinants) †† HF trial a Ref. [110] b Ref. [111] reduced by increasing the number of random walkers, N wlk , used in the MC evaluation of the HS integral inEq. (5). However, given that the required computationalexpense increases linearly with N wlk , using a brute-forceapproach that simply increases N wlk is less practical formany systems. For energy differences, correlated sam-pling provides a much cheaper alternative as it allows fora dramatic reduction in the N wlk required to achieve agiven statistical error.The errors of the cumulative averages of the IP of K inthe 6-31+G* basis (M=35) and the fixed-geometry de-protonation of methanol in the cc-pVTZ basis (M=116)are shown for different values of N wlk in Figs. 8 and9, respectively. For both reaction types, we make thefollowing observations: First, for a given N wlk , the stan-dard error is significantly lower in the correlated samplingcase. Second, the magnitude of this reduction is greaterfor smaller N wlk . Finally, while the standard errors ofthe uncorrelated runs increase as N wlk is reduced, in thecorrelated sampling runs the error is relatively insensitiveto N wlk .In light of these findings, we are in a position to un-derstand how a reduction in the standard error due tocorrelating the AFs translates into a reduction in CPU-time. Considering first the IP of K in the 6-31+G* basis,we choose a target standard error of 0.5 mHa on the QMCenergy (which defines the 99% confidence interval as thecumulative average of the energy ± i.e. within the 99% con-fidence interval) with a benchmark result obtained with5040 walkers. In fact, as few as 6 walkers produced thesame level of accuracy in some cases. The inset of Fig.8 shows that the statistical error falls below the targetat τ ∼
7. The total CPU-time required to propagate 11repeats for this length of imaginary-time is 41.2 minuteson a single 2.60 GHz Intel Xeon processor. Using thesame number of walkers without correlating the AFs, wefind that the target error is not reached even after 200
FIG. 8: Dependence of standard error on the number of randomwalkers per repeat for the IP of K in the 6-31+G* basis. Theinset highlights the errors of the uncorrelated run with 1056walkers and the correlated run with 24 walkers, compared withthe 0.5 mHa error target. Ha − . Using 1056 walkers gives rise to a standard errorthat falls below 0.5 mHa after 10 Ha − , as shown in theinset of Fig. 8, and a resulting QMC energy that is inagreement with the benchmark result. This calculationtakes 2262.6 minutes on a single processor, and we thusconclude that our correlated sampling approach reducesthe CPU-time by a factor of approximately 55.For the deprotonation of methanol we use a target errorof 1 mHa. As shown in the inset of Fig. 9, both the un-correlated run with 2400 walkers and the correlated runwith 96 walkers yield results that fall below our target er-ror at τ ∼
10. We use the fact that the total CPU-timeis proportional to the product of the number of walkersand the propagation time to estimate that correlatingthe AFs reduces the CPU-time by a factor of approxi-mately 25. Currently the total calculation, including all11 repeats, requires ∼
154 hours on a single CPU core.We perform a similar analysis for the dissociation of H · from the O-H bond of methanol. Due to the relativelylarge computational cost of using a CASSCF trial func-2 FIG. 9: Dependence of standard error on the number of randomwalkers per repeat for the fixed-geometry deprotonation ofmethanol in the cc-pVTZ basis. The inset highlights the errors ofthe uncorrelated run with 2400 walkers and the correlated runwith 96 walkers, compared with the 1 mHa error target. tion (with the same active space and truncation schemedescribed in Table III), we increase our target error to 2mHa. Figure 10 shows that the errors on the cumulativeaverages from both an uncorrelated run with 288 walkersand a correlated run with 12 walkers fall below our targeterror at τ ∼
6. Thus we deduce a speed-up factor of 24for this H abstraction reaction. On a single CPU corethis requires ∼
157 hours.
FIG. 10: Comparison of the standard errors resulting from theuse of correlated sampling with 12 walkers per repeat anduncorrelated sampling with 288 walkers per repeat to calculatethe energy difference associated with the fixed-geometry removalof H · from the O-H bond of methanol in the cc-pVTZ basis. The2 mHa error target is shown in black. As the previous sections have shown, the magnitude ofthe reduction in the statistical error, and consequentlythe CPU-time, as compared with an uncorrelated calcu-lation depends on the number of walkers used and the sizeof the basis set. In addition, one intuitively expects ourcorrelated sampling approach to work better for systems that more closely resemble each other. This is crudely thecase, yet the subtleties involved in such a claim warrantfurther discussion. Indeed, referring to the propagator in(11), while our correlated sampling method ensures thatthe ˆv operators and the AFs, x , are the same for boththe primary and secondary systems, the FBs, ¯x , as de-fined in (9) and the expectation values with respect tothe trial functions (cid:104) ˆv (cid:105) will in general be different. In thelimit that the primary and secondary systems are identi-cal, the entire propagator in (11) is identical for both sys-tems, and the statistical error in the energy difference willbe exactly and trivially zero as a result of perfect walker-pair correlation. Otherwise, the reduction in statisticalerror afforded by our correlated sampling approach be-comes less pronounced the more the trial wavefunctionsof the primary and secondary systems differ, since thethe FBs and (cid:104) ˆv (cid:105) are expectation values that depend ex-plicitly on the trial wavefunctions. It is encouraging tonote that despite the differences in trial functions corre-lating only the AFs yields such large speed-ups in CPU-time. While additionally correlating the (cid:104) ˆv (cid:105) , possibly byusing some combination of the trial functions for bothsystems, would compromise the accuracy of the phase-less approximations, future studies will explore optimalways to pair walkers such that the similarity in the FBsof walker pairs is maximized. It is encouraging that evenfor systems for which the electronic energies of the pri-mary and secondary systems differ by some 200 eV, as isthe case in the deprotonation of methanol, our correlatedsampling approach still yields dramatic efficiency gains. IV. CONCLUSIONS AND OUTLOOK
In this work we have devised a correlated samplingprotocol for the calculation of chemically relevant en-ergy differences within the exact and phaseless variantsof AFQMC. For molecules we utilize a two step strat-egy in which optimal walker-pair correlation is achievedin the ph-AFQMC description of the fixed-geometry pro-cess, while the geometry relaxation energy is calculatedwith the confines of MP2. Together with an MP2-assisted CBS extrapolation method we obtain calculatedIPs, EAs, deprotonation free-energies, and bond dissoci-ation energies that are in excellent agreement with ex-periments. Moreover, our correlated sampling approachyields large reductions in the statistical errors relativeto those obtained from uncorrelated approaches. In con-trast to uncorrelated AFQMC, where the error bars arefound to increase with system size and/or the numberof basis functions, correlating the AFs keeps the statis-tical error relatively constant as chemical complexity in-creases. In addition, our approach drastically reducesthe number of walkers required to achieve a given statis-tical error target, which translates into large reductionsin CPU-time.We utilize a “ghost” basis function strategy that en-ables the application of a correlated sampling-based ap-3proach to processes involving large energetic changes. In-deed, given that the correlated sampling scheme outlinedhere is successful for electron, proton, and H · transfer re-actions, we are optimistic that other chemical changesare within reach. In a future work we will systematicallyinvestigate chemical reactions which involve substantialchanges in geometry, including the addition/removal oflarger, more complex functional groups. Along the samelines we are optimistic about the savings that our corre-lated sampling approach may yield when basis functionsthat are independent of the nuclear coordinates, such asplane-waves, are used. We anticipate that the insensitiv-ity of the statistical error to basis set size that we observein this work may partially or totally offset the relativelylarge number of plane-waves typically required for con-vergence.Additional improvements in computational efficiencyare necessary before our correlated sampling-basedAFQMC approach can be routinely applied to the large,complex systems that we ultimately hope to investigate.In addition to incorporating orbital localization tech- niques and experimenting with more compact basis sets,we will pursue a variety of algorithmic ( e.g. the imple-mentation of mixed-precision and sparse linear algebraroutines) and hardware ( e.g. porting the code to runon GPUs) optimizations. Future targets of investigationinclude the IP of B F ,[113, 114] redox potentials andpKa’s of TM clusters and battery electrolytes,[48, 115–121] and the relative energies of low-lying states of inter-acting TM centers such as Fe-S complexes and the activesite of Photosystem II.[23, 42, 55] V. ACKNOWLEDGEMENTS
JS gratefully acknowledge Fengjie Ma, Mario Motta,Hao Shi, and James Dama for their generous help andguidance in implementing and understanding AFQMCmethods. DRR acknowledges funding from grant NSF-CHE-1464802. The work of SZ was supported by grantsde-sc0001303 and de-sc0008627. [1] K. A. Peterson, D. Feller, and D. A. Dixon, Theor.Chem. Acc. , 1 (2012).[2] R. A. Friesner, Proc. Natl. Acad. Sci. USA , 6648(2005).[3] S. Langhoff,
Quantum mechanical electronic structurecalculations with chemical accuracy , Vol. 13 (SpringerScience & Business Media, 2012).[4] C. D. Sherrill, J. Chem. Phys. , 110902 (2010).[5] J. N. Harvey, Annu. Rep. Prog. Chem., Sect. C: Phys.Chem. , 203 (2006).[6] A. Szabo and N. S. Ostlund,
Modern quantum chem-istry: introduction to advanced electronic structure the-ory (Courier Corporation, 1989).[7] G. D. Purvis III and R. J. Bartlett, J. Chem. Phys. ,1910 (1982).[8] R. J. Bartlett and M. Musia(cid:32)l, Rev. Mod. Phys. , 291(2007).[9] W. Purwanto, S. Zhang, and H. Krakauer, J. Chem.Phys. , 064302 (2015).[10] A. Dutta and C. D. Sherrill, J. Chem. Phys. , 1610(2003).[11] B. O. Roos, Advances in Chemical Physics: Ab InitioMethods in Quantum Chemistry Part 2, Volume 69 ,399 (1987).[12] J. Olsen, Int. J. Quantum Chem. , 3267 (2011).[13] M. W. Schmidt and M. S. Gordon, Annu. Rev. Phys.Chem. , 233 (1998).[14] K. Andersson, P. A. Malmqvist, B. O. Roos, A. J.Sadlej, and K. Wolinski, J. Phys. Chem. , 5483(1990).[15] K. Andersson, P.-˚A. Malmqvist, and B. O. Roos, J.Chem. Phys. , 1218 (1992).[16] G. H. Booth, A. J. W. Thom, andA. Alavi, J. Chem. Phys. , 054106 (2009),http://aip.scitation.org/doi/pdf/10.1063/1.3193710. [17] G. H. Booth, D. Cleland, A. J. W. Thom, andA. Alavi, J. Chem. Phys. , 084104 (2011),http://dx.doi.org/10.1063/1.3624383.[18] D. Cleland, G. H. Booth, and A. Alavi, J. Chem. Phys. , 024112 (2011).[19] D. Cleland, G. H. Booth, C. Overy, and A. Alavi, J.Chem. Theory Comput. , 4138 (2012).[20] K. Pierloot, Int. J. Quantum Chem. , 3291 (2011).[21] R. E. Thomas, G. H. Booth, and A. Alavi, Phys. Rev.Lett. , 033001 (2015).[22] G. K.-L. Chan and S. Sharma, Annu. Rev. Phys. Chem. , 465 (2011).[23] S. Sharma, K. Sivalingam, F. Neese, and G. K.-L. Chan,Nat. Chem. , 927 (2014).[24] Y. Kurashige, G. K.-L. Chan, and T. Yanai, Nat. Chem. , 660 (2013).[25] W. Kohn, A. D. Becke, and R. G. Parr, J. Phys. Chem. , 12974 (1996).[26] W. A. Lester Jr and B. L. Hammond, Annu. Rev. Phys.Chem. , 283 (1990).[27] B. L. Hammond, W. A. Lester Jr, and P. J. Reynolds, Monte Carlo methods in ab initio quantum chemistry ,Vol. 1 (World Scientific, 1994).[28] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Ra-jagopal, Rev. Mod. Phys. , 33 (2001).[29] S. Zhang, Emergent Phenomena in Correlated Matter (2013).[30] S. Zhang and H. Krakauer, Phys. Rev. Lett. , 136401(2003).[31] W. Al-Saidi, H. Krakauer, and S. Zhang, J. Chem.Phys. , 154110 (2006).[32] Y. Virgus, W. Purwanto, H. Krakauer, and S. Zhang,Phys. Rev. Lett. , 175502 (2014).[33] W. A. Al-Saidi, H. Krakauer, and S. Zhang, Phys. Rev.B , 075103 (2006).[34] C. Umrigar, Int. J. Quantum Chem. , 217 (1989). [35] C. Traynor and J. B. Anderson, Chem. Phys. Lett. ,389 (1988).[36] B. H. Wells, Chem. Phys. Lett. , 89 (1985).[37] J. Vrbik, D. A. Legare, and S. M. Rothstein, J. Chem.Phys. , 1221 (1990).[38] C. Filippi and C. J. Umrigar, Phys. Rev. B , R16291(2000).[39] Y. Kwon, D. M. Ceperley, and R. M. Martin, Phys.Rev. B , 1684 (1994).[40] D. R. Garmer, J. Comput. Chem. , 176 (1989).[41] D. L. Nelson, A. L. Lehninger, and M. M. Cox, Lehninger principles of biochemistry (Macmillan, 2008).[42] P. E. Siegbahn and M. R. Blomberg, Chem. Rev. ,421 (2000).[43] B. Meunier, S. P. De Visser, and S. Shaik, Chem. Rev. , 3947 (2004).[44] M. R. Palacin, Chemical Society Reviews , 2565(2009).[45] M. Gratzel, Energy resources through photochemistryand catalysis (Elsevier, 2012).[46] C. Costentin, M. Robert, and J.-M. Sav´eant, ChemicalSociety Reviews , 2423 (2013).[47] E. Alexov, E. L. Mehler, N. Baker, A. M Baptista,Y. Huang, F. Milletti, J. Erik Nielsen, D. Farrell,T. Carstensen, M. H. Olsson, J. K. Shen, J. Warwicker,S. Williams, and M. J. Word, Proteins: Struct., Funct.,Bioinf. , 3260 (2011).[48] S. V. Jerome, T. F. Hughes, and R. A. Friesner, J.Phys. Chem. B , 8008 (2014).[49] A. D. Bochevarov, M. A. Watson, J. R. Green-wood, and D. M. Philipp, J. Chem. The-ory Comput. , 6001 (2016), pMID: 27951674,http://dx.doi.org/10.1021/acs.jctc.6b00805.[50] H. A. Rajapakse, P. G. Nantermet, H. G. Selnick,J. C. Barrow, G. B. McGaughey, S. Munshi, S. R.Lindsley, M. B. Young, P. L. Ngo, M. K. Holloway,M.-T. Laid, A. S. Espesethd, X.-P. Shid, D. Colus-sid, B. Pietrakd, M.-C. Crouthameld, K. Tugushevad,Q. Huangd, M. Xud, A. J. Simond, L. Kuoc, D. J.Hazudad, S. Grahama, and J. P. Vaccaa, Bioorg. Med.Chem. Lett. , 1885 (2010).[51] D. Sprous, R. Palmer, J. Swanson, and M. Lawless,Current topics in medicinal chemistry , 619 (2010).[52] J. Cheng and M. Sprik, J. Chem. Theory Comput. ,880 (2010).[53] D. R. Gallus, R. Wagner, S. Wiemers-Meyer, M. Win-ter, and I. Cekic-Laskovic, Electrochim. Acta , 410(2015).[54] V. S. Bryantsev, Chem. Phys. Lett. , 42 (2013).[55] W. Ames, D. A. Pantazis, V. Krewald, N. Cox,J. Messinger, W. Lubitz, and F. Neese, J. Am. Chem.Soc. , 19743 (2011).[56] J. El Yazal, F. G. Prendergast, D. E. Shaw, and Y.-P.Pang, J. Am. Chem. Soc. , 11411 (2000).[57] J. Chen, Y.-F. Li, P. Sit, and A. Selloni, J. Am. Chem.Soc. , 18774 (2013).[58] G. A. Olah and A. Molnar, Hydrocarbon chemistry (John Wiley & Sons, 2003).[59] H. Wang and M. Frenklach, J. Phys. Chem. , 11465(1994).[60] M. Page and D. W. Brenner, J. Am. Chem. Soc. ,3270 (1991).[61] G. W. Burton and K. U. Ingold, Ann. N. Y. Acad. Sci. , 7 (1989). [62] J. M. Mayer, Acc. Chem. Res. , 441 (1998).[63] M. R. Blomberg, P. E. Siegbahn, S. Styring, G. T. Bab-cock, B. ˚Akermark, and P. Korall, J. Am. Chem. Soc. , 8285 (1997).[64] R. K. Grasselli, Catal. Today , 141 (1999).[65] B. B. Snider, Chem. Rev. , 339 (1996).[66] H. Basch and S. Hoz, J. Phys. Chem. A , 4416(1997).[67] M. L. Coote, J. Phys. Chem. A , 3865 (2004).[68] J. Pu and D. G. Truhlar, J. Phys. Chem. A , 773(2005).[69] E. F. Carvalho, A. N. Barauna, F. B. Machado, andO. Roberto-Neto, Int. J. Quantum Chem. , 2476(2008).[70] E. Carvalho, A. N. Barauna, F. B. Machado, andO. Roberto-Neto, Chem. Phys. Lett. , 33 (2008).[71] F. Fracchia, C. Filippi, and C. Amovilli, J. Chem. The-ory Comput. , 3453 (2013).[72] Y. Kanai and N. Takeuchi, J. Chem. Phys. , 214708(2009).[73] A. Kollias, O. Couronne, and W. Lester Jr, J. Chem.Phys. , 1357 (2004).[74] L. A. Curtiss, P. C. Redfern, K. Raghavachari, andJ. A. Pople, J. Chem. Phys. , 42 (1998).[75] W. Purwanto, H. Krakauer, Y. Virgus, and S. Zhang,J. Chem. Phys. , 164105 (2011).[76] H. F. Trotter, Proc. Am. Math. Soc. , 545 (1959).[77] M. Suzuki, Commun. Math. Phys. , 183 (1976).[78] R. Stratonovich, Dokl. Akad. Nauk SSSR , 1097(1957).[79] J. Hubbard, Phys. Rev. Lett. , 77 (1959).[80] D. R. Hamann and S. B. Fahy, Phys. Rev. B , 11352(1990).[81] H. Nguyen, Quantum Monte Carlo Calculation of theImaginary-Time Greens Function in the Hubbard Model ,Undergraduate thesis, Reed College (2014).[82] B. M. Rubenstein, S. Zhang, and D. R. Reichman,Phys. Rev. A , 053606 (2012).[83] P. J. Reynolds, D. M. Ceperley, B. J. Alder, and W. A.Lester Jr, J. Chem. Phys. , 5593 (1982).[84] H. Shi and S. Zhang, Phys. Rev. B , 125132 (2013).[85] S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev.B , 7464 (1997).[86] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R.White, D. J. Scalapino, and R. L. Sugar, Phys. Rev. B , 9301 (1990).[87] M. Troyer and U.-J. Wiese, Phys. Rev. Lett. , 170201(2005).[88] M. Motta, D. E. Galli, S. Moroni, and E. Vitali, J.Chem. Phys. , 024107 (2014).[89] W. Purwanto and S. Zhang, Phys. Rev. E , 056702(2004).[90] W. Al-Saidi, S. Zhang, and H. Krakauer, J. Chem.Phys. , 224101 (2006).[91] H. Flyvbjerg and H. G. Petersen, J. Chem. Phys. ,461 (1989).[92] H. Nguyen, H. Shi, J. Xu, and S. Zhang, Comput. Phys.Commun. , 3344 (2014).[93] M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P.Straatsma, H. J. Van Dam, D. Wang, J. Nieplocha,E. Apra, T. L. Windus, and W. A. de Jong, Comput.Phys. Commun. , 1477 (2010).[94] W. Purwanto, W. Al-Saidi, H. Krakauer, and S. Zhang,J. Chem. Phys. , 114309 (2008). [95] “PySCF,” https://github.com/sunqm/pyscf, 2014.[96] S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh,J. E. Gubernatis, and R. T. Scalettar, Phys. Rev. B , 506 (1989).[97] W. Purwanto, H. Krakauer, and S. Zhang, Phys. Rev.B , 214116 (2009).[98] M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T.Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Mat-sunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis,and J. A. Montgomery Jr, J. Comput. Chem. , 1347(1993).[99] T. H. Dunning Jr, J. Chem. Phys. , 1007 (1989).[100] R. A. Kendall, T. H. Dunning Jr, and R. J. Harrison,J. Chem. Phys. , 6796 (1992).[101] R. J. Johnson III, “NIST Computational Chem-istry Comparison and Benchmark Database,”http://cccbdb.nist.gov/, NIST Standard ReferenceDatabase Number 101, Release 17b, September 2015.[102] Q. Wang, X. Zhang, Y. Zhang, and Q. Yi, in Proceed-ings of the International Conference on High Perfor-mance Computing, Networking, Storage and Analysis ,SC ’13 (ACM, New York, NY, USA, 2013) pp. 25:1–25:12.[103] Z. Xianyi, W. Qian, and Z. Yunquan, in (2012) pp. 684–691.[104] R. B. Sidje, ACM Transactions on Mathematical Soft-ware (TOMS) , 130 (1998).[105] M. Mascagni and A. Srinivasan, ACM Transactions onMathematical Software (TOMS) , 436 (2000). [106] L. A. Curtiss, K. Raghavachari, G. W. Trucks, andJ. A. Pople, J. Chem. Phys. , 7221 (1991).[107] A. D. Becke, J. Chem. Phys. , 5648 (1993).[108] C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B ,785 (1988).[109] G. Tawa, I. Topol, S. Burt, R. Caldwell, and A. Rashin,J. Chem. Phys. , 4852 (1998).[110] G. N. Merrill and S. R. Kass, J. Phys. Chem. , 17465(1996).[111] L. A. Curtiss, L. D. Kock, and J. A. Pople, J. Chem.Phys. , 4040 (1991).[112] M. Suewattana, W. Purwanto, S. Zhang, H. Krakauer,and E. J. Walter, Phys. Rev. B , 245123 (2007).[113] B. Chan, A. J. Trevitt, S. J. Blanksby, and L. Radom,J. Phys. Chem. A , 9214 (2012).[114] Z.-H. Li and K.-N. Fan, J. Phys. Chem. A , 6659(2002).[115] D. Coskun, S. V. Jerome, and R. A. Friesner, J. Chem.Theory Comput. , 1121 (2016).[116] L. E. Roy, E. Jakubikova, M. G. Guthrie, and E. R.Batista, J. Phys. Chem. A , 6745 (2009).[117] S. J. Konezny, M. D. Doherty, O. R. Luca, R. H. Crab-tree, G. L. Soloveichik, and V. S. Batista, J. Phys.Chem. C , 6349 (2012).[118] M.-H. Baik and R. A. Friesner, J. Phys. Chem. A ,7407 (2002).[119] O. Borodin, W. Behl, and T. R. Jow, J. Phys. Chem.C , 8661 (2013).[120] N. Shao, X.-G. Sun, S. Dai, and D.-e. Jiang, J. Phys.Chem. B , 3235 (2012).[121] X. Qu, A. Jain, N. N. Rajput, L. Cheng, Y. Zhang, S. P.Ong, M. Brafman, E. Maginn, L. A. Curtiss, and K. A.Persson, Comput. Mater. Sci.103