QMCPACK: an open source ab initio quantum Monte Carlo package for the electronic structure of atoms, molecules and solids

The main production algorithms in QMCPACK use methods based on Monte Carlo sampling of electron positions in real space to produce highly accurate estimates of the many-body ground-state wavefunction $|\Phi_0\rangle$ and its associated properties. Note that QMC also has complementary orbital space-based approaches that work within second quantization, as described in section 13.

VMC and DMC are the most commonly applied real-space QMC methods. Within VMC, the simplest scheme, Monte Carlo sampling is used to obtain estimates of the energy of a trial wavefunction

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:vmc} E_{\rm VMC}=\int \Psi^* \hat{H} \Psi . \nonumber \end{align} \tag{ 1 }$$

In DMC, the ground-state wavefunction is obtained by projection of the imaginary time Schrödinger equation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{imaginary_time_schrodinger} -\frac{\partial \Psi}{\partial \beta} = \hat{H}\Psi \nonumber \end{align} \tag{ 2 }$$

to long time, where $\beta = {\rm i}t$ and has units of imaginary time. (Hartree units are used here and throughout, except where noted.) Crucial to both methods is an accurate trial or guiding wavefunction. Clearly, in VMC the trial wavefunction completely determines the accuracy and statistical efficiency of the result. In DMC it is the nodal surface of the trial wavefunction that determines the accuracy, while the overall trial wavefunction determines the statistical efficiency. The full range of supported trial wavefunctions is described in section 6.

The real-space methods use a Hamiltonian within the Born–Oppenheimer approximation (BOA):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{vmc:electron_ham} \hat{H} = -\frac{1}{2} \sum_i \nabla^2_i + \frac{1}{2} \sum_{i\neq j} \frac{1}{|{{\bf r}_i} - {{\bf r}_j}|} + \sum_{i,J} v^{eJ}({{\bf r}_i},{\bf R}_J) \nonumber \end{align} \tag{ 3 }$$

where the lower case indices and positions ${{\bf r}}_i$ refer to the electrons, and the upper case indices and positions ${\bf R}_I$ refer to the ions. In order, the terms in (3) correspond to the kinetic energy of the electrons, the potential energy of the electrons, and the potential energy due to interactions between electrons and ions. The energy contribution due to the Coulomb interactions of the atoms is constant within the BOA, and is computed by an Ewald sum. Further details are given in section 8

For technical reasons, in the following we will work with the 'importance sampled' Schrödinger equation, which can be obtained from the imaginary time Schrödinger equation by rewriting it in terms of $f({\bf r}, \beta) = \Psi_T({\bf r})\Psi({\bf r}, \beta)$.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{pmc:impsamp_schrodinger} \frac{\partial f}{\partial \tau} = \hat{L} f({\bf r},\beta) \nonumber \end{align} \tag{ 4 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle = \left[ \lambda_e \nabla \cdot (\nabla - {\bf F}({\bf r})) - (E_L({\bf r})-E_T) \right]\,f({\bf r},\beta) \nonumber \end{align} \tag{ 5 }$$

where $\Psi_T({\bf r})$ is the trial wavefunction, ${\bf F}({\bf r})=2\nabla \log \Psi_T({\bf r})$ is the 'wavefunction force', and $E_T({\bf r})=\Psi_T({\bf r}){}^{-1}\hat{H}\Psi_T({\bf r})$ is the 'local energy'. $\hat{L}$ is the 'importance-sampled Hamiltonian' operator. To help future discussion, we split $\hat{L}$ into a 'drift/diffusion' operator $\hat{K}$ and a 'branching' operator $\hat{E}$, given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{K} = \frac{1}{2}\nabla \cdot \left(\nabla - {\bf F}({\bf r})\right) \nonumber \end{align} \tag{ 6 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{E} = - \left(E_L({\bf r}) - E_T\right). \nonumber \end{align} \tag{ 7 }$$

One advantage of working with the physical Hamiltonian as opposed to an auxiliary problem (as in Kohn–Sham DFT), is that the variational theorem of quantum mechanics holds, which states that for a trial wavefunction $\Psi_T$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{vmc:var_energy} E_{\Psi_T} = \frac{\langle \Psi_T | \hat{H} | \Psi_T \rangle }{\langle \Psi_T | \Psi_T \rangle} = \frac{\langle \Psi_T | E_L(\hat{{\bf r}}) | \Psi_T \rangle }{\langle \Psi_T | \Psi_T \rangle} \geqslant E_0. \nonumber \end{align} \tag{ 8 }$$

The strict equality holds if $\Psi_T$ is the ground-state of $\hat{H}$. A corollary is that the variance of the trial wavefunction also obeys the following:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{vmc:var_variance} \sigma^2_{\Psi_T} = \frac{\langle \Psi_T | \hat{H}^2 - (E_{\Psi_T})^2 | \Psi_T \rangle}{\langle \Psi_T | \Psi_T \rangle} \geqslant 0. \nonumber \end{align} \tag{ 9 }$$

The variational principle is significant, since it gives us a well-defined metric for which wavefunctions are better or worse approximations to the ground state. This can be turned into an actual algorithm by parameterizing families of wavefunction ansatz with parameters ${\bf c}$. One can then minimize equation (8) with respect to ${\bf c}$ to obtain a best estimate for the state.

The oldest approach for dealing with the Schrödinger equation for realistic systems involves writing an approximation for the ground-state wavefunction and evaluating expectation values. There are two ingredients in this procedure: evaluating equation (8) for some set of variational parameters ${\bf c}$, and then minimizing. The optimization procedure is covered in detail in section 10.

The energy expectation value in equation (8), as well as all physical expectation values, are integrals of a form that are amenable to Metropolis Monte Carlo sampling. Thus, we can evaluate equation (8) (for example) by the following:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{\Psi_T} = \frac{\int {\rm d}{\bf r} |\Psi_T({\bf r})|^2 E_L({\bf r})}{\int {\rm d}{\bf r} |\Psi_T({\bf r})|^2} = \frac{1}{N_s} \sum_{i} E_L({\bf r}(t_i)) + \xi. \nonumber \end{align} \tag{ 10 }$$

N_s is the number of samples and ξ is a Gaussian-distributed statistical error whose variance scales like $1/\sqrt{N_s}$. We write the sample configurations as ${\bf r}(t_i)$ to emphasize that metropolis Monte Carlo generates samples sequentially via a random walk along a Markov chain. To parallelize the algorithm, multiple independent Markov chains or 'walkers' are used.

5.2.1. Trial moves.

QMCPACK supports VMC trial moves with and without drift. This means that the move ${\bf r}\to{\bf r}'$ is drawn from the transition probability distribution given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{trans_prob_alle} T({\bf r}\to{\bf r}',\tau) = \frac{1}{(4\pi\lambda\tau)^{3N/2}}\exp \left(-\frac{\left({\bf r'} -{\bf r} - 2\lambda\tau {\bf F}({\bf r})\right)^2}{4\lambda \tau}\right). \nonumber \end{align} \tag{ 11 }$$

For drift based moves, ${\bf F}({\bf r})$ is taken to be the same wavefunction force as appears in equation (4). For moves without drift, ${\bf F}({\bf r})=0$. In the absence of pathologies in the trial wavefunction, the use of the drift term is almost always more efficient.

In addition to drift or no-drift based moves, the code supports particle-by-particle or all-electron moves. All-electron moves are conceptually the simplest. One proposes the move ${\bf r}\to{\bf r}'={\bf r}+{\bf \Delta}$ by drawing the 3N_e dimensional vector ${\bf \Delta}$ from the distribution in equation (11). This move is then accepted or rejected with probability:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{accept_reject} A({\bf r}\to{\bf r}') = \min \left(1.0, \frac{|\Psi_T({\bf r}')|^2}{|\Psi_T({\bf r})|^2}\frac{T({\bf r}'\to{\bf r})}{T({\bf r}\to{\bf r}')} \right). \nonumber \end{align} \tag{ 12 }$$

In contrast, particle-by-particle moves work by iterating sequentially over all electrons. Considering an electron i at position ${\bf r}_i$. A particle-by-particle move is executed by first drawing a new position for electron i from the following probability distribution.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & T({\bf r}_0,\ldots,{\bf r}_i \to {\bf r}'_i,\ldots,{\bf r}_{N_e})=\frac{1}{(4\pi\lambda\tau)^{3/2}}\nonumber \\ &\times \exp \left(-\frac{\left({\bf r'}_i -{\bf r}_i - 2\lambda\tau {\bf F}_i({\bf r})\right)^2}{4\lambda \tau}\right). \nonumber \end{align} \tag{ 13 }$$

Then, the move is accepted or rejected using a similar acceptance probability as in equation (12). Particle-by-particle moves are typically favored over all-electron moves, due to their higher statistical and numerical efficiencies in practice. However, all-electron moves may be competitive for small systems or for sophisticated trial wavefunctions where single particle moves can not be cheaply evaluated numerically.

One can substantially improve upon the accuracy of VMC by using projector Monte Carlo methods such as DMC. The 'projector' is the formal solution of the imaginary time Schrödinger equation $\newcommand{\e}{{\rm e}} \hat{G}(\beta) = \exp(-\beta\hat{H})$, and has the very desirable property that given any trial wavefunction $|\Psi_T\rangle$ which is non-orthogonal to the ground-state wavefunction, one can obtain the ground-state $|\Phi_0\rangle$ by the following:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lim_{\beta \to \infty} \hat{G}(\beta)|\Psi_T\rangle = {\rm e}^{-\beta E_0}|\Phi_0\rangle. \nonumber \end{align} \tag{ 14 }$$

For efficiency reasons, we consider the projector $\tilde{G}({\bf r}, {\bf r}', \beta)$ associated with the importance sampled Schrödinger equation. For realistic systems, it is exceedingly rare to have exact analytic expressions for the projector. However, we can solve for the Green's function $\hat{\tilde{G}}(\tau)$ of equation (4) approximately for short times τ. Solving the drift/diffusion equations and rate equations independently in the short-time limits, one uses the symmetric Trotter formula:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \exp \left(\tau(\hat{A}+\hat{B})\right)=\exp \left(\frac{\tau}{2}\hat{B} \right)\exp \left(\tau\hat{A} \right)\exp \left(\frac{\tau}{2}\hat{B}\right)+O(\tau^2) \nonumber \end{align} \tag{ 15 }$$

to stitch these independent solutions together into an approximate solution for the importance sampled Green's function:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\bra}[1]{\bigl< #1 |} \newcommand{\ket}[1]{| #1 \bigr>} \displaystyle \label{pmc:g_shorttime} \tilde{G}({\bf r},{\bf r}',\tau) = \bra{{\bf r}}\hat{\tilde{G}}(\tau)\ket{{\bf r}'} = G_{DD}({\bf r},{\bf r}',\tau) G_{B}({\bf r},{\bf r}',\tau)+O(\tau^2). \nonumber \end{align} \tag{ 16 }$$

$G_{DD}({\bf r}, {\bf r}', \tau)$ is the Green's function for the drift/diffusion operator $\lambda_e \nabla \cdot (\nabla - {\bf F}({\bf r}))$. Assuming that ${\bf F}({\bf r})$ is slowly varying, its solution is given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{pmc:g_driftdif\!fusion} \tilde{G}_{DD}({\bf r},{\bf r'},\tau) = \frac{1}{(4\pi\lambda\tau)^{3N/2}}\exp \left(- \frac{({\bf r'}-{\bf r} - 2\lambda \tau {\bf F}({\bf r})) ^2}{4\lambda \tau} \right). \nonumber \end{align} \tag{ 17 }$$

The Green's function for the local energy operator is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{pmc:g_growth} \tilde{G}_B({\bf r},{\bf r}',\tau) = P_0 \exp \left(- \frac{1}{2}(E_L({\bf r})+E_L({\bf r}') - 2E_T)\tau\right). \nonumber \end{align} \tag{ 18 }$$

Near the nodes of $\Psi_T({\bf r})$ and near bare ions, singularities render the 'slowly-varying' approximation used in equation (17) invalid. Improved drift-diffusion projectors have been derived which have been shown to reduce the time step error [46]. QMCPACK implements drift rescaling based on proximity to the nodal surface, following the prescription in [46] for single-electron and all-electron moves. Rescaling based on proximity to bare ions is not yet implemented.

Diffusion Monte Carlo works by stochastically simulating the imaginary-time evolution of an initially prepared state $f({\bf r}, 0)=|\Psi_T({\bf r})|^2$. This is done by exploiting the mathematical correspondence between Fokker–Planck equations for the evolution of probability distributions, and Langevin equations describing the stochastic evolution of particle trajectories. To shift to a Langevin picture, we represent an initial state $f({\bf r}, 0)$ by an ensemble of N_w walkers $\lbrace {\bf r}_0(0), {\bf r}_1(0), \ldots, {\bf r}_{N_w-1}(0) \rbrace$ distributed according to $f({\bf r}, 0)$. Assume each walker also has an associated weight $w_i(\beta)$ where w_i(0)  =  1.0. We consider the action of the short-time Green's function $\tilde{G}({\bf r}, {\bf r}', \tau)$ on this distribution.

The action of the drift-diffusion propagator can be simulated with a random drift-diffusion step, given by:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{pmc:driftdif\!fusion_step} {\bf r}(\beta+\tau) = {\bf r}(\beta) + 2\lambda\tau{\bf F}({\bf r}(\beta)) + \sqrt{2\lambda\tau} \boldsymbol{\xi}. \nonumber \end{align} \tag{ 19 }$$

Here, $\boldsymbol{\xi}$ is a 3N_e dimensional Gaussian random vector with unit variance. Once a new position generated, the $G_B({\bf r}, {\bf r}', \tau)$ contribution is dealt with by updating the walker weight $w_i(\beta)$ with the following formula:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{pmc:walker_weight} w(\beta+\tau) = w(\beta)G_B({\bf r}(\beta),{\bf r}(\beta+\tau),\tau). \nonumber \end{align} \tag{ 20 }$$

Expectation values of local observables $\mathcal{A}({\bf r})$ over the distribution $f({\bf r}, \beta)$ are obtained by a weighted average:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{pmc:weighted_average} \langle\mathcal{A}\rangle_f(\beta)=\frac{\sum_{i=0}^{N_w-1} w_i(\beta)\mathcal{A}({\bf r_i(\beta))}}{\sum_{i=0}^{N_w-1} w_i(\beta)}. \nonumber \end{align} \tag{ 21 }$$

Implementing everything discussed up to this point results in 'pure diffusion' Monte Carlo. However, due to the exponential growth/decay of the walker weights with β, the efficiency of this method decays exponentially with the projection time. 'Branching' diffusion Monte Carlo circumvents this problem by implementing equation (20) stochastically through the replication/removal of walkers. For each walker i, M_i copies of the walker are made after the drift/diffusion step according to the formula $M_i = {\rm INT}(w_i(\beta+\tau)+\xi)$, where ξ is a uniform random number between 0 and 1. The weights of these M_i walkers are renormalized and the copies then proceed to the next time step. Notice that M_i  =  0 implies that the walker is killed.

To avoid a walker population explosion or collapse, practical DMC simulations adjust E_T dynamically to keep the population finite and stable. In QMCPACK this is achieved using either a variable number of walkers combined with the above population control, or via a fixed walker count scheme ('stochastic reconfiguration' [45]). Both schemes can potentially introduce a 'population control bias' that must be checked and controlled for, particularly for small populations. To minimize or check for bias, the population should be as large as possible, should be allowed to fluctuate significantly, and the simulation run for a long time.

To correctly simulate Fermionic systems and avoid a collapse of the propagated wavefunction to a Bosonic solution, the 'fixed-node approximation' is implemented. This constrains the projected solution to the nodal surface of the trial wavefunction, thereby preserving Fermionicity. Proposed moves that result in a nodal crossing are detected through the change in sign of the wavefunction and are rejected. This is usually the most significant approximation made, and requires accurate trial wavefunctions.

Finally, we note that several important modifications for production calculations: QMCPACK implements the 'small time-step error' algorithm due to Umrigar et al [46], where the drift term is modified near wavefunction nodes and effective time step introduced to improve the time step convergence of the algorithm. The recently proposed size-consistent variation [7, 11] is also implemented. This is particularly effective for computing energy differences between very different sized-systems, such as absorption energies or the formation energies of molecular crystals [11].

Reptation Monte Carlo is constructed by exploiting the Feynman path-integral formulation of Schrödinger's equation [47]. Its primary advantages over DMC are its ability to estimate observables over pure distributions in polynomial time, and lack of population bias and control issues. Consider the ground-state 'partition function':

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{Z}(\beta) = \langle \Psi_T | {\rm e}^{-\beta\hat{H}} | \Psi_T \rangle. \nonumber \end{align} \tag{ 22 }$$

First, we split ${\rm e}^{-\beta\hat{H}}$ into n segments each spanning an imaginary time $\tau=\beta/n$. After inserting n  +  1 position space resolutions of the identity and rewriting the resulting expression in terms of the importance sampled projector, we find that $\mathcal{Z}(\beta)$ can be written as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{Z}(\beta) = \int {\rm d}{\bf r}(t_0)\ldots {\rm d}{\bf r}(t_n) \mathcal{P}[X] {\rm e}^{-\sum_{i=0}^{n-1} L_{\rm DMC}({\bf r}(t_i),{\bf r}(t_{i+1})) } \nonumber \end{align} \tag{ 23 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{P}[X] = |\Psi_T({\bf r}(t_0))|^2 \tilde{G}_{DD}({\bf r}(t_0),{\bf r}(t_1)) \ldots \tilde{G}_{DD}({\bf r}(t_{n-1}),{\bf r}(t_n)) \nonumber \end{align} \tag{ 24 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle L_{\rm DMC}({\bf r},{\bf r}') = \frac{\tau}{2} \left(E_L({\bf r}) + E_L({\bf r}') -2E_T\right). \nonumber \end{align} \tag{ 25 }$$

X is shorthand for a 'path' $X=({\bf r}(t_0), \ldots, {\bf r}(t_n))$. The reader will recognize $\mathcal{Z}(\beta)$ as a path integral. L_DMC is called the 'link action', which when summed over t, gives the action for the path X. $\mathcal{P}[X]$ is the probability of a walker which is initially distributed according to $|\Psi_T({\bf r})|^2$ executing the directed random walk from ${\bf r}(t_0)$ to ${\bf r}(t_n)$ along the path X.

All reptation moves in QMCPACK use the 'bounce algorithm' [48]. For proposed moves, the improved propagators described in the DMC section are directly used in reptation. In addition to nodal drift-rescaling, we incorporate the DMC effective time step and energy filtering methods directly into the link-action, which helps to significantly reduce ergodicity problems associated with reptiles getting stuck in low energy regions of configuration space. In the event that reptiles still get stuck, the age of all reptile beads is accumulated. If a bead exceeds some specified age, the entire reptile is forced to propagate for n steps without rejection and then re-equilibrate.

Since $\mathcal{P}[X]$ cancels out of the reptile accept/reject step, any all-electron move which is a valid VMC configuration is supported in RMC. In addition to traditional all-electron moves, QMCPACK also supports reptile proposals which are built from a sequence of N_e particle-by-particle moves. Reptile moves proposed with these 'particle-by-particle' moves exhibit higher acceptance ratios than the traditional all-electron moves, and are thus favored if memory is available.

Propagation in RMC is supported for all-electron, local, and semi-local pseudopotential Hamiltonians. The fixed-node constraint is enforced by rejecting proposed node crossings and immediately bouncing.

Since the random walk of each reptile is totally independent of other reptiles, RMC is straightforwardly parallelized. In addition to generic MPI parallelization, QMCPACK 's RMC driver is able to place one reptile per OpenMP thread on shared memory systems.

For the vast majority of molecular and solid-state studies, the trial wavefunction is written as the product of an antisymmetric function and a symmetric Jastrow function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:msdj} \Psi_T=\sum_{i=1}^{M}c_{i}D^\uparrow_iD^\downarrow_i e^J, \nonumber \end{align} \tag{ 26 }$$

where the N electron trial wavefunction $\Psi_T$ is expanded in a weighted sum of products of up and down spin determinants, D. These are in turn multiplied by a real-space Jastrow factor, J. When exponentiated, this factor is nodeless and the nodes of the trial wavefunction are therefore purely determined by the determinantal parts. A single product of up-spin and down-spin determinants would correspond to a mean-field or Hartree–Fock starting point. Larger determinantal sums can be obtained, e.g. from multi-configuration self-consistent field quantum chemical calculations, CIPSI section 13.3, or be constructed based on physical or chemical reasoning. Excited states may be constructed by manipulating the occupancy of the Slater determinants in the input, e.g. to create an exciton. Wavefunctions with greater or fewer electrons than the neutral ground state may be similarly prepared to compute electronic affinities or ionization potentials.

Due to the potentially large computational cost in evaluating the trial wavefunction, QMCPACK uses previously computed data and optimized methods to avoid full recomputation wherever effective and practical. For single electron moves, QMCPACK uses the Sherman–Morrison algorithm, as described in [49]. For large calculations with thousands of electrons, the delayed update scheme of [44] is currently being implemented. For calculations with multiple determinants, QMCPACK implements the 'table method' of Clark et al [50]. This exploits the relationship between largely similar determinants to cheaply compute the determinant values while only requiring the full N² memory cost of a single determinant. This enables, e.g. molecular calculations that approach or even reach 'chemical accuracy' to be performed [51]. To date, calculations with up to O(10⁶) determinants have been performed (see section 13.3), with larger calculations clearly possible [52].

The single-particle orbitals in the Slater determinants are generally determined by another electronic structure code and imported into QMCPACK for calculations. QMCPACK has an easily extensible mechanism for adding new ways of representing single-particle orbitals. This can be particularly useful when addressing model systems or performing specialized tests. For example, QMCPACK supports specialty homogeneous electron gas and plane-wave based wavefunctions, and for work on spherical quantum dots, radial numerical functions [53]. However, by far the most common sources of orbitals are plane-wave based from Quantum Espresso (QE) [54, 55], and Gaussian based from the GAMESS code [56]. Converters from these codes are provided and can straightforwardly be extended to other methodologically related codes.

6.2.1. B-spline basis sets.

6.2.2. Gaussian basis sets.

6.2.3. Specialized basis sets.

Improvement of the nodal surface can be achieved through backflow wavefunctions, complementing the multideterminant route. The formal justification for backflow wavefunctions rests on the homogeneous electron gas and Fermi liquid theory [61]. Backflow appears promising for bulk applications [62], and has also been shown to aid in capturing dynamical correlations in molecular systems when used in conjunction with multideterminant wavefunctions [63].

Backflow wavefunctions are constructed from determinantal wavefunctions as follows. Instead of evaluating the Slater matrix $M_{ij}=\phi_j ({\bf r}_i)$ at the bare electron coordinates ${\bf r}_i$, we evaluate it at new quasiparticle coordinates $\tilde{M}_{ij}=\phi_j ({\bf q}_i)$. The 'backflow transformation' from ${\bf r}_i\to {\bf q}_i$ is defined as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf q}_{i_\alpha}={\bf r}_{i_\alpha}+\sum_{\alpha \leqslant \beta} \sum_{i_\alpha \neq j_\beta} \eta^{\alpha\beta}(|{\bf r}_{i_\alpha}-{\bf r}_{j_\beta}|)({\bf r}_{i_\alpha}-{\bf r}_{j_\beta}). \nonumber \end{align} \tag{ 27 }$$

In QMCPACK, the $\newcommand{\e}{{\rm e}} \eta^{\alpha\beta}(r)$ are short-ranged, spherically symmetric functions represented by fully optimizable B-splines. QMCPACK allows for separate optimization of same-spin, opposite-spin, and electron-ion terms. Currently, backflow is fully supported only with single determinant wavefunctions, but it can be used in both bulk and molecular systems.

The long-ranged Coulombic interactions of the electrons and ions must be handled with care in order to ensure that the potential energy does not diverge when using periodic boundary conditions. In QMCPACK, the interparticle interactions are computed using an optimized implementation [76] of the well-known strategy of decomposing the interactions into short and long ranged components, and performing sums over the former and latter in real and reciprocal space, respectively [75].

Bloch's theorem demonstrates how a finite wavefunction can be used to simulate an infinite lattice within periodic boundary conditions by incorporating the following symmetry:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Psi({\bf r}_{1} + {\bf L_{m}}, {\bf r}_{2}, \dots, {\bf r}_{N}) = {\rm e}^{{\rm i}{\bf K}\cdot{\bf L}_{m}} \Psi({\bf r}_{1}, {\bf r}_{2}, \dots, {\bf r}_{N}) \nonumber \end{align} \tag{ 43 }$$

where ${\bf K}$ is a vector in reciprocal space, ${\bf L}_{m}$ is a lattice vector of the supercell, and $\newcommand{\e}{{\rm e}} \Theta_{m} \equiv {\bf K} \cdot {\bf L}_{m}$, is the 'twist angle' [2]. For pure periodic boundary conditions (in which $\Theta=0$), systems converge slowly to their thermodynamic limit due to shell effects and quantization of momentum [79]. Therefore, to improve convergence speed and accuracy, one should average over many simulations done with different twist angles, a scheme called 'twist-averaged boundary conditions'. In QMCPACK, the averaging is done in post-processing, using e.g. qmca section 14.1 and/or Nexus section 15.

QMCPACK optimizes trial functions using an implementation of the linear method (LM) [80] that includes modifications to improve stability in the face of variables of greatly differing stiffnesses, facilitate the optimization of excited states, and reducing the memory footprint when optimizing large numbers of variational parameters. The LM is sometime referred to as 'energy minimization', although the approach is more general. The LM gets its name from the way that it employs a linear expansion of the wavefunction,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:lm_wfn} |\Phi\rangle = \sum_{y=0}^{N_{\rm v}} c_y |\Psi^y\rangle, \nonumber \end{align} \tag{ 44 }$$

where $|\Psi^y\rangle$ for $y\in\{1, 2, 3, ...\}$ is the derivative of the trial function $|\Psi\rangle$ with respect to its yth variational parameter and $\newcommand{\e}{{\rm e}} |\Psi^0\rangle\equiv |\Psi\rangle$, within an expanded energy expression,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:lm_expanded_energy} E(\vec{c}) =\frac{\langle\Phi|H|\Phi\rangle}{\langle\Phi|\Phi\rangle}. \nonumber \end{align} \tag{ 45 }$$

Using this linear approximation to how the energy changes with the variational parameters, minimizing E with respect to $\vec{c}$ can be achieved by solving the generalized eigenvalue problem

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:lm_gev} \sum_{y=0}^{N_{\rm v}} \frac{\langle\Psi^x|H|\Psi^y\rangle}{\langle\Psi|\Psi\rangle} c_y = E \sum_{y=0}^{N_{\rm v}} \frac{\langle\Psi^x|\Psi^y\rangle}{\langle\Psi|\Psi\rangle} c_y \nonumber \end{align} \tag{ 46 }$$

or, written in matrix-vector notation,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:lm_gev_mv} {\bf H} \hspace{1mm} \vec{c} = E \hspace{1mm} {\bf S} \hspace{1mm} \vec{c}, \nonumber \end{align} \tag{ 47 }$$

the matrix elements for which are evaluated by Monte Carlo integration [81, 82] in direct analogy to how VMC evaluates the energy. If one assumes the improved trial function $|\Phi\rangle$ is similar to the previous trial function $|\Psi\rangle$, which implies that the ratio $c_x/c_0$ is small for all $x\in\{1, 2, 3, ...\}$, then a reasonable approximation to $|\Phi\rangle$ can be had by replacing $\mu_x\rightarrow\mu_x+c_x/c_0$ for each variational parameter $\mu_x$ in $|\Psi\rangle$. As for other optimization methods that compute an update based on some local approximation to the target function, such as Newton–Raphson, this process is then repeated until further updates no longer lower the energy.

In practice, it is important to implement an analogue to the trust radius schemes common to Newton–Raphson in order to ensure that the solution of equation (46) does not correspond to an unreasonably long step in variable space, or, put another way, to ensure that the ratio $c_x/c_0$ is not too large. The LM optimizer in QMCPACK supports two mechanisms for preventing too-large updates: a diagonal shift α as employed in the original algorithm [80] as well as an overlap-based shift β that becomes important when parameters of greatly different stiffnesses are present. Using these shifts, the Hamiltonian matrix is modified to become

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:lm_shifted_ham} {\bf H} \rightarrow {\bf H} + \alpha {\bf A} + \beta {\bf B}, \nonumber \end{align} \tag{ 48 }$$

where ${\bf A}$ and ${\bf B}$ provide stabilization via the original and overlap shifts, respectively. As in the original method, QMCPACK uses $A_{xy}=\delta_{xy}(1-\delta_{x0})$ and the adjustable shift strength α to effectively raise the energy along each direction of change while leaving the current wavefunction $|\Psi^0\rangle=|\Psi\rangle$ unaffected.

While the original shifting scheme has been effective in many cases, it can struggle if two different variational parameters produce wavefunction derivatives of vastly different sizes. For example, imagine a two-variable wavefunction whose overlap matrix evaluates to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:lm_bad_overlap} {\bf S} = \left[ \begin{array}{@{}c c c@{}} 1 \hspace{1.5mm} & 0 & 0 \nonumber \\ 0 & 1 & 0 \nonumber \\ 0 & 0 & 10^6 \end{array} \right]. \nonumber \end{align} \tag{ 49 }$$

Performing the usual $\vec{\nu}=({\bf S}^{\frac{1}{2}})\vec{c}$ transformation to produce a standard eigenvalue problem (with β set to zero for now) gives us

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:lm_innefective_shift} {\bf S}^{-\frac{1}{2}} {\bf H} {\bf S}^{-\frac{1}{2}} \vec{\nu} + \left[ \begin{array}{@{}c c c@{}} 1 & 0 & 0 \nonumber \\ 0 & \alpha & 0 \nonumber \\ 0 & 0 & \frac{\alpha}{10^6} \nonumber \\ \end{array} \right] \vec{\nu} = E \hspace{0.3mm}\vec{\nu}. \nonumber \end{align} \tag{ 50 }$$

We see that, if we were to make α large enough to significantly penalize the second variable direction, the first direction would be penalized so much that it would essentially become a fixed parameter.

The purpose of the overlap shift is to resolve this issue by adding an energy penalty based on the norm of the part of $\vec{c}$ corresponding to directions orthogonal to the current wavefunction $|\Psi\rangle$, which would correctly penalize steps along directions of large derivative norms more than those along directions of small derivative norms. This goal is accomplished by the definitions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:lm_overlap_shift_matrices} Q_{xy} = \delta_{xy} - \delta_{x0}(1-\delta_{y0}) S_{0y} \nonumber \end{align} \tag{ 51 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_{xy} = (1 - \delta_{x0}\delta_{y0}) \left[{\bf Q}^+{\bf S} {\bf Q}\right]_{xy} \nonumber \end{align} \tag{ 52 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf B} = \left({\bf Q}^+\right)^{-1} {\bf T} {\bf Q}^{-1} \nonumber \end{align} \tag{ 53 }$$

in which ${\bf Q}$ transforms into a basis in which all update directions are orthogonal to the current wavefunction $|\Psi\rangle$ (this transformation is equivalent to that of equation (24) of [81]). ${\bf T}$ is the overlap matrix in this basis with its first element zeroed out so that the current wavefunction is not penalized. Finally, the inverses of ${\bf Q}$ and its adjoint transform us back to the basis of the original generalized eigenvalue problem so that the effect of the overlap shift cay be written in the form of equation (48). Note that, in practice, it is not necessary to construct ${\bf B}$ explicitly, as QMCPACK solves the generalized eigenvalue equation by iterative Krylov subspace expansion, during which the Krylov basis (whose first vector is always $|\Psi\rangle$) is kept orthonormal by the Gram–Schmidt procedure. In this Krylov basis, applying the overlap shift involves merely adding β to the diagonal of the subspace Hamiltonian matrix (except, of course, to the first element corresponding to $|\Psi\rangle$). This Krylov approach also has the benefit of ensuring that the overall update is orthogonal to the current wavefunction, which is related to norm-conservation and was found to be desirable by the LM's original developers [80–82].

Although like most trust-radius schemes the optimal choices for α and β are somewhat heuristic, QMCPACK automatically adjusts them after each iteration of the LM by solving for the updates generated by three different sets of shifts and retaining the shift that gave the best update, as determined by a correlated-sampling comparison of their energies on a fresh sample. For maximum efficiency in regimes where optimization is not difficult but sampling is expensive, QMCPACK retains the ability to run in a single-shift, no-second-sample mode. When running instead in multi-shift mode, we have observed that successful optimizations often result with the simple initial choice of $\alpha=\beta=1$. In principle, however, one might expect $\alpha<\beta$ to be more effective, because when the β shift is filling the role of limiting the update size, α is only needed to penalize (hopefully rare) linear dependencies between update directions that β, being overlap-based, cannot address.

QMCPACK's current LM optimization engine supports both standard energy minimization and the minimization of a recently introduced [83] excited state target function, $\langle\Psi|(\omega-\hat{H})|\Psi\rangle / \langle\Psi|(\omega-\hat{H}){}^2|\Psi\rangle$, whose global minimum is the exact energy eigenstate immediately above the targeted energy ω. Although this technology is a very recent development and will doubtless evolve in time as the science behind excited state targeting matures, we felt it important to make an early version of it available to the community. Optimization proceeds in much the same way as for a ground state, with the user specifying ω and the stabilization shifts α and β and the LM repeatedly solving generalized eigenvalue equations analogous to equation (46) to generate wavefunction updates. Additional methods for automatically selecting and updating ω have been developed [84]. For details into this targeting function and how it is optimized, we refer the reader to the original publication [83].

One important limitation of the LM comes when the number of variational parameters rises to 10 000 or more, at which point the contributions to ${\bf H}$ and ${\bf S}$ made by each Markov chain become cumbersome to store in memory, especially when running one Markov chain per core on a large parallel system in which per-core memory is limited. QMCPACK currently addresses this memory bottleneck using the blocked LM [35], a recent algorithm that separates the variable space into blocks, estimates the most important variable-change directions within each block, and then uses these directions to construct a reduced and vastly more memory efficient LM eigenvalue problem to generate an update direction in the overall variable space. Like excited state targeting, this is a new feature that can be expected to evolve in time, and has been made openly available to the community in the spirit of rapid dissemination. As of this writing, it has not been widely tested outside of the work in its original publication [35], but in time we expect to have a clearer picture of its capabilities.

QMCPACK also supports optimizing variational parameters based on not only the total energy but also variance. In certain situations, the best target object may not be the energy only but a cost function mixing both energy and variance which reduces to zero when the wavefunction is exact. The cost function can be any linear combination of energy and variance. QMCPACK picks the optimal parameter set corresponding to the minimal value of a quartic function fitting the cost function evaluated on seven shifts by correlated-sampling.

The particle number density operator is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{n}_r=\sum_{i} \delta ({\bf r}-{\bf r}_i). \nonumber \end{align} \tag{ 54 }$$

This estimator accumulates the number density on a uniform histogram grid over the simulation cell. The value obtained for a grid cell c with volume $\Omega_c$ is then the average number of particles in that cell

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle n_c=\int{} {\rm d}{\bf R}|\Psi|^2\int{_{\Omega_c}}{\rm d}{\bf r}\sum_i\delta({\bf r}-{\bf r}_i). \nonumber \end{align} \tag{ 55 }$$

When using periodic boundary conditions, the density will be collected for the cell (or supercell) defined by the simulation. When using non-periodic boundary conditions, a cell has to be defined in order to set a grid. It is then recommended to center the system (molecule) in the middle of the defined cell. The collected data is stored in HDF5 format in a .stat.h5 file. Using Nexus section 15, one can use the qdens tool to extract the data in a *.xsf format readable with visualization tools such as XCrysDens [88, 89] or VESTA [90]. Examples of density plots are shown in figure 4. Similar to the density, the spin-density estimator can be also collected for each independent spin, as shown and analyzed in [91] for magnetic states in ${\rm Ti_4O_7}$.

Zoom In Zoom Out Reset image size

The functional form of the species-resolved radial pair correlation function operator is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g_{ss'}(r) = \frac{V}{4\pi r^2N_sN_{s'}}\sum_{i_s=1}^{N_s}\sum_{j_{s'}=1}^{N_{s'}}\delta(r-|{\bf r}_{i_s}-{\bf r}_{j_{s'}}|). \nonumber \end{align} \tag{ 56 }$$

Here N_s is the number of particles of species s and $V$ is the supercell volume. If $s=s'$, then the sum is restricted so that $i_s\ne j_s$.

An estimate of $g_{ss'}(r)$ is obtained as a radial histogram with a set of N_b uniform bins of width $\delta r$. This can be expressed analytically as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{g}_{ss'}(r) = \frac{V}{4\pi r^2N_sN_{s'}}\sum_{i=1}^{N_s}\sum_{j=1}^{N_{s'}}\frac{1}{\delta r}\int_{r-\delta r/2}^{r+\delta r/2}{\rm d}r'\delta(r'-|{\bf r}_{si}-{\bf r}_{s'j}|), \nonumber \end{align} \tag{ 57 }$$

where the radial position r is restricted to reside at the bin centers $\delta r/2, 3 \delta r/2, \ldots$.

Let $\rho^e_{{\bf k}}=\sum_j {\rm e}^{{\rm i} {\bf k}\cdot{\bf r}_j^e}$ be the Fourier space electron density, with ${\bf r}^e_j$ being the coordinate of the jth electron. ${\bf k}$ is a wavevector commensurate with the simulation cell. The static electron structure factor $S({\bf k})$ can be measured at all commensurate ${\bf k}$ such that $|{\bf k}| \leqslant ({\rm LR\_DIM\_CUTOFF}) r_c$. N^e is the number of electrons, LR_DIM_CUTOFF is the optimized breakup parameter, and r_c is the Wigner–Seitz radius. It is defined as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S({\bf k}) = \frac{1}{N^e}\langle \rho^e_{-{\bf k}} \rho^e_{{\bf k}} \rangle. \nonumber \end{align} \tag{ 58 }$$

An energy density operator, $\hat{\mathcal{E}}_r$, satisfies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int {\rm d}r \hat{\mathcal{E}}_r = \hat{H}, \nonumber \end{align} \tag{ 59 }$$

where the integral is over all space and $\hat{H}$ is the Hamiltonian. In QMCPACK, the energy density is split into kinetic and potential components

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{\mathcal{E}}_r = \hat{\mathcal{T}}_r + \hat{\mathcal{V}}_r \nonumber \end{align} \tag{ 60 }$$

with each component given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{\mathcal{T}}_r &= \frac{1}{2}\sum_i\delta({\bf r}-{\bf r}_i)\hat{p}_i^2 \nonumber \\ \hat{\mathcal{V}}_r &= \sum_{i<j}\frac{\delta({\bf r}-{\bf r}_i)+\delta({\bf r}-{\bf r}_j)}{2}\hat{v}^{ee}({\bf r}_i,{\bf r}_j)\nonumber \\ &\quad + \sum_{i\ell}\frac{\delta({\bf r}-{\bf r}_i)+\delta({\bf r}-\tilde{{\bf r}}_\ell)}{2}\hat{v}^{eI}({\bf r}_i,\tilde{{\bf r}}_\ell) \nonumber \\ &\quad + \sum_{\ell< m}\frac{\delta({\bf r}-\tilde{{\bf r}}_\ell)+\delta({\bf r}-\tilde{{\bf r}}_m)}{2}\hat{v}^{II}(\tilde{{\bf r}}_\ell,\tilde{{\bf r}}_m).\nonumber \end{align} \tag{ 61 }$$

Here ${\bf r}_i$ and $\newcommand{\e}{{\rm e}} \tilde{{\bf r}}_\ell$ represent electron and ion positions, respectively, $\hat{p}_i$ is a single electron momentum operator, and $\hat{v}^{ee}({\bf r}_i, {\bf r}_j)$, $\newcommand{\e}{{\rm e}} \hat{v}^{eI}({\bf r}_i, \tilde{{\bf r}}_\ell)$, $\newcommand{\e}{{\rm e}} \hat{v}^{II}(\tilde{{\bf r}}_\ell, \tilde{{\bf r}}_m)$ are the electron–electron, electron–ion, and ion–ion pair potential operators (including non-local pseudopotentials, if present). This form of the energy density is size-consistent, i.e. the partially integrated energy density operators of well separated atoms gives the isolated Hamiltonians of the respective atoms. For periodic systems with twist averaged boundary conditions, the energy density is formally correct only for either a set of supercell k-points that correspond to real-valued wavefunctions, or a k-point set that has inversion symmetry around a k-point having a real-valued wavefunction. For more information about the energy density, see [94].

The energy density can be accumulated on piecewise uniform three dimensional grids in generalized Cartesian, cylindrical, or spherical coordinates. The energy density integrated within Voronoi volumes centered on ion positions is also available. The total particle number density is also accumulated on the same grids by the energy density estimator for convenience so that related quantities, such as the regional energy per particle, can be computed easily.

The N-body density matrix in DMC is $\newcommand{\ket}[1]{| #1 \bigr>} \newcommand{\bra}[1]{\bigl< #1 |} \newcommand{\operator}[3]{\ket{#1} #2 \bra{#3}} \hat{\rho}_N=\operator{\Psi_{T}}{}{\Psi_{FN}}$ (for VMC, substitute $\Psi_T$ for $\Psi_{FN}$). The one body reduced density matrix (1RDM) is obtained by tracing out all particle coordinates but one:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\operator}[3]{\ket{#1} #2 \bra{#3}} \newcommand{\bra}[1]{\bigl< #1 |} \newcommand{\ket}[1]{| #1 \bigr>} \displaystyle \hat{n}_1 = \sum_nTr_{R_n}\operator{\Psi_{T}}{}{\Psi_{\rm FN}}. \nonumber \end{align} \tag{ 62 }$$

In the formula above, the sum is over all electron indices and $\newcommand{\ket}[1]{| #1 \bigr>} \newcommand{\bra}[1]{\bigl< #1 |} \newcommand{\expval}[3]{\bra{#1}#2\ket{#3}} \newcommand{\e}{{\rm e}} Tr_{{\bf R}_n}(*)\equiv \int {\rm d}{\bf R}_n\expval{{\bf R}_n}{*}{{\bf R}_n}$ with ${\bf R}_n=[{\bf r}_1, ..., {\bf r}_{n-1}, {\bf r}_{n+1}, ..., {\bf r}_N]$. When the sum is restricted over spin up or down electrons, one obtains a density matrix for each spin species. The 1RDM computed by QMCPACK is partitioned in this way.

In real space, the matrix elements of the 1RDM are

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\expval}[3]{\bra{#1}#2\ket{#3}} \newcommand{\bra}[1]{\bigl< #1 |} \newcommand{\ket}[1]{| #1 \bigr>} \displaystyle n_1({\bf r},{\bf r}') = \expval{{\bf r}}{\hat{n}_1}{{\bf r}'} = \sum_n\int {\rm d}{\bf R}_n \Psi_T({\bf r},{\bf R}_n)\Psi_{FN}^*({\bf r}',{\bf R}_n). \nonumber \end{align} \tag{ 63 }$$

A more efficient and compact representation of the 1RDM is obtained by expanding in single particle orbitals, e.g. from a Hartree–Fock or DFT calculation, $\{\phi_i\}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\expval}[3]{\bra{#1}#2\ket{#3}} \newcommand{\bra}[1]{\bigl< #1 |} \newcommand{\ket}[1]{| #1 \bigr>} \displaystyle \label{eq:dm1b_direct} n_1(i,j) &= \expval{\phi_i}{\hat{n}_1}{\phi_j} \nonumber \\ &= \int {\rm d}{\bf R} \Psi_{FN}^*({\bf R})\Psi_{T}({\bf R}) \sum_n\int {\rm d}{\bf r}'_n \frac{\Psi_T({\bf r}_n',{\bf R}_n)}{\Psi_T({\bf r}_n,{\bf R}_n)}\phi_i({\bf r}_n')^* \phi_j({\bf r}_n). \nonumber \end{align} \tag{ 64 }$$

The integration over ${\bf r}'$ in equation (64) is inefficient when one is also interested in obtaining matrices involving energetic quantities, such as the energy density matrix [94] or the related and more well known Generalized Fock matrix. For this reason, we compute: [94]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle n_1(i,j) \approx \int {\rm d}{\bf R} \Psi_{FN}({\bf R})^*\Psi_T({\bf R}) \sum_n \int {\rm d}{\bf r}_n' \frac{\Psi_T({\bf r}_n',{\bf R}_n)^*}{\Psi_T({\bf r}_n,{\bf R}_n)^*}\phi_i({\bf r}_n)^* \phi_j({\bf r}_n'). \nonumber \end{align} \tag{ 65 }$$

For all-electron calculations, naïvely estimating the bare Coulomb Hellman–Feynman force with quantum Monte Carlo suffers from a fatal problem: while the expectation value of this estimator is well defined, the 1/r² divergence causes the variance to be infinite, meaning we can not obtain a meaningful error bar for this quantity. There are several schemes to circumvent this. For all-electron calculations, QMCPACK can currently calculate forces and stress using the Chiesa estimator [95] in both open and periodic boundary conditions. Implementation details and validation of forces in periodic boundary conditions can be found in [96]. In the future, pseudopotential forces will be supported, and methods to reduce the variance of existing estimators are currently being explored.

In addition to real-space QMC methods, QMCPACK also supports orbital-space QMC approaches for the study of atomic, molecular and solid-state systems. AFQMC is implemented internally, while interfaces to selected configuration interaction (SCI) methods have been developed. [97–100] The starting point of orbital-space approaches is the Hamiltonian in second quantization, typically defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{H}=\sum_{i,j} h_{ij} \hat{c}^{\dagger}_i \hat{c}_j + \frac{1}{2} \sum_{i,j,k,l} V_{ijkl} \hat{c}^{\dagger}_i \hat{c}^{\dagger}_j \hat{c}_l \hat{c}_k, \nonumber \end{align} \tag{ 66 }$$

where $\hat{c}^{\dagger}_i (\hat{c}_i)$ are the creation (annihilation) operators for spinors associated with a given single particle basis set, with associated one- and two-electron matrix elements given by h_ij and $V_{ijkl}$. The choice of the single particle basis along with the calculation of the appropriate Hamiltonian matrix elements must be performed by a separate electronic structure package. The Hamiltonian matrix elements are expected in the FCIDUMP format used by codes including Molpro [101], PySCF [102] and VASP [103–105]. The calculations are typically performed on a single particle basis defined by the solution of a Hartree–Fock or DFT calculation. Both finite and periodic calculations are possible.

The fundamental idea behind AFQMC is identical to that of DMC, namely that the propagation of many-body states in imaginary time leads to the lowest eigenstate of the Hamiltonian with non-zero overlap [27]. In contrast with DMC, AFQMC operates in the Hilbert space of non-orthogonal Slater determinants and uses the Hubbard–Stratonovich transformation [106, 107] to express the short-time approximation of the propagator as an integral over propagators that contain only one-body terms. The application of one-body propagators to walkers in the algorithm (represented by Slater determinants) leads to rotations of the corresponding Slater determinants that define a random walk, similar to the random walk in real-space followed by walkers in DMC [108]. QMCPACK implements the constrained-path algorithm of Zhang and Krakauer with the phaseless approximation [27, 109]. Similar to the algorithm of Zhang and Krakauer, we use importance sampling and force-bias to improve the sampling efficiency of the algorithm. For a complete description of the implemented algorithm, see the lecture notes on AFQMC by Zhang in the open-access book [28] and [29].

The AFQMC implementation in QMCPACK attempts to minimize the memory requirements of the calculation, while increasing the performance of the associated computations. This is done by a combination of: (1) distributed sparse representations of large data structures (e.g. two-electron integrals), (2) efficient use of shared-memory on multi-core architectures, (3) combination of efficient BLAS and sparse-BLAS routines for all major computations, and (4) an efficient distributed algorithm for walker propagation. Notice that the code is able to distribute the work associated with the propagation of a walker over many nodes, enabling access to systems with thousands of basis functions with a full ab initio representation. Both single determinant as well as multi-determinant trial wave-functions are implemented. In the case of multi-determinant expansions, both orthogonal as well as non-orthogonal expansions are efficiently implemented. For orthogonal expansions, a fast algorithm based on the Sherman–Morrison–Woodbury formula is implemented which leads to a modest increase in computing time for determinant expansions involving even many thousands of terms.

As discussed previously, a direct path towards improving the accuracy of a QMC calculation is through a better trial wavefunction. One approach is to use selected CI methods such as CIPSI (configuration interaction using a perturbative selection done iteratively), or the recently developed adaptive sampling CI (ASCI) [99] and heat bath CI (HBCI) [100]. The principle behind selected CI methods was first published in 1955 by Nesbet [97]. The first calculations on atoms were performed by Diner, Malrieu and Claverie [110] in 1967. Many advances have since been made with selected CI techniques, and it has been applied widely to atomic, molecular and periodic systems [111–118]. The method is based on an iterative process during which a wavefunction is improved at each step. During each iteration, the current wavefunction is used in conjunction with the Hamiltonian to find important contributions that will be added to the wavefunction in the next iteration. In most selected CI approaches, the importance of a contribution is determined from a many body perturbation theory estimate. A full description of CIPSI, its algorithms, and results on various systems can be found in [98, 119, 120]. A description of new improvements to selected CI techniques that have been demonstrated with ASCI and HBCI can be found in [99, 100]. The CIPSI method [98, 119–121] is implemented in the Quantum Package (QP) code [59] developed by the Caffarel group. QMCPACK does not implement CIPSI, but is able to use output from the QP code via tight integration.

In the following we use the C₂O₂H₃N molecule, figure 5, to illustrate the use of CIPSI to obtain an improved trial wavefunction. The C₂O₂H₃N molecule is part of the cycloreversion of heterocyclic rings database [122], for which the geometry was optimized with DFT using the B3LYP function in a 6–31G basis set. Orbitals are represented within the aug-ccpVTZ basis set. The energetics of this molecule are known to have a strong dependence on the choice of functional in DFT simulations [122]. Diagnostics based on coupled cluster theory (CC) with single, double, and peturbative triple excitations (CCSD(T)) [123] suggest a multireference character [124], a known problem for these techniques [125]. The multireference capability of DMC-CIPSI makes it an ideal tool for treating difficult systems with large static correlations.

Zoom In Zoom Out Reset image size

The FCI space for C₂O₂H₃N in aug-ccpVTZ is approximately 10⁸⁸ determinants. Fortunately, using all of these determinants is not necessary to converge a QMC calculation to chemical accuracy. We truncate determinants based on their magnitudes with a user defined threshold [119], which allows the wavefunction to be evaluated in QMCPACK with a cost growing as $\sqrt[]{N_{\rm det}}$, where N_det is the number of determinants. Truncation values of 10⁻³, 10⁻⁴, 10⁻⁵ and 10⁻⁶ result in wavefunctions of 239, 44 539, 541 380 and 908 128 determinants, respectively. For each truncated wavefunction we optimize one, two and three-body Jastrow factors with VMC. To isolate the improvement of the nodal surface when adding determinants from CIPSI, the coefficients of the determinants were not optimized, although this could result in a further improvement in the wavefunctions. DMC results are extrapolated to a zero time-step using time-steps of $\tau=0.001, ~0.0008$ and 0.0005.

DMC results in table 1 show that the DMC results are converged close to 0.001 Ha, or better than chemical accuracy of 1 kcal mol⁻¹, by around 500 000 determiants. Figure 6 shows the energy as a function of different single-determinant trial wavefunctions as well as multideterminant wavefunctions generated with CIPSI. The latter show a systematic improvement of the nodal surface as a function of the number of determinants. However, it is also interesting to note that in this case the single determinant B3LYP results are quite accurate, highlighting the importance of orbital selection and optimization to improve efficiency [126].

Zoom In Zoom Out Reset image size

Table 1. Energies in Hartree of C₂O₂H₃N as a function of the number of determinants N_det in the trial wavefunction obtained using CIPSI(E), CIPSI(E  +  PT2) and DMC. CIPSI(E) is the variational energy, while CIPSI(E  +  PT2) includes perturbative corrections.

Ndet	CIPSI(E)	CIPSI(E  +  PT2)	DMC
1	−281.6729	−283.0063	$-283.070(1)$
239	−281.7423	−282.9063	$-283.073(1)$
44 539	−282.0802	−282.7339	$-283.078(1)$
541 380	−282.2951	−282.6772	$-283.088(1)$
908 128	−282.4029	−282.6775	$-283.089(1)$

QMCPACK includes the qmca Python-based tool to average quantities in the output files and aid in performing statistical analysis. Given the name of an output file, qmca will compute the average of each quantity in the file. Results of separate simulations can also be aggregated, such as for different twists (twist averaging), multiple steps (autocorrelation analysis) or multiple Jastrow parameters (Jastrow optimization).

In addition to all the quantities computed by QMCPACK, qmca computes the data variance and efficiency. qmca also allows visualizing the evolution of MC quantities over the course of the simulation by a trace offering a quick picture of whether the random walk had expected behavior as in example figure 7.

Zoom In Zoom Out Reset image size

An important step before running a QMC calculation is to obtain the trial wavefunction from another electronic structure or quantum chemical code and convert it into a format readable by QMCPACK. In addition to the large set of converters available through Nexus, QMCPACK comes with two converters upon compilation. Connections to other codes will be developed on request.

convert4qmc: When compiling QMCPACK, an extra binary called convert4qmc is also created. convert4qmc manages gaussian trial wavefunctions from codes such as GAMESS [127], VSVB [128] or quantum package [59]. convert4qmc handles the conversion of single determinant, multideterminants (CASSCF, CI, CIPSI), numerical and Gaussian basis sets. The output file generated can be either in an XML or HDF5 format. convert4qmc allows the user to add multiple option to the wavefunction such as a Jastrow function (two-body, three-body or Pade), cusp conditions, or limit the number of determinants to include.

pw2qmcpack: When using a plane-wave trial wavefunction from the PWSCF code in the QE suite [54, 55], pw2qmcpack.x is used. Source code patches are included with QMCPACK to produce the pw2qmcpack.x binary for specific QE versions, necessary to collect and write the wavefunction in the correct format for QMCPACK.

Black phosphorus (BP) has enormous potential for technological applications because it is layered like graphite and can be exfoliated to form a 2D material that is naturally a semiconductor. There are some technical complications in making devices due to effects such as the degradation of the material when exposed to air or water. Ab initio calculations can help understand such processes and act as a laboratory in which to test mitigation strategies. However, the interaction between the BP layers was thought to be dispersive in nature. Non-covalent interactions are a classically difficult case to treat with density functional based approaches, while, QMC has been shown to perform as well as the most accurate quantum chemical methods while also being applicable to extended and larger systems [132–134]. We employed QMCPACK to calculate the interaction between BP layers as well as the bulk material [135]. We found that in one particular arrangement the interaction between the layers was very well reproduced by one van der Waals (vdW)-corrected DFT functional, but as soon as the stacking between the layers was changed, the functional's performance degraded.

In order to understand this phenomenon, we looked at the change in the charge density induced by the presence of nearby layers of BP. Classically, one might expect very little change in the charge density due to a vdW interaction, and this property is rigorously upheld by various vdW-corrected DFT functionals. However, we found a large charge reorganization caused by the presence of nearby BP layers and of a character that was different than that predicted by DFT.

This study took advantage of QMCPACK's ability to handle mixed boundary conditions section 9. For BP monolayers and bilayers, we used periodic boundary conditions parallel to the layers, but open boundary conditions perpendicular to them. Also, we found the interaction energies were very sensitive to finite size effects, and the ability to perform calculation on many different supercells of the material was crucial to determine their exact form. The B-spline representation was used for the orbitals in the trial wavefunctions section 6.2.1. Finally, the ability of the code to evaluate the charge density from the calculation was crucial.

Carbon can form two-dimensional sheets of sp–sp²-hybridized atoms, α-graphyne. Its existence was predicted three decades ago [136] and has recently received a great deal of attention because of its intriguing potential as a new Dirac material [137–139]. Among various available forms of a graphyne-based structures, a bilayer α-graphyne consisting of two stacked two-dimensional hexagonal lattices can be energetically stabilized in two different stacking modes (AB and Ab mode) out of a total of six available modes; note that the α-graphyne hexagons are much larger than those of graphene. While an AB mode has been predicted to exhibit electronic properties similar to those of an AB bilayer of graphene, Ab-stacked graphyne is expected to possess split Dirac cones at the Fermi level, and its gaps can be opened with and applied electric field normal to the surface [140]. Theoretical predictions of these electronic properties of bilayer graphynes were made by using first-principle DFT methods [140]. However, DFT studies failed to suggest which of the AB and Ab stacking mode is the most energetically stable one because of the approximate description of vdW forces within the non-interacting Kohn–Sham scheme. As seen in table 2, computed interlayer distances and binding energies obtained from DFT calculations are strongly scattered depending on which specific exchange-correlation functional was used.

Table 2. Equilibrium interlayer distance R₀ (Å) and binding energies E_b (meV/atom) for an Ab and AB bilayer α-graphyne using DMC and various vdW-corrected DFT functionals [141]. $\Delta E_{AB-Ab}$ indicates the binding energy difference between AB and Ab α-graphyne.

Method	α-graphyne(AB)		α-graphyne(Ab)		$\Delta E_{AB-Ab}$
	R0	Eb	R0	Eb
DFT-D2	3.25	13.4	3.37	13.6	−0.2
vdW-DF	3.47	19.8	3.64	18.5	1.3
rVV10	3.27	17.9	3.41	17.8	0.1
DMC	3.24(1)	23.2(2)	3.43(2)	22.3(3)	0.9(4)

In order to establish which stacking mode is the most stable one, DMC calculations, that include without any approximations the long-range vdW interactions, were performed to compute the equilibrium interlayer distance and binding energy for Ab and AB bilayer graphyne [141]. DMC equilibrium binding energy for the AB stacking mode (23.2(2) meV/atom) is estimated to be slightly larger than that for Ab stacking mode (22.3(3) meV/atom), which suggests that the AB-stacked bilayer is energetically favored over the Ab stacking mode. In comparing the DMC results with results from vdW-corrected DFT functionals, including the pairwise correction of Grimme (DFT-D2) [142, 143], self-consistent non-local vdW functional (vdW-DF) [144], and simplified non-local vdW-correction (rVV10) [145], it was found that those significantly underestimate the interlayer binding energies for both stacking modes (see figure 8).

Zoom In Zoom Out Reset image size

It is interesting to find which vdW-corrected DFT functional most accurately describe the weak interlayer vdW interaction in bilayer α-graphynes by comparing the DFT results to the DMC ones. As table 2 shows, there is no vdW-corrected DFT functional that achieves good accuracy for both interlayer distance and binding energy. However, the rVV10 functional produces almost identical equilibrium interlayer distances as DMC, and the vdW-DF results for interlayer binding energies and $\Delta E_{AB-Ab}$ are the closest ones to the DMC results for both stacking modes. Therefore, one may conclude that rVV10 can provide well-optimized vdW geometries for low-dimensional carbon allotropes, while vdW-DF gives the best vdW energetics among the various vdW-corrected DFT functionals investigated here.

The row 3 and row 4 transition metals are extremely useful elements for a number of applications, ranging from magnetic applications to solar cells and catalysis. The reason for their versatility is the partially filled d-shell, which gives rise to a number of oxidation and spin states. The 3d electronic states are also rather localized and give rise to relatively strong electronic correlations, especially in the transition metal oxides, such as NiO and TiO₂. This has as a consequence that standard computational approaches using density functional theory within the local density approximation (LDA) or the generalized gradient approximation (GGA) are inadequate and frequently give incorrect results.

TiO₂ occurs in a variety of polymorphs. Three of those occur naturally: rutile, anatase and brookite. Rutile is the most abundant one and is used as a white pigment as well as an opacifier and an ultraviolet radiation absorber, while anatase is the most photocatalytically active polymorph [146, 147].

In spite of the prevalence and many applications of TiO₂, it is very difficult to determine which of the polymorphs is the most stable one. Experimental studies have indicated that rutile is the most stable structure, and anatase and brookite are metastable. However, the enthalpy differences between the polymorphs are very small, of the order of 1 mHa per formula unit, making precise determinations very difficult [148–155].

Electronic structure calculations using a broad range of local and gradient corrected functionals give anatase as the lowest-energy polymorph, and thereafter brookite and rutile [156, 157], in supposed disagreement with experiments. DFT  +  U can give the apparent correct ordering of the polymorphs for a sufficiently large U value [158, 159]. Hybrid functionals can also give the correct energetic ordering of the polymorphs, but only using very large fractions of exact exchange [159].

In two recent studies [93, 160], diffusion Monte Carlo calculations were used to examine TiO₂ polymorphs with the goal of determining the energetic ordering at zero and finite temperatures. Both studies gave the result that anatase is the most stable polymorph at 0 K, with very small energy differences between rutile and brookite. Finite-temperature Helmholtz free energies were then calculated by including phonon contributions based on DFT phonon calculations. The results in both studies were very similar: the finite-temperature contributions from lattice vibrations drove the free energy of the rutile polymorph to be the smallest at temperatures above approximately 650 K (see figure 9). These two studies may finally resolve the question of the most stable polymorphs of TiO₂: at 0 K, anatase is the most stable one, while lattice vibrations drive a transition from anatase to rutile at about 650 K.

Zoom In Zoom Out Reset image size

Because of its many oxidation states, titanium can form a variety of stoichiometrically different oxides. One particular set are the Magnéli phases with the generic formula Ti$_{2n}{\rm O}_{(2n-1)}$. These form ordered crystals with diminishing band gaps with increasing n. In particular, Ti₄O₇ forms a magnetic semiconductor at temperatures below 120 K [161, 162]. DFT finds multiple low energy states, but the ordering depends on the functional used.

DMC calculations [91] found the same energetic ordering of the lowest three states as the DFT calculations by Zhong et al [163], but with larger energy separations. A detailed examination of the spin densities and local moments on the Ti⁴⁺ ions showed that the DFT methods actually gave very good representations of the total moments. However, the orbital alignments were different in DMC, especially in the FM state. This certainly will give rise to different energy differences compared to DFT as the energy differences between the states can rather accurately be attributed to Heisenberg exchange between the magnetic Ti ions, which will differ with the different orbital alignments.

VO₂ has a rich phase diagram which can be exploited in novel device applications. At ambient pressure, the unstrained VO₂ undergoes a metal to insulator transition (MIT) at $T_{\rm c}\approx 341$ K from the low temperature monoclinic M1 phase to the high temperature rutile (R) phase [164–166]. This phase transition has been studied rather extensively, in part due to the on-going debate on the driving mechanism of this MIT. The challenge in describing this material is often related to the representation of electron correlations in this strongly correlated electron material [167–173]. Commonly, experimental studies are accompanied with density functional theory (DFT) calculations for better understanding. Additionally, DFT has been used in isolation to provide insight into the mechanism of the MIT [174–177]. However, the failings of DFT in this context are also well known [178–180]. If functional development led to systematic improvement [181], this should be measurable in both the total energy and the electron density; the two properties that the exact functional must perfectly reproduce.

Recent QMC calculations of VO₂ were used as a reference in assessing various DFT formulations [182]. Supercells of 48 atoms were used to model the antiferromagnetic ground state of the M1 and R phases, and a $3\times 3\times 3$ grid was used in twist-averaging. The QMC calculations used LDA  +  U orbitals for the nodes, where the optimized U-value of 3.5 eV was obtained from DMC as the value that minimized the variational DMC energy.

In general, it was observed that the best description of energetics between the structural phases did not correspond to the best accuracy in the charge density. An accurate spin density was found to lead to a correct energetic ordering of the different magnetic states. However, local, semilocal, and meta-GGA functionals tend to erroneously favor demagnetization of the vanadium sites, which can be reconciled in terms of the self-interaction error. The metrics used also revealed the limitations in the description of correlated 3d-orbital physics present in currently available functionals. This is evident, e.g. in the density metrics shown in figure 10, where spatial variations in the electron density with respect to DMC reference are shown; the extrapolated estimator is used as the DMC reference.

Zoom In Zoom Out Reset image size

Nickel oxide is the poster-child for correlated transition metal oxides. It has a simple rocksalt structure and is at low temperatures a Type II antiferromagnetic insulator. Not only is NiO interesting as a prototypical correlated transition metal oxide, but it is typically a p-type semiconductor, one of a very small number of p-type oxides, because as-deposited NiO typically has a Ni deficiency. Stoichiometric NiO can also be p-doped with suitable monovalent elements, e.g. with Li or K. This brings the basic questions of what are the energetics of vacancies or substitutional dopants in NiO, and how do their descriptions in DFT-based calculations differ from those based on QMC?

DMC calculations [183] were performed on 64-atom supercells with a single substitutional K dopant or Ni or O vacancies. DMC results for ground state properties were in very good agreement with experiments, compared to the DFT results (see table 3). There was also a large difference in defect creation energies, with DFT energies underestimating them table 4. The optical gaps were also underestimated by DFT, as expected. DMC calculations of the optical gap were larger, in fact much larger than experimental values. Calculations of the gap using a 128-atom supercell reduced the gap, but the variance was large enough for this size cell that the reduction of the gap for the larger supercell was not statistically significant. In any case, the gap calculations indicate that finite-size corrections can be significant in DMC gap calculations and also point to the need for improved excited state methodologies in QMC calculations of solids, where the errors are currently larger than for ground states.

Testing pseudopotentials is an important part of QMC due to (i) the historic challenge of developing accurate pseudopotentials for many-body methods such as QMC and historical reuse of DFT or Hartree–Fock-derived potentials, and (ii) the importance of checking any biases unique to QMC such as the T-moves and locality approximations in DMC. Comparing atomic properties such as ionization potentials, and dimer properties such as bond lengths and binding energies obtained by DMC with experimental or quantum chemical results, provides a robust, inexpensive, and transferable test of pseudopotential quality. Recently, QMCPACK and Nexus [129] were used jointly to validate a collection of newly developed pseudopotentials for the early-row transition metals (Sc–Zn) [190].

DMC calculations were performed for five atomic charge states (ranging from neutral to 4+) and at nine transition metal–oxygen dimer bond lengths for each species. Orbitals were generated with QE using the experimental spin multiplicity within LDA [5, 191] or HSE [192, 193] for atoms and LDA for dimers. Two-body Jastrow functions were optimized at each atomic charge state and at the equilibrium geometry for dimers. Subsequent DMC calculations were performed with both T-moves and the locality approximation to assess the affect of pseudopotential localization errors. This large number of calculations is best handled using workflow software such as Nexus.

The resulting DMC atomic ionization potentials and dimer bond lengths, binding energies and vibration frequencies were compared with prior DMC results using Gaussian-based Hartree–Fock [194–196] pseudopotentials, all electron quantum chemistry results using MRCI [197, 198] and CCSD(T) [199–201] approaches, and experiment. On essential all measures, the various theoretical approaches performed similarly well compared with experiment. The current pseudopotential DMC results were within 0.2 eV of experiment on average for atomic ionization potentials and dimer binding energies, with the T-moves and locality approximation generally agreeing to 0.05 eV for these energy differences. Equilibrium bond lengths were found to be within 0.5% of the experimental values, while the more sensitive vibration frequencies agreed to around 3%. This work as well as subsequent studies have verified the quality of the new potentials for QMC studies of transition metal containing systems, but further improvements are desirable section 17.1.

The atomic core energies scale quadratically with the nuclear charge, Z, while the valence energies stay essentially constant across all chemical elements. The energy fluctuations from the core therefore dominate any QMC calculation and, indeed, the overall cost is $\propto Z^a$ with a between 5.5 and 6.5 [202, 203]. It is therefore highly desirable to construct an effective, valence-only Hamiltonian with the atomic cores removed. In fact, a related problem is encountered also in one-particle calculations with plane-waves. For that matter, even in many-body calculations with heavy atoms fully correlating the core(s) would eventually dominate the calculations regardless of the employed method.

Pseudopotentials (PPs) and the closely related effective core potentials (ECPs) have been used in condensed matter calculations as well as in basis set quantum chemical calculations for several decades [1]. At the very basic level, the atomic cores and corresponding degrees of freedom are replaced by PPs/ECPs operators that mimic the action of the core states on the valence electrons. Traditional PP/ECP constructions are based on one-particle solutions of the atom and typically involve norm conservation/shape consistency so that the pseudo-orbitals match the true original all-electron orbitals outside some appropriate core radius. Many advances for ECPs/PPs have been proposed that generalize, make more efficient, improve or more accurately reproduce the true atomic properties, e.g. [204–207]. For QMC calculations, one improvement that has been used to improve the transferability has been based on fitting Hartree–Fock or Dirac–Fock energy differences for a set of atomic and ionic excitations. Furthermore, many-body constructions have been suggested by means of reproducing the (correlated) one-particle density matrices beyond a certain radius or by improving upon DFT solution for a given atom using many-body perturbation theory.

Despite all these elaborate and sophisticated efforts, at present, PPs remain a very laborious part of QMC studies. The main reason is that QMC often reaches accuracy beyond and sometimes well beyond the accuracy of traditional and even currently most advanced PP/ECP constructions. Additional complications come from dealing with the non-locality in the QMC framework that introduces further complications and demands on quality of the trial wavefunctions. As a result, PP/ECP for every element has to be painstakingly retested and benchmarked anew and if the inaccuracies sufficiently bias the valence properties of interest, one has to go back to square one and construct a more accurate PP/ECP. Furthermore, techniques in DFT that use solid state results to improve transferability, e.g. [208], are not practical in QMC because of the computational cost of the large supercells required to converge finite size errors.

To overcome these highly technical but important barriers it is desirable to develop a new generation of PPs/ECPs based on correlated constructions from the outset, provide high accuracy that would not limit the subsequent QMC results, and also be efficiently used in QMC. In addition, one would like to have flexibility in choosing the core-valence partitioning and transferability not only at ambient (equilibrium) conditions but also at high pressures, non-equilibrium conformations and other broader set of conditions.

The goal is therefore not only to reproduce the properties of the atom within one-particle theory and then hope for best in many-body calculation. Our effort would be focused on reproducing the true many-body properties of the original system(s), i.e. atom(s), in a variety of settings. For this purpose, we plan to derive new generation of PPs/ECPs that would be based on a number of new criteria targeted to uphold its accuracy and fidelity to the true original many-body Hamiltonian. In particular,

• (i)  
  we plan to construct an atomic ECP operator that will match a subset of the many-body spectrum as close as possible to the original atomic Hamiltonian. This will involve a set of states that have the largest weights in molecular, surface, solid or other chemical settings;
• (ii)  
  include more options for multiple core-valence partitioning whenever appropriate;
• (iii)  
  express PPs/ECPs in a simple representations/forms that enable their use with multiple methods ranging from traditional DFT and plane wave-based packages to many-body approaches based on stochastic and explicit expansions in basis set methods;
• (iv)  
  try to capture all of the relevant physics that is feasible at the current state-of-the-art, e.g. impact of correlated cores together with correlated valence, describe relativity with the best available account of correlation, explicit treatment of spin–orbit effects, etc.

Clearly this is an ambitious plan that will require significant effort and time, and as such it is almost a never ending task (since it is almost always possible to slightly improve upon the previous version). Nevertheless, we believe that equally important is the adoption of new standards: many-body instead of one-particle framework, testing and benchmarking by a multitude of methods that cross-validate the quality of the PPs/ECPs, and systematic documentation and improvements so that PPs/ECPs can be used without endless retesting and with true many-body quality of the corresponding operators shown upfront. To-date, explicitly correlated ECPs have been developed for a selection of first and second row elements [209], with future developments for the rest of the periodic table underway.

These developments will be collected on a new community website, http://pseudopotentiallibrary.org, that will be used for storing the data and will also include tests and benchmarks. The PPs/ECPs will also be available as a part of the QMCPACK package.

Until now, almost all electronic structure QMC calculations have been carried out with spins of individual electrons being fixed, up or down. This is easy to justify for Hamiltonians without explicit spin operators. Such Hamiltonians commute with any spin, so that up or down orientations are conserved and the solution of the many-body problem is therefore confined to the space of spatial coordinates, and spin states are imposed as a symmetry (e.g. a triplet state). Many interesting phenomena involve the interaction between the spin and spatial degrees of freedom, such as the spin–orbit interaction. Recently, dynamic spins as quantum variables has been realized in DMC in order to treat such interactions directly, see [210] and [211]. We plan to implement dynamic spins into QMCPACK, and a brief outline of the necessary changes are outlined below.

In order to implement spin-dependent operators, spin variables need to be introduced as dynamic quantum variables. This requires a change from one-particle orbitals to one-particle spinors

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi_n({\bf r}_i,s_i) = \alpha \varphi^\uparrow_n({\bf r}_i)\chi^\uparrow(s_i) + \beta \varphi^\downarrow_n({\bf r}_i)\chi^\downarrow(s_i) \label{eqn:spinor} \nonumber \end{align} \tag{ 67 }$$

where the spatial orbitals can be different for the spin up and down channels. The spin functions $\chi^{\uparrow(\downarrow)}$ take on discrete values in the minimal spin representation, namely $\chi^\uparrow(1/2) = \chi^\downarrow(-1/2) = 1$ and $\chi^\uparrow(-1/2) = \chi^\downarrow(1/2) = 0$. The simplest antisymmetric wavefunction is then a single determinant of the one-particle spinors, rather than the product of spin-up and spin-down determinants. In order to increase the efficiency of sampling the spin degrees of freedom, the spinors will be represented using a continuous (overcomplete) representation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \chi^{\uparrow(\downarrow)}(s_i) = {\rm e}^{\pm {\rm i} s_i} \label{eqn:spin_function} \nonumber \end{align} \tag{ 68 }$$

where $s_i \in (0, 2\pi)$. These explicitly varying spins result in complex wavefunctions, for which we must abandon the fixed-node approximation that applies to real-valued wavefunctions. By writing the many-body wavefunction in terms of an amplitude and phase $\newcommand{\e}{{\rm e}} \Psi = \rho \exp\left[ i\Phi \right]$, the real part of the Schrödinger equation becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle -\frac{\partial \rho}{\partial \tau} = \left[ -\frac{1}{2}\nabla^2 + V + \frac{1}{2}|\nabla \Phi|^2 + W^{Re}\right]\rho \nonumber \end{align} \tag{ 69 }$$

where W^Re is ${\rm Re}\left[ \Psi^{-1} W \Psi \right]$ and W is any nonlocal operator such as a PP. Since we do not know the exact phase, we approximate it by the use of a trial phase using the fixed-phase approximation [71]. Since ρ is positive-definite, there is no nodal surface and the DMC algorithm is seemingly the same as for bosonic ground states with an additional potential provided by the trial phase.

Sampling of the spin variables can be achieved by introducing a spin 'kinetic energy', namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T_{s_i} = -\frac{1}{2\mu_s} \left[ \frac{\partial^2}{\partial s_i^2} + 1 \right] \label{eqn:spin_kinetic} \nonumber \end{align} \tag{ 70 }$$

such that it annihilates an arbitrary spinor and does not contribute to the total energy. However, the introduction of this operator modifies the DMC Green's function to include a diffusion and drift term for the spin variables ${\bf S} = (s_1, s_2, \ldots, s_N)$. The $\mu_s$ is a spin mass, which can be interpreted as a time step for the spin variable propagation.

Once the spins are treated as variables rather than static labels, spin-dependent Hamiltonians can be treated. For the spin–orbit integration, a generalization of the non-local operators used in QMC calculations is necessary. Relativistic quantum chemistry calculations utilize the semi-local form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle W_i = \sum\limits_\ell \sum\limits_{j=-|\ell-1/2|}^{\ell+1/2} \sum\limits_{m=-j}^{\,j}W_{\ell j}(r_i) | \ell j m_j \rangle \langle \ell j m_j | \label{eqn:rel_ecp} \nonumber \end{align} \tag{ 71 }$$

for fully relativistic PPs that include scalar-relativistic and spin–orbit effects [212]. Utilizing the spin representation introduced above, the contribution from the non-local PP is inherently complex. When evaluating the localization of the PP with the trial wavefunction ${\rm Re}\left[\Psi_T^{-1}W\Psi_T\right]$, one will encounter terms such as $\newcommand{\e}{{\rm e}} \langle \Omega s | \ell j m \rangle$, where Ω is the solid-angle and s is the spin variable for an individual electron. These terms are given as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle \Omega s | \ell, j = \ell\pm1/2, m_j \rangle = \pm \sqrt{\frac{\ell\pm m + 1/2}{2\ell+1}}Y_{\ell,m-1/2}(\Omega){\rm e}^{{\rm i}s} \nonumber \\ \quad +\sqrt{\frac{\ell \mp m + 1/2}{2\ell + 1}}Y_{\ell,m+1/2}(\Omega){\rm e}^{-{\rm i}s}. \label{eqn:ecp_eval} \nonumber \end{align} \tag{ 72 }$$

By summing over each electron and ion and taking the real part, one obtains the contribution to the local energy from the complex non-local PP in the locality approximation. As is described in [211], a generalization to the Casula T-moves [36, 78] has been constructed for complex non-local operators such as the relativistic PPs above, and will also be included.

Predicting details about future high performance computing architectures is difficult, but some general trends can be observed that should motivate application design decisions in an effort of increasing forward portability. As Moore's Law and Dennard scaling slow due to physical limitations, recent strategies for gaining performance have largely consisted of increasing the parallel capability of hardware. However, the time and power costs of moving data during program execution also limits the traditional increase in parallelism by scaling out commodity-type nodes and connecting them by a high performance interconnect. As a result, parallelism is being added both at the overall system level as before, and also within the processors by increasing thread count, and within the node by adding specialized compute accelerators such as GPUs, FPGAs, etc. Besides increasingly complex compute hierarchies, memory is also becoming hierarchical and more complex, e.g. with the addition of non-volatile memory and various types of on-package high-bandwidth memories. Although the traditional memory structure already involved multi-level caches and system RAM, this complexity was almost entirely managed by the hardware controllers and compilers. However, as the memory hierarchy deepens, it could become necessary for the application developer to more actively manage its usage to realize desired performance increases. Adapting to these changes is very important for QMC because the high computational cost often results in the methods being run on the latest supercomputers with the newest architectures.

While QMC methods have considerable natural parallelism that makes maintaining parallel efficiency fairly straightforward on most traditional hardware, the real desire is to reduce the overall time to solution for any given problem of interest. This means that simply mapping more walkers to more available threads does not scale indefinitely for our purposes, and a multifaceted strategy is needed to leverage increasingly complex modern computer architectures. The emergence of GPU-based HPC architectures around half a decade ago indicated that maintaining high performance was no longer 'business as usual'.

In order to take advantage of the massive increase in the number of lightweight threads that were offered by GPUs, the parallel strategy in QMCPACK was re-mapped from walker-based parallelism to particle-based parallelism. However, as compute accelerators get more powerful with an increasing number of available threads, and more powerful threads have their own vectorization capabilities, additional parallelism must be exploited. For example, so that multiple threads can be assigned to a single particle within each walker when working on compute intensive operations such as the inverse matrix update and Jastrow function calculations. Meanwhile, the natural parallelism present in the Monte Carlo method cannot be abandoned, as it will always be desirable to use more walkers. For this piece of the strategy to continue to be viable, better protocols for equilibration are needed to ensure that the overall time to solution can continue to be reduced.

It is common for HPC application developers to strive for a forward portable design to minimize the effort required in adapting their code to new architectures. With the trend for increasing complexity and a variety of compute and memory architecture configurations, a flexible application design is as important as ever for attaining good 'performance portability'—the ability to run efficiently on significantly different architectures, e.g. conventional processors, GPUs, and potentially even FPGAs.

QMCPACK is being actively developed to adopt to the changing computational landscape. In particular, to (i) run efficiently on systems with fewer 'fat' nodes with deeper heterogeneous compute capabilities as well as systems with a larger number of more homogeneous 'thin' nodes by supporting multiple granularities of computation, (ii) accommodating multiple types of fine-grained multi-threading and vectorization to fully utilize the processors and compute accelerators, and (iii) make effective use of the increasingly complex memory hierarchy. To facilitate this task, simplified 'miniapps' have been developed for real space QMC and AFQMC. These are distributed via https://github.com/QMCPACK. The miniapps encapsulate the main operations of each method while avoiding the full complexity of the main application and enabling use as testbeds for different programming models and middleware layers that help treat the above complexities. In the long term, we hope for a future C++ standard that provides a migration path to an effective and portable solution.