Audit cleanup, deduplicated internals, and a runnable example set

README.md

@@ -197,10 +197,7 @@ gcmc = GrandCanonicalEnsemble(
     move_selector=moves)
 ```
 
-Runnable examples: `examples/gcmc_molecule_mace.py` (O2 on Ag(111); any g2 molecule via `--molecule`) and
-`examples/re_gcmc_co_cupd_batched.py` (CO on a CuPd nanoparticle with batched
-replica exchange and a phase diagram). The notebooks in `notebooks/` walk
-through the machinery step by step.
+The notebooks in `notebooks/` walk through the machinery step by step.
 
 ## Contributing
 

benchmark/_parity_common.py

@@ -0,0 +1,116 @@
+"""Shared scaffolding for the LAMMPS parity benchmarks.
+
+``lammps_gcmc_parity.py`` (LJ units) and ``mace_gcmc_parity.py`` (metal units,
+MACE) run the same comparison against LAMMPS ``fix gcmc``, so they need the same
+four things: drive LAMMPS on a generated deck, parse its thermo output, block
+average a series, and print a pass/fail line. Both had their own copy;
+``block_stats`` was byte-identical and the rest differed only by a field width,
+a dropped docstring, and one environment variable.
+"""
+import math
+import os
+import subprocess
+
+import numpy as np
+
+
+class ZeroCalc:
+    """Zero energy, for the ideal-gas stage where only the prefactor matters."""
+
+    def get_potential_energy(self, atoms):
+        return 0.0
+
+
+def run_lammps(deck, outdir, tag, lmp, extra_env=None):
+    """Write ``deck``, run ``lmp`` on it, and return the parsed thermo output.
+
+    ``extra_env`` is merged over ``os.environ`` for the child; the MACE
+    benchmark uses it to put libtorch on ``LD_LIBRARY_PATH``.
+    """
+    deckfile = os.path.join(outdir, f'in.{tag}')
+    logfile = os.path.join(outdir, f'log.{tag}')
+    with open(deckfile, 'w') as f:
+        f.write(deck)
+    env = None
+    if extra_env:
+        env = dict(os.environ)
+        env.update(extra_env)
+    subprocess.run([lmp, '-in', os.path.abspath(deckfile),
+                    '-log', os.path.abspath(logfile), '-screen', 'none'],
+                   check=True, cwd=outdir, env=env)
+    return parse_thermo(logfile)
+
+
+def parse_thermo(logfile):
+    """Collect thermo rows from every run section. Returns dict of arrays
+    keyed by lower-cased column names (step, atoms, pe, ...)."""
+    cols, rows = None, []
+    with open(logfile) as f:
+        in_section = False
+        for line in f:
+            s = line.split()
+            if not s:
+                continue
+            if s[0] == 'Step':
+                cols = [c.lower() for c in s]
+                in_section = True
+                continue
+            if in_section:
+                if s[0] == 'Loop' or line.startswith('WARNING'):
+                    in_section = False
+                    continue
+                try:
+                    rows.append([float(v) for v in s])
+                except ValueError:
+                    in_section = False
+    if cols is None or not rows:
+        raise RuntimeError(f'no thermo data parsed from {logfile}')
+    rows = [r for r in rows if len(r) == len(cols)]
+    data = np.array(rows)
+    out = {c: data[:, i] for i, c in enumerate(cols)}
+    if 'poteng' in out:  # thermo header prints "PotEng" for pe
+        out['pe'] = out['poteng']
+    return out
+
+
+def block_stats(series, burn_frac=0.4, nblocks=15):
+    """Block-averaged mean and standard error after discarding ``burn_frac``."""
+    x = np.asarray(series, dtype=float)
+    x = x[int(len(x) * burn_frac):]
+    blocks = np.array_split(x, nblocks)
+    means = np.array([b.mean() for b in blocks])
+    return means.mean(), means.std(ddof=1) / math.sqrt(nblocks)
+
+
+def compare(label, m_mcpy, s_mcpy, m_lmp, s_lmp, nsig, width=10):
+    """Print one pass/fail line and return ``(ok, row)``.
+
+    ``ok`` is ``|difference| < nsig * sigma_combined``; ``row`` is the tuple the
+    CSV writer expects. ``width`` sizes the mean columns (metal-units energies
+    need more room than LJ ones).
+    """
+    diff = abs(m_mcpy - m_lmp)
+    scomb = math.sqrt(s_mcpy ** 2 + s_lmp ** 2)
+    ok = diff < nsig * scomb
+    sigmas = diff / scomb if scomb > 0 else float('inf')
+    print(f'    {label:6s} mcpy {m_mcpy:{width}.4f} ± {s_mcpy:.4f}   '
+          f'lammps {m_lmp:{width}.4f} ± {s_lmp:.4f}   '
+          f'|d|={diff:.4f} ({sigmas:.2f} sigma)  '
+          f'{"PASS" if ok else "FAIL"}')
+    return ok, (m_mcpy, s_mcpy, m_lmp, s_lmp, diff, scomb)
+
+
+def demo():
+    """Self-check: block_stats on a known series, compare's verdict logic."""
+    mean, err = block_stats(list(range(100)) * 2, burn_frac=0.0, nblocks=10)
+    assert abs(mean - 49.5) < 1e-9, mean
+    assert err > 0
+    ok, row = compare('demo', 1.0, 0.1, 1.05, 0.1, 2)
+    assert ok and len(row) == 6
+    ok, _ = compare('demo', 1.0, 0.01, 2.0, 0.01, 2)
+    assert not ok
+    print('ok')
+
+
+if __name__ == '__main__':
+    demo()

benchmark/lammps_gcmc_parity.py

@@ -16,12 +16,12 @@
 import argparse
 import math
 import os
-import subprocess
 
 import numpy as np
 from ase import Atoms
 from ase.calculators.lj import LennardJones
 
+from _parity_common import ZeroCalc, block_stats, compare, run_lammps
 from mcpy.cell import Cell
 from mcpy.ensembles.grand_canonical_ensemble import GrandCanonicalEnsemble
 from mcpy.moves import (DeletionMove, DisplacementMove, InsertionMove,
@@ -52,57 +52,6 @@ def mu_map(mu_mcpy, lam):
 # shared infra
 # --------------------------------------------------------------------------
 
-def run_lammps(deck, outdir, tag, lmp):
-    deckfile = os.path.join(outdir, f'in.{tag}')
-    logfile = os.path.join(outdir, f'log.{tag}')
-    with open(deckfile, 'w') as f:
-        f.write(deck)
-    subprocess.run([lmp, '-in', os.path.abspath(deckfile),
-                    '-log', os.path.abspath(logfile), '-screen', 'none'],
-                   check=True, cwd=outdir)
-    return parse_thermo(logfile)
-
-
-def parse_thermo(logfile):
-    """Collect thermo rows from every run section. Returns dict of arrays
-    keyed by lower-cased column names (step, atoms, pe, ...)."""
-    cols, rows = None, []
-    with open(logfile) as f:
-        in_section = False
-        for line in f:
-            s = line.split()
-            if not s:
-                continue
-            if s[0] == 'Step':
-                cols = [c.lower() for c in s]
-                in_section = True
-                continue
-            if in_section:
-                if s[0] == 'Loop' or line.startswith('WARNING'):
-                    in_section = False
-                    continue
-                try:
-                    rows.append([float(v) for v in s])
-                except ValueError:
-                    in_section = False
-    if cols is None or not rows:
-        raise RuntimeError(f'no thermo data parsed from {logfile}')
-    rows = [r for r in rows if len(r) == len(cols)]
-    data = np.array(rows)
-    out = {c: data[:, i] for i, c in enumerate(cols)}
-    if 'poteng' in out:  # thermo header prints "PotEng" for pe
-        out['pe'] = out['poteng']
-    return out
-
-
-def block_stats(series, burn_frac=0.4, nblocks=15):
-    x = np.asarray(series, dtype=float)
-    x = x[int(len(x) * burn_frac):]
-    blocks = np.array_split(x, nblocks)
-    means = np.array([b.mean() for b in blocks])
-    return means.mean(), means.std(ddof=1) / math.sqrt(nblocks)
-
-
 def write_csv(outdir, stage, results):
     """One row per (mu, observable): means, stderrs, diff, verdict."""
     path = os.path.join(outdir, f'stage{stage}_results.csv')
@@ -145,26 +94,10 @@ def plot_isotherms(outdir, s1_results, s2_results):
     return path
 
 
-def compare(label, m_mcpy, s_mcpy, m_lmp, s_lmp, nsig):
-    diff = abs(m_mcpy - m_lmp)
-    scomb = math.sqrt(s_mcpy ** 2 + s_lmp ** 2)
-    ok = diff < nsig * scomb
-    print(f'    {label:6s} mcpy {m_mcpy:10.4f} ± {s_mcpy:.4f}   '
-          f'lammps {m_lmp:10.4f} ± {s_lmp:.4f}   '
-          f'|d|={diff:.4f} ({diff / scomb if scomb > 0 else float("inf"):.2f} sigma)  '
-          f'{"PASS" if ok else "FAIL"}')
-    return ok, (m_mcpy, s_mcpy, m_lmp, s_lmp, diff, scomb)
-
-
 # --------------------------------------------------------------------------
 # mcpy side
 # --------------------------------------------------------------------------
 
-class ZeroCalc:
-    def get_potential_energy(self, atoms):
-        return 0.0
-
-
 class LJCalc:
     """Vectorized truncated-and-shifted LJ (sigma=eps=1), minimum image in a
     cubic box of side L. Verified against ASE LennardJones by

benchmark/mace_gcmc_parity.py

@@ -22,12 +22,12 @@
 import argparse
 import math
 import os
-import subprocess
 
 import numpy as np
 from ase import Atoms
 from ase.build import fcc111, molecule
 
+from _parity_common import ZeroCalc, block_stats, compare, run_lammps
 from mcpy.cell import Cell, CustomCell
 from mcpy.ensembles.grand_canonical_ensemble import GrandCanonicalEnsemble
 from mcpy.moves import (MoleculeDeletionMove, MoleculeInsertionMove,
@@ -45,81 +45,15 @@
         '~/miniconda3/envs/mcpy/lib/python3.14/site-packages/torch/lib'))
 
 
-# --------------------------------------------------------------------------
-# shared infra (mirrors benchmark/lammps_gcmc_parity.py)
-# --------------------------------------------------------------------------
-
-def run_lammps(deck, outdir, tag, lmp):
-    deckfile = os.path.join(outdir, f'in.{tag}')
-    logfile = os.path.join(outdir, f'log.{tag}')
-    with open(deckfile, 'w') as f:
-        f.write(deck)
-    env = dict(os.environ)
-    env['LD_LIBRARY_PATH'] = TORCH_LIB + ':' + env.get('LD_LIBRARY_PATH', '')
-    subprocess.run([lmp, '-in', os.path.abspath(deckfile),
-                    '-log', os.path.abspath(logfile), '-screen', 'none'],
-                   check=True, cwd=outdir, env=env)
-    return parse_thermo(logfile)
-
-
-def parse_thermo(logfile):
-    cols, rows = None, []
-    with open(logfile) as f:
-        in_section = False
-        for line in f:
-            s = line.split()
-            if not s:
-                continue
-            if s[0] == 'Step':
-                cols = [c.lower() for c in s]
-                in_section = True
-                continue
-            if in_section:
-                if s[0] == 'Loop' or line.startswith('WARNING'):
-                    in_section = False
-                    continue
-                try:
-                    rows.append([float(v) for v in s])
-                except ValueError:
-                    in_section = False
-    if cols is None or not rows:
-        raise RuntimeError(f'no thermo data parsed from {logfile}')
-    rows = [r for r in rows if len(r) == len(cols)]
-    data = np.array(rows)
-    out = {c: data[:, i] for i, c in enumerate(cols)}
-    if 'poteng' in out:
-        out['pe'] = out['poteng']
-    return out
-
-
-def block_stats(series, burn_frac=0.4, nblocks=15):
-    x = np.asarray(series, dtype=float)
-    x = x[int(len(x) * burn_frac):]
-    blocks = np.array_split(x, nblocks)
-    means = np.array([b.mean() for b in blocks])
-    return means.mean(), means.std(ddof=1) / math.sqrt(nblocks)
-
-
-def compare(label, m_mcpy, s_mcpy, m_lmp, s_lmp, nsig):
-    diff = abs(m_mcpy - m_lmp)
-    scomb = math.sqrt(s_mcpy ** 2 + s_lmp ** 2)
-    ok = diff < nsig * scomb
-    print(f'    {label:6s} mcpy {m_mcpy:12.4f} ± {s_mcpy:.4f}   '
-          f'lammps {m_lmp:12.4f} ± {s_lmp:.4f}   '
-          f'|d|={diff:.4f} ({diff / scomb if scomb > 0 else float("inf"):.2f} sigma)  '
-          f'{"PASS" if ok else "FAIL"}')
-    return ok
+def _libtorch_env():
+    """LAMMPS pair_style mace needs libtorch on the loader path."""
+    return {'LD_LIBRARY_PATH': TORCH_LIB + ':' + os.environ.get('LD_LIBRARY_PATH', '')}
 
 
 # --------------------------------------------------------------------------
 # mcpy side
 # --------------------------------------------------------------------------
 
-class ZeroCalc:
-    def get_potential_energy(self, atoms):
-        return 0.0
-
-
 class MaceSinglePoint:
     """Bare single-point MACE energy (no relaxation), float32 on GPU."""
 
@@ -224,7 +158,7 @@ def gate(outdir, lmp, model, py_calc):
 thermo_style custom step atoms pe
 run 0
 """
-    th = run_lammps(deck, outdir, 'gate', lmp)
+    th = run_lammps(deck, outdir, 'gate', lmp, extra_env=_libtorch_env())
     e_l = th['pe'][-1]
     diff = abs(e_py - e_l)
     ok = diff < 5e-4
@@ -276,7 +210,7 @@ def stage1(outdir, lmp, steps_mcpy, lmp_steps):
 thermo_style custom step atoms pe
 run {lmp_steps}
 """
-    th = run_lammps(deck, outdir, 'stage1', lmp)
+    th = run_lammps(deck, outdir, 'stage1', lmp, extra_env=_libtorch_env())
     ml, sl = block_stats(th['atoms'] / 2.0)
     ok_lmp = abs(ml - n_exact) < 2 * sl
     print(f'  lammps <N_mol> = {ml:.3f} ± {sl:.3f}  '
@@ -338,12 +272,12 @@ def stage2(outdir, lmp, model, py_calc, steps_mcpy, lmp_attempts, deltas):
 thermo_style custom step atoms pe
 run {nsteps}
 """
-        th = run_lammps(deck, outdir, f'stage2_d{delta}', lmp)
+        th = run_lammps(deck, outdir, f'stage2_d{delta}', lmp, extra_env=_libtorch_env())
         mln, sln = block_stats((th['atoms'] - n_slab) / 2.0)
         mle, sle = block_stats(th['pe'])
 
-        okn = compare('<Nmol>', mn, sn, mln, sln, 2)
-        oke = compare('<PE>', me, se, mle, sle, 3)
+        okn, _ = compare('<Nmol>', mn, sn, mln, sln, 2, width=12)
+        oke, _ = compare('<PE>', me, se, mle, sle, 3, width=12)
         all_ok &= okn and oke
     print(f'  Stage II: {"PASS" if all_ok else "FAIL"}')
     return all_ok

docs/cells.rst

@@ -166,8 +166,6 @@ Key parameters:
    )
    dome_cell.calculate_volume(atoms)
 
-A worked end-to-end run is in :doc:`examples/gcmc_dome_supported`.
-
 
 Choosing the right cell
 -----------------------

docs/conf.py

@@ -20,7 +20,17 @@
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
 
-extensions = ['sphinx.ext.autodoc', 'sphinx.ext.todo', 'sphinx.ext.viewcode']
+extensions = ['sphinx.ext.autodoc', 'sphinx.ext.napoleon',
+              'sphinx.ext.todo', 'sphinx.ext.viewcode']
+
+# The optional backends are not installed on the docs builder, and importing a
+# real torch here is both slow and, in a misconfigured env, fatal (SIGBUS).
+# Mocking lets autodoc read the wrappers' signatures and docstrings anyway.
+autodoc_mock_imports = ['torch', 'mace', 'nvalchemi', 'cuequivariance', 'mpi4py']
+autodoc_member_order = 'bysource'
+# Parameters are documented in __init__, so render both docstrings.
+autoclass_content = 'both'
+autodoc_typehints = 'description'
 
 templates_path = ['_templates']
 exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']