Skip to content

Audit cleanup, deduplicated internals, and a runnable example set - #32

Merged
farrisric merged 20 commits into
mainfrom
research/metric-capture
Aug 24, 2026
Merged

farrisric merged 20 commits into
mainfrom
research/metric-capture

Conversation

@farrisric

Copy link
Copy Markdown
Owner

20 commits, 74 files, +2331/-5317. Three strands, in the order they were done.

Repo-wide over-engineering audit and its follow-through

docs/research/ponytail-audit.md records the audit; the commits after it act on
it. Dead code removed (mcpy/moves/go_moves.py, mcpy/utils/utils.py,
unreachable MoveSelector methods, dead parameters), the hand-written API
reference replaced by sphinx.ext.autodoc, and the docs pruned to the examples
that actually exist.

Deduplicated internals

One accept/reject step shared between the serial and batched GCMC loops, one
base class in place of the Alchemi _common module, shared cell-region tests,
shared ensemble bookkeeping, and shared LAMMPS parity scaffolding. Behaviour is
unchanged; tests/test_ag111_gcmc_golden.py pins the GPU path end to end
(opt-in via MCPY_GPU_TESTS=1, since it needs a GPU and a local checkpoint).

An example set a new user can actually run

The two existing examples both needed a GPU and the nvalchemi backend, so a new
user could run nothing without one. Now:

Example Shows Needs
00_hello_gcmc.py the GCMC loop: cell, moves, mu, output files ASE only
01_canonical_mc.py canonical MC and basin hopping ASE only
02_molecule_gcmc.py molecular insertion/deletion/displacement, molecule_id ASE only
03_gcmc_surface_mace.py production GCMC, O and Ag both exchanged GPU
04_replica_exchange_batched.py batched replica exchange on one GPU GPU

The CPU examples use Lennard-Jones and EMT deliberately. EMT has no usable
oxygen chemistry (O on Ag(111) either piles up without bound or never adsorbs,
0.1 eV apart in delta_mu), so the GCMC starters use the LJ system the LAMMPS
parity benchmark already validates, where the expected averages are known.

Two fixes fell out of writing them:

  • example 03 exchanged only O, so its lattice could never reconstruct into an
    oxide. Ag exchange is restored, one cell per species, and the mu_Ag the
    script already derived now drives moves.
  • the batched replica-exchange example built its cells, selector and ensemble
    with no seeds, which silently made a seeded run unrepeatable.

Both GPU examples now read as constants at the top instead of an argparse
block, matching the new ones. The IQTC cluster tutorial is dropped: it was
site-specific SLURM instructions that had to be kept in step with the example
flags it submitted.

Verification

  • flake8 mcpy/ examples/ tests/ clean
  • pytest tests/ 202 passed, 14 skipped
  • examples 00-02 run to completion on CPU; 03 and 04 run on a 5090 against a
    real MACE checkpoint
  • sphinx -b html docs builds with no warnings

The batch position array is float32, so assigning the float64 frozen
coordinates into it truncated them: the "restored" FixAtoms rows came back
about 5e-7 A off their pinned positions at slab-sized coordinates, in x and y,
growing with the coordinate magnitude. It did not accumulate across steps
(each restore re-truncates an already-truncated value) and it flipped no
acceptance decision in the Ag(111) fixture, but every Alchemi-relaxed frozen
substrate was slightly off its constraint.

Widen to float64 before restoring, in both write-back paths (the single
structure and the batched harvest). Frozen drift goes from 5.4e-7 A to
2.0e-15 A. Found by tests/test_ag111_gcmc_golden.py on its first run, which is
also the first test ever to execute this package.

Also folds the duplicated inline FixAtoms loop into _fixed_indices.
Fourteen end-to-end assertions on a 4x4x3 Ag(111) slab with the bottom layer
fixed, run through AlchemiFCalculator with the MACE small-density-agnesi
checkpoint. Two fixtures: single-replica GCMC, and a three-rung mu ladder
through BatchedReplicaExchange.

This is the first test to execute mcpy/calculators at all (the package sat at
0% coverage), and the first to exercise the Alchemi FIRE path, FreezeAtomsHook,
the batched swap rule and the batched forward pass.

Everything is seeded, so the accept/reject sequence is exact rather than
statistical: same seeds reproduce every discrete decision, with positions
agreeing to ~4e-6 A and energies to ~1e-4 eV (float32 GPU sums are not bitwise
reproducible, the decisions are).

Parameters are the calibrated set from the old examples/gcmc_custom_cell.py:
a=4.165, the O-Ag distance 2.068 used both as the exclusion radius and as the
region depth below the top layer, height 7, min_insert 0.5, and the -0.176 bulk
Ag correction. delta_mu_O = -0.5 is the balance point where both insertion and
deletion accept; past about -0.7 nothing adsorbs and the test would silently
stop covering insertion, which test_both_move_types_accept guards.

The ladder test asserts coverage is monotonic in mu (O = 8, 5, 2 for
delta_mu = -0.45, -0.50, -0.55). That is the assertion a degenerate ladder
fails while every other number still looks plausible: the all-accept mode
averages the rungs together and flattens exactly this gradient.

Opt-in, not CI: needs a GPU, nvalchemi and a local checkpoint, and takes about
six minutes. Gated on MCPY_GPU_TESTS rather than an import probe, because
importing a misconfigured torch raises SIGBUS, which no try/except can catch
and which would take the whole collection down.
Forty-three findings with paths and line numbers: dead code, duplication,
hand-rolled stdlib, and config nobody sets. About -2,700 lines available, no
dependency removable. Also lists what looked bloated but earns its lines, so a
later pass does not "simplify" a guard that encodes a real past bug.

Two hot spots worth naming here: the serial and batched GCMC accept/reject
loops are the same sixty-five lines kept in step by hand, and so are the single
structure and batched Alchemi relaxation paths (which have already drifted).
The suggested order puts the dead-code deletions first and the shared-helper
extractions last, behind the golden test.

Findings are listed only; nothing is applied.
Promotes the calibrated Ag(111)/O script out of the ignored scratch/ directory
into examples/gcmc_ag111.py, so the setup behind the golden test is tracked. It
carries the calibrated geometry, derives its reference chemical potentials from
the running potential instead of hardcoding them, and seeds all five RNGs
(cells, both moves, the selector and the ensemble) so a run with --seed is
actually repeatable; previously --seed reached only the two move RNGs.

Every other example moves to a local directory pending a decision on what to
ship. Note that eleven of them are still referenced by docs/ or README.md, so
those references now point at files that are not in the tree.
Deletes the five example pages that documented parked scripts and drops them
from the examples index; the four surviving pages carry their code inline and
never named a file. Retargets the live references: the installation page and
the cluster tutorial now point at the two tracked scripts, and the molecular
adsorbates page and the README keep only the notebook pointer.

Two cross-references broke with the deleted pages (docs/cells.rst pointing at
the dome walk-through, and the NVT example pointing at the replica-exchange
page); both are removed or redirected, and the Sphinx build is warning-free.

CHANGELOG.md, docs/superpowers/ and docs/research/ keep their references on
purpose: they are historical records of when those scripts existed. The
notebooks still reference three parked scripts and are left untouched, since
editing executed notebook JSON is a separate job.
Findings 1, 7, 11 and 15, all zero-behaviour-change:

- mcpy/moves/go_moves.py, whole file. Three of its classes were stale
  copy-based twins of the exported moves; the other four had no caller, no
  export and no test.
- mcpy/calculators/mace_calculator.py. No code in the repo used it:
  MACE_F_Calculator builds the upstream mace.calculators.MACECalculator, not
  this one, which hand-rolled the AtomicData and batch dicts instead.
- Six unused helpers in mcpy/utils/utils.py (find_surface_indices,
  sphere_volume, total_volume, total_volume_with_overlap, overlap_volume,
  get_volume) and the module logger they were the only users of. The analytic
  overlap volume was superseded by the cells' Monte Carlo free-volume sampling.
  get_p_at_support stays: it is a real helper, currently used only by parked
  examples.
- Duplicate get_species() in SphericalCell and CustomCell, plus the
  species_radii reassignment those two and DomeCell repeated after
  Cell.__init__ had already stored it. The audit claimed DomeCell also
  duplicated get_species; it did not.

Two tests went with them: the overlap_volume pair, which existed only to keep
that function alive, and one assertion in test_logging that required
mcpy.utils.utils to expose a module logger, which it no longer needs because
nothing in it logs.

Sphinx build is warning-free and the torch-free suite passes (200 tests).
Finding 13. Signatures that promised something the code never did:

- BaseEnsemble(units_type=, user_tag=). units_type appeared only in the
  signature and was never read; user_tag was stored and never read.
  GrandCanonicalEnsemble keeps its own units_type, which is real: it builds
  SetUnits from it. It just stops forwarding it to a base class that discards
  it.
- CanonicalEnsemble(units_type='metal'). Accepted and ignored: the class
  computes beta from ase.units.kB directly, so passing units_type='LJ' silently
  got you metal units. Removing it turns a wrong answer into a TypeError.
- RandomNumberGenerator(warm_up=0). The docstring already said the Mersenne
  Twister does not need warm-up and the default was 0; no caller set it.
- MoveSelector.get_operator and .calculate_volumes: no caller anywhere.
- MoveSelector.get_acceptance_ration, the deprecated typo'd alias, and the one
  test that existed to assert it still warned. get_acceptance_ratio stays.

Reference docs follow: units_type moves from the shared BaseEnsemble parameter
list to the GrandCanonicalEnsemble section where it belongs.

Two of these are caller-visible changes (CanonicalEnsemble's units_type and the
typo'd alias). No caller in the repo passed either.

Sphinx build warning-free, 199 tests pass.
Finding 2. The five reference pages restated by hand what the docstrings
already say, so they drifted: the BaseEnsemble signature there still listed
units_type and user_tag until the previous commit removed them, and
MACECalculator had its own section until it was deleted. Replaced the
hand-typed signature and parameter blocks with autoclass/automodule, keeping
the narrative intros that carry information no docstring holds (the two cell
predicates, the move return contract, which backends import conditionally).

688 lines of reference prose become 208 lines of directives, while the rendered
pages get larger, not smaller: full method documentation, inherited signatures
and [source] links, all generated from the code so they cannot drift again.

Enabling this needed three things in conf.py: autodoc_mock_imports for the
optional backends (importing a real torch on the docs builder is unnecessary,
and in a misconfigured env fatal, SIGBUS rather than ImportError),
autoclass_content='both' because the parameters are documented in __init__, and
napoleon for the Args: blocks.

Also fixes the SetUnits docstring, whose "Parameters:" block was malformed RST
(unindented items followed by an indented continuation) and which autodoc
surfaced as a docutils error the moment it was rendered. It is now a napoleon
Args: block, and documents temperature, which it had never mentioned.

Sphinx build is warning-free and every previously documented class still
appears on its page.
Finding X1. Both examples called configure_logging() between their imports,
which forced a `# noqa: E402` on every mcpy import below it. The call belongs
in main(): configure() is idempotent, and library modules only ever call
getLogger, so nothing needs logging configured at import time. Thirteen noqa
comments go away and the imports sort normally.

Two unrelated E402 sources, found while checking the first:

- tests/test_molecule_moves.py grew nine sectional mid-file imports. Hoisted
  into one sorted block; the sections still read fine without them.
- tests/test_phase_diagram.py must call matplotlib.use('Agg') before anything
  imports pyplot, so its E402 is correct code. Marked with noqa and a comment
  saying why, rather than "fixed".

tests/test_calculator_gcmc_equivalence.py keeps its noqa markers: they sit
below pytest.importorskip calls, which have to run before the imports they
guard.

flake8 is now clean across mcpy/, tests/ and examples/. 199 tests pass, and
both examples still respond to --help in the alchemi env.
Finding B1. The two parity scripts each carried their own run_lammps,
parse_thermo, block_stats, compare and ZeroCalc; mace_gcmc_parity.py even said
so in a comment ("shared infra (mirrors benchmark/lammps_gcmc_parity.py)").
block_stats was byte-identical, and the rest differed only by a field width, a
dropped docstring, and one environment variable.

They now live in benchmark/_parity_common.py, with the two real differences as
parameters: compare(..., width=) sizes the mean columns, since metal-units
energies need more room than LJ ones, and run_lammps(..., extra_env=) carries
the libtorch loader path the MACE deck needs. The MACE side's compare returned
a bare bool while the LJ side returned (ok, row); the shared one returns the
tuple and the MACE call sites unpack it.

The module has a demo() self-check under a __main__ guard, since parse_thermo
and block_stats are the kind of parsing and averaging that fails quietly.
Neither benchmark had any test before.
Finding E1, the audit's biggest duplication. do_gcmc_step and
BatchedReplicaExchange._batched_single_move were the same sixty-five lines,
comments included: the arrays+constraints snapshot, the "is False or is None"
sentinel with its identity-check explanation, the "returned a different Atoms
object" guard, the de Broglie count fallback, and the whole accept branch
(wrap, n_atoms, E_old, counter, cell volumes, _record_minimum). Keeping them in
step was manual, and the repo carried a paired regression test for each loop
precisely because they could drift.

Split into GrandCanonicalEnsemble._propose() and ._commit_or_rollback(), with
_restore() shared by both. The split is what the batched loop needs anyway:
every replica must have proposed before their energies can go into one batched
forward pass, and each replica's MoveSelector still holds its own pending
volume and exchange count when it is scored.

One behaviour change, deliberate: the batched loop had a defensive branch for a
move returning a non-tuple ("result if isinstance(result, tuple) else
(result, 0, None)"). The serial loop never had it and every move returns a
three-tuple, so it is gone; such a move now raises instead of being silently
read as delta_particles=0.

test_batched_rejected_deletion_restores_constraint used to hand-build a
SimpleNamespace replica, which stubbed out exactly the code the batched loop now
delegates to. It uses a real GrandCanonicalEnsemble replica instead, so it tests
the shared path rather than a mock of it.

Verified against the GPU golden test: 14/14 pass unchanged, both the
single-replica fixture and the three-rung mu ladder, so every accept/reject
decision and swap outcome is identical. 199 torch-free tests pass.
Findings A1, A2, A3, A6, A7. A5 was already applied with the float64 fix.

A1, the one that mattered: get_potential_energy re-implemented the batched
relaxation for a single structure, and the two had drifted apart. This path
called opt.run, which steps the whole batch until every graph has converged,
while the batched path retires each graph at its own first convergence. One
structure is a batch of one, so it now delegates: get_potential_energies
([atoms])[0].

A2: the bootstrap ritual (pre-allocate the tensors compute() writes into, build
the neighbor list, compute initial forces, zero the frozen rows) appeared three
times with the same three-line comment. Now _prepare_batch and _bootstrap.

A3: run_md was byte-identical in both calculators, docstring included. Now a
_MDMixin in the shared module. The energy_only guard reads through getattr, so
the F variant, which has no such attribute and rejects a forces-less model in
its constructor instead, inherits it harmlessly.

A6: _make_batch was exactly _make_multi_batch([atoms]); one function now takes
either.

A7: _per_graph_energies carried six lines of comment describing a
scatter-reduced fallback above a branch that only raises.

A4 is NOT applied: the audit called the fire2 optimizer uncalled, but its only
caller then lived in the ignored scratch/ directory. It is now the default in
examples/gcmc_ag111.py and what the golden test uses, so the registry stays.

Verified three ways: the GPU golden test passes 14/14 with unchanged goldens,
so the single-graph compacted path is equivalent to opt.run as expected;
benchmark/internal/verify_compact_parity.py passes (max |dE| 0.0006 eV against
the serial path, fresh max|F| within fmax, FixAtoms rows bit-identical); and the
199 torch-free tests pass. That guard needed its import updated for A6.
Findings E2 to E8, minus two the audit got wrong.

E2: the seven-column replica table was typed four times (header and rows, in
each replica-exchange driver), and ReplicaExchange printed its whole parameter
banner twice, once through logger.info and once to the outfile. Now re_header,
re_row, RE_RULE and re_banner. The two drivers label their first column
differently (MPI numbers ranks, the batched driver numbers replicas), so that
column is a parameter and both outputs are byte-identical to before.

E3: CanonicalEnsemble kept the current energy in three places at once:
_current_energy, atoms.info['key_value_pairs']['potential_energy'] (an ASE-GA
convention nothing else here reads), and lowest_energy beside the base class's
_best_energy. relax() now returns the energy instead of stashing it in
atoms.info, trial_step compares against _current_energy, and lowest_energy is
gone in favour of the base _best_energy. set_state no longer has to reconstruct
the info dict so a swapped configuration carries a baseline.

E4: _write_global_minimum existed twice; only picking the winner differed (MPI
gather vs local min), so the open/write/except tail moves to
base_ensemble.write_global_minimum.

E5: _write_initial_row re-typed write_outfile's format string to emit N/A
placeholders. Both call _row now.

E6: consolidate_logging reached for the replicas' logger by hardcoded dotted
path. It asks BaseEnsemble.__module__ instead, so moving the file cannot
silently break the silencing. The audit suggested dropping the feature and
logging at DEBUG, which would change what every serial run prints; not worth it
for a string.

E7: BaseEnsemble.__del__, a third cleanup path after finalize_run and __exit__,
which swallowed every exception.

E8: CanonicalEnsemble(optimizer=None, move_selector=None) were optional in the
signature and mandatory in fact, either default reaching a TypeError inside
relax() or do_mutation(). Now required. Callers all pass them by keyword.

NOT applied: the audit also called _exchange_pairs a one-caller helper. It has
a test (test_re_exchange_pairs_alternate_by_offset), which the audit missed by
grepping only mcpy/; it stays.

199 torch-free tests pass, the GPU golden test passes 14/14 with unchanged
goldens, and the Sphinx build is warning-free.
Findings 6, 9, 16 and 18.

6: the free-volume sampler's covered-point test (one cKDTree per distinct
radius, nearest atom within it) was written three times. Now Cell._covered_mask,
with positions and radii overridable for the periodic-image expansion CustomCell
needs. The cells still sample their own points: that is where they genuinely
differ, and it is what consumes their RNG, so the draw order is untouched.

9: get_atoms_specie_inside_cell repeated the same symbol-mask and str/list
normalisation in four cells, differing only in the region predicate. One base
implementation now calls a per-cell _inside_mask(atoms); the box cell answers
"everywhere" and each region cell keeps its own test, comments included.

16: the metal-unit constants were written out by hand (Planck, Boltzmann,
amu->kg, sqrt(e)) at 6 to 7 significant digits. They come from ase.units now,
which is already a hard dependency. The values shift by about 1e-6 relative,
which moves Lambda^3 by 4e-6, so the acceptance prefactor changes in its sixth
digit; the golden fixture's decisions are unaffected.

18: the weighted move pick was cumsum + searchsorted + an index clamp. It is
random.choices(cum_weights=...), which bisects the same cumulative weights
against one random() draw. Verified equivalent over 20k draws: identical
selections and identical final RNG state, so no run changes.

Three findings in this group are NOT applied, on evidence:

- 17, replacing the hand-built quaternion rotation with
  scipy.spatial.transform.Rotation.random: that takes a numpy Generator and
  rejects a python Random ("SeedSequence expects int or sequence of ints"),
  while the move draws from its own seeded RandomNumberGenerator. Applying it
  would need a second generator per move, so a single --seed would no longer
  determine the run. The twelve hand-written lines buy reproducibility.

- 10, replacing write_xyz with ase.io extxyz: ase writes a different comment
  line (Properties= and pbc= present, energy last, different float formatting),
  so every trajectory this package has ever produced would change shape. That
  is a user-visible output format change for 40 lines, and the writer has
  tests. Not worth it.

- 20, retiring NullCell by making BaseMove.cell optional: NullCell is exported,
  tested, and reads as a deliberate null object; replacing it spreads None
  checks through five moves to save 35 lines.

199 torch-free tests pass, the GPU golden test passes 14/14 with unchanged
goldens.
_alchemi_common.py exported twelve underscore-prefixed names, of which only
four were used by both calculators (_load_model, _make_batch,
_per_graph_energies, _MDMixin). Five belonged solely to AlchemiFCalculator and
one solely to AlchemiCalculator. "Common" was a location, not a concept, and
things drifted into it: a private module that two sibling modules and two
benchmark scripts all imported from.

What the two calculators actually shared was state and a template, which free
functions cannot absorb. Their __init__ bodies agreed on 19 of 29/39 lines (the
whole signature plus every assignment) and their get_potential_energies on 18
(the entire chunking loop, differing only in what happens per chunk). Both
survived the previous dedupe pass for exactly that reason.

So: one module, mcpy/calculators/alchemi.py, with _AlchemiBase holding the
shared state, the batching and neighbour-list helpers, the bootstrap, run_md,
get_potential_energy and the chunked get_potential_energies template. Each
subclass now implements only _evaluate_chunk: one forward pass, or a batched
FIRE relaxation. 726 lines across three files become 646 in one, and the 26
lines of cross-module private imports are gone.

The subclasses keep their full explicit signatures rather than *args/**kwargs:
the parameters are the public API, positional callers must keep working, and
autodoc has to be able to render them.

Two guards kept asymmetric on purpose, since a base class makes it tempting to
unify them: AlchemiCalculator refuses energy_only with a shared pre-loaded
wrapper and overrides run_md to reject a forces-less model, while
AlchemiFCalculator rejects one in its constructor instead.

Verified: GPU golden test 14/14 with unchanged goldens, verify_compact_parity
PASS (max |dE| 0.0003 eV, fresh max|F| within fmax), 199 torch-free tests,
Sphinx warning-free. The reference page now documents the base class too.

benchmark/internal/verify_compact_parity.py and verify_compile_parity.py import
the moved names; both updated, though that directory is untracked.
The two existing examples both need a GPU and the nvalchemi backend, so a
new user could not run anything without one. These three run on any CPU
with nothing but ASE installed, and cover the parts of the library the GPU
examples do not: the canonical ensemble, basin hopping, and the molecular
moves.

They deliberately use toy potentials. EMT has no usable oxygen chemistry
(O on Ag(111) either piles up without bound or never adsorbs, depending on
delta_mu by 0.1 eV), so the GCMC examples use the Lennard-Jones system that
the LAMMPS parity benchmark already validates, where the expected averages
are known. EMT is kept for the canonical example, where Cu and Au are
parameters it actually has.
`gcmc_ag111.py` and `re_gcmc_batched.py` named their system and their
implementation, not what a reader would learn from them, and they carried no
position in the reading order the new CPU examples introduced.

  gcmc_ag111.py     -> 03_gcmc_surface_mace.py
  re_gcmc_batched.py -> 04_replica_exchange_batched.py

Also drops two references to `examples/gcmc.py` and `examples/re_gcmc.py`,
which no longer exist.
…examples

The Ag(111) example exchanged only O, so the lattice it started with was the
lattice it ended with: O could decorate the surface but never reconstruct it
into an oxide. The two-species setup that examples/gcmc.py carried before it
was pruned is restored here, one cell per species (an incoming Ag needs a
full Ag-Ag distance where an incoming O needs one O-Ag bond), with the mu_Ag
the script already derived now actually driving moves. EXCHANGE_AG = False
falls back to the O-only setup the golden test pins.

Both GPU examples now read top to bottom as constants instead of an
argparse block: an example is read more often than it is invoked, and the
flags were duplicating every default twice. The cluster tutorial submits
them without flags accordingly, changing into the results directory since
outputs now land in the working directory.

Also seeds the replica-exchange cells, selector and ensemble. They were the
one unseeded path left in the examples, which silently made a --seed run
unrepeatable.
Site-specific install and SLURM instructions for one university cluster,
which every mcpy release has to keep in step with the examples it submits
even though nobody outside that cluster can follow it. Its only remaining
tutorial neighbour is an :orphan: redirect stub, so the Tutorials toctree
goes with it.
@farrisric
farrisric merged commit b43445e into main Aug 24, 2026
4 checks passed
@farrisric
farrisric deleted the research/metric-capture branch August 24, 2026 13:38
farrisric added a commit that referenced this pull request Aug 24, 2026
The GitHub Pages docs job builds with -W, and it does not install the package
either, so it failed with 37 "No module named 'mcpy'" warnings on the merge of
#32. Same cause as the Read the Docs fix in the previous commit, different
builder.

It also ran only on pushes to main, which is why no pull request could ever
catch it: the one place the docs are compiled with warnings as errors never
saw a PR. It now runs on pull requests too, with the Pages upload, the deploy
job and the Read the Docs trigger gated to non-PR events so nothing publishes
from a branch.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant