API Reference

Core module: fundamental types & utilities.

Data representation (current): positions :: Vector{SVector{D,T}} velocities :: Vector{SVector{D,T}} forces(R) :: Vector{SVector{D,T}} (or (F,U) with return_potential=true)

Matrix-based NxD layouts are deprecated.

struct CubicBox{T<:Real}

Simple cubic periodic box with side length L.

ParticleSystem{Dims,T_float}

A lightweight container for particle states.

Fields

• positions::Vector{SVector{Dims,T_float}}: length-N vector of positions (each an SVector)
• velocities::Vector{SVector{Dims,T_float}}: length-N vector of velocities
• masses::Vector{T_float}: length-N vector of particle masses

kinetic_energy(system::ParticleSystem) -> T_float

Total kinetic energy: ∑ ½ mᵢ ‖vᵢ‖².

minimum_image(Δ, box::CubicBox)
minimum_image(Δ, L)

Apply the minimum-image convention to the displacement vector Δ for the periodic box (wraps components to (-L/2, L/2]) and return the new vector.

potential_energy(system::ParticleSystem{Dims,T_float}, forces::Function) -> T_float

Try to obtain total potential energy using forces(...; return_potential=true) convention.

velocity_verlet!(system::ParticleSystem{Dims,T_float}, forces::Function, dt::T_float)

Advance system by one timestep dt using the velocity Verlet integrator.

wrap_positions!(R, box::CubicBox)

Wrap particle positions R into the primary periodic image of box in-place. This is useful before building neighbor lists or measuring displacements.

The resulting positions will be in the range (-L/2, L/2].

CellGrid{IT<:Integer, T<:Real}

Linked-list cell grid for cubic periodic boxes. The domain is split into a uniform dims = (nx,ny,nz) grid of cubic cells with an intrusive linked list (heads, next) storing particle indices in each cell. This enables O(N) rebinning at fixed density and powers the O(N) neighbor build.

Fields

• L::T: cubic box length.
• cell_size::T: effective uniform cell width actually used for binning.
• dims::NTuple{3,IT}: number of cells along each axis (each ≥ 1).
• heads::Vector{IT}: length nx*ny*nz; head index per cell (0 sentinel = empty).
• next::Vector{IT}: length N; linked list “next” pointer per particle (0 = end).

NeighborList{IT<:Integer,T<:Real}

Compressed-sparse-row (CSR) neighbor list with snapshot of reference positions.

Fields

• cutoff::T: physical cutoff distance (rcut).
• skin::T: extra buffer distance beyond cutoff.
• pairs::Vector{IT}: concatenated neighbor indices (CSR-style).
• offsets::Vector{IT}: length N+1; neighbors of particle i are stored in pairs[offsets[i]:(offsets[i+1]-1)].
• ref_positions::Matrix{T}: snapshot of positions at last rebuild, used for displacement checks.

Use build_neighborlist to construct, and maybe_rebuild! to keep it updated.

build_cellgrid(R, box; cell_size)

Create a new CellGrid sized for cell_size and bin positions R (N×3). Returns a populated grid with linked lists set for particle indices 1..N.

build_neighborlist(R, box; cutoff, skin=0.3) -> NeighborList

Construct a neighbor list for positions R inside a CubicBox. Pairs are included if they are within cutoff + skin under minimum-image distance.

This is an O(N²) operation but performed infrequently. The returned NeighborList uses a CSR (compressed sparse row) layout for compact storage and fast iteration.

build_neighborlist_cells(R, box; cutoff, skin=0.3, grid=nothing) -> NeighborList-like

Construct a half neighbor list (each pair appears exactly once with j>i) using an O(N) cell-linked grid build. Returns an object with .offsets and .pairs fields compatible with typical CSR traversal.

lj_forces(R, box, nlist; ϵ=1.0, σ=1.0, rcut=NaN, shift=false, return_potential=false)

Compute Lennard–Jones forces using a prebuilt NeighborList. This is more efficient than the brute-force O(N²) kernel, scaling ~O(N) for fixed density.

Arguments

• R::Vector{SVector{Dims,T}}: positions.
• box::CubicBox: periodic cubic box.
• nlist::NeighborList: neighbor list built with at least the requested cutoff.

Keywords

• ϵ: LJ well depth (default 1.0).
• σ: LJ length scale (default 1.0).
• rcut: optional cutoff to use for the force calculation (defaults to nlist.cutoff).
• shift::Bool: if true, subtract potential at cutoff to make it continuous.
• return_potential::Bool: if true, also return total potential energy.

Returns

• F::Vector{SVector{Dims,T}}: force on each particle (same shape as R).
• U::Float64 (if return_potential=true).

Example

box = CubicBox(8.0)
R = randn(SVector{3, Float64}, 10)
wrap_positions!(R, box)
nlist = build_neighborlist(R, box; cutoff=2.5, skin=0.4)
F, U = lj_forces(R, box, nlist; rcut=2.5, return_potential=true)

maybe_rebuild!(nlist, R, box) -> Bool

Check whether the neighbor list nlist is still valid given new positions R. If the maximum displacement since the last build is greater than skin/2, rebuild the list and return true. Otherwise leave it unchanged and return false.

rebin!(grid, R, box) -> grid

Reset the grid's heads and next and bin the positions R into cells according to the current cell_size and periodic cubic box box.

struct DistanceConstraints{T_int,T_float}

Immutable set of pairwise distance constraints for SHAKE/RATTLE.

apply_rattle!(ps, cons)

Correct velocities to satisfy velocity constraints.

apply_shake!(ps, cons, dt)

Iteratively correct positions to satisfy constraints.

constraint_residuals(ps::ParticleSystem, cons::DistanceConstraints) -> (; maxC, rmsC, maxCd, rmsCd)

Compute the constraint residuals for a system subject to pairwise distance constraints.

For each constraint l with atoms (i, j) and target distance r0_l:

• Position residual: C_l = ||r_i - r_j||^2 - r0_l^2
• Velocity residual: Ċ_l = 2 (r_i - r_j) ⋅ (v_i - v_j)

Returns a named tuple with maxC, rmsC, maxCd, rmsCd.

Example

ps = ParticleSystem([SVector(0.0, 0, 0), SVector(1.0, 0, 0)], [SVector(0.0, 0, 0), SVector(0.0, 0, 0)], [1.0, 1.0])
cons = DistanceConstraints([(1,2)], [1.0])
constraint_residuals(ps, cons)

remove_com_motion!(ps; which=:velocity)

Remove center-of-mass motion from velocities/positions.

velocity_verlet_shake_rattle!(ps, forces, dt, cons)

Constrained velocity Verlet step with SHAKE/RATTLE.

• degreesoffreedom(ps) -> Int

Return the effective translational degrees of freedom of the system. Currently this is simply N * D for N particles in D spatial dimensions.

This function is a hook for future extensions (e.g., constraints, rigid bodies, or COM removal) that reduce the effective DoF used in temperature calculations.

degrees_of_freedom(ps; constraints=nothing, remove_com=false) -> Int

Effective translational DoF, reduced by number of constraints and optionally COM removal. This is the canonical method. A no-keyword fallback is provided for backward compatibility.

• instantaneous_temperature(ps; kB=1.0) -> Float64

Compute the instantaneous scalar temperature via equipartition:

KE = 0.5 * Σᵢ mᵢ ||vᵢ||² T = 2 * KE / (kB * dof)

where dof = degrees_of_freedom(ps). For now, dof = N * D (translational DoF only).

Keywords

• kB::Real = 1.0: Boltzmann constant in your unit system.

Returns

• Float64: instantaneous temperature.

Notes

• For small N, expect large fluctuations. Compare time averages to the target T when thermostatting.

instantaneous_temperature(ps; kB=1.0) -> Float64

Compute instantaneous temperature via equipartition: T = 2 * KE / (kB * dof).

• langevinbaoab!(ps::ParticleSystem, forces, dt; γ, T, kB=1.0, rng=Random.defaultrng())

Advance the system by one NVT step using the Langevin BAOAB integrator.

This scheme integrates the Langevin SDE

m dv = [F(r) - γ m v] dt + √(2γ m kB T) dW, dr = v dt

via the B → A → O → A → B splitting with an exact Ornstein–Uhlenbeck (O) step. It is widely used for robust thermostatting and near-optimal configurational sampling. When γ == 0, the O-step is the identity and BAOAB reduces to velocity–Verlet (up to roundoff).

Arguments

• ps::ParticleSystem: System with positions::Vector{SVector{Dims,T}}, velocities::Vector{SVector{Dims,T}}, and masses::Vector{T}.
• forces: A callable F(R)::Vector{SVector{Dims,T}} returning forces evaluated at positions R.
• dt::Real: Time step (same time units as used for γ).

Keywords

• γ::Real: Friction coefficient [1/time].
• T::Real: Target temperature. Must be consistent with kB (so that kB*T has units of energy).
• kB::Real = 1.0: Boltzmann constant in your unit system.
• rng::AbstractRNG = Random.default_rng(): RNG used only in the stochastic O-step for reproducibility.

Details

The O-step uses the exact OU update: c = exp(-γ*dt) sᵢ = sqrt((1 - c^2) * kB*T / mᵢ) # per-particle, per-component scale vᵢ[k] ← c * vᵢ[k] + sᵢ * ξ, where ξ ~ 𝒩(0,1) For numerical stability when γ*dt ≪ 1, implementations should use the linearized forms c ≈ 1 - γ*dt and 1 - c^2 ≈ 2γ*dt to avoid catastrophic cancellation.

Performance tips

• Avoid allocations in the inner loop. Reuse force buffers if you add multi-step drivers.
• Use @views for in-place operations on ps.positions/ps.velocities.
• Pass an explicit rng to ensure reproducibility across runs and Julia sessions.

Pitfalls

• Units matter: ensure γ is in 1/time, and kB*T is in energy, consistent with r, v, and m.
• Very large γ*dt yields overdamped dynamics (OK for sampling, not for kinetics).
• With heterogeneous masses, compute sᵢ per particle; do not share a single scale across species.
• For tiny systems, instantaneous temperature fluctuates strongly — average over time for diagnostics.

Example

```julia

Pseudocode showing typical usage (assumes your package exports ParticleSystem and BAOAB):

using Random N, D = 64, 3 ps = ParticleSystem(zeros(N, D), randn(N, D), ones(N)) dt, γ, T = 0.005, 1.0, 1.5 rng = MersenneTwister(123) forces(R) = -0.1 .* R # harmonic "bath"

for _ in 1:2_000

• langevin_baoab!(ps, forces, dt; γ=γ, T=T, kB=1.0, rng=rng) end ```

langevin_baoab!(ps, forces, dt; γ, T, kB=1.0, rng=Random.default_rng())

Advance one step with the Langevin BAOAB integrator: B(half) → A(half) → O(OU) → A(half) → B(half). With γ=0, this reduces to velocity-Verlet (deterministic).

Arguments

• ps : particle system (positions N×D, velocities N×D, masses length N)
• forces : callable F = forces(positions)::N×D
• dt : time step

Keywords

• γ : friction (1/time)
• T : target temperature
• kB : Boltzmann constant (default 1.0)
• rng: AbstractRNG for reproducibility

langevin_baoab_constrained!(ps::ParticleSystem, forces, dt, cons::DistanceConstraints;
                            γ, T, kB=1.0, rng=Random.default_rng())

Advance one constrained NVT BAOAB step with SHAKE/RATTLE projections.

Splitting and projections:

1. B (half kick) → apply_rattle! (velocities)
2. A (half drift) → apply_shake! (positions)
3. O (OU stochastic) → apply_rattle! (velocities)
4. A (half drift) → apply_shake! (positions)
5. Recompute forces
6. B (half kick) → apply_rattle! (velocities)

Notes:

• Uses exp(-γ*dt) for the OU decay. Noise variance is (1 - c^2) * kB*T / m_i per component.
• With constraints, prefer this method over the unconstrained langevin_baoab!.

• velocity_rescale!(ps, T; kB=1.0)

Deterministically rescale all velocities to match a target temperature T. This is a one-shot utility often used for quick pre-equilibration or initializing an NVT run. It is not a thermostat and does not generate the correct kinetic energy distribution by itself.

Keywords

• kB::Real = 1.0: Boltzmann constant.

Notes

Let T₀ = instantaneous_temperature(ps; kB) and λ = √(T / max(T₀, eps())). This function applies vᵢ ← λ * vᵢ in-place.

velocity_rescale!(ps, T; kB=1.0)

Deterministically rescale velocities to match target temperature T.