Fix twelve findings from the 2026-08-17 review

docs/replica_exchange_ladder_spacing.rst

@@ -79,7 +79,9 @@ and neighbouring rungs exchange when their distributions overlap, i.e. when
 
    \Delta\mu(\mu) = \frac{1}{\sqrt{\beta\,dN/d\mu}} ,
    \qquad
-   n_\text{rungs} = \int \sqrt{\beta\,\frac{dN}{d\mu}}\; d\mu .
+   n_\text{gaps} = \int \sqrt{\beta\,\frac{dN}{d\mu}}\; d\mu ,
+   \qquad
+   n_\text{rungs} = n_\text{gaps} + 1 .
 
 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
@@ -118,14 +120,23 @@ is the disease being diagnosed.
        """
        beta = 1.0 / (8.617333e-5 * temperature)
        mu, n_ads = np.asarray(mu, float), np.asarray(n_ads, float)
+       order = np.argsort(mu)             # desorption scans arrive descending
+       mu, n_ads = mu[order], n_ads[order]
        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))
+       # Clip noise-induced negative dN/dmu to zero: a NaN here would
+       # propagate into ``cumulative`` and crash np.interp cryptically.
+       density = np.sqrt(beta * np.maximum(np.interp(fine, midpoints, slope),
+                                           0.0))
        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)
+       # linspace, not arange: the ladder must reach mu_max, so the last gap
+       # shrinks rather than the top rung being dropped.
+       targets = np.linspace(0.0, cumulative[-1],
+                             int(np.ceil(cumulative[-1])) + 1)
+       return np.interp(targets, cumulative, fine)
 
 
 Calibration

mcpy/calculators/_alchemi_common.py

@@ -48,6 +48,15 @@ def _load_model(
     head: Union[str, int, None] = None,
 ) -> MACEWrapper:
     if isinstance(checkpoint, MACEWrapper):
+        if head is not None:
+            raise ValueError(
+                'head= cannot be applied to a pre-loaded MACEWrapper: the '
+                'head is baked in at wrap time, so it would be silently '
+                'ignored and a multihead model would evaluate on head 0 '
+                '(usually the pretrain head). Load from the .model path '
+                'with head=, or wrap the raw model in _HeadMACEWrapper '
+                'yourself before sharing it.'
+            )
         return checkpoint
     if compile_model:
         # torch >= 2.12.1 donates compiled-backward buffers that the FIRE

mcpy/calculators/alchemi_calculator.py

@@ -56,7 +56,7 @@ def __init__(
     ) -> 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 —
+            # 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
@@ -65,8 +65,8 @@ def __init__(
                 "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."
+                "wrapper. Pass the checkpoint path instead, so this "
+                "calculator loads its own model."
             )
         self.device = device
         self.dtype = dtype
@@ -163,6 +163,15 @@ def run_md(
         model is loaded. ``friction`` is in 1/fs, ``dt`` in fs, ``temperature``
         in K. ``FixAtoms`` constraints are honored.
         """
+        if self.energy_only:
+            # The integrator would otherwise die steps deep inside nvalchemi
+            # ("NVTLangevin requires forces...") without naming the cause.
+            raise ValueError(
+                'run_md needs forces, but this calculator was built with '
+                'energy_only=True (forces are stripped from the model '
+                'outputs). Build a separate calculator without energy_only '
+                'for MD.'
+            )
         _run_langevin_md(
             self.model, self._nl_config, atoms,
             temperature=temperature, friction=friction, dt=dt, steps=steps,

mcpy/calculators/alchemi_f_calculator.py

@@ -97,6 +97,16 @@ def __init__(
         self._optimizer_cls = _ALCHEMI_OPTIMIZERS[optimizer]
         self.optimizer_name = optimizer
         self.model = _load_model(checkpoint, device, dtype, enable_cueq, compile_model, head)
+        if 'forces' not in self.model.model_config.active_outputs:
+            # A model reused from an energy_only AlchemiCalculator (its
+            # ``.model`` attribute) has forces stripped; FIRE would only fail
+            # steps deep inside nvalchemi without naming the cause.
+            raise ValueError(
+                "the provided model has 'forces' disabled in its active "
+                'outputs (an energy_only AlchemiCalculator strips them), but '
+                'FIRE relaxation needs forces. Load a fresh wrapper for this '
+                'calculator instead of sharing an energy-only one.'
+            )
         self._nl_config = self.model.model_config.neighbor_config
 
     def get_potential_energy(self, atoms: Atoms) -> float:

mcpy/cell/base_cell.py

@@ -39,3 +39,19 @@ def get_volume(self):
         :return: Volume of the cell.
         """
         raise NotImplementedError("This method should be implemented in subclasses.")
+
+    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 ``get_atoms_specie_inside_cell``, and the
+        predicate the molecule moves use. It is separate from
+        ``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. Defined here on the ABC so cells written
+        against the original single-predicate contract (only
+        ``is_point_inside``) keep working; cells whose two regions coincide
+        (the box, the sphere, the dome) inherit it unchanged.
+        """
+        return self.is_point_inside(point)

mcpy/cell/cell.py

@@ -95,17 +95,3 @@ 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)

mcpy/ensembles/batched_replica_exchange.py

@@ -75,6 +75,26 @@ def __init__(
                 gcmc_factory(mu=mu, rank=i) for i, mu in enumerate(self.mus)
             ]
 
+        # Batching shares one process across replicas, so nothing isolates
+        # them but the factory contract above. A factory written for the MPI
+        # class (module-level atoms/moves, safe under one-process-per-rank)
+        # silently corrupts every replica here; refuse it up front.
+        shared = [
+            attr for attr in ('atoms', 'move_selector')
+            if len({id(getattr(r, attr, None)) for r in self.replicas})
+            < self.n_replicas
+        ]
+        cell_ids = [id(c) for r in self.replicas for c in getattr(r, 'cells', [])]
+        if len(set(cell_ids)) < len(cell_ids):
+            shared.append('cells')
+        if shared:
+            raise ValueError(
+                f"replicas share {', '.join(shared)}; gcmc_factory must build "
+                'a fresh ensemble with its own atoms, cells and move_selector '
+                'per replica (see the module docstring). Shared objects mix '
+                'the replicas\' volumes, counters and RNG streams silently.'
+            )
+
         if not hasattr(calculator, 'get_potential_energies'):
             raise TypeError(
                 "calculator must implement get_potential_energies(atoms_list). "
@@ -92,6 +112,11 @@ def __init__(
 
         self.exchange_attempts = [0] * self.n_replicas
         self.exchange_successes = [0] * self.n_replicas
+        # Midpoint snapshot of the tallies, for the second-half per-rung
+        # degeneracy check (early free swaps inflate the cumulative rate
+        # forever -- docs/replica_exchange_ladder_spacing.rst).
+        self._mid_attempts = None
+        self._mid_successes = None
 
         self._outfile = outfile
         self._outfile_handle = None
@@ -120,6 +145,9 @@ def run(self) -> None:
 
         try:
             for step in range(self.gcmc_steps):
+                if step == self.gcmc_steps // 2:
+                    self._mid_attempts = list(self.exchange_attempts)
+                    self._mid_successes = list(self.exchange_successes)
                 self._batched_gcmc_step()
 
                 if step > 0 and step % self.exchange_interval == 0:
@@ -306,6 +334,12 @@ def _accept_swap(self, i: int, j: int) -> bool:
         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.
+
+        At unequal temperatures the grand-canonical weight also carries the
+        de Broglie factors Λ_s^(-3 N_s), and Λ depends on T, so swapping
+        configurations with different particle counts shifts the joint weight
+        by a further 3 Σ_s (N_X,s - N_Y,s) ln(Λ_i,s / Λ_j,s). The term
+        vanishes on a μ ladder (shared T) and in LJ units (Λ = 1).
         """
         ri, rj = self.replicas[i], self.replicas[j]
         beta_i, beta_j = ri.units.beta, rj.units.beta
@@ -315,6 +349,17 @@ def _accept_swap(self, i: int, j: int) -> bool:
         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)
+        if beta_i != beta_j:
+            delta += self._de_broglie_cross_term(ri, rj)
+        if not np.isfinite(delta):
+            # A NaN would silently reject every draw from here on (the old
+            # min(1, exp(delta)) form silently *accepted*); a dead rung must
+            # not be quieter than a diverged relaxation.
+            self.logger.warning(
+                'swap %d<->%d: non-finite acceptance exponent (delta=%r); '
+                'rejecting. A NaN/inf replica energy usually means a '
+                'diverged relaxation.', i, j, delta)
+            return False
         # 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(
@@ -324,6 +369,23 @@ def _accept_swap(self, i: int, j: int) -> bool:
         )
         return self.rng.get_uniform() < p
 
+    @staticmethod
+    def _de_broglie_cross_term(ri, rj) -> float:
+        """3 Σ_s (N_X,s - N_Y,s) ln(Λ_i,s / Λ_j,s) for configs X (slot i) and
+        Y (slot j). Species are counted exactly as the replicas' own
+        acceptance counts them (molecules for molecular μ keys, atoms
+        otherwise), via ``GrandCanonicalEnsemble._species_count``. A replica
+        without μ (e.g. canonical) exchanges nothing with a reservoir and
+        contributes no Λ factors."""
+        term = 0.0
+        for specie in getattr(ri, '_mu', None) or {}:
+            dn = (ri._species_count(ri.atoms, specie)
+                  - rj._species_count(rj.atoms, specie))
+            if dn:
+                term += 3.0 * dn * np.log(
+                    ri.units.lambda_dbs[specie] / rj.units.lambda_dbs[specie])
+        return term
+
     @staticmethod
     def _grand_potential(r) -> float:
         """Φ = E - Σ_s μ_s N_s for a grand-canonical replica; bare E without μ.
@@ -400,8 +462,17 @@ def _warn_if_ladder_degenerate(self) -> None:
         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.
+
+        A ladder can also die at individual rungs while the others stay
+        alive, which the whole-run tally cannot show (those early free swaps
+        inflate it forever); ``_warn_dead_rungs`` covers that with the
+        second-half per-rung test from
+        docs/replica_exchange_ladder_spacing.rst.
         """
-        attempts, successes = sum(self.exchange_attempts), sum(self.exchange_successes)
+        # Both slots of a pair are tallied per attempt, so the slot sums are
+        # exactly twice the number of swap attempts; report attempts.
+        attempts = sum(self.exchange_attempts) // 2
+        successes = sum(self.exchange_successes) // 2
         if attempts < _SWAP_DEGENERACY_MIN_ATTEMPTS:
             return
         if successes == attempts:
@@ -412,12 +483,32 @@ def _warn_if_ladder_degenerate(self) -> None:
             reason = ('no swap was ever accepted -- rungs are too far apart '
                       'to exchange, so each replica sampled in isolation')
         else:
+            self._warn_dead_rungs()
             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)
+            'replica-exchange ladder looks degenerate: %d/%d swap attempts '
+            'accepted. Widen or tighten the %s spacing: %s',
+            successes, attempts,
+            'mu' if self.mus is not None else 'temperature', reason)
+
+    def _warn_dead_rungs(self) -> None:
+        """Per-rung second-half check: a rung whose swaps all fail after the
+        midpoint is dead even when the cumulative tally looks mixed (early
+        equilibration swaps accept freely and inflate it permanently)."""
+        if self._mid_attempts is None:
+            return
+        for k, (att0, suc0) in enumerate(
+                zip(self._mid_attempts, self._mid_successes)):
+            att = self.exchange_attempts[k] - att0
+            suc = self.exchange_successes[k] - suc0
+            # Half the whole-run floor: a rung sees roughly half the rounds.
+            if att >= _SWAP_DEGENERACY_MIN_ATTEMPTS // 2 and suc == 0:
+                self.logger.warning(
+                    'replica-exchange rung %d accepted no swaps in the second '
+                    'half of the run (0/%d): the ladder is dead there; '
+                    'respace the %s ladder around that rung '
+                    '(docs/replica_exchange_ladder_spacing.rst).',
+                    k, att, 'mu' if self.mus is not None else 'temperature')
 
     def _log_summary(self) -> None:
         """Final consolidated summary: per-replica move acceptances."""

mcpy/ensembles/canonical_ensemble.py

@@ -95,8 +95,14 @@ 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"]
-        # Step count and exchange statistics stay with the slot, not the
-        # config -- see GrandCanonicalEnsemble.set_state.
+        # Restored only when present, so the restart round-trip is lossless;
+        # ReplicaExchange strips these slot-bound keys before a swap -- see
+        # GrandCanonicalEnsemble.set_state.
+        self._step = state.get("step", self._step)
+        self.exchange_attempts = state.get(
+            "exchange_attempts", self.exchange_attempts)
+        self.exchange_successes = state.get(
+            "exchange_successes", self.exchange_successes)
 
     def _acceptance_condition(self, potential_diff: float) -> bool:
         if potential_diff <= 0:

mcpy/ensembles/grand_canonical_ensemble.py

@@ -89,6 +89,7 @@ def get_state(self) -> Dict[str, any]:
             "mu": self._mu,
             "temperature": self._temperature,
             "beta": self.units.beta,
+            "lambda_dbs": self.units.lambda_dbs,
             "step": self._step,
             "exchange_attempts": self.exchange_attempts,
             "exchange_successes": self.exchange_successes,
@@ -98,11 +99,16 @@ def set_state(self, state: Dict[str, any]) -> None:
         self.atoms = state["atoms"]
         self.E_old = state["energy"]
         self.n_atoms = state["n_atoms"]
-        # 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.
+        # Step count and exchange statistics are restored only when the dict
+        # carries them, so the get_state()/set_state() restart round-trip is
+        # lossless. ``ReplicaExchange`` strips them before applying a swap:
+        # they describe the *slot*, not the config, and trading them between
+        # ranks made the per-rank "Accepted Exchange (%)" column meaningless.
+        self._step = state.get("step", self._step)
+        self.exchange_attempts = state.get(
+            "exchange_attempts", self.exchange_attempts)
+        self.exchange_successes = state.get(
+            "exchange_successes", self.exchange_successes)
         # 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
@@ -169,20 +175,23 @@ def write_outfile(self, step: int = None, energy: float = None) -> None:
         except (OSError, AttributeError):
             self.logger.exception("Error writing to file %s", self._outfile)
 
+    def _species_count(self, atoms: Atoms, specie: str) -> int:
+        """Particle count for one μ species: molecules for molecular keys,
+        atoms otherwise. The one counting convention shared by
+        ``_minimum_score`` and the replica-exchange swap criterion."""
+        if specie in self.units.molecules:
+            template = self.units.molecules[specie]
+            return len(find_molecules(
+                atoms, sorted(template.get_chemical_symbols())))
+        return atoms.get_chemical_symbols().count(specie)
+
     def _minimum_score(self, atoms: Atoms, energy: float) -> float:
         """Grand potential Ω = E − Σ μ_i N_i. Comparing raw E across moves
         that change N is not meaningful in the grand canonical ensemble.
         Molecular species count molecules, atomic species count atoms."""
         score = energy
-        symbols = atoms.get_chemical_symbols()
         for specie, mu in self._mu.items():
-            if specie in self.units.molecules:
-                template = self.units.molecules[specie]
-                n = len(find_molecules(
-                    atoms, sorted(template.get_chemical_symbols())))
-            else:
-                n = symbols.count(specie)
-            score -= mu * n
+            score -= mu * self._species_count(atoms, specie)
         return score
 
     def _acceptance_condition(self,