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. Every ensemble stores the atoms, cells, and calculator, writes the trajectory and log at the configured intervals, and tracks the running minimum (see Running minima and basin-hopping output).

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).\]

temperature=0 gives \(\beta = \infty\), a greedy descent that accepts only moves that do not raise the energy. Each trial is relaxed with the configured optimizer for at most steps steps (default 1000); a warning is logged when a relaxation stops at the cap before reaching fmax.

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 – the cells whose free volume the ensemble refreshes after every accepted move. The ensemble adds the cell of every insertion move to this list automatically, because the insertion acceptance reads the cached volume of its cell. You can still pass other cells here,
move_selector – the weighted collection of trial moves.

The constructor checks the setup before the first energy call and raises ValueError when:

an entry of species is not an element symbol (molecular species go in molecules),
a move exchanges a species that is not declared or has no entry in mu,
a species is inserted and deleted with different probabilities per trial in the same cell (see Moves),
a molecule_id group matches no molecule template (see Simulating molecular adsorbates).

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. For atomic moves, \(N\) is the total atom count before the move: every atom of every species, the fixed substrate included. This departs from the textbook per-species count on purpose, to keep the group’s earlier results reproducible (docs/adr/0001-gcmc-de-broglie-total-atom-count.md). Molecular moves use the textbook count, the number of molecules of the exchanged species in the cell. GCMC acceptance: the de Broglie particle count records both conventions. The thermal de Broglie wavelength is

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

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 import GrandCanonicalEnsemble

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

The current configuration is appended to traj_file (extended XYZ, with energy= and Lattice= on the comment line) every trajectory_write_interval steps, whether or not the last trial was accepted.
The header of outfile lists the run parameters, including the mcpy version (plus the git describe output when mcpy runs from a git checkout).
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.

Particle count against MC trial for a Lennard-Jones GCMC run. — The signature of a converged GCMC run: the particle count climbs from the empty box and settles into a fluctuating steady state around the LAMMPS `fix gcmc` average for the same conditions (Lennard-Jones fluid, \(\mu^* = -2.4\), 30 000 trials, from the executed GCMC-basics notebook).

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: the current configuration every trajectory_write_interval steps, whether or not the last trial was accepted, plus the starting configuration as frame 0. 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. A new ensemble starts the file anew; a later run() of the same ensemble continues it.
'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, each replica’s best configuration is collected; the MPI variant gathers per-rank bests to rank 0, and the batched variant reads its in-process replicas.

If all replicas share one \(\mu\), their scores \(\Omega\) are comparable. The lowest-score configuration is written as a single XYZ frame.
If the replicas have different \(\mu\), each score uses its own replica’s \(\mu\), so no common score exists. This is true on a \(\mu\) ladder, and on a temperature ladder where gcmc_factory sets a different \(\mu\) on each rung. The file then holds one frame per rung, in rung order. Each frame is tagged rank=<i> and carries that rung’s own-\(\mu\) score.

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 – the maximum number of random points used in the estimate. The first estimate draws all of them. Later estimates draw only enough for about 0.5% relative noise at the free fraction of the previous estimate (see Cells).

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.

At each exchange attempt, rank 0 draws one random number. It selects the even pairing (0-1, 2-3, …) or the odd pairing (1-2, 3-4, …), so every replica mixes with both neighbours. Each pair accepts or rejects its swap with the rule in Swap rule. After an accepted swap, the two ranks exchange configurations and energies. Each rank keeps its own temperature and chemical potentials.

Launch one MPI rank per replica, e.g. mpirun -n 4 python -m mpi4py my_run.py for a four-rung ladder; Replica exchange on one GPU or over MPI has a complete script. python -m mpi4py makes an uncaught exception on one rank abort all ranks.

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;
the factory must seed every rank independently: the ensemble, the move selector, every move and every cell. Derive each seed from rank. If two ranks start any random stream in the same state, the constructor raises ValueError on every rank, because common random numbers couple the replicas and bias the swap rule;
at the end of the run, rank 0 collects each replica’s running best (via comm.gather) and writes it to global_minimum_file (default "global_minimum.xyz"; pass None to disable): one frame on a temperature ladder, one frame per rung on a \(\mu\) ladder. See Running minima and basin-hopping output;
an exception on any rank aborts the whole MPI job (comm.Abort), because the other ranks would otherwise wait forever in a collective call.

Swap rule

Slot \(k\) gives a configuration \(Z\) the weight

\[w_k(Z) \propto e^{-\beta_k \Phi_Z^{(k)}} \prod_s \Lambda_s(T_k)^{-3 N_{Z,s}}, \qquad \Phi_Z^{(k)} = E_Z - \sum_s \mu_s^{(k)} N_{Z,s}.\]

\(\Phi_Z^{(k)}\) is the grand potential of \(Z\) with the chemical potentials of slot \(k\). \(N_{Z,s}\) is the number of particles of species \(s\) in \(Z\). A canonical replica has no \(\mu\) and no de Broglie factors, so it uses \(\Phi = E\). A swap moves configuration \(X\) from slot \(i\) to slot \(j\), and configuration \(Y\) from slot \(j\) to slot \(i\). The swap is accepted with \(P_{\mathrm{exch}} = \min(1, e^{\Delta})\), where

\[\Delta = \beta_i \Phi_X^{(i)} + \beta_j \Phi_Y^{(j)} - \beta_i \Phi_Y^{(i)} - \beta_j \Phi_X^{(j)} + 3 \sum_s (N_{X,s} - N_{Y,s}) \ln\frac{\Lambda_s(T_i)}{\Lambda_s(T_j)}.\]

The last term comes from the de Broglie factors. Because \(\Lambda \propto T^{-1/2}\), it equals \(\tfrac{3}{2} \sum_s (N_{X,s} - N_{Y,s}) \ln(T_j / T_i)\). It is zero when the two slots have the same temperature, and in LJ units, where \(\Lambda = 1\). This is the hyper-parallel tempering rule of Yan and de Pablo (doi:10.1063/1.480282).

The same rule applies to every ladder:

On a temperature ladder with one shared \(\mu\), the rule becomes \(\Delta = (\beta_j - \beta_i)(\Phi_Y - \Phi_X)\) plus the de Broglie term. Do not use bare energies here, because the replicas can have different \(N\).
On a \(\mu\) ladder at one temperature, the rule becomes \(\Delta = \beta \sum_s (\mu_s^{(i)} - \mu_s^{(j)})(N_{Y,s} - N_{X,s})\). The energies cancel.
If gcmc_factory sets a different \(\mu\) on each temperature rung, for example a \(\mu(T)\) with a thermal correction, the cross terms keep the rule exact.

A non-finite \(\Delta\), for example from a NaN replica energy, rejects the swap and logs a warning.

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. Shared objects or equal random streams between two replicas raise ValueError.

Swaps use the same rule as ReplicaExchange (see Swap rule).

Output layout mirrors ReplicaExchange: per-replica outfile/traj_file/ minima_file (typically tagged with rank) plus a single global_minimum_file written after the run completes (see Running minima and basin-hopping output).

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.