diff --git a/examples/gcmc_big_alchemi.py b/examples/gcmc_big_alchemi.py
index 7ade460..57527b4 100644
--- a/examples/gcmc_big_alchemi.py
+++ b/examples/gcmc_big_alchemi.py
@@ -39,6 +39,7 @@
 from mcpy.calculators import AlchemiFCalculator  # noqa: E402
 from mcpy.cell import DomeCell  # noqa: E402
 from mcpy.utils.utils import get_p_at_support  # noqa: E402
+from mcpy.utils import derive_mu_bulk, derive_mu_gas  # noqa: E402
 
 
 def parse_args():
@@ -47,8 +48,10 @@ def parse_args():
     p.add_argument('--support', default='Al2O3.poscar', help='POSCAR of the support')
     p.add_argument('--T', type=float, default=500.0, help='Temperature (K)')
     p.add_argument('--steps', type=int, default=300, help='Number of GCMC steps')
-    p.add_argument('--mu-Ag', type=float, default=-2.99, help='Chemical potential of Ag (eV)')
-    p.add_argument('--mu-O', type=float, default=-4.91, help='Reference mu_O (eV)')
+    p.add_argument('--mu-Ag', type=float, default=None,
+                   help='Chemical potential of Ag (eV); default: derived from the potential')
+    p.add_argument('--mu-O', type=float, default=None,
+                   help='Reference mu_O (eV); default: E(O2)/2 derived from the potential')
     p.add_argument('--delta-mu-O', type=float, default=-0.5, help='Shift applied to mu_O (eV)')
     p.add_argument('--relax-steps', type=int, default=100, help='Max FIRE steps per energy eval')
     p.add_argument('--fmax', type=float, default=0.05, help='FIRE force threshold (eV/A)')
@@ -103,7 +106,9 @@ def main():
          InsertionMove(scell, species=species, min_insert=0.5, seed=seed_ins)],
     )
 
-    mus = {'Ag': args.mu_Ag, 'O': args.mu_O + args.delta_mu_O}
+    mu_ag = args.mu_Ag if args.mu_Ag is not None else derive_mu_bulk(calculator, 'Ag')
+    mu_o = args.mu_O if args.mu_O is not None else derive_mu_gas(calculator)
+    mus = {'Ag': mu_ag, 'O': mu_o + args.delta_mu_O}
 
     tag = f'{atoms.get_chemical_formula()}_dmu_{args.delta_mu_O}_{args.optimizer}'
     gcmc = GrandCanonicalEnsemble(
diff --git a/examples/gcmc_nano.py b/examples/gcmc_nano.py
index 7f3f590..5752957 100644
--- a/examples/gcmc_nano.py
+++ b/examples/gcmc_nano.py
@@ -14,14 +14,17 @@
 from mcpy.moves.move_selector import MoveSelector  # noqa: E402
 from mcpy.ensembles.grand_canonical_ensemble import GrandCanonicalEnsemble  # noqa: E402
 from mcpy.cell import SphericalCell  # noqa: E402
+from mcpy.utils import derive_mu_bulk, derive_mu_gas  # noqa: E402
 
 
 def parse_args():
     p = argparse.ArgumentParser(description=__doc__)
     p.add_argument('--T', type=float, default=500.0, help='Temperature (K)')
     p.add_argument('--steps', type=int, default=1_000_000, help='Number of GCMC steps')
-    p.add_argument('--mu-Ag', type=float, default=-2.99, help='Chemical potential of Ag (eV)')
-    p.add_argument('--mu-O', type=float, default=-4.91, help='Reference mu_O (eV)')
+    p.add_argument('--mu-Ag', type=float, default=None,
+                   help='Chemical potential of Ag (eV); default: derived from the potential')
+    p.add_argument('--mu-O', type=float, default=None,
+                   help='Reference mu_O (eV); default: E(O2)/2 derived from the potential')
     p.add_argument('--delta-mu-O', type=float, default=-0.5, help='Shift applied to mu_O (eV)')
     p.add_argument('--seed', type=int, default=None, help='Master seed (random if unset)')
     p.add_argument('--device', default='cuda', help='Torch device for MACE')
@@ -50,7 +53,9 @@ def main():
          InsertionMove(scell, species=species, min_insert=0.5, seed=seed_ins)],
     )
 
-    mus = {'Ag': args.mu_Ag, 'O': args.mu_O + args.delta_mu_O}
+    mu_ag = args.mu_Ag if args.mu_Ag is not None else derive_mu_bulk(calculator, 'Ag')
+    mu_o = args.mu_O if args.mu_O is not None else derive_mu_gas(calculator)
+    mus = {'Ag': mu_ag, 'O': mu_o + args.delta_mu_O}
 
     tag = f'{atoms.get_chemical_formula()}_dmu_{args.delta_mu_O}'
     gcmc = GrandCanonicalEnsemble(
diff --git a/examples/gcmc_nano_alchemi.py b/examples/gcmc_nano_alchemi.py
index 3276860..8b61cf5 100644
--- a/examples/gcmc_nano_alchemi.py
+++ b/examples/gcmc_nano_alchemi.py
@@ -36,6 +36,7 @@
 from mcpy.ensembles.grand_canonical_ensemble import GrandCanonicalEnsemble  # noqa: E402
 from mcpy.calculators import AlchemiCalculator, AlchemiFCalculator  # noqa: E402
 from mcpy.cell import SphericalCell  # noqa: E402
+from mcpy.utils import derive_mu_bulk, derive_mu_gas  # noqa: E402
 
 
 def parse_args():
@@ -57,8 +58,10 @@ def parse_args():
                    help='Max FIRE steps per energy evaluation [relax mode only]')
     p.add_argument('--T', type=float, default=500.0, help='Temperature (K)')
     p.add_argument('--steps', type=int, default=1_000_000, help='Number of GCMC steps')
-    p.add_argument('--mu-Ag', type=float, default=-2.99, help='Chemical potential of Ag (eV)')
-    p.add_argument('--mu-O', type=float, default=-4.91, help='Reference mu_O (eV)')
+    p.add_argument('--mu-Ag', type=float, default=None,
+                   help='Chemical potential of Ag (eV); default: derived from the potential')
+    p.add_argument('--mu-O', type=float, default=None,
+                   help='Reference mu_O (eV); default: E(O2)/2 derived from the potential')
     p.add_argument('--delta-mu-O', type=float, default=-0.5, help='Shift applied to mu_O (eV)')
     p.add_argument('--min-insert', type=float, default=0.5,
                    help='Min insertion distance to existing atoms (Å)')
@@ -112,7 +115,9 @@ def main():
          InsertionMove(scell, species=species, min_insert=args.min_insert, seed=seed_ins)],
     )
 
-    mus = {'Ag': args.mu_Ag, 'O': args.mu_O + args.delta_mu_O}
+    mu_ag = args.mu_Ag if args.mu_Ag is not None else derive_mu_bulk(calculator, 'Ag')
+    mu_o = args.mu_O if args.mu_O is not None else derive_mu_gas(calculator)
+    mus = {'Ag': mu_ag, 'O': mu_o + args.delta_mu_O}
 
     mode = 'nrelax' if args.no_relax else f'fire_fmax{args.fmax}'
     tag = f'{atoms.get_chemical_formula()}_dmu_{args.delta_mu_O}_{mode}'
diff --git a/examples/gcmc_nano_supported.py b/examples/gcmc_nano_supported.py
index 53fa1db..947f260 100644
--- a/examples/gcmc_nano_supported.py
+++ b/examples/gcmc_nano_supported.py
@@ -17,6 +17,7 @@
 from mcpy.ensembles.grand_canonical_ensemble import GrandCanonicalEnsemble  # noqa: E402
 from mcpy.cell import CustomCell  # noqa: E402
 from mcpy.utils.utils import get_p_at_support  # noqa: E402
+from mcpy.utils import derive_mu_bulk, derive_mu_gas  # noqa: E402
 
 
 def parse_args():
@@ -24,8 +25,10 @@ def parse_args():
     p.add_argument('--support', default='Al2O3.poscar', help='POSCAR of the support')
     p.add_argument('--T', type=float, default=500.0, help='Temperature (K)')
     p.add_argument('--steps', type=int, default=1_000_000, help='Number of GCMC steps')
-    p.add_argument('--mu-Ag', type=float, default=-2.99, help='Chemical potential of Ag (eV)')
-    p.add_argument('--mu-O', type=float, default=-4.91, help='Reference mu_O (eV)')
+    p.add_argument('--mu-Ag', type=float, default=None,
+                   help='Chemical potential of Ag (eV); default: derived from the potential')
+    p.add_argument('--mu-O', type=float, default=None,
+                   help='Reference mu_O (eV); default: E(O2)/2 derived from the potential')
     p.add_argument('--delta-mu-O', type=float, default=-0.5, help='Shift applied to mu_O (eV)')
     p.add_argument('--seed', type=int, default=None, help='Master seed (random if unset)')
     p.add_argument('--device', default='cuda', help='Torch device for MACE')
@@ -63,7 +66,9 @@ def main():
          InsertionMove(scell, species=species, min_insert=0.5, seed=seed_ins)],
     )
 
-    mus = {'Ag': args.mu_Ag, 'O': args.mu_O + args.delta_mu_O}
+    mu_ag = args.mu_Ag if args.mu_Ag is not None else derive_mu_bulk(calculator, 'Ag')
+    mu_o = args.mu_O if args.mu_O is not None else derive_mu_gas(calculator)
+    mus = {'Ag': mu_ag, 'O': mu_o + args.delta_mu_O}
 
     tag = f'{atoms.get_chemical_formula()}_dmu_{args.delta_mu_O}'
     gcmc = GrandCanonicalEnsemble(
diff --git a/examples/re_gcmc.py b/examples/re_gcmc.py
index c55a129..c1966d9 100644
--- a/examples/re_gcmc.py
+++ b/examples/re_gcmc.py
@@ -21,6 +21,7 @@
 from mcpy.ensembles.grand_canonical_ensemble import GrandCanonicalEnsemble  # noqa: E402
 from mcpy.cell import CustomCell as Cell  # noqa: E402
 from mcpy.ensembles import ReplicaExchange  # noqa: E402
+from mcpy.utils import derive_mu_bulk, derive_mu_gas  # noqa: E402
 
 
 def parse_args():
@@ -31,8 +32,10 @@ def parse_args():
     p.add_argument('--gcmc-steps', type=int, default=200, help='GCMC steps between exchanges')
     p.add_argument('--exchange-interval', type=int, default=10, help='Steps between exchanges')
     p.add_argument('--write-interval', type=int, default=1, help='Outfile/traj write interval')
-    p.add_argument('--mu-Ag', type=float, default=-2.99, help='Chemical potential of Ag (eV)')
-    p.add_argument('--mu-O', type=float, default=-4.91, help='Reference mu_O (eV)')
+    p.add_argument('--mu-Ag', type=float, default=None,
+                   help='Chemical potential of Ag (eV); default: derived from the potential')
+    p.add_argument('--mu-O', type=float, default=None,
+                   help='Reference mu_O (eV); default: E(O2)/2 derived from the potential')
     p.add_argument('--delta-mu-O', type=float, default=-0.5, help='Shift applied to mu_O (eV)')
     p.add_argument('--rel-max-steps', type=int, default=40, help='LBFGS relaxation steps')
     p.add_argument('--rel-fmax', type=float, default=0.1, help='LBFGS force convergence (eV/Å)')
@@ -65,7 +68,9 @@ def main():
                                 steps=args.rel_max_steps, fmax=args.rel_fmax)
 
     species = ['Ag', 'O']
-    mus = {'Ag': args.mu_Ag, 'O': args.mu_O + args.delta_mu_O}
+    mu_ag = args.mu_Ag if args.mu_Ag is not None else derive_mu_bulk(calculator, 'Ag')
+    mu_o = args.mu_O if args.mu_O is not None else derive_mu_gas(calculator)
+    mus = {'Ag': mu_ag, 'O': mu_o + args.delta_mu_O}
 
     move_selector = MoveSelector(
         [25, 25, 25, 25],
diff --git a/examples/re_gcmc_batched.py b/examples/re_gcmc_batched.py
index 6b3b951..7dc4ae5 100644
--- a/examples/re_gcmc_batched.py
+++ b/examples/re_gcmc_batched.py
@@ -1,7 +1,7 @@
 """Batched replica-exchange GCMC on a single GPU.
 
 Single-process variant of ``re_gcmc.py``. All replicas live in one Python
-process and share one ``AlchemiCalculator`` instance; energies for trial
+process and share one ``AlchemiFCalculator`` instance; energies for trial
 moves are evaluated in a single batched forward pass per MC step.
 
 Run::
@@ -25,27 +25,30 @@
 from mcpy.moves import DeletionMove, InsertionMove  # noqa: E402
 from mcpy.moves.move_selector import MoveSelector  # noqa: E402
 from mcpy.ensembles.grand_canonical_ensemble import GrandCanonicalEnsemble  # noqa: E402
-from mcpy.calculators import AlchemiCalculator  # noqa: E402
+from mcpy.calculators import AlchemiFCalculator  # noqa: E402
 from mcpy.cell import CustomCell as Cell  # noqa: E402
 from mcpy.ensembles import BatchedReplicaExchange  # noqa: E402
+from mcpy.utils import derive_mu_bulk, derive_mu_gas  # noqa: E402
 
 
 def parse_args():
     p = argparse.ArgumentParser(description=__doc__,
                                 formatter_class=argparse.RawDescriptionHelpFormatter)
     p.add_argument('--temperatures', type=float, nargs='+',
-                   default=[250, 300, 350, 400, 450, 500],
+                   default=[300, 400, 500, 600, 700, 800],
                    help='Replica temperatures (K)')
     p.add_argument('--gcmc-steps', type=int, default=200)
     p.add_argument('--exchange-interval', type=int, default=10)
     p.add_argument('--write-interval', type=int, default=1)
-    p.add_argument('--checkpoint', default='medium-mpa-0')
+    p.add_argument('--checkpoint', default='mace-small-density-agnesi-stress.model')
     p.add_argument('--device', default='cuda')
     p.add_argument('--no-cueq', action='store_true')
     p.add_argument('--no-compile', action='store_true')
-    p.add_argument('--mu-Ag', type=float, default=-2.99)
-    p.add_argument('--mu-O', type=float, default=-4.91)
-    p.add_argument('--delta-mu-O', type=float, default=-0.5)
+    p.add_argument('--mu-Ag', type=float, default=None,
+                   help='Chemical potential of Ag (eV); default: derived from the potential')
+    p.add_argument('--mu-O', type=float, default=None,
+                   help='Reference mu_O (eV); default: E(O2)/2 derived from the potential')
+    p.add_argument('--delta-mu-O', type=float, default=-0.3)
     p.add_argument('--seed', type=int, default=None)
     p.add_argument('--outdir', default='.')
     return p.parse_args()
@@ -68,7 +71,7 @@ def main():
     base_atoms.set_constraint(FixAtoms(indices=bottom_layer))
 
     # ONE model, shared across replicas. Batched eval reuses this calculator.
-    calculator = AlchemiCalculator(
+    calculator = AlchemiFCalculator(
         checkpoint=args.checkpoint,
         device=args.device,
         enable_cueq=not args.no_cueq,
@@ -76,20 +79,22 @@ def main():
     )
 
     species = ['Ag', 'O']
-    mus = {'Ag': args.mu_Ag, 'O': args.mu_O + args.delta_mu_O}
+    mu_ag = args.mu_Ag if args.mu_Ag is not None else derive_mu_bulk(calculator, 'Ag')
+    mu_o = args.mu_O if args.mu_O is not None else derive_mu_gas(calculator)
+    mus = {'Ag': mu_ag, 'O': mu_o + args.delta_mu_O}
 
     def gcmc_factory(T, rank):
         atoms = base_atoms.copy()
         atoms.set_constraint(FixAtoms(indices=bottom_layer))
 
-        cell_ag_ag = Cell(atoms, custom_height=5.5, bottom_z=12.8 - 2.11,
+        cell_ag_ag = Cell(atoms, custom_height=7, bottom_z=12.8 - 2.11,
                           species_radii={'Ag': 2.75, 'O': 0})
-        cell_ag_o = Cell(atoms, custom_height=5.5, bottom_z=12.8 - 2.11,
+        cell_ag_o = Cell(atoms, custom_height=7, bottom_z=12.8 - 2.11,
                          species_radii={'Ag': 2.11, 'O': 0})
 
         s = move_seeds[4 * rank:4 * (rank + 1)]
         move_selector = MoveSelector(
-            [25, 25, 25, 25],
+            [5, 5, 5, 5],
             [DeletionMove(cell_ag_ag, species=['Ag'], seed=s[0]),
              DeletionMove(cell_ag_o, species=['O'], seed=s[1]),
              InsertionMove(cell_ag_ag, species=['Ag'], min_insert=0.5, seed=s[2]),
diff --git a/examples/re_gcmc_batched_chunked.py b/examples/re_gcmc_batched_chunked.py
index c1ade2e..845aa3c 100644
--- a/examples/re_gcmc_batched_chunked.py
+++ b/examples/re_gcmc_batched_chunked.py
@@ -50,6 +50,7 @@
 from mcpy.calculators import AlchemiFCalculator  # noqa: E402
 from mcpy.cell import SphericalCell  # noqa: E402
 from mcpy.ensembles import BatchedReplicaExchange  # noqa: E402
+from mcpy.utils import derive_mu_bulk, derive_mu_gas  # noqa: E402
 
 
 def parse_args():
@@ -72,8 +73,10 @@ def parse_args():
     p.add_argument('--no-compile', action='store_true')
     p.add_argument('--fmax', type=float, default=0.05)
     p.add_argument('--relax-steps', type=int, default=500)
-    p.add_argument('--mu-Ag', type=float, default=-2.99)
-    p.add_argument('--mu-O', type=float, default=-4.91)
+    p.add_argument('--mu-Ag', type=float, default=None,
+                   help='Chemical potential of Ag (eV); default: derived from the potential')
+    p.add_argument('--mu-O', type=float, default=None,
+                   help='Reference mu_O (eV); default: E(O2)/2 derived from the potential')
     p.add_argument('--delta-mu-O', type=float, default=-0.5)
     p.add_argument('--min-insert', type=float, default=0.5)
     p.add_argument('--vacuum', type=float, default=3.0)
@@ -111,7 +114,9 @@ def main():
     )
 
     species = ['O']
-    mus = {'Ag': args.mu_Ag, 'O': args.mu_O + args.delta_mu_O}
+    mu_ag = args.mu_Ag if args.mu_Ag is not None else derive_mu_bulk(calculator, 'Ag')
+    mu_o = args.mu_O if args.mu_O is not None else derive_mu_gas(calculator)
+    mus = {'Ag': mu_ag, 'O': mu_o + args.delta_mu_O}
 
     def gcmc_factory(T, rank):
         atoms = base_atoms.copy()
diff --git a/examples/re_gcmc_co_cupd_batched.py b/examples/re_gcmc_co_cupd_batched.py
index f444a6f..a820bd9 100644
--- a/examples/re_gcmc_co_cupd_batched.py
+++ b/examples/re_gcmc_co_cupd_batched.py
@@ -53,8 +53,8 @@ def parse_args():
                         '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)
+    p.add_argument('--rel-steps', type=int, default=100)
+    p.add_argument('--rel-fmax', type=float, default=0.05)
     p.add_argument('--min-insert', type=float, default=1.3)
     p.add_argument('--mol-disp-weight', type=int, default=0,
                    help='Weight of the rigid CO translate+rotate move '
diff --git a/examples/re_gcmc_nano_alchemi.py b/examples/re_gcmc_nano_alchemi.py
index 98acb35..c4e43d9 100644
--- a/examples/re_gcmc_nano_alchemi.py
+++ b/examples/re_gcmc_nano_alchemi.py
@@ -38,6 +38,7 @@
 from mcpy.calculators import AlchemiFCalculator  # noqa: E402
 from mcpy.cell import SphericalCell  # noqa: E402
 from mcpy.ensembles import BatchedReplicaExchange  # noqa: E402
+from mcpy.utils import derive_mu_bulk, derive_mu_gas  # noqa: E402
 
 
 def parse_args():
@@ -60,8 +61,10 @@ def parse_args():
                    help='FIRE force convergence threshold (eV/Å)')
     p.add_argument('--relax-steps', type=int, default=500,
                    help='Max FIRE steps per batched energy evaluation')
-    p.add_argument('--mu-Ag', type=float, default=-2.99, help='Chemical potential of Ag (eV)')
-    p.add_argument('--mu-O', type=float, default=-4.91, help='Reference mu_O (eV)')
+    p.add_argument('--mu-Ag', type=float, default=None,
+                   help='Chemical potential of Ag (eV); default: derived from the potential')
+    p.add_argument('--mu-O', type=float, default=None,
+                   help='Reference mu_O (eV); default: E(O2)/2 derived from the potential')
     p.add_argument('--delta-mu-O', type=float, default=-0.5, help='Shift applied to mu_O (eV)')
     p.add_argument('--min-insert', type=float, default=0.5,
                    help='Min insertion distance to existing atoms (Å)')
@@ -83,7 +86,7 @@ def main():
     master_seed = int(seeds[-1])
 
     # Shared template -- each factory call copies it and builds its own cell/moves.
-    base_atoms = Octahedron('Ag', 6, 1)
+    base_atoms = Octahedron('Ag', 9, 3)
 
     # ONE model, shared across replicas. Batched FIRE relax reuses this calculator.
     calculator = AlchemiFCalculator(
@@ -96,7 +99,9 @@ def main():
     )
 
     species = ['O']
-    mus = {'Ag': args.mu_Ag, 'O': args.mu_O + args.delta_mu_O}
+    mu_ag = args.mu_Ag if args.mu_Ag is not None else derive_mu_bulk(calculator, 'Ag')
+    mu_o = args.mu_O if args.mu_O is not None else derive_mu_gas(calculator)
+    mus = {'Ag': mu_ag, 'O': mu_o + args.delta_mu_O}
 
     def gcmc_factory(T, rank):
         atoms = base_atoms.copy()
@@ -106,7 +111,7 @@ def gcmc_factory(T, rank):
 
         s = move_seeds[2 * rank:2 * (rank + 1)]
         move_selector = MoveSelector(
-            [1, 1],
+            [25, 25],
             [DeletionMove(scell, species=species, seed=s[0]),
              InsertionMove(scell, species=species, min_insert=args.min_insert, seed=s[1])],
         )
diff --git a/mcpy/calculators/_alchemi_common.py b/mcpy/calculators/_alchemi_common.py
index 0002d1f..b11f7c3 100644
--- a/mcpy/calculators/_alchemi_common.py
+++ b/mcpy/calculators/_alchemi_common.py
@@ -16,7 +16,7 @@
 from nvalchemi.hooks import DynamicsContext
 from nvalchemi.dynamics import NVTLangevin, initialize_velocities
 from nvalchemi.dynamics.base import DynamicsStage
-from nvalchemi.dynamics.hooks import FreezeAtomsHook
+from nvalchemi.dynamics.hooks import FreezeAtomsHook, NaNDetectorHook
 
 
 class _HeadMACEWrapper(MACEWrapper):
@@ -58,12 +58,11 @@ def _load_model(
                 'yourself before sharing it.'
             )
         return checkpoint
-    if compile_model:
-        # torch >= 2.12.1 donates compiled-backward buffers that the FIRE
-        # force backward still needs; GCMC runs crash within a few steps
-        # ("compiled with non-empty donated buffers"). Zero measured cost.
-        from torch._functorch import config as _functorch_config
-        _functorch_config.donated_buffer = False
+    # The donated_buffer=False workaround (b333c8e) is gone: the crash it
+    # papered over came from the wrapper running in train mode
+    # (create_graph=True made the FIRE force backward need the buffers the
+    # compiled backward donates); with wrapper.eval() below, donation is safe
+    # (verified compile+cuEq, alias and local paths, varying-N GCMC).
     # Local .model file: load directly. MACEWrapper.from_checkpoint treats the
     # string as a download alias, so it cannot open local paths. ``head``
     # selects a multihead model's head by name or index.
@@ -96,7 +95,14 @@ def _load_model(
             for param in raw.parameters():
                 param.requires_grad = False
             raw = torch.compile(raw)
-        return MACEWrapper(raw) if idx is None else _HeadMACEWrapper(raw, idx)
+        wrapper = MACEWrapper(raw) if idx is None else _HeadMACEWrapper(raw, idx)
+        # MACEWrapper.forward passes training=self.training down to MACE, where
+        # it becomes create_graph/retain_graph on the force autograd. nn.Module
+        # defaults to train mode, and from_checkpoint's wrapper.eval() (alias
+        # path below) never runs here — without this every force pass builds a
+        # retained second-order graph.
+        wrapper.eval()
+        return wrapper
     if head is not None:
         raise ValueError(
             "head= is only supported when loading a local .model path; "
@@ -214,7 +220,8 @@ def _run_langevin_md(
         random_seed=seed,
     )
 
-    freeze_hooks = _freeze_hook_for(batch, _fixed_indices(atoms))
+    fixed = _fixed_indices(atoms)
+    freeze_hooks = _freeze_hook_for(batch, fixed)
     nl_hook = NeighborListHook(nl_config, max_neighbors=max_neighbors)
     opt = NVTLangevin(
         model=model,
@@ -223,13 +230,20 @@ def _run_langevin_md(
         friction=friction,
         random_seed=seed,
         n_steps=steps,
-        hooks=freeze_hooks,
+        # NaN forces/energy would otherwise propagate silently into the
+        # proposal the ensemble scores; the hook raises at the offending step.
+        hooks=(freeze_hooks or []) + [NaNDetectorHook()],
     )
     opt.register_hook(nl_hook, stage=DynamicsStage.BEFORE_COMPUTE)
 
     # Bootstrap: build NL and compute initial forces before the MD loop.
     _build_nl(batch, nl_hook)
     opt.compute(batch)
+    if fixed:
+        # FreezeAtomsHook zeros frozen forces only at AFTER_POST_UPDATE; the
+        # bootstrap compute() runs outside the hook loop, so step 0 would
+        # otherwise kick and displace the frozen rows once.
+        batch.forces[fixed] = 0.0
 
     opt.run(batch)
     _write_back_positions(atoms, batch)
diff --git a/mcpy/calculators/alchemi_calculator.py b/mcpy/calculators/alchemi_calculator.py
index 79bc357..3150f4c 100644
--- a/mcpy/calculators/alchemi_calculator.py
+++ b/mcpy/calculators/alchemi_calculator.py
@@ -82,6 +82,16 @@ def __init__(
             # fp32 rounding. See docs/superpowers/specs/2026-06-09-...-design.md.
             self.model.model_config.active_outputs.discard('forces')
 
+    def _forward(self, batch):
+        """Model forward. The wrapper manages positions grad itself when forces
+        are active; with energy_only nothing needs grad, and no_grad keeps the
+        forward graph-free even when model parameters still require grad
+        (compile_model=False local loads)."""
+        if self.energy_only:
+            with torch.no_grad():
+                return self.model(batch)
+        return self.model(batch)
+
     def get_potential_energy(self, atoms: Atoms) -> float:
         """
         Single forward pass — no geometry relaxation.
@@ -97,11 +107,9 @@ def get_potential_energy(self, atoms: Atoms) -> float:
             Potential energy in eV.
         """
         batch = _make_batch(atoms, self.device, self.dtype)
-        # positions must have grad enabled — MACEWrapper always computes forces via autograd
-        batch.positions.requires_grad_(True)
         nl_hook = NeighborListHook(self._nl_config, max_neighbors=self.max_neighbors)
         _build_nl(batch, nl_hook)
-        out = self.model(batch)
+        out = self._forward(batch)
         return float(out['energy'].sum().item())
 
     def get_potential_energies(
@@ -140,10 +148,9 @@ def get_potential_energies(
         for start, stop in chunk_ranges(len(atoms_list), cs):
             chunk = atoms_list[start:stop]
             batch = _make_multi_batch(chunk, self.device, self.dtype)
-            batch.positions.requires_grad_(True)
             nl_hook = NeighborListHook(self._nl_config, max_neighbors=self.max_neighbors)
             _build_nl(batch, nl_hook)
-            model_out = self.model(batch)
+            model_out = self._forward(batch)
             out.append(_per_graph_energies(model_out['energy'], len(chunk)))
         return np.concatenate(out)
 
diff --git a/mcpy/calculators/alchemi_f_calculator.py b/mcpy/calculators/alchemi_f_calculator.py
index 3029805..ff011bf 100644
--- a/mcpy/calculators/alchemi_f_calculator.py
+++ b/mcpy/calculators/alchemi_f_calculator.py
@@ -11,6 +11,7 @@
 from nvalchemi.hooks.neighbor_list import NeighborListHook
 from nvalchemi.dynamics import FIRE as AlchemiFIRE, FIRE2 as AlchemiFIRE2, ConvergenceHook
 from nvalchemi.dynamics.base import DynamicsStage
+from nvalchemi.dynamics.hooks import NaNDetectorHook
 
 from ..utils.chunking import chunk_ranges
 from ._alchemi_common import (
@@ -56,8 +57,11 @@ class AlchemiFCalculator:
     compile_model : bool
         torch.compile. Best to pre-warm in __init__ if reusing.
     dt : float
-        FIRE initial timestep (default 1.0). Matches ASE FIRE's dtmax;
-        benchmarks show dt=1.0 converges in ~half the steps vs dt=0.1.
+        FIRE initial timestep in fs (default 1.0, the value of ASE FIRE's
+        dtmax; benchmarks show dt=1.0 converges in ~half the steps vs 0.1).
+        Note nvalchemi grows the adaptive timestep up to dt_max = 10*dt,
+        unlike ASE's 1 fs cap; per-step displacement stays clamped by
+        maxstep=0.2 A either way.
     optimizer : str
         'fire' (default, classic FIRE) or 'fire2' (Guénolé et al variant —
         typically converges in fewer steps).
@@ -131,7 +135,8 @@ def get_potential_energy(self, atoms: Atoms) -> float:
         batch.forces = torch.zeros_like(batch.positions)
         batch.energy = torch.zeros(1, 1, device=self.device, dtype=self.dtype)
 
-        freeze_hooks = _freeze_hook_for(batch, _fixed_indices(atoms))
+        fixed = _fixed_indices(atoms)
+        freeze_hooks = _freeze_hook_for(batch, fixed)
 
         nl_hook = NeighborListHook(self._nl_config, max_neighbors=self.max_neighbors)
         opt = self._optimizer_cls(
@@ -139,17 +144,30 @@ def get_potential_energy(self, atoms: Atoms) -> float:
             dt=self.dt,
             convergence_hook=ConvergenceHook.from_fmax(self.fmax),
             n_steps=self.steps,
-            hooks=freeze_hooks,
+            # NaN forces/energy would otherwise flow silently into the GCMC
+            # acceptance decision; the hook raises at the offending step.
+            hooks=(freeze_hooks or []) + [NaNDetectorHook()],
         )
         opt.register_hook(nl_hook, stage=DynamicsStage.BEFORE_COMPUTE)
 
         # Bootstrap: build NL and compute initial forces before the FIRE loop
         _build_nl(batch, nl_hook)
         opt.compute(batch)
+        if fixed:
+            # FreezeAtomsHook zeros frozen forces only at AFTER_POST_UPDATE;
+            # the bootstrap compute() runs outside the hook loop, so step 0
+            # would otherwise displace the frozen rows once.
+            batch.forces[fixed] = 0.0
 
         opt.run(batch)
         self.last_relax_steps = int(opt.step_count)
         self.total_relax_steps += self.last_relax_steps
+        if self.last_relax_steps >= self.steps:
+            # Known bias trap: a too-tight cap returns truncated relaxations
+            # whose energies skew the GCMC acceptance.
+            logger.warning('FIRE reached the %d-step cap (fmax=%.3g eV/A); '
+                           'returned energy may be unconverged',
+                           self.steps, self.fmax)
         logger.debug("FIRE relaxation: %d/%d steps (fmax=%.3g eV/A)",
                      self.last_relax_steps, self.steps, self.fmax)
         _write_back_positions(atoms, batch)
@@ -177,12 +195,14 @@ def _relax_batch(self, atoms_list: List[Atoms]) -> tuple:
             dt=self.dt,
             convergence_hook=ConvergenceHook.from_fmax(self.fmax),
             n_steps=self.steps,
-            hooks=freeze_hooks,
+            hooks=(freeze_hooks or []) + [NaNDetectorHook()],
         )
         opt.register_hook(nl_hook, stage=DynamicsStage.BEFORE_COMPUTE)
 
         _build_nl(batch, nl_hook)
         opt.compute(batch)
+        if fixed:
+            batch.forces[fixed] = 0.0  # see get_potential_energy bootstrap
 
         return self._run_compacted(opt, batch, atoms_list)
 
@@ -224,6 +244,10 @@ def _run_compacted(self, opt, batch, atoms_list: List[Atoms]) -> tuple:
                     # or the next hook context builds an out-of-bounds mask.
                     opt._last_converged = None
             if alive.size:  # step cap hit: harvest the stragglers as-is
+                logger.warning('%d/%d graphs reached the %d-step FIRE cap '
+                               '(fmax=%.3g eV/A); their energies may be '
+                               'unconverged', alive.size, len(atoms_list),
+                               self.steps, self.fmax)
                 self._harvest(batch, np.arange(alive.size), alive, energies,
                               atoms_list)
         finally:
diff --git a/mcpy/cell/base_cell.py b/mcpy/cell/base_cell.py
index a30afe6..09f74cc 100644
--- a/mcpy/cell/base_cell.py
+++ b/mcpy/cell/base_cell.py
@@ -1,5 +1,8 @@
+import logging
 from abc import ABC, abstractmethod
 
+logger = logging.getLogger(__name__)
+
 
 class BaseCell(ABC):
     """
@@ -40,6 +43,22 @@ def get_volume(self):
         """
         raise NotImplementedError("This method should be implemented in subclasses.")
 
+    def _clamp_free_volume(self, free_volume, region_volume):
+        """Floor an MC free-volume estimate at one sample point's worth of
+        the region. Zero covered-everything estimates only prove the free
+        volume is below the sampler's resolution, and a literal 0.0 turns
+        the deletion prefactor N*Lambda^3/V into inf (auto-accept) or nan
+        (silent reject) in ``SetUnits.de_broglie_deletion``."""
+        floor = region_volume / self.mc_sample_points
+        if free_volume < floor:
+            logger.warning(
+                'Estimated free volume %.3g A^3 is below the sampler '
+                'resolution (%d points over %.3g A^3); flooring to %.3g A^3. '
+                'The region is effectively full -- consider enlarging it.',
+                free_volume, self.mc_sample_points, region_volume, floor)
+            return floor
+        return free_volume
+
     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).
diff --git a/mcpy/cell/custom_cell.py b/mcpy/cell/custom_cell.py
index d796e88..1b0bf96 100644
--- a/mcpy/cell/custom_cell.py
+++ b/mcpy/cell/custom_cell.py
@@ -87,7 +87,8 @@ def calculate_volume(self, atoms) -> float:
             covered |= np.isfinite(dists)
 
         occupied_fraction = float(np.count_nonzero(covered)) / self.mc_sample_points
-        self.volume = self.cell_volume * (1.0 - occupied_fraction)
+        self.volume = self._clamp_free_volume(
+            self.cell_volume * (1.0 - occupied_fraction), self.cell_volume)
 
     def _periodic_images(self, atoms):
         """Positions (cell frame) of the atoms plus their -1/0/+1 periodic
diff --git a/mcpy/cell/dome_cell.py b/mcpy/cell/dome_cell.py
index 2a3990a..3795089 100644
--- a/mcpy/cell/dome_cell.py
+++ b/mcpy/cell/dome_cell.py
@@ -134,4 +134,5 @@ def calculate_volume(self, atoms):
             covered |= np.isfinite(dists)
 
         free_fraction = float(np.count_nonzero(~covered)) / self.mc_sample_points
-        self.volume = free_fraction * self.sphere_volume
+        self.volume = self._clamp_free_volume(
+            free_fraction * self.sphere_volume, self.sphere_volume)
diff --git a/mcpy/cell/spherical_cell.py b/mcpy/cell/spherical_cell.py
index e1bbb03..a7af832 100644
--- a/mcpy/cell/spherical_cell.py
+++ b/mcpy/cell/spherical_cell.py
@@ -138,4 +138,5 @@ def calculate_volume(self, atoms):
             covered |= np.isfinite(dists)
 
         free_fraction = float(np.count_nonzero(~covered)) / self.mc_sample_points
-        self.volume = free_fraction * self.sphere_volume
+        self.volume = self._clamp_free_volume(
+            free_fraction * self.sphere_volume, self.sphere_volume)
diff --git a/mcpy/moves/alchemi_brownian_move.py b/mcpy/moves/alchemi_brownian_move.py
index 03f5a6b..2272001 100644
--- a/mcpy/moves/alchemi_brownian_move.py
+++ b/mcpy/moves/alchemi_brownian_move.py
@@ -15,6 +15,12 @@ class AlchemiBrownianMove(BaseMove):
     This move is pure Python; the nvalchemi dependency enters only through
     ``calculator.run_md`` at call time.
 
+    Positions are not wrapped back into the box after the MD trajectory.
+    This is safe: the periodic ``Cell`` selects deletion candidates by
+    species only (no position test), MACE handles PBC through neighbor-list
+    shifts, and region cells (DomeCell, CustomCell) treat out-of-region
+    drift exactly as they do for any displacement move.
+
     Parameters
     ----------
     calculator : AlchemiCalculator | AlchemiFCalculator
diff --git a/mcpy/utils/__init__.py b/mcpy/utils/__init__.py
index 3aae1f1..40d32ee 100644
--- a/mcpy/utils/__init__.py
+++ b/mcpy/utils/__init__.py
@@ -1,8 +1,11 @@
 from .random_number_generator import RandomNumberGenerator
+from .reference_mu import derive_mu_bulk, derive_mu_gas
 
 __all__ = [
     'RandomNumberGenerator',
     'analyze_phase_diagram_results',
+    'derive_mu_bulk',
+    'derive_mu_gas',
     'plot_phase_diagram',
 ]
 
diff --git a/mcpy/utils/reference_mu.py b/mcpy/utils/reference_mu.py
new file mode 100644
index 0000000..7af488a
--- /dev/null
+++ b/mcpy/utils/reference_mu.py
@@ -0,0 +1,59 @@
+"""Derive reference chemical potentials from the potential at run time.
+
+Reference energies (E(O2)/2, bulk metal energy per atom) differ between
+checkpoints, so hardcoded mu values silently change meaning whenever the
+potential is swapped. Deriving them with the run's own calculator keeps
+``delta_mu`` the only physical knob.
+"""
+import logging
+
+import numpy as np
+from ase.build import bulk, molecule
+
+logger = logging.getLogger(__name__)
+
+
+def derive_mu_gas(calculator, species: str = 'O2', box: float = 20.0) -> float:
+    """Half the molecular gas energy, e.g. mu_O = E(O2)/2.
+
+    The molecule is placed in a ``box`` A cubic cell; a relaxing calculator
+    (BaseCalculator, AlchemiFCalculator, ...) also relaxes its geometry.
+    """
+    mol = molecule(species)
+    mol.set_cell([box, box, box])
+    mol.center()
+    mol.set_pbc(True)
+    mu = float(calculator.get_potential_energy(mol)) / len(mol)
+    logger.info('Derived mu(%s/%d) = %.4f eV from the potential',
+                species, len(mol), mu)
+    return mu
+
+
+def derive_mu_bulk(calculator, symbol: str = 'Ag',
+                   crystalstructure: str = 'fcc', a: float = 4.1592,
+                   scan: float = 0.04, n_points: int = 7) -> float:
+    """Bulk energy per atom at the potential's own lattice constant.
+
+    Scans ``n_points`` lattice constants in ``a * [1-scan, 1+scan]``,
+    fits a parabola and evaluates its minimum, so the reference does not
+    inherit the lattice constant of whatever potential ``a`` came from.
+    """
+    scales = np.linspace(1.0 - scan, 1.0 + scan, n_points)
+    energies = []
+    for s in scales:
+        atoms = bulk(symbol, crystalstructure, a=a * s, cubic=True)
+        energies.append(float(calculator.get_potential_energy(atoms)) / len(atoms))
+    coeffs = np.polyfit(scales, energies, 2)
+    s_min = -coeffs[1] / (2.0 * coeffs[0])
+    if coeffs[0] <= 0.0 or not scales[0] <= s_min <= scales[-1]:
+        logger.warning(
+            'Bulk %s energy minimum fell outside the scanned lattice '
+            'constants %.3f-%.3f A (fit minimum at scale %.3f); using the '
+            'lowest sampled point. Pass a closer ``a``.',
+            symbol, a * scales[0], a * scales[-1], s_min)
+        mu = float(np.min(energies))
+    else:
+        mu = float(np.polyval(coeffs, s_min))
+        logger.info('Derived mu(%s) = %.4f eV at a = %.4f A from the potential',
+                    symbol, mu, a * s_min)
+    return mu
diff --git a/tests/test_reference_mu.py b/tests/test_reference_mu.py
new file mode 100644
index 0000000..eb5b6a9
--- /dev/null
+++ b/tests/test_reference_mu.py
@@ -0,0 +1,50 @@
+"""
+Tests for runtime reference-mu derivation. Uses a stub calculator so the
+suite stays free of torch/mace, matching the other tests.
+
+Run with: python -m pytest tests/test_reference_mu.py -v
+"""
+import numpy as np
+import pytest
+
+from mcpy.utils import derive_mu_bulk, derive_mu_gas
+
+
+class StubGasCalculator:
+    def get_potential_energy(self, atoms):
+        assert atoms.get_chemical_formula() == 'O2'
+        assert atoms.pbc.all()
+        return -9.8
+
+
+class StubBulkCalculator:
+    """Per-atom energy parabolic in the lattice constant, minimum at a_min."""
+
+    def __init__(self, a_min, e_min):
+        self.a_min = a_min
+        self.e_min = e_min
+
+    def get_potential_energy(self, atoms):
+        a = atoms.cell.lengths()[0]  # cubic=True conventional cell
+        return len(atoms) * (self.e_min + 0.5 * (a - self.a_min) ** 2)
+
+
+def test_derive_mu_gas_is_half_molecule_energy():
+    assert derive_mu_gas(StubGasCalculator()) == pytest.approx(-4.9)
+
+
+def test_derive_mu_bulk_finds_off_guess_minimum():
+    # Minimum deliberately off the guess but inside the +-4% scan window.
+    calc = StubBulkCalculator(a_min=4.25, e_min=-2.8)
+    assert derive_mu_bulk(calc, 'Ag', a=4.1592) == pytest.approx(-2.8, abs=1e-4)
+
+
+def test_derive_mu_bulk_warns_outside_scan(caplog):
+    # Minimum far outside the scan window: fall back to lowest sample.
+    calc = StubBulkCalculator(a_min=5.5, e_min=-2.8)
+    with caplog.at_level('WARNING'):
+        mu = derive_mu_bulk(calc, 'Ag', a=4.1592)
+    scan_max = 4.1592 * 1.04
+    assert mu == pytest.approx(-2.8 + 0.5 * (scan_max - 5.5) ** 2, abs=1e-4)
+    assert any('outside the scanned' in r.message for r in caplog.records)
+    assert np.isfinite(mu)
diff --git a/tests/test_zero_free_volume.py b/tests/test_zero_free_volume.py
new file mode 100644
index 0000000..3ad2a6f
--- /dev/null
+++ b/tests/test_zero_free_volume.py
@@ -0,0 +1,63 @@
+"""
+Tests that a fully-covered region never reports a free volume of exactly
+zero. The MC estimator cannot distinguish "no free volume" from "less free
+volume than one sample point resolves", and a literal 0.0 propagates into
+``SetUnits.de_broglie_deletion`` as a division by zero: the deletion
+prefactor becomes ``inf`` (every deletion auto-accepts) or ``inf * 0 = nan``
+(every deletion silently rejects). The estimate is therefore floored at one
+sample point's worth of the region volume, with a warning.
+
+Run with: python -m pytest tests/test_zero_free_volume.py -v
+"""
+import logging
+
+import numpy as np
+import pytest
+from ase import Atoms
+from ase.build import fcc111
+
+from mcpy.cell import CustomCell, DomeCell, SphericalCell
+
+
+def test_custom_cell_full_region_floors_volume(caplog):
+    atoms = fcc111('Ag', a=4.1592, size=(4, 4, 3), periodic=True, vacuum=8)
+    # A radius this large covers every sample point in the 5.5 A region.
+    cell = CustomCell(atoms, custom_height=5.5, bottom_z=12.8 - 2.11,
+                      species_radii={'Ag': 50.0, 'O': 0},
+                      mc_sample_points=1000, seed=0)
+    with caplog.at_level(logging.WARNING):
+        cell.calculate_volume(atoms)
+    assert cell.get_volume() == pytest.approx(cell.cell_volume / 1000)
+    assert cell.get_volume() > 0.0
+    assert any('free volume' in r.message for r in caplog.records)
+
+
+def test_spherical_cell_full_region_floors_volume():
+    atoms = Atoms('Ag', positions=[[0.0, 0.0, 0.0]])
+    cell = SphericalCell(atoms, vacuum=2.0, species_radii={'Ag': 50.0},
+                         mc_sample_points=1000, seed=0)
+    cell.calculate_volume(atoms)
+    assert cell.get_volume() == pytest.approx(cell.sphere_volume / 1000)
+    assert cell.get_volume() > 0.0
+
+
+def test_dome_cell_full_region_floors_volume():
+    atoms = Atoms('Ag', positions=[[0.0, 0.0, 5.0]],
+                  cell=[10.0, 10.0, 10.0], pbc=False)
+    cell = DomeCell(atoms, particle_species='Ag', bottom_z=4.0, vacuum=2.0,
+                    species_radii={'Ag': 50.0}, mc_sample_points=1000, seed=0)
+    cell.calculate_volume(atoms)
+    assert cell.get_volume() == pytest.approx(cell.sphere_volume / 1000)
+    assert cell.get_volume() > 0.0
+
+
+def test_partial_coverage_is_not_floored():
+    """A region with real free volume must keep its MC estimate."""
+    atoms = fcc111('Ag', a=4.1592, size=(4, 4, 3), periodic=True, vacuum=8)
+    cell = CustomCell(atoms, custom_height=5.5, bottom_z=12.8 - 2.11,
+                      species_radii={'Ag': 2.11, 'O': 0},
+                      mc_sample_points=20_000, seed=0)
+    cell.calculate_volume(atoms)
+    assert cell.get_volume() > 10 * cell.cell_volume / 20_000
+    assert cell.get_volume() < cell.cell_volume
+    assert np.isfinite(cell.get_volume())