Review follow-ups: runtime reference mu, free-volume floor, article example settings

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(

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(

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}'

examples/gcmc_nano_supported.py

@@ -17,15 +17,18 @@
 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():
     p = argparse.ArgumentParser(description=__doc__)
     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(

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],

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,28 +71,30 @@ 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,
         compile_model=not args.no_compile,
     )
 
     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]),

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()

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 '

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])],
         )