diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 21dcb90..fddfb1e 100644 --- a/.github/workflows/tests.yml +++ b/.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/ diff --git a/CHANGELOG.md b/CHANGELOG.md index 69d96a3..cdf226d 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/docs/gcmc_acceptance_convention.rst b/docs/gcmc_acceptance_convention.rst index e9d8516..0798eb9 100644 --- a/docs/gcmc_acceptance_convention.rst +++ b/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 diff --git a/docs/index.rst b/docs/index.rst index 137c59e..7bdaf81 100644 --- a/docs/index.rst +++ b/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 diff --git a/docs/reference/calculators.rst b/docs/reference/calculators.rst index b1e2855..03c3740 100644 --- a/docs/reference/calculators.rst +++ b/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 diff --git a/docs/reference/cells.rst b/docs/reference/cells.rst index 6534f83..4e966e1 100644 --- a/docs/reference/cells.rst +++ b/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 ---------- diff --git a/docs/replica_exchange_ladder_spacing.rst b/docs/replica_exchange_ladder_spacing.rst new file mode 100644 index 0000000..54ea9d4 --- /dev/null +++ b/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. diff --git a/examples/re_gcmc_co_cupd_batched.py b/examples/re_gcmc_co_cupd_batched.py index d604b42..f444a6f 100644 --- a/examples/re_gcmc_co_cupd_batched.py +++ b/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( diff --git a/mcpy/calculators/alchemi_calculator.py b/mcpy/calculators/alchemi_calculator.py index 2b6c4ee..e9cebd2 100644 --- a/mcpy/calculators/alchemi_calculator.py +++ b/mcpy/calculators/alchemi_calculator.py @@ -54,6 +54,20 @@ def __init__( energy_only: bool = False, head: Union[str, int, None] = None, ) -> None: + if energy_only and isinstance(checkpoint, MACEWrapper): + # energy_only mutates the wrapper's own model_config (below), and a + # pre-loaded wrapper is shared by every calculator built from it — + # including an AlchemiFCalculator whose FIRE relaxation needs the + # forces this would switch off. There is no way to scope the change + # to one calculator, so refuse instead of breaking the other one + # silently and order-dependently. + raise ValueError( + "energy_only=True cannot be combined with a pre-loaded " + "MACEWrapper: dropping 'forces' from its active outputs would " + "also disable forces for every other calculator sharing that " + "wrapper. Pass the checkpoint path so this calculator loads " + "its own model, or set energy_only on all of the sharers." + ) self.device = device self.dtype = dtype self.max_neighbors = max_neighbors diff --git a/mcpy/cell/cell.py b/mcpy/cell/cell.py index 7d172d6..581a785 100644 --- a/mcpy/cell/cell.py +++ b/mcpy/cell/cell.py @@ -66,6 +66,28 @@ def get_species(self): """ return list(self.species_radii.keys()) + def _radii_for(self, atoms): + """Per-atom exclusion radii for the free-volume samplers. + + Raises a message that names the missing symbols instead of the bare + ``KeyError`` a dict lookup would throw from inside the sampler. The + check lives here rather than in ``__init__`` because GCMC inserts + species that are absent when the cell is built. + + :return: ndarray of radii, one per atom, in ``atoms`` order. + """ + symbols = atoms.get_chemical_symbols() + missing = sorted(set(symbols) - set(self.species_radii)) + if missing: + raise ValueError( + f'{type(self).__name__}.species_radii has no radius for ' + f'{missing}. Every species the system can contain needs one, ' + f'including species inserted during the run; got ' + f'{sorted(self.species_radii)}.' + ) + return np.fromiter((self.species_radii[s] for s in symbols), + dtype=float, count=len(symbols)) + def is_point_inside(self, point): """The box cell spans the whole periodic cell: every point is inside. @@ -73,3 +95,17 @@ def is_point_inside(self, point): any cell type (the region cells implement a real test). """ return True + + def is_point_exchangeable(self, point): + """Whether a molecule whose center of mass sits at ``point`` may be + exchanged with the reservoir (counted, deleted, displaced). + + The point counterpart of :meth:`get_atoms_specie_inside_cell`, and the + predicate the molecule moves use. It is separate from + :meth:`is_point_inside` because a region may deliberately accept + molecules it would never *propose* — :class:`CustomCell` drops the z + upper bound so a desorbed molecule stays deletable instead of + accumulating forever. Cells whose two regions coincide (the box, the + sphere, the dome) inherit this delegation unchanged. + """ + return self.is_point_inside(point) diff --git a/mcpy/cell/custom_cell.py b/mcpy/cell/custom_cell.py index a17374a..d796e88 100644 --- a/mcpy/cell/custom_cell.py +++ b/mcpy/cell/custom_cell.py @@ -75,12 +75,7 @@ def calculate_volume(self, atoms) -> float: cart_coords = frac_coords @ self.dimensions # in cell frame positions = self._periodic_images(atoms) - symbols = atoms.get_chemical_symbols() - radii = np.fromiter( - (self.species_radii[s] for s in symbols), - dtype=float, count=n_atoms, - ) - radii = np.tile(radii, len(positions) // n_atoms) + radii = np.tile(self._radii_for(atoms), len(positions) // n_atoms) covered = np.zeros(self.mc_sample_points, dtype=bool) for r in np.unique(radii): @@ -144,6 +139,21 @@ def is_point_inside(self, point): frac = (point - self.offset) @ self._dim_inv return bool(np.all((frac >= 0.0) & (frac < 1.0))) + def is_point_exchangeable(self, point): + """Same asymmetry as :meth:`get_atoms_specie_inside_cell`: inside the + xy footprint and at or above the cell floor, with the z upper bound + dropped. + + Without this a molecule that desorbs above the cell top would stop + being a deletion candidate *and* would drop out of the per-species + count in the acceptance factor ``V/((N+1)Λ³)`` -- the runaway + insertion mode described in docs/gcmc_acceptance_convention.rst. + Molecules below the floor stay excluded, like atoms. + """ + frac = (point - self.offset) @ self._dim_inv + return bool(np.all(frac[:2] >= 0.0) and np.all(frac[:2] < 1.0) + and frac[2] >= 0.0) + def get_species(self): """ Get the species present in the custom cell. diff --git a/mcpy/cell/dome_cell.py b/mcpy/cell/dome_cell.py index 1a5b71c..2a3990a 100644 --- a/mcpy/cell/dome_cell.py +++ b/mcpy/cell/dome_cell.py @@ -121,11 +121,7 @@ def calculate_volume(self, atoms): self.volume = dome_fraction * self.sphere_volume return - symbols = atoms.get_chemical_symbols() - radii = np.fromiter( - (self.species_radii[s] for s in symbols), - dtype=float, count=n_atoms, - ) + radii = self._radii_for(atoms) positions = atoms.positions covered = np.zeros(len(pts), dtype=bool) diff --git a/mcpy/cell/spherical_cell.py b/mcpy/cell/spherical_cell.py index 75c919f..e1bbb03 100644 --- a/mcpy/cell/spherical_cell.py +++ b/mcpy/cell/spherical_cell.py @@ -121,11 +121,7 @@ def calculate_volume(self, atoms): pts = self._sample_sphere_points(self.mc_sample_points) - symbols = atoms.get_chemical_symbols() - radii = np.fromiter( - (self.species_radii[s] for s in symbols), - dtype=float, count=n_atoms, - ) + radii = self._radii_for(atoms) positions = atoms.positions covered = np.zeros(self.mc_sample_points, dtype=bool) diff --git a/mcpy/ensembles/batched_replica_exchange.py b/mcpy/ensembles/batched_replica_exchange.py index 3c2d117..b5e4ae0 100644 --- a/mcpy/ensembles/batched_replica_exchange.py +++ b/mcpy/ensembles/batched_replica_exchange.py @@ -36,6 +36,9 @@ logger = logging.getLogger(__name__) +# Below this many attempts an all-accept / all-reject tally says nothing. +_SWAP_DEGENERACY_MIN_ATTEMPTS = 20 + class BatchedReplicaExchange: def __init__( @@ -224,9 +227,12 @@ def _batched_single_move(self, active: List[int]) -> None: result if isinstance(result, tuple) else (result, 0, None) ) if atoms_new is False or atoms_new is None: - # Move couldn't propose — atoms unchanged, snapshot harmless. - # Identity check, not truthiness: an empty Atoms (last atom - # deleted) is falsy but is a real proposal. + # Move couldn't propose. Identity check, not truthiness: an + # empty Atoms (last atom deleted) is falsy but is a real + # proposal. Restore rather than trust the sentinel to mean + # "nothing was touched" (see GrandCanonicalEnsemble). + r.atoms.arrays = snapshots[i] + r.atoms.set_constraint(constraint_snapshots[i]) continue if atoms_new is not r.atoms: raise RuntimeError( @@ -284,21 +290,37 @@ def _attempt_exchanges(self) -> None: self.exchange_successes[j] += 1 def _accept_swap(self, i: int, j: int) -> bool: - """Temperature-RE Metropolis: P = exp((β_j - β_i)(Φ_j - Φ_i)). + """Replica-exchange Metropolis for a temperature *or* a μ ladder. Replicas here are grand-canonical (fixed μ, fluctuating N), so configs - at different temperatures must be compared through the grand potential - Φ = E - Σ_s μ_s N_s rather than bare energy. With no chemical potential - this reduces to the standard energy-only swap. + must be compared through the grand potential Φ = E - Σ_s μ_s N_s rather + than bare energy. Exchanging config X (slot i) with config Y (slot j) + changes the joint weight by + + ln(w'/w) = β_i Φ_X^(i) + β_j Φ_Y^(j) - β_i Φ_Y^(i) - β_j Φ_X^(j) + + where Φ_Z^(k) is configuration Z scored with slot k's chemical + potentials. The cross terms are what makes this work on both ladders: + with a shared μ it collapses to the familiar (β_j - β_i)(Φ_j - Φ_i), + but on a μ ladder (one shared temperature) that shortened form is + identically zero and would accept every swap, destroying the ladder. + With shared β it reduces to β Σ_s (μ_i,s - μ_j,s)(N_Y,s - N_X,s), the + rule ``ReplicaExchange._exchange_prob_mu`` implements for MPI. """ ri, rj = self.replicas[i], self.replicas[j] beta_i, beta_j = ri.units.beta, rj.units.beta - phi_i, phi_j = self._grand_potential(ri), self._grand_potential(rj) - delta = (beta_j - beta_i) * (phi_j - phi_i) - p = min(1.0, float(np.exp(delta))) + phi_ii = self._grand_potential(ri) # config X, μ of i + phi_jj = self._grand_potential(rj) # config Y, μ of j + phi_ji = ri._minimum_score(rj.atoms, rj.E_old) # config Y, μ of i + phi_ij = rj._minimum_score(ri.atoms, ri.E_old) # config X, μ of j + delta = (beta_i * phi_ii + beta_j * phi_jj + - beta_i * phi_ji - beta_j * phi_ij) + # exp() overflows to inf for a large positive delta; short-circuit. + p = 1.0 if delta >= 0.0 else float(np.exp(delta)) self.logger.debug( - "swap %d<->%d: beta_i=%.3e beta_j=%.3e Phi_i=%.3f Phi_j=%.3f delta=%.3f p=%.3f", - i, j, beta_i, beta_j, phi_i, phi_j, delta, p, + "swap %d<->%d: beta_i=%.3e beta_j=%.3e Phi_ii=%.3f Phi_jj=%.3f " + "Phi_ji=%.3f Phi_ij=%.3f delta=%.3f p=%.3f", + i, j, beta_i, beta_j, phi_ii, phi_jj, phi_ji, phi_ij, delta, p, ) return self.rng.get_uniform() < p @@ -364,8 +386,42 @@ def _log_status(self, step: int) -> None: self.logger.info('RE %d/%d | N: %s | E(eV): %s | swap acc %s', step, self.gcmc_steps, n, e, swap) + def _warn_if_ladder_degenerate(self) -> None: + """Flag a ladder that stopped behaving like one. + + A ladder accepting *every* swap is not exchanging information between + rungs, it is averaging them together; one accepting *none* leaves each + replica sampling in isolation. Neither shows up as an error, and the + all-accept case in particular yields a plausible-looking but + μ-independent isotherm, so it needs saying out loud. + + Checked once on the whole-run tally rather than live: a ladder whose + replicas have not differentiated yet (identical starting configs, an + adsorbate that has not begun to adsorb) legitimately accepts every + early swap, and a warning that fires during equilibration would be + noise. The trade is that this is a post-mortem. + """ + attempts, successes = sum(self.exchange_attempts), sum(self.exchange_successes) + if attempts < _SWAP_DEGENERACY_MIN_ATTEMPTS: + return + if successes == attempts: + reason = ('every swap was accepted -- neighbouring rungs are ' + 'indistinguishable, so the ladder averages over them ' + 'and per-rung results lose their T/mu dependence') + elif successes == 0: + reason = ('no swap was ever accepted -- rungs are too far apart ' + 'to exchange, so each replica sampled in isolation') + else: + return + self.logger.warning( + 'replica-exchange ladder looks degenerate: %d/%d swaps accepted. ' + 'Widen or tighten the %s spacing.', + successes, attempts, 'mu' if self.mus is not None else 'temperature') + self.logger.warning(' %s', reason) + def _log_summary(self) -> None: """Final consolidated summary: per-replica move acceptances.""" + self._warn_if_ladder_degenerate() for i, r in enumerate(self.replicas): ratios = ' '.join( f'{x:.0%}' if x == x else 'n/a' diff --git a/mcpy/ensembles/canonical_ensemble.py b/mcpy/ensembles/canonical_ensemble.py index 77f1a9a..4465bac 100644 --- a/mcpy/ensembles/canonical_ensemble.py +++ b/mcpy/ensembles/canonical_ensemble.py @@ -95,12 +95,8 @@ def set_state(self, state): # config swap, without relying on the pickled .info surviving MPI. self.atoms.info.setdefault("key_value_pairs", {}) self.atoms.info["key_value_pairs"]["potential_energy"] = state["energy"] - if "step" in state: - self._step = state["step"] - if "exchange_attempts" in state: - self.exchange_attempts = state["exchange_attempts"] - if "exchange_successes" in state: - self.exchange_successes = state["exchange_successes"] + # Step count and exchange statistics stay with the slot, not the + # config -- see GrandCanonicalEnsemble.set_state. def _acceptance_condition(self, potential_diff: float) -> bool: if potential_diff <= 0: diff --git a/mcpy/ensembles/grand_canonical_ensemble.py b/mcpy/ensembles/grand_canonical_ensemble.py index ca8345b..fac1b6d 100644 --- a/mcpy/ensembles/grand_canonical_ensemble.py +++ b/mcpy/ensembles/grand_canonical_ensemble.py @@ -98,13 +98,11 @@ def set_state(self, state: Dict[str, any]) -> None: self.atoms = state["atoms"] self.E_old = state["energy"] self.n_atoms = state["n_atoms"] - # Restore optional bookkeeping if present (used by restart, not by RE). - if "step" in state: - self._step = state["step"] - if "exchange_attempts" in state: - self.exchange_attempts = state["exchange_attempts"] - if "exchange_successes" in state: - self.exchange_successes = state["exchange_successes"] + # Only the configuration travels. Step count and exchange statistics + # describe the *slot*, not the config: ``ReplicaExchange`` passes a + # full ``get_state()`` dict here on every accepted swap, so restoring + # them would trade the two ranks' swap tallies back and forth and make + # the per-rank "Accepted Exchange (%)" column meaningless. # The configuration changed under the cells: refresh their free # volumes now, or the next insertion/deletion acceptance would use # the previous configuration's volume (GCMC only recalculates on @@ -143,12 +141,17 @@ def write_outfile(self, step: int = None, energy: float = None) -> None: :class:`BaseEnsemble` signature but ignored — GCMC always logs its own ``_step``/``E_old`` so the row matches the sampler state. """ - if self._outfile is None or self._outfile_handle is None: - return if self._last_logged_step == self._step: return # already wrote this step (e.g. finalize_run after a triggered write) acceptance_ratios = self.move_selector.interval_ratios() + # Reset before the disabled-outfile bail-out: the counters are + # per-write-interval, so skipping the reset when there is nothing to + # write turns interval_ratios() into total_ratios() for the rest of + # the run. self.move_selector.reset_counters() + if self._outfile is None or self._outfile_handle is None: + self._last_logged_step = self._step + return ratio_str = ", ".join( f"{r * 100:.1f}%" if not np.isnan(r) else "N/A" for r in acceptance_ratios @@ -246,11 +249,18 @@ def do_gcmc_step(self) -> None: atoms_new, delta_particles, species = self.move_selector.do_trial_move(atoms) if atoms_new is False or atoms_new is None: - # Move couldn't be proposed (e.g. empty cell). The move did - # not mutate ``atoms``; MoveSelector already recorded the - # failure so it won't depress the acceptance ratio. Identity - # check, not truthiness: an empty Atoms (last atom deleted) - # is falsy but is a real proposal that must be scored. + # Move couldn't be proposed (e.g. empty cell). MoveSelector + # already recorded the failure so it won't depress the + # acceptance ratio. Identity check, not truthiness: an empty + # Atoms (last atom deleted) is falsy but is a real proposal + # that must be scored. + # + # Restore the snapshot rather than trusting the sentinel to + # mean "nothing was touched": a move that mutates before + # bailing would otherwise leave the configuration changed + # while ``E_old`` still describes the previous one. + atoms.arrays = saved_arrays + atoms.set_constraint(saved_constraints) continue if atoms_new is not atoms: diff --git a/mcpy/moves/molecule_displacement_move.py b/mcpy/moves/molecule_displacement_move.py index c154421..500e7d0 100644 --- a/mcpy/moves/molecule_displacement_move.py +++ b/mcpy/moves/molecule_displacement_move.py @@ -61,11 +61,14 @@ def do_trial_move(self, atoms) -> Atoms: # Reject displacements that carry the COM out of the region: the # reverse proposal would have zero probability (the molecule stops # being a candidate), so crossing would be a one-way door that - # strands molecules outside the grand-canonical bookkeeping. + # strands molecules outside the grand-canonical bookkeeping. The + # region tested is the exchangeable one -- the same predicate that + # decides candidacy in ``find_molecules`` -- so the door stays + # two-way for every cell type. new_com = com + shift if np.any(np.asarray(atoms.pbc)): new_com = wrap_positions([new_com], atoms.cell, pbc=atoms.pbc)[0] - if not self.cell.is_point_inside(new_com): + if not self.cell.is_point_exchangeable(new_com): return False, 0, self.name if self.max_angle is None: rotation = random_rotation_matrix(self.rng) diff --git a/mcpy/moves/molecule_utils.py b/mcpy/moves/molecule_utils.py index 8f172a2..2f2a176 100644 --- a/mcpy/moves/molecule_utils.py +++ b/mcpy/moves/molecule_utils.py @@ -42,8 +42,10 @@ def find_molecules(atoms, template_symbols, cell=None): ``template_symbols`` is the sorted list of the template's chemical symbols. When ``cell`` is given, only molecules whose center of mass - lies inside it (``cell.is_point_inside``) are returned; ``None`` skips - the spatial filter (used by the grand-potential bookkeeping). + the cell counts as exchangeable (``cell.is_point_exchangeable``) are + returned; ``None`` skips the spatial filter (used by the grand-potential + bookkeeping). The exchangeable region is not always the proposal region: + see :meth:`mcpy.cell.Cell.is_point_exchangeable`. """ ids = atoms.arrays.get('molecule_id') if ids is None: @@ -56,7 +58,8 @@ def find_molecules(atoms, template_symbols, cell=None): members = np.where(ids == mid)[0] if sorted(symbols[members]) != list(template_symbols): continue - if cell is not None and not cell.is_point_inside(molecule_com(atoms, members)): + if cell is not None and not cell.is_point_exchangeable( + molecule_com(atoms, members)): continue groups.append(members) return groups diff --git a/mcpy/moves/move_selector.py b/mcpy/moves/move_selector.py index 290502d..b4db54e 100644 --- a/mcpy/moves/move_selector.py +++ b/mcpy/moves/move_selector.py @@ -1,9 +1,26 @@ +import re import warnings import numpy as np from mcpy.utils import RandomNumberGenerator +def _abbreviate(class_name): + """Short but still distinct label for a move class. + + A plain three-character slice collapses ``MoleculeInsertionMove``, + ``MoleculeDeletionMove`` and ``MoleculeDisplacementMove`` to the same + 'Mol', which makes the outfile's acceptance-ratio header unreadable for + molecular runs. Split the CamelCase name, drop the trailing 'Move', and + keep three letters per remaining word: 'MolIns', 'MolDel', 'MolDis'. + Single-word moves are unchanged ('Ins', 'Del', 'Dis', 'Per'). + """ + words = re.findall('[A-Z][a-z0-9]*', class_name) or [class_name] + if len(words) > 1 and words[-1] == 'Move': + words = words[:-1] + return ''.join(word[:3] for word in words) + + class MoveSelector: """Class used to randomly select a Monte Carlo Move from a list of Moves. @@ -35,7 +52,7 @@ class MoveSelector: def __init__(self, probabilities, move_list, seed=None, n_moves=None): assert len(probabilities) == len(move_list) self.move_list = move_list - self.move_list_names = [move.__class__.__name__[:3] for move in move_list] + self.move_list_names = [_abbreviate(move.__class__.__name__) for move in move_list] if n_moves is None: n_moves = max(1, round(sum(probabilities))) if n_moves < 1: diff --git a/mcpy/moves/permutation_move.py b/mcpy/moves/permutation_move.py index 6bfda42..61ebb7b 100644 --- a/mcpy/moves/permutation_move.py +++ b/mcpy/moves/permutation_move.py @@ -33,12 +33,18 @@ def do_trial_move(self, atoms): """ nums = atoms.arrays['numbers'] symbols = np.asarray(atoms.get_chemical_symbols()) + # Resolve the usable species once, up front. A swap preserves every + # species count, so a species absent now stays absent for the whole + # trial: bailing out from inside the loop would leave the earlier + # swaps applied while returning the "couldn't propose" sentinel, which + # the ensembles read as "the atoms were not touched". + present = [s for s in self.species if np.any(symbols == s)] + if len(present) < 2: + return False, 0, 'X' for _ in range(self.n_swaps): - species_pair = self.rng.random.sample(self.species, 2) + species_pair = self.rng.random.sample(present, 2) indices_a = np.where(symbols == species_pair[0])[0] indices_b = np.where(symbols == species_pair[1])[0] - if len(indices_a) == 0 or len(indices_b) == 0: - return False, 0, 'X' i = int(self.rng.random.choice(indices_a)) j = int(self.rng.random.choice(indices_b)) nums[i], nums[j] = nums[j], nums[i] diff --git a/paper/paper.bib b/paper/paper.bib index b044d56..07f52f3 100644 --- a/paper/paper.bib +++ b/paper/paper.bib @@ -111,3 +111,24 @@ @article{reuter2001 year = {2001}, doi = {10.1103/PhysRevB.65.035406} } + +@inproceedings{telari2026csi, + author = {Telari, Emanuele and Farris, Riccardo and Quinlivan Dom\'{i}nguez, Jon and Neyman, Konstantin M. and Bruix, Albert}, + title = {Modelling Nanoparticle Transformations in Reactive Environments}, + booktitle = {Book of Abstracts, CSI 2026 -- Cluster-Surface Interactions}, + address = {Paris, France}, + pages = {26}, + year = {2026}, + note = {Conference abstract, 1--4 June 2026}, + url = {https://csi2026.sciencesconf.org/data/pages/Book_Abstracts.pdf} +} + +@inproceedings{bruix2026ai4am, + author = {Bruix, Albert}, + title = {Accelerating the Structural and Chemical Characterization of Nanostructured Materials under Reaction Conditions with {ML}-guided Grand Canonical Global Optimization}, + booktitle = {Book of Abstracts, Artificial Intelligence for Advanced Materials ({AI4AM} 2026)}, + address = {Madrid, Spain}, + year = {2026}, + note = {Conference abstract, 19--21 May 2026}, + url = {https://phantomsfoundation.com/AI4AM/2026/Abstracts/AI4AM2026_Bruix_Albert_58.pdf} +} diff --git a/paper/paper.md b/paper/paper.md index 860bdcf..9a9fd64 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -180,7 +180,11 @@ and surfaces under reactive atmospheres at near-DFT accuracy. The replica-exchange GCMC workflow was built during doctoral research at the Universitat de Barcelona [@farris2025thesis] and underpins a forthcoming manuscript on the oxidation thermodynamics of silver nanoparticles, with further application to additional -catalytic metal/oxide systems. \autoref{fig:phase} shows +catalytic metal/oxide systems. Simulations run with `mcpy` have contributed to +the group's grand canonical modelling of nanostructured catalysts under reaction +conditions, presented at recent international conferences on cluster science and +on machine learning for materials [@telari2026csi; @bruix2026ai4am]. +\autoref{fig:phase} shows a representative result from this workflow: the hydrogenation phase diagram of a 201-atom PdAg nanoparticle, in which `mcpy` recovers the progressive hydrogen loading of the particle as the chemical potential increases. By coupling diff --git a/paper/paper.pdf b/paper/paper.pdf index dd93b40..14a254b 100644 Binary files a/paper/paper.pdf and b/paper/paper.pdf differ diff --git a/pyproject.toml b/pyproject.toml index dff563c..2d41b0f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "mcpy" -version = "1.3.0" +version = "1.4.0" description = "A python package to run atomistic Monte Carlo simulations" readme = "README.md" requires-python = ">=3.11" diff --git a/tests/test_audit_regressions.py b/tests/test_audit_regressions.py index 20ab53a..ac2bfa5 100644 --- a/tests/test_audit_regressions.py +++ b/tests/test_audit_regressions.py @@ -398,3 +398,137 @@ def test_import_does_not_mutate_matplotlib_rcparams(): "matplotlib.rcParams['figure.dpi']\n" ) subprocess.run([sys.executable, '-c', code], check=True) + + +# -------------------------------------------------------------------------- +# 2026-07-27 review +# -------------------------------------------------------------------------- + +# A move that mutates and then reports failure must not corrupt the sampler +# (bug: the ensembles trusted the sentinel to mean "atoms untouched" and +# skipped the rollback, leaving the config out of step with E_old) + +class _DirtyBailingMove(BaseMove): + """Contract violator: mutates the atoms, then returns the ``False`` + 'could not propose' sentinel.""" + + def __init__(self, seed=1): + super().__init__(NullCell(), ['X'], seed) + + def do_trial_move(self, atoms): + atoms.positions += 1.0 + return False, 0, 'X' + + +def _h2(): + return Atoms('H2', positions=[[0.0, 0.0, 0.0], [2.0, 0.0, 0.0]], + cell=[10.0, 10.0, 10.0]) + + +def test_gcmc_rolls_back_a_move_that_mutates_then_reports_failure(): + atoms = _h2() + before = atoms.positions.copy() + g = _gcmc(atoms, [NullCell()], MoveSelector([1], [_DirtyBailingMove()]), + {'H': 0.0}) + energy_before = g.E_old + g.do_gcmc_step() + np.testing.assert_allclose(g.atoms.positions, before) + assert g.E_old == energy_before + + +def test_batched_re_rolls_back_a_move_that_mutates_then_reports_failure(): + from mcpy.ensembles.batched_replica_exchange import BatchedReplicaExchange + + atoms = _h2() + before = atoms.positions.copy() + replica = _gcmc(atoms, [NullCell()], + MoveSelector([1], [_DirtyBailingMove()]), {'H': 0.0}) + + class _BatchCalc(StubCalc): + def get_potential_energies(self, atoms_list): + return np.array([self.get_potential_energy(a) for a in atoms_list]) + + pt = BatchedReplicaExchange.__new__(BatchedReplicaExchange) + pt.replicas = [replica] + pt.calculator = _BatchCalc() + pt._batched_single_move([0]) + np.testing.assert_allclose(replica.atoms.positions, before) + + +# Replica exchange swaps configurations, not slot bookkeeping (bug: get_state +# emitted step/exchange counters and set_state restored them, so accepted +# swaps traded the two ranks' tallies) + +def test_gcmc_set_state_keeps_local_step_and_exchange_counters(): + g = _gcmc(_h2(), [NullCell()], MoveSelector([1], [_DirtyBailingMove()]), + {'H': 0.0}) + g._step, g.exchange_attempts, g.exchange_successes = 40, 8, 3 + + partner = _gcmc(_h2(), [NullCell()], + MoveSelector([1], [_DirtyBailingMove()]), {'H': 0.0}) + partner._step, partner.exchange_attempts, partner.exchange_successes = 40, 8, 7 + g.set_state(partner.get_state()) + + assert (g._step, g.exchange_attempts, g.exchange_successes) == (40, 8, 3) + assert g.atoms is partner.atoms # the configuration still travels + + +def test_canonical_set_state_keeps_local_step_and_exchange_counters(tmp_path): + mc = _canonical(tmp_path) + mc._step, mc.exchange_attempts, mc.exchange_successes = 12, 4, 1 + mc.set_state({'atoms': _h2(), 'energy': -1.0, 'step': 99, + 'exchange_attempts': 40, 'exchange_successes': 39}) + assert (mc._step, mc.exchange_attempts, mc.exchange_successes) == (12, 4, 1) + assert mc._current_energy == -1.0 + + +# Per-interval acceptance counters are cleared even with the outfile disabled +# (bug: write_outfile returned before reset_counters, so interval_ratios() +# silently degraded into total_ratios()) + +def test_interval_ratios_reset_with_outfile_disabled(): + atoms = _h2() + g = _gcmc(atoms, [NullCell()], MoveSelector([1], [DisplacementMove( + species=['H'], seed=5, max_displacement=0.1)]), {'H': 0.0}) + assert g._outfile is None + g.do_gcmc_step() + assert g.move_selector.move_counter == [1] + g.write_outfile() + assert g.move_selector.move_counter == [0] + assert g.move_selector.move_counter_total == [1] # cumulative untouched + + +# Distinct outfile labels for the molecule moves (bug: a 3-character slice +# collapsed every Molecule*Move to 'Mol') + +def test_move_labels_distinguish_the_molecule_moves(): + from mcpy.moves.move_selector import _abbreviate + + assert _abbreviate('MoleculeInsertionMove') == 'MolIns' + assert _abbreviate('MoleculeDeletionMove') == 'MolDel' + assert _abbreviate('MoleculeDisplacementMove') == 'MolDis' + labels = [_abbreviate(n) for n in ('MoleculeInsertionMove', + 'MoleculeDeletionMove', + 'MoleculeDisplacementMove')] + assert len(set(labels)) == 3 + + +def test_move_labels_unchanged_for_single_word_moves(): + """Atomic runs keep the labels their existing outfiles were written with.""" + from mcpy.moves.move_selector import _abbreviate + + for name in ('InsertionMove', 'DeletionMove', 'DisplacementMove', + 'PermutationMove', 'ShakeMove', 'BrownianMove'): + assert _abbreviate(name) == name[:3] + + +# A missing exclusion radius names the species (bug: bare KeyError from inside +# the free-volume sampler) + +def test_missing_species_radius_raises_a_named_error(): + atoms = Atoms('Cu2O', positions=[[0, 0, 0], [2, 0, 0], [0, 2, 0]], + cell=[10.0, 10.0, 10.0]) + cell = CustomCell(atoms, custom_height=5.0, bottom_z=0.0, + species_radii={'Cu': 2.4}, mc_sample_points=64, seed=0) + with pytest.raises(ValueError, match=r"no radius for \['O'\]"): + cell.calculate_volume(atoms) diff --git a/tests/test_batched_re_acceptance.py b/tests/test_batched_re_acceptance.py index 8824dba..479b93d 100644 --- a/tests/test_batched_re_acceptance.py +++ b/tests/test_batched_re_acceptance.py @@ -7,6 +7,9 @@ """ import logging +import numpy as np +import pytest + from mcpy.ensembles.batched_replica_exchange import BatchedReplicaExchange @@ -85,6 +88,123 @@ def test_accept_swap_uses_grand_potential_not_bare_energy(): assert re._accept_swap(0, 1) is True +# -------------------------------------------------------------------------- +# mu-ladder swaps (bug: the temperature-only criterion collapses to delta=0 +# when both replicas share a temperature, accepting every swap unconditionally) +# -------------------------------------------------------------------------- + +_BETA_300K = 1.0 / (8.617333e-5 * 300.0) + + +def _mu_ladder(u): + """Two same-temperature replicas differing only in mu and N.""" + re = _bare_re() + ri = _FakeReplica(beta=_BETA_300K, energy=-100.0, mu={'O': -5.0}, + symbols=['O'] * 10) + rj = _FakeReplica(beta=_BETA_300K, energy=-100.0, mu={'O': -3.0}, + symbols=['O'] * 60) + re.replicas = [ri, rj] + re.rng = _FakeRng(u) + return re + + +def test_mu_ladder_swap_is_not_always_accepted(): + # beta_i == beta_j, so (beta_j - beta_i)(Phi_j - Phi_i) is identically 0 + # and the old rule returned p = 1 for every draw. The cross terms give + # exp(beta (mu_i - mu_j)(N_j - N_i)) = exp(-3868) here: always reject. + assert not _mu_ladder(u=1e-12)._accept_swap(0, 1) + + +def _partial_mu_ladder(): + """mu-ladder tuned to an acceptance well inside (0, 1), so that both the + 'just below p' and the 'just above p' draw are legal uniforms -- a p + pinned at 1 would let the buggy always-accept rule pass either way. + Energies differ too: they must cancel out of a same-temperature swap.""" + beta, mu_i, n_i, n_j = _BETA_300K, -5.0, 10, 12 + mu_j = mu_i + 0.5 / (beta * (n_j - n_i)) # -> delta = -0.5 + re = _bare_re() + re.replicas = [ + _FakeReplica(beta=beta, energy=-100.0, mu={'O': mu_i}, + symbols=['O'] * n_i), + _FakeReplica(beta=beta, energy=-42.0, mu={'O': mu_j}, + symbols=['O'] * n_j), + ] + p = float(np.exp(beta * (mu_i - mu_j) * (n_j - n_i))) + assert p == pytest.approx(np.exp(-0.5)) + return re, p + + +def test_mu_ladder_swap_matches_analytic_probability(): + re, p = _partial_mu_ladder() + re.rng = _FakeRng(p * 0.99) + assert re._accept_swap(0, 1) + re.rng = _FakeRng(p * 1.01) + assert not re._accept_swap(0, 1) + + +def test_mu_ladder_swap_ignores_the_energy_difference(): + """At one temperature the two configs' energies cancel exactly, so the + swap probability must not move when either replica's energy changes.""" + re, p = _partial_mu_ladder() + re.rng = _FakeRng(p * 0.99) + re.replicas[0].E_old += 37.0 + re.replicas[1].E_old -= 12.5 + assert re._accept_swap(0, 1) + re.rng = _FakeRng(p * 1.01) + assert not re._accept_swap(0, 1) + + +def test_mu_ladder_swap_toward_favoured_replica_always_accepted(): + # Reverse the population imbalance: the same expression is now positive, + # so the swap is downhill in the joint weight and must always be taken. + re = _bare_re() + re.replicas = [ + _FakeReplica(beta=_BETA_300K, energy=-100.0, mu={'O': -3.0}, + symbols=['O'] * 10), + _FakeReplica(beta=_BETA_300K, energy=-100.0, mu={'O': -5.0}, + symbols=['O'] * 60), + ] + re.rng = _FakeRng(1.0 - 1e-12) + assert re._accept_swap(0, 1) # p clamped to 1 without an exp() overflow + + +def test_swap_probability_is_symmetric_in_slot_order(): + """Both orderings must agree, or the pairing loop's (i, j) choice would + silently bias the ladder.""" + re, p = _partial_mu_ladder() + for u in (p * 0.99, p * 1.01): + re.rng = _FakeRng(u) + assert re._accept_swap(0, 1) == re._accept_swap(1, 0) + + +def test_temperature_ladder_swap_unchanged_by_the_mu_fix(): + """With a shared mu the general form must reduce exactly to the old + (beta_j - beta_i)(Phi_j - Phi_i).""" + ri = _FakeReplica(beta=1.0, energy=10.0, mu={'Ag': 1.0}, symbols=['Ag'] * 5) + rj = _FakeReplica(beta=0.5, energy=20.0, mu={'Ag': 1.0}, symbols=['Ag'] * 18) + re = _bare_re() + re.replicas = [ri, rj] + # Phi_i = 5, Phi_j = 2 -> delta = (0.5 - 1.0)(2 - 5) = 1.5 -> p = 1. + legacy_delta = (rj.units.beta - ri.units.beta) * ( + re._grand_potential(rj) - re._grand_potential(ri)) + assert legacy_delta == pytest.approx(1.5) + re.rng = _FakeRng(1.0 - 1e-12) + assert re._accept_swap(0, 1) + + +def test_temperature_ladder_uphill_swap_matches_legacy_probability(): + ri = _FakeReplica(beta=1.0, energy=10.0, mu={'Ag': 1.0}, symbols=['Ag'] * 18) + rj = _FakeReplica(beta=0.5, energy=20.0, mu={'Ag': 1.0}, symbols=['Ag'] * 5) + re = _bare_re() + re.replicas = [ri, rj] + # Phi_i = -8, Phi_j = 15 -> delta = (0.5 - 1.0)(15 - -8) = -11.5. + p = float(np.exp(-11.5)) + re.rng = _FakeRng(p * 0.99) + assert re._accept_swap(0, 1) + re.rng = _FakeRng(p * 1.01) + assert not re._accept_swap(0, 1) + + def test_consolidated_status_line(caplog): """One console line covers all replicas; per-replica detail stays in files.""" import logging @@ -107,3 +227,58 @@ def test_consolidated_status_line(caplog): assert len(caplog.records) == 1 msg = caplog.records[0].getMessage() assert 'RE 30/100' in msg and '42' in msg and '44' in msg and '25%' in msg + + +# -------------------------------------------------------------------------- +# Degenerate-ladder warning. A 100%-accept ladder sat unremarked in the +# "Accepted Exchange (%)" column of six consecutive production runs before the +# mu-ladder bug above was found; an all-accept tally must not stay silent. +# -------------------------------------------------------------------------- + +def _tally_re(attempts, successes, mus=None): + pt = BatchedReplicaExchange.__new__(BatchedReplicaExchange) + pt.logger = logging.getLogger('mcpy.ensembles.batched_replica_exchange') + pt.exchange_attempts = attempts + pt.exchange_successes = successes + pt.mus = mus + return pt + + +@pytest.mark.parametrize('successes, expect', [ + ([40, 40], True), # every swap accepted + ([0, 0], True), # no swap ever accepted + ([12, 12], False), # healthy partial acceptance +]) +def test_degenerate_ladder_warning(caplog, successes, expect): + pt = _tally_re([40, 40], successes) + with caplog.at_level(logging.WARNING, + logger='mcpy.ensembles.batched_replica_exchange'): + pt._warn_if_ladder_degenerate() + warned = any(r.levelno == logging.WARNING for r in caplog.records) + assert warned is expect + + +def test_degenerate_ladder_warning_needs_enough_attempts(caplog): + """Two all-accept swaps prove nothing; don't cry wolf.""" + pt = _tally_re([1, 1], [1, 1]) + with caplog.at_level(logging.WARNING, + logger='mcpy.ensembles.batched_replica_exchange'): + pt._warn_if_ladder_degenerate() + assert not caplog.records + + +def test_degenerate_ladder_warning_names_the_ladder_kind(caplog): + pt = _tally_re([40, 40], [40, 40], mus=[{'CO': -1.0}, {'CO': -0.5}]) + with caplog.at_level(logging.WARNING, + logger='mcpy.ensembles.batched_replica_exchange'): + pt._warn_if_ladder_degenerate() + text = ' '.join(r.getMessage() for r in caplog.records) + assert '80/80' in text and 'mu spacing' in text + + caplog.clear() + pt = _tally_re([40, 40], [40, 40], mus=None) + with caplog.at_level(logging.WARNING, + logger='mcpy.ensembles.batched_replica_exchange'): + pt._warn_if_ladder_degenerate() + assert 'temperature spacing' in ' '.join( + r.getMessage() for r in caplog.records) diff --git a/tests/test_compound_moves.py b/tests/test_compound_moves.py index 44a9f44..e1447ce 100644 --- a/tests/test_compound_moves.py +++ b/tests/test_compound_moves.py @@ -78,3 +78,46 @@ def test_displacement_n_steps_exceeds_movable_raises(): move = DisplacementMove(species=['Au'], seed=1, n_steps=5) with pytest.raises(ValueError): move.do_trial_move(atoms) + + +# -------------------------------------------------------------------------- +# A trial must be all-or-nothing (bug: a multi-swap trial could apply k swaps +# and then report the "couldn't propose" sentinel, which the ensembles read as +# "the atoms were not touched" -- so the config changed while E_old did not) +# -------------------------------------------------------------------------- + +def test_permutation_absent_species_does_not_abort_a_multi_swap_trial(): + """Declaring a species that is not in the system must not turn a later + iteration into a mid-trial bail-out. Species counts are swap-invariant, so + the usable pair is resolved once, up front.""" + for seed in range(200): + atoms = _balanced_alloy() + before = atoms.get_atomic_numbers().copy() + move = PermutationMove(species=['Au', 'Pt', 'Ag'], seed=seed, n_swaps=3) + result, delta, _ = move.do_trial_move(atoms) + assert result is atoms, f'seed={seed} bailed out mid-trial' + assert delta == 0 + assert not np.array_equal(atoms.get_atomic_numbers(), before) + + +def test_permutation_without_two_present_species_is_a_clean_no_op(): + """The one case that genuinely cannot propose must leave the atoms + byte-identical, since the ensembles skip the rollback on that path.""" + atoms = Atoms('Au4', positions=[[0, 0, 0], [2, 0, 0], [0, 2, 0], [0, 0, 2]]) + before = atoms.get_atomic_numbers().copy() + move = PermutationMove(species=['Au', 'Pt'], seed=1, n_swaps=3) + result, delta, name = move.do_trial_move(atoms) + assert result is False + assert (delta, name) == (0, 'X') + np.testing.assert_array_equal(atoms.get_atomic_numbers(), before) + + +def test_permutation_stream_unchanged_when_every_species_is_present(): + """Filtering to the present species must be a no-op for a healthy setup, + or previously published runs would no longer reproduce.""" + results = [] + for species in (['Au', 'Pt'], ['Au', 'Pt']): + atoms = _balanced_alloy() + PermutationMove(species=species, seed=11, n_swaps=5).do_trial_move(atoms) + results.append(atoms.get_atomic_numbers().copy()) + np.testing.assert_array_equal(*results) diff --git a/tests/test_custom_cell_region.py b/tests/test_custom_cell_region.py index dfca434..f751e37 100644 --- a/tests/test_custom_cell_region.py +++ b/tests/test_custom_cell_region.py @@ -64,3 +64,84 @@ def test_subsurface_oxygen_is_excluded_from_deletion_candidates(): assert total_o == 2 # both O are present in the structure assert len(counted) == 1 # only the above-floor O is a deletion candidate + + +# -------------------------------------------------------------------------- +# Exchangeable region vs proposal region (bug: molecules that desorbed above +# the cell top stopped being deletion candidates AND dropped out of the +# per-species de Broglie count, driving runaway insertion) +# -------------------------------------------------------------------------- + +def _point_at(cell, frac_z): + """A point in the middle of the footprint at fractional height ``frac_z`` + (``> 1`` is above the cell top).""" + return np.array([0.5, 0.5, frac_z]) @ cell.dimensions + cell.offset + + +def test_point_above_the_cell_top_is_exchangeable_but_not_inside(): + """The two predicates are deliberately different: ``is_point_inside`` + bounds the *proposal* region, ``is_point_exchangeable`` the region the + reservoir can take molecules back from -- the same asymmetry + ``get_atoms_specie_inside_cell`` already applies to single atoms.""" + _, cell = make_cell() + escaped = _point_at(cell, 1.4) + assert not cell.is_point_inside(escaped) + assert cell.is_point_exchangeable(escaped) + + +def test_point_below_the_cell_floor_is_neither(): + _, cell = make_cell() + buried = _point_at(cell, -0.2) + assert not cell.is_point_inside(buried) + assert not cell.is_point_exchangeable(buried) + + +def test_point_outside_the_xy_footprint_is_neither(): + _, cell = make_cell() + outside = np.array([1.4, 0.5, 0.5]) @ cell.dimensions + cell.offset + assert not cell.is_point_inside(outside) + assert not cell.is_point_exchangeable(outside) + + +def test_desorbed_molecule_stays_a_deletion_candidate(): + """A CO whose center of mass drifts above the cell top must remain + findable, exactly as its individual C and O atoms already do. Otherwise it + can never be deleted and its absence from ``last_exchange_count`` inflates + V/((N+1)Lambda^3) -- the runaway insertion mode documented in + docs/gcmc_acceptance_convention.rst. + """ + from ase import Atoms + + from mcpy.moves.molecule_utils import find_molecules + + atoms, cell = make_cell() + atoms.new_array('molecule_id', np.full(len(atoms), -1, dtype=int)) + for mol_id, frac_z in enumerate((0.5, 1.4)): # inside, and above the top + frag = Atoms('CO', positions=[[0, 0, 0], [0, 0, 1.13]]) + frag.positions += _point_at(cell, frac_z) + frag.new_array('molecule_id', np.full(len(frag), mol_id, dtype=int)) + atoms += frag + + template = sorted(['C', 'O']) + assert len(find_molecules(atoms, template)) == 2 + assert len(find_molecules(atoms, template, cell)) == 2 + # The atomic path has always seen all four member atoms; the molecular + # path must not disagree with it. + assert len(cell.get_atoms_specie_inside_cell(atoms, ['C', 'O'])) == 4 + + +def test_buried_molecule_is_still_excluded(): + """The floor exclusion is intentional (buried species are kept), so the + fix above must not open that side of the region too.""" + from ase import Atoms + + from mcpy.moves.molecule_utils import find_molecules + + atoms, cell = make_cell() + atoms.new_array('molecule_id', np.full(len(atoms), -1, dtype=int)) + frag = Atoms('CO', positions=[[0, 0, 0], [0, 0, 1.13]]) + frag.positions += _point_at(cell, -0.5) + frag.new_array('molecule_id', np.zeros(len(frag), dtype=int)) + atoms += frag + + assert find_molecules(atoms, sorted(['C', 'O']), cell) == [] diff --git a/tests/test_molecule_moves.py b/tests/test_molecule_moves.py index 3eaf3dc..633a0c5 100644 --- a/tests/test_molecule_moves.py +++ b/tests/test_molecule_moves.py @@ -86,7 +86,7 @@ def test_find_molecules_missing_array_and_cell_filter(): assert find_molecules(plain, ['H', 'H']) == [] class _NowhereCell: - def is_point_inside(self, point): + def is_point_exchangeable(self, point): return False atoms = _water_box() @@ -448,13 +448,13 @@ def test_write_xyz_molecule_id_roundtrip(tmp_path): def test_find_molecules_queries_com_point(): - """The point handed to cell.is_point_inside must be the molecule's + """The point handed to cell.is_point_exchangeable must be the molecule's center of mass, not e.g. the first member's position.""" atoms = _water_box() recorded = {} class _RecordingCell: - def is_point_inside(self, point): + def is_point_exchangeable(self, point): recorded['point'] = point return True