---
abstract: 'We present a way to include non local potentials in the standard Diffusion Monte Carlo method without using the locality approximation. We define a stochastic projection based on a fixed node effective Hamiltonian, whose lowest energy is an upper bound of the true ground state energy, even in the presence of non local operators in the Hamiltonian. The variational property of the resulting algorithm provides a stable diffusion process, even in the case of divergent non local potentials, like the hard-core pseudopotentials. It turns out that the modification required to improve the standard Diffusion Monte Carlo algorithm is simple.'
author:
- Michele Casula
title: Beyond the locality approximation in the standard diffusion Monte Carlo method
---

Diffusion Monte Carlo (DMC) is one of the most successful methods to compute the ground state properties of quantum systems. Although the fixed node (FN) approximation is needed to cure the infamous sign problem for fermions, the accuracy of the DMC framework has yielded many benchmark results[@foulkesreview]. However, when the DMC method is applied to “ab initio” realistic Hamiltonians, its computational cost scales $\propto Z^{6.5}$, where $Z$ is the atomic number[@ceperley86]. Therefore, the use of pseudopotentials is necessary to make those calculations feasible.

Since the pseudopotentials are usually non local, the “locality approximation” is made besides the FN, by replacing the true Hamiltonian $H$ with an *effective* one $H^{\mathrm{eff}}$, which reads[@mitas]: $$H^{\mathrm{eff}} = K + V_{\mathrm{loc}} + \frac{\int dx^\prime \langle
      x^\prime | V_{\mathrm{non~loc}} | x \rangle \Psi_T(x^\prime)}{ \Psi_T(x)}, 
\label{H_locality}$$ where $K$ is the kinetic operator, $V_{\mathrm{loc}}$ is the local potential, and the last term in Eq. \[H\_locality\] is the non local potential localized by means of the trial wave function $\Psi_T$. The projection is then realized by iteratively applying the operator $G=\exp(-\tau (H^{\mathrm{eff}} - E_{\mathrm{eff}}))$ to $\Psi_T$ in order to filter out its high energy components. The localized potential enters in the branching part (birth and death process) of the algorithm, while the usual FN constraint is employed to limit the diffusion process within the nodal pockets of $\Psi_T$, and avoid the fermionic sign problem. Thus $E_{\mathrm{eff}}$ is the FN ground state energy of $H^{\mathrm{eff}}$, computed during the sampling of the mixed distribution $\Psi_{\mathrm{eff}} \Psi_T$: $$E_{\mathrm{eff}}=\frac{\langle \Psi_{\mathrm{eff}} | H^{\mathrm{eff}} | \Psi_T
  \rangle}{\langle \Psi_{\mathrm{eff}} | \Psi_T \rangle} = 
\frac{\langle \Psi_{\mathrm{eff}} | H | \Psi_T
  \rangle}{\langle \Psi_{\mathrm{eff}} | \Psi_T \rangle} = E_{MA}.
\label{E_MA_identity}$$ $E_{MA}$ is the *mixed average* of $H$, and the above identity holds because $H^{\mathrm{eff}} \Psi_T/\Psi_T = H \Psi_T/\Psi_T$. Since $\Psi_{\mathrm{eff}}$ is the FN ground state of $H^{\mathrm{eff}}$, which differs from $H$, $E_{MA}$ is no longer equal to the variational FN energy of $H$, defined as: $$E_{FN}=\langle \Psi_{\mathrm{eff}} | H | \Psi_{\mathrm{eff}} \rangle /
\langle \Psi_{\mathrm{eff}} | \Psi_{\mathrm{eff}} \rangle.
\label{E_FN}$$ Therefore, in contrast with the case of local Hamiltonians, $E_{MA}$ calculated with the locality approximation does not in general give an upper bound to the ground state energy of $H$ (variational principle).

In a previous work[@lrdmc], we introduced the Lattice Regularized Diffusion Monte Carlo algorithm (LRDMC), which provides an upper bound for the true ground state energy and allows estimate $E_{FN}$, even in the case of non local potentials. In this paper we propose an extension of the standard DMC framework that gives the same results as the LRDMC method, after a proper modification of the DMC propagator.

We start by considering the importance sampling Green function $$G(x^\prime \leftarrow x, \tau) = \frac{\Psi_T(x^\prime)}{\Psi_T(x)} \langle
x^\prime | e^{-\tau (H - E_T)} | x \rangle,  
\label{green_dmc}$$ where $E_T$ is an energy offset, $\tau$ the time step, and $x$ a vector of particle coordinates. In the diffusion Monte Carlo method, $G(x^\prime
\leftarrow x, \tau)$ is iteratively applied to $\Psi_T^2$, in order to sample stochastically the mixed distribution $\Phi(x,t)=\Psi_T(x)
\Psi(x,t)$, $\Psi(x,t)$ converging to the lowest possible state in energy. To rewrite $G(x^\prime \leftarrow x,\tau)$ (Eq. \[green\_dmc\]) in a practical way, it is necessary to resort to the Trotter break up, which is exact in the limit of $\tau
\rightarrow 0$. Here we split the Hamiltonian into local and non local operators, and we end up with the following expression for the Green function: $$G(x^\prime \leftarrow x,\tau) \simeq \int dx'' ~ T_{x^\prime,x''}(\tau) ~
G_{DMC}(x'' \leftarrow x, \tau),
\label{new_green_function}$$ where $G_{DMC}(x^\prime \leftarrow x, \tau)$ is the usual DMC propagator[@foulkesreview], $$\frac{1}{(2 \pi \tau)^{\frac{3N}{2}}} 
\exp\left[-\frac{(x^\prime - x - \tau v(x))^2}{2 \tau}\right]
e^{-\tau(E^{\mathrm{loc}}_L(x^\prime)- E_T)},
\label{DMC_green_function}$$ and $T_{x^\prime,x}(\tau)$ is the matrix containing the non local potential, $$\frac{\Psi_T(x^\prime)}{\Psi_T(x)} \langle x^\prime | e^{-\tau
V_{\mathrm{non~loc}}}| x \rangle 
  \simeq   \delta_{x^\prime,x} - \tau V_{x^\prime,x}.$$ In the above Eqs. $N$ is the total number of particles, $v(x)=\nabla \ln
|\Psi_T(x)| $ the drift velocity, $E^{\mathrm{loc}}_L(x)=(K+V_{\mathrm{loc}})\Psi_T(x)/\Psi_T(x)$ the contribution to the local energy coming from the local operators, and $V_{x^\prime,x}=\frac{\Psi_T(x^\prime)}{\Psi_T(x)} \langle
x^\prime | V_{\mathrm{non~loc}} | x \rangle$. The final form of $G_{DMC}$ has been obtained by further splitting the Hamiltonian into the kinetic and potential part, while the exponential of the non local potential in $T$ has been linearized up to order $\tau$.

If the case of pseudopotentials, the number of non-zero matrix elements $V_{x^\prime,x}$ will be *finite*, once a quadrature rule with a discrete mesh of points is applied to evaluate the projection over the angular components of the pseudopotential[@fahy; @mitas]. Therefore, the process in $G(x^\prime \leftarrow x,\tau)$ driven by $T_{x^\prime,x}(\tau)$ can be calculated using a heat bath algorithm, since $T_{x^\prime,x}(\tau)
/ \sum_{x''} T_{x'',x}(\tau)$ can be seen as a transition probability, and it can be computed *a priori* for all possible new coordinates $x^\prime$. We notice that the matrix elements $T_{x^\prime,x}(\tau)$ are easily evaluated in a standard DMC algorithm, since $V_{x^\prime,x}$ are already computed to calculate the localized pseudopotential in Eq. \[H\_locality\]: $$\frac{\int dx^\prime \langle x^\prime | V_{\mathrm{non~
      loc}} | x \rangle \Psi_T(x^\prime)}{ \Psi_T(x)} = \sum_{x^\prime} V_{x^\prime,x}. 
\label{potential_locality}$$ At variance with the locality approximation, $ V_{x^\prime,x}$ contribute now to move the particles, according to the transition matrix $T$ ($T$-moves).

An important limitation of this idea is given by the sign problem. Indeed both $\frac{\Psi_T(x^\prime)}{\Psi_T(x)}$ and $\langle
x^\prime | V_{\mathrm{non~loc}} | x \rangle$ can change sign, which should be included in the weights, but this yields averages with exponentially increasing noise. A solution is to apply the FN approximation not only to $G_{DMC}$ but also to $T$, which becomes: $$T^{FN}_{x^\prime,x}(\tau) = \delta_{x,x^\prime} - \tau V^-_{x^\prime,x},
\label{T_FN}$$ where we defined $V^\pm_{x^\prime,x}=1/2 (V_{x^\prime,x} \pm | V_{x^\prime,x} |)$. In practice, we keep only those matrix elements which give a positive $T_{x^\prime,x}(\tau)$. Moreover, we add to the *diagonal* potential the so called “sign flip term”, i.e. the sum over the discarded matrix elements $V^+_{x^\prime,x}$. Therefore, the local potential becomes $$V^{\textrm{eff}}(x)=V_{\mathrm{loc}}(x) + \sum_{x^\prime} V^+_{x^\prime,x}.$$ This is equivalent to work with a new effective FN Hamiltonian $$\begin{aligned}
\label{H_effective}
H^{\textrm{eff}}_{x,x} & = & K + V^{\textrm{eff}}(x) \\
H^{\textrm{eff}}_{x^\prime,x} & = & \langle x^\prime | V_{\mathrm{non~loc}} | x \rangle  
\textrm{~~~if $V_{x^\prime,x}<0$}.\nonumber \end{aligned}$$

In contrast to the effective Hamiltonian of the locality approximation written in Eq. \[H\_locality\], the ground state energy $E_{\mathrm{eff}}(=E_{MA})$ of the above $H^{\textrm{eff}}$ is an upper bound for the ground state energy of the true $H$. As shown in Ref.  for the Lattice Green function Monte Carlo, this variational property is due to the sign flip term (*positive* contribution) added to the local potential, *and* the $T$-moves driven by the off diagonal matrix elements $V^-_{x^\prime,x}$. Instead, in the locality approximation also $V^-_{x^\prime,x}$ is summed in the diagonal part (Eq. \[potential\_locality\]), and this leads to an *attractive* potential, which cannot provide a variational property for $E_{MA}$. Moreover, we found that the negative divergences of the fully localized potential on the nodes of $\Psi_T$ are responsible in some case (e.g. see Fig. \[noise\]) for numerical instabilities in the locality approximation, which disappear once $H^{\textrm{eff}}$ in Eq. \[H\_effective\] is used together with the $T^{FN}$-moves. Indeed, whenever $V^-_{x^\prime,x}$ is large, it pushes the walker away from the attractive regions of the localized potential, and protects the sampling from divergences in the weights.

Once a $T^{FN}$-move is generated according to the transition probability $T^{FN}_{x^\prime,x}(\tau) / \sum_{x''}
T^{FN}_{x'',x}(\tau)$, the walker should acquire the weight $w_T(x,\tau)=\sum_{x''}
T^{FN}_{x'',x}(\tau)$ due to the normalization of the $T^{FN}$ matrix. This weight can be recast as an exponential form valid up to order $\tau$, $$w_T=1-\tau \sum_{x^\prime} V^-_{x^\prime,x} \simeq \exp\left[-\tau \sum_{x^\prime}
  V^-_{x^\prime,x}\right].
\label{T_weight}$$ Thus the overall weight $w(x,\tau)$ of $G(x^\prime \leftarrow x,\tau)$ will be $$w(x,\tau) = w_{DMC}~ w_T = \exp \left[-\tau (E_L(x) - E_T)
  \right], 
\label{final_weight}$$ where $w_{DMC}$ is the weight of $G_{DMC}$ for the effective Hamiltonian ($V_{\mathrm{loc}}$ replaced by $V^{\textrm{eff}}$), and $E_L(x)=H^{\mathrm{eff}} \Psi_T/\Psi_T = H \Psi_T/\Psi_T$ is the local energy. Notice that a non-symmetric branching factor has been included in $G_{DMC}$ (Eq. \[DMC\_green\_function\]). When we use the exponential form in Eq. \[T\_weight\], and consequently the weight in Eq. \[final\_weight\], the time step error is usually smaller than that obtained with the linear form. This can be understood in the limit of perfect importance sampling. Indeed, if $\Psi_T$ is close to the ground state of $H$, the weight in Eq. \[final\_weight\] is almost constant, since the variance of $E_L(x)$ is small, and the time step bias is reduced.

The proposed DMC scheme for fixed node Hamiltonians with non local potentials is the following: (i) perform a diffusion-drift move according to $G_\textrm{diff}(x^\prime \leftarrow
x, \tau)= \exp\left[-(x^\prime - x - \tau v(x))^2 /2 \tau\right]/(2 \pi
\tau)^{\frac{3N}{2}} $ *as is done* in the standard DMC algorithm, and accept or reject this move according to the probability $$min\left[1,\frac{G_\textrm{diff}(x \leftarrow x^\prime, \tau)
      \Psi_T^2(x^\prime)}{G_\textrm{diff}(x^\prime \leftarrow x, \tau)
      \Psi_T^2(x)}\right];
\label{acc_rej}$$ (ii) weight the walker with the factor $\exp \left[-\tau (E_L(x^\prime) - E_T) \right]$; (iii) displace the walker a second time, with a $T$-move selected according to the transition probability $p(x'' \leftarrow x^\prime, \tau)=
T^{FN}_{x'',x^\prime}(\tau) / \sum_y T^{FN}_{y,x^\prime}(\tau)$, computed *a priori* for all possible new $x''$. The branching process will be the same as in the usual DMC algorithm. In practice, only the $T$-move is the new step, which is performed *after* weighting the walker[@time_step_error]. Although we perform an acceptance/rejection step (Eq. \[acc\_rej\]), which has been shown to reduce the time step error[@cyrus_dmc] in $G_\textrm{DMC}$, the algorithm does not satisfy exactly the detailed balance except in the limit of $\tau \rightarrow 0$, due to the break up of $G$ into $G_\textrm{DMC}$ and $T^{FN}$ (Eq. \[new\_green\_function\]), and the use of a non symmetric branching factor in Eq. \[DMC\_green\_function\].

In order to estimate the variational FN energy $E_{FN}$ (Eq. \[E\_FN\]), and study the quality of the locality approximation, we introduce a more general effective Hamiltonian[@lrdmc] $H^{\alpha,\gamma}$, $$\begin{aligned}
\label{H_last}
H^{\alpha,\gamma}_{x,x} & = & K + V_{\mathrm{loc}}(x) + (1+\gamma) \sum_{x^\prime}
  V^+_{x^\prime,x} \nonumber \\
 & & + \alpha (1+\gamma) \sum_{x^\prime}  V^-_{x^\prime,x} \\
H^{\alpha,\gamma}_{x^\prime,x} & = & -\gamma ~ \langle x^\prime |
V_{\mathrm{non~loc}} | x \rangle  \textrm{~~~~~~~~~~~~~~~~~if $V_{x^\prime,x}>0$} \nonumber \\
H^{\alpha,\gamma}_{x^\prime,x} & = & (1-\alpha(1+\gamma))~ \langle x^\prime |
V_{\mathrm{non~loc}} | x \rangle  \textrm{~~if $V_{x^\prime,x}<0$}, \nonumber \end{aligned}$$ where $0 \le \alpha \le 1$ and $0 \le \gamma \le 1/\alpha -1$ are two external parameters. In order to sample the Green function $G(x^\prime \leftarrow x, \tau)$ for $H^{\alpha,\gamma}$, it is sufficient to modify the matrix $T_{x^\prime,x}(\tau)$, which becomes $$T^{\alpha,\gamma}_{x^\prime,x} = \left \{
\begin{array}{ll}
1     & \textrm{if $x = x^\prime$} \label{T_matrix_new} \\
\tau ~ \gamma ~  V^+_{x^\prime,x}  & \textrm{if $V_{x^\prime,x}>0$} \\ 
- \tau ~ (1-\alpha(1+\gamma)) ~ V^-_{x^\prime,x} & \textrm{if $V_{x^\prime,x}<0$}. 
\end{array} 
\right.$$ The ground state $E(\alpha,\gamma)$ of $H^{\alpha,\gamma}$ is equal to $E_{MA}(\alpha,\gamma)$ (Eq. \[E\_MA\_identity\]), since $H^{\alpha,\gamma} \Psi_T/\Psi_T = H\Psi_T/\Psi_T$ by construction. The Hamiltonian in Eq. \[H\_effective\] is recovered with $\alpha=0$ and $\gamma=0$, while the Hamiltonian of the locality approximation (Eq. \[H\_locality\]) is obtained with $\alpha=1$ and $\gamma=0$. Therefore, $H^{\alpha,\gamma}$ can interpolate between these two extremes, but the variational principle for $E_{MA}(\alpha,\gamma)$ is not guaranteed as soon as $\alpha \ne 0$, since the attractive term $\alpha (1+\gamma) \sum_{x^\prime}  V^-_{x^\prime,x}$ is added to the diagonal potential. However by means of $H^{\alpha,\gamma}$ one can estimate the value of $E_{FN}(\alpha,\gamma)$ (Eq. \[E\_FN\]), which is *variational* for every $\alpha$ and $\gamma$, since it is the expectation value of the true $H$ on the ground state of $H^{\alpha,\gamma}$. Indeed $H=H^{\alpha,\gamma}-(1+\gamma) \partial_\gamma H^{\alpha,\gamma}$, and the Hellmann-Feynman theorem leads to the relation $$E_{FN}(\alpha,\gamma)=E(\alpha,\gamma)-(1+\gamma)~ \partial_\gamma E(\alpha,\gamma). 
\label{evaluate_E_FN}$$ One can show [@sandro_effective] that, for a given value of $\alpha$, the lowest $E_{FN}(\alpha,\gamma)$ is obtained for $\gamma=0$. Therefore, in order to find the best variational estimate of the ground state of $H$, it is enough to calculate the expression in Eq. \[evaluate\_E\_FN\] with $\gamma=0$. In this way one can check which $\alpha$ provides the best variational state for $H$. The derivative $\partial_\gamma E(\alpha,0)$ can be computed with either finite differences or correlated sampling. In both cases, one should keep in mind that $\gamma < 1/\alpha -1$, to guarantee the positivity of the $T^{\alpha,\gamma}$ matrix (Eq. \[T\_matrix\_new\]), and so calculating $E_{FN}(\alpha,0)$ becomes harder as $\alpha$ gets closer to $1$.

Here we present the application of the method to the $Si$ and $C$ pseudoatoms. We computed $E_{MA}(\alpha,0)$ and $E_{FN}(\alpha,0)$ for $\alpha=0,0.5,0.9$, and the DMC energy with the locality approximation, which corresponds to $E_{MA}(1,0)$. With an aim to quantify the locality error, and the correction provided by the effective Hamiltonian $H^{\alpha,\gamma}$, we used three different trial wave functions (with no Jastrow, a two-body, and a three-body (electron-electron-ion) Jastrow factor respectively), sharing the same determinantal part, and hence the same nodes. In this way, the FN error can be separated from the effect of the locality approximation, which causes a dependence of the DMC energy on the *shape* of $\Psi_T$.

For the $Si$ atom we used an $s-p$ norm-conserving Hartree-Fock (HF) pseudopotential, which is soft and has been generated using the Vanderbilt construction[@vanderbilt]. The determinantal part of $\Psi_T$ is a HF wave function with a $6s6p/1s1p$ Gaussian basis set. The 2-body Jastrow is from Ref. , while the 3-body Jastrow factor is from Ref. . Both of the Jastrow factors have been optimized using an energy minimization procedure[@sorella_hessian]. The results are reported in Fig. \[Si\_data\_alpha\]. The variational $E_{FN}(\alpha,0)$ improves going from $\alpha=0$ to $\alpha=0.9$, i.e. approaching the locality approximation. It means that at least for this soft pseudopotential the locality approximation ($\alpha=1,\gamma=0$) gives a ground state which is a good variational wave function for $H$. Notice however that the standard DMC energies $E_{MA}(1,0)$ have a sizable locality error, while $E_{FN}$ with $\alpha=0.9$ depends only slightly on the shape of the trial wave function. A similar result was obtained with the LRDMC method for the same pseudoatom[@lrdmc].

For the $C$ atom we chose to work with an SBK pseudopotential[@SBK], which is extremely hard, since it diverges like $1/r^2$ in the $s$ channel, and $1/r$ in its local component. The Slater part of $\Psi_T$ is an antisymmetrized geminal power (AGP) wave function[@casula_atoms; @casula_wf] with a $2s2p$ Gaussian basis set, optimized in the presence of the 3-body Jastrow factor by minimizing its variational energy[@sorella_hessian]. The determinantal part has been kept fixed in the other two $\Psi_T$’s, which differ only by their Jastrow factors. The results are plotted in Fig. \[C\_data\_alpha\]. Here the locality approximation is very poor, as it leads to non variational $E_{MA}$. The spikes in Fig. \[noise\], coming from regions of the configuration space where the *effective* potential is attractive, are surely responsible of the non variational results. Surprisingly, $\Psi_T$ without Jastrow, which has a higher energy, leads to much more stable DMC simulations. The locality approximation, which relies on the quality of the shape of $\Psi_T$ in the core, performs poorly with this hard-core pseudopotential, since it is difficult to find the optimal shape of $\Psi_T$ in the core region, due to the divergence of the non local pseudopotential. Indeed $E_{FN}(\alpha,0)$ is higher for $\alpha=0.9$, being the worst for the 3-body Jastrow factor. On the other hand, the best variational $E_{FN}(\alpha,0)$ is obtained for $\alpha=0$, irrespective of the form of the Jastrow factor.

To summarize, we have described a scheme to treat non local potentials within the standard DMC method. We have extended the DMC formalism to handle a generic Hamiltonian with discrete off-diagonal matrix elements and the fixed node approximation. Only a simple modification of the standard algorithm is required to include the $T$-moves generated according to the non local potentials. By using an effective Hamiltonian approach, we showed that it is possible to have stable simulations, even in the case of divergent hard-core pseudopotentials, and obtain variational results. A similar effective Hamiltonian has been successfully used in the LRDMC method. The difference is in the kinetic part, which is discretized in the lattice regularized approach. The LRDMC and the DMC methods have the same efficiency for small $Z$, although it is possible to have a gain in the LRDMC efficiency by an *ad hoc* choice of the kinetic parameters, particularly for heavier elements. We conclude, by noting that the same Green function presented here can be used in the Reptation Quantum Monte Carlo method[@baroni99].

This work was supported by the NSF grant DMR-0404853. We thank S. Sorella, C. Umrigar, D. M. Ceperley, S. Chiesa, and J. Kim for useful discussions.

[99]{} W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. [**73**]{}, 33 (2001). D. M. Ceperley, J. Stat. Phys. [**43**]{}, 815 (1986); A. Ma, N. D. Drummond, M. D. Towler, and R. J. Needs, Phys. Rev. E [**71**]{}, 066704 (2005). L. Mitáš, E. L. Shirley and D. M. Ceperley, J. Chem. Phys. [**95**]{}, 3467 (1991). M. Casula, C. Filippi, and S. Sorella, Phys. Rev. Lett. [ **95**]{}, 100201 (2005). S. Fahy, X. W. Wang and Steven G. Louie, Phys. Rev. B [**42**]{}, 3503 (1990). D. F. B. ten Haaf, H. J. M. van Bemmel, J. M. J. van Leeuwen, W. van Saarloos, and D. M. Ceperley, , 13039 (1995). Since the normalization of $T^\textrm{FN}_{x',x}$ (Eq. \[T\_weight\]) depends on $x$, the best choice is to weight the walker before the $T$-move. P. J. Reynolds, *et al.*, J. Chem. Phys. [**77**]{}, 5593 (1982); C. J. Umrigar, *et al.*, J. Chem. Phys. [**99**]{}, 2865 (1993). S. Sorella, cond-mat/0201388. D. Vanderbilt, Phys. Rev. B [**32**]{}, 8412 (1985). C. Filippi and C. J. Umrigar, J. Chem. Phys. [**105**]{}, 213 (1996). M. Casula, C. Attaccalite, and S. Sorella, J. Chem. Phys. [**121**]{}, 7710 (2004). S. Sorella, Phys. Rev. B [**71**]{}, 241103 (2005). W. J. Stevens, H. Basch, and M. Krauss, J. Chem. Phys. [**81**]{}, 6026 (1984). M. Casula, and S. Sorella, J. Chem. Phys. [**119**]{}, 6500 (2003). S. Baroni, and S. Moroni, Phys. Rev. Lett. [**82**]{}, 4745 (1999).
