Auxiliary-Field Quantum Monte Carlo#

Overview#

The AFQMC method is an orbitally-based many-body method which is formulated in terms of a generic, interacting second-quantized Hamiltonian. We refer the reader to one of the review articles on the method [MZM20, PZ04, Zha13, ZK03] for a detailed description of the algorithm. Here, we provide an overview of AFQMC with a focus on how the formalism and the inputs to the SAFIRE code relate to each other.

AFQMC uses imaginary-time projection to generate a many-body wavefunction, \(| \Psi \rangle\), starting from an initial wavefunction, \(| \Psi_I \rangle\), which has non-zero overlap with the target wavefunction - i.e. \(\langle \Psi | \Psi_I \rangle \neq 0\)

(11)#\[\lim_{\beta \to \infty} e^{-\beta \hat{H}} | \Psi_I \rangle \rightarrow | \Psi \rangle\]

where \(\beta\) is the total projection time, and \(\hat{H}\) is the interacting Hamiltonian which, in its most general form, is given by

(12)#\[\begin{split}H =& H_0 + \hat{H}_1 + \hat{H}_2 \\ =& E_0 + \sum\limits_{il\sigma\sigma'} h^{\sigma \sigma'}_{il} \hat{c}^\dagger_{i\sigma}\hat{c}_{l\sigma'} + \sum\limits_{ijkl\sigma\sigma'\sigma'\sigma} v_{ijkl}^{\sigma\sigma'\sigma'\sigma} \hat{c}^\dagger_{i\sigma} \hat{c}^\dagger_{j\sigma'}\hat{c}_{k\sigma'}\hat{c}_{l\sigma}\end{split}\]

and \(\hat{c}^\dagger_{i\sigma}\)/ \(\hat{c}_{i\sigma}\) are fermionic creation/annihilation operators which create/annihilate a particle in orbital \(i\) with spin \(\sigma\), \(H_0\) contains all constant contributions to the Hamiltonian, \(\hat{H}_1\) contains all one-body Hamiltonian terms (for example, the kinetic energy in ab initio calculations, or the hoping matrix in lattice models ), and \(\hat{H}_2\) contains all two-body interactions such as the Coulomb interaction in ab initio calculations or the Hubbard U in lattice models. The matrix elements of \(\hat{H}_1\) are given generically as \(h_{ij}^{\sigma\sigma'}\) and the details of which terms are included depends on the specific system of interest. Similarly for the matrix elements of \(\hat{H}_2\) which are generally \(v_{ijkl}^{\sigma\sigma'\sigma'\sigma}\). In specific calculations, the general form of the Hamiltonian simplifies considerably especially for \(\hat{H}_2\). The different forms of Hamiltonian implement in SAFIRE are described in the The Hamiltonian file formats section.

Imaginary-time projection#

Thouless’ theorem [Tho60] states that the operation of the exponential of a one-body operator on a Slater determinant is simply another Slater determinant. However, the application of the imaginary-time projection operator is non-trivial due to the interaction term in the hamiltonian. To handle the interaction term, the projection operator is mapped onto a high-dimensional integral where the integrand consists of exponentials of one-body operators only. This achieved as follows. First, the imaginary-time projection is divided into \(N_{steps}\) imaginary-time steps of size \(\tau\) which allows a Trotter–Suzuki (TS) [Suz76, Tro59] decomposition to be applied such that the propagator now has the form,

(13)#\[e^{-\tau \hat{H}} \approx e^{-\tau \hat{H}_1/2} e^{-\tau \hat{H}_2} e^{-\tau \hat{H}_1 / 2} + \mathcal(O)(\tau^3)\]

where the accuracy of the TS decomposition can be improved by using a smaller time step. Next, the interaction term is rewritten as a quadratic form of one-body operators,

(14)#\[\hat{H}_2 = \sum_\lambda (\hat{v}^\gamma)^2,\]

where \(\hat{v}^\gamma\) are one-body operators. Next, \(exp[-\tau \hat{H}_2]\) can be rewritten via a Hubbard-Stratonovich (HS) transformation as

(15)#\[e^{-\tau \sum_\lambda (\hat{v}^\gamma)^2} = \int d\mathbf{\sigma} P(\mathbf{\sigma}) e^{-\tau \mathbf{\sigma} \cdot \hat{\mathbf{v}}},\]

where \(\mathbf{\sigma}\) is a vector of auxiliary fields, \(\hat{\mathbf{v}}\) is the vector of one-body operators \(\hat{v}^\gamma\), and \(P(\mathbf{\sigma})\) is a normal distribution function. Combining the TS decomposition and the HS transformation, we arrive at

(16)#\[\begin{split}e^{-\tau \hat{H} } \approx& \int d\mathbf{\sigma} P(\mathbf{\sigma}) e^{-\tau \hat{H}_1/2} e^{-\tau \mathbf{\sigma} \cdot \hat{\mathbf{v}}} e^{-\tau \hat{H}_1/2} , \\ \approx& \int d\mathbf{\sigma} P(\mathbf{\sigma}) B(\mathbf{\sigma})\end{split}\]

where \(B(\sigma) = exp^{-\tau \hat{H}_1/2} exp^{-\tau \mathbf{\sigma} \cdot \hat{\mathbf{v}}} exp^{-\tau \hat{H}_1/2}\). Additionally, the many-body wavefunction is represented, at time step \(s\), in an over-complete basis of non-orthogonal Slater determinants, \(| \Phi^s_n \rangle\), as,

(17)#\[| \Psi^s \rangle \doteq \sum_n W^s_n | \Phi^s_n \rangle,\]

where \(W^s_n\) is a weight which is accumulated during the projection.

Finally, the many-body wavefunction can be sampled as

(18)#\[e^{-\tau \hat{H} } | \Psi^s \rangle \approx \int d\mathbf{\sigma} P(\mathbf{\sigma}) B(\mathbf{\sigma}) | \Psi^{s-1} \rangle\]

by sampling auxiliary-field configurations from \(P(\mathbf{\sigma})\), constructing \(B(\mathbf{\sigma})\) with the sampled set of \(\mathbf{\sigma}\), and evaluating \(B(\mathbf{\sigma}) | \Phi^s_n \rangle\) for each walker. Observables are computed on the fly to achieve low-order polynomial scaling with system size. For example, the energy is evaluated as

\[E = \frac{\langle \Psi_\mathrm{T} | \hat{H} | \Psi^s \rangle}{\langle \Psi_\mathrm{T} | \Psi^s \rangle} \approx \frac{1}{\sum_n W^s_n} \sum_n W^s_n \frac{\langle \Psi_\mathrm{T} | \hat{H} | \Phi^s_n \rangle}{\langle \Psi_\mathrm{T} | \Phi^s_n \rangle},\]

Importance sampling#

In practice, importance sampling is used to reduce variance. The auxiliary fields from the Hubbard-Stratonovich transformation are shifted by the so-called “force bias”, \(\bar{\mathbf{\sigma}}\) and the walkers are re-weighted according to an importance function, \(I(\mathbf{\sigma},\bar{\mathbf{\sigma}},\Phi^{s-1}_n)\), at each projection step. Formally, the force bias is given by

(19)#\[\bar{\mathbf{\sigma}}=-\sqrt{\tau}\frac{\langle\Psi_\mathrm{T}|\hat{\mathbf{v}}|\Phi^s_n \rangle}{\langle\Psi_\mathrm{T}|\Phi^s_n \rangle},\]

where \(|\Psi_\mathrm{T} \rangle\) is the trial wavefunction which is essentially a guess of the true many-body wavefunction. The walkers weights are updated at each projection step based on the importance function as

(20)#\[W^s_n = I(\mathbf{\sigma},\bar{\mathbf{\sigma}},\Phi^{s-1}_n) W^{s-1}_n,\]

where \(I(\mathbf{\sigma},\bar{\mathbf{\sigma}},\Phi^{s-1}_n)\) is given by

(21)#\[I(\mathbf{\sigma},\bar{\mathbf{\sigma}},\Phi) \equiv \frac{\langle\Psi_\mathrm{T}|\hat{B}(\mathbf{\sigma} - \bar{\mathbf{\sigma}})|\Phi \rangle} {\langle\Psi_\mathrm{T}|\Phi \rangle} \textrm{Exp}\left(\bar{\mathbf{\sigma}} \cdot \mathbf{\sigma} - \frac{\bar{\mathbf{\sigma}} \cdot \bar{\mathbf{\sigma}}}{2}\right)\]

For lattice models, some of the details are slightly different since the form of the interactions is simpler.

The Sign/Phase Problem#

As formulated above, the imaginary-time projection suffers from an exponential decrease in the signal-to-noise ratio which arises from the invariance of any observable to an overall phase of the many-body wavefunction, \(e^{-i \theta} | \Psi \rangle\). This is a manifestation of the generic fermionic sign problem that arises in all quantum Monte Carlo approaches. The projector, \(B(\mathbf{\sigma})\), is typically complex-valued and introduces an arbitrary phase to the Slater determinant random walkers. As the projection proceeds, the random walkers accumulate random phases, leading to a finite density of walkers in the complex plane defined by \(\langle \Psi_\mathrm{T} | \Phi^s_n \rangle\). This leads to large fluctuations in the weights that ultimately diverge. This is known as the phase problem in the AFQMC literature [ZK03]. For lattice models, it is possible to choose a HS transformation such that \(B(\mathbf{\sigma})\) is real-valued. In this case, the random walkers accumulate random signs — i.e., \(\theta \in \{0, \pi\}\) — leading to the fermionic sign problem.

The sign/phase problem is controlled using importance sampling, as shown above, and by imposing a constraint based on the trial wavefunction, \(| \Psi_\mathrm{T} \rangle\). For the phase problem, the so-called “phaseless” approximation [ZK03] is used. At each step, walkers are individually projected onto an evolving line in the complex plane defined by \(\langle \Psi_\mathrm{T} | \Phi^s_n \rangle\) by multiplying the walker weights by \(\max\{0, \cos(\Delta \theta)\}\), where \(\Delta \theta = \arg\left[\frac{\langle\Psi_\mathrm{T}|\hat{B}(\mathbf{\sigma} - \bar{\mathbf{\sigma}})|\Phi^s_n \rangle}{\langle\Psi_\mathrm{T}|\Phi^s_n \rangle}\right]\). For the sign problem, the “constrained path” approximation [ZCG97] is used, where walkers are eliminated if \(\langle \Psi_\mathrm{T} | \Phi^s_n \rangle < 0\). This can be viewed as a special case of the phaseless approximation. In both cases, the constraint introduces a bias that can be reduced based on the quality of the trial wavefunction.

Observables#

Observables are evaluated stochastically as

(22)#\[\langle \hat{O} \rangle_\mathrm{Mixed} = \frac{1}{\sum_k W_{n,k}} \sum_k W_{n,k} \frac{\langle \Psi_\mathrm{T} | \hat{O} | \Phi_{n,k} \rangle }{\langle \Psi_\mathrm{T} | \Phi_{n,k} \rangle}\]

where \(n\) is the projection step index, \(| \Phi_{n,k} \rangle\) are Slater determinant random walkers with weight \(W_{n,k}\) (from importance sampling), and \(| \Psi_\mathrm{T} \rangle\) is the trial wavefunction. The mixed estimator provides an unbiased estimate for any observable which commutes with the Hamiltonian. For the AFQMC energy, \(\hat{O} = \hat{H}\), and trivially \([ \hat{O}, \hat{H} ] = 0\); therefore, the mixed estimator is an unbiased choice for the total ground state energy.

On the other hand, it is often interesting to compute observables which do not commute with the Hamiltonian. For example, the one-body reduced density matrix, from which many properties can be computed, does not commute with the Hamiltonian. In these cases, a back-propagated estimator is necessary to mitigate the bias that would be incurred by using a mixed-estimator. The back-propagated estimator has the form,

(23)#\[\langle \hat{O} \rangle_\mathrm{BP} = \frac{1}{\sum_k W_{s+m,k}} \sum_k W_{s+m,k} \frac{\langle \tilde{\Phi}_{m,k} | \hat{O} | \Phi_{s,k} \rangle }{\langle \tilde{\Phi}_{m,k} |\Phi_{s,k}\rangle}\]

where \(| \Phi_{s,k} \rangle\) are the usual forward-projected Slater determinant random walkers, and \(| \tilde{\Phi}_{m,k} \rangle\) are the back-propagated walkers given by,

\[| \tilde{\Phi}_{m,k} \rangle = \hat{B}^\dagger( (x - \bar{x})_{s,k} ) ... \hat{B}^\dagger( (x - \bar{x})_{s+m-1,k} ) | \Psi_\mathrm{T} \rangle.\]

The index \(s\) corresponds to the current forward projection step, and \(m\) is the back-propagated step index. We note that each random walker has a corresponding back-propagated partner which share the same path in auxiliary-field space. It is important to note that the two Slater determinants are distinct from each other since the propagator \(B(x - \bar{x})\) is applied differently to each. i.e. the propagator is applied to different Slater determinants when moving in the forward-, and backward-directions.

Glossary of Symbols and Notation#

Indices#

  • \(i,j,k,l\) are orbital basis indices. In quantum chemistry or in solids, these typically refer to spatial orbitals \(\phi_i (\vec{r})\). In lattice models, these refer to a combined lattice site and band index.

  • \(\sigma, \sigma', ...\) are spin indices. We use additional “primes” to introduce distinct spin indices.

  • \(n\) is the walker index for sums over random walkers.

  • \(s\) is the AFQMC projection time step index.

Hamiltonian#

  • \(H_0\), the constant energy contributions to the Hamiltonian. For example, in ab initio calculations performed in the Born-Oppenheimer approximation, \(H_0\) includes the nuclear repulsion energy.

  • \(\hat{H}_1\), a generic “one-body” Hamiltonian consisting of all one-body terms in the Hamiltonian with the form \(\hat{H}_1 = \sum_{ij} h_{ij} \hat{c}^\dagger_i \hat{c}_j\)

  • \(\hat{H}_2\), a generic “two-body” interaction consisting of all two-body terms in the Hamiltonian with the form \(\hat{H}_2 = \sum_{ijkl} v_{ijkl} \hat{c}^\dagger_i \hat{c}^\dagger_j \hat{c}_k \hat{c}_l\)

Wavefunctions#

  • \(|\Psi\rangle\) , a generic many-body wavefunction, typically (but not necessarily!) consisting of a linear combination of Slater determinants.

  • \(|\Psi_I\rangle\) , the initial wavefunction

  • \(| \Phi \rangle\) , a single Slater determinant. Subscript and/or superscript indices may also apply.

AFQMC parameters#

  • \(\beta\), total projection time of the AFQMC imaginary-time projection

  • \(\tau\), imaginary-time projection step size

bibliography#

[ADVFerre+09]

Francesco Aquilante, Luca De Vico, Nicolas Ferré, Giovanni Ghigo, Per-åke Malmqvist, Pavel Neogrády, Thomas Bondo Pedersen, Michal Pitoňák, Markus Reiher, Björn O. Roos, Luis Serrano-Andrés, Miroslav Urban, Valera Veryazov, and Roland Lindh. MOLCAS 7: The Next Generation. J. Comput. Chem., 31(1):224, 2009. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/jcc.21318.

[BL77]

Nelson H. F. Beebe and Jan Linderberg. Simplifications in the generation and transformation of two-electron integrals in molecular calculations. Int. J. Quantum Chem., 12(4):683, 1977. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/qua.560120408.

[KdMerasP03]

Henrik Koch, Alfredo Sánchez de Merás, and Thomas Bondo Pedersen. Reduced scaling in electronic structure calculations using Cholesky decompositions. The Journal of Chemical Physics, 118(21):9481, 2003. URL: https://doi.org/10.1063/1.1578621.

[MZM20]

Fionn D. Malone, Shuai Zhang, and Miguel A. Morales. Accelerating auxiliary-field quantum monte carlo simulations of solids with graphical processing units. Journal of Chemical Theory and Computation, 16(7):4286–4297, 2020. URL: https://doi.org/10.1021/acs.jctc.0c00262, doi:10.1021/acs.jctc.0c00262.

[MZ18]

Mario Motta and Shiwei Zhang. Ab initio computations of molecular systems by the auxiliary-field quantum Monte Carlo method. WIREs Computational Molecular Science, 8(5):e1364, 2018. doi:10.1002/wcms.1364.

[MZC19]

Mario Motta, Shiwei Zhang, and Garnet Kin-Lic Chan. Hamiltonian symmetries in auxiliary-field quantum Monte Carlo calculations for electronic structure. Physical Review B, 100:045127, July 2019. URL: https://link.aps.org/doi/10.1103/PhysRevB.100.045127, doi:10.1103/PhysRevB.100.045127.

[PKVZ11]

Wirawan Purwanto, Henry Krakauer, Yudistira Virgus, and Shiwei Zhang. Assessing weak hydrogen binding on Ca+ centers: An accurate many-body study with large basis sets. The Journal of Chemical Physics, 135(16):164105, 2011. URL: https://doi.org/10.1063/1.3654002, doi:10.1063/1.3654002.

[PZ04]

Wirawan Purwanto and Shiwei Zhang. Quantum monte carlo method for the ground state of many-boson systems. Phys. Rev. E, 70:056702, November 2004. doi:10.1103/PhysRevE.70.056702.

[PZK13]

Wirawan Purwanto, Shiwei Zhang, and Henry Krakauer. Frozen-Orbital and Downfolding Calculations with Auxiliary-Field Quantum Monte Carlo. Journal of Chemical Theory and Computation, 9(11):4825–4833, 2013. URL: https://doi.org/10.1021/ct4006486.

[Suz76]

M. Suzuki. Generalized trotter's formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems. Communications in Mathematical Physics, 51:183–190, 1976. URL: https://link.springer.com/article/10.1007/BF01609348, doi:https://doi.org/10.1007/BF01609348.

[Tho60]

D.J. Thouless. Stability conditions and nuclear rotations in the hartree-fock theory. Nuclear Physics, 21:225–232, 1960. URL: https://www.sciencedirect.com/science/article/pii/0029558260900481, doi:https://doi.org/10.1016/0029-5582(60)90048-1.

[Tro59]

H. F. Trotter. On the product of semi-groups of operators. Proceedings of the American Mathematical Society, 10:545–551, 1959. URL: https://www.ams.org/journals/proc/1959-010-04/S0002-9939-1959-0108732-6/, doi:https://doi.org/10.1090/S0002-9939-1959-0108732-6.

[Zha13]

Shiwei Zhang. Auxiliary-field quantum monte carlo for correlated electron systems. Modeling and Simulation, 3:, 2013. URL: http://hdl.handle.net/2128/5389, doi:.

[ZCG97]

Shiwei Zhang, J. Carlson, and J. E. Gubernatis. Constrained path Monte Carlo method for fermion ground states. Phys. Rev. B, 55(12):7464–7477, March 1997. doi:10.1103/PhysRevB.55.7464.

[ZK03] (1,2,3)

Shiwei Zhang and Henry Krakauer. Quantum monte carlo method using phase-free random walks with slater determinants. Phys. Rev. Lett., 90:136401, April 2003. doi:10.1103/PhysRevLett.90.136401.