Rejection-free Monte Carlo#

spheropack.rejection_free implements the event-driven Monte Carlo method of Peters and de With [PdW12]. It samples the canonical distribution \(\exp(-U/kT)\) of particles with a pair potential, without ever rejecting a move.

The method#

All particles move along straight lines, \(\mathrm{d}\mathbf x_i/\mathrm{d}s = \mathbf v_i\), where \(s\) plays the role of time. For every pair potential \(U_\alpha\) separately, the increase of \(U_\alpha\) along the path is accumulated; decreases are free. The probability that no reflection due to \(U_\alpha\) has happened yet is

\[ P_\mathrm{no\,refl}(s) = \exp\left[-\beta\int_{s_0}^{s}\max\left(\frac{\mathrm{d}U_\alpha}{\mathrm{d}s'}, 0\right)\mathrm{d}s'\right], \]

so a pair reflects where the accumulated increase equals \(-kT\ln u\) with \(u\) uniform in \((0, 1]\). A reflection is an elastic collision of two equal masses along the line of centres, which leaves \(\sum_i |\mathbf v_i|^2\) unchanged. Positions sampled at equidistant times \(s_n = n\,\Delta s\) follow the canonical distribution. The velocities are not physical momenta: they only set the directions and the pace of the moves.

For an isotropic pair potential the separation along a straight path first decreases to the closest approach and then increases. With at most one minimum \(r_m\) of \(U(r)\), the uphill parts are the inward path inside \(r_m\) (repulsion) and the outward path between \(r_m\) and the cutoff (attraction), so the reflection point follows from inverting \(U(r)\) on one branch and solving a quadratic for the time.

As in the packing engine, every particle keeps one predicted event in a priority queue. Because a pair can be predicted by either particle, and again later, the random number \(u\) of a pair is a counter-based hash of the two particle indices, the number of velocity changes of each, and the periodic image, and the accumulated increase is counted from the start of the pair’s current straight segment. All predictions of the same pair therefore agree; independent draws would double the reflection rate.

Usage#

import spheropack as sp
from spheropack import rejection_free as rf

system = rf.System(n=1000, density=0.317, potential=rf.LennardJones(cutoff=2.5), kT=1.085, seed=1)
system.run(100.0)  # equilibrate
g = rf.radial_distribution(system, n_samples=300, interval=2.0, r_max=5.0)
for x in system.samples(100, interval=1.0):  # positions every 1.0
    ...
system.snapshot().to_xyz("lj.xyz")  # analysis and output as for packings

Potentials: LennardJones() (truncated and shifted), WCA(), DPD(), SoftSpheres() (harmonic, Hertzian) and HardSpheres(). The box must be periodic with edges of at least twice the cutoff. For strongly repulsive potentials start from a packing, rf.System(sp.pack(...).positions, box=..., ...), rather than from random positions.

Event chains#

The paper’s second implementation moves one particle at a time (rf.System(..., method="event_chain")). The active particle moves along a coordinate axis; for every pair with it the uphill part of the potential along the path is accumulated as above, and at a reflection the active particle stops and the partner moves on in the same direction (a lift). A chain ends after a displacement chain_length, and the next chain starts at a random particle. No event queue is needed, which makes this variant much faster per unit of displacement.

By default the chains are irreversible: the direction cycles through \(+x, +y, +z\). This breaks detailed balance but satisfies global balance, so the canonical distribution is still sampled, and the motion is ballistic rather than diffusive, which speeds up the decorrelation considerably. With irreversible=False every chain gets a random axis and sign (detailed balance). Near the critical point, as in the Lennard-Jones example of the paper, the reversible chains equilibrate the long-range density fluctuations slowly; use the default there.

In both methods run(t) moves every particle by t on average (the event-chain method performs chains with a total displacement \(N t\)), so sampling intervals are comparable.

Validation#

  • Two particles. The distance distribution of two particles in a periodic box is \(r^{D-1}\exp(-U/kT)\); the tests compare the sampled histogram and the bound fraction for every potential in 2D and 3D, and check that a 20% error in the temperature would be detected.

  • Lennard-Jones fluid at \(\rho = 0.317\), \(T = 1.085\), \(r_c = 2.5\sigma\) (the state point of the paper): \(g(r)\) agrees with the paper’s long Metropolis runs within 0.002 on average for both methods (see the notebook).

  • Methods against each other at \(T = 2\): the mean potential energy per particle of the collision method and of irreversible and reversible event chains agree within about two standard errors (0.1%).

  • DPD liquid at \(\rho = 3\), \(a = 25\,kT\): \(g(r)\) agrees with the paper’s data.

Not (yet) included#

Potentials with several minima, tabulated potentials, and containers with walls.