Variational Monte Carlo

Variational Monte Carlo is implemented in the casino.vmc.VMC class. It samples the probability density \(|\Psi(\mathbf{R})|^2\) of a trial wave function with the Metropolis algorithm and estimates the expectation value of the Hamiltonian as the mean of the local energy over the resulting chain. The Hamiltonian enters only through that average: the walk itself is driven by \(|\Psi|^2\) alone.

Three sampling modes are named by vmc_method, of which two are implemented:

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

• DBDS (vmc_method : 2) — one spin determinant is displaced at a time. Not implemented, and refused by the input check;

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

The single free parameter of the walk is the step size, set by dtvmc. The formulas below are written in \(\tau\), the standard deviation of one displacement component, so that dtvmc, which is by convention the variance of that component, is \(\tau^2\). The step is too small when successive configurations are strongly correlated and too large when almost every move is rejected, and the efficiency of the walk is flat enough between those extremes that any acceptance from roughly 30 % to 70 % costs little [23]; opt_dtvmc : T tunes it to 50 % in optimize_vmc_step, starting from the closed-form estimate in approximate_step_size. The rest of this page derives that estimate, because the derivation also fixes what the optimizer can and cannot be expected to find, and where the remaining empirical constant comes from.

Notation

Write \(\Lambda(\mathbf{R}) = \ln|\Psi(\mathbf{R})|\) for the logarithm of the trial function, \(\mathbf{F} = \nabla\Lambda\) for the drift vector and \(\mathsf{H} = \nabla\nabla\Lambda\) for its Hessian, all in the full \(3N_e\)-dimensional configuration space.

Two averages appear below and must not be confused:

Symbol	Averaged over	Held fixed
\(\mathbb{E}_\Delta[\,\cdot\,|\,\mathbf{R}]\)	the proposed displacement \(\boldsymbol{\Delta}\)	the configuration \(\mathbf{R}\), so \(\Lambda\), \(\mathbf{F}\), \(\mathsf{H}\) are constants
\(\mathbb{E}_\mathbf{R}[\,\cdot\,]\)	the chain, with density \(|\Psi|^2 / \int |\Psi|^2\)	nothing

The full average is \(\mathbb{E} = \mathbb{E}_\mathbf{R}\,\mathbb{E}_\Delta\), innermost first, and that is the order the algorithm performs: the walker sits at \(\mathbf{R}\) and draws \(\boldsymbol{\Delta}\), and only the chain distributes \(\mathbf{R}\) itself. The order matters, because the acceptance probability is a non-linear function of the configuration and the two stages do not commute. The derivation below commutes them anyway — it averages the local kinetic energies over the chain and only then applies the non-linear acceptance function, rather than the reverse — and the whole of the last section is the price of doing so. That price is paid by the estimate of the step size and by nothing else: the walk itself never commutes anything and samples \(|\Psi|^2\) exactly at any step size.

Three local kinetic energies are used. They are functions of \(\mathbf{R}\), hence constants with respect to \(\mathbb{E}_\Delta\):

\[\begin{split}\begin{aligned} T_L(\mathbf{R}) &= -\frac{1}{2}\frac{\nabla^2\Psi}{\Psi} \\ T_D(\mathbf{R}) &= \frac{1}{2}\,\mathbf{F}\cdot\mathbf{F} = \frac{1}{2}\sum_i |\nabla_i\Lambda|^2 \\ T_M(\mathbf{R}) &= \frac{1}{2}\left(T_L + T_D\right) \end{aligned}\end{split}\]

The identity \(\nabla^2\Psi/\Psi = \nabla^2\Lambda + |\nabla\Lambda|^2\) ties them to the Laplacian of the logarithm,

\[\nabla^2\Lambda = -2\,(T_L + T_D) = -4\,T_M\]

These are the three kinetic-energy estimators Casino reports: \(T_D\) is FISQ, \(T_M\) is TI, and \(T_L = 2\,T_M - T_D\) is KEI, the one entering the total energy. In Pycasino \(T_L\) is what wfn_kinetic_energy returns, and \(T_D\) is computed inside it as \(\mathbf{F}\cdot\mathbf{F}/2\).

The local energy \(E_L = T_L + V\) plays no part in what follows, and neither does the potential. The Metropolis criterion contains \(|\Psi|^2\) only. The kinetic energies appear as combinations of derivatives of \(\Lambda\), which is an identity, not physics; \(V\) is not differentiated anywhere.

Finally, \(\Phi\) and \(\varphi\) denote the standard normal distribution function and density. Every inverse of \(\Phi\) below is written as an inverse error function, \(\Phi^{-1}(p) = \sqrt{2}\,\mathrm{erfinv}(2p - 1)\), the form the code uses, so that a target acceptance \(A\) enters through \(\mathrm{erfinv}(1 - A)\) rather than through \(\Phi^{-1}(1 - A/2)\), and

\[s \equiv \sqrt{2}\,\mathrm{erfinv}(1/2) = 0.6744898, \qquad \varphi(s) = 0.3177766\]

What is accepted

CBCS displaces all electrons at once, \(\mathbf{R}' = \mathbf{R} + \boldsymbol{\Delta}\), where the \(3N_e\) components \(\delta_k\) are independent Gaussians of standard deviation \(\tau\). The proposal is symmetric, so the Metropolis ratio reduces to the density ratio alone:

\[\begin{split}\begin{aligned} a(\mathbf{R} \to \mathbf{R}') &= \min(1, e^X) \\ X = \ln\frac{|\Psi(\mathbf{R}')|^2}{|\Psi(\mathbf{R})|^2} &= 2\left[\Lambda(\mathbf{R}') - \Lambda(\mathbf{R})\right] \end{aligned}\end{split}\]

Everything reduces to the statistics of the single scalar \(X\). The mean acceptance is a two-stage average, in the order the algorithm performs it — the walker sits at \(\mathbf{R}\) and draws \(\boldsymbol{\Delta}\), and only then is \(\mathbf{R}\) itself distributed according to \(|\Psi|^2\):

\[\begin{split}\begin{aligned} A(\tau) &= \mathbb{E}_\mathbf{R}\left[\alpha(\mathbf{R})\right] \\ \alpha(\mathbf{R}) &= \mathbb{E}_\Delta\left[\min(1, e^X)\,\middle|\,\mathbf{R}\right] \end{aligned}\end{split}\]

Stage I: averaging over the proposal

Expanding to second order in the displacement,

\[X = 2\,\boldsymbol{\Delta}\cdot\mathbf{F} + \boldsymbol{\Delta}^{\mathsf T}\mathsf{H}\,\boldsymbol{\Delta} + O(\Delta^3)\]

The moments of the Gaussian proposal are \(\mathbb{E}_\Delta[\delta_k] = 0\), \(\mathbb{E}_\Delta[\delta_k \delta_l] = \tau^2\,\delta_{kl}\) and \(\mathbb{E}_\Delta[\delta_k \delta_l \delta_m] = 0\). Since \(\mathbf{F}\) and \(\mathsf{H}\) are fixed by \(\mathbf{R}\), they come out of the average as constants:

\[\mathbb{E}_\Delta\left[2\,\boldsymbol{\Delta}\cdot\mathbf{F}\,\middle|\,\mathbf{R}\right] = 2\sum_k F_k\,\mathbb{E}_\Delta[\delta_k] = 0\]

The linear term dies from the symmetry of the proposal, not from any property of the drift. The quadratic term survives as a trace,

\[\mathbb{E}_\Delta\left[\boldsymbol{\Delta}^{\mathsf T}\mathsf{H}\boldsymbol{\Delta} \,\middle|\,\mathbf{R}\right] = \sum_{kl} \mathsf{H}_{kl}\,\mathbb{E}_\Delta[\delta_k \delta_l] = \tau^2\,\mathrm{Tr}\,\mathsf{H} = \tau^2\,\nabla^2\Lambda\]

and in the variance the cross term vanishes because it involves only third moments of a symmetric distribution:

\[\begin{split}\begin{aligned} \mathrm{Var}_\Delta[X|\mathbf{R}] &= \mathrm{Var}_\Delta\!\left(2\boldsymbol{\Delta}\cdot\mathbf{F}\right) + \underbrace{\mathrm{Var}_\Delta\!\left(\boldsymbol{\Delta}^{\mathsf T}\mathsf{H} \boldsymbol{\Delta}\right)}_{O(\tau^4)} + \underbrace{2\,\mathrm{Cov}_\Delta}_{0} \\ &= 4\,\tau^2\,\mathbf{F}\cdot\mathbf{F} + O(\tau^4) \end{aligned}\end{split}\]

The mean is therefore of second order in the step size while the variance survives at first non-vanishing order, and both are kinetic energies:

\[\begin{split}\begin{aligned} m(\mathbf{R}) &\equiv \mathbb{E}_\Delta[X|\mathbf{R}] = -4\,\tau^2\,T_M(\mathbf{R}) \\ s^2(\mathbf{R}) &\equiv \mathrm{Var}_\Delta[X|\mathbf{R}] = 8\,\tau^2\,T_D(\mathbf{R}) \end{aligned}\end{split}\]

The asymmetry is essential: the mean sees the combination \(T_L + T_D\), the variance sees \(T_D\) alone, and the Laplacian form does not enter the variance at all.

At fixed \(\mathbf{R}\) the leading part of \(X\) is a linear function of a Gaussian vector and is therefore exactly normal, whatever the number of components that move. This is why the proposal is Gaussian: with a uniform one, normality would rest on a central limit theorem over the \(3N_e\) summands, which EBES, moving three at a time, does not have. With that,

\[\begin{split}\begin{aligned} \alpha(\mathbf{R}) &= g\big(m(\mathbf{R}), s(\mathbf{R})\big) \\ g(m, s) &= \Phi\!\left(\frac{m}{s}\right) + e^{m + s^2/2}\,\Phi\!\left(-s - \frac{m}{s}\right) \end{aligned}\end{split}\]

which is the standard integral

\[\mathbb{E}[\min(1,e^X)] = P(X \ge 0) + \int_{-\infty}^{0} e^x\,\mathcal{N}(x; m, s^2)\,dx\]

Note that pointwise

\[m(\mathbf{R}) + \tfrac{1}{2}s^2(\mathbf{R}) = 2\,\tau^2\left[T_D(\mathbf{R}) - T_L(\mathbf{R})\right] \neq 0\]

so \(g\) cannot be simplified at this stage. It becomes simplifiable only after stage II.

Stage II: averaging along the chain

Integrating by parts against \(|\Psi|^2\), with the boundary term vanishing for a bound state,

\[\mathbb{E}_\mathbf{R}[T_D] = \frac{1}{2}\int |\nabla\Psi|^2 = -\frac{1}{2}\int \Psi \nabla^2\Psi = \mathbb{E}_\mathbf{R}[T_L]\]

so the three estimators are unbiased for the same quantity and differ only in variance:

\[\mathbb{E}_\mathbf{R}[T_L] = \mathbb{E}_\mathbf{R}[T_D] = \mathbb{E}_\mathbf{R}[T_M] \equiv \langle T \rangle\]

By the laws of total expectation and total variance,

\[\begin{split}\begin{aligned} \mu \equiv \mathbb{E}[X] &= \mathbb{E}_\mathbf{R}[m] = -4\,\tau^2 \langle T \rangle \\ \sigma^2 \equiv \mathrm{Var}[X] &= \mathbb{E}_\mathbf{R}[s^2] + \mathrm{Var}_\mathbf{R}[m] = 8\,\tau^2 \langle T \rangle + O(\tau^4) \end{aligned}\end{split}\]

the discarded piece being \(16\,\tau^4\,\mathrm{Var}_\mathbf{R}(T_M)\). Hence \(\mu = -\sigma^2/2\), which is not an assumption but a consequence: stationarity of the chain requires \(\mathbb{E}[e^X] = 1\), that is \(\mu + \sigma^2/2 = 0\) for normal \(X\), and the expansion reproduces it by itself. It appears only after stage II — stage I does not give it.

The second relation is a sum rule:

\[\mathrm{Var}[X] = 8\,\tau^2 \langle T \rangle\]

and nothing else enters. Neither the geometry of the molecule, nor the distribution of nuclear charge, nor the Jastrow factor and backflow appear other than through the single scalar \(\langle T \rangle\). The reason is that the sum rule is a full trace,

\[\mathbb{E}_\mathbf{R}\left[\sum_i |\nabla_i \Lambda|^2\right] = \mathrm{Tr}\,\mathbb{E}_\mathbf{R}\left[\nabla_\alpha \Lambda \nabla_\beta \Lambda\right] = 2 \langle T \rangle\]

and an isotropic step sees the trace only. Every multipole of the density cancels identically, so the law is blind to shape: linear, planar and cage-like molecules of the same composition obey it alike. Anisotropy would enter only through a step that treats the Cartesian directions differently and so probes the traceless part of \(Q_{\alpha\beta} = \mathbb{E}_\mathbf{R}[\sum_i \nabla_\alpha \Lambda \nabla_\beta \Lambda]\); no such step is implemented.

Interchanging the two stages

The exact acceptance is \(A = \mathbb{E}_\mathbf{R}[g(m(\mathbf{R}), s(\mathbf{R}))]\). Replacing \(T_D(\mathbf{R})\) and \(T_M(\mathbf{R})\) by their common mean \(\langle T \rangle\) before applying \(g\) — that is, commuting \(\mathbb{E}_\mathbf{R}\) with a non-linear function — sends \(m \to -\sigma^2/2\) and \(s \to \sigma\), so the exponential becomes unity and both arguments collapse to \(-\sigma/2\):

\[A(\sigma) \approx 2\,\Phi(-\sigma/2)\]

a one-parameter curve: acceptance depends on the step size through \(\sigma\) only. Setting \(A = 1/2\) gives \(\sigma = 2s = 1.3489795\) and, with the sum rule,

\[\tau_{50}\,\sqrt{\langle T \rangle} = \frac{s}{\sqrt{2}} = \mathrm{erfinv}(1/2) = 0.4769362\]

For a general target \(A\) the constant is \(\mathrm{erfinv}(1 - A)\); at the Roberts–Gelman–Gilks optimal acceptance 0.234 [22] it would be 0.8416. There is nothing fundamental about 0.477: it is fixed by three conventions — the step quoted as a standard deviation rather than as a variance, a 50 % target, and atomic units. Nor is there any reason for it to equal one, the prefactor being dimensional, \([\tau] = \sqrt{1/\mathrm{energy}}\).

Only the second moment of the proposal enters the sum rule, so the step size is quoted as that moment. This is what dtvmc is, and both codes draw each component from a Gaussian of variance dtvmc, so \(\mathtt{dtvmc} = \tau^2\) and the two agree on

\[\mathtt{dtvmc} \cdot \langle T \rangle = \left[\mathrm{erfinv}(1/2)\right]^2 = 0.2274682\]

which is the form to compare codes in, and the one used below.

approximate_step_size

To use the result as an initial guess, \(\langle T \rangle\) is needed before any sampling has been done. The virial theorem gives \(\langle T \rangle = |E|\), for an ion no less than for a neutral atom since the net charge enters nowhere, and the Thomas–Fermi expansion with its Scott and Dirac corrections [24] [25] supplies the energy with no fitted constant at all:

\[\langle T \rangle \approx \sum_a \left[ 0.7687\,Z_a^{7/3} - 0.5\,Z_a^2 + 0.2699\,Z_a^{5/3} \right]\]

which reproduces Hartree–Fock energies to 1.8 % rms. The leading term alone is much worse than it looks — 4.2 % rms, with the effective coefficient \(|E| / \sum_a Z_a^{7/3}\) drifting from 0.54 to 0.64 across the periodic table — and that drift is exactly the Scott term.

Neutrality is assumed by the expansion rather than by the virial theorem: its leading term minimizes the Thomas–Fermi functional subject to \(\int\!\rho = Z\). The number of electrons deliberately does not appear — substituting it for \(Z_a\) would put \(\mathrm{Be}^{2+}\) 80 % low in \(\langle T \rangle\) and 123 % high in \(\tau\) — so what the formula returns for an ion is the neutral atom carrying the same nucleus: \(-1.0\,\%\) on \(\mathrm{Li}^+\) and \(+4.7\,\%\) on \(\mathrm{Be}^{2+}\), \(+0.5\,\%\) and \(-2.3\,\%\) in \(\tau\). The two errors partly cancel, the expansion running low on a neutral atom while a cation sits below its neutral parent.

The expansion does extend to ions, and the extension is worth stating because it is not the obvious one. Thomas–Fermi theory scales as \(E = -Z^{7/3}e(\lambda)\) in \(\lambda = N/Z\) [24], the Scott term survives untouched — it comes from \(r \sim 1/Z\), where no valence electron screens anything — and only the exchange part of the \(Z^{5/3}\) coefficient follows the density. But there is no term linear in the net charge, because the neutral Thomas–Fermi atom has zero chemical potential: \(\partial E/\partial N\) vanishes at \(N = Z\), hence \(e'(1) = 0\). Solving the equation with the ion boundary condition confirms it, \(e(1) - e(\lambda)\) falling as \((q/Z)^{2.45}\) at \(q/Z = 0.02\) and approaching the classical \((q/Z)^{7/3}\) from above. Carried out exactly the correction works: it moves \(\mathrm{Li}^+\) to \(-4.3\,\%\) and \(\mathrm{Be}^{2+}\) to \(-4.1\,\%\), which is where the same expansion already sits on neutral lithium and beryllium, \(-3.6\,\%\) and \(-2.3\,\%\). The ion-specific error is genuinely removed and what is left is the series being 2 to 4 % low at \(Z\) of 3 and 4. That is also why it is not used here: the residual it exposes is the one that was cancelling against it, so nothing is gained in \(\tau\), while the price is a numerical solution of the Thomas–Fermi equation and a per-atom partition of the electron count that a molecular ion does not admit.

Writing \(T_a\) for the share of \(\langle T \rangle\) carried by nucleus \(a\), the remaining correction is expressed through its participation ratio,

\[\begin{split}\begin{aligned} n_\mathrm{nuc} &= \frac{\left(\sum_a T_a\right)^2}{\sum_a T_a^2} \\ \tau &= \mathrm{erfinv}(1/2)\, \frac{1 + 0.080 / n_\mathrm{nuc}^{0.82}}{\sqrt{\langle T \rangle}} \end{aligned}\end{split}\]

which is \(1\) for a single heavy atom and the number of equivalent nuclei for a symmetric molecule. Everything here follows from the sum rule and from tabulated atomic physics except the \(0.080\) and the exponent, the only fitted numbers in the estimate; the next section explains what they stand for. Measured on the systems in examples/time_step/CBCS the 50 % acceptance point lands within 0.3 % rms of this guess, and 0.8 % at worst, for atoms, ions, hydrides and hydrocarbons alike.

The other two sampling modes take the same estimate through their own sum rule, since only the block of coordinates that moves changes:

\[\tau(\mathrm{EBES}) = \sqrt{N_e}\;\tau, \qquad \tau(\mathrm{DBDS}) = \sqrt{2}\;\tau\]

Both are exact statements about \(\mathrm{Var}[X]\), and in dtvmc they read \(N_e\) and \(2\). What the EBES factor is knowingly missing is not a failure of the Gaussian law but the spread of the drift between electrons, which the sum rule averages over; that is left to optimize_vmc_step, whose fixed point does not depend on the model.

Cost of the interchange

Suppose the trial function were exact. The zero-variance property fixes \(E_L\) at the eigenvalue and with it \(T_L = E - V\) pointwise, so the better the wave function, the more rigidly the Laplacian estimator is tied to the potential — which diverges at every nucleus. Perfection does not remove the fluctuation the interchange discards; it prescribes it.

What the interchange needs is not an accurate wave function but a featureless one, \(\mathrm{Var}(T_D) = 0\), a log-density of constant slope, and it is again the exact wave function that forbids it. The Kato cusp fixes \(|\nabla_i\Lambda| \to Z_a\) as electron \(i\) reaches nucleus \(a\), the asymptotic decay \(\Psi \sim e^{-\sqrt{2I}\,r}\), with \(I\) the ionization energy the departing electron leaves against, fixes it at \(\sqrt{2I}\) far out, so each electron’s share of \(T_D\) swings between \(Z_a^2/2\) and \(I\) — on neon, 50 against 0.79 — whatever is done to the ansatz. The correction below is therefore a property of the system and not a defect of the trial function, which is why nothing about the latter enters it. Measured on the campaign, the coefficient of variation \(\mathrm{CV}(T_D) = \sqrt{\mathrm{Var}_\mathbf{R}(T_D)}\,/\,\langle T \rangle\) grows with the nuclear charge, 0.10 on helium, 2.5 on neon, 3.1 on krypton, and falls when the same charge is spread over more nuclei, 1.3 on \(\mathrm{C_2H_2}\) against 0.5 on \(\mathrm{C_6H_6}\).

The same asymmetry is reported from the diffusion side by Umrigar, Nightingale and Runge [26], who observe that the acceptance of a drift–diffusion move is considerably lower for electrons close to a nucleus than for electrons far from any, and who raise it by replacing the Gaussian proposal near a nucleus with one carrying the cusp. Their remedy does not transfer: it repairs a Green function that is an approximation, whereas the Metropolis walk samples \(|\Psi|^2\) exactly at every step size, so here the same nonuniformity costs efficiency and not accuracy. The observation does transfer, and it is the one being quantified above.

The interchange is thus the only uncontrolled step of the derivation. Its sign follows without any calculation: \(g\) is convex in \(T_D\), so by Jensen’s inequality the spread of \(T_D\) raises the acceptance and a larger step is needed. The correction is positive by construction and is exactly the Jensen gap of commuting \(\mathbb{E}_\mathbf{R}\) with \(g\).

In terms of the dimensionless fluctuations \(u = T_D/\langle T\rangle - 1\) and \(v = T_L/\langle T\rangle - 1\), expanding the exact acceptance to second order and solving \(A = 1/2\) gives

\[\frac{\delta\,\tau}{\tau} = c_{DD}\,\mathbb{E}_\mathbf{R}[u^2] + c_{DL}\,\mathbb{E}_\mathbf{R}[uv] + c_{LL}\,\mathbb{E}_\mathbf{R}[v^2]\]

\[\begin{split}\begin{aligned} c_{DD} &= \frac{1}{8} - \frac{s^2}{4} + \frac{s^3}{16\varphi(s)} = +0.071617 \\ c_{DL} &= \frac{s^2}{2} - \frac{s^3}{8\varphi(s)} = +0.106766 \\ c_{LL} &= -\frac{s^2}{8} + \frac{s^3}{16\varphi(s)} = +0.003484 \end{aligned}\end{split}\]

with \(c_{DD} + c_{DL} + c_{LL} = (1 + s^2)/8\), \(c_{DD} - c_{LL} = (1 - s^2)/8\) and \(c_{DL} + 2c_{LL} = s^2/4\); the terms in \(s^3/\varphi\) cancel in the sum. Treating \(\nabla^2\Lambda\) as non-fluctuating, that is setting \(T_L(\mathbf{R}) = T_D(\mathbf{R})\) pointwise, collapses this to a single term,

\[\frac{\delta\,\tau}{\tau} = \frac{1 + s^2}{8}\,\frac{\mathrm{Var}_\mathbf{R}(T_D)}{\langle T \rangle^2} = 0.18187\,\mathrm{CV}^2(T_D)\]

and this is the form the fitted \(1 + 0.080/n_\mathrm{nuc}^{0.82}\) was read off: \(T_D\) is dominated by whichever electron is currently near a nucleus, since \(|\nabla\Lambda| \sim Z\) there, so \(n\) equivalent nuclei contribute \(n\) independent terms and \(\mathrm{CV}^2(T_D)\) falls as \(1/n\). The exponent comes out at \(0.82\) rather than the \(1\) that argument predicts.

The campaign says the identification is not that direct. Of the two factors the measured correction separates into, the Jensen gap is the second one, the shape of \(X\), and there the collapsed formula is exact where \(\mathrm{CV}(T_D)\) is small and out of range where it is not: helium, at \(\mathrm{CV} = 0.10\), is predicted at \(+0.2\,\%\) and measured at \(0.0\,\%\), while krypton, at \(3.1\), is predicted at \(+175\,\%\) and measured at \(+16\,\%\). The first factor, the \(O(\tau^2)\) deficit of the sum rule itself, is not a Jensen gap at all — it survives at \(\mathrm{Var}(T_D) = 0\), being a property of how \(\nabla\Lambda\) turns under a finite displacement rather than of how it fluctuates — and on the light atoms it carries the whole correction. The constant is as flat in \(Z\) as it is because those two cancel, so the mechanism above should be read as accounting for one of them and not for the fit.

This one-constant form should not be trusted quantitatively. Rewritten in the variables Casino prints, \(t = T_M/\langle T\rangle - 1 = (u + v)/2\), the same expression reads

\[\frac{\delta\,\tau}{\tau} = -0.03167\,\mathbb{E}_\mathbf{R}[u^2] + 0.19959\,\mathbb{E}_\mathbf{R}[ut] + 0.01395\,\mathbb{E}_\mathbf{R}[t^2]\]

The coefficient of \(\mathrm{Var}(T_D)\) has turned negative and the covariance dominates, because \(T_M = -\nabla^2\Lambda/4\) diverges as \(Z/2r\) at a nucleus — from \(\Lambda \approx -Zr\) follows \(\nabla^2\Lambda = -2Z/r\) — whereas \(T_D\) stays bounded there, \(|\nabla\Lambda| \to Z\). The variances of the estimators are not comparable: on helium the two error bars differ by an order of magnitude, so \(\mathrm{Var}(T_M)/\mathrm{Var}(T_D)\) is of order 70 and \(\mathrm{CV}(T_L)\) is of order unity, where a second-order expansion in \(v\) no longer applies. A single fitted constant absorbs \(\mathbb{E}_\mathbf{R}[uv]\), \(\mathbb{E}_\mathbf{R}[v^2]\) and the truncation of the series alike.

The interchange is avoidable in principle

The exact acceptance is the average of a known function of two quantities that are already evaluated at every configuration, so accumulating \(T_M\) and \(T_D\) over the equilibration walk — which runs at \(|\Psi|^2\) before the step is optimized — and solving \(A(\tau) = 1/2\) numerically would remove the Thomas–Fermi estimate, the virial theorem and the fitted constant together, and would account for the Jastrow factor, backflow, ions and pseudopotentials for free, since the function actually being sampled is the one being measured.

Range of validity

Of stage I:

1. the expansion of \(X\) is truncated at second order in \(\boldsymbol{\Delta}\), and at the 50 % point \(\tau^2 \langle T\rangle \approx 0.23\), so the step is not small. Measured, the sum_rule column at that point runs from 0.85 on helium to 1.15 on krypton, so what is truncated is a 15 % effect in either direction;

2. \(\Psi\) is not analytic at a nuclear cusp, which is precisely where \(T_L\) diverges, so the Taylor expansion formally fails there;

3. the leading term is exactly normal, the proposal being Gaussian, so nothing here rests on a central limit theorem. What is left non-normal at fixed \(\mathbf{R}\) is the quadratic term of point 1, and it enters at order \(\tau^4\).

Of stage II:

4. the \(\mathrm{Var}_\mathbf{R}[m]\) contribution to \(\mathrm{Var}[X]\) is dropped at order \(\tau^4\);

5. the interchange is controlled only to second order and requires small \(\mathrm{CV}(u)\) and \(\mathrm{CV}(v)\); the second is not satisfied in practice.

Outside both stages:

6. the virial theorem holds for an eigenstate, or for a trial function optimized with respect to a uniform scaling of all coordinates; the sum rule itself holds for any smooth \(\Psi\).

Pseudopotentials and the electron gas

Neither stage of the derivation ever differentiates the potential: everything is built from \(\nabla\Lambda\) and \(\nabla\nabla\Lambda\), and \(V\) enters only the local energy, which the walk never sees. The sum rule \(\mathrm{Var}[X] = 8\,\tau^2\langle T \rangle\) therefore holds verbatim for a pseudopotential, local or non-local, and for a periodic system. What breaks is only point 6 above, the route from the Hamiltonian to a number for \(\langle T \rangle\):

• the virial theorem needs a potential homogeneous of degree \(-1\) under a uniform scaling. A pseudopotential is not, so \(\langle T \rangle \neq |E|\); in the Casino runs under examples/ppotential_HF the ratio to the valence energy runs from 0.62 on carbon to 0.81 on neon, with no plateau to extrapolate to;

• Thomas–Fermi is a semiclassical theory of many electrons and says nothing about the handful of valence electrons a pseudopotential leaves behind. Fed the pseudo charge it overshoots fourfold on carbon — \(13.2\) against a measured \(3.3\) — because it places those four electrons in a \(1s\) shell.

The valence electrons are hydrogenic instead. With \(Z_a\) the pseudo charge of nucleus \(a\), \(n_a\) the principal quantum number of its row, and Slater screening by the other valence electrons,

\[\langle T \rangle \approx \sum_a \frac{Z_a (0.65\,Z_a + 0.35)^2}{2\,n_a^2}\]

This is exact on hydrogen and 11 % to 55 % high across the second row, a pseudo-orbital being smoother than the hydrogenic one it is modelled on, which halves into a 5 % to 20 % low estimate of \(\tau\) — well inside the basin of optimize_vmc_step.

That the sum rule itself survives can be read off existing Casino output, in the proposal-independent form. Taking every CBCS run under examples and multiplying the optimized dtvmc by the KEI the same run reports:

Runs	Number	\(\langle T \rangle\) (au)	\(\mathtt{dtvmc}\cdot\langle T \rangle\)
all-electron	184	2.4 … 2781	0.261 ± 15 %
pseudopotential	72	0.49 … 57	0.293 ± 14 %

The product is constant over a factor of 1160 in \(\langle T \rangle\), with and without a pseudopotential — which is the claim being tested, and the only one these runs settle. Both medians sit above the \(0.2274682\) of the law by the Jensen gap, \(+7\,\%\) and \(+14\,\%\) in \(\tau\) — the all-electron figure compatible with the \(8.0\,\%\) fitted above, the pseudopotential one larger as expected, since removing the core leaves only three to eight electrons to average \(T_D\) over. The campaign lands in the same place from its own side, \(\mathtt{dtvmc}\cdot\langle T\rangle\) measuring 0.266 over the atoms and ions and 0.258 over all seventeen all-electron systems. The \(15\,\%\) scatter, however, is not physics but the tolerance of Casino’s own step-size optimizer, which stops as soon as the acceptance is near enough to 50 %, so nothing finer — the \(n_\mathrm{nuc}\) collapse in particular — can be extracted from these runs. That is what the campaign described in optimize_vmc_step is for. There is also no pseudopotential molecule to test it on: \(\mathrm{B_2H_6}\) under examples/ppotential_HF carries a pseudopotential on hydrogen only, and pseudo-hydrogen is hydrogen.

The homogeneous electron gas is the opposite limit and the cleanest test the law admits. \(\langle T \rangle\) is known in closed form, without a virial theorem or a fit,

\[\langle T \rangle = \frac{3}{10} \left(\frac{6\pi^2}{\Omega}\right)^{2/3} \left(N_\uparrow^{5/3} + N_\downarrow^{5/3}\right) = \frac{1.10495\,N_e}{r_s^2} \quad (\text{unpolarized})\]

so the prediction carries no adjustable constant at all: \(\mathtt{dtvmc} = 0.2058\,r_s^2 / N_e\) in CBCS, and \(N_e\) times that, \(0.2058\,r_s^2\), in EBES, where the kinetic energy enters per electron and the size of the system drops out.

One condition attaches to that closed form. What the sum rule needs is \(\langle T_D \rangle\), and \(\Lambda = \ln|\Psi|\) knows the modulus only, so for a complex wave function the physical kinetic energy exceeds it by the current term \(\langle |\nabla\theta|^2 \rangle / 2\), which \(|\Psi|^2\) cannot see. At the \(\Gamma\) point with a closed shell the determinant is real up to a global phase and the two coincide, which is the case the number above describes. Under twisted boundary conditions — the standard route to finite-size corrections — or with an open shell, the determinant carries a current, \(3 k_F^2 N_e / 10\) is then larger than what the sum rule is asking for, and the step comes out too small. Nothing in the derivation breaks: it is the closed-form input that stops being the right quantity, and the free-particle expression should be read as the modulus-only kinetic energy wherever the two part company.

Translational invariance costs nothing here, which is worth saying because it looks as though it should. Displacing every electron by the same vector leaves \(|\Psi|^2\) exactly unchanged, so of the \(3N_e\) components a CBCS move proposes, three do not enter \(X\) at all, while EBES moves one electron and cannot produce such a displacement in the first place. Neither formula changes, because neither counts degrees of freedom: both follow from \(\sum_k F_k^2\), and the invariance is precisely the statement \(\sum_i \nabla_i \Lambda = 0\), which puts no weight in that subspace to begin with. The bookkeeping is automatic. What the free subspace does cost is efficiency, at order \(1/N_e\): the centre of mass random-walks with the full step and its displacement enters the diffusion criterion below without ever having affected the acceptance. No atom or molecule has this subspace, the clamped nuclei supplying the missing force.

The Jensen gap should be far smaller than in any atom, since there is no nuclear cusp for \(T_D\) to be dominated by and \(\mathrm{CV}(T_D)\) falls as \(1/\sqrt{N_e}\), so the interchange that costs 4.5 % on a single atom should cost almost nothing here. Pycasino has no homogeneous-gas wave function, so this is a prediction rather than a measurement; atom_kinetic_energy returns an empty array for a system without nuclei and the guess falls back to unity, which optimize_vmc_step recovers from without help.

The gas is also the only system for which a published estimate of the step size exists to check against. Section II.C of [23] sets the RMS distance an electron diffuses in one EBES move equal to \(r_s\), assumes the acceptance is 50 %, and obtains \(\mathtt{dtvmc} = 2 r_s^2/3 = 0.667\,r_s^2\) — a factor 3.24 above the \(0.2058\,r_s^2\) the sum rule gives at that same 50 %. The estimate does not produce the acceptance it assumes, which is one reason its authors call it far from optimal. But the 50 % target is a convention too, and the alternative Casino implements as OPT_DTVMC 2 is to maximize the diffusion constant, the mean square displacement per attempted move, which is proportional to \(\mathtt{dtvmc}\cdot A\). Writing \(u = \sigma/2\), that is the maximum of \(u^2\,\Phi(-u)\), stationary where

\[2\,\Phi(-u) = u\,\varphi(u), \qquad u = 1.1906, \qquad A = 0.2338\]

so the diffusion-optimal dtvmc is \((u/s)^2 = 3.12\) times the one the 50 % rule sets, and the acceptance it wants is the Roberts–Gelman–Gilks 0.234 [22] — obtained here from the sum rule and a definition of efficiency rather than from a scaling limit. Carried through, \(\mathtt{dtvmc} = 0.641\,r_s^2\), within 4 % of \(2 r_s^2/3\). The published estimate is therefore not the 50 % step it is derived as but, to that accuracy, the diffusion-optimal one: taking \(r_s\) as the length overshoots the 50 % step by 3.24, the target the criterion should have used overshoots it by 3.12, and the two cancel. Both figures inherit the reservation of the next section: \(A(\sigma)\) describes a single electron and the gas has only one kind of electron, so here it should hold — it is in an atom, where a core and a valence electron are the same move at wildly different drift, that it fails.

EBES

Only stage I changes. With \(\mathbf{R}\) fixed and the electron \(i\) chosen, the sum runs over 3 components rather than \(3N_e\)

\[\mathrm{Var}_\Delta[X | \mathbf{R}, i] = 4\,\tau^2\,|\nabla_i \Lambda|^2\]

Averaging over the uniform choice of electron as well gives \(\mathbb{E}_i \mathbb{E}_\mathbf{R}[|\nabla_i\Lambda|^2] = 2\langle T\rangle / N_e\), so

\[\mathrm{Var}[X] = 8\,\tau^2\,\frac{\langle T \rangle}{N_e}, \qquad \tau_{50}(\mathrm{EBES}) = \sqrt{N_e}\;\tau_{50}(\mathrm{CBCS})\]

The EBES step is therefore set by the kinetic energy per electron, an intensive quantity: at fixed composition it does not depend on the size of the system, and for neutral atoms it scales as \(Z^{-2/3}\).

This mode is why the proposal is Gaussian. Only three components move, so a uniform proposal would leave \(X\) a weighted Irwin–Hall variable with no limit for a central limit theorem to be taken in, and the drift of a single electron near a nucleus is nearly radial, which drives the participation ratio \((\sum_k F_k^2)^2 / \sum_k F_k^4\) from 3 towards 1 — one scaled uniform variable, bounded support, excess kurtosis \(-6/5\). How much that mattered was estimated by integrating the cube against \(\int\min(1,e^x)\), which gave a step \(3.0\,\%\) from the Gaussian answer at a participation ratio of 1 and under \(0.5\,\%\) at 2 or 3 at the 50 % target, two symmetric distributions of equal variance agreeing in the body.

Both proposals have since been run over the same thirty-one systems, and the estimate was an order of magnitude low. Taking the width of \(X\) at its measured variance — the second factor of the decomposition below, which is exactly the quantity that integral was after — the Gaussian gives \(1.000\) on every system with a single occupied shell, \(0.999\) to \(1.011\) in CBCS and \(0.9995\) to \(1.013\) in EBES, while the cube gives \(0.956\) in CBCS and \(0.918\) in EBES. The deficit doubles when \(3N_e\) components are replaced by three, which is the missing central limit theorem measured directly, and it moves the 50 % step by \(2.3\,\%\) and \(8.3\,\%\) respectively — not by the \(0.5\,\%\) argued above. The shape itself is where it should be: the cube’s excess kurtosis at the smallest step is \(-0.719\) on pseudo-hydrogen, against the \(-1.2\,\mathbb{E}[\sum_k n_k^4] = -0.72\) of three components and an isotropic drift, and \(-0.337\) on helium where six components share. What the integral got wrong is how much a small deviation from normality is worth at the 50 % point. And the deviation is not even a constant: it runs from \(-0.72\) to zero along the acceptance grid, so no fitted correction could have absorbed it.

A Gaussian displacement removes the whole question, \(X\) being exactly normal at leading order however few components move, and the bounded support that a cube offers buys nothing here.

What the measurement says

The acceptance curve of all thirty-one calibration systems has been measured in this mode, on the same grid and the same conventions as the CBCS set, at \(10^6\) steps per point. The cube campaign it replaces is kept beside it in examples/time_step/EBES.cube, at \(10^5\).

The sum rule is confirmed by its sharpest consequence: at fixed composition the step does not depend on the size of the system. C₂H₂, C₄H₄ and C₆H₆, from 14 to 42 electrons, take dtvmc of 0.17307, 0.17295 and 0.17115, equal to \(1.1\,\%\).

The Gaussian acceptance law on top of it does not hold, and splitting correction into its two factors says which half fails where.

Where every electron is equivalent — pseudo-hydrogen, helium, the two-electron ions — the second factor is 1.000 to 1.013. \(X\) is Gaussian, with nothing left of the \(-0.33\) the cube used to put there. The step is nevertheless 1.13 to 1.15 times the \(\sqrt{N_e}\) law, and all of that is the first factor, the \(O(\tau^2)\) deficit of the sum rule at a step whose per-electron variance is twice the CBCS one. The prediction that the Gaussian alone would take the one-shell excess to unity was half right: the shape part went to 1.000, the finite-step part was never the proposal’s to remove.

Everywhere else the second factor is 1.94 to 2.29 — beryllium 1.94, nitrogen 2.19, neon 2.04, argon 2.19, krypton 2.18, CH₄ 2.25, O₃ 2.16, C₆H₆ 2.22 — while the pseudoatoms, carrying a valence shell and no core, sit at 1.21 to 1.25. What it responds to is the presence of a core, not the number of shells and not the nuclear charge. The trend the raw correction shows — about 2.1 across the second period against 2.56 on argon and 2.82 on krypton — is the other factor, which runs from 0.92 on beryllium to 1.30 on krypton.

The reason is that the acceptance is a nonlinear functional of \(\mathrm{Var}[X|\mathbf{R},i] \propto |\nabla_i\Lambda|^2\), whereas the sum rule constrains only its mean. At a step where the average electron would sit at 50 %, a core electron is rejected almost always and a valence one accepted almost always, so the 50 % is made up by the electrons whose drift is below the mean, and the step that achieves it is correspondingly larger. In those terms the second factor squared, 3.8 on beryllium and between 4.2 and 5.2 everywhere else, is the ratio of \(\mathbb{E}_i[|\nabla_i\Lambda|^2]\) to the value carried by the electron sitting at the median of the acceptance; that it varies so little over a factor of 21 in \(\langle T\rangle/N_e\) is not explained by anything here and deserves a look of its own. The same mixture is what the kurtosis column reports: 2.0 to 5.3 on the all-electron systems with a core, 4.6 to 6.1 on the pseudoatoms, and 0.00 on the one-shell systems, where the Gaussian proposal leaves nothing at all.

The step this implies is the root of \(\mathbb{E}_i[2\Phi(-\sigma_i/2)] = 1/2\) over the spread of \(|\nabla_i\Lambda|^2\) across electrons rather than over its mean, which requires no fitted constant, only the per-electron drift. That replacement has not been made: approximate_step_size still returns the \(\sqrt{N_e}\) of the leading order, and optimize_vmc_step carries the rest in its three iterations.

DBDS

vmc_method : 2 selects determinant-by-determinant sampling, in which one spin determinant is displaced at a time — all \(N_\uparrow\) up-spin electrons together, then all \(N_\downarrow\) down-spin ones, each with its own accept/reject step. It sits between the two implemented modes: EBES with a block size of one electron, CBCS with a block size of \(N_e\).

It is not implemented. casino.vmc.VMC.random_step() dispatches methods 1 and 3 only and returns False for anything else, so the walker would never move and the run would silently produce the initial configuration repeated vmc_nstep times. The input check therefore refuses the keyword outright; what is left of the mode is the branch of approximate_step_size returning the \(\sqrt{2}\) of the sum rule below, and the sum rule itself.

Casino removed its own Method 2 as well: it “did not offer any advantage over the other methods and was hard to support”, and the one feature worth keeping — using acceptance probabilities in the accumulation — was moved into Method 3.

The sum rule says what could be expected of it. A move of the \(\sigma\) determinant is the derivation of stage I with the sum restricted to \(3N_\sigma\) components, so

\[\mathrm{Var}[X] = 8\,\tau^2 \langle T_\sigma \rangle, \qquad \langle T_\sigma \rangle = \frac{1}{2}\,\mathbb{E}_\mathbf{R} \left[\sum_{i \in \sigma} |\nabla_i\Lambda|^2\right]\]

For a closed-shell system the two spin channels carry half the kinetic energy each, so \(\tau_{50}(\mathrm{DBDS}) = \sqrt{2}\,\tau_{50}(\mathrm{CBCS})\) — a factor \(\sqrt{2}\), against \(\sqrt{N_e}\) for EBES. Three limitations follow.

• The gain in step size is bounded by \(\sqrt{2}\) whatever the system, because there are only ever two determinants. It does not grow with \(N_e\), unlike the EBES gain.

• An open-shell system cannot be served by one dtvmc. The optimal step of each channel is \(0.477/\sqrt{\langle T_\sigma \rangle}\), and the two differ whenever \(N_\uparrow \neq N_\downarrow\); a single step size puts one determinant off its 50 % target. Neither the input keyword nor optimize_vmc_step has room for two values.

• The saving DBDS aims at is already taken by EBES. Its point is that moving one determinant requires updating only that determinant’s inverse, not both. EBES keeps the inverse of every slater matrix in the walker and updates it by a rank-one formula when a single electron moves, which is the same machinery at a block size of one, so DBDS would add a coarser version of what is there.

The one theoretical argument that used to be in its favour has expired. With a uniform proposal DBDS kept the central limit theorem that EBES lost, \(X\) being a sum of \(3N_\sigma\) terms rather than 3, so the Gaussian form of \(g\) remained usable where EBES had to do without it. A Gaussian proposal gives that to EBES for nothing, and what is left is a block size of \(N_\sigma\) buying \(\sqrt{2}\).

acceptance_ratio

Measures the fraction of accepted moves at the current step size, counted over the proposals themselves: one electron in EBES, the whole configuration in CBCS. That is the quantity the 50 % rule is about in either mode, and both are then read off the same counter.

Reading it off the walk instead — the fraction of steps in which the configuration changed at all — works in CBCS and saturates in EBES, where a step is a sweep over the electrons. The fraction then measured is \(1 - (1 - a)^{N_e}\) for a per-electron probability \(a\), and inverting that expression is exact but useless: at the 50 % target Kr changes at least one of its 36 electrons with probability \(1 - 2^{-36}\), so every sample of it reads exactly one whatever the step size, the inverse returns one, and optimize_vmc_step is left with no scale to work from — it multiplies the step by the ten its clipped acceptance asks for, overshoots, and the three iterations end wherever the last overshoot left them.

optimize_vmc_step

Since \(\sigma \propto \tau\), one measurement already fixes the whole curve. Inverting \(A = 2\Phi(-\sigma/2)\) at the measured acceptance and reading it back at \(A = 1/2\) gives

\[\tau \;\longleftarrow\; \tau\, \frac{\mathrm{erfinv}(1/2)}{\mathrm{erfinv}\!\left(1 - A(\tau)\right)}\]

which lands on the target in a single shot whenever the Gaussian law is exact. Two properties make this safe to iterate rather than fit:

• its fixed point is \(A = 1/2\) for any monotone acceptance curve, since the two inverse error functions cancel identically there. The law sets the rate of convergence, not the answer, so the map cannot converge to the wrong step size where the Gaussian fails — EBES, a one-electron system, a badly optimized wave function;

• it is a genuine Newton step in the variable the law is linear in, so it converges from far away. Replayed on the measured curves in examples/time_step/CBCS, a 40 % error becomes 2 %, then 0.2 %, then 0.02 %.

Three iterations are used. They cost a third of the eleven-point scan and empirical two-parameter fit they replace, they start from the current step size rather than resetting to the guess — so a re-optimization after a wave-function update is a single correction — and the residual is set by the noise of the acceptance, \(\delta\tau/\tau = \delta A / 2\varphi(s)\,s = 2.33\,\delta A\), rather than by the model.

casino.pycasino.Casino.vmc_step_graph() measures the shape of the curve for the same reason it is not fitted here. Its grid is laid out equally in acceptance — nineteen targets from 0.95 to 0.05 — through the same inverse, \(\tau = \tau_{50}\,\mathrm{erfinv}(1 - A)/\mathrm{erfinv}(1/2)\), anchored on the \(\tau_{50}\) the optimizer has just measured. A grid in step size instead spends most of its points in the saturated tails, where the acceptance carries no information about \(\tau_{50}\), and spends them at system-dependent places, so nothing is comparable between files; a grid in acceptance gives every system — atom, molecule, pseudoatom, gas — the same points in the only variable the law knows about. The step sizes are written out in atomic units, so the data stay unnormalized by the law they are used to test.

What the columns test

The acceptance on its own cannot say which half of the law has failed, since it is produced by the sum rule and the Gaussian shape acting together. The moments of \(X\) can, because the sum rule is a statement about the second moment alone and says nothing about shape, while detailed balance constrains \(X\) whatever its shape is. Each grid point is therefore recorded as six numbers taken from one and the same walk, arranged so that no two of them test the same statement.

acc_ratio

the acceptance, averaged as \(\langle \min(1, e^X) \rangle\) over the proposals rather than counted from the accept/reject decisions that follow them. The same mean with a smaller variance, the coin flip having been integrated out analytically.

correction

the step size actually used, divided by the one the law would prescribe for the acceptance just measured,

\[\text{correction} = \tau \left/ \left[ \frac{\mathrm{erfinv}\!\left(1 - A\right)}{\sqrt{\langle T\rangle}}\,\right]\right.\]

One means the law is exact at that point. This is the factor approximate_step_size carries as \(1 + 0.080/n_\mathrm{nuc}^{0.82}\), expressed in the same units, so the file offers the number to be fitted rather than an energy to be converted into one; the value in the 50 % row is that constant measured on that system. In EBES the denominator carries \(\sqrt{N_e}\) besides, the sum rule of a one-electron move seeing the kinetic energy per electron.

sum_rule

\(\mathrm{Var}(X) \big/ 8\,\tau^2\langle T\rangle\), the sum rule by itself. It assumes nothing about the shape of the distribution and never refers to the acceptance. One in the limit \(\tau \to 0\), where the rule is exact; what it loses at finite \(\tau\) is the \(O(\tau^2)\) term that the leading order drops.

gaussian

\(\langle X \rangle + \mathrm{Var}(X)/2\). A stationary walk with a symmetric proposal has \(\langle e^X\rangle = 1\) exactly, and the cumulant expansion of that identity, \(\kappa_1 + \kappa_2/2 + \kappa_3/6 + \kappa_4/24 + \dots = 0\), makes this column equal to \(-(\kappa_3/6 + \kappa_4/24 + \dots)\) — every cumulant above the second and nothing besides. It vanishes if and only if \(X\) is Gaussian, and it does so as \(\tau^3\), being led by the third cumulant of a quantity linear in the displacement. Divided by \(\mathrm{Var}(X)^{3/2}\) it becomes a skewness, which is the form in which systems may be compared; undivided it grows along the grid for dimensional reasons alone.

exp_mean

\(\langle e^X \rangle\), which is one for any correct stationary walk with a symmetric proposal. It tests the implementation rather than the physics, and costs nothing to report. Rare large \(X\) dominate the estimator, so it is informative only where the acceptance is not tiny.

kurtosis

the excess kurtosis of \(X\). Being dimensionless already, it has a finite \(\tau \to 0\) limit, and that limit is the shape of \(X\) at leading order. One mechanism sets it: the variance of \(X\) conditioned on the position is proportional to the local \(T_D\), so \(X\) is a scale mixture over the walk and is leptokurtic. The column therefore has a definite sign, and what it measures is the spread of the local kinetic energy — near zero for an atom with one occupied shell, large where a core and a valence shell differ by orders of magnitude in \(|\nabla\ln\Psi|^2\). A negative value is now a defect rather than physics; the cube files kept beside the campaign carry the \(-6/5\) per component underneath, which competed with the mixture and won wherever the mixture was weak — \(-0.72\) on pseudo-hydrogen in EBES, three components and an isotropic drift, against \(0.06\) for the same system with a Gaussian. The variance of this estimator is set by the eighth moment, so the small step end of the column is a sign and not a number.

The middle two columns multiply back into the first:

\[\text{correction} = \frac{1}{\sqrt{\text{sum\_rule}}}\;\cdot\; \frac{\sqrt{\mathrm{Var}(X)}}{2\sqrt{2}\,\mathrm{erfinv}\!\left(1 - A\right)}\]

which is an identity and not an approximation. The first factor uses \(\langle T\rangle\) and never the acceptance; the second uses the acceptance and never \(\langle T\rangle\). The failure of the law therefore separates cleanly into the part where the sum rule is breaking down and the part where \(X\) is not Gaussian, and the second factor — the true width of \(X\) against the width a Gaussian would need to reproduce the same acceptance — is above one exactly when \(X\) is leptokurtic.

The separation matters because the two factors are individually system-dependent and move in opposite directions. Across the all-electron atoms, at the 50 % row itself, both are monotone in \(Z\) and each spans some 17 %: the first falls from 1.087 on helium through 0.987 on neon to 0.931 on krypton, the second rises from 1.000 through 1.095 to 1.162, and their product — the number the initial guess is fitted to — stays inside 1.0 %, from 1.0754 to 1.0861. On the light atoms the sum rule is undershooting and \(X\) is exactly Gaussian; on the heavy ones the sum rule is overshooting and \(X\) is strongly leptokurtic; the fitted constant is what is left when the two cancel. At 5 % acceptance the cancellation is no longer complete, the factors spanning 1.056 to 1.323 and 1.020 to 1.176 with the product ranging over 8.6 %. The near-universality of the raw correction curve is therefore a cancellation between two effects, not a single law, and quoting it without the decomposition would hide that.

The level of the whole test is fixed by \(\langle T \rangle\), which the same run measures on the same walk and writes into the header, so no fitted constant enters anywhere.

Both estimators are measured, with error bars, and the two questions to ask of them come apart. On variance neither dominates: \(\mathrm{Var}(T_L)/\mathrm{Var}(T_D)\) measures 199 on helium, 0.1 on krypton, 0.2 on neon and 0.03 on pseudo-carbon. On whether the error bar can be believed, \(T_D\) wins everywhere. Between \(10^4\) and \(10^5\) samples the mean of \(T_L\) drifts from 2.590 to 2.732 on helium and from 2557 to 2796 on krypton, many times its own quoted bar, while \(T_D\) moves from 2.840 to 2.853 and from 2908 to 2754, each time within its. \(T_L\) has a heavy right tail — it is the unbounded one — and a sample standard deviation does not describe it.

So \(T_D\) is the reference, its larger bar on heavy atoms being an honest one, and \(T_L\) is reported beside it because the two share a mean by parts: their agreement within the quoted bars is the convergence check, and where it fails the file says so instead of quietly reporting one number.

vmc_corr_graph

Everything above fixes the step size and stops there. It says nothing about what the step is worth, which is the other half of

\[\mathcal{E} = \frac{1}{\mathrm{Var}(E_L)\; n_\mathrm{corr}\; T_\mathrm{iter}}\]

casino.pycasino.Casino.vmc_corr_graph() measures the missing factors on the grid vmc_step_graph already uses, so the two campaigns line up row for row. Each point is one unthinned walk, from which come the integrated autocorrelation times of the local energy, of a soft coordinate \(\sum_i r_i^2\) and of \(T_D\), the diffusion constant as a measured mean square displacement rather than as \(\mathtt{dtvmc}\cdot A\), the variance of \(E_L\), and the two wall times a Metropolis step is built from — one move and one local energy. CBCS only: in EBES a proposal moves one electron, and at the step that gives the average electron 50 % a core electron is accepted once in thousands of sweeps, so a correlation time measured there describes a coordinate that never moved rather than a slow mode.

Two prerequisites, both learned the hard way. The cusp correction must be on: an uncorrected Gaussian orbital leaves \(E_L \sim -Z/r\) at the nucleus, whose variance is finite but whose fourth moment is not, so the sample variance never settles — measured cuspless it swings by a factor of 400 between rows of one system, and every efficiency drawn from it is noise. And the timings must be taken as minima over the rows rather than averaged: a move costs the same whatever the step size, so the spread in those two columns is scheduler noise, which is one-sided.

The correlation time and the diffusion constant peak in different places

The central measurement is that \(n_\mathrm{corr}(E_L)\) has a minimum near 40 % to 55 % acceptance and rises again at larger steps, while \(D\) grows until roughly 20 % and then flattens. On helium, \(10^5\) steps per point:

acc   0.95  0.90  0.84  0.79  0.74  0.69  0.64  0.59  0.54  0.50  0.46  0.41  0.37  0.32  0.29  0.24  0.19  0.15  0.10
tau    268    74    32    19    16    10   8.7   8.2   6.4   8.1   8.6    18    14    17    21    14    39    39    69

and over the thirty-one systems measured, \(A\) at the minimum of \(n_\mathrm{corr}\) spans 0.23 to 0.63 while \(A\) at the maximum of \(D\) spans 0.14 to 0.27. Two mechanisms separate them, and neither involves a per-electron acceptance — CBCS proposes one displacement for the whole configuration and accepts or rejects it whole.

A rejected move repeats the sample. At acceptance \(A\) a state is held for \(1/A\) steps on average, so even if every accepted move landed somewhere statistically independent,

\[n_\mathrm{corr} \ge \frac{2}{A} - 1\]

This floor is pure loss for a bounded observable, and it is invisible to \(D\), which counts distance covered and has the zeros of the rejected moves already averaged into it: one rare long jump carries as much displacement as many short ones, and the local energy is not improved by the trade.

The repetition is not uniform over configuration space. The acceptance at \(\mathbf{R}\) is governed by \(T_D(\mathbf{R})\) through \(\mathrm{Var}(X|\mathbf{R}) = 8\,\mathtt{dtvmc}\,T_D(\mathbf{R})\), and where an electron is near a nucleus \(T_D\) is orders of magnitude above its mean, so the whole configuration is rejected almost every time. Those are exactly the configurations that carry the variance of \(E_L\), so the walk sticks in the tail of the very distribution being averaged, and the contrast between sticky and slippery regions sharpens exponentially with the step. Measured on helium the ratio of \(n_\mathrm{corr}\) to the floor above is 1.9 to 4.6 across the curve, never below one.

Which optimum is the efficiency depends on the decorrelation period

The two wall times differ by a fixed factor, \(T_\mathrm{energy}/T_\mathrm{move} = 4.7\) to 6.1 over the whole set, and the cost of one independent sample takes two different forms depending on whether moves may be thinned:

\[\text{no thinning:}\quad n_\mathrm{corr}\,\bigl(T_\mathrm{move} + A\,T_\mathrm{energy}\bigr) \qquad \text{optimal } p:\quad n_\mathrm{corr}\,T_\mathrm{move} + T_\mathrm{energy}\]

In the first, the energy is paid on every stored configuration, and the factor \(A\) appears because a rejected move leaves the configuration alone and needs no recomputation. A larger step is then attractive twice over — it lowers \(n_\mathrm{corr}\) down to the minimum, and it makes the expensive part happen less often — and the second reason carries the optimum past that minimum, down to where \(D\) is maximal. In the second, the energy is paid once per decorrelation, its share no longer depends on the step, and only \(n_\mathrm{corr}\) is left to minimize: the optimum sits at the minimum of \(n_\mathrm{corr}\) itself.

Measured, with the efficiency reconstructed from the columns at \(p = 1\) and at the \(p\) that minimizes the cost:

System	\(N_e\)	\(A\) at min \(n_\mathrm{corr}\)	\(A_\mathrm{opt}(p)\)	\(A\) at max \(D\)	\(A_\mathrm{opt}(1)\)	\(p\) at 50 %
CH₄	10	0.46	0.46	0.19	0.14	14
C₂H₂	14	0.35	0.45	0.22	0.16	16
Ar	18	0.40	0.32	0.18	0.18	18
O₃	24	0.35	0.35	0.25	0.25	23
C₄H₄	28	0.31	0.27	0.27	0.27	31
Kr	36	0.35	0.40	0.21	0.12	25

The left pair tracks and the right pair tracks, which is what makes the account above a mechanism rather than a story.

Which of the four combinations then gives the smallest error bar for a given amount of computer time? Taking the best achievable — thinning, at the step size that mode prefers — as one, the medians over the thirty-one systems are

	target 50 %	target max \(D\)
optimal \(p\)	0.84	0.66
no thinning	0.35	0.46

Thinning at the 50 % target wins outright, and never falls below 0.57 on any all-electron system in the set. Maximizing \(D\) while still thinning is worse everywhere it differs, because the thinning has already removed the only reason to prefer the larger step. Refusing to thin costs a factor of about two in efficiency whatever the target, which is a 1.17 to 1.57 times wider error bar for the same time — median 1.37, or 1.9 times the computer time for the same bar — and it is only within that already-lost regime that maximizing \(D\) beats the 50 % rule. So the recommendation to maximize \(D\) is sound for a calculation that stores every move, and is superseded by choosing \(p\).

The unit of that table is itself the maximum of a flat and noisy curve, so it is biased upwards and 0.84 is if anything an underestimate of what the default settings achieve.

How long the decorrelation loop should be

The \(p\) that minimizes the cost runs from 7 on the two-electron ions to 33 on benzene, inside the 8 … 36 that [23] reports for CBCS, and it is what casino.pycasino.Casino.optimize_decorr_period() estimates from the production block. Over all thirty-one systems it is described by

\[p_\mathrm{opt} = 1.35\,\sqrt{n_\mathrm{corr}\;T_\mathrm{energy}/T_\mathrm{move}}\]

to 18 %, where the obvious alternative \(p \propto n_\mathrm{corr}\) scatters by a factor of 2.3 across the same data, from 0.46 on benzene to 1.08 on helium. Both factors are needed and each explains one trend: the loop lengthens with system size through the correlation time, and it lengthens again wherever the energy is expensive relative to a move.

Against that, vmc_decorr_period defaults to 10 in Pycasino under CBCS, where Casino’s default is 3. Measured at the 50 % step on the thirty-one systems, a period of 3 costs a median 1.69 times the computer time of the optimum for the same error bar, up to 2.70 on fluorine with a pseudopotential, and a period of 1 costs 2.41; a flat 10 lands within 6 % of the optimum in the median and within 27 % everywhere. Carried to the production wave functions with the \(T_\mathrm{energy}/T_\mathrm{move}\) measured above, the same comparison against a period of 3 gives 1.44 for neon with a Jastrow, 2.02 with backflow, 2.54 for pseudo-nitrogen with a Jastrow and 3.20 with backflow. The default is worth more than any other single number in this chapter, and unlike the step size it costs nothing to get approximately right.

Under EBES the same measurement gives 4, and the default follows it back to Casino’s 3. Two things change at once and both shorten the loop. A sweep proposes a move for every electron in turn, so it moves something almost every time — the measured fraction is 1 to five figures on neon and on krypton — and the local energy is therefore paid on every stored configuration whatever the period, leaving the loop to save on \(p\,T_\mathrm{move}\) alone. And the correlation time it is working against is three to four times shorter to begin with, 3 to 8 sweeps at the optimum where CBCS needs 18 to 47 configurations. What remains is a shallow optimum at 4: a period of 3 costs 1 to 7 % over it, a period of 1 costs 1.5 to 1.8, and the 10 that is right for CBCS costs a quarter to a half.

That second factor is where the pseudopotentials part company with everything else. \(T_\mathrm{energy}/T_\mathrm{move}\) is 4.7 to 6.1 on every all-electron system, Gaussian or Slater, light or heavy — a remarkably stable number, since both times scale with \(N_e\) in the same way — but 11.4 to 18.1 on the pseudoatoms, where the non-local quadrature is evaluated on a grid of points per electron and the move is not. Neon is the clean comparison: all-electron it wants \(p = 16\) on ten electrons, pseudised it wants \(p = 23\) on eight. The consequence is practical: thinning matters most exactly where the local energy is most expensive, and refusing it costs 0.19 … 0.35 of the achievable efficiency on the pseudoatoms against 0.21 … 0.60 on the all-electron systems. Since pseudopotentials are the norm for the extended systems that [23] recommends maximizing \(D\) for, this is not a corner of the parameter space.

What the correlation times measure

For a diffusive walk the correlation time of an observable is a squared length divided by a diffusion constant,

\[n_\mathrm{corr}(f)\,D = L_f^2 = \frac{\mathrm{Var}(f)}{\langle|\nabla f|^2\rangle}\]

exact for an Ornstein–Uhlenbeck process and a variational bound in general. The product is the quantity to compare between systems, and it splits the same way the step size does: \(D\) is sampling and \(L_f\) is chemistry. For the local energy, over the all-electron atoms from helium to krypton in both bases at 50 % acceptance,

\[L_E \propto Z^{-0.81}, \qquad L_E\,Z = 1.36 \pm 0.28, \qquad n_\mathrm{corr} \propto Z^{0.71}\]

that is, the length scale of the core, as expected of an observable whose fluctuations are collected there, and the two exponents are consistent through \(D \propto \mathtt{dtvmc} \propto 1/\langle T\rangle \propto Z^{-7/3}\). The pseudoatoms sit an order of magnitude away, \(L_E\) of 0.36 to 1.3 bohr against 0.05 to 0.15 on the heavy all-electron atoms: with the core removed the local energy varies on the valence scale instead. The soft coordinate returns \(L_{r^2}\) of 1.0 to 2.0 bohr on the light systems, the size of the electron cloud, and nothing usable on argon and krypton, where \(n_\mathrm{corr}(r^2)\) reaches \(10^4\) and \(10^5\) steps per point cannot resolve it.

The Gaussian and Slater versions of the same system are the control on all of this, and they agree: over the seven pairs measured, \(n_\mathrm{corr}\), \(L_E^2\) and \(p_\mathrm{opt}\) differ by 5 % to 25 % with no systematic sign. What is measured is a property of the wave function, not of how the orbitals are expanded.

One system has no interior minimum at all. Pseudo-hydrogen carries one electron and a nearly exact wave function — \(\mathrm{Var}(E_L) = 0.0011\) — and its \(n_\mathrm{corr}\) rises monotonically with the step, 1.9 to 23, so the apparent optimum at 84 % acceptance is the edge of the grid rather than a measurement. That is the second mechanism above failing to exist: with one electron there is no configuration space for the acceptance to be inhomogeneous over.

Three further reservations on data taken at \(10^5\) steps per point. Rows above 90 % acceptance are unusable — C₄H₄ reports \(n_\mathrm{corr} = 1764 \pm 1048\) there and the optimal period saturates its cap of 100. A single close encounter still contaminates one row in twenty, visible as a variance three to four times its neighbours’, and it depresses the correlation time of the same row. And \(A_\mathrm{opt}\) itself is the argument of a maximum of a flat curve, so it moves by \(\pm 0.1\) between repeat runs: the ratios above are the robust quantities, not the position of the optimum.

References

[22] (1,2)

G. O. Roberts, A. Gelman, and W. R. Gilks, Weak convergence and optimal scaling of random walk Metropolis algorithms, Ann. Appl. Probab. 7, 110 (1997).

[23] (1,2,3,4)

R. M. Lee, G. J. Conduit, N. Nemec, P. López Ríos, and N. D. Drummond, Strategies for improving the efficiency of quantum Monte Carlo calculations, Phys. Rev. E 83, 066706 (2011).

[24] (1,2)

E. H. Lieb, Thomas-Fermi and related theories of atoms and molecules, Rev. Mod. Phys. 53, 603 (1981).

[25]

J. M. C. Scott, The binding energy of the Thomas-Fermi atom, Philos. Mag. 43, 859 (1952).

[26]

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).