This page solves the Schrödinger equation of a few quantum systems with the variational Monte Carlo (VMC) method: the density |ψ|² of a trial wave function is sampled with the Metropolis algorithm and the energy is estimated as the average of the local energy. Each system is implemented in four languages (Javascript, PHP, Python, C++) with one iron rule: given the same seed, the four backends produce numerically identical results — same random number generator (xorshift64*), same draw order, same stochastic trajectory. The timing comparison thus becomes a controlled experiment: only the language changes (~2 orders of magnitude between C++ and Python). The Javascript backend runs directly in the browser (Web Worker); the others run on the server. Beyond the multi-language suite there are the many-electron systems with a general Slater determinant (Be, Ne, N₂, CH₄) and the general engine with external wave function: the chemistry lives in a JSON file generated once with Hartree-Fock (PySCF), the C++ engine is unique and provides both VMC and fixed-node DMC, the top rung of the HF → VMC → DMC ladder. This tab collects its theoretical companion, with full derivations — from the Schrödinger equation to the variational principle, from Metropolis sampling to fixed-node DMC — and closes with the bibliography.

1. The Schrödinger equation

The starting point is the stationary Schrödinger equation [1] $$ \hat H \,\Psi(\mathbf R) = E \,\Psi(\mathbf R), $$ where \( \mathbf R = (\mathbf r_1, \dots, \mathbf r_N) \) collects the coordinates of the \( N \) electrons. The whole suite works in atomic units (\( \hbar = m_e = e = 4\pi\varepsilon_0 = 1 \)): lengths in bohr (\( a_0 \simeq 0.529\,\text{Å} \)) and energies in hartree (\( 1\,\text{Ha} \simeq 27.211\,\text{eV} \)). In the Born–Oppenheimer approximation the nuclei are fixed and the electronic Hamiltonian reads $$ \hat H = -\frac{1}{2}\sum_{i=1}^{N} \nabla_i^2 \;-\; \sum_{i,A} \frac{Z_A}{|\mathbf r_i - \mathbf R_A|} \;+\; \sum_{i<j} \frac{1}{r_{ij}} \;+\; \sum_{A<B} \frac{Z_A Z_B}{R_{AB}}, $$ with \( r_{ij} = |\mathbf r_i - \mathbf r_j| \). The electron–electron repulsion makes the problem non-separable: the wave function lives in \( 3N \) dimensions and a grid with \( M \) points per dimension would cost \( M^{3N} \) evaluations. This is the curse of dimensionality, and it is why one samples with the Monte Carlo method [10, 13]: the statistical error scales as \( 1/\sqrt{n} \) regardless of dimension.

1.1 Stationary states and unitary evolution (the visualizer)

The rest of the page solves the stationary equation; the “Hydrogen visualizer” tab instead shows the time-dependent equation for the only atom with a closed analytic solution. The hydrogen eigenstates are $$ \psi_{nlm}(\mathbf r) = R_{nl}(r)\, Y_l^m(\theta,\varphi), \qquad E_n = -\frac{1}{2 n^2}\ \text{Ha} = -\frac{13.6}{n^2}\ \text{eV}, $$ with \( R_{nl} \) the radial part (Laguerre polynomials) and \( Y_l^m \) spherical harmonics. A single eigenstate is stationary: \( |\psi|^2 \) does not depend on time. A superposition of them instead evolves unitarily, $$ \Psi(\mathbf r,t) = c_A\,\psi_A\,e^{-iE_A t} + c_B\,\psi_B\,e^{-iE_B t}, $$ and the density \( |\Psi(t)|^2 \) oscillates (“beats”) with angular frequency \( \Delta E = |E_A - E_B| \) and period $$ T = \frac{2\pi}{\Delta E} $$ (atomic units; the visualizer reports it in fs). The colour of each point encodes the local phase \( \arg\Psi \), and the cloud is sampled from \( |\Psi|^2 \) with the same Metropolis algorithm (§2) that, made quantitative, becomes the VMC: it is the visual bridge between stochastic sampling and the computational method of the other tabs.

2. Variational principle and variational Monte Carlo (VMC)

For any normalisable trial function \( \Psi_T \) the variational principle holds: $$ E_V \;=\; \frac{\langle \Psi_T | \hat H | \Psi_T \rangle} {\langle \Psi_T | \Psi_T \rangle} \;\ge\; E_0 , $$ with equality if and only if \( \Psi_T \) is the ground state. The key step of VMC [3, 10] is rewriting the functional as a weighted average over the density \( |\Psi_T|^2 \): $$ E_V = \int \pi(\mathbf R)\, E_L(\mathbf R)\, d\mathbf R, \qquad \pi(\mathbf R) = \frac{|\Psi_T(\mathbf R)|^2}{\int |\Psi_T|^2\, d\mathbf R}, \qquad E_L(\mathbf R) = \frac{\hat H \Psi_T(\mathbf R)}{\Psi_T(\mathbf R)} , $$ where \( E_L \) is the local energy. The average is estimated by sampling \( \pi(\mathbf R) \) with the Metropolis algorithm [2]: from \( \mathbf R \) one proposes \( \mathbf R' = \mathbf R + \delta \) (uniform displacement of half-width step) and accepts with probability $$ A(\mathbf R \to \mathbf R') = \min\!\left(1,\; \left|\frac{\Psi_T(\mathbf R')}{\Psi_T(\mathbf R)}\right|^2\right). $$ The resulting chain satisfies detailed balance \( \pi(\mathbf R)\,P(\mathbf R\to\mathbf R') = \pi(\mathbf R')\,P(\mathbf R'\to\mathbf R) \) and converges to \( \pi \) with no need for the normalisation constant. The first steps (burn_in) are discarded; subsequent samples are serially correlated, so the statistical error is estimated with block averages: the run is split into \( N_b \) blocks (default 100), and $$ \sigma_{E}^2 \simeq \frac{1}{N_b(N_b-1)} \sum_{b=1}^{N_b} \left( \bar E_b - \bar E \right)^2 . $$

The zero-variance principle also holds: if \( \Psi_T \) is an exact eigenfunction, \( E_L(\mathbf R) = E_0 \) is constant and the variance \( \sigma^2(E_L) \) vanishes. Besides the energy, the variance thus measures the quality of the trial function — and it is the absolute correctness test used by Hooke's atom (§4.4).

3. Jastrow factor and Kato cusp conditions

The trial functions of the suite have the product form \( \Psi_T = \Phi \cdot J \), where \( \Phi \) fixes the orbital structure (products or determinants) and \( J = e^{\sum_{i<j} u(r_{ij})} \) is the Jastrow factor [4], here in Padé form: $$ u(r) = \frac{k\, r}{1 + \beta r} . $$ Its job is to correlate electron pairs explicitly and, in particular, to cure the divergences of the Coulomb potential. When two charged particles coincide, \( 1/r \to \infty \): for \( E_L \) to stay finite, the kinetic term must produce an equal and opposite divergence. This yields the Kato cusp conditions [5]: $$ \left.\frac{1}{\Psi}\frac{\partial \Psi}{\partial r_{iA}}\right|_{r_{iA}=0} = -Z_A \qquad\text{(electron–nucleus)}, $$ $$ \left.\frac{1}{\Psi}\frac{\partial \Psi}{\partial r_{ij}}\right|_{r_{ij}=0} = \frac{1}{2} \;\;\text{(antiparallel spins)}, \qquad \frac{1}{4} \;\;\text{(parallel spins)} , $$ where the reduced parallel-spin value follows from antisymmetry (the wave function already vanishes linearly at coalescence). In the Padé form the slope at the origin is \( u'(0) = k \): choosing \( k = 1/2 \) or \( 1/4 \) per pair enforces the exact electron–electron cusp, while \( \beta \) remains the variational parameter setting its range. Electron–nucleus cusps are instead enforced by the orbital exponents (or, in the general engine, by the cusp correction, §8).

4. The six systems of the multi-language suite

4.1 Helium

Padé–Jastrow trial function with exact-cusp \( 1s \) orbitals (\( Z = 2 \)): $$ \Psi_T = e^{-2 r_1}\, e^{-2 r_2}\, \exp\!\left[ \frac{r_{12}}{2\,(1+\alpha\, r_{12})} \right] . $$ The local energy is available in closed form (derivation in Thijssen [12]): $$ E_L = -4 + \frac{(\hat{\mathbf r}_1 - \hat{\mathbf r}_2)\cdot(\mathbf r_1 - \mathbf r_2)} {r_{12}\,(1+\alpha r_{12})^2} - \frac{1}{r_{12}\,(1+\alpha r_{12})^3} - \frac{1}{4\,(1+\alpha r_{12})^4} + \frac{1}{r_{12}} . $$ With \( \alpha = 0.15 \) one gets \( E = -2.87756 \pm 0.00061 \) Ha, against the exact value \( -2.90372 \) Ha [15].

4.2 H₂⁺ ion

One electron and two centres at distance \( R \): LCAO \( \phi = e^{-a r_A} + e^{-a r_B} \). The exact cusp on both nuclei requires accounting, at coalescence on A, for the contribution of the orbital centred on B: the condition gives the fixed-point equation $$ a = 1 + e^{-aR} , $$ solved by iteration. The total energy (with \( 1/R \)) has its minimum around \( R \simeq 2.0\text{–}2.2 \) bohr with \( E \simeq -0.576 \) Ha, qualitatively correct for a minimal LCAO (one-electron bonding is demonstrated).

4.3 Molecular hydrogen H₂

Two electrons in the same bonding orbital, for each \( \phi(i) = e^{-a r_{iA}} + e^{-a r_{iB}} \) with the same fixed-point cusp, plus the Jastrow with \( k = 1/2 \) (antiparallel spins). With \( L_i = \nabla_i^2 \phi_i / \phi_i \) and \( \mathbf F_i = \nabla_i \phi_i / \phi_i \), the local kinetic energy is $$ \text{KE} = -\frac{1}{2}\Big[ L_1 + L_2 + 2u'' + \frac{4u'}{r_{12}} + 2u'\,(\mathbf F_1 - \mathbf F_2)\cdot \hat{\mathbf r}_{12} + 2u'^2 \Big] , $$ where the \( u \) terms come from \( \nabla^2 e^{u} / e^{u} = u'' + 2u'/r + u'^2 \) applied to both electrons and from the cross product of orbital and Jastrow gradients. With \( \beta = 0.7 \) the minimum falls at \( R = 1.40 \) bohr with \( E = -1.1513 \pm 0.0009 \) Ha: 98% of the exact value \( -1.1745 \) Ha of Kołos and Wolniewicz [16].

4.4 Hooke's atom (a zero-variance case)

Two electrons in a harmonic trap with Coulomb repulsion: $$ \hat H = -\frac{1}{2}(\nabla_1^2 + \nabla_2^2) + \frac{\omega^2}{2}(r_1^2 + r_2^2) + \frac{1}{r_{12}} . $$ With \( \Psi_T = e^{-\omega (r_1^2+r_2^2)/2}\, e^{u(r_{12})} \) the local energy is $$ E_L = 3\omega - u'' - \frac{2u'}{r_{12}} + \omega\, u' r_{12} - u'^2 + \frac{1}{r_{12}} . $$ Taut [9] showed that for particular \( \omega \) the solution is analytic: for \( \omega = 1/2 \) the ground state has correlation factor \( J = 1 + r_{12}/2 \) and energy exactly \( E = 2 \) Ha. In exact mode the suite reproduces \( E = 2.0000000000 \) with zero variance: every configuration yields the same local energy, the most severe correctness test a VMC code can pass (§2).

4.5 H₃⁺ ion

Two electrons and three protons at the vertices of an equilateral triangle of side \( d \) in the \( xy \) plane (circumradius \( d/\sqrt3 \)). Three-centre LCAO with cusp from the fixed point \( a = 1 + 2e^{-ad} \) (at coalescence on one nucleus the other two contribute), Jastrow with \( k=1/2 \), nuclear repulsion \( 3/d \). Minimum at \( d \simeq 1.65\text{–}1.9 \) bohr with \( E \simeq -1.248 \) Ha (93% of the reference \( -1.3438 \) Ha [17]), below the H₂ energy: the stability of the ion — one of the most abundant in the interstellar medium — is demonstrated.

4.6 Lithium hydride (LiH)

Four electrons (2↑, 2↓): antisymmetry is enforced by separate \( 2\times2 \) Slater determinants per spin, \( \Phi = D^{\uparrow} D^{\downarrow} \), with orbitals $$ \phi_1 = e^{-a_1 s}, \qquad \phi_2 = e^{-a_2 t} + c\, e^{-a_3 s} , $$ where \( s \) and \( t \) are the distances from Li and H. The determinant is invariant under \( \phi_2 \to \phi_2 + \lambda\, \phi_1 \) (column operations), a property used as a test. The Jastrow is pairwise with spin-dependent cusp (\( k = 1/4 \) parallel, \( 1/2 \) antiparallel, §3). At the defaults \( E = -7.988 \pm 0.003 \) Ha, essentially the Hartree–Fock limit (\( -7.987 \)); the exact non-relativistic value is \( \simeq -8.070 \) Ha: the difference is the correlation energy that a single determinant cannot capture.

5. General Slater determinant (Be, Ne, N₂, CH₄)

For many-electron systems the suite uses a general \( N\times N \) determinant machine per spin species: given the Slater matrix $$ A_{ij} = \phi_j(\mathbf r_i) , $$ the modulus of the determinant is evaluated as \( \ln|\det A| \) via LU decomposition with pivoting (stable also near the nodes) and the inverse by Gauss–Jordan. Differentiating \( \ln D \) is inverse-matrix algebra: $$ \nabla_i \ln D = \sum_j (A^{-1})_{ji}\, \nabla \phi_j(\mathbf r_i), \qquad \frac{\nabla_i^2 D}{D} = \sum_j (A^{-1})_{ji}\, \nabla^2 \phi_j(\mathbf r_i) . $$ With the pairwise Jastrow, the local kinetic energy per electron assembles as $$ \text{KE}_i = T_i + \sum_{j \ne i}\Big( u'' + \frac{2u'}{r_{ij}} \Big) + 2\, \mathbf d_i \cdot \mathbf j_i + |\mathbf j_i|^2 , $$ where \( T_i \) is the pure determinantal term, \( \mathbf d_i = \nabla_i \ln D \) and \( \mathbf j_i = \nabla_i \ln J \). The code structure is identical for the four systems; only the orbitals \( \phi_j \) change (Slater-type with variational exponents, molecular symmetries fixed: D_∞h for N₂, T_d for CH₄). The engines also report the Hartree–Fock and exact references for comparison.

6. The general engine with external wave function

The architectural leap of Stage C is separating the engine from the chemistry: the wave function becomes a JSON file generated once by a restricted Hartree–Fock (RHF) calculation with PySCF [19], and the C++ engine is unique: $$ \Psi_T = D^{\uparrow}[\phi_k]\; D^{\downarrow}[\phi_k]\; J, \qquad \phi_k = \sum_{\mu} C_{k\mu}\, \chi_\mu , $$ where the \( \chi_\mu \) are contracted Gaussian atomic orbitals of \( s \) and \( p \) type: $$ \chi_\mu(\mathbf r) = P_\mu(\mathbf u) \sum_i c_i\, e^{-a_i u^2}, \qquad \mathbf u = \mathbf r - \mathbf R_c , $$ with \( P_\mu = 1 \) (type \( s \)) or \( u_x, u_y, u_z \) (type \( p \)). Defining \( S_n = \sum_i c_i\, a_i^n\, e^{-a_i u^2} \), the derivatives needed by the local energy are in closed form: $$ s:\;\; \chi = S_0, \quad \nabla\chi = -2 S_1 \mathbf u, \quad \nabla^2\chi = 4 S_2 u^2 - 6 S_1 ; $$ $$ p:\;\; \chi = u_m S_0, \quad \nabla\chi = -2 S_1 u_m \mathbf u + S_0\, \hat{\mathbf e}_m, \quad \nabla^2\chi = u_m (4 S_2 u^2 - 10 S_1) . $$ Validation is three-tiered, strongest to weakest: (i) pointwise: the engine's local energy matches, to \( \sim 10^{-6} \), the numerical Laplacian of \( \ln\Psi \) computed independently; (ii) statistical: with the Jastrow off (\( \beta = 0 \)) and cusps disabled, VMC reproduces the PySCF RHF energy within error bars; (iii) physical: with the Jastrow on, the energy drops below HF while staying above the exact value (variational bound, §2).

6.1 The molecule set

The library comprises 18 wave functions: the 7 baseline ones (H₂O and CH₄ in STO-3G and 6-31G, ethane, ethanol, benzene) plus 11 added to cover, each, a distinct chemical concept. All are closed-shell single-determinant RHF; unless stated otherwise the basis is STO-3G.

Molecule	Concept illustrated
F2	The emblematic DMC case: at Hartree–Fock level F2 comes out almost unbound; VMC and DMC recover the correct binding energy — a live demonstration of why correlation matters.
C2H4 (ethene)	C=C double bond (sp2 hybridisation).
C2H2 (acetylene)	C≡C triple bond (sp); together with ethane (single C–C, sp3) it completes the sp3→sp2→sp series.
(H2O)2 (dimer)	Hydrogen bond: a weak interaction (~5 kcal/mol) that HF underestimates and DMC recovers.
NH3 (ammonia)	Lone pair and umbrella inversion.
H2CO (formaldehyde)	C=O carbonyl group.
CH3OH (methanol)	Lighter companion to ethanol, for series comparison.
CO2	Linear molecule with cumulated double bonds (D∞h symmetry).
HF (hydrogen fluoride)	Strongly polar diatomic, with a well-known "exact" reference.
C5H5N (pyridine)	Heteroaromatic: like benzene but with a nitrogen in the ring — same computational cost, one extra heteroatom.
C10H8 (naphthalene)	Extends the aromatic series beyond benzene (68 electrons): a testbed for the method's scalability.

6.2 What this engine cannot generate

The trial function is a single closed-shell RHF determinant times the Jastrow, all-electron with cusp correction. This choice — simple and robust — leaves out four families of systems:

• Multireference / biradical (C₂, O₃, many bond-breaking species): a single determinant gives qualitatively wrong nodes and therefore a large fixed-node error (§10). A multi-determinant wave function would be needed.
• Open-shell (radicals, O₂ in its triplet ground state): the generator enforces zero total spin (RHF); ROHF/UHF orbitals and different-sized determinants for the two spins would be required.
• Bases with d or f functions (6-31G*, cc-pVDZ, …): the engine only evaluates \( s \)- and \( p \)-type orbitals (§6).
• Heavy elements (third row onward, transition metals): all-electron, the deep cores blow up cost and local energy; pseudopotentials would be needed, and are not implemented.

In all these cases the route is not "add a molecule" but extend the wf format and the engine (multi-determinant, ROHF/UHF, pseudopotentials): a far larger undertaking than simply exporting a JSON.

7. Single-electron moves and the Sherman–Morrison formula

Moving all electrons at once and redoing the LU would cost \( O(N^3) \) per move. The engine moves one electron at a time: displacing \( \mathbf r_i \to \mathbf r_i' \) changes a single row of \( A \), and the determinant ratio is read off the maintained inverse with no recomputation: $$ \rho = \frac{\det A'}{\det A} = \sum_j \phi_j(\mathbf r_i')\,(A^{-1})_{ji} . $$ If the move is accepted, the inverse is updated in \( O(N^2) \) with the Sherman–Morrison formula [18]: $$ S_k = \sum_l \phi_l(\mathbf r')\,(A^{-1})_{lk} \qquad (S_i = \rho), $$ $$ (A^{-1})'_{jk} = (A^{-1})_{jk} - \frac{(A^{-1})_{ji}\, S_k}{\rho} \;\; (k \ne i), \qquad (A^{-1})'_{ji} = \frac{(A^{-1})_{ji}}{\rho} . $$ A "sample" becomes a sweep (one proposed move per electron) plus one energy evaluation from the maintained state, with per-electron caches of molecular-orbital values, gradients and Laplacians: no decomposition in the hot loop. For numerical hygiene the inverse is rebuilt from scratch every refresh sweeps (default 250). Measured efficiency gain (\( \sigma^2 \times t \) metric) on H₂O: 19× over all-electron moves.

8. Cusp correction for Gaussian basis sets

Gaussians have zero derivative at the nucleus: they cannot satisfy the Kato cusp (§3), so \( E_L \) diverges like \( -Z/r \) near the nuclei and the variance blows up. The remedy is the scheme of Ma, Towler, Drummond and Needs [11] (here in a simplified version): within a radius \( r_c \) around each nucleus, the \( s \) part of each molecular orbital is replaced by $$ R(r) = \operatorname{sgn} \cdot e^{\,p(r)}, \qquad p(r) = \alpha_0 + \alpha_1 r + \alpha_2 r^2 + \alpha_3 r^3 + \alpha_4 r^4 , $$ with five constraints: the exact cusp \( p'(0) = -Z_{\text{eff}} \), the \( C^2 \) matching at \( r_c \) (continuity of \( p, p', p'' \)) and the central value \( \alpha_0 \). The delicate point — discovered at the expense of the H₂O variance, which initially got worse — is that \( \alpha_0 \) must not be pinned to the HF value: it must be left free and optimised by minimising the deviation of the orbital local energy $$ E_L^{\text{orb}}(r) = -\frac{1}{2}\left( p'' + p'^2 + \frac{2p'}{r} \right) - \frac{Z_{\text{eff}}}{r} $$ inside the sphere, so that the polynomial stays tame (on the oxygen core \( \alpha_4 \) goes from \( +2\cdot10^4 \) to \( -78 \)). With the full recipe the H₂O variance collapses from 966 to 39.5 Ha² (24×). The fitting happens in the exporter; the engine reads the ready-made coefficients from the JSON. Cusp correction and electron–nucleus Jastrow would impose the cusp twice: the engine rejects the combination.

9. Variational optimisation of β with correlated sampling

The Jastrow parameter \( \beta \) must be optimised: the variational principle guarantees that the minimum of \( E_V(\beta) \) is the best choice. Re-running a VMC for each \( \beta \) would be expensive and noisy (different minima due to statistical noise); one uses correlated sampling instead: sample once at \( \beta_{\text{ref}} \) and reweight $$ E(\beta) \simeq \frac{\sum_n w_n E_L^{(\beta)}(\mathbf R_n)}{\sum_n w_n}, \qquad w_n = \left| \frac{\Psi_\beta(\mathbf R_n)}{\Psi_{\beta_{\text{ref}}}(\mathbf R_n)} \right|^2 = e^{\,2[U_\beta(\mathbf R_n) - U_{\beta_{\text{ref}}}(\mathbf R_n)]} , $$ where the determinantal parts cancel in the ratio: the \( T \) and \( \mathbf d \) terms already computed are reused and only the Jastrow terms change — the whole \( E(\beta) \) curve comes out of a single pass. The effective-sample fraction \( w_{\text{eff}} = (\sum w_n)^2 / (n \sum w_n^2) \) signals when the overlap between the two distributions degrades. Two practical devices: the reference is moved iteratively onto the argmin (with parabolic refinement), and the weight exponent is recentred on the first sample — in large systems \( \Delta U \) has a huge constant part that saturates the numerical clamp and would make all weights equal (a bogus \( w_{\text{eff}} = 1 \)), but weights are defined up to a constant which cancels in the ratio. A non-obvious result: the better the basis, the more harm small \( \beta \) does (the long-range Jastrow distorts the one-electron density), so the optimal \( \beta \) grows with basis quality (H₂O: 7.51 in STO-3G, 4.66 in 6-31G).

10. Fixed-node diffusion Monte Carlo (DMC)

VMC is limited by the form of \( \Psi_T \). DMC [6, 7, 20] projects out the excitations by solving the Schrödinger equation in imaginary time \( \tau = it \): $$ -\frac{\partial \Phi}{\partial \tau} = (\hat H - E_T)\,\Phi \qquad\Longrightarrow\qquad \Phi(\tau) = e^{-(\hat H - E_T)\tau}\, \Phi(0) \;\xrightarrow[\tau\to\infty]{}\; c_0\, \Psi_0 , $$ because each eigencomponent decays as \( e^{-(E_n - E_T)\tau} \) and the slowest to die is the ground state. Sampling \( \Phi \) directly is inefficient (and impossible where it changes sign); one samples the mixed distribution \( f = \Psi_T\, \Phi \), which satisfies $$ -\frac{\partial f}{\partial \tau} = -\frac{1}{2}\nabla^2 f + \nabla\cdot\left[ \mathbf v(\mathbf R)\, f \right] + \left[ E_L(\mathbf R) - E_T \right] f , \qquad \mathbf v = \nabla \ln |\Psi_T| . $$ The three terms are read as stochastic processes on a population of walkers: Gaussian diffusion of width \( \sqrt{\tau} \), drift along \( \mathbf v \) (pushing walkers where \( \Psi_T \) is large), and branching with weight $$ w = \exp\!\left[ -\tau \left( \frac{E_L(\mathbf R) + E_L(\mathbf R')}{2} - E_T \right) \right] , $$ realised as \( \lfloor w + u \rfloor \) copies (\( u \) uniform). The reference \( E_T \) is updated to stabilise the population (default 400 walkers) and the DMC energy is the weighted average of \( E_L \). A Metropolis acceptance on the approximate Green's function corrects the finite-time-step error.

Fixed-node. For fermions \( \Phi \) must change sign, but \( f = \Psi_T \Phi \) must stay positive: one imposes that \( \Phi \) has the same nodes as \( \Psi_T \) [20, 7], rejecting node-crossing moves (the determinant sign is tracked through the LU pivots). The FN-DMC energy is exact for the given nodes and remains an upper bound: the residual error is the nodal error alone, as seen by comparing H₂O with STO-3G nodes (\( -76.311 \)) and 6-31G nodes (\( -76.426 \pm 0.008 \), 12 mHa from the exact \( -76.438 \)).

Stability: the limited drift. Near the nodes \( \mathbf v \sim 1/d \) diverges: a walker proposes huge jumps, is always rejected and, stuck with a huge \( E_L \), clones itself until the population collapses (the pathology observed at small \( \tau \)). The cure by Umrigar, Nightingale and Runge [8] replaces the drift with the limited version $$ \tilde{\mathbf v} = \mathbf v\; \frac{-1 + \sqrt{1 + 2\tau |\mathbf v|^2}}{\tau |\mathbf v|^2} , $$ which tends to \( \mathbf v \) for \( \tau |\mathbf v|^2 \to 0 \) and stays \( O(\sqrt{2/\tau}) \) for \( |\mathbf v| \to \infty \), used consistently in the proposal and in both Green's functions; in the branching \( \tau \) is replaced by \( \tau_{\text{eff}} = \tau \cdot \langle \text{acceptance} \rangle \) and the \( E_L \) tails are cut at \( \pm 2/\sqrt{\tau} \). After the cure, \( \tau = 0.005 \) and \( 0.0025 \) agree (the time-step error is under control) and the acceptance grows monotonically as \( \tau \) decreases.

11. Cross-language numerical identity and the random number generator

The iron rule of the suite — identical results in C++, Python, Javascript and PHP at equal seed — rests on the xorshift64* generator of Marsaglia and Vigna [14]: $$ s \mathrel{\oplus{=}} s \gg 12, \qquad s \mathrel{\oplus{=}} s \ll 25, \qquad s \mathrel{\oplus{=}} s \gg 27, \qquad \text{out} = s \times \text{0x2545F4914F6CDD1D} , $$ with the uniform in \( [0,1) \) built from the top 53 bits divided by \( 2^{53} \) (exactly a double's mantissa). The traps lie in languages without native 64-bit integers: Python masks with & 0xFFFFFFFFFFFFFFFF at every step; Javascript keeps the state as a pair of 32-bit words and multiplies by 16-bit limbs; PHP emulates the unsigned logical shift. Same draw order everywhere ⇒ same stochastic trajectory, down to the MD5 of the walker list. In the parallel engine the seeds of chains and walkers are derived with splitmix64 from (seed, step, index): the result is bit-identical with any number of threads — the same property that keeps the future GPU port reproducible (counter-based scheme à la Philox, as in QMCPACK [21]).

12. Bibliography

1. E. Schrödinger, Quantisierung als Eigenwertproblem (Erste Mitteilung), Annalen der Physik 384 (79), 361–376 (1926). doi:10.1002/andp.19263840404
2. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, Equation of State Calculations by Fast Computing Machines, J. Chem. Phys. 21, 1087–1092 (1953). doi:10.1063/1.1699114
3. W. L. McMillan, Ground State of Liquid He⁴, Phys. Rev. 138, A442–A451 (1965). doi:10.1103/PhysRev.138.A442
4. R. Jastrow, Many-Body Problem with Strong Forces, Phys. Rev. 98, 1479–1484 (1955). doi:10.1103/PhysRev.98.1479
5. T. Kato, On the eigenfunctions of many-particle systems in quantum mechanics, Commun. Pure Appl. Math. 10, 151–177 (1957). doi:10.1002/cpa.3160100201
6. D. M. Ceperley, B. J. Alder, Ground State of the Electron Gas by a Stochastic Method, Phys. Rev. Lett. 45, 566–569 (1980). doi:10.1103/PhysRevLett.45.566
7. P. J. Reynolds, D. M. Ceperley, B. J. Alder, W. A. Lester Jr., Fixed-node quantum Monte Carlo for molecules, J. Chem. Phys. 77, 5593–5603 (1982). doi:10.1063/1.443766
8. C. J. Umrigar, M. P. Nightingale, K. J. Runge, A diffusion Monte Carlo algorithm with very small time-step errors, J. Chem. Phys. 99, 2865–2890 (1993). doi:10.1063/1.465195
9. M. Taut, Two electrons in an external oscillator potential: Particular analytic solutions of a Coulomb correlation problem, Phys. Rev. A 48, 3561–3566 (1993). doi:10.1103/PhysRevA.48.3561
10. W. M. C. Foulkes, L. Mitas, R. J. Needs, G. Rajagopal, Quantum Monte Carlo simulations of solids, Rev. Mod. Phys. 73, 33–83 (2001). doi:10.1103/RevModPhys.73.33
11. A. Ma, M. D. Towler, N. D. Drummond, R. J. Needs, Scheme for adding electron–nucleus cusps to Gaussian orbitals, J. Chem. Phys. 122, 224322 (2005). doi:10.1063/1.1940588
12. J. M. Thijssen, Computational Physics, 2nd ed., Cambridge University Press (2007). ISBN 978-0-521-83346-2.
13. B. L. Hammond, W. A. Lester Jr., P. J. Reynolds, Monte Carlo Methods in Ab Initio Quantum Chemistry, World Scientific, Singapore (1994). doi:10.1142/1170
14. G. Marsaglia, Xorshift RNGs, J. Stat. Softw. 8(14), 1–6 (2003). doi:10.18637/jss.v008.i14; S. Vigna, An experimental exploration of Marsaglia's xorshift generators, scrambled, ACM Trans. Math. Softw. 42(4), 30 (2016). doi:10.1145/2845077
15. C. L. Pekeris, 1¹S and 2³S States of Helium, Phys. Rev. 115, 1216–1221 (1959). doi:10.1103/PhysRev.115.1216
16. W. Kołos, L. Wolniewicz, Potential-Energy Curves for the X¹Σg⁺, b³Σu⁺, and C¹Πu States of the Hydrogen Molecule, J. Chem. Phys. 43, 2429–2441 (1965). doi:10.1063/1.1697142
17. W. Cencek, J. Rychlewski, R. Jaquet, W. Kutzelnigg, Sub-microhartree accuracy potential energy surface for H₃⁺ … I. Calculation of the potential points, J. Chem. Phys. 108, 2831–2836 (1998). doi:10.1063/1.475702
18. J. Sherman, W. J. Morrison, Adjustment of an Inverse Matrix Corresponding to a Change in One Element of a Given Matrix, Ann. Math. Statist. 21, 124–127 (1950). doi:10.1214/aoms/1177729893
19. Q. Sun et al., PySCF: the Python-based simulations of chemistry framework, WIREs Comput. Mol. Sci. 8, e1340 (2018). doi:10.1002/wcms.1340
20. J. B. Anderson, A random-walk simulation of the Schrödinger equation: H₃⁺, J. Chem. Phys. 63, 1499–1503 (1975). doi:10.1063/1.431514
21. J. Kim et al., QMCPACK: an open source ab initio quantum Monte Carlo package…, J. Phys.: Condens. Matter 30, 195901 (2018). doi:10.1088/1361-648X/aab9c3

Note on numbering: in-text citations use [1]=Schrödinger, [2]=Metropolis, [3]=McMillan, [4]=Jastrow, [5]=Kato, [6]=Ceperley–Alder, [7]=Reynolds et al., [8]=Umrigar et al., [9]=Taut, [10]=Foulkes et al., [11]=Ma et al., [12]=Thijssen, [13]=Hammond et al., [14]=Marsaglia/Vigna, [15]=Pekeris, [16]=Kołos–Wolniewicz, [17]=Cencek et al., [18]=Sherman–Morrison, [19]=PySCF, [20]=Anderson, [21]=QMCPACK.