C++ API#

The core is header-only (include/spheropack/) and can be used without Python.

namespace spheropack#

Typedefs

template<int D> using Vec = std::array<double, D>#

Fixed-size vector of D doubles.

Template Parameters:: D – dimension (2 or 3)

Enums

enum class Status#

Reason why LSPacking::run() returned.

Values:

enumerator running#: run() has not finished (initial value)

enumerator target#: radii reached a_i * target_scale

enumerator pressure#: reduced pressure exceeded max_pressure (jammed or crystallised)

enumerator jammed#: jamming protocol reached jammed_pressure (median over spheres)

enumerator stall#: radii grew by less than stall_tol over stall_window collisions per particle

enumerator collisions#: max_collisions reached

enumerator time_limit#: max_time (simulation time) reached

enumerator timeout#: wall-clock timeout reached

enumerator interrupted#: the interrupt callback returned true (e.g. Ctrl-C)

enumerator no_events#: no future events (all particles at rest)

enumerator box_limit#: the largest sphere touches its own periodic image

enum class CollisionRule#

Collision rule for growing spheres.

Here \(u_n\) is the normal relative velocity before contact (negative when approaching) and the growth speed is the rate at which the contact distance grows.

Values:

enumerator legacy#

Normal relative velocity after contact: max(-u_n, growth speed + margin).

Elastic when the approach is faster than the growth, perfectly inelastic in the frame of the growing surfaces otherwise (the legacy rule).

enumerator elastic_growing#

Elastic in the frame of the growing surfaces: u_n’ = 2 * growth speed - u_n (Donev, Torquato and Stillinger).

Heats the system; the thermostat removes it.

Functions

inline double ball_measure(int m, double r) noexcept#

Measure of a ball of radius r in m dimensions (m = 1, 2, 3).

Returns:: \(2r\), \(\pi r^2\) or \(\tfrac{4}{3}\pi r^3\)

template<int D> constexpr double ball_volume(double r) noexcept#

Volume (D=3) or area (D=2) of a ball of radius r.

Template Parameters:: D – dimension (2 or 3)
Parameters:: r – radius
Returns:: \(\pi r^2\) for D=2, \(\tfrac{4}{3}\pi r^3\) for D=3

inline const char *to_string(Status s) noexcept#

Name of a status as used in the Python interface.

Parameters:: s – status value
Returns:: the enumerator name as a string, e.g. "time_limit"; "unknown" for invalid values

template<int D, class F> void for_each_pair(const std::vector<Vec<D>> &x, const Box<D> &box, double cutoff, F &&f)#

Calls f(i, j, r) once for every pair i < j and every periodic image of j with |r| < cutoff, where r = x_j + shift * L - x_i.

Positions must lie inside the box. On periodic axes the cutoff may exceed half the box length but not the box length itself; every image within the cutoff is then reported separately. Images of a particle with itself are not reported.

Template Parameters:

D – dimension (2 or 3)

Parameters:

x – positions, inside the box
box – container
cutoff – largest distance of interest
f – callable f(std::uint32_t i, std::uint32_t j, const Vec<D>& r)

inline double first_contact(double A, double B, double C) noexcept#

Earliest \(\tau \ge 0\) at which the quadratic \(f(\tau) = A\tau^2 + 2B\tau + C\) reaches zero while decreasing.

Here \(f\) is the squared distance of a pair minus its squared contact distance, so the root is the first contact. \(C \le 0\) means contact (or round-off overlap) now. The roots use the cancellation-free form \(q = -(B + \mathrm{sign}(B)\sqrt{B^2 - AC})\), with roots \(q/A\) and \(C/q\). With growing particles \(A < 0\) is legal.

Parameters:

A – coefficient of \(\tau^2\)
B – half the coefficient of \(\tau\) ( \(B<0\): pair closing now)
C – value at \(\tau = 0\)

Returns:

contact time, 0 if already in contact, or never if the pair does not meet

inline double first_contact_linear(double gap, double rate) noexcept#

Earliest \(\tau \ge 0\) at which \(\mathrm{gap} + \mathrm{rate}\,\tau\) reaches zero (flat wall, gap linear in time).

Parameters:

gap – current distance to the wall
rate – rate of change of the gap

Returns:

contact time (0 if the gap is already non-positive and closing), or never if the gap does not decrease

inline double reflection_time(double r2, double b, double u2, double du, const PairPotential &pot) noexcept#

Time \(\tau \ge 0\) along a straight relative path at which a pair reflects.

The pair separation is \(\mathbf r(\tau) = \mathbf r_0 + \mathbf u\,\tau\). Along the inward leg (up to the closest approach) the potential rises inside \(r_m\), along the outward leg it rises between \(r_m\) and \(r_c\). The reflection happens where the accumulated rise equals du.

Parameters:

r2 – \(|\mathbf r_0|^2\)
b – \(\mathbf r_0 \cdot \mathbf u\)
u2 – \(|\mathbf u|^2\)
du – drawn energy, \(-kT \ln u\)
pot – the pair potential

Returns:

\(\tau\), or never if the pair does not reflect before it separates beyond the cutoff

inline bool may_reflect(double r2, double b, double u2, const PairPotential &pot) noexcept#

Whether a straight relative path has any uphill part for pot (a cheap necessary condition for reflection_time() to be finite).

Parameters:

r2 – \(|\mathbf r_0|^2\)
b – \(\mathbf r_0 \cdot \mathbf u\)
u2 – \(|\mathbf u|^2\)
pot – the pair potential

constexpr std::uint64_t splitmix64(std::uint64_t &state) noexcept#

One step of SplitMix64.

Parameters:: state – generator state, advanced by the golden-ratio increment
Returns:: a well-mixed 64-bit value

constexpr std::uint64_t mix64(std::uint64_t key) noexcept#

Stateless 64-bit mix of a key, used for counter-based random numbers.

Parameters:: key – input value (for example a counter or an index)
Returns:: one SplitMix64 output for the state key

constexpr double to_unit_interval(std::uint64_t bits) noexcept#

Uniform double in \([0, 1)\) from the 53 high bits of a 64-bit integer.

Parameters:: bits – raw 64-bit random value
Returns:: (bits >> 11) * 2^-53, exactly representable

template<int D> inline double dot(const Vec<D> &a, const Vec<D> &b) noexcept#

Euclidean inner product.

Template Parameters:: D – dimension
Returns:: \(\sum_k a_k b_k\)

template<int D> inline double norm2(const Vec<D> &a) noexcept#

Squared Euclidean norm.

Template Parameters:: D – dimension
Returns:: \(\sum_k a_k^2\)

Variables

constexpr double never = std::numeric_limits<double>::infinity()#: Event time meaning “no contact” (positive infinity).

template<int D> struct BallWall#

#include <box.hpp>

Curved wall: the region \(\sum_{k \in A} (x_k - c_k)^2 \le R^2\) over a set of axes \(A\).

Over two axes of a 3D box it is a cylinder, over all axes a spherical container (3D) or a disk (2D).

Template Parameters:: D – dimension (2 or 3)

Public Members

bool active = false#: false: no curved wall

std::array<bool, D> axes = {}#: axes the ball constrains

Vec<D> center = {}#: centre (only the components on axes are used)

double radius = 0.0#: radius \(R\)

template<int D> struct Box#

#include <box.hpp>

Axis-aligned box \([0, L_0) \times \dots \times [0, L_{D-1})\) with a boundary type per axis, and an optional curved wall inside it.

Along a periodic axis the box repeats. Along a non-periodic axis that the ball does not constrain, flat walls bound the box at 0 and \(L_k\). Along the axes of the ball, the ball wall bounds the spheres; it must lie inside the box.

Template Parameters:: D – dimension (2 or 3)

Public Functions

inline bool flat_walls(int k) const noexcept#

Returns:: true if axis k is bounded by flat walls at 0 and \(L_k\)

inline void validate() const#

Check the geometry.

Throws:: std::invalid_argument – if an edge length is not positive (or NaN), or if the ball wall is periodic along one of its axes or does not fit in the box

inline double volume() const noexcept#

Returns:: volume (D=3) or area (D=2) available to the spheres: the box, or the ball over its axes times the box edges along the other axes

inline double ball_distance2(const Vec<D> &x) const noexcept#: Squared distance from the ball centre, over the ball axes only.

Public Members

Vec<D> L = {}#: edge lengths; the box spans [0, L_k) along axis k

std::array<bool, D> periodic = {}#: per axis; false: walls (flat, or the ball)

BallWall<D> ball = {}#: optional curved wall

template<int D> class CellGrid#

#include <cell_grid.hpp>

Cell list with O(1) insertion, removal and move of particles.

Particles are stored in one doubly linked list per cell. Particle indices and cell indices are 32-bit; none marks the end of a list and a particle in no cell.

Template Parameters:: D – dimension (2 or 3)

Public Types

using Index = std::array<int, D>#: Integer cell coordinates, one per axis.

Public Functions

inline void build(const Box<D> &box, double min_width, std::size_t max_cells, int stencil = 1)#

Choose the cell layout for a box and empty all cells.

Cells are at least min_width wide along every axis; the total number of cells is capped near max_cells by widening cells (never by narrowing them). Call clear_particles afterwards, before inserting.

Parameters:

box – container (copied; its periodicity is used by for_each_neighbor)
min_width – lower bound on the cell width along every axis, typically the largest interaction distance; a value <= 0 requests the finest grid allowed
max_cells – approximate upper bound on the total number of cells
stencil – stencil radius m in cells: neighbours are searched in offsets -m..m per axis, so cells may be as narrow as the interaction distance divided by m

inline void clear_particles(std::size_t n_particles)#

Remove all particles and size the per-particle arrays.

Parameters:: n_particles – number of particles that will be inserted

inline double min_width() const noexcept#

Returns:: smallest cell width over all axes

inline int n(int k) const noexcept#

Returns:: number of cells along axis k

inline int stencil() const noexcept#

Returns:: the stencil radius m

inline double width(int k) const noexcept#

Returns:: cell width along axis k

inline std::size_t num_cells() const noexcept#

Returns:: total number of cells

inline Index coords_of(const Vec<D> &x) const noexcept#

Cell coordinates containing a position.

Parameters:: x – position; coordinates outside the box are clamped to the edge cells
Returns:: integer cell coordinates

inline std::int32_t flat(const Index &c) const noexcept#

Flatten cell coordinates (axis 0 varies fastest).

Parameters:: c – cell coordinates, within range
Returns:: cell index

inline Index coords(std::int32_t f) const noexcept#

Inverse of flat.

Parameters:: f – cell index
Returns:: cell coordinates

inline std::int32_t cell_of(std::uint32_t i) const noexcept#

Parameters:: i – particle index
Returns:: cell holding particle i, or none if it is not inserted

inline void insert(std::uint32_t i, std::int32_t cell) noexcept#

Put particle i at the head of the list of cell.

Parameters:

i – particle index; must not currently be in a cell
cell – destination cell index

inline void remove(std::uint32_t i) noexcept#

Take particle i out of its cell; afterwards cell_of(i) == none.

Parameters:: i – particle index; must currently be in a cell

inline void move(std::uint32_t i, std::int32_t cell) noexcept#

Move particle i to another cell (remove followed by insert).

Parameters:

i – particle index; must currently be in a cell
cell – destination cell index

template<class F> inline void for_each_neighbor_layer(std::int32_t cell, int axis, int dir, F &&f) const#

Like for_each_neighbor(), but only the layer of the stencil at offset dir (+1 or -1) along axis: the cells that become adjacent when a particle moves into cell across that face (offset m along axis).

Callers fall back to a full scan when there are fewer than 2m + 1 cells along axis.

Parameters:

cell – centre cell of the stencil
axis – axis of the move
dir – direction of the move along axis (+1 or -1)
f – callable f(std::uint32_t j, const std::array<int, D>& shift)

template<class F> inline void for_each_neighbor(std::int32_t cell, F &&f) const#

Visit every particle in the \((2m+1)^D\) stencil of cells around cell.

Calls f(j, shift) for each particle j found, where the image of j to use is \(x_j + \mathrm{shift}\cdot L\) (componentwise integers). Every stencil offset is a distinct periodic image. On periodic axes with fewer than 2m + 1 cells the stencil visits the same cell with different shifts, which is what makes small periodic boxes correct. On non-periodic axes stencil cells outside the grid are skipped. Includes the particle itself with shift 0.

Parameters:

cell – centre cell of the stencil
f – callable f(std::uint32_t j, const std::array<int, D>& shift); the grid must not be modified from inside f

Public Static Attributes

static constexpr std::int32_t none = -1#: Sentinel for “no particle” / “no cell”.

template<int D> class EventChainMC#

#include <rejection_free.hpp>

Straight event-chain variant of the rejection-free method (Peters and de With 2012, section “straight event-chain collision”; Bernard, Krauth and Wilson 2009 for hard spheres).

One particle moves at a time, along a coordinate axis. For every pair with the moving particle the uphill part of the pair potential along the path is accumulated; at a reflection the moving particle stops and the partner continues with the same displacement direction (a lift). A chain ends after a total displacement chain_length; the next chain starts at a random particle. In the irreversible mode (the default) the direction cycles through +x, +y, +z; this breaks detailed balance but satisfies global balance, so the canonical distribution is sampled. In the reversible mode each chain gets a random axis and sign.

Template Parameters:: D – dimension (2 or 3)

Public Functions

inline EventChainMC(const Box<D> &box, const std::vector<Vec<D>> &positions, PairPotential potential, double kT, std::uint64_t seed, bool irreversible = true)#

Set up particles at the given positions.

Parameters:

box – periodic box (all axes periodic), edges at least twice the cutoff
positions – initial positions
potential – pair potential (prepared here)
kT – temperature in energy units
seed – random seed
irreversible – cycle the directions +x, +y, +z (true) or draw random axes and signs

inline const EventChainStats &run(double displacement, double chain_length, const std::function<bool()> &interrupted = {})#

Run chains until the total displacement reaches displacement.

Parameters:

displacement – total displacement of all chains together
chain_length – displacement per chain
interrupted – polled regularly; returning true stops the run early

Returns:

the accumulated statistics

inline std::vector<Vec<D>> positions() const#

Returns:: current positions, wrapped into the box

inline double potential_energy() const#

Returns:: total potential energy

inline const EventChainStats &stats() const noexcept#

Returns:: the accumulated statistics

struct EventChainStats#

#include <rejection_free.hpp>

Run statistics of an event-chain simulation.

Public Members

double displacement = 0.0#: total displacement of all chains

std::uint64_t n_lifts = 0#: number of lifts (the motion passed to another particle)

std::uint64_t n_chains = 0#: number of chains

double wall_time = 0.0#: wall-clock seconds spent in run()

bool interrupted = false#: the last run() was stopped by the interrupt callback

class EventHeap#

#include <event_queue.hpp>

Indexed binary min-heap with one key per particle.

Particle i always owns slot i; its key can be changed in O(log n). The top is the particle with the smallest key (the next event). Keys equal to infinity mean “no event”.

Public Functions

inline void reset(std::uint32_t n)#

Resize to n particles and set every key to infinity.

Parameters:: n – number of particles

inline void assign(const std::vector<double> &keys)#

Set all keys at once (O(n) heapify).

Parameters:: keys – one key per particle; must have the size given to reset

inline void update(std::uint32_t i, double t)#

Change the key of particle i and restore the heap order, O(log n).

Parameters:

i – particle index
t – new key (event time)

inline std::uint32_t top() const noexcept#

Returns:: particle with the smallest key; undefined if empty()

inline double top_key() const noexcept#

Returns:: smallest key; undefined if empty()

inline double key(std::uint32_t i) const noexcept#

Parameters:: i – particle index
Returns:: current key of particle i

inline bool empty() const noexcept#

Returns:: true if the heap holds no particles

struct HistoryRecord#

#include <ls_packing.hpp>

One record per pressure window.

Public Members

double time#: simulation time at the end of the window

double scale#: s at the end of the window

double reduced_pressure#: mean over spheres

double median_pressure#: median over spheres

double kT#: time-averaged temperature over the window

std::uint64_t n_collisions#: collisions since the start of the run

int phase#: 0 growing, 1 relaxing (fixed radii), 2 growth step

template<int D> class LSPacking#

#include <ls_packing.hpp>

Event-driven Lubachevsky-Stillinger packing of polydisperse growing hard spheres.

The spheres move ballistically with unit mass at temperature kT = 1 and collide elastically while their radii \(a_i s(t)\) grow at a constant rate. The temperature is restored to 1 at every synchronisation while the spheres grow. Pressure is measured in windows of sync_interval collisions per particle as the reduced pressure \(Z = PV/(NkT) = 1 + \sum \mathbf r\cdot\mathbf J / (2 \int KE\, dt)\), from the virial of the collision impulses \(\mathbf J\) and the time integral of the kinetic energy. The container may have periodic axes, flat walls and a spherical wall (see Box).

Template Parameters:: D – dimension (2 or 3)

Public Functions

inline LSPacking(const Box<D> &box, std::vector<double> radii, const PackingOptions &opt)#

Creates a packing with random positions at zero radius and random velocities.

The growth rate is set so that the mean diameter grows at opt.growth_rate thermal speeds. The generator is seeded with opt.seed.

Parameters:

box – container; validated here
radii – relative radii a_i; the radii at scale s are a_i * s
opt – run parameters

Throws:

std::invalid_argument – if the box is invalid, radii is empty or has too many entries, a radius is not positive and finite, growth_rate < 0 or target_scale <= 0

inline void randomize_positions()#: Places the particles uniformly at random (inside the ball wall, if any) and sets the scale to s = 0, so that the radii are zero.

inline void randomize_velocities()#: Draws Maxwellian velocities with zero total momentum, scaled to kT = 1 exactly.

inline void set_positions(const std::vector<Vec<D>> &x, double scale)#

Starts from given positions at scale s (radii a_i * s); they must not overlap.

Resets the periodic image counters. The positions are not checked for overlaps.

Parameters:

x – positions, inside the box
scale – initial scale s, at least 0

Throws:

std::invalid_argument – if the number of positions differs from the number of particles or scale < 0

inline void set_velocities(const std::vector<Vec<D>> &v)#

Sets the velocities, replacing the random ones.

The temperature is rescaled to kT = 1 at the first synchronisation of a growing run.

Parameters:: v – velocities, one per particle
Throws:: std::invalid_argument – if the number of velocities differs from the number of particles

inline const PackingStats &run(const std::function<bool()> &interrupted = {})#

Runs until a stopping criterion is met.

The criteria are the target scale, the pressure limits, the jamming protocol, stall detection, the collision, simulation time and wall-clock limits, the interrupt callback, the absence of events and the box limit; see Status. On return the radii are shrunk uniformly if round-off left overlaps (PackingStats::shrink_factor).

Parameters:

interrupted – optional callback polled regularly; the run stops with Status::interrupted when it returns true

Throws:

std::invalid_argument – if the growth rate is zero and the initial scale is zero
std::runtime_error – if the box is smaller than the largest particle diameter

Returns:

the statistics of the run, valid until the next call of run()

inline std::size_t size() const noexcept#

Returns:: the number of particles

inline const std::vector<HistoryRecord> &history() const noexcept#

Returns:: one record per pressure window of the last run

inline std::vector<double> sphere_pressure() const#

Per-sphere reduced pressure of the last completed window: 1 + N w_i / (2 int KE dt).

Without walls its mean over spheres is the reduced pressure; wall contacts add to the spheres that touch a wall. Close to 1 for rattlers.

Returns:: one value per sphere; empty before the first window has closed

inline void set_full_rescan(bool on) noexcept#

Testing aid: re-predict with a full neighbour scan after every cell crossing instead of scanning only the new cell layer.

The trajectory must be the same.

inline const PackingStats &stats() const noexcept#

Returns:: the statistics of the last run (status Status::running before the first run)

inline const Box<D> &box() const noexcept#

Returns:: the container

inline std::vector<Vec<D>> positions() const#

Positions wrapped into the box; valid after run().

Returns:: one position per particle

inline std::vector<Vec<D>> velocities() const#

Returns:: one velocity per particle

inline std::vector<double> radii() const#

Final radii a_i * stats().scale; valid after run().

Returns:: one radius per particle

struct PackingOptions#

#include <ls_packing.hpp>

Parameters of a run.

Every stopping criterion is off when set to its default.

Lengths are in the units of the box; never is positive infinity (see predict.hpp).

Jamming protocol

Off when jammed_pressure is infinite: grow until the median reduced pressure reaches jam_start_pressure, then alternate relaxation at fixed radii (jam_relax_windows pressure windows) with growth steps of relative size jam_step / Z_median, until Z_median after relaxation reaches jammed_pressure.

double jammed_pressure = never #: median reduced pressure at which the packing counts as jammed

double jam_start_pressure = 1e3#: median reduced pressure at which the protocol starts

int jam_relax_windows = 4#: pressure windows per relaxation at fixed radii

double jam_step = 0.2#: fraction of the remaining relative growth closed per step

Public Members

double target_scale = 1.0#: Final radii are a_i * target_scale; infinity: no target.

double growth_rate = 0.02#: Growth speed of the mean diameter in units of the thermal speed sqrt(kT/m).

double max_pressure = never #: Reduced pressure PV/(NkT) at which growth stops (mean over a measurement window).

double stall_tol = 0.0#: 0: off. Stop when s grew by less than this (relative)

double stall_window = 100.0#: over this many collisions per particle

std::uint64_t max_collisions = 0#: total collisions at which to stop; 0: unlimited

CollisionRule rule = CollisionRule::legacy #: collision rule

double max_time = never #: simulation time

double timeout = never #: wall-clock seconds

double sync_interval = 10.0#: collisions per particle between synchronisations

double separation_margin = 1e-9#: After a collision the pair separates at least this much faster than the contact distance grows (thermal speed units); keeps the gap opening in floating point.

std::uint64_t seed = 0#: seed of the random number generator

struct PackingStats#

#include <ls_packing.hpp>

Summary of a finished run, returned by LSPacking::run().

Public Members

Status status = Status::running #: why the run ended

double scale = 0.0#: final s: radii are a_i * scale

double density = 0.0#: volume (area) fraction

double reduced_pressure = 0.0#: last measured PV/(NkT) (mean over spheres); 0 if never measured

double median_pressure = 0.0#: median over spheres of the per-sphere reduced pressure

std::uint64_t n_collisions = 0#: number of collisions that changed velocities

std::uint64_t n_events = 0#: number of processed events, including cell crossings

double sim_time = 0.0#: simulation time

double wall_time = 0.0#: wall-clock seconds

double min_gap_ratio = 0.0#: min over contacts of distance / contact distance - 1, before shrink

double shrink_factor = 1.0#: uniform radius factor applied to remove round-off overlaps

double max_contact_error = 0.0#: max over collisions of |distance / contact distance - 1|

struct PairPotential#

#include <rejection_free.hpp>

Isotropic pair potential with at most one minimum, zero at and beyond the cutoff.

U decreases on \((0, r_m]\) and increases on \([r_m, r_c]\) to \(U(r_c) = 0\). Purely repulsive potentials have \(r_m = r_c\). Hard spheres have an infinite step at \(\sigma\).

Public Types

enum class Kind#

Potential families.

Values:

enumerator soft#: \(U = \tfrac{\epsilon}{\alpha}(1 - r/\sigma)^\alpha\) for \(r < \sigma\) (DPD: \(\alpha = 2\))

enumerator lennard_jones#: \(4\epsilon[(\sigma/r)^{12} - (\sigma/r)^6] - U(r_c)\) for \(r < r_c\)

enumerator hard#: infinite for \(r < \sigma\), zero otherwise

Public Functions

inline void prepare()#: Check parameters and derive the shift and the minimum.

inline double r_min() const noexcept#

Returns:: radius of the minimum \(r_m\)

inline double u_min() const noexcept#

Returns:: \(U(r_m)\), the lowest value (0 for purely repulsive potentials)

inline double value(double r) const noexcept#

Returns:: \(U(r)\)

inline double inner_radius(double u) const noexcept#: Radius on the repulsive branch \(r \le r_m\) where \(U(r) = u\), for \(u \ge U(r_m)\).

inline double outer_radius(double u) const noexcept#: Radius on the attractive branch \(r_m \le r \le r_c\) where \(U(r) = u\), for \(U(r_m) \le u \le 0\).

Public Members

Kind kind = Kind::soft #: potential family

double epsilon = 1.0#: energy scale

double sigma = 1.0#: length scale (soft and hard: range / diameter)

double cutoff = 1.0#: \(r_c\) (soft and hard: equal to sigma)

double alpha = 2.0#: exponent of the soft potential

template<int D> class RejectionFreeMC#

#include <rejection_free.hpp>

Rejection-free event-driven Monte Carlo in a periodic box.

Template Parameters:: D – dimension (2 or 3)

Public Functions

inline RejectionFreeMC(const Box<D> &box, const std::vector<Vec<D>> &positions, PairPotential potential, double kT, std::uint64_t seed)#

Set up particles at the given positions with Gaussian move velocities.

Parameters:

box – periodic box (all axes periodic)
positions – initial positions inside the box
potential – pair potential (prepared here)
kT – temperature in energy units
seed – random seed (velocities and reflection energies)

inline void set_velocities(const std::vector<Vec<D>> &v)#: Replace the move velocities (for example to redraw them).

inline void redraw_velocities()#: Redraw all move velocities from a standard normal distribution.

inline const RejectionFreeStats &run(double duration, const std::function<bool()> &interrupted = {})#

Advance the simulation by duration (time = contour length along the moves).

Parameters:

duration – simulation time to advance
interrupted – polled regularly; returning true stops the run early

Returns:

the accumulated statistics

inline double time() const noexcept#

Returns:: simulation time since the start

inline std::vector<Vec<D>> positions()#

Returns:: positions at the current time, wrapped into the box

inline std::vector<Vec<D>> velocities() const#

Returns:: current move velocities

inline double potential_energy()#

Returns:: total potential energy at the current time

inline void set_full_rescan(bool on) noexcept#

Testing aid: re-predict with a full neighbour scan after every cell crossing instead of scanning only the new cell layer.

The trajectory must be the same.

inline const RejectionFreeStats &stats() const noexcept#

Returns:: the accumulated run statistics

inline std::size_t size() const noexcept#

Returns:: the number of particles

struct RejectionFreeStats#

#include <rejection_free.hpp>

Run statistics of a rejection-free simulation.

Public Members

double time = 0.0#: simulation time (“contour length” of the moves)

std::uint64_t n_reflections = 0#: number of pair reflections

std::uint64_t n_events = 0#: processed events (reflections, cell crossings, stale predictions)

double wall_time = 0.0#: wall-clock seconds spent in run()

bool interrupted = false#: the last run() was stopped by the interrupt callback

class Rng#

#include <rng.hpp>

xoshiro256++ pseudo-random generator (Blackman and Vigna, 2019).

Satisfies the interface of a C++ uniform random bit generator (result_type, min, max, operator()), and provides portable uniform, exponential and normal variates. The state is 256 bits, seeded through SplitMix64.

Public Types

using result_type = std::uint64_t#: Type of the raw integer output.

Public Functions

inline explicit Rng(std::uint64_t seed = 0) noexcept#

Construct and seed the generator.

Parameters:: seed – any 64-bit seed; the four state words are drawn from SplitMix64

inline void reseed(std::uint64_t seed) noexcept#

Reset the stream to the one defined by seed and drop any cached normal variate.

Parameters:: seed – any 64-bit seed

inline result_type operator()() noexcept#

Next raw 64-bit value of the xoshiro256++ stream.

Returns:: uniformly distributed 64-bit integer

inline double uniform() noexcept#

Returns:: uniform variate in \([0, 1)\)

inline double uniform_open0() noexcept#

Returns:: uniform variate in \((0, 1]\), safe as argument of log

inline double exponential() noexcept#

Returns:: standard exponential variate (rate 1), \(-\ln u\) with \(u \in (0,1]\)

inline double normal() noexcept#

Standard normal variate by the Box-Muller transform.

Each transform yields two independent values; the second is cached and returned by the next call.

Returns:: normal variate with zero mean and unit variance

inline std::array<std::uint64_t, 4> state() const noexcept#

Returns:: copy of the four 64-bit state words (the cached normal is not included)

Public Static Functions

static inline constexpr result_type min() noexcept#: Smallest raw output (0).

static inline constexpr result_type max() noexcept#: Largest raw output (all bits set).