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:

- :class:`CanonicalEnsemble` -- fixed composition (NVT) Metropolis sampling.
- :class:`GrandCanonicalEnsemble` -- variable composition at fixed :math:`(\mu, V, T)`.
- :class:`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 :math:`i` with energy :math:`E_i`, the canonical Boltzmann weight is

.. math::

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

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

.. math::

   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 :math:`\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 :math:`(N \rightarrow N-1)`:

.. math::

   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 :math:`(N \rightarrow N+1)`:

.. math::

   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 :math:`z = e^{\beta \mu}` is the species activity, and

.. math::

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

is the thermal de Broglie wavelength. The geometric volume :math:`V` of the textbook
formulation is replaced by the accessible free volume :math:`V_{\mathrm{free}}` returned by
the configured cell -- see :ref:`free-volume`.

.. code-block:: python

   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:

Free volume and `species_radii`
-------------------------------

In the textbook GCMC equations, :math:`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 :math:`V_{\mathrm{free}}` by Monte Carlo sampling and that
quantity is used in the acceptance rules above.

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

.. math::

   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}

The occupied fraction and the resulting free volume are

.. math::

   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 :doc:`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 :math:`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

.. math::

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

and an analogous expression involving :math:`\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
---------------------------

- :class:`CanonicalEnsemble` -- fixed-composition sampling, structural relaxation, or
  thermal annealing.
- :class:`GrandCanonicalEnsemble` -- variable-composition adsorption/desorption studies at
  fixed :math:`(\mu, T)`, e.g. oxidation phase diagrams.
- :class:`ReplicaExchange` -- when a single GCMC trajectory is prone to getting trapped in
  a single basin and broader thermodynamic coverage is needed.