From 5b0bb68fee4417c7a9aee568421496a94ae48adf Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Fri, 7 Aug 2026 23:50:10 +0900
Subject: [PATCH 01/27] Add variational Monte Carlo solver

---
 Project.toml      |   3 +
 docs/Project.toml |   3 +
 docs/make.jl      |   1 +
 docs/src/VMC.md   |  62 ++++++++++++++
 src/TwoBody.jl    |   1 +
 src/VMC.jl        | 205 ++++++++++++++++++++++++++++++++++++++++++++++
 test/Project.toml |   4 +
 test/VMC.jl       |  80 ++++++++++++++++++
 test/runtests.jl  |   3 +-
 9 files changed, 361 insertions(+), 1 deletion(-)
 create mode 100644 docs/src/VMC.md
 create mode 100644 src/VMC.jl
 create mode 100644 test/VMC.jl

diff --git a/Project.toml b/Project.toml
index 6b23157..70490c5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,9 +6,11 @@ version = "0.0.9"
 [deps]
 ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d"
 FiniteDifferenceMatrices = "a7a66f33-e7b8-47af-b618-f9b5bea05f3d"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Subscripts = "2b7f82d5-8785-4f63-971e-f18ddbeb808e"
@@ -16,6 +18,7 @@ Subscripts = "2b7f82d5-8785-4f63-971e-f18ddbeb808e"
 [compat]
 ArnoldiMethod = "0.4.0"
 FiniteDifferenceMatrices = "0.1.0"
+ForwardDiff = "0.10, 1"
 Optim = "1.9.4"
 SpecialFunctions = "2.3.1"
 Subscripts = "0.1.3"
diff --git a/docs/Project.toml b/docs/Project.toml
index 106e3c9..1a19e6e 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -3,3 +3,6 @@ Antique = "be6e5d0e-34a5-4c8f-af83-e1b5389203d8"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 TwoBody = "a92d7657-722c-45a6-9d18-9da4c8a753b6"
+
+[compat]
+Antique = "0.13"
diff --git a/docs/make.jl b/docs/make.jl
index 57353b7..c8cbee8 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -20,6 +20,7 @@ makedocs(;
         "Hamiltonian" => "Hamiltonian.md",
         "Rayleigh-Ritz Method" => "Rayleigh-Ritz.md",
         "Finite Difference Method" => "FDM.md",
+        "Variational Monte Carlo" => "VMC.md",
         "API reference" => "API.md",
     ],
 )
diff --git a/docs/src/VMC.md b/docs/src/VMC.md
new file mode 100644
index 0000000..68d3c65
--- /dev/null
+++ b/docs/src/VMC.md
@@ -0,0 +1,62 @@
+```@meta
+CurrentModule = TwoBody
+```
+
+# Variational Monte Carlo
+
+Variational Monte Carlo estimates the energy of a trial wavefunction from the
+local energy
+
+```math
+E_\mathrm{loc}(\mathbf{r}) =
+\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
+```
+
+at positions sampled from ``|\psi(\mathbf{r})|^2``. The Laplacian in a
+non-relativistic kinetic term is evaluated by automatic differentiation with
+ForwardDiff.jl.
+
+## Usage
+
+The following example uses a one-Gaussian trial wavefunction for the hydrogen
+atom. A nonzero initial position avoids starting exactly at the Coulomb
+singularity.
+
+```julia
+using Random
+using TwoBody
+
+H = Hamiltonian(
+  NonRelativisticKinetic(ℏ=1, m=1),
+  CoulombPotential(coefficient=-1),
+)
+
+α = 0.2829
+ψ(r) = exp(-α * sum(abs2, r))
+
+method = VariationalMonteCarlo(
+  n_steps=100_000,
+  burn_in=1_000,
+  δ=0.5,
+  r₀=[1.0, 0.0, 0.0],
+)
+
+result = solve(H, ψ, method; rng=MersenneTwister(123))
+result.E
+```
+
+The result also contains the sample variance, a naive standard error, the
+acceptance rate, local energies, and sampled positions. Because
+successive Markov-chain samples are correlated, use batching or an autocorrelation
+analysis when a rigorous uncertainty estimate is required.
+
+Non-finite local energies at isolated singular points are excluded from the
+average and counted in `result.n_discarded`.
+
+## API reference
+
+```@docs; canonical=false
+TwoBody.VariationalMonteCarlo
+TwoBody.local_energy
+TwoBody.solve(hamiltonian::Hamiltonian, wavefunction::Function, method::VariationalMonteCarlo)
+```
diff --git a/src/TwoBody.jl b/src/TwoBody.jl
index a455922..68eabda 100644
--- a/src/TwoBody.jl
+++ b/src/TwoBody.jl
@@ -9,5 +9,6 @@ include("./Basis.jl")
 # Solvers
 include("./Rayleigh-Ritz.jl")
 include("./FDM.jl")
+include("./VMC.jl")
 
 end
diff --git a/src/VMC.jl b/src/VMC.jl
new file mode 100644
index 0000000..9867ebf
--- /dev/null
+++ b/src/VMC.jl
@@ -0,0 +1,205 @@
+export VariationalMonteCarlo, local_energy
+
+import ForwardDiff
+import LinearAlgebra
+import Random
+
+struct VariationalMonteCarlo{T<:AbstractFloat,V<:AbstractVector{T}}
+  n_steps::Int
+  burn_in::Int
+  thinning::Int
+  δ::T
+  r₀::V
+
+  function VariationalMonteCarlo(;
+    n_steps::Int=10^5,
+    burn_in::Int=10^3,
+    thinning::Int=1,
+    δ::Real=0.5,
+    r₀::AbstractVector{<:Real}=[1.0, 0.0, 0.0],
+  )
+    0 < n_steps || throw(ArgumentError("n_steps must be positive"))
+    0 ≤ burn_in || throw(ArgumentError("burn_in must be nonnegative"))
+    0 < thinning || throw(ArgumentError("thinning must be positive"))
+    isfinite(δ) && 0 < δ || throw(ArgumentError("δ must be positive and finite"))
+    isempty(r₀) && throw(ArgumentError("r₀ must contain at least one coordinate"))
+    all(isfinite, r₀) || throw(ArgumentError("r₀ must contain only finite coordinates"))
+
+    position = float.(collect(r₀))
+    step_size = convert(eltype(position), δ)
+    new{eltype(position),typeof(position)}(n_steps, burn_in, thinning, step_size, position)
+  end
+end
+
+Base.string(method::VariationalMonteCarlo) =
+  "VariationalMonteCarlo(" *
+  join(["$(symbol)=$(getproperty(method, symbol))" for symbol in fieldnames(typeof(method))], ", ") *
+  ")"
+Base.show(io::IO, method::VariationalMonteCarlo) = print(io, Base.string(method))
+
+function _laplacian(wavefunction::Function, position::AbstractVector{<:Real})
+  return LinearAlgebra.tr(ForwardDiff.hessian(wavefunction, position))
+end
+
+function _local_energy(term::Laplacian, wavefunction::Function, position, ψ)
+  return term.coefficient * _laplacian(wavefunction, position) / ψ
+end
+
+function _local_energy(term::NonRelativisticKinetic, wavefunction::Function, position, ψ)
+  return -term.ℏ^2 / (2 * term.m) * _laplacian(wavefunction, position) / ψ
+end
+
+_local_energy(term::RestEnergy, wavefunction::Function, position, ψ) = term.m * term.c^2
+_local_energy(term::PotentialTerm, wavefunction::Function, position, ψ) =
+  V(term, LinearAlgebra.norm(position))
+
+function _local_energy(term::KineticTerm, wavefunction::Function, position, ψ)
+  throw(ArgumentError("$(typeof(term)) is not supported by VariationalMonteCarlo"))
+end
+
+function local_energy(
+  hamiltonian::Hamiltonian,
+  wavefunction::Function,
+  position::AbstractVector{<:Real},
+)
+  ψ = wavefunction(position)
+  return sum(_local_energy(term, wavefunction, position, ψ) for term in hamiltonian.terms)
+end
+
+function _probability_density(wavefunction::Function, position)
+  probability = abs2(wavefunction(position))
+  probability isa Real || throw(ArgumentError("abs2(wavefunction(r)) must be real"))
+  return probability
+end
+
+function _metropolis_samples(
+  wavefunction::Function,
+  method::VariationalMonteCarlo,
+  rng::Random.AbstractRNG,
+)
+  position = copy(method.r₀)
+  probability = _probability_density(wavefunction, position)
+  isfinite(probability) && 0 < probability ||
+    throw(ArgumentError("the probability density at r₀ must be positive and finite"))
+
+  samples = Matrix{eltype(position)}(undef, length(position), method.n_steps)
+  accepted = 0
+  attempted = 0
+  stored = 0
+  n_transitions = method.burn_in + method.n_steps * method.thinning
+  half = eltype(position)(1 // 2)
+
+  for step in 1:n_transitions
+    proposal = position .+ method.δ .* (rand(rng, eltype(position), length(position)) .- half)
+    proposed_probability = _probability_density(wavefunction, proposal)
+    attempted += 1
+
+    if isfinite(proposed_probability) && 0 ≤ proposed_probability
+      acceptance_probability = min(one(probability), proposed_probability / probability)
+      if rand(rng) < acceptance_probability
+        position = proposal
+        probability = proposed_probability
+        accepted += 1
+      end
+    end
+
+    if method.burn_in < step && (step - method.burn_in) % method.thinning == 0
+      stored += 1
+      samples[:, stored] = position
+    end
+  end
+
+  return samples, accepted / attempted
+end
+
+function _isfinite_number(value::Number)
+  return isfinite(real(value)) && isfinite(imag(value))
+end
+
+function solve(
+  hamiltonian::Hamiltonian,
+  wavefunction::Function,
+  method::VariationalMonteCarlo;
+  rng::Random.AbstractRNG=Random.default_rng(),
+  info::Int=1,
+)
+  samples, acceptance_rate = _metropolis_samples(wavefunction, method, rng)
+  energies = [local_energy(hamiltonian, wavefunction, samples[:, i]) for i in axes(samples, 2)]
+  finite_energies = filter(_isfinite_number, energies)
+  isempty(finite_energies) &&
+    throw(ArgumentError("all local energies are non-finite; check the wavefunction and initial position"))
+
+  energy = sum(finite_energies) / length(finite_energies)
+  variance = if length(finite_energies) == 1
+    zero(abs2(first(finite_energies)))
+  else
+    sum(abs2(value - energy) for value in finite_energies) / (length(finite_energies) - 1)
+  end
+  standard_error = sqrt(variance / length(finite_energies))
+
+  result = (
+    hamiltonian=hamiltonian,
+    method=method,
+    E=energy,
+    variance=variance,
+    standard_error=standard_error,
+    acceptance_rate=acceptance_rate,
+    n_samples=length(finite_energies),
+    n_discarded=length(energies) - length(finite_energies),
+    local_energies=energies,
+    samples=samples,
+  )
+
+  if 0 < info
+    println("\n# method\n")
+    println(method)
+    println("\n# energy\n")
+    println("E = $(result.E) ± $(result.standard_error)")
+    println("\n# sampling\n")
+    println("acceptance rate = $(result.acceptance_rate)")
+    println("finite local energies = $(result.n_samples) / $(method.n_steps)")
+  end
+
+  return result
+end
+
+@doc raw"""
+`VariationalMonteCarlo(n_steps=10^5, burn_in=10^3, thinning=1, δ=0.5, r₀=[1.0, 0.0, 0.0])`
+
+Options for variational Monte Carlo with a symmetric, uniform Metropolis proposal.
+The sampler targets ``|\psi(\mathbf{r})|^2``. `n_steps` is the number of retained
+samples, `burn_in` is the number of discarded initial transitions, `thinning` is
+the number of transitions between retained samples, `δ` is the proposal-box
+width, and `r₀` is the initial position.
+""" VariationalMonteCarlo
+
+@doc raw"""
+`local_energy(hamiltonian, wavefunction, position)`
+
+Evaluate
+
+```math
+E_\mathrm{loc}(\mathbf{r}) =
+\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}.
+```
+
+The Laplacian of a real-valued trial wavefunction is evaluated with
+`ForwardDiff.hessian`. Non-relativistic kinetic, Laplacian, rest-energy, and
+potential terms with a defined `V` method are supported.
+""" local_energy
+
+@doc raw"""
+`solve(hamiltonian, wavefunction, method::VariationalMonteCarlo; rng=Random.default_rng(), info=1)`
+
+Estimate the variational energy by averaging the local energy over Metropolis
+samples from ``|\psi|^2``. A seeded random-number generator can be supplied with
+`rng`. Non-finite local energies, which can occur at a measure-zero singularity
+such as the origin of a Coulomb potential, are excluded and reported as
+`n_discarded` in the result.
+""" solve(
+  hamiltonian::Hamiltonian,
+  wavefunction::Function,
+  method::VariationalMonteCarlo;
+  rng::Random.AbstractRNG=Random.default_rng(),
+  info::Int=1,
+)
diff --git a/test/Project.toml b/test/Project.toml
index 8c54d18..5ac7edd 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -3,5 +3,9 @@ Antique = "be6e5d0e-34a5-4c8f-af83-e1b5389203d8"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[compat]
+Antique = "0.13"
diff --git a/test/VMC.jl b/test/VMC.jl
new file mode 100644
index 0000000..be520bc
--- /dev/null
+++ b/test/VMC.jl
@@ -0,0 +1,80 @@
+@testset "VMC.jl" begin
+  using Random
+
+  @testset "method validation" begin
+    @test_throws ArgumentError VariationalMonteCarlo(n_steps=0)
+    @test_throws ArgumentError VariationalMonteCarlo(burn_in=-1)
+    @test_throws ArgumentError VariationalMonteCarlo(thinning=0)
+    @test_throws ArgumentError VariationalMonteCarlo(δ=0)
+    @test_throws ArgumentError VariationalMonteCarlo(δ=Inf)
+    @test_throws ArgumentError VariationalMonteCarlo(r₀=Float64[])
+  end
+
+  @testset "local energy" begin
+    H = Hamiltonian(
+      NonRelativisticKinetic(ℏ=1, m=1),
+      PowerLawPotential(coefficient=1 / 2, exponent=2),
+    )
+    ψ(r) = exp(-sum(abs2, r) / 2)
+
+    @test local_energy(H, ψ, [0.2, -0.3, 0.4]) ≈ 1.5 atol=1e-12
+    @test_throws ArgumentError local_energy(
+      Hamiltonian(RelativisticKinetic()),
+      ψ,
+      [0.2, -0.3, 0.4],
+    )
+  end
+
+  @testset "harmonic oscillator" begin
+    H = Hamiltonian(
+      NonRelativisticKinetic(ℏ=1, m=1),
+      PowerLawPotential(coefficient=1 / 2, exponent=2),
+    )
+    ψ(r) = exp(-sum(abs2, r) / 2)
+    method = VariationalMonteCarlo(
+      n_steps=200,
+      burn_in=20,
+      thinning=2,
+      δ=1.0,
+      r₀=[1.0, 0.0, 0.0],
+    )
+    result = solve(H, ψ, method; rng=MersenneTwister(123), info=0)
+
+    @test result.E ≈ 1.5 atol=1e-12
+    @test result.n_samples == method.n_steps
+    @test result.n_discarded == 0
+    @test size(result.samples) == (3, method.n_steps)
+    @test length(result.local_energies) == method.n_steps
+    @test 0 < result.acceptance_rate < 1
+    mean_radius_squared =
+      sum(sum(abs2, result.samples[:, i]) for i in axes(result.samples, 2)) / method.n_steps
+    @test mean_radius_squared ≈ 1.5 atol=0.5
+
+    singular_H = Hamiltonian(
+      NonRelativisticKinetic(ℏ=1, m=1),
+      FunctionPotential(f=r -> r < 1 ? Inf : r^2 / 2),
+    )
+    singular_result = solve(
+      singular_H,
+      ψ,
+      method;
+      rng=MersenneTwister(123),
+      info=0,
+    )
+    @test 0 < singular_result.n_discarded < method.n_steps
+  end
+
+  @testset "Coulomb singularity" begin
+    H = Hamiltonian(
+      NonRelativisticKinetic(ℏ=1, m=1),
+      CoulombPotential(coefficient=-1),
+    )
+    α = 0.2829
+    ψ(r) = exp(-α * sum(abs2, r))
+    method = VariationalMonteCarlo(n_steps=10_000, burn_in=2_000, δ=2.0)
+    result = solve(H, ψ, method; rng=MersenneTwister(123), info=0)
+
+    @test isfinite(result.E)
+    @test result.E ≈ -0.42441317922678223 atol=0.03
+  end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 89abff9..847a49b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,4 +10,5 @@ using ForwardDiff
   include("Basis.jl")
   include("Rayleigh-Ritz.jl")
   include("FDM.jl")
-end
\ No newline at end of file
+  include("VMC.jl")
+end

From 6f49c090ef323cecbc3c46d8a27f9ea8962ae6ef Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Sat, 8 Aug 2026 02:33:26 +0900
Subject: [PATCH 02/27] Add VMC bibliography

---
 docs/src/VMC.md | 26 ++++++++++++++++++++++----
 1 file changed, 22 insertions(+), 4 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 68d3c65..cb2959b 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -12,9 +12,11 @@ E_\mathrm{loc}(\mathbf{r}) =
 \frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
 ```
 
-at positions sampled from ``|\psi(\mathbf{r})|^2``. The Laplacian in a
-non-relativistic kinetic term is evaluated by automatic differentiation with
-ForwardDiff.jl.
+at positions sampled from ``|\psi(\mathbf{r})|^2``. This local-energy
+formulation and its use in VMC are reviewed by Foulkes *et al.* [1]. Samples are
+generated with a symmetric random-walk proposal and the Metropolis acceptance
+rule [2]. The Laplacian in a non-relativistic kinetic term is evaluated by
+forward-mode automatic differentiation with ForwardDiff.jl [4].
 
 ## Usage
 
@@ -48,11 +50,27 @@ result.E
 The result also contains the sample variance, a naive standard error, the
 acceptance rate, local energies, and sampled positions. Because
 successive Markov-chain samples are correlated, use batching or an autocorrelation
-analysis when a rigorous uncertainty estimate is required.
+analysis when a rigorous uncertainty estimate is required; the blocking method
+of Flyvbjerg and Petersen [3] is one standard approach.
 
 Non-finite local energies at isolated singular points are excluded from the
 average and counted in `result.n_discarded`.
 
+## Bibliography
+
+1. W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal,
+   ["Quantum Monte Carlo simulations of solids"](https://doi.org/10.1103/RevModPhys.73.33),
+   *Reviews of Modern Physics* **73**, 33–83 (2001).
+2. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller,
+   ["Equation of State Calculations by Fast Computing Machines"](https://doi.org/10.1063/1.1699114),
+   *The Journal of Chemical Physics* **21**, 1087–1092 (1953).
+3. H. Flyvbjerg and H. G. Petersen,
+   ["Error estimates on averages of correlated data"](https://doi.org/10.1063/1.457480),
+   *The Journal of Chemical Physics* **91**, 461–466 (1989).
+4. J. Revels, M. Lubin, and T. Papamarkou,
+   ["Forward-Mode Automatic Differentiation in Julia"](https://arxiv.org/abs/1607.07892),
+   arXiv:1607.07892 (2016).
+
 ## API reference
 
 ```@docs; canonical=false

From bf4442a46866bd3f2cd9b4e10c0cd077e39a8d4d Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Sat, 8 Aug 2026 03:03:01 +0900
Subject: [PATCH 03/27] Show computed VMC energy in docs

---
 docs/src/VMC.md | 26 ++++++++++++++++++++------
 1 file changed, 20 insertions(+), 6 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index cb2959b..87b5d82 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -24,7 +24,7 @@ The following example uses a one-Gaussian trial wavefunction for the hydrogen
 atom. A nonzero initial position avoids starting exactly at the Coulomb
 singularity.
 
-```julia
+```@example vmc-hydrogen
 using Random
 using TwoBody
 
@@ -37,16 +37,30 @@ H = Hamiltonian(
 ψ(r) = exp(-α * sum(abs2, r))
 
 method = VariationalMonteCarlo(
-  n_steps=100_000,
-  burn_in=1_000,
-  δ=0.5,
+  n_steps=2_000,
+  burn_in=500,
+  δ=2.0,
   r₀=[1.0, 0.0, 0.0],
 )
 
-result = solve(H, ψ, method; rng=MersenneTwister(123))
-result.E
+result = solve(H, ψ, method; rng=MersenneTwister(123), info=0)
+gaussian_expectation = 3α / 2 - 2sqrt(2α / π)
+
+(;
+  VMC_energy=round(result.E; digits=4),
+  gaussian_expectation=round(gaussian_expectation; digits=4),
+  exact_energy=-0.5,
+  acceptance_rate=round(result.acceptance_rate; digits=3),
+)
 ```
 
+The sampled VMC energy is close to the analytical expectation value of the same
+Gaussian trial wavefunction. That expectation value is above the exact hydrogen
+ground-state energy, as required by the variational principle. Increasing the
+number of Monte Carlo samples reduces statistical noise but does not remove the
+variational bias caused by the restricted one-Gaussian trial wavefunction. A
+finite-sample VMC estimate can fluctuate to either side of its expectation value.
+
 The result also contains the sample variance, a naive standard error, the
 acceptance rate, local energies, and sampled positions. Because
 successive Markov-chain samples are correlated, use batching or an autocorrelation

From 56b918f318cec96085e052bc71cc8fb93eaf09f0 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Sat, 8 Aug 2026 06:50:46 +0900
Subject: [PATCH 04/27] Label sections in VMC example

---
 docs/src/VMC.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 87b5d82..d4c781a 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -28,14 +28,17 @@ singularity.
 using Random
 using TwoBody
 
+# Hamiltonian
 H = Hamiltonian(
   NonRelativisticKinetic(ℏ=1, m=1),
   CoulombPotential(coefficient=-1),
 )
 
+# Trial wave function
 α = 0.2829
 ψ(r) = exp(-α * sum(abs2, r))
 
+# VMC options
 method = VariationalMonteCarlo(
   n_steps=2_000,
   burn_in=500,
@@ -43,6 +46,7 @@ method = VariationalMonteCarlo(
   r₀=[1.0, 0.0, 0.0],
 )
 
+# Solve
 result = solve(H, ψ, method; rng=MersenneTwister(123), info=0)
 gaussian_expectation = 3α / 2 - 2sqrt(2α / π)
 

From a8bd0f88be9d63ed09cd797fec181ed133d88f2b Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Sat, 8 Aug 2026 07:15:25 +0900
Subject: [PATCH 05/27] Support multiple VMC walkers

---
 docs/src/VMC.md | 38 ++++++++++++++++++++
 src/VMC.jl      | 93 ++++++++++++++++++++++++++++++-------------------
 test/VMC.jl     | 29 +++++++++++++++
 3 files changed, 125 insertions(+), 35 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index d4c781a..9843440 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -74,6 +74,42 @@ of Flyvbjerg and Petersen [3] is one standard approach.
 Non-finite local energies at isolated singular points are excluded from the
 average and counted in `result.n_discarded`.
 
+## Multiple walkers and equilibration
+
+Thijssen [5, Table 12.1] describes calculations with 300 walkers, 12,000
+attempted displacements per walker, and the first 2,000 states of each walker
+discarded for equilibration. With `thinning=1`, the same sampling counts are
+specified by retaining the remaining 10,000 states per walker:
+
+```julia
+method = VariationalMonteCarlo(
+  n_walkers=300,
+  n_steps=10_000,
+  burn_in=2_000,
+  thinning=1,
+  δ=2.0,
+  r₀=[1.0, 0.0, 0.0],
+)
+
+result = solve(H, ψ, method; rng=MersenneTwister(123), info=0)
+
+result.E        # expectation value of the energy
+result.variance # sample variance of the retained local energies
+result.n_attempted          # 3_600_000 attempted displacements
+result.n_burn_in_discarded  # 600_000 states discarded for equilibration
+```
+
+Each walker performs `burn_in + n_steps * thinning` transitions, so the
+configuration above attempts 12,000 displacements per walker and retains
+3,000,000 samples in total. Samples and local energies are stored consecutively
+by walker; for example, they can be grouped as
+`reshape(result.local_energies, method.n_steps, method.n_walkers)`. Walkers are
+advanced sequentially using separate random draws from the supplied random
+number generator.
+
+`result.n_discarded` has a different meaning: it counts retained samples whose
+local energy was non-finite and therefore excluded from the energy statistics.
+
 ## Bibliography
 
 1. W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal,
@@ -88,6 +124,8 @@ average and counted in `result.n_discarded`.
 4. J. Revels, M. Lubin, and T. Papamarkou,
    ["Forward-Mode Automatic Differentiation in Julia"](https://arxiv.org/abs/1607.07892),
    arXiv:1607.07892 (2016).
+5. J. M. ティッセン 著, 松田和典, 道廣嘉隆, 谷村吉隆, 高須昌子, 吉江友照 訳,
+   『計算物理学』, 丸善出版 (2012).
 
 ## API reference
 
diff --git a/src/VMC.jl b/src/VMC.jl
index 9867ebf..ae4ed58 100644
--- a/src/VMC.jl
+++ b/src/VMC.jl
@@ -8,6 +8,7 @@ struct VariationalMonteCarlo{T<:AbstractFloat,V<:AbstractVector{T}}
   n_steps::Int
   burn_in::Int
   thinning::Int
+  n_walkers::Int
   δ::T
   r₀::V
 
@@ -15,19 +16,28 @@ struct VariationalMonteCarlo{T<:AbstractFloat,V<:AbstractVector{T}}
     n_steps::Int=10^5,
     burn_in::Int=10^3,
     thinning::Int=1,
+    n_walkers::Int=1,
     δ::Real=0.5,
     r₀::AbstractVector{<:Real}=[1.0, 0.0, 0.0],
   )
     0 < n_steps || throw(ArgumentError("n_steps must be positive"))
     0 ≤ burn_in || throw(ArgumentError("burn_in must be nonnegative"))
     0 < thinning || throw(ArgumentError("thinning must be positive"))
+    0 < n_walkers || throw(ArgumentError("n_walkers must be positive"))
     isfinite(δ) && 0 < δ || throw(ArgumentError("δ must be positive and finite"))
     isempty(r₀) && throw(ArgumentError("r₀ must contain at least one coordinate"))
     all(isfinite, r₀) || throw(ArgumentError("r₀ must contain only finite coordinates"))
 
     position = float.(collect(r₀))
     step_size = convert(eltype(position), δ)
-    new{eltype(position),typeof(position)}(n_steps, burn_in, thinning, step_size, position)
+    new{eltype(position),typeof(position)}(
+      n_steps,
+      burn_in,
+      thinning,
+      n_walkers,
+      step_size,
+      position,
+    )
   end
 end
 
@@ -77,39 +87,43 @@ function _metropolis_samples(
   method::VariationalMonteCarlo,
   rng::Random.AbstractRNG,
 )
-  position = copy(method.r₀)
-  probability = _probability_density(wavefunction, position)
-  isfinite(probability) && 0 < probability ||
-    throw(ArgumentError("the probability density at r₀ must be positive and finite"))
-
-  samples = Matrix{eltype(position)}(undef, length(position), method.n_steps)
+  total_samples = method.n_walkers * method.n_steps
+  samples = Matrix{eltype(method.r₀)}(undef, length(method.r₀), total_samples)
   accepted = 0
   attempted = 0
   stored = 0
   n_transitions = method.burn_in + method.n_steps * method.thinning
-  half = eltype(position)(1 // 2)
-
-  for step in 1:n_transitions
-    proposal = position .+ method.δ .* (rand(rng, eltype(position), length(position)) .- half)
-    proposed_probability = _probability_density(wavefunction, proposal)
-    attempted += 1
-
-    if isfinite(proposed_probability) && 0 ≤ proposed_probability
-      acceptance_probability = min(one(probability), proposed_probability / probability)
-      if rand(rng) < acceptance_probability
-        position = proposal
-        probability = proposed_probability
-        accepted += 1
+  half = eltype(method.r₀)(1 // 2)
+
+  for _ in 1:method.n_walkers
+    position = copy(method.r₀)
+    probability = _probability_density(wavefunction, position)
+    isfinite(probability) && 0 < probability ||
+      throw(ArgumentError("the probability density at r₀ must be positive and finite"))
+
+    for step in 1:n_transitions
+      proposal =
+        position .+ method.δ .* (rand(rng, eltype(position), length(position)) .- half)
+      proposed_probability = _probability_density(wavefunction, proposal)
+      attempted += 1
+
+      if isfinite(proposed_probability) && 0 ≤ proposed_probability
+        acceptance_probability = min(one(probability), proposed_probability / probability)
+        if rand(rng) < acceptance_probability
+          position = proposal
+          probability = proposed_probability
+          accepted += 1
+        end
       end
-    end
 
-    if method.burn_in < step && (step - method.burn_in) % method.thinning == 0
-      stored += 1
-      samples[:, stored] = position
+      if method.burn_in < step && (step - method.burn_in) % method.thinning == 0
+        stored += 1
+        samples[:, stored] = position
+      end
     end
   end
 
-  return samples, accepted / attempted
+  return samples, accepted, attempted
 end
 
 function _isfinite_number(value::Number)
@@ -123,7 +137,8 @@ function solve(
   rng::Random.AbstractRNG=Random.default_rng(),
   info::Int=1,
 )
-  samples, acceptance_rate = _metropolis_samples(wavefunction, method, rng)
+  samples, n_accepted, n_attempted = _metropolis_samples(wavefunction, method, rng)
+  acceptance_rate = n_accepted / n_attempted
   energies = [local_energy(hamiltonian, wavefunction, samples[:, i]) for i in axes(samples, 2)]
   finite_energies = filter(_isfinite_number, energies)
   isempty(finite_energies) &&
@@ -144,6 +159,9 @@ function solve(
     variance=variance,
     standard_error=standard_error,
     acceptance_rate=acceptance_rate,
+    n_accepted=n_accepted,
+    n_attempted=n_attempted,
+    n_burn_in_discarded=method.n_walkers * method.burn_in,
     n_samples=length(finite_energies),
     n_discarded=length(energies) - length(finite_energies),
     local_energies=energies,
@@ -157,20 +175,22 @@ function solve(
     println("E = $(result.E) ± $(result.standard_error)")
     println("\n# sampling\n")
     println("acceptance rate = $(result.acceptance_rate)")
-    println("finite local energies = $(result.n_samples) / $(method.n_steps)")
+    println("finite local energies = $(result.n_samples) / $(method.n_steps * method.n_walkers)")
   end
 
   return result
 end
 
 @doc raw"""
-`VariationalMonteCarlo(n_steps=10^5, burn_in=10^3, thinning=1, δ=0.5, r₀=[1.0, 0.0, 0.0])`
+`VariationalMonteCarlo(n_steps=10^5, burn_in=10^3, thinning=1, n_walkers=1, δ=0.5, r₀=[1.0, 0.0, 0.0])`
 
 Options for variational Monte Carlo with a symmetric, uniform Metropolis proposal.
 The sampler targets ``|\psi(\mathbf{r})|^2``. `n_steps` is the number of retained
-samples, `burn_in` is the number of discarded initial transitions, `thinning` is
-the number of transitions between retained samples, `δ` is the proposal-box
-width, and `r₀` is the initial position.
+samples per walker, `burn_in` is the number of initial transitions discarded from
+each walker, `thinning` is the number of transitions between retained samples,
+`n_walkers` is the number of Markov chains, `δ` is the proposal-box width, and
+`r₀` is the initial position of every walker. Each walker performs
+`burn_in + n_steps * thinning` transitions.
 """ VariationalMonteCarlo
 
 @doc raw"""
@@ -192,10 +212,13 @@ potential terms with a defined `V` method are supported.
 `solve(hamiltonian, wavefunction, method::VariationalMonteCarlo; rng=Random.default_rng(), info=1)`
 
 Estimate the variational energy by averaging the local energy over Metropolis
-samples from ``|\psi|^2``. A seeded random-number generator can be supplied with
-`rng`. Non-finite local energies, which can occur at a measure-zero singularity
-such as the origin of a Coulomb potential, are excluded and reported as
-`n_discarded` in the result.
+samples from ``|\psi|^2``. Samples from multiple walkers are stored consecutively
+in the columns of `result.samples`. A seeded random-number generator can be
+supplied with `rng`. Non-finite local energies, which can occur at a measure-zero
+singularity such as the origin of a Coulomb potential, are excluded and reported
+as `n_discarded` in the result. The result also reports `n_attempted`, `n_accepted`,
+and the number of equilibration transitions discarded across all walkers as
+`n_burn_in_discarded`.
 """ solve(
   hamiltonian::Hamiltonian,
   wavefunction::Function,
diff --git a/test/VMC.jl b/test/VMC.jl
index be520bc..ac146d9 100644
--- a/test/VMC.jl
+++ b/test/VMC.jl
@@ -5,6 +5,7 @@
     @test_throws ArgumentError VariationalMonteCarlo(n_steps=0)
     @test_throws ArgumentError VariationalMonteCarlo(burn_in=-1)
     @test_throws ArgumentError VariationalMonteCarlo(thinning=0)
+    @test_throws ArgumentError VariationalMonteCarlo(n_walkers=0)
     @test_throws ArgumentError VariationalMonteCarlo(δ=0)
     @test_throws ArgumentError VariationalMonteCarlo(δ=Inf)
     @test_throws ArgumentError VariationalMonteCarlo(r₀=Float64[])
@@ -64,6 +65,34 @@
     @test 0 < singular_result.n_discarded < method.n_steps
   end
 
+  @testset "multiple walkers" begin
+    H = Hamiltonian(
+      NonRelativisticKinetic(ℏ=1, m=1),
+      PowerLawPotential(coefficient=1 / 2, exponent=2),
+    )
+    ψ(r) = exp(-sum(abs2, r) / 2)
+    method = VariationalMonteCarlo(
+      n_steps=50,
+      burn_in=10,
+      n_walkers=3,
+      δ=1.0,
+      r₀=[1.0, 0.0, 0.0],
+    )
+    result = solve(H, ψ, method; rng=MersenneTwister(123), info=0)
+    total_samples = method.n_walkers * method.n_steps
+    transitions_per_walker = method.burn_in + method.n_steps * method.thinning
+
+    @test result.E ≈ 1.5 atol=1e-12
+    @test result.n_samples == total_samples
+    @test result.n_discarded == 0
+    @test size(result.samples) == (3, total_samples)
+    @test length(result.local_energies) == total_samples
+    @test result.n_attempted == method.n_walkers * transitions_per_walker
+    @test result.n_accepted == round(Int, result.acceptance_rate * result.n_attempted)
+    @test result.n_burn_in_discarded == method.n_walkers * method.burn_in
+    @test 0 < result.acceptance_rate < 1
+  end
+
   @testset "Coulomb singularity" begin
     H = Hamiltonian(
       NonRelativisticKinetic(ℏ=1, m=1),

From 75da1b0261c83573ae703eaebf36c857df97173c Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Sat, 8 Aug 2026 07:19:37 +0900
Subject: [PATCH 06/27] Unify VMC energy equations

---
 docs/src/VMC.md | 38 +++++++++++++++++++++++++++++---------
 1 file changed, 29 insertions(+), 9 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 9843440..ecc33d8 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -4,19 +4,39 @@ CurrentModule = TwoBody
 
 # Variational Monte Carlo
 
-Variational Monte Carlo estimates the energy of a trial wavefunction from the
-local energy
+Combining the variational energy expectation value with the probability density
+and local energy gives the single expression
 
 ```math
-E_\mathrm{loc}(\mathbf{r}) =
-\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
+\begin{aligned}
+\langle E \rangle
+&= \frac{
+  \displaystyle \int \mathrm{d}\mathbf{r}\,
+  \psi^*(\mathbf{r})\hat{H}\psi(\mathbf{r})
+}{
+  \displaystyle \int \mathrm{d}\mathbf{r}\,
+  |\psi(\mathbf{r})|^2
+} \\
+&= \int \mathrm{d}\mathbf{r}\,
+\underbrace{
+  \frac{|\psi(\mathbf{r})|^2}
+  {\displaystyle \int \mathrm{d}\mathbf{r}'\,|\psi(\mathbf{r}')|^2}
+}_{P(\mathbf{r})}
+\underbrace{
+  \frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
+}_{E_\mathrm{loc}(\mathbf{r})}
+= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})E_\mathrm{loc}(\mathbf{r})
+\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i),
+\qquad \mathbf{r}_i \sim P.
+\end{aligned}
 ```
 
-at positions sampled from ``|\psi(\mathbf{r})|^2``. This local-energy
-formulation and its use in VMC are reviewed by Foulkes *et al.* [1]. Samples are
-generated with a symmetric random-walk proposal and the Metropolis acceptance
-rule [2]. The Laplacian in a non-relativistic kinetic term is evaluated by
-forward-mode automatic differentiation with ForwardDiff.jl [4].
+Thus VMC estimates ``\langle E\rangle`` by averaging the local energy over
+positions sampled from the normalized density ``P(\mathbf{r})``. This
+local-energy formulation and its use in VMC are reviewed by Foulkes *et al.*
+[1]. Samples are generated with a symmetric random-walk proposal and the
+Metropolis acceptance rule [2]. The Laplacian in a non-relativistic kinetic term
+is evaluated by forward-mode automatic differentiation with ForwardDiff.jl [4].
 
 ## Usage
 

From 95e1f446449e269a8ec6ab6867e8dbfcf29896bd Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Sat, 8 Aug 2026 07:23:59 +0900
Subject: [PATCH 07/27] Use Thijssen hydrogen trial function

---
 docs/src/VMC.md | 34 +++++++++++++++++++---------------
 1 file changed, 19 insertions(+), 15 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index ecc33d8..fa88627 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -40,9 +40,9 @@ is evaluated by forward-mode automatic differentiation with ForwardDiff.jl [4].
 
 ## Usage
 
-The following example uses a one-Gaussian trial wavefunction for the hydrogen
-atom. A nonzero initial position avoids starting exactly at the Coulomb
-singularity.
+The following example uses the hydrogen trial wavefunction
+``\psi(r)=\exp(-0.8r)`` from Thijssen [5, Table 12.1]. A nonzero initial
+position avoids starting exactly at the Coulomb singularity.
 
 ```@example vmc-hydrogen
 using Random
@@ -55,35 +55,39 @@ H = Hamiltonian(
 )
 
 # Trial wave function
-α = 0.2829
-ψ(r) = exp(-α * sum(abs2, r))
+α = 0.8
+ψ(r) = exp(-α * sqrt(sum(abs2, r)))
 
 # VMC options
 method = VariationalMonteCarlo(
-  n_steps=2_000,
-  burn_in=500,
+  n_walkers=20,
+  n_steps=100,
+  burn_in=100,
   δ=2.0,
   r₀=[1.0, 0.0, 0.0],
 )
 
 # Solve
 result = solve(H, ψ, method; rng=MersenneTwister(123), info=0)
-gaussian_expectation = 3α / 2 - 2sqrt(2α / π)
+analytical_expectation = α^2 / 2 - α
 
 (;
   VMC_energy=round(result.E; digits=4),
-  gaussian_expectation=round(gaussian_expectation; digits=4),
+  analytical_expectation=round(analytical_expectation; digits=4),
+  thijssen_table_12_1="-0.4813(6)",
   exact_energy=-0.5,
   acceptance_rate=round(result.acceptance_rate; digits=3),
 )
 ```
 
-The sampled VMC energy is close to the analytical expectation value of the same
-Gaussian trial wavefunction. That expectation value is above the exact hydrogen
-ground-state energy, as required by the variational principle. Increasing the
-number of Monte Carlo samples reduces statistical noise but does not remove the
-variational bias caused by the restricted one-Gaussian trial wavefunction. A
-finite-sample VMC estimate can fluctuate to either side of its expectation value.
+For ``\psi(r)=\exp(-\alpha r)``, the analytical expectation value in atomic
+units is ``\alpha^2/2-\alpha=-0.48`` at ``\alpha=0.8``. Thijssen [5, Table 12.1]
+reports ``-0.4813(6)`` using the larger sampling counts described below. The
+expectation value is above the exact hydrogen ground-state energy, as required
+by the variational principle. Increasing the number of Monte Carlo samples
+reduces statistical noise but does not remove the variational bias of the trial
+wavefunction. A finite-sample VMC estimate can fluctuate to either side of its
+expectation value.
 
 The result also contains the sample variance, a naive standard error, the
 acceptance rate, local energies, and sampled positions. Because

From 129742bd5d5c289e423fa493a352f68bcf90b695 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Sat, 8 Aug 2026 07:41:35 +0900
Subject: [PATCH 08/27] Simplify VMC example output

---
 docs/src/VMC.md | 13 ++++---------
 1 file changed, 4 insertions(+), 9 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index fa88627..7427014 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -69,15 +69,10 @@ method = VariationalMonteCarlo(
 
 # Solve
 result = solve(H, ψ, method; rng=MersenneTwister(123), info=0)
-analytical_expectation = α^2 / 2 - α
-
-(;
-  VMC_energy=round(result.E; digits=4),
-  analytical_expectation=round(analytical_expectation; digits=4),
-  thijssen_table_12_1="-0.4813(6)",
-  exact_energy=-0.5,
-  acceptance_rate=round(result.acceptance_rate; digits=3),
-)
+
+# Display
+println("This work: $(result.E)")
+println("Reference: -0.4813(6)")
 ```
 
 For ``\psi(r)=\exp(-\alpha r)``, the analytical expectation value in atomic

From 0c9a935a02253153b7ec7bac0f6711d2cc0a25c7 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Sat, 8 Aug 2026 07:45:43 +0900
Subject: [PATCH 09/27] Make VMC defaults reproducible and quiet

---
 docs/src/VMC.md |  5 ++---
 src/VMC.jl      | 18 ++++++++++--------
 test/VMC.jl     |  5 +++++
 3 files changed, 17 insertions(+), 11 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 7427014..f6d3594 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -45,7 +45,6 @@ The following example uses the hydrogen trial wavefunction
 position avoids starting exactly at the Coulomb singularity.
 
 ```@example vmc-hydrogen
-using Random
 using TwoBody
 
 # Hamiltonian
@@ -68,7 +67,7 @@ method = VariationalMonteCarlo(
 )
 
 # Solve
-result = solve(H, ψ, method; rng=MersenneTwister(123), info=0)
+result = solve(H, ψ, method)
 
 # Display
 println("This work: $(result.E)")
@@ -110,7 +109,7 @@ method = VariationalMonteCarlo(
   r₀=[1.0, 0.0, 0.0],
 )
 
-result = solve(H, ψ, method; rng=MersenneTwister(123), info=0)
+result = solve(H, ψ, method)
 
 result.E        # expectation value of the energy
 result.variance # sample variance of the retained local energies
diff --git a/src/VMC.jl b/src/VMC.jl
index ae4ed58..cd0b3f8 100644
--- a/src/VMC.jl
+++ b/src/VMC.jl
@@ -134,8 +134,8 @@ function solve(
   hamiltonian::Hamiltonian,
   wavefunction::Function,
   method::VariationalMonteCarlo;
-  rng::Random.AbstractRNG=Random.default_rng(),
-  info::Int=1,
+  rng::Random.AbstractRNG=Random.MersenneTwister(123),
+  info::Int=0,
 )
   samples, n_accepted, n_attempted = _metropolis_samples(wavefunction, method, rng)
   acceptance_rate = n_accepted / n_attempted
@@ -209,20 +209,22 @@ potential terms with a defined `V` method are supported.
 """ local_energy
 
 @doc raw"""
-`solve(hamiltonian, wavefunction, method::VariationalMonteCarlo; rng=Random.default_rng(), info=1)`
+`solve(hamiltonian, wavefunction, method::VariationalMonteCarlo; rng=Random.MersenneTwister(123), info=0)`
 
 Estimate the variational energy by averaging the local energy over Metropolis
 samples from ``|\psi|^2``. Samples from multiple walkers are stored consecutively
 in the columns of `result.samples`. A seeded random-number generator can be
-supplied with `rng`. Non-finite local energies, which can occur at a measure-zero
-singularity such as the origin of a Coulomb potential, are excluded and reported
-as `n_discarded` in the result. The result also reports `n_attempted`, `n_accepted`,
+supplied with `rng`; by default, a new `MersenneTwister(123)` is used for a
+reproducible calculation. Set `info` to a positive value to print a result
+summary. Non-finite local energies, which can occur at a measure-zero singularity
+such as the origin of a Coulomb potential, are excluded and reported as
+`n_discarded` in the result. The result also reports `n_attempted`, `n_accepted`,
 and the number of equilibration transitions discarded across all walkers as
 `n_burn_in_discarded`.
 """ solve(
   hamiltonian::Hamiltonian,
   wavefunction::Function,
   method::VariationalMonteCarlo;
-  rng::Random.AbstractRNG=Random.default_rng(),
-  info::Int=1,
+  rng::Random.AbstractRNG=Random.MersenneTwister(123),
+  info::Int=0,
 )
diff --git a/test/VMC.jl b/test/VMC.jl
index ac146d9..6d9558c 100644
--- a/test/VMC.jl
+++ b/test/VMC.jl
@@ -91,6 +91,11 @@
     @test result.n_accepted == round(Int, result.acceptance_rate * result.n_attempted)
     @test result.n_burn_in_discarded == method.n_walkers * method.burn_in
     @test 0 < result.acceptance_rate < 1
+
+    default_result = solve(H, ψ, method)
+    repeated_default_result = solve(H, ψ, method)
+    @test default_result.samples == repeated_default_result.samples
+    @test default_result.E == repeated_default_result.E
   end
 
   @testset "Coulomb singularity" begin

From 91e5bc6b68b34449a76d2040e257eb90305a68c5 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 00:06:29 +0900
Subject: [PATCH 10/27] Break VMC equation across steps

---
 docs/src/VMC.md | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index f6d3594..9695374 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -25,8 +25,9 @@ and local energy gives the single expression
 \underbrace{
   \frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
 }_{E_\mathrm{loc}(\mathbf{r})}
-= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})E_\mathrm{loc}(\mathbf{r})
-\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i),
+\\
+&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})E_\mathrm{loc}(\mathbf{r}) \\
+&\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i),
 \qquad \mathbf{r}_i \sim P.
 \end{aligned}
 ```

From 54302e4f80b01d3949dafbe3cde38214b10af7a0 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 00:33:37 +0900
Subject: [PATCH 11/27] Show normalization step in VMC equation

---
 docs/src/VMC.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 9695374..f89a13d 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -17,6 +17,10 @@ and local energy gives the single expression
   \displaystyle \int \mathrm{d}\mathbf{r}\,
   |\psi(\mathbf{r})|^2
 } \\
+&= \frac{1}{Z}\int \mathrm{d}\mathbf{r}\,
+\psi^*(\mathbf{r})\hat{H}\psi(\mathbf{r}),
+\qquad
+Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2 \\
 &= \int \mathrm{d}\mathbf{r}\,
 \underbrace{
   \frac{|\psi(\mathbf{r})|^2}

From 79644c151996ae6c8a49833019cc4b698730515e Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 01:23:43 +0900
Subject: [PATCH 12/27] Clarify VMC probability substitution

---
 docs/src/VMC.md | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index f89a13d..407a9a6 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -20,19 +20,19 @@ and local energy gives the single expression
 &= \frac{1}{Z}\int \mathrm{d}\mathbf{r}\,
 \psi^*(\mathbf{r})\hat{H}\psi(\mathbf{r}),
 \qquad
-Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2 \\
-&= \int \mathrm{d}\mathbf{r}\,
-\underbrace{
-  \frac{|\psi(\mathbf{r})|^2}
-  {\displaystyle \int \mathrm{d}\mathbf{r}'\,|\psi(\mathbf{r}')|^2}
-}_{P(\mathbf{r})}
-\underbrace{
-  \frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
-}_{E_\mathrm{loc}(\mathbf{r})}
-\\
-&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})E_\mathrm{loc}(\mathbf{r}) \\
+\left(Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2\right) \\
+&= \frac{1}{Z}\int \mathrm{d}\mathbf{r}\,
+|\psi(\mathbf{r})|^2
+\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})} \\
+&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})E_\mathrm{loc}(\mathbf{r}),
+\qquad
+\left(
+P(\mathbf{r}) = \frac{|\psi(\mathbf{r})|^2}{Z},\quad
+E_\mathrm{loc}(\mathbf{r}) =
+\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
+\right) \\
 &\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i),
-\qquad \mathbf{r}_i \sim P.
+\qquad \left(\mathbf{r}_i \sim P\right).
 \end{aligned}
 ```
 

From a860de0b920e44128d2eb87e8801c7b251912084 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 01:36:38 +0900
Subject: [PATCH 13/27] Restore annotations in VMC equation

---
 docs/src/VMC.md | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 407a9a6..5ef3f5e 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -23,14 +23,17 @@ and local energy gives the single expression
 \left(Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2\right) \\
 &= \frac{1}{Z}\int \mathrm{d}\mathbf{r}\,
 |\psi(\mathbf{r})|^2
+\cdot
 \frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})} \\
-&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})E_\mathrm{loc}(\mathbf{r}),
-\qquad
-\left(
-P(\mathbf{r}) = \frac{|\psi(\mathbf{r})|^2}{Z},\quad
-E_\mathrm{loc}(\mathbf{r}) =
+&= \int \mathrm{d}\mathbf{r}\,
+\underbrace{
+\frac{|\psi(\mathbf{r})|^2}{Z}
+}_{P(\mathbf{r})}
+\cdot
+\underbrace{
 \frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
-\right) \\
+}_{E_\mathrm{loc}(\mathbf{r})} \\
+&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})\cdot E_\mathrm{loc}(\mathbf{r}) \\
 &\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i),
 \qquad \left(\mathbf{r}_i \sim P\right).
 \end{aligned}

From 72aeb57156268c0bfe3e8e478c08275f67ffbd24 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 01:39:12 +0900
Subject: [PATCH 14/27] Explain wavefunction modulus identity

---
 docs/src/VMC.md | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 5ef3f5e..873ea99 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -24,7 +24,10 @@ and local energy gives the single expression
 &= \frac{1}{Z}\int \mathrm{d}\mathbf{r}\,
 |\psi(\mathbf{r})|^2
 \cdot
-\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})} \\
+\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})},
+\qquad
+\left(|\psi(\mathbf{r})|^2 =
+\psi^*(\mathbf{r})\psi(\mathbf{r})\right) \\
 &= \int \mathrm{d}\mathbf{r}\,
 \underbrace{
 \frac{|\psi(\mathbf{r})|^2}{Z}

From e7f2982661251db3e09ab7e110f3038c67a1277a Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 01:50:47 +0900
Subject: [PATCH 15/27] Merge VMC probability definition step

---
 docs/src/VMC.md | 14 +++++---------
 1 file changed, 5 insertions(+), 9 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 873ea99..95cb9f0 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -21,21 +21,17 @@ and local energy gives the single expression
 \psi^*(\mathbf{r})\hat{H}\psi(\mathbf{r}),
 \qquad
 \left(Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2\right) \\
-&= \frac{1}{Z}\int \mathrm{d}\mathbf{r}\,
-|\psi(\mathbf{r})|^2
-\cdot
-\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})},
-\qquad
-\left(|\psi(\mathbf{r})|^2 =
-\psi^*(\mathbf{r})\psi(\mathbf{r})\right) \\
 &= \int \mathrm{d}\mathbf{r}\,
 \underbrace{
-\frac{|\psi(\mathbf{r})|^2}{Z}
+\frac{1}{Z}|\psi(\mathbf{r})|^2
 }_{P(\mathbf{r})}
 \cdot
 \underbrace{
 \frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
-}_{E_\mathrm{loc}(\mathbf{r})} \\
+}_{E_\mathrm{loc}(\mathbf{r})},
+\qquad
+\left(|\psi(\mathbf{r})|^2 =
+\psi^*(\mathbf{r})\psi(\mathbf{r})\right) \\
 &= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})\cdot E_\mathrm{loc}(\mathbf{r}) \\
 &\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i),
 \qquad \left(\mathbf{r}_i \sim P\right).

From 8488df16f5d4599b796d21861840733d1e63fa49 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 01:54:09 +0900
Subject: [PATCH 16/27] Use consistent VMC equation definitions

---
 docs/src/VMC.md | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 95cb9f0..1da994a 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -22,17 +22,19 @@ and local energy gives the single expression
 \qquad
 \left(Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2\right) \\
 &= \int \mathrm{d}\mathbf{r}\,
-\underbrace{
 \frac{1}{Z}|\psi(\mathbf{r})|^2
-}_{P(\mathbf{r})}
 \cdot
-\underbrace{
-\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
-}_{E_\mathrm{loc}(\mathbf{r})},
+\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})},
 \qquad
 \left(|\psi(\mathbf{r})|^2 =
 \psi^*(\mathbf{r})\psi(\mathbf{r})\right) \\
-&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})\cdot E_\mathrm{loc}(\mathbf{r}) \\
+&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})\cdot E_\mathrm{loc}(\mathbf{r}),
+\qquad
+\left(
+P(\mathbf{r}) = \frac{1}{Z}|\psi(\mathbf{r})|^2,\quad
+E_\mathrm{loc}(\mathbf{r}) =
+\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
+\right) \\
 &\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i),
 \qquad \left(\mathbf{r}_i \sim P\right).
 \end{aligned}

From 15fec55ef5358679b8677919f97233c6570960bb Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 01:59:23 +0900
Subject: [PATCH 17/27] Align VMC equation annotations

---
 docs/src/VMC.md | 24 ++++++++++--------------
 1 file changed, 10 insertions(+), 14 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 1da994a..d5d3bad 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -18,25 +18,21 @@ and local energy gives the single expression
   |\psi(\mathbf{r})|^2
 } \\
 &= \frac{1}{Z}\int \mathrm{d}\mathbf{r}\,
-\psi^*(\mathbf{r})\hat{H}\psi(\mathbf{r}),
-\qquad
-\left(Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2\right) \\
+\psi^*(\mathbf{r})\hat{H}\psi(\mathbf{r})
+&& Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2 \\
 &= \int \mathrm{d}\mathbf{r}\,
 \frac{1}{Z}|\psi(\mathbf{r})|^2
 \cdot
-\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})},
-\qquad
-\left(|\psi(\mathbf{r})|^2 =
-\psi^*(\mathbf{r})\psi(\mathbf{r})\right) \\
-&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})\cdot E_\mathrm{loc}(\mathbf{r}),
-\qquad
-\left(
-P(\mathbf{r}) = \frac{1}{Z}|\psi(\mathbf{r})|^2,\quad
+\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
+&& |\psi(\mathbf{r})|^2 =
+\psi^*(\mathbf{r})\psi(\mathbf{r}) \\
+&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})\cdot E_\mathrm{loc}(\mathbf{r})
+&& P(\mathbf{r}) = \frac{1}{Z}|\psi(\mathbf{r})|^2,\quad
 E_\mathrm{loc}(\mathbf{r}) =
 \frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
-\right) \\
-&\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i),
-\qquad \left(\mathbf{r}_i \sim P\right).
+\\
+&\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i)
+&& \mathbf{r}_i \sim P.
 \end{aligned}
 ```
 

From 3e162312641458d639beed56a5177a7e6aa79d20 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 10:34:13 +0900
Subject: [PATCH 18/27] Clarify VMC equation annotations

---
 docs/src/VMC.md | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index d5d3bad..9e5dfec 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -18,21 +18,21 @@ and local energy gives the single expression
   |\psi(\mathbf{r})|^2
 } \\
 &= \frac{1}{Z}\int \mathrm{d}\mathbf{r}\,
-\psi^*(\mathbf{r})\hat{H}\psi(\mathbf{r})
-&& Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2 \\
+\psi^*(\mathbf{r})\hat{H}\psi(\mathbf{r}),
+&& \text{where}\quad Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2 \\
 &= \int \mathrm{d}\mathbf{r}\,
 \frac{1}{Z}|\psi(\mathbf{r})|^2
 \cdot
-\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
-&& |\psi(\mathbf{r})|^2 =
+\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})},
+&& \text{using}\quad |\psi(\mathbf{r})|^2 =
 \psi^*(\mathbf{r})\psi(\mathbf{r}) \\
-&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})\cdot E_\mathrm{loc}(\mathbf{r})
-&& P(\mathbf{r}) = \frac{1}{Z}|\psi(\mathbf{r})|^2,\quad
+&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})\cdot E_\mathrm{loc}(\mathbf{r}),
+&& \text{where}\quad P(\mathbf{r}) = \frac{1}{Z}|\psi(\mathbf{r})|^2,\quad
 E_\mathrm{loc}(\mathbf{r}) =
 \frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
 \\
-&\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i)
-&& \mathbf{r}_i \sim P.
+&\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i),
+&& \text{with}\quad \mathbf{r}_i \sim P.
 \end{aligned}
 ```
 

From 443b04c85795ea157f7dc3cbe759bdf237b3cfdf Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 10:39:55 +0900
Subject: [PATCH 19/27] Use ellipses for VMC equation annotations

---
 docs/src/VMC.md | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 9e5dfec..06ac876 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -18,21 +18,21 @@ and local energy gives the single expression
   |\psi(\mathbf{r})|^2
 } \\
 &= \frac{1}{Z}\int \mathrm{d}\mathbf{r}\,
-\psi^*(\mathbf{r})\hat{H}\psi(\mathbf{r}),
-&& \text{where}\quad Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2 \\
+\psi^*(\mathbf{r})\hat{H}\psi(\mathbf{r})
+&& \qquad\cdots\qquad Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2 \\
 &= \int \mathrm{d}\mathbf{r}\,
 \frac{1}{Z}|\psi(\mathbf{r})|^2
 \cdot
-\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})},
-&& \text{using}\quad |\psi(\mathbf{r})|^2 =
+\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
+&& \qquad\cdots\qquad |\psi(\mathbf{r})|^2 =
 \psi^*(\mathbf{r})\psi(\mathbf{r}) \\
-&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})\cdot E_\mathrm{loc}(\mathbf{r}),
-&& \text{where}\quad P(\mathbf{r}) = \frac{1}{Z}|\psi(\mathbf{r})|^2,\quad
+&= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})\cdot E_\mathrm{loc}(\mathbf{r})
+&& \qquad\cdots\qquad P(\mathbf{r}) = \frac{1}{Z}|\psi(\mathbf{r})|^2,\quad
 E_\mathrm{loc}(\mathbf{r}) =
 \frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
 \\
-&\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i),
-&& \text{with}\quad \mathbf{r}_i \sim P.
+&\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i)
+&& \qquad\cdots\qquad \mathbf{r}_i \sim P.
 \end{aligned}
 ```
 

From 3b20a4306c61345a76701ce83615e504b45f617b Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 10:46:46 +0900
Subject: [PATCH 20/27] Use inline style for VMC annotations

---
 docs/src/VMC.md | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 06ac876..45bc635 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -19,20 +19,22 @@ and local energy gives the single expression
 } \\
 &= \frac{1}{Z}\int \mathrm{d}\mathbf{r}\,
 \psi^*(\mathbf{r})\hat{H}\psi(\mathbf{r})
-&& \qquad\cdots\qquad Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2 \\
+&& \qquad\cdots\qquad {\textstyle
+Z = \int \mathrm{d}\mathbf{r}\,|\psi(\mathbf{r})|^2} \\
 &= \int \mathrm{d}\mathbf{r}\,
 \frac{1}{Z}|\psi(\mathbf{r})|^2
 \cdot
 \frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
-&& \qquad\cdots\qquad |\psi(\mathbf{r})|^2 =
-\psi^*(\mathbf{r})\psi(\mathbf{r}) \\
+&& \qquad\cdots\qquad {\textstyle
+|\psi(\mathbf{r})|^2 = \psi^*(\mathbf{r})\psi(\mathbf{r})} \\
 &= \int \mathrm{d}\mathbf{r}\,P(\mathbf{r})\cdot E_\mathrm{loc}(\mathbf{r})
-&& \qquad\cdots\qquad P(\mathbf{r}) = \frac{1}{Z}|\psi(\mathbf{r})|^2,\quad
+&& \qquad\cdots\qquad {\textstyle
+P(\mathbf{r}) = \frac{1}{Z}|\psi(\mathbf{r})|^2,\quad
 E_\mathrm{loc}(\mathbf{r}) =
-\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}
+\frac{\hat{H}\psi(\mathbf{r})}{\psi(\mathbf{r})}}
 \\
 &\approx \frac{1}{N}\sum_{i=1}^{N}E_\mathrm{loc}(\mathbf{r}_i)
-&& \qquad\cdots\qquad \mathbf{r}_i \sim P.
+&& \qquad\cdots\qquad {\textstyle \mathbf{r}_i \sim P.}
 \end{aligned}
 ```
 

From b193edd5e8c218fed27115f1cd6a91825e9e6d43 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 10:54:08 +0900
Subject: [PATCH 21/27] Consolidate VMC introduction

---
 docs/src/VMC.md | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 45bc635..f6e5c9a 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -4,8 +4,7 @@ CurrentModule = TwoBody
 
 # Variational Monte Carlo
 
-Combining the variational energy expectation value with the probability density
-and local energy gives the single expression
+Combining the variational energy expectation value with the probability density and local energy gives a single expression in which VMC estimates ``\langle E\rangle`` by averaging the local energy over positions sampled from the normalized density ``P(\mathbf{r})``:
 
 ```math
 \begin{aligned}
@@ -38,9 +37,7 @@ E_\mathrm{loc}(\mathbf{r}) =
 \end{aligned}
 ```
 
-Thus VMC estimates ``\langle E\rangle`` by averaging the local energy over
-positions sampled from the normalized density ``P(\mathbf{r})``. This
-local-energy formulation and its use in VMC are reviewed by Foulkes *et al.*
+This local-energy formulation and its use in VMC are reviewed by Foulkes *et al.*
 [1]. Samples are generated with a symmetric random-walk proposal and the
 Metropolis acceptance rule [2]. The Laplacian in a non-relativistic kinetic term
 is evaluated by forward-mode automatic differentiation with ForwardDiff.jl [4].

From 67dcef3814811e4b712416275438f7998cbe3599 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 11:04:35 +0900
Subject: [PATCH 22/27] Consolidate VMC usage documentation

---
 docs/src/VMC.md | 65 ++++++-------------------------------------------
 1 file changed, 8 insertions(+), 57 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index f6e5c9a..78df29d 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -44,9 +44,7 @@ is evaluated by forward-mode automatic differentiation with ForwardDiff.jl [4].
 
 ## Usage
 
-The following example uses the hydrogen trial wavefunction
-``\psi(r)=\exp(-0.8r)`` from Thijssen [5, Table 12.1]. A nonzero initial
-position avoids starting exactly at the Coulomb singularity.
+The following example uses the hydrogen trial wavefunction ``\psi(r)=\exp(-0.8r)`` and sampling conditions from Thijssen [5, Table 12.1]: 300 walkers, 12,000 attempted displacements per walker, and the first 2,000 states of each walker discarded for equilibration. Since `n_steps` counts retained samples, the remaining 10,000 states are specified with `burn_in=2_000` and `thinning=1`. A nonzero initial position avoids the Coulomb singularity.
 
 ```@example vmc-hydrogen
 using TwoBody
@@ -63,9 +61,10 @@ H = Hamiltonian(
 
 # VMC options
 method = VariationalMonteCarlo(
-  n_walkers=20,
-  n_steps=100,
-  burn_in=100,
+  n_walkers=300,
+  n_steps=10_000,
+  burn_in=2_000,
+  thinning=1,
   δ=2.0,
   r₀=[1.0, 0.0, 0.0],
 )
@@ -78,59 +77,11 @@ println("This work: $(result.E)")
 println("Reference: -0.4813(6)")
 ```
 
-For ``\psi(r)=\exp(-\alpha r)``, the analytical expectation value in atomic
-units is ``\alpha^2/2-\alpha=-0.48`` at ``\alpha=0.8``. Thijssen [5, Table 12.1]
-reports ``-0.4813(6)`` using the larger sampling counts described below. The
-expectation value is above the exact hydrogen ground-state energy, as required
-by the variational principle. Increasing the number of Monte Carlo samples
-reduces statistical noise but does not remove the variational bias of the trial
-wavefunction. A finite-sample VMC estimate can fluctuate to either side of its
-expectation value.
-
-The result also contains the sample variance, a naive standard error, the
-acceptance rate, local energies, and sampled positions. Because
-successive Markov-chain samples are correlated, use batching or an autocorrelation
-analysis when a rigorous uncertainty estimate is required; the blocking method
-of Flyvbjerg and Petersen [3] is one standard approach.
-
-Non-finite local energies at isolated singular points are excluded from the
-average and counted in `result.n_discarded`.
-
-## Multiple walkers and equilibration
-
-Thijssen [5, Table 12.1] describes calculations with 300 walkers, 12,000
-attempted displacements per walker, and the first 2,000 states of each walker
-discarded for equilibration. With `thinning=1`, the same sampling counts are
-specified by retaining the remaining 10,000 states per walker:
-
-```julia
-method = VariationalMonteCarlo(
-  n_walkers=300,
-  n_steps=10_000,
-  burn_in=2_000,
-  thinning=1,
-  δ=2.0,
-  r₀=[1.0, 0.0, 0.0],
-)
-
-result = solve(H, ψ, method)
-
-result.E        # expectation value of the energy
-result.variance # sample variance of the retained local energies
-result.n_attempted          # 3_600_000 attempted displacements
-result.n_burn_in_discarded  # 600_000 states discarded for equilibration
-```
+For ``\psi(r)=\exp(-\alpha r)``, the analytical expectation value in atomic units is ``\alpha^2/2-\alpha=-0.48`` at ``\alpha=0.8``; Thijssen reports ``-0.4813(6)``. The expectation value obeys the variational bound, while a finite-sample estimate can fluctuate around it. More samples reduce statistical noise but not the trial wavefunction's variational bias.
 
-Each walker performs `burn_in + n_steps * thinning` transitions, so the
-configuration above attempts 12,000 displacements per walker and retains
-3,000,000 samples in total. Samples and local energies are stored consecutively
-by walker; for example, they can be grouped as
-`reshape(result.local_energies, method.n_steps, method.n_walkers)`. Walkers are
-advanced sequentially using separate random draws from the supplied random
-number generator.
+Each walker performs `burn_in + n_steps * thinning` transitions, giving 3,600,000 attempted displacements, 600,000 discarded burn-in states, and 3,000,000 retained samples. These counts are available as `result.n_attempted`, `result.n_burn_in_discarded`, and `result.n_samples`. The result also contains the sample variance, a naive standard error, the acceptance rate, local energies, and positions; walker-wise local energies can be obtained with `reshape(result.local_energies, method.n_steps, method.n_walkers)`.
 
-`result.n_discarded` has a different meaning: it counts retained samples whose
-local energy was non-finite and therefore excluded from the energy statistics.
+Because successive Markov-chain samples are correlated, rigorous uncertainty estimates require batching or autocorrelation analysis; the blocking method of Flyvbjerg and Petersen [3] is one standard approach. Retained samples with non-finite local energies are excluded from the statistics and counted separately in `result.n_discarded`.
 
 ## Bibliography
 

From 8b6f68a787c905590f9982e6c5c6a0ae30bfa345 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 12:27:01 +0900
Subject: [PATCH 23/27] Update Antique compatibility to 0.15

---
 docs/Project.toml | 2 +-
 test/Project.toml | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 1a19e6e..e5fe655 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -5,4 +5,4 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 TwoBody = "a92d7657-722c-45a6-9d18-9da4c8a753b6"
 
 [compat]
-Antique = "0.13"
+Antique = "0.15"
diff --git a/test/Project.toml b/test/Project.toml
index 5ac7edd..b1befd1 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -8,4 +8,4 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-Antique = "0.13"
+Antique = "0.15"

From eeaf1d5ab897577c1ceb98673c494afabffe2f16 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 13:58:02 +0900
Subject: [PATCH 24/27] Deduplicate VMC guide and API docs

---
 docs/src/VMC.md | 12 +++++-------
 src/VMC.jl      | 21 ++++++++++++---------
 2 files changed, 17 insertions(+), 16 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 78df29d..9eb57ee 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -4,7 +4,7 @@ CurrentModule = TwoBody
 
 # Variational Monte Carlo
 
-Combining the variational energy expectation value with the probability density and local energy gives a single expression in which VMC estimates ``\langle E\rangle`` by averaging the local energy over positions sampled from the normalized density ``P(\mathbf{r})``:
+VMC estimates ``\langle E\rangle`` by averaging the local energy over positions sampled from the normalized density ``P(\mathbf{r})``:
 
 ```math
 \begin{aligned}
@@ -37,14 +37,14 @@ E_\mathrm{loc}(\mathbf{r}) =
 \end{aligned}
 ```
 
-This local-energy formulation and its use in VMC are reviewed by Foulkes *et al.*
+This local-energy formulation is reviewed by Foulkes *et al.*
 [1]. Samples are generated with a symmetric random-walk proposal and the
 Metropolis acceptance rule [2]. The Laplacian in a non-relativistic kinetic term
 is evaluated by forward-mode automatic differentiation with ForwardDiff.jl [4].
 
 ## Usage
 
-The following example uses the hydrogen trial wavefunction ``\psi(r)=\exp(-0.8r)`` and sampling conditions from Thijssen [5, Table 12.1]: 300 walkers, 12,000 attempted displacements per walker, and the first 2,000 states of each walker discarded for equilibration. Since `n_steps` counts retained samples, the remaining 10,000 states are specified with `burn_in=2_000` and `thinning=1`. A nonzero initial position avoids the Coulomb singularity.
+The following example uses the hydrogen trial wavefunction ``\psi(r)=\exp(-0.8r)`` and sampling conditions from Thijssen [5, Table 12.1]: 300 walkers, 12,000 attempted displacements per walker, and the first 2,000 states of each walker discarded for equilibration (`burn_in=2_000`; `n_steps=10_000` counts the retained samples). A nonzero initial position avoids starting exactly at the Coulomb singularity.
 
 ```@example vmc-hydrogen
 using TwoBody
@@ -77,11 +77,9 @@ println("This work: $(result.E)")
 println("Reference: -0.4813(6)")
 ```
 
-For ``\psi(r)=\exp(-\alpha r)``, the analytical expectation value in atomic units is ``\alpha^2/2-\alpha=-0.48`` at ``\alpha=0.8``; Thijssen reports ``-0.4813(6)``. The expectation value obeys the variational bound, while a finite-sample estimate can fluctuate around it. More samples reduce statistical noise but not the trial wavefunction's variational bias.
+For ``\psi(r)=\exp(-\alpha r)``, the analytical expectation value in atomic units is ``\alpha^2/2-\alpha=-0.48`` at ``\alpha=0.8``; Thijssen reports ``-0.4813(6)``. A finite-sample estimate fluctuates around the analytical expectation value; more samples reduce statistical noise but not the trial wavefunction's variational bias.
 
-Each walker performs `burn_in + n_steps * thinning` transitions, giving 3,600,000 attempted displacements, 600,000 discarded burn-in states, and 3,000,000 retained samples. These counts are available as `result.n_attempted`, `result.n_burn_in_discarded`, and `result.n_samples`. The result also contains the sample variance, a naive standard error, the acceptance rate, local energies, and positions; walker-wise local energies can be obtained with `reshape(result.local_energies, method.n_steps, method.n_walkers)`.
-
-Because successive Markov-chain samples are correlated, rigorous uncertainty estimates require batching or autocorrelation analysis; the blocking method of Flyvbjerg and Petersen [3] is one standard approach. Retained samples with non-finite local energies are excluded from the statistics and counted separately in `result.n_discarded`.
+Because successive Markov-chain samples are correlated, rigorous uncertainty estimates require batching or autocorrelation analysis; the blocking method of Flyvbjerg and Petersen [3] is one standard approach.
 
 ## Bibliography
 
diff --git a/src/VMC.jl b/src/VMC.jl
index cd0b3f8..b3199d6 100644
--- a/src/VMC.jl
+++ b/src/VMC.jl
@@ -212,15 +212,18 @@ potential terms with a defined `V` method are supported.
 `solve(hamiltonian, wavefunction, method::VariationalMonteCarlo; rng=Random.MersenneTwister(123), info=0)`
 
 Estimate the variational energy by averaging the local energy over Metropolis
-samples from ``|\psi|^2``. Samples from multiple walkers are stored consecutively
-in the columns of `result.samples`. A seeded random-number generator can be
-supplied with `rng`; by default, a new `MersenneTwister(123)` is used for a
-reproducible calculation. Set `info` to a positive value to print a result
-summary. Non-finite local energies, which can occur at a measure-zero singularity
-such as the origin of a Coulomb potential, are excluded and reported as
-`n_discarded` in the result. The result also reports `n_attempted`, `n_accepted`,
-and the number of equilibration transitions discarded across all walkers as
-`n_burn_in_discarded`.
+samples from ``|\psi|^2``. A seeded random-number generator can be supplied with
+`rng`; by default, a new `MersenneTwister(123)` is used for a reproducible
+calculation. Set `info` to a positive value to print a result summary.
+
+The returned named tuple contains `E`, `variance`, the naive `standard_error`,
+`acceptance_rate`, the sampling counts `n_accepted`, `n_attempted`,
+`n_burn_in_discarded`, `n_samples`, and `n_discarded`, as well as
+`local_energies` and `samples`. Data from multiple walkers are stored
+consecutively, with positions in the columns of `samples`. Non-finite local
+energies, which can occur at a measure-zero singularity such as the origin of a
+Coulomb potential, are excluded from the energy statistics and counted in
+`n_discarded`.
 """ solve(
   hamiltonian::Hamiltonian,
   wavefunction::Function,

From 794ace6c3d029e50ab2d5e4fdf03935d4d265d6d Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 14:04:52 +0900
Subject: [PATCH 25/27] Link VMC usage to API details

---
 docs/src/VMC.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 9eb57ee..a005de9 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -79,6 +79,8 @@ println("Reference: -0.4813(6)")
 
 For ``\psi(r)=\exp(-\alpha r)``, the analytical expectation value in atomic units is ``\alpha^2/2-\alpha=-0.48`` at ``\alpha=0.8``; Thijssen reports ``-0.4813(6)``. A finite-sample estimate fluctuates around the analytical expectation value; more samples reduce statistical noise but not the trial wavefunction's variational bias.
 
+`solve` also returns sampling diagnostics and retained data; see the [API reference](#API-reference) for details.
+
 Because successive Markov-chain samples are correlated, rigorous uncertainty estimates require batching or autocorrelation analysis; the blocking method of Flyvbjerg and Petersen [3] is one standard approach.
 
 ## Bibliography

From 96af95091b6caef1c025f5de04ec9ef1a7a89dbf Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 14:23:06 +0900
Subject: [PATCH 26/27] Show exact VMC energy in example

---
 docs/src/VMC.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index a005de9..7b605a9 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -75,6 +75,7 @@ result = solve(H, ψ, method)
 # Display
 println("This work: $(result.E)")
 println("Reference: -0.4813(6)")
+println("Exact    : -0.4800")
 ```
 
 For ``\psi(r)=\exp(-\alpha r)``, the analytical expectation value in atomic units is ``\alpha^2/2-\alpha=-0.48`` at ``\alpha=0.8``; Thijssen reports ``-0.4813(6)``. A finite-sample estimate fluctuates around the analytical expectation value; more samples reduce statistical noise but not the trial wavefunction's variational bias.

From d0b18df210182d7159c51a538276da21323dd431 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Mon, 10 Aug 2026 17:13:55 +0900
Subject: [PATCH 27/27] Apply VMC review feedback

---
 docs/src/VMC.md | 8 +++++---
 src/VMC.jl      | 3 ++-
 test/VMC.jl     | 3 +++
 3 files changed, 10 insertions(+), 4 deletions(-)

diff --git a/docs/src/VMC.md b/docs/src/VMC.md
index 7b605a9..f07c22f 100644
--- a/docs/src/VMC.md
+++ b/docs/src/VMC.md
@@ -82,7 +82,7 @@ For ``\psi(r)=\exp(-\alpha r)``, the analytical expectation value in atomic unit
 
 `solve` also returns sampling diagnostics and retained data; see the [API reference](#API-reference) for details.
 
-Because successive Markov-chain samples are correlated, rigorous uncertainty estimates require batching or autocorrelation analysis; the blocking method of Flyvbjerg and Petersen [3] is one standard approach.
+Because successive Markov-chain samples are correlated, the naive `standard_error` reported by `solve` typically underestimates the true uncertainty; rigorous estimates require batching or autocorrelation analysis, such as the blocking method of Flyvbjerg and Petersen [3].
 
 ## Bibliography
 
@@ -98,8 +98,10 @@ Because successive Markov-chain samples are correlated, rigorous uncertainty est
 4. J. Revels, M. Lubin, and T. Papamarkou,
    ["Forward-Mode Automatic Differentiation in Julia"](https://arxiv.org/abs/1607.07892),
    arXiv:1607.07892 (2016).
-5. J. M. ティッセン 著, 松田和典, 道廣嘉隆, 谷村吉隆, 高須昌子, 吉江友照 訳,
-   『計算物理学』, 丸善出版 (2012).
+5. J. M. Thijssen,
+   [*Computational Physics*, 2nd ed.](https://doi.org/10.1017/CBO9781139171397),
+   Cambridge University Press (2007);
+   邦訳: 松田和典, 道廣嘉隆, 谷村吉隆, 高須昌子, 吉江友照 訳, 『計算物理学』, 丸善出版 (2012).
 
 ## API reference
 
diff --git a/src/VMC.jl b/src/VMC.jl
index b3199d6..1071baf 100644
--- a/src/VMC.jl
+++ b/src/VMC.jl
@@ -216,7 +216,8 @@ samples from ``|\psi|^2``. A seeded random-number generator can be supplied with
 `rng`; by default, a new `MersenneTwister(123)` is used for a reproducible
 calculation. Set `info` to a positive value to print a result summary.
 
-The returned named tuple contains `E`, `variance`, the naive `standard_error`,
+The returned named tuple echoes the input `hamiltonian` and `method` and
+contains `E`, `variance`, the naive `standard_error`,
 `acceptance_rate`, the sampling counts `n_accepted`, `n_attempted`,
 `n_burn_in_discarded`, `n_samples`, and `n_discarded`, as well as
 `local_energies` and `samples`. Data from multiple walkers are stored
diff --git a/test/VMC.jl b/test/VMC.jl
index 6d9558c..9079d43 100644
--- a/test/VMC.jl
+++ b/test/VMC.jl
@@ -47,6 +47,7 @@
     @test size(result.samples) == (3, method.n_steps)
     @test length(result.local_energies) == method.n_steps
     @test 0 < result.acceptance_rate < 1
+    @test result.n_attempted == method.burn_in + method.n_steps * method.thinning
     mean_radius_squared =
       sum(sum(abs2, result.samples[:, i]) for i in axes(result.samples, 2)) / method.n_steps
     @test mean_radius_squared ≈ 1.5 atol=0.5
@@ -63,6 +64,8 @@
       info=0,
     )
     @test 0 < singular_result.n_discarded < method.n_steps
+    @test isfinite(singular_result.E)
+    @test singular_result.n_samples + singular_result.n_discarded == method.n_steps
   end
 
   @testset "multiple walkers" begin