Ensembles

An ensemble in mcpy is the object that owns the Monte Carlo loop. It holds the atomic configuration, the energy calculator, the set of cells used for free-volume estimation, and the MoveSelector that proposes trial configurations. It then evaluates the appropriate acceptance rule and writes the trajectory and log.

Four ensembles are available:

• CanonicalEnsemble – fixed composition (NVT) Metropolis sampling.

• GrandCanonicalEnsemble – variable composition at fixed \((\mu, V, T)\).

• ReplicaExchange – an MPI wrapper that runs several GCMC replicas at different temperatures or chemical potentials and periodically attempts exchanges between them.

• BatchedReplicaExchange – single-process variant that batches every replica’s trial-move energy into a single forward pass (one GPU context, one kernel launch per layer).

All ensembles share a basin-hopping-style running minimum output: whenever a strictly lower-score configuration is observed, it is appended to (or overwrites) a dedicated minima_file independent of the per-step sampling trajectory.

The remainder of this page describes each one, the acceptance criteria they implement, the free-volume estimator used by GCMC insertion/deletion, and the running-minimum output.

BaseEnsemble

Shared infrastructure used by every concrete ensemble:

• storage of atoms, cells, and the ASE-compatible calculator,

• single-point energy evaluation (optionally preceded by a local relaxation),

• trajectory and log writers driven by trajectory_write_interval and outfile_write_interval,

• initialization, finalization, and step-counter management,

• a running-minimum hook (_record_minimum) that snapshots the best-so-far configuration and optionally writes it to minima_file – see Running minima and basin-hopping output.

BaseEnsemble is abstract – you instantiate one of the concrete subclasses below.

CanonicalEnsemble

Samples the canonical (NVT) distribution at fixed composition. Trial moves are proposed by the configured MoveSelector and accepted with the Metropolis criterion.

For a configuration \(i\) with energy \(E_i\), the canonical Boltzmann weight is

\[P_i = \frac{e^{-\beta E_i}}{Z}, \qquad Z = \sum_j e^{-\beta E_j},\]

with \(\beta = 1/(k_B T)\). The acceptance probability for a trial move \(i \rightarrow j\) reduces to the symmetric Metropolis form

\[P_{ij}^{\mathrm{acc}} = \min\!\left(1,\; e^{-\beta (E_j - E_i)}\right).\]

Typical use: structural optimization or thermal sampling at fixed stoichiometry.

GrandCanonicalEnsemble

Samples adsorption/desorption equilibria where the number of selected species fluctuates in contact with a reservoir at fixed chemical potentials.

In a grand-canonical run, the user fixes:

• temperature – sets \(\beta = 1/(k_B T)\),

• mu – a dictionary mapping each chemical species to its reservoir chemical potential,

• species – the chemical symbols that the move set is allowed to insert or delete,

• cells – one or more cell objects defining the insertion region and providing the accessible (free) volume,

• move_selector – the weighted collection of trial moves.

Acceptance rules

For a number-conserving move (e.g. displacement, permutation), GrandCanonicalEnsemble falls back to the Metropolis criterion above. For number-changing moves, the de-Broglie factors appear:

Deletion \((N \rightarrow N-1)\):

\[P_{ij}^{\,N \rightarrow N-1} = \min\!\left(1,\; \frac{N\,\Lambda^3}{z\,V_{\mathrm{free}}} \,e^{-\beta (E_j - E_i)}\right).\]

Insertion \((N \rightarrow N+1)\):

\[P_{ij}^{\,N \rightarrow N+1} = \min\!\left(1,\; \frac{z\,V_{\mathrm{free}}}{(N+1)\,\Lambda^3} \,e^{-\beta (E_j - E_i)}\right).\]

Here \(z = e^{\beta \mu}\) is the species activity, and

\[\Lambda = \frac{h}{\sqrt{2\pi m\, k_B T}}\]

is the thermal de Broglie wavelength. The geometric volume \(V\) of the textbook formulation is replaced by the accessible free volume \(V_{\mathrm{free}}\) returned by the configured cell – see Free volume and species_radii.

from mcpy.ensembles.grand_canonical_ensemble import GrandCanonicalEnsemble
from mcpy.moves import MoveSelector

gcmc = GrandCanonicalEnsemble(
    atoms=atoms,
    cells=[cell],
    calculator=calculator,
    mu={"O": -5.0, "Ag": -2.9},
    units_type="metal",
    species=["O"],
    temperature=500,
    move_selector=move_selector,
)
gcmc.run(steps=10000)

Outputs

• A trajectory is appended to traj_file (extended XYZ, with energy= and Lattice= on the comment line) every trajectory_write_interval steps.

• A human-readable log is written to outfile every outfile_write_interval steps, containing the step index, particle count, current energy, and the per-interval acceptance ratio of each registered move. The interval counters are reset after every write; cumulative ratios are reported at the end of the run.

Setting a unit system

units_type selects the constants the acceptance rules use. Pass 'metal' for eV, Å, and K, or 'LJ' for reduced Lennard-Jones units. The choice must match the units of the energies and lengths the calculator returns. A MACE or DFT calculator works in metal units; a reduced Lennard-Jones potential works in LJ units.

In 'metal' units mcpy uses physical constants and reads each species’ mass from ASE. The inverse temperature is

\[\beta = \frac{1}{k_B T}, \qquad k_B = 8.617333\times10^{-5}\ \mathrm{eV/K},\]

and the thermal de Broglie wavelength \(\Lambda\) of each species follows from its mass and the temperature through the expression given in the acceptance rules above, evaluated in Å.

In 'LJ' units every constant equals one. The Boltzmann and Planck constants are 1, all masses are 1, and \(\Lambda = 1\) for every species, so \(\beta = 1/T^*\) with \(T^*\) the reduced temperature. The de Broglie factors then drop out of the insertion and deletion rules, leaving their reduced-unit forms.

These constants live in the SetUnits helper, which each ensemble builds from units_type, temperature, and species. It computes \(\beta\) and the per-species \(\Lambda\), and exposes the two factors, de_broglie_insertion and de_broglie_deletion, that enter the number-changing acceptance rules above. See Utilities for its signature.

Running minima and basin-hopping output

When the trial move is followed by a local relaxation, the Metropolis loop is, by construction, a basin-hopping sampler: each accepted move corresponds to a local minimum of the potential energy surface. Two output channels coexist:

• traj_file (extended XYZ) – the sampling trajectory, written every trajectory_write_interval accepted steps. Use this for thermodynamic averages and visualisation.

• minima_file (extended XYZ, off by default) – the basin-hopping output. A frame is written only when the score of the just-accepted configuration is strictly lower than every previous score observed in the run.

Either channel can be disabled independently by passing None. The four useful combinations:

traj_file	minima_file	What you get
'traj.xyz'	None	Standard sampling (default).
None	'minima.xyz'	Pure basin hopping – only improving minima are persisted.
'traj.xyz'	'minima.xyz'	Both sampling and the basin-hopping history.
None	None	No XYZ output (energies still go to outfile).

minima_mode controls the policy of the minima file:

• 'a' (default) – append every new minimum, building a history of improving structures. The last frame of the file is always the running best.

• 'w' – overwrite, keeping a single frame with the current best configuration.

Each minima frame carries step= and score= in its comment line, so the improvement trace can be reconstructed from the file alone.

Score function

The quantity compared across configurations is the score, returned by _minimum_score. Defaults:

• CanonicalEnsemble – the potential energy \(E\).

• GrandCanonicalEnsemble – the grand potential

  \[\Omega = E - \sum_i \mu_i N_i,\]

  which is the meaningful comparison across moves that change particle counts. Subclasses may override _minimum_score to introduce custom criteria (e.g. order parameters).

Global minimum across replicas

Both ReplicaExchange and BatchedReplicaExchange accept a global_minimum_file argument. After all replicas finalise, the lowest-score configuration observed across the whole ladder is written as a single XYZ frame to that file. The MPI variant gathers per-rank bests to rank 0; the batched variant minimises over the in-process replicas. The per-replica minima_file history is independent of this and continues to live in the per-rank outputs (e.g. minima_rank_3.xyz).

Free volume and species_radii

In the textbook GCMC equations, \(V\) is the geometric volume of the simulation box. In practice, dense atomistic systems leave only a fraction of that volume actually accessible to a new particle. mcpy follows the hybrid scheme of Senftle et al.: every cell object estimates an accessible volume \(V_{\mathrm{free}}\) by Monte Carlo sampling and that quantity is used in the acceptance rules above.

Let \(V_{\mathrm{cell}}\) be the geometric volume of the insertion region and \(N_{\mathrm{MC}}\) the number of random sample points \(\mathbf{x}_k\). Each sample point is classified by an indicator that tests whether it falls inside the exclusion sphere of any atom \(a\):

\[\begin{split}I_k = \begin{cases} 1, & \exists\, a\ \text{such that}\ \|\mathbf{x}_k-\mathbf{r}_a\|^2 \le r_{\mathrm{species}(a)}^2,\\ 0, & \text{otherwise.} \end{cases}\end{split}\]

The occupied fraction and the resulting free volume are

\[f_{\mathrm{occ}} = \frac{1}{N_{\mathrm{MC}}}\sum_{k=1}^{N_{\mathrm{MC}}} I_k, \qquad V_{\mathrm{free}} = V_{\mathrm{cell}}\,(1 - f_{\mathrm{occ}}).\]

Two parameters control this estimator on every cell:

• species_radii – an element-wise mapping from chemical species to an exclusion radius. Physically meaningful values are critical: if the radii are too small, GCMC will spend most of its time proposing overlapping insertions; if they are too large, the accessible volume collapses and insertions become impossible. See Calibrating species_radii for a reproducible workflow to calibrate them from short relaxation trials.

• mc_sample_points – number of random points used in the estimate. More points reduce statistical noise on \(V_{\mathrm{free}}\) at the cost of a single up-front sweep per configuration change.

With or without relaxation

GCMC workflows in mcpy are typically run in one of two styles:

• No relaxation – propose, evaluate the trial energy, accept/reject. Fast per step, but acceptance is dominated by overlap rejections in dense regions.

• With relaxation – propose, run a short local optimization on the trial structure, then accept/reject using the relaxed energy. This is the workflow used in the reference application: relaxation absorbs unphysical overlaps and lets the system reach stable acceptance behaviour in roughly an order of magnitude fewer steps.

In both cases the free-volume estimate enters the insertion/deletion acceptance, so good species_radii calibration pays off either way.

ReplicaExchange

ReplicaExchange parallelises sampling by running several GCMC replicas at different temperatures or different chemical potentials, one per MPI rank, and periodically attempting to swap their configurations. It requires mpi4py and a user-supplied factory that constructs the rank-local GrandCanonicalEnsemble.

The exchange between two neighbouring replicas in temperature ladders is accepted with the standard parallel-tempering criterion

\[P_{\mathrm{exch}} = \min\!\left(1,\; e^{(\beta_2 - \beta_1)(E_2 - E_1)}\right),\]

and an analogous expression involving \(\beta\mu \Delta N\) is used when replicas differ in chemical potentials. After acceptance, replicas swap their configurations and energies; partner selection alternates between even and odd ranks on consecutive exchange attempts to ensure every replica eventually mixes with both neighbours.

Notes:

• the ensemble factory must produce per-rank output filenames (e.g. by including rank in the trajectory, log, and minima paths) so that ranks do not race on the same files;

• one MPI rank corresponds to one replica;

• at the end of the run, rank 0 collects each replica’s running best (via comm.gather) and writes the global minimum as a single XYZ frame to global_minimum_file (default "global_minimum.xyz"; pass None to disable). See Running minima and basin-hopping output.

BatchedReplicaExchange

BatchedReplicaExchange runs every replica in a single Python process and evaluates the trial-move energy of all replicas in one batched calculator call. The accept/reject step is still sequential per replica – the batching is over replicas at the same logical step. This is the right choice on a single GPU, where MPI-based replica exchange would either oversubscribe one CUDA context or serialise through it.

Requirements:

• the calculator must implement get_potential_energies(atoms_list) -> ndarray – the optional AlchemiCalculator does;

• gcmc_factory(T=..., rank=i) (or mu=..., rank=i for chemical-potential ladders) must return a fresh, independently-seeded GrandCanonicalEnsemble with its own cells and move selector.

Output layout mirrors ReplicaExchange: per-replica outfile/traj_file/ minima_file (typically tagged with rank) plus a single global_minimum_file written from min(replicas, key=best_score) after the run completes.

Choosing the right ensemble

• CanonicalEnsemble – fixed-composition sampling, structural relaxation, or thermal annealing.

• GrandCanonicalEnsemble – variable-composition adsorption/desorption studies at fixed \((\mu, T)\), e.g. oxidation phase diagrams.

• ReplicaExchange – multi-node MPI runs, when a single GCMC trajectory is prone to getting trapped in a single basin and broader thermodynamic coverage is needed.

• BatchedReplicaExchange – single-GPU runs, where batching trial energies across replicas is faster than launching one CUDA context per replica.