Hamiltonian symmetries in auxiliary-field quantum Monte Carlo calculations for electronic structure

TL;DR

This paper introduces a method to incorporate Hamiltonian symmetries into auxiliary-field quantum Monte Carlo calculations, significantly reducing computational costs especially for crystalline solids with many k points.

Contribution

The authors develop a formalism to include Abelian symmetries in AFQMC, leading to computational cost reductions and enabling more efficient simulations of large systems.

Findings

Cost reduction by factor of N_k^{-1} in AFQMC steps

Effective for molecular and crystalline systems

Enables calculations approaching the thermodynamic limit

Abstract

We describe how to incorporate symmetries of the Hamiltonian into auxiliary-field quantum Monte Carlo calculations (AFQMC). Focusing on the case of Abelian symmetries, we show that the computational cost of most steps of an AFQMC calculation is reduced by $N_k^{-1}$, where $N_k$ is the number of irreducible representations of the symmetry group. We apply the formalism to a molecular system as well as to several crystalline solids. In the latter case, the lattice translational group provides increasing savings as the number of k points is increased, which is important in enabling calculations that approach the thermodynamic limit. The extension to non-Abelian symmetries is briefly discussed.

Tables

Table 1. Table 1: Comparison of the computational cost of the main operations in an AFQMC calculation, without (left) and with the use of symmetry (right). In most cases, a cost reduction by a factor of N k − 1 superscript subscript 𝑁 𝑘 1 N_{k}^{-1} is seen. The three terms under “storage of integrals” refer to the one- and two-body parts of H ^ ^ 𝐻 \hat{H} and the trial wavefunction respectively. The two terms under “overlap, Θ Θ \Theta matrix” refer to the construction and inversion of the overlap matrix and construction of the Θ Θ \Theta matrix respectively. The two terms under “propagation (potential)” refer to construction and application of the 𝒜 𝒜 \mathcal{A} matrix. The two terms under “local energy calculation” refer to one- and two-body parts of the local energy.

operation	without symmetry	with symmetry	savings
storage of integrals	$M^2+M^2{N}_{\gamma }+MN$	$m^2{N}_k+m^2{n}_{\gamma }{N}_k^2+mn{N}_k$	$N_k^{-1}$
no. of auxiliary fields	$N_{\gamma }$	$n_{\gamma }{N}_k$	none
force bias calculation	$MN{N}_{\gamma }$	$mn{n}_{\gamma }{N}_k$	$N_k^{-1}$
overlap, $\Theta$ matrix	$M{N}^2+N^3+M^{2}N$	$m{N}^2+N^3+M^{2}N$	$N_k^{-1}$ (partial)
propagation (kinetic)	$M^{2}N$	$m^{2}N{N}_k$	$N_k^{-1}$
propagation (potential)	$M^2{N}_{\gamma }+M^{2}N$	$m^2{n}_{\gamma }{N}_k+M^{2}N$	$N_k^{-1}$ (partial)
local energy calculation	$NM+M{N}^2{N}_{\gamma }$	$mn{N}_k+m{n}^2{n}_{\gamma }{N}_k^3$	$N_k^{-1}$

Equations

$$\begin{array}{rl}\hat{H}& =\sum _{pq}{h}_{pq}{\hat{a}}_p^{†}{\hat{a}}_q^{}+\frac{1}{2}\sum _{pqrs}(pr|qs){\hat{a}}_p^{†}{\hat{a}}_q^{†}{\hat{a}}_s^{}{\hat{a}}_r^{},\\ h_{pq}& =\int d\mathbf{r}{\phi }_p^{\ast }(\mathbf{r})\left(-\frac{1}{2}{\nabla }^2-\sum _a\frac{Z_a}{|\mathbf{r}-{\mathbf{R}}_a|}\right){\phi }_q^{}(\mathbf{r}),\\ (pr|qs)& =\int d\mathbf{r}d{\mathbf{r}}^{\prime }{\phi }_p^{\ast }(\mathbf{r}){\phi }_r^{}\left(\mathbf{r}\right)\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{\phi }_q^{\ast }({\mathbf{r}}^{\prime }){\phi }_s^{}\left({\mathbf{r}}^{\prime }\right),\end{array}$$

$$[\hat{s},\hat{H}]=0.$$

$$\begin{array}{rl}\hat{\Gamma }(s)|{\phi }_p⟩& =\sum _r\Gamma (s{)}_{pr}|{\phi }_r⟩,\end{array}$$

$$\mathcal{S}\simeq {\mathbb{Z}}_{n_0}\times \cdots \times {\mathbb{Z}}_{n_{r-1}}\equiv {\mathbb{Z}}_{\mathcal{S}},\prod _{i=0}^{r-1}{n}_i=|\mathcal{S}|,$$

$$\mathbf{s}+\mathbf{t}=(s_0+t_0\text{mbox}mod{n}_0\dots s_{r-1}+t_{r-1}\text{mbox}mod{n}_{r-1}).$$

$${\hat{\Pi }}_{\mathbf{k}}=\sum _{\mathbf{s}}\frac{e^{-2\pi i\mathbf{k}\cdot \mathbf{s}}}{|\mathcal{S}|}\hat{\Gamma }(\mathbf{s}),\mathbf{k}\in {\mathbb{Z}}_{\mathcal{S}},\mathbf{k}\cdot \mathbf{s}=\sum _{i=0}^{r-1}\frac{k_i{s}_i}{n_i},$$

$${\hat{\Pi }}_{\mathbf{k}}{\hat{\Pi }}_{{\mathbf{k}}^{\prime }}={\delta }_{\mathbf{k}{\mathbf{k}}^{\prime }}{\hat{\Pi }}_{\mathbf{k}},\sum _{\mathbf{k}}{\hat{\Pi }}_{\mathbf{k}}=\hat{\mathbb{I}}.$$

$${\hat{\Pi }}_{\mathbf{k}}|{\tilde{\phi }}_{p{\mathbf{k}}_p}⟩={\delta }_{\mathbf{k},{\mathbf{k}}_p}|{\tilde{\phi }}_{p{\mathbf{k}}_p}⟩.$$

$$\begin{array}{rl}& \hat{H}=E_0+\sum _{\begin{array}{c}\mathbf{k}\\ pq\end{array}}{h}_{pq}(\mathbf{k}){\hat{a}}_{p\mathbf{k}}^{†}{\hat{a}}_{q\mathbf{k}}^{}\\ +& \sum _{\begin{array}{c}{\mathbf{k}}_p{\mathbf{k}}_r{\mathbf{k}}_q{\mathbf{k}}_s\\ prqs\end{array}}^{\ast }\frac{(p{\mathbf{k}}_p,r{\mathbf{k}}_r|q{\mathbf{k}}_q,s{\mathbf{k}}_s)}{2}{\hat{a}}_{p{\mathbf{k}}_p}^{†}{\hat{a}}_{q{\mathbf{k}}_q}^{†}{\hat{a}}_{s{\mathbf{k}}_s}^{}{\hat{a}}_{r{\mathbf{k}}_r}^{}\end{array}$$

$$\begin{array}{rl}\hat{H}& =E_0+\sum _{\begin{array}{c}\mathbf{k}\\ pq\end{array}}{h}_{pq}(\mathbf{k}){\hat{a}}_{p\mathbf{k}}^{†}{\hat{a}}_{q\mathbf{k}}^{}\\ & +\sum _{\begin{array}{c}\mathbf{Q}{\mathbf{k}}_r{\mathbf{k}}_s\\ prqs\end{array}}\frac{(p{\mathbf{k}}_r+\mathbf{Q},r{\mathbf{k}}_r|q{\mathbf{k}}_s-\mathbf{Q},s{\mathbf{k}}_s)}{2}\\ & {\hat{a}}_{p{\mathbf{k}}_r+\mathbf{Q}}^{†}{\hat{a}}_{q{\mathbf{k}}_s-\mathbf{Q}}^{†}{\hat{a}}_{s{\mathbf{k}}_s}^{}{\hat{a}}_{r{\mathbf{k}}_r}^{}\end{array}$$

$$|{\Psi }_0⟩\propto \underset{n\to \infty }{\lim }{e}^{-n\Delta \tau (\hat{H}-E_0)}|{\Phi }_I⟩,$$

$$\hat{H}-E_0={\hat{H}}_1-\frac{1}{2}\sum _{\gamma }{\hat{v}}_{\gamma }^2,$$

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

$$\hat{B}(\mathbf{x})=e^{-\frac{\Delta \tau }{2}{\hat{H}}_1}{e}^{\sqrt{\Delta \tau }{\sum }_{\gamma }{x}_{\gamma }{\hat{v}}_{\gamma }}{e}^{-\frac{\Delta \tau }{2}{\hat{H}}_1}$$

$$\begin{array}{rl}|{\Phi }_n⟩& =e^{-n\Delta \tau (\hat{H}-E_0)}|{\Phi }_I⟩=\int \prod _{l=0}^{n-1}d{\mathbf{x}}_{l}p({\mathbf{x}}_l)\hat{B}({\mathbf{x}}_l)|{\Phi }_I⟩.\end{array}$$

$$\begin{array}{rl}& \frac{⟨{\Psi }_T|\hat{H}|{\Phi }_n⟩}{⟨{\Psi }_T|{\Phi }_n⟩}=\\ =& \frac{\int {\prod }_{l=0}^{n-1}d{\mathbf{x}}_{l}p({\mathbf{x}}_l)⟨{\Psi }_T|\hat{H}|{\prod }_{l=0}^{n-1}\hat{B}({\mathbf{x}}_l)|{\Phi }_I⟩}{\int {\prod }_{l=0}^{n-1}d{\mathbf{x}}_{l}p({\mathbf{x}}_l)⟨{\Psi }_T|{\prod }_{l=0}^{n-1}\hat{B}({\mathbf{x}}_l)|{\Phi }_I⟩}=\\ =& \frac{\int d\mathbf{X}p(\mathbf{X})W(\mathbf{X}){\mathcal{E}}_{loc}({\Phi }_n(\mathbf{X}))}{\int d\mathbf{X}p(\mathbf{X})W(\mathbf{X})},\end{array}$$

$$\begin{array}{rl}W(\mathbf{X})& =⟨{\Psi }_T|\prod _{l=0}^{n-1}\hat{B}({\mathbf{x}}_l)|{\Psi }_I⟩\equiv ⟨{\Psi }_T⟩{\Phi }_n(\mathbf{X}),\\ {\mathcal{E}}_{loc}({\Phi }_n(\mathbf{X}))& =\frac{⟨{\Psi }_T|\hat{H}|{\Phi }_n(\mathbf{X})⟩}{⟨{\Psi }_T|{\Phi }_n(\mathbf{X})⟩},\end{array}$$

$$\frac{⟨{\Psi }_T|\hat{H}|{\Psi }_0⟩}{⟨{\Psi }_T|{\Psi }_0⟩}\simeq \frac{{\sum }_i{W}_i{\mathcal{E}}_{loc}({\Phi }_i)}{{\sum }_i{W}_i}.$$

$$\hat{H}-E_0={\hat{H}}_1^{\prime }-\frac{1}{2}\sum _{\gamma }{\left({\hat{v}}_{\gamma }-⟨{\hat{v}}_{\gamma }{⟩}_T\right)}^2,$$

$$\begin{split}&\hat{B}({\bf{x}})\to\hat{B}^{\prime}({\bf{x}}-\overline{{\bf{x}}})=e^{-\frac{\Delta\tau}{2}\hat{H}^{\prime}_{1}}e^{\sqrt{\Delta\tau}\sum_{\gamma}(x_{\gamma}-\overline{x}_{\gamma})\hat{v}^{\prime}_{\gamma}}e^{-\frac{\Delta\tau}{2}\hat{H}^{\prime}_{1}}\,\,,\\ &\prod_{l=0}^{n-1}\hat{B}({\bf{x}}_{l})\ket{\Phi_{I}}\to\prod_{l=0}^{n-1}\hat{B}^{\prime}({\bf{x}}_{l}-\overline{{\bf{x}}}_{l})\ket{\Phi_{I}}\equiv\ket{\Phi^{\prime}_{n}({\bf{X}},\overline{{\bf{X}}})}\,\,,\\ &W^{\prime}({\bf{X}})\to\prod_{l=0}^{n-1}I\Big{(}{\bf{x}}_{l},\overline{{\bf{x}}}_{l};\Phi^{\prime}_{l}({\bf{X}},\overline{{\bf{X}}})\Big{)}\,\,,\\ &I\big{(}{\bf{x}},\overline{{\bf{x}}};\Phi\big{)}=\frac{p({\bf{x}}-\overline{{\bf{x}}})}{p({\bf{x}})}\frac{\langle\Psi_{T}|\hat{B}^{\prime}({\bf{x}}_{l}-\overline{{\bf{x}}}_{l})\ket{\Phi}}{\langle\Psi_{T}\ket{\Phi}}\,\,,\end{split}$$

$${\overline{\mathbf{x}}}_{\gamma }=-\sqrt{\Delta \tau }\frac{⟨{\Psi }_T|{\hat{v}}_{\gamma }^{\prime }|\Phi ⟩}{⟨{\Psi }_T|\Phi ⟩}\equiv -\sqrt{\Delta \tau }⟨{\hat{v}}_{\gamma }^{\prime }⟩$$

$$\begin{split}I\big{(}{\bf{x}},\overline{{\bf{x}}};\Phi\big{)}&\simeq e^{-\Delta\tau\,\mbox{Re}\,\left(\mathcal{E}_{loc}(\Phi)-E_{0}\right)}\,\max\left(0,\cos(\Delta\theta)\right)\quad,\\ \Delta\theta&=\mbox{Arg}\,\frac{\langle\Psi_{T}|\hat{B}^{\prime}({\bf{x}}-\overline{{\bf{x}}})\ket{\Phi}}{\langle\Psi_{T}\ket{\Phi}}\quad.\end{split}$$

$$\frac{(p{\mathbf{k}}_r+\mathbf{Q},r{\mathbf{k}}_r|q{\mathbf{k}}_s-\mathbf{Q},s{\mathbf{k}}_s)}{2}=\sum _{\gamma }{L}_{p{\mathbf{k}}_r+\mathbf{Q},r{\mathbf{k}}_r}^{\gamma ,\mathbf{Q}}{L}_{q{\mathbf{k}}_s-\mathbf{Q},s{\mathbf{k}}_s}^{\gamma ,-\mathbf{Q}},$$

$$\begin{array}{rl}\hat{H}-E_0& =\sum _{\begin{array}{c}\mathbf{k}\\ pq\end{array}}{\tilde{h}}_{pq}(\mathbf{k}){\hat{a}}_{p\mathbf{k}}^{†}{\hat{a}}_{q\mathbf{k}}^{}-\frac{1}{2}\left[\sum _{\gamma \mathbf{Q}}{\left(\frac{i{\hat{L}}_{\gamma ,\mathbf{Q}}+i{\hat{L}}_{\gamma ,-\mathbf{Q}}}{\sqrt{2}}\right)}^2+{\left(\frac{{\hat{L}}_{\gamma ,\mathbf{Q}}-{\hat{L}}_{\gamma ,-\mathbf{Q}}}{\sqrt{2}}\right)}^2\right],\end{array}$$

$${\hat{L}}_{\gamma ,\mathbf{Q}}=\sum _{\begin{array}{c}{\mathbf{k}}_r\\ rp\end{array}}{L}_{p{\mathbf{k}}_r+\mathbf{Q},r{\mathbf{k}}_r}^{\gamma ,\mathbf{Q}}{\hat{a}}_{p{\mathbf{k}}_r+\mathbf{Q}}^{†}{\hat{a}}_{r{\mathbf{k}}_r}^{}.$$

$$\begin{split}\ket{\Psi_{T}}&=\prod_{i{\bf{k}}_{i}}\hat{a}^{\dagger}_{\psi_{i{\bf{k}}_{i}}}\ket{\emptyset}\quad,\\ \hat{a}^{\dagger}_{\psi_{i{\bf{k}}_{i}}}&=\sum_{r}\Big{(}{\Psi_{T}}({\bf{k}}_{i})\Big{)}_{ri}\,\hat{a}^{\dagger}_{r{\bf{k}}_{r}}\quad,\\ \end{split}$$

$$\begin{array}{rl}{\rho }_{r{\mathbf{k}}_r,p{\mathbf{k}}_p}& =⟨{\Psi }_T⟩{\hat{a}}_{p{\mathbf{k}}_p}^{†}{\hat{a}}_{r{\mathbf{k}}_r}^{}{\Psi }_T\\ & ={\delta }_{{\mathbf{k}}_p{\mathbf{k}}_r}{\left[{\Psi }_T({\mathbf{k}}_r){\Psi }_T^{†}({\mathbf{k}}_r)\right]}_{rp}\end{array}$$

$$\begin{array}{r}⟨{\Psi }_T|{\hat{L}}_{\gamma ,\mathbf{Q}}|{\Psi }_T⟩={\delta }_{\mathbf{Q},0}{\ell }_{\gamma }.\end{array}$$

$${\hat{L}}_{\gamma ,\mathbf{Q}}^{\prime }={\hat{L}}_{\gamma ,\mathbf{Q}}-\frac{{\delta }_{\mathbf{Q},0}{\ell }_{\gamma }}{N}\hat{N},$$

$$\begin{array}{rl}{\hat{v}}_{\gamma \mathbf{Q},1}^{\prime }& =\frac{i{\hat{L}}_{\gamma ,\mathbf{Q}}^{\prime }+i{\hat{L}}_{\gamma ,-\mathbf{Q}}^{\prime }}{\sqrt{2}}\\ {\hat{v}}_{\gamma \mathbf{Q},2}^{\prime }& =\frac{{\hat{L}}_{\gamma ,\mathbf{Q}}^{\prime }-{\hat{L}}_{\gamma ,-\mathbf{Q}}^{\prime }}{\sqrt{2}}\end{array}$$

Full text

I Introduction

The basic task of electronic structure (ES) theory is to solve the time-independent Schrödinger equation $\hat{H}\ket{\Psi_{\mu}}=E_{\mu}\ket{\Psi_{\mu}}$, to determine the eigenvalues $E_{\mu}$ and eigenstates $\ket{\Psi_{\mu}}$ of a Hamiltonian $\hat{H}$. When relativistic effects and nuclear motion are neglected, the second-quantized Born-Oppenheimer Hamiltonian operator has the form Born and Oppenheimer (1927); Ziman (1965); Szabo and Ostlund (1989),

[TABLE]

expressed here in atomic units, using a basis of $M$ spin-orbitals $\{\varphi_{p}\}_{p=1}^{M}$ for the one-electron Hilbert space. In Eq. (1), $a$ labels nuclei with positions ${\bf R}_{a}$, and ${\bf r},{\bf r}^{\prime}$ denote electronic coordinates.

Many ES methods take advantage of Hamiltonian symmetries Dupuis and King (1977); Szabo and Ostlund (1989); Stanton et al. (1991); Chan and Head-Gordon (2002); McClain et al. (2017); Sun et al. to improve the efficiency of calculations. In the present work, we give an instructional account of how to incorporate symmetry into the auxiliary-field quantum Monte Carlo (AFQMC) method Blankenbecler et al. (1981); Sugiyama and Koonin (1986); Zhang et al. (1997); Rom et al. (1998); Baer et al. (1998); Zhang and Krakauer (2003); Al-Saidi et al. (2006); Zhang (2013); Motta and Zhang (2018a) through symmetry adapted orbitals Dupuis and King (1977); Szabo and Ostlund (1989). We illustrate the formalism through the use of reflection symmetry in molecules and lattice translational symmetry in periodic solids. In the case of lattice translational symmetry, the symmetry orbitals become crystalline orbitals, and thus the AFQMC calculation is carried out in a similar framework to other recent implementations of quantum chemical many-body methods in crystals Evarestov (2013); Dovesi et al. (2014); Booth et al. (2016); McClain et al. (2017); Sun et al. (2017, ); Gruber et al. (2018). This use of translational symmetry should be distinguished from twist averaging  Lin et al. (2001); Ma et al. (2015); Zhang et al. (2018); Malone et al. (2019), as incorporating larger translational groups reduces many-body size effects, rather than only the one-body size effects captured by twist averaging.

The remainder of the paper is structured as follows. In Section II we briefly review the connection between Hamiltonian symmetries, matrix element sparsity, and the AFQMC formalism. We then describe in detail how to use Hamiltonian symmetries to decrease the cost of the various operations involved in an AFQMC calculation. In Section IV we apply the formalism to a test molecule, using reflection symmetries, and crystalline solids, with increasing sizes of the translational group. Conclusions are drawn in Section V. Further implementation details are provided in the Appendices.

II Background

A transformation operator $\hat{s}$ is a symmetry operator for a Hamiltonian $\hat{H}$ if the latter is invariant under the transformation $\hat{s}$,

[TABLE]

The set of operators $\hat{s}$ such that (2) holds, forms a group $\mathcal{S}$ under composition, termed the symmetry group of $\hat{H}$. Here, we will focus on Abelian symmetries, i.e. we will assume $[\hat{s}_{1},\hat{s}_{2}]=0$ for all $\hat{s}_{1},\hat{s}_{2}\in\mathcal{S}$.

II.1 Hamiltonian symmetries and sparsity

In the one-electron Hilbert space, the action of symmetry transformations is captured by one-body operators

[TABLE]

where $\Gamma(s)_{pr}$ is a $M$-dimensional matrix representation of $\mathcal{S}$.

To make the abstract group transformations $\hat{s}$ more concrete, and amenable to storage and numerical manipulation, it is useful to employ the structure theorem for finitely generated Abelian groups Hungerford (1980); Dummit and Foote (2004); Rudin (1962), which states that $\mathcal{S}$ is isomorphic to a direct product

[TABLE]

of cyclic groups $\mathbb{Z}_{n_{i}}$, of orders $n_{i}$ multiplying to the number $|\mathcal{S}|$ of elements of $\mathcal{S}$. In what follows, symmetries $\hat{s}\in\mathcal{S}$ will be labeled with strings ${\bf{s}}\in\mathbb{Z}_{\mathcal{S}}$, and sums of such strings will be understood to be modulo the orders $n_{i}$ of the cyclic groups,

[TABLE]

As detailed in Appendix A.1, from the properties of the discrete Fourier transform Rudin (1962); Jozsa (1998), and the relation $\hat{\Gamma}({\bf{s}})\hat{\Gamma}({\bf{t}})=\hat{\Gamma}({\bf{s}}+{\bf{t}})$, the operators

[TABLE]

form a complete set of orthogonal projectors,

[TABLE]

Applying the projectors $\hat{\Pi}_{{\bf{k}}}$ to the basis functions, and orthonormalizing the resulting vectors, yields an orthonormal basis of symmetry-adapted orbitals $\{\tilde{\varphi}_{p{\bf{k}}_{p}}\}$. Here, and in the remainder of the present work, ${\bf{k}}_{p}\in\mathbb{Z}_{\mathcal{S}}$ denotes an irreducible representation or irrep, and $p$ a group of $m_{{\bf{k}}_{p}}$ orbitals labelled by the irrep ${\bf{k}}_{p}$. The numbers $m_{{\bf{k}}_{p}}$ sum to $M$, and

[TABLE]

The number of irreps will be denoted $N_{k}=|\mathcal{S}|$. Since the orbitals $\tilde{\varphi}_{p{\bf{k}}_{p}}$ are eigenfunctions of the projectors $\hat{\Pi}_{{\bf{k}}}$, and the latter commute with the one- and two-body parts of the Hamiltonian, the matrix elements of $\hat{H}$ are sparse, as revealed by the expression

[TABLE]

where the $*$ over the summation denotes the constraint ${\bf{k}}_{p}+{\bf{k}}_{q}={\bf{k}}_{r}+{\bf{k}}_{s}$ (with equality holding in the modular arithmetic of $\mathbb{Z}_{\mathcal{S}}$). (9) can be rewritten in terms of the transfer parameter ${\bf{Q}}={\bf{k}}_{p}-{\bf{k}}_{r}={\bf{k}}_{s}-{\bf{k}}_{q}$, leading to

[TABLE]

The structure of the Hamiltonian operator is illustrated in Figure 1. In the forthcoming sections, after providing a brief description of the AFQMC method, we will discuss in detail how the form of the Hamiltonian in (10) leads to savings in many of the operations in the method.

II.2 The AFQMC method

The AFQMC method Zhang and Krakauer (2003); Purwanto and Zhang (2004); Zhang (2013); Motta and Zhang (2018a) expresses the many-body ground state $\Psi_{0}$ of a Hamiltonian $\hat{H}$ through an imaginary time evolution,

[TABLE]

where $\Delta\tau$ (the time step) is chosen small and the initial state $\Phi_{I}$, which should not be orthogonal to $\Psi_{0}$, is often a single Slater determinant. To sample the many-body propagator, we rewrite the Hamiltonian as

[TABLE]

where $\hat{H}_{1}$, $\hat{v}_{\gamma}$ are one-body operators. Then, using a Hubbard-Stratonovich transformation Hubbard (1959); Stratonovich (1958), the short-imaginary-time propagator is

[TABLE]

where

[TABLE]

is a one-body propagator that is a function of the multi-dimensional vector ${\bf{x}}$, and $p({\bf{x}})$ is the standard normal probability distribution. AFQMC thus represents the many-body wave function as a superposition of non-orthogonal Slater determinants,

[TABLE]

The ground-state expectation value of $\hat{H}$ is obtained as

[TABLE]

where $\Psi_{T}$ is a second many-body state, called the trial wavefunction. In Eq. (LABEL:eq:esta), the overlap and local energy,

[TABLE]

are defined on a path ${\bf{X}}=({\bf{x}}_{n-1}\dots\bf{x}_{0})$ of auxiliary fields at each time slice up to $n-1$. The expectation value is computed over a collection of Monte Carlo (MC) samples (labeled by $i$) as

[TABLE]

and the stochastically sampled determinants $\Phi_{i}$ are called walkers in the AFQMC literature.

Because the propagator $\hat{B}({\bf{x}})$ contains stochastically fluctuating fields, the MC sampling will lead to complex overlaps $W_{i}$, which causes the variance of this estimator to grow exponentially with the number of time-steps $n$. This phase problem can be controlled by an approximate gauge condition known as the phaseless approximation Zhang and Krakauer (2003); Purwanto and Zhang (2004); Zhang (2013); Motta and Zhang (2018a), which we summarize below:

mean-field background subtraction:

the expectation values $\langle\hat{v}_{\gamma}\rangle_{T}=\langle\Psi_{T}|\hat{v}_{\gamma}|\Psi_{T}\rangle$ are computed and the Hamiltonian is rewritten as

[TABLE] 2. 2.

importance sampling transformation:

the Hubbard-Stratonovich is defined up to a shift ${\bf{x}}\to{\bf{x}}-\overline{{\bf{x}}}$, $\overline{{\bf{x}}}\in\mathbb{C}$, and this additional freedom is exploited to rewrite the estimator (LABEL:eq:esta) with the replacements

[TABLE]

where $\hat{v}^{\prime}_{\gamma}=\hat{v}_{\gamma}-\langle\hat{v}_{\gamma}\rangle_{T}$. One makes the choice

[TABLE]

to minimize fluctuations in the importance function $I({\bf{x}},\overline{{\bf{x}}};\Phi)$ to leading order in $\Delta\tau$. 3. 3.

Real local energy and cosine approximations:

the importance function is approximated as

[TABLE]

Steps 1, 2 are simply re-parametrizations which reduce fluctuations in the estimators of physical properties when $\overline{\bf{x}}_{\gamma}$ is real. When it is complex, these transformations ensure that the gauge variation is minimized, enabling to remove the sign problem in step 3. The bias resulting from the approximations in step 3 can be reduced by improving the trial wavefunction. AFQMC has been successfully applied to multiple lattice models of correlated electrons LeBlanc et al. (2015); Qin et al. (2016); Zheng et al. (2017) and real materials Purwanto et al. (2015); Motta et al. (2017); Shee et al. (2019) achieving accuracies competitive with other high-level wavefunction methods. There are many additional algorithmic improvements and extensions, for example to compute properties Motta and Zhang (2017, 2018b); Shee et al. (2017), and to improve the efficiency of the method Motta and Zhang (2018a); Motta et al. (2018).

III Method

We now illustrate how to account for Hamiltonian symmetries in the above procedure. Additional implementation details appear in the Appendices. In what follows, we consider a symmetry group composed of $N_{k}$ Abelian symmetries, with irreps labelled by ${\bf{k}}$. We use the symbols $M$,$N$,$N_{\gamma}$ to denote the total number of basis functions, particles and auxiliary fields respectively. Correspondingly, we use $m_{\bf{k}}$, $n_{\bf{k}}$ and $n_{\gamma,{\bf{k}}}$ to denote the number of basis functions, particles and auxiliary fields labelled by the irrep ${\bf{k}}$, and $m$, $n$, $n_{\gamma}$ to denote the average numbers of basis functions, particles and auxiliary fields per irrep, $m=\sum_{\bf{k}}m_{\bf{k}}/N_{k}$ etc. Indices $prqs$, $ij$, $\gamma$ run over basis functions, particles and auxiliary fields labelled by a specific irrep ${\bf{k}}$, respectively.

III.1 Hamiltonian representation

As expressed in (12), in AFQMC we must express the two-body part of the Hamiltonian as a sum of squares of one-body operators. To achieve this, one commonly relies on a density fitting (DF) Whitten (1973); Hohenstein and Sherrill (2010) or Cholesky decomposition (CD) Beebe and Linderberg (1977); Koch et al. (2003); Aquilante et al. (2010); Purwanto et al. (2011) of the electron repulsion integrals, $(pr|qs)\simeq\sum_{\gamma=1}^{N_{\gamma}}L_{pr}^{\gamma}L_{qs}^{\gamma}$ where $N_{\gamma}$ is the number of components. As detailed in Appendix A.2, in the presence of symmetries, such a decomposition becomes

[TABLE]

where, as illustrated in Figure 1, components $\gamma$ are labelled by irreps ${\bf{Q}}$, $-{\bf{Q}}$, the number of which $n_{\gamma,{\bf{Q}}}$ sums to $N_{\gamma}$. Then (as detailed in Appendix B) we interchange the creation and destruction operators in (10) and use (23) to rewrite the Hamiltonian as a sum of squares,

[TABLE]

where

[TABLE]

The total number of auxiliary fields is thus $2N_{\gamma}$, where the correspondence with the label $\gamma$ in (12) is $\gamma\to(\gamma,{\bf{Q}},1)$, $(\gamma,{\bf{Q}},2)$ and the auxiliary field operators are $\hat{v}_{\gamma,{\bf{Q}},1}=(i\hat{L}_{\gamma,{\bf{Q}}}+i\hat{L}_{\gamma,-{\bf{Q}}})/\sqrt{2}$ and $\hat{v}_{\gamma,{\bf{Q}},2}=(\hat{L}_{\gamma,{\bf{Q}}}-\hat{L}_{\gamma,-{\bf{Q}}})/\sqrt{2}$ respectively. A simple technique to reduce the number of auxiliary fields to $N_{\gamma}$ is described in Appendix B.3.

Importantly, the correction to the one-body part of the Hamiltonian resulting from the interchange of creation and destruction operators in (10) does not mix orbitals labeled by different irreps (i.e. it is block-diagonal). Additional details are given in Appendix B.1.

It is clear that the Hamiltonian integrals, such as $(p{\bf{k}}_{r}{\bf{+Q}},r{\bf{k}}_{r}|q{\bf{k}}_{s}{\bf{-Q}},s{\bf{k}}_{s})$, require $1/N_{k}$ less storage than without symmetry. The coefficients defining $\hat{L}_{\gamma,{\bf{Q}}}$ also show a $1/N_{k}$ reduction in storage compared to without symmetry.

III.2 Mean-field wavefunction and background subtraction

We assume that the trial wavefunction is a single determinant. If the single determinant wavefunction transforms as an irrep of $\mathcal{S}$, then its orbitals $\psi_{i{{\bf{k}}_{i}}}$ can also be labelled by irreps, thus

[TABLE]

where the numbers $n_{\bf{k}}$ of particles in each orbital irrep satisfy $\sum_{{\bf{k}}}n_{\bf{k}}=N$ and the coefficient matrix ${\Psi_{T}}({\bf{k}}_{i})$ is blocked by symmetry. In contrast, a generic Slater determinant (such as the walkers in an AFQMC calculation) is parametrized by a dense $M\times N$ matrix $\Phi_{r{\bf{k}}_{r},i{\bf{k}}_{i}}$. An important property of (26) is that the one-body density matrix

[TABLE]

is also block-diagonal. The mean-field expectation values of the operators $\hat{L}_{\gamma,{\bf{Q}}}$ thus read

[TABLE]

Defining the operators

[TABLE]

which by construction have zero average over $\Psi_{T}$, we can obtain operators with the mean-field background subtracted, as in (19) and detailed in Appendix B.2,

[TABLE]

In this step, the trial wavefunction requires storing $\sum_{\bf{k}}m_{\bf{k}}n_{\bf{k}}\simeq mnN_{k}$ coefficients, a reduction of $1/N_{k}$ compared to without symmetry. Similarly, the block structure of the trial wavefunction and density matrix means that computing the mean-field density matrix and subsequent background term also involves a $1/N_{k}$ reduction in the number of operations.

III.3 Overlap calculation

The overlap matrix between the trial wavefunction and a walker

[TABLE]

is $\braket{\Psi_{T}}{\Phi}=\det(\Omega)$, with

[TABLE]

Due to the symmetry block structure of $\Psi_{T}$, computing $\Omega$ requires $\mathcal{O}(mN^{2})$ operations, even though $\Omega$ is in general dense. Compared to the corresponding cost of $\mathcal{O}(MN^{2})$ in a calculation that does not exploit symmetry, this is more efficient by a factor of $1/N_{k}$.

III.4 Force bias calculation

Using symmetry, the force bias calculation involves the quantities $\langle\hat{v}^{\prime}_{\gamma}\rangle$. These are easily related to

[TABLE]

with ${{\bf{Q}}}={\bf{k}}_{p}-{\bf{k}}_{r}$, which are calculated as

[TABLE]

where ${\bf{Q}}={\bf{k}}_{p}-{\bf{k}}_{r}$ and $\Theta=\Omega^{-1}\Phi$. $\mathcal{L}$ can be precomputed, storing $mnn_{\gamma}N_{k}^{2}$ complex numbers. The computational cost to obtain the force bias with precomputation is $\mathcal{O}(mnn_{\gamma}N_{k}^{2})$. Compared with the cost of a calculation that does not exploit symmetry ($MNN_{\gamma}$ operations) this is more efficient by a factor $1/N_{k}$. Note, however, that while we can use symmetry when generating the matrix $\Omega$, as discussed in Section III.3, the cost of computing its determinant and its inverse is not reduced by symmetry, because it is in general dense. There is thus only a partial gain due to symmetry for the force bias calculation.

III.5 Walker propagation

With symmetry, the small-imaginary time propagator in (13) takes the form

[TABLE]

for operators $\hat{A}({\bf{x}})$, defined as a linear combination of $x_{{\bf{Q}}\gamma 1}\hat{v}^{\prime}_{\gamma,{\bf{Q}},1}$, $x_{{\bf{Q}}\gamma 2}\hat{v}^{\prime}_{\gamma,{\bf{Q}},2}$, as derived in Appendix B.3. Applying the exponential of $\hat{H}^{\prime}_{1}$ to a (walker) Slater determinant can be carried out efficiently. In fact, since $\hat{H}^{\prime}_{1}$ is symmetric and thus does not mix irreps, it has the form $\hat{H}^{\prime}_{1}=\sum_{{\bf{k}}pq}h^{\prime}_{pq}({\bf{k}})\,\hat{a}^{\dagger}_{p{\bf{k}}}\hat{a}^{\phantom{\dagger}}_{q{\bf{k}}}$. Thus if the walker $\ket{\Phi}$ is parametrized by the matrix $\Phi_{p{\bf{k}}_{p},i{\bf{k}}_{i}}$, its image $\ket{\Phi^{\prime}}=e^{-\frac{\Delta\tau}{2}\hat{H}_{1}^{\prime}}\ket{\Phi}$ from the one-body propagator is parametrized by the matrix

[TABLE]

This can be computed using $\mathcal{O}(m^{2}NN_{k})$ operations, as compared to $\mathcal{O}(M^{2}N)$ operations in a calculation that does not use symmetry. To propagate the interacting part of $\hat{H}$, one possible strategy is to construct the matrix $\mathcal{A}_{p{\bf{k}}_{p},r{\bf{k}}_{r}}$ associated with the operator $\hat{A}({\bf{x}})$ and then apply $e^{\mathcal{A}}$ to $\Phi$. The cost to construct $\mathcal{A}$ is $\mathcal{O}(m^{2}n_{\gamma}N_{k}^{2})$ (it is a linear combination of block sparse matrices), a reduction of $1/N_{k}$ compared to without symmetry. However, since in general $\mathcal{A}$ lacks any sparsity properties, the cost of applying its exponential to $\Phi$ requires $\mathcal{O}(M^{2}N)$ operations. Note that, in principle, sparsity could be exploited by applying individual terms of $e^{\hat{A}}$, such as $e^{\sqrt{\Delta\tau}ix_{{\bf{Q}}\gamma 1}\hat{v}^{\prime}_{\gamma,{\bf{Q}},1}}$, where such an exponential is further expanded as a power series with each term transforming as an irrep and corresponding to sparse matrix multiplication onto $\Phi$. However, this is only a savings for symmetry groups where the number of irreps far exceeds the number of terms in the power series. We have not observed savings with this second strategy in this work. Thus we only report calculations where we apply the full $\mathcal{A}$ matrix, with only a partial (i.e. limited to the generation of $\mathcal{A}$) reduction in cost in this step from symmetry.

III.6 Local energy calculation

For the one-body part of the local energy, the form of the estimator and the speedup over implementations without symmetry are very similar to the case of the force bias. Indeed, an analogous calculation leads to the expressions

[TABLE]

For the considerably more expensive two-body part, we use the generalized Wick’s theorem Wick (1950); Balian and Brezin (1969) to obtain

[TABLE]

where $ij$ are associated with particles labelled by the irreps ${\bf{k}}_{p}={\bf{k}}_{r}+{\bf{Q}}$, ${\bf{k}}_{q}={\bf{k}}_{s}-{\bf{Q}}$ respectively, and the tensor $f$ is defined as

[TABLE]

As seen, the cost of the procedure is $\mathcal{O}(mn^{2}n_{\gamma}N_{k}^{3})$ operations to generate the tensor $f$, and $\mathcal{O}(n^{2}n_{\gamma}N_{k}^{3})$ to perform the final contraction. Compared with a calculation without symmetry, this is more efficient by a factor of $1/N_{k}$.

III.7 Summary

The acceleration achieved in AFQMC due to the use of symmetries is summarized in Table 1. Most computational steps are accelerated by a factor of $1/N_{k}$, and storage is reduced by $1/N_{k}$ as well. In a standard mean-field calculation, symmetries lead to an acceleration by a factor of $1/N_{k}^{2}$, due to symmetries in both the Hamiltonian as well as the wavefunction. In AFQMC, the acceleration is limited to $1/N_{k}$ because the walkers do not transform as irreps of $\mathcal{S}$, as individual components of the Hubbard-Stratonovich such as $\hat{L}_{\gamma,{\bf{Q}}}$ all transform as different irreps of $\mathcal{S}$.

IV Results

We now present some illustrative calculations using symmetry in AFQMC calculations for a molecule and for crystalline systems. Restricted Hartree-Fock (RHF), density functional theory (DFT), Møller-Plesset perturbation theory (MP2) and coupled-cluster with singles and doubles (CCSD) calculations were performed with the PySCF package Sun et al. .

Molecular calculations were all-electron calculations using the cc-pVTZ basis Dunning (1989); Woon and Dunning (1993). The auxiliary field decomposition was performed using Cholesky decomposition, and the RHF state was used as a trial wavefunction in the AFQMC calculations.

In the crystal calculations presented below, core electrons were replaced by norm-conserving GTH Padé pseudopotentials Goedecker et al. (1996); Hartwigsen et al. (1998); McClain et al. (2017). Hamiltonian matrix elements were computed with the PySCF program Sun et al. using the GTH series of Gaussian bases Hutter et al. (2014). Gaussian density fitting was used to treat the electron-electron interaction and to obtain the auxiliary field decomposition Sun et al. (2017). RHF energies are reported using the leading finite size correction for the ${\bf{G}}=0$ contribution to the Hartree-Fock exchange (exxdiv=ewald) and total energies from other methods were obtained by adding this RHF energy to the respective correlation energies, with integrals computed omitting the ${\bf{G}}=0$ term (exxdiv=None) McClain et al. (2017); Sun et al. . The RHF state was used as a trial wavefunction in the AFQMC calculations.

IV.1 Molecular systems and point group symmetry

As a simple test-case, we first consider a molecular system, SF6, where we use the reflection group symmetry. The reflection group is isomorphic to $\mathbb{Z}_{2}^{3}$, where $3$ is the number of reflection planes, giving 8 irreps in the group.

In Fig. 2 we show the equation of state $E(R)$ of the molecule (where $R$ is the S-F bond length) using AFQMC, DFT with the B3LYP functional, RHF and CCSD. As can be seen, the AFQMC calculations performed with and without reflection symmetry (red circles and dark red crosses in Fig. 2) yield identical results to within statistical error. CCSD and AFQMC yield potential energy surfaces in good agreement with each other.

IV.2 Crystalline solids

The computational saving from symmetries is especially important in systems with a large symmetry group, and crystalline solids form one such example. Consider a crystal with a primitive cell with lattice vectors ${\bf{a}}_{0}$, ${\bf{a}}_{1}$, ${\bf{a}}_{2}$. We use translational-symmetry-adapted (crystalline) Gaussian atomic orbitals McClain et al. (2017) (AOs) as a symmetry basis. Starting from a set of Gaussian AOs in the primitive cell $\varphi_{\mu}$, these can be written as

[TABLE]

where ${\bf{i}}$ is an integer vector $(i_{0},i_{1},i_{2})$ denoting $\varphi_{\mu}$ translated from the primitive cell by lattice vector $\sum_{r=0}^{2}i_{r}{\bf{a}}_{r}$, and ${\bf{k}}$ has the form $\sum_{r=0}^{2}\frac{k_{r}}{N_{r}}{\bf{b}}_{r}$, where $k_{r}$ is an integer vector with $0\leq k_{r}<N_{r}$ and ${\bf{b}}_{r}$ are the reciprocal lattice vectors. This choice of $k_{r}$ is equivalent to sampling the Brillouin zone with a mesh of $N_{0}\times N_{1}\times N_{2}$ wave-vectors including the $\Gamma$ point (origin). Note the above basis representation spans the same Hilbert space as a $N_{0}\times N_{1}\times N_{2}$ supercell calculation with ($\Gamma$ point) periodic boundary conditions Evarestov (2013).

The Hamiltonian symmetries are the lattice translations, corresponding to integer multiples $\sum_{r=0}^{2}s_{r}{\bf{a}}_{r}$ of the ${\bf{a}}_{0}$, ${\bf{a}}_{1}$, ${\bf{a}}_{2}$ vectors. Under such translations, the basis transforms as

[TABLE]

thus, the translation group is isomorphic to $\mathbb{Z}_{N_{0}}\times\mathbb{Z}_{N_{1}}\times\mathbb{Z}_{N_{2}}$.

To demonstrate the symmetry-adapted AFQMC using the lattice translation group we first compute the equilibrium lattice constants of C diamond and Si FCC in Figure 3, using a 2$\times$2$\times$2 $k$-point mesh, at the GTH-DZV level. Here we find that AFQMC is in good agreement with CCSD using the same $k$-point mesh, and significantly improves on RHF and MP2. This trend can be seen in the potential energy surfaces in Figure 3, as well as in the corresponding equilibrium lattice constants.

Using translational symmetry, we can further consider larger symmetry groups in order to extrapolate to the thermodynamic limit (TDL). Note that increasing the size of the translational symmetry group yields the same result as a calculation with increased supercell size, but with much reduced cost.

We illustrate the extrapolation of results to the thermodynamic limit in Figure 4, using C diamond as a test system. RHF and correlation energies were computed for 2$\times$2$\times$2, 2$\times$2$\times$3, 3$\times$3$\times$3, 3$\times$3$\times$4, 4$\times$4$\times$3 and 4$\times$4$\times$4 meshes of $k$-points, using the GTH-DZV basis. In the upper panel of Figure 4 we show the equation of state $E_{TDL}(R)$ extrapolated to the thermodynamic limit (with the minimum value $E_{TDL}(R_{eq})$ subtracted) from RHF, MP2 and AFQMC.

We extrapolate RHF total energies and AFQMC, MP2 correlation energies (per cell) to the TDL using power-law Ansatz $E(N_{k})=\alpha+\beta N_{k}^{\mu}$, with $\mu=-1$ for RHF and AFQMC Kwee et al. (2008) and $\mu=-\frac{1}{3}$ for MP2 McClain et al. (2017). Extrapolation of RHF (correlation) energies is carried out using data for all but the smallest two (the smallest) $k$-point meshes. In the lower panel of Figure 4, we illustrate the extrapolation of MP2 (left) and AFQMC (right) correlation energies at the representative bondlength $R=3.6\,\mbox{\AA}$.

Fitting the TDL curves to the Morse potential Ansatz $E(R)=E_{0}+\Delta E\,(1-e^{-\alpha(R-R_{eq})})^{2}$ gives an equilibrium bondlength of $R_{eq,AFQMC}=3.575(1)\,\mbox{\AA}$. For the 2$\times$2$\times$2 supercell, the same procedure yields $R_{eq,AFQMC}=3.657(1)\,\mbox{\AA}$, thus TDL extrapolation significantly shortens the AFQMC equilibrium bondlength. For reference, the experimental bondlength is $R=3.553\,\mbox{\AA}$, corrected for zero-point vibrational effects Schimka et al. (2011); for a more faithful comparison, a larger basis set should be used Zhang et al. (2018).

In Figure 5 we carry out a similar calculation for 2D hexagonal boron nitride. RHF energies and AFQMC correlation energies (inset) were computed for 4$\times$4$\times$1 meshes of $k$-points. Total energies are shown in the upper panel, measured from the minimum value, $E_{min}=E(R_{eq})$. The AFQMC equilibrium bondlengths are $R_{eq}=1.5133(9)$, $1.4613(16)$, $1.4478(5)$, $1.4455(12)$ Å for GTH-SZV, GTH-DZV, GTH-DZVP and GTH-TZVP respectively; for reference, the reported experimental equilibrium bondlength is $R_{eq}=1.45\mbox{\AA}$ 111 https://github.com/cryos/avogadro/blob/master/

crystals/nitrides/BN.cif.

So far, we have illustrated the use of symmetry when calculating total energies and lattice constants. We now briefly show that symmetry adaptation can be used when computing arbitrary ground-state properties in AFQMC, such as the electron density, within the back-propagation algorithm Zhang et al. (1997); Purwanto and Zhang (2004); Motta and Zhang (2017).

The electron density is computed by contracting the spin-summed one-body density matrix with the basis orbitals $\varphi_{p{\bf{k}}_{p}}({\bf{x}})$, evaluated on a mesh of points ${\bf{x}}$ along the lattice plane,

[TABLE]

The one-body density matrix is evaluated using the back-propagation algorithm as

[TABLE]

Here, $|\Psi_{i}\rangle$ is a stochastically sampled Slater determinant sampled at imaginary time $n\,\Delta\tau$, $W_{i}$ its future weight at some time $(n+m)\,\Delta\tau$ and $\langle\Phi_{i}|$ is obtained back-propagating (i.e. propagating as a bra or linear functional, rather than a ket or vector) $\langle\Psi_{T}|$ along the segment of the future path of $|\Psi_{i}\rangle$, sampled during the time interval between $n\,\Delta\tau$ and $(n+m)\,\Delta\tau$ Motta and Zhang (2017).

In Figure 6, we compute the electronic density of two low-dimensional materials, 2D hexagonal BN and graphene, within the GTH-DZV basis. The electron density illustrates the different nature of the two materials: while in graphene the density is distributed uniformly around C atoms, in BN there is a net concentration of electrons around N atoms, consistent with the charge-transfer nature of the material.

IV.3 Timings

In Figure 7 we compare the timings of the standard and symmetry-adapted AFQMC implementations. Timings were performed on a cluster with Intel E5-2680, 2.4 GHz CPUs. In the various panels, the times for force bias and local energy evaluation, Hubbard-Stratonovich operator construction and walker propagation, the most expensive steps of an AFQMC calculation, are shown for standard AFQMC calculations of BN at the GTH-DZV level, using supercells of increasingly large size $N_{s}$, and symmetry-adapted calculations using $k$-point meshes of increasingly large size $N_{k}$. In a standard calculation, the local energy evaluation scales as $N_{s}^{4}$ and all the other subroutines as $N_{s}^{3}$. In a symmetry-adapted calculation, the local energy evaluation scales as $N_{k}^{3}$ and all the other subroutines scale as $N_{k}^{2}$, confirming the reduction in scaling by one power arising from symmetry adaptation. In the current implementation, the lower scaling comes at the cost of an increased prefactor, so that crossover between the two strategies occurs around $N_{s}=N_{k}\simeq 10$ for the local energy evaluation and $N_{s}=N_{k}\simeq 20$ for all other subroutines.

V Conclusions

In this work, we presented a formalism to perform AFQMC calculations that take advantage of Abelian Hamiltonian symmetries. We described how within a symmetry adapted orbital basis, the matrix elements of the Hamiltonian operator acquire block sparsity, which, when combined with a trial state that transforms as an irrep of the symmetry group and Hubbard-Stratonovich fields that also transform as irreps of the symmetry group, it is possible to reduce the cost and memory of the main steps in the AFQMC calculation by a factor of $N_{k}^{-1}$, where $N_{k}$ is the order of the group.

Extending this formalism to non-Abelian symmetries is straightforward. The only difference arises because irreps of non-Abelian groups need not be one-dimensional. Thus products of objects that transform as irreps (such as $\varphi^{*}_{p{\bf{k}}_{p}}({\bf{r}})\,\varphi_{r{\bf{k}}_{r}}({\bf{r}})$) no longer simply transform as a single irrep ${\bf{k}}_{r}-{\bf{k}}_{p}$, but correspond to a linear combination of objects, each transforming according to potentially different irreps. Nonetheless, all quantities that are block sparse in the current algorithm remain block sparse in the non-Abelian generalization, and a similar speedup of $\mathcal{O}(N_{k}^{-1})$ will be achieved as it is observed here.

As we showed in our demonstration calculations, the use of Abelian symmetries is particularly beneficial in the context of the large translational group associated with crystalline calculations. Thus we believe the present work will be particularly important in accelerating AFQMC calculations in realistic materials, and in particular, in removing finite size effects and in extrapolations to the thermodynamic limit.

VI Acknowledgments

M. M. acknowledges Qiming Sun and James McClain for assistance and discussions regarding ES calculations for crystalline solids. This work was supported by the US Department of Energy, Office of Science (via Grant No. SC0019390 to G. K.-L. C.). S. Z. acknowledges support from DOE (Grant No. DE-SC0001303). Additional software developments for Hamiltonian symmetries implemented in PySCF were supported by US NSF (Grant No. 1657286). Computations were carried out on facilities supported by the US Department of Energy, National Energy Research Scientific Computing Center (NERSC), on facilities supported by the Scientific Computing Core at the Flatiron Institute, on the Pauling cluster at the California Institute of Technology, and on the Storm and SciClone Clusters at the College of William and Mary. The Flatiron Institute is a division of the Simons Foundation.

Appendix A Additional theoretical details

A.1 Properties of the $\hat{\Pi}_{\bf{k}}$ operators

The relation

[TABLE]

readily implies that the operators $\hat{\Pi}_{\bf{k}}$ are orthogonal projectors,

[TABLE]

Completeness holds, since

[TABLE]

and the neutral element ${\bf{0}}$ of $\mathbb{Z}_{\mathcal{S}}$ is mapped onto the neutral element $\hat{\Gamma}({\bf{0}})=\hat{\mathbb{I}}$ of $\mathcal{S}$. Finally,

[TABLE]

Hamiltonian sparsity easily follows from the fact that symmetry-adapted orbitals are eigenfunctions of projection operators, and that projection operators commute with the one-body and two-body parts of the Hamiltonian. Indeed,

[TABLE]

and, since the projectors onto symmetry adapted orbitals in the two-particle Hilbert space become

[TABLE]

one has

[TABLE]

A.2 Density fitting and Cholesky decomposition

In this Section we show how the structure (23) emerges when the electron-electron interaction is treated within the density fitting (DF) or Cholesky (CD) decomposition. Within DF, the electron repulsion integral is approximated by density fitting with an auxiliary basis of atom-centered Gaussian atomic orbitals $\{\chi_{\gamma}\}_{\gamma=1}^{N_{\gamma}}$,

[TABLE]

where $S_{\gamma\delta}=\langle\chi_{\gamma}|\chi_{\delta}\rangle$. The action of the symmetry group on the auxiliary basis $\{\chi_{\gamma}\}_{\gamma}$ is captured by a family of operators

[TABLE]

so that, by following the procedure outlined in Section II.1, one can produce an orthonormal basis of symmetry-adapted auxiliary basis functions $\tilde{\chi}_{\gamma{\bf{Q}}}$ (with overlap matrix equal to the identity). The electron repulsion integral reads, in the symmetry-adapted molecular and auxiliary bases,

[TABLE]

where summation is restricted to auxiliary basis functions belonging to the irrep labelled by ${\bf{Q}}$ for the pair $(p{\bf{k}}_{r}+{\bf{Q}},r{\bf{k}}_{r})$ and by $-{\bf{Q}}$ for the pair $(q{\bf{k}}_{s}-{\bf{Q}},s{\bf{k}}_{s})$ respectively.

Performing a Cholesky decomposition of the electron repulsion integral,

[TABLE]

may not lead to the form (23), i.e. the tensor $L$ may not be sparse. The desired structure can be extracted performing a SVD of the rank-three tensor $L^{\gamma}_{q{\bf{k}}_{q},s{\bf{k}}_{s}}=\sum_{\mu}U^{\mu}_{q{\bf{k}}_{q},s{\bf{k}}_{s}}\sigma_{\mu}V^{\mu\gamma}$ After the SVD is taken, the ERI reads

[TABLE]

and the tensor $U$ is non-zero only for certain values of the index $\mu$, that depend only on the difference ${\bf{k}}_{p}-{\bf{k}}_{r}={\bf{Q}}$. Indices $\mu$ can thus be parametrized as pairs $\gamma,{\bf{Q}}$, and the ERI takes the desired form in Sec. II.1.

Appendix B Additional algorithmic details

B.1 Interaction as squares of one-body operators

Starting from (10), we interchange the creation and destruction operators,

[TABLE]

and inserting this equation in (10), obtaining

[TABLE]

with $\hat{L}_{\gamma,{\bf{Q}}}$ as in Eq. (25). To obtain a representation as a sum of squares of one-body operators, we observe that

[TABLE]

B.2 Mean-field background subtraction

The mean-field background subtraction requires replacing the operators $\hat{L}_{{\bf{Q}},\gamma}$ with $\hat{L}^{\prime}_{{\bf{Q}},\gamma}$ in (24). This leads to

[TABLE]

where $\ell$ is defined as in Eq. (28). The operator (59) has the same form as in (19) with $\hat{H}^{\prime}_{1}=\sum_{\bf{k}}h^{\prime}_{pq}({\bf{k}})\hat{a}^{\dagger}_{p{\bf{k}}}\hat{a}^{\phantom{\dagger}}_{q{\bf{k}}}$,

[TABLE]

and $\hat{v}^{\prime}_{\gamma}$ as detailed in the main text.

B.3 Reducing the number of auxiliary fields

by Lagrangian partition

In (24), the Hamiltonian was expressed as

[TABLE]

clearly leading to $2N_{\gamma}$ auxiliary fields. This is more than in a calculation that does not incorporate symmetry. To reduce the number of auxiliary fields, we partition the irreps ${\bf{Q}}$ of the symmetry group into three sets:

• •

the set $\mathcal{P}_{0}=\{{\bf{Q}}:{\bf{Q}}=-{\bf{Q}}\}$ of ${\bf{Q}}$ coinciding with their inverse

• •

any subset $\mathcal{P}_{+}\subset\mathbb{Z}_{\mathcal{S}}-\mathcal{P}_{0}$ such that, if ${\bf{Q}}\in\mathcal{P}_{+}$, then $-{\bf{Q}}\notin\mathcal{P}_{+}$

• •

$\mathcal{P}_{-}=\mathbb{Z}_{\mathcal{S}}-\mathcal{P}_{0}-\mathcal{P}_{+}$

According to Lagrange’s theorem Hungerford (1980); Dummit and Foote (2004); Rudin (1962), the set $\mathcal{P}_{0}$ contains elements other than ${\bf{0}}$ in, and only in, groups with even order $|\mathcal{S}|$. Then clearly one has

[TABLE]

Now, since $|\mathcal{P}_{0}|+2|\mathcal{P}_{+}|=|\mathcal{S}|$, the interaction part of the Hamiltonian has been reduced to a sum of $N_{\gamma}$ squares of one-body operators, the same as in a calculation that does not enforce symmetries.

With this representation of the Hamiltonian, the small-imaginary-time propagator has the form (35) with

[TABLE]

This leads immediately to the form of the matrix $\mathcal{A}$ associated with $\hat{A}({\bf{x}})$.