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.

Three ensembles are available:

CanonicalEnsemble – fixed composition (NVT) Metropolis sampling.
GrandCanonicalEnsemble – variable composition at fixed \((\mu, V, T)\).
ReplicaExchange – a wrapper that runs several GCMC replicas at different temperatures or chemical potentials and periodically attempts exchanges between them.

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

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.

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.

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 and log paths) so that ranks do not race on the same files;
one MPI rank corresponds to one replica.

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 – when a single GCMC trajectory is prone to getting trapped in a single basin and broader thermodynamic coverage is needed.