Metropolis Monte Carlo: mathematics and assumptions#

Metropolis sampling moves an existing non-overlapping packing at fixed particle count, radii, composition, and box volume. Its postulate is the canonical equilibrium distribution for rigid hard spheres. An RSA packing can supply the starting centres, but equilibration must be assessed before interpreting samples as equilibrium measurements. See Monte-Carlo packings for API usage and RSA Monte Carlo: mathematics and assumptions for irreversible deposition.

Hard-sphere potential and equilibrium measure#

For species radii \(a_i,a_j\), the additive hard-sphere potential is

\[\begin{split}u_{ij}(r)=\begin{cases} +\infty,&r<a_i+a_j,\\ 0,&r\geq a_i+a_j. \end{cases}\end{split}\]

For a configuration \(\mathbf{X}\), the potential energy \(U\) sums these pair interactions and any wall exclusions. Its configurational probability density is

\[\pi(\mathbf{X})=Z^{-1}\exp[-\beta U(\mathbf{X})],\qquad \beta=(k_BT)^{-1}.\]

\(Z\) normalizes probability over allowed configurations. Consequently, all valid configurations have the same density, and overlapping states have zero probability. Temperature cancels from this purely hard-sphere acceptance rule. Packing fraction stays fixed because particle radii and volume do.

Symmetric proposals and acceptance#

PackLab proposes one displacement per particle per sweep, in particle index order. Each component is sampled independently from \(\mathrm{Uniform}[-d,d]\), where \(d\) is maximum_displacement. Thus \(d\) is a coordinate bound, not a bound on the vector length, whose maximum is \(\sqrt{3}d\).

A symmetric proposal has equal forward and reverse probability density. The single-particle Metropolis rule is then

\[P_{\mathrm{accept}}(\mathbf{X}\to\mathbf{X}')= \min\{1,\exp[-\beta(U(\mathbf{X}')-U(\mathbf{X}))]\}.\]

Starting from a valid packing, this is one if the new centre satisfies walls and all pair exclusions, and zero otherwise. A rejected move retains the previous centre. Rejected states must still contribute to time-indexed measurements; measuring only accepted moves changes the sampled distribution. Each particle update obeys detailed balance (equal equilibrium probability flow forward and backward); their sequential sweep preserves the equilibrium measure, although the composed sweep need not itself obey detailed balance.

Reaching equilibrium also requires adequate exploration, or mixing. Local moves can mix very slowly in crowded states. A valid configuration or high acceptance rate alone does not establish equilibration. Small proposals are often accepted but explore slowly; large proposals can be mostly rejected.

Periodic geometry and accumulated displacement#

A periodic image is a particle translated by integer box lengths. The minimum-image convention finds the nearest such copy for overlap checks; its formula and a worked example are in RSA Monte Carlo: mathematics and assumptions. Spheres must satisfy \(2a\leq\min(L_x,L_y,L_z)\) to avoid their own periodic images. Stored centres are wrapped into the primary box after accepted periodic moves.

Wrapped centres are suitable for pair distances, but lose the number of times a particle has crossed a face. A winding counts these traversals. Unwrapped displacement retains every accepted proposal vector:

\[\mathbf{D}_n(s)=\sum_{\ell\leq s} I_{n\ell}\,\delta\mathbf{x}_{n\ell},\qquad \mathrm{MSD}(s;s_0)=\frac{1}{N}\sum_n \left\|\mathbf{D}_n(s)-\mathbf{D}_n(s_0)\right\|^2.\]

\(I_{n\ell}\) is one for an accepted move and zero for a rejection. simulator.cumulative_displacements stores \(\mathbf{D}_n\) as a length quantity of shape (N, 3) and clears it on reset. run_metropolis_diagnostics measures mean squared displacement (MSD) relative to the start of that call, even when the simulator has already run. Tracking the actual proposal preserves travel even across several box lengths in one accepted move. Minimum-image subtraction of two endpoints would lose this information and artificially bound the MSD.

MSD here measures movement in Monte Carlo sweeps. Metropolis proposals are numerical sampling steps, not Newtonian or Brownian time evolution. A physical diffusion coefficient requires a separately justified dynamical model. Because MSD relative to the beginning of the run typically changes with time, it is a mobility diagnostic rather than a stationary equilibrium observable. By default the diagnostic report contains MSD and move statistics, and its observable and uncertainty fields are None. Supply an observable callable to request uncertainty estimates for a scalar quantity chosen to be approximately stationary after burn-in. Choosing an observable does not itself establish stationarity or equilibration.

Burn-in, autocorrelation, and uncertainty#

Burn-in is an initial interval excluded from measurements because the chain still remembers its starting state. Stationarity means the measured observable’s distribution no longer changes with sampling time. Check several observables, different initial states, and longer runs; no fixed burn-in count guarantees stationarity.

For a stationary scalar observable \(A_t\), define the lag correlation and integrated autocorrelation time (in retained-sample intervals) by

\[\rho_A(\ell)=\frac{\operatorname{Cov}(A_t,A_{t+\ell})} {\operatorname{Var}(A_t)},\qquad \tau_{\mathrm{int}}=\frac12+\sum_{\ell=1}^{\infty}\rho_A(\ell),\qquad N_{\mathrm{eff}}\simeq\frac{M}{2\tau_{\mathrm{int}}}.\]

\(M\) is the number of retained observations, and \(N_{\mathrm{eff}}\) is the effective sample size for that observable. It is not the particle count. PackLab truncates the correlation sum at its first non-positive lag, so estimates can be unreliable for short or complex trajectories. Block averaging divides observations into complete blocks; if block means are approximately independent, the standard error is their sample standard deviation divided by the square root of the block count. PackLab’s interval uses a normal quantile; few blocks can give poor coverage. If every retained value is identical, autocorrelation and uncertainty are reported as unavailable: the trace alone cannot distinguish a truly constant quantity from a chain that has failed to explore its distribution.

For example, after constructing a periodic simulator, the cosine component of its longest-wavelength density mode can be monitored as follows:

import numpy as np
from PackLab import monte_carlo

box_length_m = simulator.domain.length_x.to("meter").magnitude

def density_mode(configuration):
    x = configuration.positions.to("meter").magnitude[:, 0]
    return float(np.cos(2 * np.pi * x / box_length_m).mean())

report = monte_carlo.run_metropolis_diagnostics(
    simulator, number_of_sweeps=10_000, sample_interval=10,
    burn_in_sweeps=2_000, observable=density_mode,
    observable_name="density_mode_x", block_size=50,
)
print(report.effective_sample_size, report.confidence_interval)

This is an illustration, not an established burn-in or block size: inspect the observable, compare independent chains, and check block-size sensitivity. Statistics use only regularly spaced post-burn-in samples; a final shorter interval remains in the recorded trajectory but is excluded from them. The report does not prove equilibration. Periodic boundaries reduce wall effects but preserve finite-box effects; compare larger boxes and independent seeded trajectories before drawing bulk conclusions. Equilibrated Metropolis and PY comparisons test an approximation; RSA and PY can differ for physical reasons even with matched inputs.

Metropolis Monte Carlo: mathematics and assumptions#

Hard-sphere potential and equilibrium measure#

Symmetric proposals and acceptance#

Periodic geometry and accumulated displacement#

Burn-in, autocorrelation, and uncertainty#

Further reading#