Diffusion Monte Carlo

Diffusion Monte Carlo is implemented in the casino.dmc.DMC class. It projects out the ground state from a trial wave function by evolving an ensemble of walkers (electron configurations) in imaginary time:

\[\Phi_0 = \lim_{t \to \infty} e^{-t(\hat H - E_T)} \Psi\]

where \(E_T\) is the reference (trial) energy. With importance sampling the walkers are distributed according to the mixed distribution \(f(\mathbf{r}, t) = \Psi(\mathbf{r})\Phi(\mathbf{r}, t)\), which satisfies the Fokker–Planck-like equation with drift, diffusion and branching terms. For fermions the nodal surface of \(\Phi\) is constrained to that of \(\Psi\) (the fixed-node approximation [17]): moves that change the sign of \(\Psi\) are rejected.

For a short time step \(\tau\) the importance-sampled Green’s function factorizes into a drift-diffusion part and a branching part:

\[G(\mathbf{r} \to \mathbf{r}', \tau) \approx \underbrace{(2\pi\tau)^{-3N_e/2} \exp\left[-\frac{(\mathbf{r}' - \mathbf{r} - \tau \mathbf{v}(\mathbf{r}))^2}{2\tau}\right]}_{G_D} \cdot \underbrace{\exp\left[\frac{\tau (S(\mathbf{r}) + S(\mathbf{r}'))}{2}\right]}_{G_B}\]

where \(\mathbf{v} = \nabla \ln \Psi\) is the drift velocity and \(S\) is the branching energy. One DMC step consists of a drift-diffusion move for every walker, a branching (birth/death) step, and a T-move step if nonlocal pseudopotentials are present.

Summary of Methods

DMC class has the following methods:

Method	Description
random_walk	perform a number of DMC steps, return the history of \(E_{best}\)
random_step	one DMC step: drift-diffusion, branching, T-move, update of \(E_T\)
drift_diffusion	EBES (method 1) or CBCS (method 2) drift-diffusion move of all walkers
limiting_velocity	drift velocity with UNR limiting applied
alimit_vector	position-dependent limiting parameter \(a(\mathbf{r}_i)\)
branching_energy	ZSGMA-limited branching energy \(S(\mathbf{r})\)
branching	stochastic birth/death of walkers according to their weights
t_move	Casula T-move for nonlocal pseudopotentials
load_balancing	redistribute walkers evenly across MPI ranks

drift_diffusion

Each walker is propagated by the Langevin step

\[\mathbf{r}' = \mathbf{r} + \tau \bar{\mathbf{v}}(\mathbf{r}) + \boldsymbol{\eta}, \qquad \boldsymbol{\eta} \sim \mathcal{N}(0, \tau)\]

where \(\bar{\mathbf{v}}\) is the limited drift velocity. The short-time Green’s function is not exact, so a Metropolis accept/reject step enforces detailed balance:

\[p = \min\left(1, \frac{G_D(\mathbf{r}' \to \mathbf{r}) |\Psi(\mathbf{r}')|^2} {G_D(\mathbf{r} \to \mathbf{r}') |\Psi(\mathbf{r})|^2}\right)\]

Moves crossing the nodal surface (\(\mathrm{sign}\,\Psi(\mathbf{r}') \neq \mathrm{sign}\,\Psi(\mathbf{r})\)) are rejected. To prevent persistent (stuck) walkers, the acceptance probability of a walker that has not moved for more than 50 steps (CBCS) or 20 steps (EBES) is multiplied by \(1.1^{\text{age}}\).

Two sampling modes are implemented:

• CBCS (dmc_method : 2) — all electrons are moved at once, single accept/reject.

• EBES (dmc_method : 1) — electrons are moved one at a time, each with its own accept/reject.

Rejected moves slow the diffusion down, which is accounted for by the effective time step used in branching [18]:

\[\tau_{eff} = \tau \frac{\langle p |\Delta \mathbf{r}|^2 \rangle}{\langle |\Delta \mathbf{r}|^2 \rangle}\]

where \(\Delta \mathbf{r}\) is the diffusive part of a proposed electron displacement, that is the displacement measured from the drifted position, and the averages run over the proposals for all electrons of all walkers, including those rejected at the nodal surface.

Modifications for bare nuclei (nucleus_gf_mods : T, all-electron systems only [18]): the drift of an electron at distance \(z\) from the nearest nucleus is decomposed into components parallel and perpendicular to \(\hat{\mathbf{z}}\), the unit vector from the nucleus to the electron. The drifted position is

\[z'' = \max(z + v_z \tau, 0), \qquad \mathbf{r}_{drift} = z'' \left(\hat{\mathbf{z}} + \frac{2 \mathbf{v}_\rho \tau}{z + z''}\right) + \mathbf{R}_I\]

so that the drift never carries an electron through the nucleus. The new position is then sampled from a mixture of a Gaussian centred at \(\mathbf{r}_{drift}\) and an exponential centred at the nucleus:

\[q = \frac{1}{2}\,\mathrm{erfc}\left(\frac{z + v_z \tau}{\sqrt{2\tau}}\right), \qquad (1 - q)\, \mathcal{N}(\mathbf{r}_{drift}, \tau) + q\, \frac{\zeta^3}{\pi} e^{-2\zeta |\mathbf{r} - \mathbf{R}_I|}, \qquad \zeta = \sqrt{Z^2 + 1/\tau}\]

which mimics the cusp in the Green’s function near the nucleus.

limiting_velocity

Close to the nodal surface the drift velocity diverges, causing large time-step errors. The velocity is smoothly limited following Umrigar, Nightingale and Runge [18]:

\[\bar{\mathbf{v}}_i = \mathbf{v}_i\,\frac{-1 + \sqrt{1 + 2 a v_i^2 \tau}}{a v_i^2 \tau}\]

where \(v_i^2 = |\mathbf{v}_i|^2\) is the squared drift of electron \(i\): every electron is limited by its own velocity, not by that of the whole configuration. For all-electron systems with nucleus_gf_mods : T the parameter \(a\) is position-dependent (alimit_vector):

\[a(\mathbf{r}_i) = \frac{1 + \hat{\mathbf{v}}_i \cdot \hat{\mathbf{z}}_i}{2} + \frac{(Z z_i)^2}{10\,(4 + (Z z_i)^2)}\]

with \(z_i\) the distance to the nearest nucleus; otherwise the constant alimit keyword (default 0.5) is used.

branching_energy

The branching energy of a walker is

\[S(\mathbf{r}) = (E_T - E_{best}) + \bar E_{cut}\,\frac{|\bar{\mathbf{v}}|}{|\mathbf{v}|}\]

where the ratio of the limited to the unlimited drift velocity suppresses branching near the nodal surface [18], and the deviation of the local energy from the current best estimate is limited according to Zen, Sorella, Gillan, Michaelides and Alfè [19]:

\[\bar E_{cut} = \mathrm{sign}(E_{best} - E_L)\, \min\left(|E_{best} - E_L|,\ 0.2\sqrt{N_e/\tau}\right)\]

This scheme makes the time-step error size-consistent: the energy cutoff grows with the number of electrons \(N_e\).

branching

After the drift-diffusion move each walker carries the weight

\[W = \exp\left[\tau_{eff}\left(\left(1 - \frac{p}{2}\right) S(\mathbf{r}) + \frac{p}{2} S(\mathbf{r}')\right)\right]\]

the reweighting factor of UNR [18] averaged over accepting and rejecting the proposed move. Averaging the exponent rather than applying \(\exp[\tau_{eff}(S(\mathbf{r}) + S(\mathbf{r}'))/2]\) to accepted moves and \(\exp[\tau_{eff} S(\mathbf{r})]\) to rejected ones gives the same \(\tau \to 0\) limit with a smaller time-step error, but it needs the local energy at \(\mathbf{r}'\) even when the move is rejected. A move rejected at the nodal surface has \(p = 0\) and the expression reduces to \(W = \exp[\tau_{eff} S(\mathbf{r})]\). In EBES every electron carries its own acceptance probability and the walker has no single \(p\), so there the symmetric average \(W = \exp[\tau_{eff}(S(\mathbf{r}) + S(\mathbf{r}'))/2]\) over the whole electron sweep is used instead. The walker is then replaced by \(M = \lfloor W + u \rfloor\) copies, where \(u\) is a uniform random number on \([0, 1]\); walkers with \(M = 0\) are killed.

The reference energy \(E_T\) is updated after every step to keep the total walker population close to the target (dmc_target_weight) [18]:

\[E_T = E_{best} - \ln\left(\frac{W_{tot}}{W_{target}}\right) \frac{\tau}{\tau_{eff}}\]

where \(E_{best}\) is the population average of the local energy over all MPI ranks.

t_move

The nonlocal pseudopotential energy is evaluated with the localization approximation, which makes the DMC energy non-variational and sensitive to the quality of \(\Psi\). The T-move scheme of Casula [20] (use_tmove : T) restores the upper-bound property by turning the sign-problem-free part of the nonlocal operator into an additional walker move: after branching, each electron of each walker is offered a heat-bath move to one of the angular quadrature points \(\mathbf{r}_q\) (see nonlocal_potential) with probability proportional to

\[-\tau\, V_{nl}(\mathbf{r}_q)\,\frac{\Psi(\mathbf{r}_q)}{\Psi(\mathbf{r})}\]

restricted to the points where this quantity is positive; with the complementary probability the walker stays in place. Electrons are processed one after another (the size-consistent version of the algorithm [21]).

load_balancing

Branching makes the walker population fluctuate independently on every MPI rank. Every 500 steps the walkers are redistributed pairwise between ranks so that each rank carries the same number of walkers. The load-balancing efficiency

\[\eta = \frac{\bar N_{walkers}}{\max_{rank} N_{walkers}}\]

is accumulated and reported at the end of the run.

random_step

One DMC step chains the stages described above:

1. drift_diffusion — propagate and accept/reject every walker;

2. branching — duplicate/kill walkers according to their weights;

3. t_move — T-move for pseudopotential systems;

4. update \(E_{best}\), \(\tau_{eff}\) and \(E_T\) (MPI-reduced over all ranks).

random_walk

Runs steps DMC steps, calling load_balancing every 500 steps, and returns the per-step history of \(E_{best}\) from which the DMC energy and its statistical error are estimated by the reblocking analysis.

References

[17]

P. J. Reynolds, D. M. Ceperley, B. J. Alder, and W. A. Lester, Fixed-node quantum Monte Carlo for molecules, J. Chem. Phys. 77, 5593 (1982).

[18] (1,2,3,4,5,6)

C. J. Umrigar, M. P. Nightingale, and K. J. Runge, A diffusion Monte Carlo algorithm with very small time-step errors, J. Chem. Phys. 99, 2865 (1993).

[19]

A. Zen, S. Sorella, M. J. Gillan, A. Michaelides, and D. Alfè, Boosting the accuracy and speed of quantum Monte Carlo: Size consistency and time step, Phys. Rev. B 93, 241118(R) (2016).

[20]

M. Casula, Beyond the locality approximation in the standard diffusion Monte Carlo method, Phys. Rev. B 74, 161102(R) (2006).

[21]

M. Casula, S. Moroni, S. Sorella, and C. Filippi, Size-consistent variational approaches to nonlocal pseudopotentials: Standard and lattice regularized diffusion Monte Carlo methods revisited, J. Chem. Phys. 132, 154113 (2010).