Lattice Regularized Diffusion Monte Carlo

• LRDMC is a projector quantum Monte Carlo technique that refines guiding wave functions by regularizing the Hamiltonian on a lattice.
• The method projects a prepared trial state toward the fixed-node ground state, controlling discretization bias through lattice spacing a with quadratic error reduction.
• Applications in atomization energies, reactive oxygen species, and all-electron systems demonstrate LRDMC’s effectiveness in evaluating nodal quality and achieving continuum extrapolation.

Lattice Regularized Diffusion Monte Carlo (LRDMC) is a projector quantum Monte Carlo method in which a continuum many-body Hamiltonian is replaced by a lattice-regularized Hamiltonian $\hat{\mathcal H}^a$ satisfying $\hat{\mathcal H}^a \to \hat{\mathcal H}$ as the lattice spacing $a \to 0$, and the resulting Green’s-function dynamics are used to project a guiding wave function toward the fixed-node ground state compatible with its nodal surface. In the literature considered here, LRDMC occupies a specific position within the VMC–DMC hierarchy: VMC prepares and optimizes the trial state, whereas LRDMC performs the beyond-variational projection without altering the trial nodes. The method is described as variational and size-consistent by construction for non-local pseudopotentials, and as a practical realization of diffusion Monte Carlo within a Green’s Function Monte Carlo framework (Casula et al., 2010, Nakano et al., 16 Aug 2025, Raghav et al., 2022).

1. Conceptual role within projector quantum Monte Carlo

The basic logic of LRDMC is the same as that of fixed-node diffusion Monte Carlo: a trial or guiding wave function $\Psi_G$ is projected in imaginary time toward the lowest-energy state compatible with its sign structure. What distinguishes LRDMC from conventional short-time DMC is the discretization strategy. Instead of introducing only a time-step approximation to the propagator, LRDMC regularizes the Hamiltonian itself on a lattice and controls the bias through the lattice spacing $a$, with the continuum limit recovered by extrapolating $a \to 0$ (Nakano et al., 16 Aug 2025, Raghav et al., 2022).

This role is especially clear in application papers. In atomization-energy benchmarks on the G2 set, VMC is used first to optimize flexible many-body ansätze, including nodal parameters, and LRDMC is then applied as the final fixed-node projection step. In acenes, multiple variationally optimized singlet ansätze are compared on equal footing by further projecting each one with fixed-node LRDMC, so that lower LRDMC energy can be interpreted as evidence of better nodes among the tested trial states. In reactive oxygen species and in all-electron Na$_2$, LRDMC is similarly used as the more accurate projector stage after VMC optimization rather than as an autonomous black-box procedure (Raghav et al., 2022, Dupuy et al., 2018, Zen et al., 2014, Nakano et al., 2019).

Within this hierarchy, a recurring theme is that LRDMC is not primarily used to compensate for a poor trial function. Rather, it refines a carefully prepared one. This is why the method is repeatedly discussed together with wave-function optimization, nodal quality, and error cancellation in energy differences, rather than as a stand-alone stochastic solver.

2. Lattice-regularized Hamiltonian and continuum limit

In the formal development revisiting standard and lattice-regularized DMC for non-local pseudopotentials, the lattice Hamiltonian is written as

$\mathcal H^a = \mathcal K^a + \mathcal V^a,$

with

$$\mathcal K^a= -\sum_{i=1}^N \Delta_i^a(\theta_i,\phi_i)/2 + \mathcal V^{NL},$$

and

$$\mathcal V_i^a({\bf R}) = \nu ({\bf r}_i) - \frac{(\Delta_i-\Delta_i^a)\Psi_T({\bf R})}{2\Psi_T({\bf R})}, \qquad \mathcal V^a = \sum_{i=1}^N \mathcal V_i^a + \mathcal V_{ee} + \mathcal V_{nn}.$$

Here the discretized Laplacian $\hat{\mathcal H}^a \to \hat{\mathcal H}$0 is evaluated in a randomly rotated Cartesian frame and satisfies $\hat{\mathcal H}^a \to \hat{\mathcal H}$1. The resulting Hamiltonian obeys

$\hat{\mathcal H}^a \to \hat{\mathcal H}$2

so the leading lattice-space error is quadratic in $\hat{\mathcal H}^a \to \hat{\mathcal H}$3 (Casula et al., 2010).

A distinctive feature of this construction is that the correction operator satisfies

$\hat{\mathcal H}^a \to \hat{\mathcal H}$4

which yields the asymptotic energy structure

$\hat{\mathcal H}^a \to \hat{\mathcal H}$5

The practical implication is that the coefficient of the $\hat{\mathcal H}^a \to \hat{\mathcal H}$6 error becomes small when the trial function is good. In numerical work this motivates polynomial extrapolations such as

$\hat{\mathcal H}^a \to \hat{\mathcal H}$7

which are explicitly used in later application papers for continuum estimates (Casula et al., 2010, Nakano et al., 2019).

In the more recent GFMC-oriented formulation, the discretized kinetic operator is expressed through a finite-difference Laplacian, and the local potential is regularized so that the local energy is preserved: $\hat{\mathcal H}^a \to \hat{\mathcal H}$8 That preservation is presented as important for accelerating convergence to the continuum limit. The same framework also underlies later load-balanced formulations, where the lattice-regularized Hamiltonian remains the starting point and the residual discretization bias is still described as $\hat{\mathcal H}^a \to \hat{\mathcal H}$9 (Nakano et al., 16 Aug 2025).

3. Fixed-node constraint, guiding functions, and nodal quality

Across the benchmark literature, the fixed-node approximation is the central interpretive principle of LRDMC. The projected state is improved relative to the original trial function, but only within the nodal manifold supplied by the guiding wave function. Consequently, lower LRDMC energy among comparable trial states is used as a diagnostic of better nodal surfaces rather than merely better amplitudes (Dupuy et al., 2018, Nakano et al., 2021).

For this reason, the structure of the trial wave function is not ancillary. The common pattern is a Jastrow-multiplied antisymmetric ansatz,

$a \to 0$0

or, in JAGP notation,

$a \to 0$1

Because the Jastrow factor is positive, it does not change the nodes; the nodal surface is determined entirely by the antisymmetric component. Application papers therefore distinguish sharply between standard single-determinant forms such as JSD or JDFT and more flexible pairing-based forms such as AGP, JAGP, or JsAGPs (Zen et al., 2014, Raghav et al., 2022).

The AGP family replaces the occupied-orbital determinant by an antisymmetrized pairing state. In the G2 study, the singlet pairing function is expanded in atomic orbitals and the AGP wave function is constructed as a determinant of the corresponding pairing matrix; for open-shell systems, unpaired orbitals are appended to make the matrix square. The authors emphasize that AGP is multiconfigurational in nature but in practice is as efficient as a single determinant ansatz, and that this added variational freedom matters for LRDMC because it alters the nodal surface on which the projection is carried out (Raghav et al., 2022).

This nodal emphasis is explicit in acenes. There, RVB/JAGP, OSS-oriented JDD, and closed-shell JSD trial states are all optimized at VMC level and then projected by fixed-node LRDMC. The result is not just a variational ordering but a nodal ordering: the RVB/JAGP wave function is reported to have always a lower variational energy and better nodes than the OSS/JDD state for all molecular species considered (Dupuy et al., 2018).

4. Numerical realizations and algorithmic evolution

The LRDMC literature also contains a clear line of algorithmic refinement. In the 2010 revision of variational DMC and LRDMC for non-local pseudopotentials, two concrete improvements were introduced for LRDMC: an improved lattice Hamiltonian regularization with a controlled $a \to 0$2 error and a single-lattice randomization strategy. Earlier implementations had used two incommensurate lattice spacings; the revised scheme instead uses a single spacing $a \to 0$3 while randomizing the orientation of the Cartesian axes for each electron update. This reduces the number of off-diagonal kinetic Green-function elements from $a \to 0$4 to $a \to 0$5, giving a factor-of-two speedup in full-core calculations and a speedup factor $a \to 0$6 with pseudopotentials (Casula et al., 2010).

For all-electron heavy atoms, a different optimization was later proposed through a generalized double-grid LRDMC algorithm. The central idea is to use a small lattice spacing near nuclei and a larger one in the valence region, controlled by a position-dependent function $a \to 0$7. In the improved parameterization, the core radius is chosen as

$a \to 0$8

and the coarse-to-fine spacing ratio scales as

$a \to 0$9

Benchmark timings on all-electron atoms gave $\Psi_G$0 scaling for single-grid LRDMC and $\Psi_G$1 for the double-grid algorithm, with speedups growing from $\Psi_G$2 for He to $\Psi_G$3 for Xe. For benzene, the reported double-grid bias relative to single-grid at the same fine spacing was $\Psi_G$4 mHa, with a speedup of $\Psi_G$5 (Nakano et al., 2019).

A more recent development targets massively parallel hardware. In the conventional continuous-time LRDMC algorithm, walkers propagated over a fixed imaginary-time interval require a random number of off-diagonal moves before branching, so the synchronization overhead grows as $\Psi_G$6. The load-balanced reformulation removes this source of imbalance by sampling only the off-diagonal part of the lattice Hamiltonian and propagating every walker for the same fixed number $\Psi_G$7 of moves between branching steps. On the Leonardo supercomputer with NVIDIA A100 GPUs, this variant maintained approximately $\Psi_G$8 weak-scaling parallel efficiency up to $\Psi_G$9 GPUs, corresponding to $a$0, and gave a direct speedup of $a$1 relative to the conventional algorithm with the same number of walkers (Nakano et al., 16 Aug 2025).

5. Energetics and electronic-structure benchmarks

Representative applications show LRDMC being used chiefly for energetics that are sensitive to nodal quality and error cancellation. The benchmark domains below summarize the pattern seen across molecular applications (Raghav et al., 2022, Dupuy et al., 2018, Nakano et al., 2019, Zen et al., 2014).

Benchmark domain	Use of LRDMC	Representative outcome
G2 atomization energies	Final fixed-node projection after VMC optimization of JDFT and JsAGPs	MAD $a$2 kcal/mol; chemical accuracy for 26 molecules
Oligoacenes $a$3–$a$4	Compare JSD, JDD, and RVB/JAGP nodes	RVB about 5 kcal/mol below JSD; JDD about 2 kcal/mol below JSD
All-electron Na$a$5	Continuum-extrapolated PES from optimized JAGP	$a$6 mHa, $a$7 Å, $a$8 cm$a$9
Reactive oxygen species	Selective energetic refinement on top of VMC/JAGP	Triplet O$a \to 0$0 $a \to 0$1 eV; NO EA $a \to 0$2 eV

In the G2 study, the LRDMC workflow is explicitly post-variational: equilibrium geometries are taken from prior benchmarks, flexible many-body ansätze are optimized at VMC level, and LRDMC is then extrapolated to zero lattice spacing. The central numerical result is a reduction of the atomization-energy mean absolute deviation from about $a \to 0$3 kcal/mol for JDFT to about $a \to 0$4 kcal/mol for JsAGPs, with chemical accuracy achieved for $a \to 0$5 molecules. For the N$a \to 0$6 test case, JsAGPs-LRDMC gives an atomization energy of $a \to 0$7 eV against an experimental and estimated exact value of $a \to 0$8 eV (Raghav et al., 2022).

In oligoacenes, LRDMC serves chiefly as a nodal arbiter. The direct competition between JSD, JDD, and RVB/JAGP singlets preserves the same energetic hierarchy after projection as at VMC level: RVB remains systematically lower than JDD, and both lie below JSD. The average energy gain relative to JSD is reported as about $a \to 0$9 kcal/mol for JDD and about $_2$0 kcal/mol for RVB. Separately, JSD-LRDMC adiabatic singlet–triplet gaps from anthracene to nonacene remain negative, from $_2$1 to $_2$2 kcal/mol, so the singlet stays below the triplet throughout the series (Dupuy et al., 2018).

The Na$_2$3 study illustrates a different use case: all-electron weak binding. VMC-JAGP reproduces the qualitative shape of the potential-energy surface but underbinds strongly, whereas LRDMC on top of the optimized JAGP guide yields $_2$4 mHa, $_2$5 Å, and $_2$6 cm$_2$7, in very good agreement with the experimental references quoted in the paper. The same study also shows the dependence on guiding nodes: at $_2$8 Å, the reported LRDMC binding energies are $_2$9 mHa for JDFT, $\mathcal H^a = \mathcal K^a + \mathcal V^a,$0 mHa for JSD, and $\mathcal H^a = \mathcal K^a + \mathcal V^a,$1 mHa for JAGP (Nakano et al., 2019).

For reactive oxygen species, LRDMC is used selectively in “crucial cases” where VMC/JAGP energetics require refinement. Reported improvements include triplet O$\mathcal H^a = \mathcal K^a + \mathcal V^a,$2 dissociation energy $\mathcal H^a = \mathcal K^a + \mathcal V^a,$3 eV, singlet O$\mathcal H^a = \mathcal K^a + \mathcal V^a,$4 $\mathcal H^a = \mathcal K^a + \mathcal V^a,$5 eV, superoxide $\mathcal H^a = \mathcal K^a + \mathcal V^a,$6 eV, and NO electron affinity $\mathcal H^a = \mathcal K^a + \mathcal V^a,$7 eV. In hydrotrioxyl, the HO–OO dissociation energy changes qualitatively from the wrong sign at VMC, $\mathcal H^a = \mathcal K^a + \mathcal V^a,$8 kcal/mol for trans-HOOO, to $\mathcal H^a = \mathcal K^a + \mathcal V^a,$9 kcal/mol at LRDMC, although this remains far below the $$\mathcal K^a= -\sum_{i=1}^N \Delta_i^a(\theta_i,\phi_i)/2 + \mathcal V^{NL},$$0 kcal/mol CASPT2(19,15) estimate cited in the same study (Zen et al., 2014).

6. Forces and structural observables

Although LRDMC is most often used for total energies, it has also been extended to ionic forces. In this setting the central quantity is the exact derivative

$$\mathcal K^a= -\sum_{i=1}^N \Delta_i^a(\theta_i,\phi_i)/2 + \mathcal V^{NL},$$1

which the force paper decomposes into Hellmann–Feynman, Pulay, and variational contributions. The difficulty specific to DMC/LRDMC is that the projected fixed-node state $$\mathcal K^a= -\sum_{i=1}^N \Delta_i^a(\theta_i,\phi_i)/2 + \mathcal V^{NL},$$2 is not known in closed form, so $$\mathcal K^a= -\sum_{i=1}^N \Delta_i^a(\theta_i,\phi_i)/2 + \mathcal V^{NL},$$3 must be approximated. Two approximations are benchmarked: the Reynolds (RE) approximation, which sets $$\mathcal K^a= -\sum_{i=1}^N \Delta_i^a(\theta_i,\phi_i)/2 + \mathcal V^{NL},$$4, and the variational-drift (VD) approximation, which uses a trajectory-based history term (Nakano et al., 2021).

On six all-electron dimers, LRDMC potential-energy surfaces were reported to improve significantly on VMC. The mean absolute errors with respect to experiment for equilibrium bond lengths and harmonic frequencies are

$$\mathcal K^a= -\sum_{i=1}^N \Delta_i^a(\theta_i,\phi_i)/2 + \mathcal V^{NL},$$5

compared with VMC values

$$\mathcal K^a= -\sum_{i=1}^N \Delta_i^a(\theta_i,\phi_i)/2 + \mathcal V^{NL},$$6

When force-derived parameters are compared with the LRDMC energy-derived reference, RE-LRDMC forces improve on VMC forces: $$\mathcal K^a= -\sum_{i=1}^N \Delta_i^a(\theta_i,\phi_i)/2 + \mathcal V^{NL},$$7 whereas

$$\mathcal K^a= -\sum_{i=1}^N \Delta_i^a(\theta_i,\phi_i)/2 + \mathcal V^{NL},$$8

A hybrid estimator gives $$\mathcal K^a= -\sum_{i=1}^N \Delta_i^a(\theta_i,\phi_i)/2 + \mathcal V^{NL},$$9 Å and $$\mathcal V_i^a({\bf R}) = \nu ({\bf r}_i) - \frac{(\Delta_i-\Delta_i^a)\Psi_T({\bf R})}{2\Psi_T({\bf R})}, \qquad \mathcal V^a = \sum_{i=1}^N \mathcal V_i^a + \mathcal V_{ee} + \mathcal V_{nn}.$$0 cm$$\mathcal V_i^a({\bf R}) = \nu ({\bf r}_i) - \frac{(\Delta_i-\Delta_i^a)\Psi_T({\bf R})}{2\Psi_T({\bf R})}, \qquad \mathcal V^a = \sum_{i=1}^N \mathcal V_i^a + \mathcal V_{ee} + \mathcal V_{nn}.$$1, while VD-LRDMC does not show systematic improvement, yielding $$\mathcal V_i^a({\bf R}) = \nu ({\bf r}_i) - \frac{(\Delta_i-\Delta_i^a)\Psi_T({\bf R})}{2\Psi_T({\bf R})}, \qquad \mathcal V^a = \sum_{i=1}^N \mathcal V_i^a + \mathcal V_{ee} + \mathcal V_{nn}.$$2 Å and $$\mathcal V_i^a({\bf R}) = \nu ({\bf r}_i) - \frac{(\Delta_i-\Delta_i^a)\Psi_T({\bf R})}{2\Psi_T({\bf R})}, \qquad \mathcal V^a = \sum_{i=1}^N \mathcal V_i^a + \mathcal V_{ee} + \mathcal V_{nn}.$$3 cm$$\mathcal V_i^a({\bf R}) = \nu ({\bf r}_i) - \frac{(\Delta_i-\Delta_i^a)\Psi_T({\bf R})}{2\Psi_T({\bf R})}, \qquad \mathcal V^a = \sum_{i=1}^N \mathcal V_i^a + \mathcal V_{ee} + \mathcal V_{nn}.$$4 (Nakano et al., 2021).

A major practical result of that study is the effect of the space-warp coordinate transformation (SWCT) on scaling. Without SWCT, the ratio of force variance to energy variance scales approximately as $$\mathcal V_i^a({\bf R}) = \nu ({\bf r}_i) - \frac{(\Delta_i-\Delta_i^a)\Psi_T({\bf R})}{2\Psi_T({\bf R})}, \qquad \mathcal V^a = \sum_{i=1}^N \mathcal V_i^a + \mathcal V_{ee} + \mathcal V_{nn}.$$5 for LRDMC; with SWCT, the fitted form is $$\mathcal V_i^a({\bf R}) = \nu ({\bf r}_i) - \frac{(\Delta_i-\Delta_i^a)\Psi_T({\bf R})}{2\Psi_T({\bf R})}, \qquad \mathcal V^a = \sum_{i=1}^N \mathcal V_i^a + \mathcal V_{ee} + \mathcal V_{nn}.$$6, effectively independent of $$\mathcal V_i^a({\bf R}) = \nu ({\bf r}_i) - \frac{(\Delta_i-\Delta_i^a)\Psi_T({\bf R})}{2\Psi_T({\bf R})}, \qquad \mathcal V^a = \sum_{i=1}^N \mathcal V_i^a + \mathcal V_{ee} + \mathcal V_{nn}.$$7. The paper therefore concludes that, once SWCT is used, ionic-force evaluation does not worsen the practical atomic-number scaling relative to energy calculations (Nakano et al., 2021).

7. Limitations, misconceptions, and unresolved issues

Several methodological cautions recur across the literature. The most fundamental is that fixed-node LRDMC is only exact if the trial nodes are exact. Lower LRDMC energy therefore identifies better nodes only among the ansätze actually tested. It does not remove the fixed-node constraint itself, and it does not guarantee improved energy differences even when all total energies are lowered (Raghav et al., 2022, Dupuy et al., 2018).

This distinction is visible in atomization benchmarks. In the G2 study, JsAGPs gives lower DMC energies for all atoms and molecules studied, but the authors explicitly note that improved total energies do not always translate into better atomization energies because the relevant quantity depends on cancellation between atomic and molecular fixed-node errors. Si$$\mathcal V_i^a({\bf R}) = \nu ({\bf r}_i) - \frac{(\Delta_i-\Delta_i^a)\Psi_T({\bf R})}{2\Psi_T({\bf R})}, \qquad \mathcal V^a = \sum_{i=1}^N \mathcal V_i^a + \mathcal V_{ee} + \mathcal V_{nn}.$$8H$$\mathcal V_i^a({\bf R}) = \nu ({\bf r}_i) - \frac{(\Delta_i-\Delta_i^a)\Psi_T({\bf R})}{2\Psi_T({\bf R})}, \qquad \mathcal V^a = \sum_{i=1}^N \mathcal V_i^a + \mathcal V_{ee} + \mathcal V_{nn}.$$9 and CO$\hat{\mathcal H}^a \to \hat{\mathcal H}$00 are cited as cases where JsAGPs atomization energies are worse than JDFT despite better individual total energies (Raghav et al., 2022).

A related misconception is that a lower VMC energy or a larger basis necessarily implies better nodes. The reactive-oxygen-species study shows the opposite for the O$\hat{\mathcal H}^a \to \hat{\mathcal H}$01 singlet–triplet gap: basis C gives the best LRDMC result, while larger bases E and F worsen it, even though they are variationally more flexible. The same paper states explicitly that variationally better wave functions at the variational level have not necessarily a better nodal surface (Zen et al., 2014).

Ansatz dependence imposes further limitations. The G2 benchmark emphasizes that AGP is not systematically improvable in the same sense as selected multideterminant wave functions, and that optimization of its many nonlinear parameters can suffer from local minima or instability, with no universal guidelines for constraining or initializing the optimization. The reactive-oxygen-species study adds a separate issue: JAGP is not size-consistent for dissociation into fragments with spin $\hat{\mathcal H}^a \to \hat{\mathcal H}$02 unless the total spin equals the sum of fragment spins. In triplet O$\hat{\mathcal H}^a \to \hat{\mathcal H}$03, this nodal defect persists after projection, leaving the LRDMC dissociation curve about $\hat{\mathcal H}^a \to \hat{\mathcal H}$04 eV too high at $\hat{\mathcal H}^a \to \hat{\mathcal H}$05 Å and about $\hat{\mathcal H}^a \to \hat{\mathcal H}$06 eV too high at $\hat{\mathcal H}^a \to \hat{\mathcal H}$07 Å (Raghav et al., 2022, Zen et al., 2014).

Force calculations introduce their own unresolved problems. The ionic-force benchmark neglects the variational response term in practice, identifying this omission as the source of a self-consistency error that is more severe in LRDMC than in VMC. The same paper also reports that the VD estimator can become prohibitively noisy in all-electron work: for SiH at $\hat{\mathcal H}^a \to \hat{\mathcal H}$08 Å, the ratio $\hat{\mathcal H}^a \to \hat{\mathcal H}$09 is about $\hat{\mathcal H}^a \to \hat{\mathcal H}$10 for VD but only about $\hat{\mathcal H}^a \to \hat{\mathcal H}$11 for RE, leading the authors to call for a new regularization technique (Nakano et al., 2021).

Finally, application papers vary widely in how much of the LRDMC machinery they expose. Some report only the practical setup and defer the full formalism to earlier method papers. In acenes, a single lattice spacing $\hat{\mathcal H}^a \to \hat{\mathcal H}$12 with random orientation is stated to give converged energy differences, and no multi-$\hat{\mathcal H}^a \to \hat{\mathcal H}$13 extrapolation is reported. In the G2 benchmark, the authors state that JsAGPs was followed by LRDMC projection and extrapolation to zero lattice spacing, but they do not provide the actual set of lattice spacings or the fitting formula in the main text (Dupuy et al., 2018, Raghav et al., 2022).

Taken together, these studies portray LRDMC as a mature fixed-node projector method whose strengths are variational control, favorable continuum extrapolation, and strong compatibility with optimized many-body trial states, but whose ultimate accuracy remains inseparable from the nodal quality, size consistency, and optimization stability of the guiding ansatz.