Fix/review 2026 07 27

.github/workflows/tests.yml

@@ -34,3 +34,23 @@ jobs:
 
       - name: Run tests
         run: pytest tests/ -v
+
+  lint:
+    runs-on: ubuntu-latest
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: actions/setup-python@v5
+        with:
+          python-version: '3.12'
+          cache: 'pip'
+
+      - name: Install flake8
+        run: |
+          python -m pip install --upgrade pip
+          pip install flake8
+
+      # Same command CLAUDE.md documents; config lives in .flake8.
+      - name: Lint
+        run: flake8 mcpy/

CHANGELOG.md

@@ -4,6 +4,27 @@ All notable changes to this project are documented here. The format is based on
 [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project
 adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [1.4.0] - 2026-07-28
+
+### Fixed
+- **Batched replica exchange accepted every chemical-potential swap.** `BatchedReplicaExchange._accept_swap` implemented only the temperature-ladder criterion `(beta_j - beta_i)(Phi_j - Phi_i)`, which is identically zero when the replicas share a temperature: every mu-ladder swap was accepted with p = 1, configurations random-walked freely across the ladder, and no replica sampled its own mu. **Results produced by `BatchedReplicaExchange(..., mus=[...])` in 1.3.0 or earlier are invalid and need re-running**; temperature ladders were always correct, as was the MPI `ReplicaExchange`, which selected its mu criterion all along. The replacement carries the cross terms, `beta_i Phi_X^(i) + beta_j Phi_Y^(j) - beta_i Phi_Y^(i) - beta_j Phi_X^(j)`, where `Phi_Z^(k)` is a configuration scored with slot `k`'s chemical potentials; it reduces to the previous expression for a shared mu, to `beta (mu_i - mu_j)(N_j - N_i)` for a shared temperature, and additionally covers a joint (T, mu) ladder. The failure was not a statistical degradation: an EMT Ag(111)/Au ladder gives `N_Au = [0, 0, 0, 14]` with the fix and `[0, 0, 0, 0]` without it, because unconditional swapping averages the rungs together and erases the mu dependence entirely.
+- **`PermutationMove` could mutate the configuration and then report that it had not.** With `n_swaps > 1`, drawing a species absent from the system in a later iteration returned the "could not propose" sentinel after earlier swaps had already been applied; both GCMC loops read that sentinel as "the atoms were not touched" and skipped the rollback, leaving the configuration out of step with the stored energy for the rest of the run. The usable species are now resolved once, before any mutation (species counts are swap-invariant, so a species absent at the start is absent for the whole trial), and both ensembles restore their pre-trial snapshot on the sentinel path rather than trusting the contract.
+- Molecules whose center of mass drifted above a `CustomCell`'s top stopped being deletion candidates and dropped out of the per-species de Broglie count, inflating `V/((N+1)Lambda^3)` into the runaway insertion mode documented in `docs/gcmc_acceptance_convention.rst`. Molecule candidacy now goes through the new `is_point_exchangeable` predicate, which applies the same dropped z upper bound that `get_atoms_specie_inside_cell` already applied to single atoms; molecules below the cell floor stay excluded, as atoms do.
+- `AlchemiCalculator(energy_only=True)` discarded `'forces'` from a `model_config` that a pre-loaded `MACEWrapper` shares with every calculator built from it, silently disabling FIRE relaxation in an `AlchemiFCalculator` depending on construction order. The combination now raises with the workaround in the message.
+- `GrandCanonicalEnsemble.set_state` and `CanonicalEnsemble.set_state` restored the step count and exchange statistics from the incoming state. `ReplicaExchange` passes a full `get_state()` dict on every accepted swap, so the two ranks traded their swap tallies and the per-rank "Accepted Exchange (%)" column was meaningless. Only the configuration travels now.
+- Per-interval acceptance ratios were never cleared when the outfile was disabled, silently degrading `interval_ratios()` into `total_ratios()` for the rest of the run.
+- A missing `species_radii` entry raised a bare `KeyError` from inside the free-volume sampler; the error now names the cell and the missing species, and notes that species inserted during the run need radii too.
+
+### Added
+- `Cell.is_point_exchangeable(point)`: the point counterpart of `get_atoms_specie_inside_cell`, used by the molecule moves to decide which molecules the reservoir may take back. It defaults to `is_point_inside`, so the box, spherical and dome cells are unchanged; `CustomCell` overrides it. `MoleculeDisplacementMove` tests the same predicate for its region guard, so candidacy and displacement always agree on the region.
+- `BatchedReplicaExchange` warns at the end of a run when the whole-run swap tally is all-accept or all-reject, the two ways a ladder stops being a ladder. Checked on the whole-run tally rather than live, because replicas that have not differentiated yet legitimately accept every early swap.
+- `--chunk-size` on `examples/re_gcmc_co_cupd_batched.py`: peak GPU memory follows the largest chunk rather than the replica count, which is what lets a correctly spaced ladder run at all (an acceptance-equalized ladder for that system needs ~29 rungs, i.e. ~13k atoms in the relax batch, several times the whole-batch ceiling).
+- `docs/replica_exchange_ladder_spacing.rst`: how to detect a dead ladder (read second-half swap acceptance, never the cumulative column, which early free swaps inflate permanently), why uniform mu spacing cannot work across a coverage range (`dN/dmu` grows with coverage while `dmu` does not), and the `dmu = 1/sqrt(beta dN/dmu)` spacing rule with a reference implementation. Worked example on CO/Cu375Pd30 at 400 K: 29 rungs instead of 5 gives second-half acceptance of 18-65% (median 40%) with no dead pair, against three of four pairs at exactly zero accepted swaps on the uniform ladder, and GPU utilisation of 94-98% instead of 43%.
+- CI runs `flake8 mcpy/`, the lint command the contributor docs already specified but the workflow never invoked.
+
+### Changed
+- The outfile acceptance-ratio header abbreviates the molecule moves distinctly (`MolIns`, `MolDel`, `MolDis`) instead of collapsing all three to an indistinguishable `Mol`. Single-word move labels (`Ins`, `Del`, `Dis`, `Per`, `Sha`, `Bro`) are unchanged, so outfiles from atomic runs stay comparable across versions.
+
 ## [1.3.0] - 2026-07-08
 
 ### Added
@@ -87,6 +108,8 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 Initial public release.
 
+[1.4.0]: https://github.com/farrisric/mcpy/compare/v1.3.0...v1.4.0
+[1.3.0]: https://github.com/farrisric/mcpy/compare/v1.2.0...v1.3.0
 [1.2.0]: https://github.com/farrisric/mcpy/compare/v1.1.0...v1.2.0
 [1.1.0]: https://github.com/farrisric/mcpy/compare/v1.0.0...v1.1.0
 [1.0.0]: https://github.com/farrisric/mcpy/releases/tag/v1.0.0

docs/gcmc_acceptance_convention.rst

@@ -152,6 +152,21 @@ Atomic moves are unaffected: for them ``get_exchange_count()`` returns
 ``None`` and the ensembles fall back to the total-atom count documented
 above.
 
+Because molecular moves *do* use the per-species count, they are exposed to
+the runaway mode described in the previous section: any escape route out of
+the counted region undercounts ``N`` and inflates
+:math:`V/((N+1)\Lambda^3)`. Molecule candidacy therefore goes through
+``cell.is_point_exchangeable(com)``, not ``cell.is_point_inside(com)``. The
+two differ for :class:`mcpy.cell.CustomCell`, whose exchangeable region drops
+the z upper bound exactly as ``get_atoms_specie_inside_cell`` already does for
+single atoms, so a molecule that desorbs above the cell top stays both
+countable and deletable. The floor exclusion is kept on both paths: molecules
+that sink below ``bottom_z`` remain excluded, and the buried-species
+irreversibility noted above applies to them unchanged.
+:class:`mcpy.moves.MoleculeDisplacementMove` tests the same predicate for its
+one-way-door guard, so the candidacy region and the region a molecule may be
+displaced within always agree.
+
 The textbook form above holds for ``min_insert=None``. When ``min_insert``
 is set, ``MoleculeInsertionMove`` retries the random position/orientation
 draw (up to 1000 times) against the cell's ``species_radii`` atoms until it

docs/index.rst

@@ -61,6 +61,7 @@ the hybrid GCMC method it implements (see :doc:`bibliography`).
    molecular_adsorbates
    phase_diagrams
    gcmc_acceptance_convention
+   replica_exchange_ladder_spacing
 
 .. toctree::
    :caption: Tutorials

docs/reference/calculators.rst

@@ -89,7 +89,11 @@ chunk_size=None)``, and ``run_md(...)``.
   ``get_potential_energies``; caps peak GPU memory at one chunk. ``None``
   evaluates the whole batch in one pass. See :doc:`../calculators`.
 - ``energy_only`` (bool): drop force computation (no autograd graph) for a
-  ~12% memory saving. Energy is unchanged; forces are unavailable.
+  ~12% memory saving. Energy is unchanged; forces are unavailable. The switch
+  lives on the model's own config, so it cannot be scoped to one calculator:
+  combining it with a pre-loaded ``MACEWrapper`` (which other calculators may
+  share, including a relaxing ``AlchemiFCalculator``) raises ``ValueError``.
+  Pass a checkpoint path instead so this calculator loads its own model.
 
 
 AlchemiFCalculator

docs/reference/cells.rst

@@ -32,9 +32,18 @@ The full ASE simulation box as the active region. Random points are sampled
 uniformly in fractional coordinates. Its volume is the fixed box volume, so
 ``calculate_volume`` does no free-volume sampling.
 
-- ``species_radii`` (dict, optional): per-element exclusion radii.
+- ``species_radii`` (dict, optional): per-element exclusion radii. Every
+  species the system can hold needs one, including species inserted during the
+  run; a missing entry raises ``ValueError`` from ``calculate_volume``.
 - ``seed`` (int, optional): cell-local RNG seed.
 
+Every cell exposes two point-membership predicates. ``is_point_inside(point)``
+bounds the *proposal* region -- the one ``get_random_point`` samples.
+``is_point_exchangeable(point)`` bounds the region the reservoir may take
+molecules back from, and is what the molecule moves use to pick candidates.
+They coincide for every cell except ``CustomCell``; see
+:doc:`../gcmc_acceptance_convention`.
+
 
 CustomCell
 ----------

docs/replica_exchange_ladder_spacing.rst

@@ -0,0 +1,192 @@
+Spacing a replica-exchange ladder
+=================================
+
+A replica-exchange run can look perfectly healthy and still be doing nothing.
+This note records how to tell, and how to space a chemical-potential ladder so
+that it works. It is the ladder-design companion to
+:doc:`gcmc_acceptance_convention`, and it exists because the failure modes here
+are silent: no error, no warning, and output that reads as plausible physics.
+
+.. contents::
+   :local:
+   :depth: 1
+
+
+Read second-half acceptance, never the cumulative column
+--------------------------------------------------------
+
+Every replica starts from the same configuration, so early swaps are free: two
+rungs holding identical states always exchange. Those free successes are
+recorded permanently in the cumulative tally, which therefore overstates the
+ladder's health for the rest of the run -- badly, and for a long time.
+
+A five-rung CO/CuPd run reported cumulative per-slot acceptances of 28.6, 19.6,
+5.9, 2.0 and 0.0 %. That reads as a ladder that mixes well at one end and
+poorly at the other. In the run's *second half*, three of the four pairs
+accepted **exactly zero** swaps: only pair 0-1 was still functioning, and the
+other four replicas were independent single-:math:`\mu` chains.
+
+The test needs no extra instrumentation. Attempts accumulate linearly, so
+between the run's midpoint and its end the attempt count doubles. If a
+cumulative percentage *halves exactly* over that span, the numerator never
+changed and there were no successes at all in the second half:
+
+.. code-block:: text
+
+   second-half rate = 2 * cum_end - cum_mid
+
+Slots ``0`` and ``n-1`` each belong to exactly one pair, so their columns read
+those two pairs directly; interior slots average the pair on either side.
+
+``BatchedReplicaExchange`` warns at the end of a run when the whole-run tally is
+all-accept or all-reject, which catches a catastrophically mis-spaced ladder.
+It cannot catch the case above, where a ladder starts healthy and dies as the
+replicas differentiate. For that, run the halving test.
+
+
+Why uniform spacing cannot work
+-------------------------------
+
+For a :math:`\mu` ladder at a single temperature the swap exponent is
+
+.. math::
+
+   \ln \frac{w'}{w} = \beta\,(\mu_i - \mu_j)(N_j - N_i)
+                    = -\beta\,\Delta\mu\,\Delta N .
+
+:math:`\Delta N` between neighbouring rungs is not a constant: it grows with
+coverage, because :math:`dN/d\mu` does. A spacing chosen where the surface is
+bare is therefore far too coarse once the surface fills. In the run above, at
+400 K with :math:`\Delta\mu = 0.2` eV and :math:`\Delta N \approx 16`, the
+exponent reaches :math:`29 \times 0.2 \times 16 \approx 93`, i.e.
+:math:`p \sim 10^{-40}`. No amount of sampling recovers that pair.
+
+
+The spacing rule
+----------------
+
+Grand-canonical fluctuation-dissipation gives the width of each rung's
+:math:`N` distribution,
+
+.. math::
+
+   \sigma_N^2 = \frac{1}{\beta}\frac{dN}{d\mu} ,
+
+and neighbouring rungs exchange when their distributions overlap, i.e. when
+:math:`\beta\,\Delta\mu\,\sigma_N \sim 1`. Substituting:
+
+.. math::
+
+   \Delta\mu(\mu) = \frac{1}{\sqrt{\beta\,dN/d\mu}} ,
+   \qquad
+   n_\text{rungs} = \int \sqrt{\beta\,\frac{dN}{d\mu}}\; d\mu .
+
+Use :math:`dN/d\mu` from a measured isotherm, **not** the per-rung
+:math:`\sigma_N` observed in an unconverged run: a stuck rung reports a
+:math:`\sigma_N` far below equilibrium and a drifting one far above, because the
+"fluctuation" is drift.
+
+Histogram reweighting is not an alternative route to this. Reweighting requires
+overlap between adjacent rungs, and a mis-spaced ladder has none -- that absence
+is the disease being diagnosed.
+
+.. code-block:: python
+
+   import numpy as np
+
+   def acceptance_equalized_mu_ladder(mu, n_ads, temperature, mu_min, mu_max):
+       """Rung positions for a chemical-potential ladder of even acceptance.
+
+       Places rungs where the cumulative integral of sqrt(beta dN/dmu) crosses
+       successive integers, so every neighbouring pair has comparable
+       distribution overlap and therefore comparable swap acceptance.
+
+       Parameters
+       ----------
+       mu, n_ads : array_like
+           A measured isotherm: chemical potentials and mean adsorbate counts.
+       temperature : float
+           Ladder temperature in K (all rungs share it).
+       mu_min, mu_max : float
+           Range to span.
+
+       Returns
+       -------
+       ndarray
+           Rung chemical potentials, ascending. Spacing widens where the
+           isotherm is flat and tightens where it is steep.
+       """
+       beta = 1.0 / (8.617333e-5 * temperature)
+       mu, n_ads = np.asarray(mu, float), np.asarray(n_ads, float)
+       slope = np.diff(n_ads) / np.diff(mu)          # dN/dmu per interval
+       midpoints = 0.5 * (mu[1:] + mu[:-1])
+
+       fine = np.linspace(mu_min, mu_max, 2001)
+       density = np.sqrt(beta * np.interp(fine, midpoints, slope))
+       cumulative = np.concatenate([[0.0], np.cumsum(
+           0.5 * (density[1:] + density[:-1]) * np.diff(fine))])
+       return np.interp(np.arange(0.0, cumulative[-1], 1.0), cumulative, fine)
+
+
+Calibration
+-----------
+
+The rule above is derived from the exponent at the *mean* :math:`\Delta N`, but
+the acceptance that matters is :math:`\langle \min(1, \cdot) \rangle` over both
+rungs' fluctuations. A favourable fluctuation accepts outright, so the realised
+acceptance is higher than the point estimate suggests -- for the CO/CuPd system
+the rule aimed at ~25 % and delivered ~40 %.
+
+Treat the rule as a safe upper bound on rung count, then trim once a run has
+measured the real acceptance. A ladder derived from an unconverged isotherm is
+also a lower bound on the rungs eventually needed, since :math:`dN/d\mu`
+steepens as coverage converges. Those two biases push in opposite directions;
+re-derive from each run rather than trusting one calculation.
+
+
+Worked example: CO on a CuPd nanoparticle
+-----------------------------------------
+
+Cu\ :sub:`375`\ Pd\ :sub:`30` at 400 K, :math:`\Delta\mu` spanning -1.8 to -1.0
+eV, ``AlchemiFCalculator``.
+
+.. list-table::
+   :header-rows: 1
+
+   * - Interval (eV)
+     - :math:`dN/d\mu`
+     - Spacing
+     - Rungs
+   * - -1.8 .. -1.6
+     - 8.5
+     - 0.064 eV
+     - 3
+   * - -1.6 .. -1.4
+     - 45.1
+     - 0.028 eV
+     - 7
+   * - -1.4 .. -1.2
+     - 54.5
+     - 0.025 eV
+     - 8
+   * - -1.2 .. -1.0
+     - 88.9
+     - 0.020 eV
+     - 10
+
+29 rungs, against the 5 that a uniform 0.2 eV ladder would use. Measured
+second-half acceptance: minimum 18 %, median 40 %, maximum 65 %, with no dead
+pair -- and a smooth, strictly monotonic isotherm from 0.1 to 57 CO.
+
+Two practical consequences of going wide:
+
+- **Chunking becomes mandatory.** 29 replicas of ~460 atoms is ~13k atoms,
+  several times the whole-batch relaxation ceiling. Pass ``chunk_size`` so peak
+  memory follows the largest chunk rather than the replica count.
+- **It is faster, not slower, in absolute terms.** A narrow ladder leaves the
+  GPU idle. Five replicas ran at ~43 % utilisation; 29 replicas at
+  ``chunk_size=8`` reached 94-98 % and 15.6 GB of 32.6 GB.
+
+Once the ladder is healthy, convergence becomes the binding constraint rather
+than the sampler. Chain runs with ``--init-dir`` pointing at the previous
+output and watch the tail-half drift per rung, not just the acceptance.

examples/re_gcmc_co_cupd_batched.py

@@ -46,9 +46,12 @@ def parse_args():
     p.add_argument('--gcmc-steps', type=int, default=80,
                    help='GCMC steps per replica')
     p.add_argument('--exchange-interval', type=int, default=8,
-                   help='Steps between replica-exchange attempts; use >= 50 '
-                        'for per-mu coverage curves (small values over-mix '
-                        'the ladder)')
+                   help='Steps between replica-exchange attempts. Frequent '
+                        'exchange is fine: a correct swap rule leaves each '
+                        "rung's distribution intact. (The old advice to use "
+                        '>= 50 here was working around the mu-ladder swap '
+                        'bug, which accepted every swap and so did flatten '
+                        'the coverage curve the more often it ran.)')
     p.add_argument('--checkpoint', default='medium-mpa-0')
     p.add_argument('--rel-steps', type=int, default=30)
     p.add_argument('--rel-fmax', type=float, default=0.1)
@@ -61,6 +64,13 @@ def parse_args():
     p.add_argument('--mol-disp-angle', type=float, default=None,
                    help='Max rotation angle of the rigid move (rad); '
                         'None = full uniform rotation')
+    p.add_argument('--chunk-size', type=int, default=None,
+                   help='Relax at most this many replicas per forward pass. '
+                        'Peak GPU memory follows the largest chunk rather '
+                        'than the replica count, which is what makes a wide '
+                        'ladder fit: acceptance-equalized spacing needs ~29 '
+                        'rungs here, and 29 x ~460 atoms is several times the '
+                        'whole-batch ceiling. None relaxes the whole batch.')
     p.add_argument('--no-compile', action='store_true')
     p.add_argument('--seed', type=int, default=7,
                    help='Master seed for moves, RE, and the Pd placement')
@@ -108,6 +118,7 @@ def main():
         steps=args.rel_steps,
         fmax=args.rel_fmax,
         compile_model=not args.no_compile,
+        chunk_size=args.chunk_size,
     )
 
     e_co = calculator.get_potential_energy(