Chemical Transformations Approaching Chemical Accuracy via Correlated Sampling in Auxiliary-Field Quantum Monte Carlo

TL;DR

This paper introduces a correlated sampling method for Auxiliary-Field Quantum Monte Carlo that significantly improves efficiency in calculating energy differences, enabling accurate chemical property predictions with reduced computational cost.

Contribution

The authors develop a correlated sampling approach for AFQMC that enhances efficiency and accuracy in computing energy differences relevant to chemical transformations.

Findings

Correlated sampling accelerates AFQMC calculations by reducing statistical errors.

The method achieves chemical accuracy in ionization potentials, free energies, and bond dissociation energies.

Speed-up factors of 24 to 55 times are observed for various chemical properties.

Abstract

The exact and phaseless variants of Auxiliary-Field Quantum Monte Carlo (AFQMC) have been shown to be capable of producing accurate ground-state energies for a wide variety of systems including those which exhibit substantial electron correlation effects. The computational cost of performing these calculations has to date been relatively high, impeding many important applications of these approaches. Here we present a correlated sampling methodology for AFQMC which relies on error cancellation to dramatically accelerate the calculation of energy differences of relevance to chemical transformations. In particular, we show that our correlated sampling-based AFQMC approach is capable of calculating redox properties, deprotonation free-energies, and hydrogen abstraction energies in an efficient manner without sacrificing accuracy. We validate the computational protocol by calculating the ionization potentials and electron affinities of the atoms contained 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, to compute the ionization potential, deprotonation free-energy, and the O-H bond disocciation energy of methanol, all to within chemical accuracy. We show that the efficiency of correlated sampling relative to uncorrelated calculations increases with system and basis set size, and that correlated sampling 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 the K 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 the door to the accurate description of chemical processes in complex systems.

Figures

Tables

Table 1. Table 1: Experimental IPs and the deviations of various calculated results (theory - experiment) for atoms in the G2 Test Set in eV. QMC statistical errors in the two right-most digits are shown in parenthesis. QMC calculations using single-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.	$\Delta$ph-HF/QMC	$\Delta$FP	$\Delta$ph-CAS/QMC	$\Delta$G2	$\Delta$B3LYP
B	8.2980	-0.156(10)*	-0.0012(50)	-0.0162(26)a	0.10	-0.44
C	11.2603	0.0342(90)*		0.0045(43)b	0.08	-0.29
N	14.5341	0.1214(61)*	0.0050(84)	0.0100(38)c	0.06	-0.14
O	13.6181	-0.0830(22)*		-0.0360(24)d	0.08	-0.55
F	17.4228	0.0010(35)		0.0015(46)e	0.03	-0.34
Ne	21.5645	0.0775(51)		0.0159(34)f	-0.05	-0.21

Table 2. Table 2: Experimental EAs and the deviations of various calculated results (theory - experiment) for atoms in the G2 Test Set in eV. QMC statistical errors in the two right-most digits are shown in parenthesis. QMC calculations using single-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.	$\Delta$ph-HF/QMC	$\Delta$FP	$\Delta$ph-CAS/QMC	$\Delta$G2	$\Delta$B3LYP
B	0.2797	-0.0422(31)		-0.0090(30)a	0.09	-0.18
C	1.2621	0.0762(84)*		0.0031(43)b	0.07	-0.11
O	1.4620	0.0505(86)		0.0327(53)c	0.06	-0.22
F	3.4013	0.222(18)*	-0.033(20)	0.0474(49)d	-0.08	-0.13

Table 3. Table 3: Adiabatic Reaction Energies for Methanol. Experimental values and the deviations of the correlated sampling-based ph-AFQMC and G2 results (theory - experiment) in eV. QMC statistical errors in the two right-most digits are shown in parenthesis. All QMC calculations use 192 walkers per repeat.

Reaction	Expt.	$\Delta$ph-QMC	$\Delta$G2
Ionization Potential	10.84	0.034(27)$†$	-0.11
Deprotonation Free-Energy	16.2695	-0.005(21)$††$	0.0484a
O-H Bond Dissociation Energy	4.5359	0.039(14)$†$	0.0166b

Equations

$$\hat{H}=\sum _{ij}^M{T}_{ij}\sum _{\sigma }{c}_{i\sigma }^{†}{c}_{j\sigma }+\frac{1}{2}\sum _{ijkl}^M{V}_{ijkl}\sum _{\sigma ,\tau }{c}_{i\sigma }^{†}{c}_{j\tau }^{†}{c}_{l\tau }{c}_{k\sigma },$$

$${\hat{H}}_2=-\frac{1}{2}\sum _{\alpha }({\hat{v}}_{\alpha }-⟨{\hat{v}}_{\alpha }⟩{)}^2.$$

$$e^{-\Delta \tau \hat{H}}=e^{-\Delta \tau {\hat{H}}_1/2}{e}^{-\Delta \tau {\hat{H}}_2}{e}^{-\Delta \tau {\hat{H}}_1/2}+O(\Delta {\tau }^3).$$

$$e^{\frac{1}{2}\Delta\tau(\hat{v}_{\alpha}-\langle\hat{v}_{\alpha}\rangle)^{2}}=\int_{-\infty}^{\infty}dx_{\alpha}\bigg{(}\frac{e^{-\frac{1}{2}x_{\alpha}^{2}}}{\sqrt{2\pi}}\bigg{)}e^{\sqrt{\Delta\tau}x_{\alpha}(\hat{v}_{\alpha}-\langle\hat{v}_{\alpha}\rangle)},$$

$$e^{-\Delta \tau \hat{H}}=\int d\mathbf{x}P(\mathbf{x})\hat{B}(\mathbf{x}),$$

$$|\Phi ⟩=\sum _k\frac{w_k|{\varphi }_k⟩}{⟨{\varphi }_T|{\varphi }_k⟩},$$

$$\frac{⟨{\varphi }_T|\hat{H}|\Phi ⟩}{⟨{\varphi }_T|\Phi ⟩}=\frac{{\sum }_k{w}_k{E}_L({\varphi }_k)}{{\sum }_k{w}_k},$$

$$\frac{\partial}{\partial x_{\alpha}}\bigg{[}\frac{\langle\phi_{T}|e^{\sqrt{\Delta\tau}(x_{\alpha}-\bar{x}_{\alpha})(\hat{v}_{\alpha}-\langle\hat{v}_{\alpha}\rangle)}|\phi_{k}\rangle}{\langle\phi_{T}|\phi_{k}\rangle}e^{-\bar{x}_{\alpha}^{2}/2+x_{\alpha}\bar{x}}\bigg{]}_{x_{\alpha}=0}=0.$$

$$\bar{x}_{\alpha}=-\sqrt{\Delta\tau}\big{[}\bar{v}_{\alpha}-\langle\hat{v}_{\alpha}\rangle\big{]},$$

$$E_L(\varphi )\equiv Re\{\frac{⟨{\varphi }_T|\hat{H}|\varphi ⟩}{⟨{\varphi }_T|\varphi ⟩}\}.$$

$$\hat{B}(\mathbf{x}-\bar{\mathbf{x}})=e^{-\Delta \tau {\hat{H}}_1/2}{e}^{\sqrt{\Delta \tau }(\mathbf{x}-\bar{\mathbf{x}})(\hat{\mathbf{v}}-⟨\hat{\mathbf{v}}⟩)}{e}^{-\Delta \tau {\hat{H}}_1/2},$$

$$V_{ijkl}=\int d{\mathbf{r}}_{1}d{\mathbf{r}}_2{\varphi }_i({\mathbf{r}}_1){\varphi }_j({\mathbf{r}}_2)\frac{1}{r_{12}}{\varphi }_k({\mathbf{r}}_1){\varphi }_l({\mathbf{r}}_2).$$

Full text

I Introduction

The pursuit of chemical accuracy (as defined by errors not exceeding 1 kcal/mol) in the ab initio computation of the energetic properties of generic many-electron systems is a long-standing goal that has yet to be reached.Peterson et al. (2012); Friesner (2005); Langhoff (2012); Sherrill (2010); Harvey (2006) Traditional methods such as full configuration interaction (FCI)Szabo and Ostlund (1989) and single-reference coupled cluster with single, double, and perturbative triple excitations (CCSD(T))Purvis III and Bartlett (1982); Bartlett and Musiał (2007) scale exponentially and with the seventh power of the system size, respectively, and the latter can break down for systems which exhibit sufficiently strong electron correlation.Purwanto et al. (2015); Dutta and Sherrill (2003) Multi-reference methods such as CASSCF,Roos (1987); Olsen (2011); Schmidt and Gordon (1998) often supplemented with second order perturbation theory,Andersson et al. (1990, 1992) and FCI-QMCBooth et al. (2009, 2011); Cleland et al. (2011, 2012) have been shown to yield benchmark-quality results even for strongly correlated electronic systems, an import class of problems that includes many transition metal (TM)-containing systems.Pierloot (2011); Thomas et al. (2015) However, while a judicious choice of the active space can make such calculations tractable for small systems, exponential scaling prohibits the use of these methods in studying most realistic systems of interest in biology, materials science, and chemical catalysis. As a result, with the exception of isolated studies using more accurate, specialized methods,Chan and Sharma (2011); Sharma et al. (2014); Kurashige et al. (2013) these systems can only be investigated feasibly with less accurate but more economical approaches such as Density Functional Theory (DFT).Kohn et al. (1996) Given the plethora of investigations using this tool to generate hypotheses concerning reaction mechanisms, equilibrium structures, etc., the importance of having a systematically improvable ab initio method which is both accurate and feasible becomes readily apparent.

Among the class of Quantum Monte Carlo (QMC) methods used in electronic structure theory,Lester Jr and Hammond (1990); Hammond et al. (1994); Foulkes et al. (2001) Auxiliary-Field QMC, and in particular its “phaseless” variant (ph-AFQMC)Zhang (2013); Zhang and Krakauer (2003) which controls the sign problem at the cost of introducing a bias which in principle can be systematically reduced, has produced state-of-the-art, benchmark-quality results for systems such as first- and second- row $d$ elements,Al-Saidi et al. (2006a) the chromium dimer,Purwanto et al. (2015) cobalt adatoms on graphene,Virgus et al. (2014) and various transition metal oxides,Al-Saidi et al. (2006b) while exhibiting low-polynomial scaling ($M^{4}$ for Gaussian basis sets). However, the widespread use of ph-AFQMC in quantum chemistry has not yet taken place. One major reason is undoubtedly its relatively high computational cost, as the favorable scaling is masked by a large prefactor.

In this work we present an approach to greatly reducing this prefactor which involves the use of correlated sampling for a particular class of important processes. The general idea is that for sufficiently similar systems, energy differences are expected to converge more rapidly, i.e. with smaller error bars, than total energies when the errors or statistical fluctuations in the calculations are biased in the same direction. Indeed, error cancellation is largely responsible for the success of many approximate methods such as DFT in the computation of energy differences. Correlated sampling has previously been adapted to reduce the statistical errors in QMC approaches via the use of the same set of configurations sampled for both the primary and secondary systems.Lester Jr and Hammond (1990) This technique is often referred to in the literature as “Differential” QMC, and the details of its implementation can vary depending on the type of QMC being used. The potential energy curves of H2 and BH have been calculated using correlated sampling with Variational QMC,Umrigar (1989) and similarly that of the $H_{3}$ cation with Differential Green’s Function QMC.Traynor and Anderson (1988) The latter method has been used to compute the dipole moment of LiH, Wells (1985) and to calculate infinitesimal energy differences from which forces and various polarizabilities have been obtained.Vrbik et al. (1990) This idea has also been extended to Diffusion Monte Carlo, which has been used to compute forces and potential energy surfaces for the first row diatomics.Filippi and Umrigar (2000) Correlated sampling has also been used to calculate energy differences between ground and excited states of the same Hamiltonian, as illustrated by a Variational QMC study of particle-hole excitations in the two-dimensional electron gas.Kwon et al. (1994) In addition, the concept has been extended to enable concerted propagation of a system with different time steps, in order to extrapolate the Trotter error in Differential Diffusion QMC.Garmer (1989)

Correlated sampling is, in fact, also well-suited to model the energetics of myriad chemical reactions, since only energy differences, as opposed to total energies, are relevant. In this paper we present a novel correlated sampling-based AFQMC approach, and show that it is capable of computing ground-state energy differences corresponding to redox, deprotonation, and hydrogen abstraction reactions to a given statistical error in a fraction of the time previously required, without any loss of accuracy. Redox reactions (often involving TMs) abound, for example, in metabolic and photosynthetic processes,Nelson et al. (2008); Siegbahn and Blomberg (2000); Meunier et al. (2004) battery chemistry,Palacin (2009) and catalysis (e.g. CO2 and O2 reduction).Gratzel (2012); Costentin et al. (2013) A reliable ab initio method to calculate deprotonation free energies would provide an improvement upon existing computational approaches to determining pKa’s and protonation states,Alexov et al. (2011); Jerome et al. (2014); Bochevarov et al. (2016) which would have significant ramifications for drug discovery,Rajapakse et al. (2010); Sprous et al. (2010) materials science,Cheng and Sprik (2010); Gallus et al. (2015); Bryantsev (2013) and the structural determination of biological complexes.Ames et al. (2011); El Yazal et al. (2000) In the context of chemical catalysis, proton removal is known to be the rate-limiting step in many important reactions, e.g. oxygen reduction on the surface of TiO2.Chen et al. (2013) Hydrogen abstraction reactions are ubiquitous and play a major role in combustion and the oxidation of hydrocarbons,Olah and Molnar (2003); Wang and Frenklach (1994) diamond growth via chemical vapor deposition,Page and Brenner (1991) biochemical processes involving e.g. the antioxidant vitamin E Burton and Ingold (1989) and various metalloenzymes,Mayer (1998); Blomberg et al. (1997) industrial processes,Grasselli (1999) and organic synthesis.Snider (1996) AFQMC would be a useful benchmark for previous ab initio studiesBlomberg et al. (1997); Basch and Hoz (1997); Coote (2004); Pu and Truhlar (2005); Carvalho et al. (2008a, b); Fracchia et al. (2013); Kanai and Takeuchi (2009); Kollias et al. (2004) in predicting thermodynamic properties of this difficult class of chemical reactions. Thus the class of applications we target is large and important. We highlight the fact that the cost of our correlated sampling approach, relative to the uncorrelated method, decreases with increasing system and basis set size, opening the door to the treatment of realistic large chemical systems with correlated sampling-based AFQMC in the near future.

This work is organized as follows: Section II.A will review the exact and phaseless variants of AFQMC, while Section II.B will present our correlated sampling approach. We justify this approach for modeling molecular systems in Section II.C, and Section II.D will disclose further computational details. Section III.A will present calculations of the ionization potentials (IPs) and electron affinities (EAs) of the 1st row atoms in the G2 Ion Test Set,Curtiss et al. (1998) while Section III.B will consider adiabatic molecular properties, taking the IP, deprotonation energy, and O-H bond dissociation energy of methanol as examples. The efficacy of correlated sampling as a function of both basis set size and number of random walkers will be explored in Sections III.C. The latter subsection will provide an assessment of the reduction in CPU-time afforded by the use of correlated sampling. In Section IV we conclude, emphasizing opportunities for further gains in computational efficiency that may be possible and future targets of investigation.

II Methods

II.1 Overview of AFQMC - Exact Method and the Phaseless Constraint

The excited states of a many-body state $|\Phi\rangle$ can be projected out via imaginary time propagation, i.e. $\lim_{N\rightarrow\infty}(e^{\Delta\tau(E_{0}-\hat{H})})^{N}|\Phi\rangle\rightarrow|\Phi_{0}\rangle$, with $\langle\Phi_{0}|\Phi\rangle\neq 0$ and the general electronic Hamiltonian

[TABLE]

where $M$ is the size of the orthonormal one-particle basis, and $c_{i\sigma}^{\dagger}$ and $c_{i\sigma}$ are the second-quantized fermionic creation and annihilation operators with particle and spin labels. The two-body matrix elements, $V_{ijkl}$, can be expressed in terms of Cholesky vectors as $V_{ijkl}=\sum_{\alpha}L_{ik}^{\alpha}L_{jl}^{\alpha}$.Purwanto et al. (2011) Defining the one-body operator $\hat{v}_{\alpha}\equiv i\sum_{ik}L_{ik}^{\alpha}\sum_{\sigma}c_{i\sigma}^{\dagger}c_{k\sigma}$, and subtracting the expectation value with respect to the trial wavefunction $\langle\hat{v}_{\alpha}\rangle$ from $\hat{v}_{\alpha}$, the Hamiltonian can be written as the sum of all one-body operators, $\hat{H}_{1}$, plus the following two-body operator

[TABLE]

Use of the Trotter-Suzuki decompositionTrotter (1959); Suzuki (1976) gives

[TABLE]

The exponential terms involving $\hat{H}_{2}$ may be decomposed using a Hubbard-Stratonovich (HS) transformationStratonovich (1957); Hubbard (1959)

[TABLE]

which expresses the exponential of a two-body operator as the exponential of a one-body operator integrated over auxiliary-fields (AFs). This transformation allows for practical propagation in terms of the Thouless theorem, which states that the application of an exponential of a one-body operator on a Slater determinant produces another Slater determinant,Hamann and Fahy (1990) which can be implemented via a simple matrix multiplication.Nguyen (2014); Rubenstein et al. (2012)

The propagator (3) now takes the form of a multi-dimensional integral

[TABLE]

where $\mathbf{x}=(x_{1},x_{2},\ldots,x_{\alpha})$, $P(\mathbf{x})$ is a normal distribution with unit variance, and $\hat{B}(\mathbf{x})=e^{-\Delta\tau\hat{H}_{1}}e^{\sqrt{\Delta\tau}\mathbf{x}\cdot(\mathbf{\hat{v}}-\mathbf{\langle\hat{v}\rangle})}e^{-\Delta\tau\hat{H}_{1}}$.

The integral in (5) may be approximated using a Monte Carlo scheme, with walkers whose propagation in the space of Slater determinants is guided by the complex importance function $\langle\phi_{T}|\phi\rangle$ which is proportional to the walker weights. The representation of the total wavefunction is thus a weighted sum over walker determinants

[TABLE]

yielding essentially a multi-reference description. The energy is calculated at intervals using the mixed-estimator

[TABLE]

where the “local energy” is given by $E_{L}(\phi_{k})=\frac{\langle\phi_{T}|\hat{H}|\phi_{k}\rangle}{\langle\phi_{T}|\phi_{k}\rangle}$, in analogy to Diffusion QMC.Reynolds et al. (1982)

The method as described above is called the “Free Projection” (FP) approach.Shi and Zhang (2013) While it is formally exact and can yield excellent results for small system sizes, this method suffers from the “phase problem,” which is a generalization of the Fermionic “sign problem” to the complex plane.Zhang et al. (1997); Zhang (2013); Loh et al. (1990); Troyer and Wiese (2005) For the standard Coulomb interaction the $\hat{v}_{\alpha}$ operators are purely imaginary, and each application of $e^{\sqrt{\Delta\tau}\mathbf{x}\cdot(\mathbf{\hat{v}}-\mathbf{\langle\hat{v}\rangle})}$ can be thought of as a rotation of the Slater determinant $|\phi\rangle$, causing an evolution of the overlap $\langle\phi_{T}|\phi\rangle$ in the complex plane. Over the course of the random walk, a determinant accumulates a phase $e^{i\theta}$, and the infinitely many possible values of $\theta\ \epsilon\ [0,2\pi)$ result in the possibility of infinitely many indistinguishable determinants. Furthermore, over the course of the propagation, walkers will populate the origin where $\langle\phi_{T}|\phi\rangle=0$, and subsequent propagation yields only noise from signal cancellation effects and divergences in the weights and local energies.

The ph-AFQMC employs a multifaceted strategy to control this problem. The weights are initialized to a positive real constant, and after each propagation step are projected back onto the real axis, i.e. the rotated weights are multiplied by $max\{0,cos(\Delta\theta)\}$, where we have defined the phase of the overlap ratio $\Delta\theta\equiv Im\{\ln\frac{\langle\phi_{T}|\phi^{(\tau+1)}\rangle}{\langle\phi_{T}|\phi^{(\tau)}\rangle}\}$.Motta et al. (2014) For this phase projection to work, the AFs are shifted by a force bias (FB) $\mathbf{\bar{x}}$,Zhang and Krakauer (2003); Purwanto and Zhang (2004) the optimal choice of which is obtained by minimizing the fluctuations of the weights with respect to the AFs at their average value:

[TABLE]

This is, in essence, a stationary-phase approximation. Expanding the expression inside the brackets to $O(\sqrt{\Delta\tau})$ and taking the derivative gives

[TABLE]

where $\bar{v}_{\alpha}\equiv\frac{\langle\phi_{T}|\hat{v}_{\alpha}|\phi\rangle}{\langle\phi_{T}|\phi\rangle}$. The introduction of the FB does not add any additional approximations, as the integration variable in (4) is merely shifted by a constant, yet it is crucial for two reasons. First, it diverges when the “nodal surface” (as defined in the complex plane of overlaps) is approached, pushing the walker away from the origin. Second, since $\hat{v}_{\alpha}$ is complex, $Im[{\mathbf{\bar{x}}}]$ reduces the amount of physical information discarded in the phase projection. Similarly, the subtraction of $\langle\hat{v}\rangle$ from $\hat{v}$ in the propagator also greatly reduces the severity of this projection, as the smaller diagonal matrix elements of the resulting propagator cause milder rotations of the phases of the orbitals.Al-Saidi et al. (2006c)

The choice of FB allows the weight factor which multiplies the previous weight after a propagation step to be written as $W(\phi)=e^{-\Delta\tau E_{L}(\phi)}$. The second approximation in ph-AFQMC, which is much milder than the first, takes the real part of the local energy in the weight above:

[TABLE]

The severity of the phaseless constraint in ph-AFQMC can be reduced with the use of more accurate trial wavefunctions, especially those with the correct symmetry properties.Shi and Zhang (2013) This can be seen from (10), for when $|\phi_{T}\rangle=|\Phi_{0}\rangle$, the local energy which determines the weights and energy measurements is a real constant (equal to the exact ground state energy).

II.2 Correlated Sampling Methods for AFQMC

Given that the statistical error in an AFQMC calculation arises solely from the MC evaluation of the integral over AFs, a natural way to implement correlated sampling in the calculation of an energy difference is to pair the walkers of the two systems, and use the same set of AFs to propagate each pair of walkers. To be precise, the correlation is established on the level of each term in (2), for each value of $\alpha$. This, in essence, matches Cholesky vectors of the same iteration, and becomes more effective the more similar the interactions are in the two systems.

Sampling in this way, however, requires one to relinquish optimal importance sampling, which is usually implemented via a population control (PC) scheme in which walkers with large (small) weights are replicated (stochastically purged), since performing independent PC for each system separately would quickly destroy the walker-pair correlation. Alternatively, one can implement PC with respect to the weights of the primary system for both systems. This, however, will only be effective if the two systems are essentially identical. In the absence of a PC scheme, the noise from the accumulation and persistence of divergent walkers inevitably grows with propagation time. We find that the immediate reduction in statistical error following the correlation of the AFs, augmented as needed by what we call the “preliminary equilibration scheme” below, allows converged averages to be obtained at short and intermediate projection times. In light of the fact that the stability of the random walks at long times, as afforded by a branching scheme, appears to be only marginally relevant when our correlated sampling approach is used, we simply choose not to implement PC when the AFs are correlated. In the event that a walker’s weight becomes zero or negative (in the latter case the phase projection sets the weight to zero), the walker is no longer propagated and the random number stream is updated if necessary such that the correlation of the other walker-pairs is unaffected.

The data shown in blue in Fig. 1 illustrates features associated with the typical correlated sampling protocol used in this work. In this example, the propagation in imaginary time (measured in Ha*-1*) is performed using correlated AFs and repeated 11 times using different random number seeds. The mean energy difference and associated standard error at each $\tau$ point is computed among the repeats (see left plot). We then choose the imaginary-time at which the energy difference is seen to stabilize (in this case at $\tau\sim 4$), and for each repeat calculate the cumulative average at each $\tau$, which represents the running average of the energy measurements taken after the end of the equilibration period up to the given value of $\tau$. To obtain the final result, we compute the mean and standard error of the cumulative averages among the repeats (see right plot), and choose the value corresponding to the $\tau$ at which the standard error reaches a minimum (if there are multiple equivalent minima the one occurring earliest is chosen) or falls below a target error level. We note that our choice of 11 repeats is arbitrary; a larger number could be used to reduce the standard error as necessary. Clearly in this example correlating the AFs for $\tau\geq 0$ drastically reduces the error relative to the uncorrelated runs shown in red, and the resulting IP agrees with that obtained from independent $\tau=80$ ph-AFQMC runs of the neutral and cationic species. This benchmark, indicated by the solid black line, was found to have a negligible standard error of 0.2 mHa after employing a reblocking analysis which corrects for auto-correlation.Flyvbjerg and Petersen (1989)

As it is often the case that a very small population size is sufficient to achieve a desired statistical error via the correlated sampling approach, we choose not to employ PC in the uncorrelated comparisons since this would result in a bias that typically goes as 1/$N_{wlk}$.Nguyen et al. (2014) A second reason is that over the relatively short imaginary-time scales relevant for the correlated sampling method, PC is expected to have little effect on the uncorrelated comparisons as the walker weights usually do not stray far from unity. This is confirmed by the similarity of the error curves plotted in the insets of Fig. 1, corresponding to uncorrelated runs with (dotted green) and without (red) PC. The former is obtained by using a large enough population size such that the bias from the PC algorithm is negligible (360 walkers per repeat), and rescaling the resulting standard error by $\sqrt{\frac{360}{12}}$.

If one or both of the comparative systems requires a long equilibration time, which can be the case for, e.g., initial populations which are severely spin-contaminated or for strongly correlated systems in which the “guiding” trial function poorly describes the true ground state, then walker pair correlation can be lost prior to convergence, and the associated noise growth can make measurements impossible. The simplest way to overcome this problem is to use a better trial function, e.g. a multi-determinant CASSCF wavefunction instead of the single Hartree-Fock (HF) determinant. Alternatively, we have devised what we will refer to as the “preliminary equilibration scheme” (PES). First, one of the two systems is equilibrated using PC for the required interval, then the walkers of the secondary system are initialized with the resulting determinants of the equilibrated primary system, with weights scaled by $\frac{\langle\phi_{T}^{secondary}|\phi\rangle}{\langle\phi_{T}^{primary}|\phi\rangle}$ (and any resulting phases projected) to reflect the appropriate importance sampling. Finally, both systems are propagated using correlated sampling with their respective Hamiltonians without PC for a short period of time after which energy measurements are collected. A schematic of this procedure is shown in Fig. 2. This protocol is repeated with different random number seeds to obtain satisfactory statistics on the energy difference between the ground states. We note that a similar scheme was published many years ago by Traynor et al.Traynor and Anderson (1988)

Inspection of the ph-AFQMC propagator

[TABLE]

demonstrates that stricter approaches to correlated sampling exist. In one such alternative to correlating only the AFs, each pair of walkers is propagated using both the same AFs and FBs. However taking the simple average of the FBs of the primary and secondary systems was found to cause an early onset of noise even for marginally different systems. This issue likely arises from divergences in the local energies and weights which can occur more frequently due to the use of sub-optimal FBs, a fact which follows from the discussion centered around Eq. (9). Separately, we choose not to correlate the $\langle\mathbf{\hat{v}}\rangle$ since any gain in sampling efficiency would come at the cost of an increased bias from the phaseless constraint.Al-Saidi et al. (2006c)

Finally, we note that our method of correlated sampling is not limited to the phaseless version of AFQMC, and can be used in exactly the same manner for FP calculations. In fact, using the same AFs in the latter case results in a relatively stricter form of correlated sampling, since the FB is not present in the post-HS propagator used in FP.

II.3 Utilizing Optimal Correlation in Molecular Applications

Having justified our choice to correlate only the AFs, we claim that maximal correlation between a pair of walkers is achieved when the primary and secondary systems use the same set of basis functions. To see this, recall the definition of the two-electron matrix elements in the Hamiltonian (1):

[TABLE]

When $\{\phi\}_{primary}=\{\phi\}_{secondary}$, the $V_{ijkl}$ and, in turn, the one-body operators $\mathbf{\hat{v}}$ in the propagators (11) for both systems will be identical. This ideal condition is obviously satisfied when the two systems are an atom or molecule with the same geometry but e.g. different charges, since atomic basis functions such as Gaussian-type orbitals are usually centered on the positions of the nuclei. In fact, for the calculation of vertical IPs and EAs, exactly the same $T_{ij}$ and $V_{ijkl}$ elements are used in the imaginary-time propagation of both systems. For adiabatic redox processes, however, the ground-state geometries of the neutral and ionized species are, in general, different. In most cases the differences in the $V_{ijkl}$ are slight, and correlated sampling is still found to be effective, albeit with reduced efficiency. For example, Fig. 3 compares the correlated and uncorrelated standard errors corresponding to the vertical and adiabatic IPs of methanol. While the errors from the uncorrelated runs are roughly similar in magnitude, that from the correlated vertical case is significantly smaller and more constant compared to that from the correlated adiabatic case. This suggests the use of a two-step process to compute adiabatic energy differences, in which the fixed-geometry transition is calculated with correlated sampling-based AFQMC, and the geometry relaxation energy is obtained from a lower level of theory. Further details will be presented in Section III.B.

In the same spirit, we have devised a protocol which enables the calculation of the change in energy corresponding to the removal of a proton while the rest of the geometry remains fixed, in which the optimal condition mentioned above is realized with the use of so-called “ghost” basis functions. In this scheme, the deprotonated species uses the same set of basis functions as the acid, i.e. the position of the removed proton still serves as a center for hydrogen basis functions but not as a center of nuclear charge. As a result, the number of basis functions remains the same for the primary and secondary systems as required by our correlated sampling method. Moreover the $V_{ijkl}$ and thus the $\mathbf{\hat{v}}$ are, by construction, identical. We note that even though the $T_{ij}$ now differ 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 evaluation of the integral over AFs. Fig. 4 illustrates the efficacy of this procedure in a calculation of the deprotonation energy for methanol, showing a drastic reduction in the statistical errors when the AFs are correlated.

The “ghost” basis function strategy can also be directly applied effectively to hydrogen abstraction reactions. Fig. 5 demonstrates the error reduction afforded in a calculation of the MeOH $\rightarrow$ MeO$\cdot$ reaction energy.

Thus we have presented a correlated sampling approach which has the potential to reduce the statistical error for redox, deprotonation, and hydrogen abstraction reactions. The extension of this protocol to obtain adiabatic energy differences will be described in Section III.B.

II.4 Computational Details

One-electron and overlap integrals, Cholesky vectors, restricted open-shell and unrestricted HF trial wavefunctions were all obtained from a modified version of NWChem.Valiev et al. (2010); Purwanto et al. (2011) The maximum residual error in the Cholesky decomposition was chosen to be $1\times 10^{-6}$ Ha. This cutoff has a negligible effect on the QMC energies (i.e. orders of magnitude smaller than the statistical error bars), and is utilized to decrease the number of AFs required, which in turn decreases both the error from the random noise and the computational cost of the calculations.

For calculations using a single-determinant trial wavefunction with unrestricted reference, spin-contamination can occur due to the presence of higher multiplicity spin configurations in the trial function, which can prolong equilibration times and reduce accuracy. We use a spin-projection techniquePurwanto et al. (2008) in AFQMC to minimize such effects. In this scheme, the walkers are initialized with a restricted open-shell HF determinant, which is an eigenfunction of $\hat{S}^{2}$, such that propagation with $e^{\sqrt{\Delta\tau}\mathbf{x}\cdot(\mathbf{\hat{v}}-\langle\mathbf{\hat{v}}\rangle)}$ preserves spin symmetry. The unrestricted HF determinant, while not an eigenfunction of the total spin operator, is known in most cases to provide a better description of the ground-state energy, and is therefore used to implement the phaseless constraint and to estimate the energy.

Multi-determinant trial states were obtained from CASSCF calculations performed with PySCF.pys The resulting expansions of the wavefunctions in CI space are truncated such that determinants associated with coefficients below a specified threshold are discarded. In what follows, CASSCF($X$,$Y$) will denote an active space with $X$ electrons and $Y$ orbitals; for cases in which an energy difference is calculated we use the notation $X\equiv N_{neutral}/N_{ion}$ and $Y\equiv M_{neutral}/M_{ion}$ to denote the numbers of electrons and orbitals in the neutral and charged systems, respectively.

In all AFQMC calculations the spin -up and -down sectors of the walker determinants are separately orthonormalized via a modified Gram-Schmidt procedureWhite et al. (1989) every 0.05 Ha*-1* to preserve the anti-symmetry of the walker determinants. When required, we use a PC algorithm in which the number of walkers is fixed throughout the entire calculation.Nguyen et al. (2014) PC (when used) and energy measurements are performed every 0.1 Ha*-1*.

Use of the “hybrid” formulationPurwanto et al. (2009) of ph-AFQMC allows for the calculation of the local energy, the bottleneck in the algorithm’s scaling, at intervals rather than at every step since $E_{L}$ is not required to compute the weight factors, which are now calculated explicitly as $\frac{\langle\phi_{T}|\phi^{(\tau+1)}\rangle}{\langle\phi_{T}|\phi^{(\tau)}\rangle}e^{\mathbf{x}\cdot\mathbf{\bar{x}}-\mathbf{\bar{x}}^{2}/2}$. We chose to implement the hybrid method primarily because it is much faster than the local energy method, and offers more flexibility when devising correlated sampling strategies.

The comparison of QMC energies with experimental results in general requires extrapolations to mitigate errors associated with the use of a finite time step and basis set. The Trotter error, due to the decomposition in (3) and three bounding conditions,Purwanto et al. (2009) can be estimated from a linear extrapolation of the energy differences from three independent simulations at $\Delta\tau=0.02,\ 0.01,\ \text{and}\ 0.005$ Ha*-1*. For the IPs and EAs of the G2 atomic test set, the energies resulting from the 0.01 Ha*-1* time step were found in all cases to be well within 1 mHa of the $\Delta\tau$-extrapolated value. Hence, we use only the $\Delta\tau=0.01$ value to produce the results in Section III. Extrapolation to the complete basis set (CBS) limit was performed following the second order Møller-Plesset perturbation theory (MP2)-assisted protocol detailed in Ref. Purwanto et al. (2011) AFQMC calculations were performed at $x=$ 3 (where $x$ is the cardinal number of the basis set); UHF calculations were run at $x=$ 2-5 and MP2 calculations at $x=$ 3, 4 using GAMESSSchmidt et al. (1993) or, for the reactions which utilize “ghost” basis functions, NWChem. IPs were computed using the cc-pVxZ basis,Dunning Jr (1989) while EAs and the deprotonation energies to be compared with experimental data were computed using the aug-cc-pVxZ basis sets.Kendall et al. (1992)

QMC statistical error bars were propagated through the data analysis procedure. In assessing the statistical error on the deviation of QMC from experiment, we assume that there are negligible uncertainties associated with 1) the experimental measurements, 2) the exponential (HF) and linear (MP2) fits in the CBS extrapolation procedure where we found that the former error is an order of magnitude or more smaller than the QMC error in representative cases, and 3) the scaling factor for the zero-point energies. With regard to 1), experimental uncertainties can be found in the NIST database,Johnson III and are ignored in order to isolate the statistical error associated with the QMC measurements, since the quantities presented in this work are differences between calculated energies and experimental measurements.

Our Fortran 90 code is linked with OpenBLASWang et al. (2013); Xianyi et al. (2012) and Expokit.Sidje (1998) Random numbers were generated with the 48 bit Linear Congruential Generator with Prime Addend, as implemented in SPRNG5.Mascagni and Srinivasan (2000)

III Results

III.1 Atomic IPs and EAs

In this section, we use our correlated sampling-based AFQMC approach to calculate the IPs and EAs of the 1st row atoms included in the G2 Ion Test set, which have experimental uncertainties of less than 0.05 eV.Curtiss et al. (1998) The deviations from experiment within calculated AFQMC results are presented in Tables 1 and 2 alongside the deviations resulting from G2 theory,Curtiss et al. (1991a) and DFT with the B3LYP exchange-correlation functionalBecke (1993); Lee et al. (1988) in the 6-311+G(3df,2p) basis.Curtiss et al. (1998)

For the calculations which used the unrestricted HF determinant as the trial function, PES was used for the IPs of B, C, N, and O with equilibration times of 35, 20, 15, and 15 $Ha^{-1}$, respectively, and similarly for the EAs of C and F with 15 and 10 $Ha^{-1}$, respectively. For the remaining species, equilibration was facile and thus the AFs were correlated from the beginning of the imaginary-time propagation. Notable deviations from the experimental values are found for the IPs of B and N, and the EA of F when the UHF state is used as the trial function (nearly identical errors have been previously reported within ph-AFQMCAl-Saidi et al. (2006c)). These discrepancies are resolved in the correlated sampling-based FP results. Since the energies from FP are not biased by the phaseless constraint and thus insensitive to the trial function used, (any) small residual errors can be attributed to the MP2-assisted CBS extrapolation scheme.

The long equilibration times and inaccuracies encountered in the above cases are manifestations of the fact that single-determinant trial functions obtained from mean-field calculations are generally not well-suited to describe the open-shell systems involved in redox reactions. For instance in the IP of B, comparing FP and ph-AFQMC calculations for both the neutral and cationic species exposes the fact that the error in the IP stems from inaccuracies in the computation of the total energy of the neutral B atom, which has a single unpaired electron in the triply-degenerate $p$ orbital manifold. The use of trial wavefunctions with proper symmetry properties has previously been shown to lead to improved accuracy within the phaseless approximation,Shi and Zhang (2013) so we now consider trial functions which are eigenfunctions of $\hat{S}^{2}$ for use with the phaseless constraint. While employing a restricted open-shell HF trial affords no appreciable gain in accuracy, the use of a truncated CASSCF trial with a very modest active space size and a small number of determinants is sufficient to eliminate most of the error for the IPs of both B and N.

The accurate description of the EA of F proves to be more demanding. It is possible that the phase problem is particularly problematic for the case of F- due to the small atomic radius which may lead to very strong electronic correlations. Such sizable correlations manifest themselves in the algorithm as a large imaginary component of the propagator, which in turn leads to the elimination of crucial physical information during the projection to the real axis if the trial function does not adequately provide the gauge information on the Slater determinants of the ground state. In such cases, excitations into a large number of virtual orbitals will make significant contributions to the correlation energy, and therefore a large active space is required to generate the CASSCF trial function so that the resulting phase projections are sufficiently benign.

In general, the implementation and efficacy of the correlated sampling scheme presented in this work remain unaltered in the case of a multi-determinant trial function since walker determinants are propagated in exactly the same manner. We do find a reduction in the required equilibration times for all multi-determinant ph-AFQMC calculations that we have performed, relative to single-determinant calculations of the same systems, rendering the use of PES unnecessary. In addition, the use of multi-determinant trial functions results in drastically smaller error bars, even when more than 4x fewer walkers are used. Increasing the quality of the CASSCF trial function is a promising way to systematically reduce the error from the phaseless constraint, and we find that in all cases the use of a reasonable number of determinants produces redox energies within the maximum experimental error of 0.05 eV for the G2 Test Set.

III.2 The Case of Methanol

Motivated by the arguments set forth in Section II.C, we now describe the details of a composite method for calculating the adiabatic IPs, deprotonation free-energies, and hydrogen-dissociation energies of molecular systems, and illustrate the accuracy of our approach on the case of methanol. We utilize a stepwise process consisting of: 1) the fixed-geometry process calculated with correlated sampling-based ph-AFQMC, and 2) a geometry-relaxation step calculated within MP2. Optimal correlation can be achieved in 1) since the same set of basis functions is used, while in 2) we expect large error cancellation resulting from the fact that the initial and final ground-state geometries are typically very similar in redox, deprotonation, and hydrogen abstraction reactions. The calculations are performed in a triple-zeta basis with additional diffuse basis functions for all species involved in the deprotonation reaction. The resulting energy differences from these two steps are added together, and the endpoints are extrapolated to the CBS limit.

For all reaction types, zero-point energies are calculated at the level of HF/6-31G* and scaled by a factor of 0.899 to account for anharmonicity and the known shortcomings of HF theory, following the G2 protocol. Deprotonation free-energy results incorporate the value of -6.28 kcal/mol as the free-energy of a proton at 298 K,Tawa et al. (1998) and we use the exact ground-state energy of H$\cdot$ (-0.5 Ha) to calculate the bond dissociation enthalpy of the O-H bond.

Table 3 illustrates the accuracy of our correlated sampling protocol with respect to experimental results for all three reaction types. The quality of our IP and deprotonation free-energy results surpass that of the more costly G2 method, while the O-H bond dissociation result has an error of comparable magnitude. We note that the single-determinant trial wavefunction is sufficient to produce a near-exact deprotonation free-energy, which we attribute to the fact that in this reaction type both the protonated and deprotonated species are closed-shell. Moreover, our correlated sampling-based ph-AFQMC calculations, using an extremely modest number of walkers, predicts the energy differences corresponding to all of these reaction types to within chemical accuracy, which is not the case for the G2 method.

III.3 Basis Set Size, Number of Random Walkers, and CPU-time Reduction

III.3.1 Basis Set Size

The statistical noise of individual AFQMC runs is expected to increase with the number of AFs, $\alpha$, yet the precise scaling is subtle. Each AF contributes additional noise. On the other hand, the magnitude of the contribution is moderated by a degree which depends on the form of the interaction and the decomposition procedure leading to Eq. (5). For example, the error bar is seen to change little beyond a modest cutoff in plane-wave AFQMC.Suewattana et al. (2007) In principle $\alpha$ grows as $M^{2}$, but our (rather conservative) truncation of the Cholesky decomposition via the cutoff mentioned in Section II.D results in $\alpha\sim 10M$. As a result, we expect the statistical error in an uncorrelated calculation to increase with $M$ before saturating toward the CBS limit. This could be problematic given that large basis sets containing functions associated with high angular momenta are frequently required to accurately describe the correlation energy,Dunning Jr (1989) and also given the fact that most interesting applications involve large systems.

Here we show that the magnitude of the reduction in statistical error enabled by the use of correlated sampling grows with $M$, such that the error is nearly independent of $M$. Fig. 6 shows the standard errors resulting from calculations of the IP of the K in the 6-31G*, 6-31+G*, and 6-311+G* basis sets (which consist of 23, 35, and 45 basis functions, respectively); Fig. 7 illustrates the same effect for fixed-geometry deprotonation reactions of water, methanol, and ethanol (58, 116, and 174 basis functions). For both reaction types, while the errors from the uncorrelated calculations increase significantly with $M$, we find that those resulting from correlated sampling remained roughly constant. We anticipate that this finding will hold in general, and will be of crucial importance in future applications of correlated sampling-based AFQMC to larger molecules.

III.3.2 Reduction in the Number of Random Walkers and CPU-Time

In an uncorrelated QMC calculation for a fixed length of propagation time, the resulting standard error can be reduced by increasing the number of random walkers, $N_{wlk}$, used in the MC evaluation of the HS integral in Eq. (5). However, given that the required computational expense increases linearly with $N_{wlk}$, using a brute-force approach that simply increases $N_{wlk}$ is less practical for many systems. For energy differences, correlated sampling provides a much cheaper alternative as it allows for a dramatic reduction in the $N_{wlk}$ required to achieve a given statistical error.

The errors of the cumulative averages of the IP of K in the 6-31+G* basis (M=35) and the fixed-geometry deprotonation of methanol in the cc-pVTZ basis (M=116) are shown for different values of $N_{wlk}$ in Figs. 8 and 9, respectively. For both reaction types, we make the following observations: First, for a given $N_{wlk}$, the standard error is significantly lower in the correlated sampling case. Second, the magnitude of this reduction is greater for smaller $N_{wlk}$. Finally, while the standard errors of the uncorrelated runs increase as $N_{wlk}$ is reduced, in the correlated sampling runs the error is relatively insensitive to $N_{wlk}$.

In light of these findings, we are in a position to understand how a reduction in the standard error due to correlating 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 QMC energy (which defines the 99$\%$ confidence interval as the cumulative average of the energy $\pm$ 1 kcal/mol), and compare the total CPU-time required to achieve this via correlated and uncorrelated approaches. In the former case, we find that using only 24 walkers in each of the 11 repeat calculations is sufficient to achieve the target error and a resulting energy in agreement (i.e. within the 99$\%$ confidence interval) with a benchmark result obtained with 5040 walkers. In fact, as few as 6 walkers produced the same level of accuracy in some cases. The inset of Fig. 8 shows that the statistical error falls below the target at $\tau\sim 7$. The total CPU-time required to propagate 11 repeats for this length of imaginary-time is 41.2 minutes on a single 2.60 GHz Intel Xeon processor. Using the same number of walkers without correlating the AFs, we find that the target error is not reached even after 200 $Ha^{-1}$. Using 1056 walkers gives rise to a standard error that falls below 0.5 mHa after 10 $Ha^{-1}$, as shown in the inset of Fig. 8, and a resulting QMC energy that is in agreement with the benchmark result. This calculation takes 2262.6 minutes on a single processor, and we thus conclude that our correlated sampling approach reduces the CPU-time by a factor of approximately 55.

For the deprotonation of methanol we use a target error of 1 mHa. As shown in the inset of Fig. 9, both the uncorrelated run with 2400 walkers and the correlated run with 96 walkers yield results that fall below our target error at $\tau\sim 10$. We use the fact that the total CPU-time is proportional to the product of the number of walkers and the propagation time to estimate that correlating the AFs reduces the CPU-time by a factor of approximately 25. Currently the total calculation, including all 11 repeats, requires $\sim$154 hours on a single CPU core. We perform a similar analysis for the dissociation of H$\cdot$ from the O-H bond of methanol. Due to the relatively large computational cost of using a CASSCF trial function (with the same active space and truncation scheme described in Table III), we increase our target error to 2 mHa. Figure 10 shows that the errors on the cumulative averages from both an uncorrelated run with 288 walkers and a correlated run with 12 walkers fall below our target error at $\tau\sim 6$. Thus we deduce a speed-up factor of 24 for this H abstraction reaction. On a single CPU core this requires $\sim$157 hours.

As the previous sections have shown, the magnitude of the reduction in the statistical error, and consequently the CPU-time, as compared with an uncorrelated calculation depends on the number of walkers used and the size of the basis set. In addition, one intuitively expects our correlated sampling approach to work better for systems that more closely resemble each other. This is crudely the case, yet the subtleties involved in such a claim warrant further discussion. Indeed, referring to the propagator in (11), while our correlated sampling method ensures that the $\mathbf{\hat{v}}$ operators and the AFs, $\mathbf{x}$, are the same for both the primary and secondary systems, the FBs, $\mathbf{\bar{x}}$, as defined in (9) and the expectation values with respect to the trial functions $\langle\mathbf{\hat{v}}\rangle$ will in general be different. In the limit that the primary and secondary systems are identical, the entire propagator in (11) is identical for both systems, and the statistical error in the energy difference will be exactly and trivially zero as a result of perfect walker-pair correlation. Otherwise, the reduction in statistical error afforded by our correlated sampling approach becomes less pronounced the more the trial wavefunctions of the primary and secondary systems differ, since the the FBs and $\langle\mathbf{\hat{v}}\rangle$ are expectation values that depend explicitly on the trial wavefunctions. It is encouraging to note that despite the differences in trial functions correlating only the AFs yields such large speed-ups in CPU-time. While additionally correlating the $\langle\mathbf{\hat{v}}\rangle$, possibly by using some combination of the trial functions for both systems, would compromise the accuracy of the phaseless approximations, future studies will explore optimal ways to pair walkers such that the similarity in the FBs of walker pairs is maximized. It is encouraging that even for systems for which the electronic energies of the primary and secondary systems differ by some 200 eV, as is the case in the deprotonation of methanol, our correlated sampling approach still yields dramatic efficiency gains.

IV Conclusions and Outlook

In this work we have devised a correlated sampling protocol for the calculation of chemically relevant energy differences within the exact and phaseless variants of AFQMC. For molecules we utilize a two step strategy in which optimal walker-pair correlation is achieved in the ph-AFQMC description of the fixed-geometry process, while the geometry relaxation energy is calculated with the confines of MP2. Together with an MP2-assisted CBS extrapolation method we obtain calculated IPs, EAs, deprotonation free-energies, and bond dissociation energies that are in excellent agreement with experiments. Moreover, our correlated sampling approach yields large reductions in the statistical errors relative to those obtained from uncorrelated approaches. In contrast to uncorrelated AFQMC, where the error bars are found to increase with system size and/or the number of basis functions, correlating the AFs keeps the statistical error relatively constant as chemical complexity increases. In addition, our approach drastically reduces the number of walkers required to achieve a given statistical error target, which translates into large reductions in CPU-time.

We utilize a “ghost” basis function strategy that enables the application of a correlated sampling-based approach to processes involving large energetic changes. Indeed, given that the correlated sampling scheme outlined here is successful for electron, proton, and H$\cdot$ transfer reactions, we are optimistic that other chemical changes are within reach. In a future work we will systematically investigate chemical reactions which involve substantial changes in geometry, including the addition/removal of larger, more complex functional groups. Along the same lines we are optimistic about the savings that our correlated sampling approach may yield when basis functions that are independent of the nuclear coordinates, such as plane-waves, are used. We anticipate that the insensitivity of the statistical error to basis set size that we observe in this work may partially or totally offset the relatively large number of plane-waves typically required for convergence.

Additional improvements in computational efficiency are necessary before our correlated sampling-based AFQMC approach can be routinely applied to the large, complex systems that we ultimately hope to investigate. In addition to incorporating orbital localization techniques and experimenting with more compact basis sets, we will pursue a variety of algorithmic (e.g. the implementation of mixed-precision and sparse linear algebra routines) and hardware (e.g. porting the code to run on GPUs) optimizations. Future targets of investigation include the IP of B2F4,Chan et al. (2012); Li and Fan (2002) redox potentials and pKa’s of TM clusters and battery electrolytes,Coskun et al. (2016); Jerome et al. (2014); Roy et al. (2009); Konezny et al. (2012); Baik and Friesner (2002); Borodin et al. (2013); Shao et al. (2012); Qu et al. (2015) and the relative energies of low-lying states of interacting TM centers such as Fe-S complexes and the active site of Photosystem II.Sharma et al. (2014); Ames et al. (2011); Siegbahn and Blomberg (2000)

V Acknowledgements

JS gratefully acknowledge Fengjie Ma, Mario Motta, Hao Shi, and James Dama for their generous help and guidance in implementing and understanding AFQMC methods. DRR acknowledges funding from grant NSF-CHE-1464802. The work of SZ was supported by grants DE-SC0001303 and DE-SC0008627.