2026-10-03: the MBE step, phase 1 (the seamm_mbe library)#

The design is ~/Work/SEAMM/MBE_correction_step_design.rst (workspace root). Its “Decisions (2026-10-03)” section is binding, and its “Phases” section is the plan. The science context is the ~/Sites lab notebook record mlff-training/2026-09-27_vasp-periodic-route and the ~/Sites memory vasp-periodic-prototype. The design session (“design”) reviews each phase before its PR.

Phase 1 scope#

Phase 1 is a pure-Python library (numpy, and molsystem for the adapter only), extracted from the prototype’s gen_stage2_frame.py and assemble_frame.py (TinkerCliffs /projects/seamm/psaxe/periodic/tools/). It covers:

  • molecule typing by formula and topology, with a catalog for the campaign’s molecules, per-type charges and the designated atoms;

  • fragment enumeration under PBC, with canonical keys and the prototype’s names, general in the order;

  • the mixed low-level assignment;

  • the increment algebra for energy, forces and the origin-independent virial;

  • the assembled labels, with energy offsets, atomic and molecular pressure, per-body breakdown and QC;

  • the unit and sign conventions.

Phase 2 (mbe_step) and phase 3 (the adapters) need Paul’s OK.

The oracle#

The regression targets in the design document are those of the pilot frame (frames/pilot_opls0997 on TinkerCliffs). The original assemble_frame.py reproduces them exactly on the Mac from the pilot’s own outputs. The oracle Python was the seamm-psi4 conda env, which has dftd4 4.2.0.

The notebook’s data (~/Sites/.../data/vasp_results.json.gz and stage2_orca.tar.gz) are the earlier stage-2/2h runs of the same geometry. They give P_mol +1796.35 and P_MBE +831.78 atm and are not the oracle. The “science” session corrected the design document’s wording accordingly.

The fixture (tests/data/water64_pilot.npz, 569 KB) and its provenance are described in tests/data/README.md.

Results#

  • Enumeration is identical to the prototype’s frags.json: names in the same order, molecules, images, coordinates (largest difference 6e-17 Å), which pairs are outer pairs, which fragments go to VASP, and each triple’s sub-pairs. The counts are 64 monomers, 391 pairs < 4.5 Å, 268 outer pairs and 559 connected triples.

  • Calculations: 218 periodic fragments + the cell = 219 VASP inputs; high + molecular on all 1282 fragments = 2564 ORCA inputs. The molecular level is the closure of the molecular-assigned fragments, which is everything for this frame.

  • Level assignment: periodic 64 monomers + 154 pairs; molecular 237 pairs + 559 triples.

  • Numbers, all to rounding error (≤ 1e-9 relative):

    • P_MBE +832.381 atm, P_D4(r2SCAN) −2322.500 atm;

    • per-body P −1595.82 / +3062.44 / −634.24 atm;

    • per-body ΔE +5.961 (2-body) and −1.194 (3-body) kJ/mol per water;

    • REF_energy −203.32968869 eV, largest increment net force 1.2186 meV/Å.

  • One deliberate difference: the prototype converted VASP’s kB to atm with 986.923; seamm_mbe.units uses the exact 1e8/101325 = 986.92327. The VASP term, P_atom and P_mol therefore differ by 0.017 atm: P_mol +1795.932 vs +1795.949, P_atom −63009.665 vs −63009.649. The test asserts that the ratio is exactly the two constants’.

  • Speed: enumerating the 64-water frame takes 0.08 s.

Decisions and deviations from the plan#

Agreed with “design”, 2026-10-03:

  1. Validity bound. The cell’s smallest perpendicular width must exceed R + D, where R and D are the radius and diameter a selected fragment can have under its rule:

    • pairs and “compact”: 2c;

    • “hub” and connected triples: 3c;

    • connected 4-body: 5c.

    Contact criteria add 2 r_max per bond. The derivation is in SelectionRules.check(). It refuses rather than warns. For order ≥ 3 the enumerator also refuses two different selected fragments of the same molecules. For pairs it does not need to, since pair names carry non-minimum images. With the pair cutoff beyond L/2 and check() bypassed, both images would be selected. “design” proposed 3c for “compact”; the derivation gives 2c, since all members are within c of each other.

  2. Names are the prototype’s. Selected triples and 4-bodies carry no images, which is unique by the bound. Auxiliary fragments of order ≥ 3 carry their images (_x...) in a periodic cell and nothing in a cluster. Pair images beyond ±1 are spelled out (_x2_0_-1) instead of the prototype’s sign letters, which would have collided.

  3. Molecular level: the closure of the molecular-assigned fragments, with calculations(molecular_everywhere=True) for the consistency check.

  4. Ladders: each increment is built from sub-fragments at its own level, the mixed scheme. The plan said assign_levels would enforce that a periodic fragment’s sub-fragments are periodic. That is not needed and is dropped: each increment is self-consistent within its own ladder, whatever level its sub-fragments are assigned in the sum.

  5. Charges: a structure’s formal charges, when present and not all zero, override the catalog. A disagreement is an error naming the molecule. The parity check follows seamm_bsse’s wording.

  6. Stress sign: the batch contract says only “GPa, as the program gives it”, and seamm_mdi documents the MDI pressure convention. So every external stress needs an explicit convention (“pressure” or “stress”), and to_seamm writes σ = −P. For phase 3: vasp_step.analyze_task must state its stress convention.

  7. No dftd4: the D4 add-on is a CellTerm (or added to the periodic fragment results) supplied by the caller. The route-B twin belongs to the Dispersion step; its pilot numbers are kept in the expected-values file for that step.

Bugs found by the unit tests#

  • A bond given as (i, j) with i > j was dropped from the molecule’s graph, which mistyped it. The molsystem adapter always gives i < j, so the pilot regression could not see it.

  • mbe_correction on fragments without levels reported nothing missing. It now asks for assign_levels first.

Review (2026-10-03)#

“design” and a review subagent read every module and verified each finding with throwaway scripts. The verdict was to approve for the PR after these fixes, all made, with tests:

  1. Collision guard. The same-molecules guard ran on every connected placement before the rule filter, so “hub” and “compact” were held to the connected bound. Four Ar 2.9 Å apart in a 9.1 Å cell, hub 4-body, passed check() but the enumerator refused it. The guard now applies to the placements that pass the rule. The test fails on the old code.

  2. Packaging. The wheel installed a top-level tests: added include = ["seamm_mbe*"]. The sdist lacked tests/data and docs: fixed in MANIFEST.in.

  3. Masses. The table stops at Ba and was looked up eagerly. Masses are now lazy, falling back to molsystem’s table for heavier elements. molsystem’s own midpoint values (O 15.9995) would move P_mol, so the conventional table stays first.

  4. Cells and inputs.

    • A singular cell is refused.

    • An empty periodic system passes check().

    • assign_levels needs real booleans: 1 means 1 Å.

  5. Cutoff tables. Contradictory entries are refused: the same pair twice in __post_init__, and equally specific matches at lookup.

  6. Bond-less salts. A bond-less salt from molsystem is accepted. The halide ions and noble gases count as monatomic, and the catalog gained Na⁺, K⁺, F⁻, Cl⁻, Br⁻ and I⁻.

  7. Increments. increments() on names not closed under sub-fragments gives a clear error.

  8. Images. A molecule bonded to its own periodic image is refused. The minimum_image docstring no longer claims exactness for long vectors in non-reduced cells.

The review also verified, so there is no need to re-audit:

  • triclinic cells against brute force;

  • the R + D bound, including compact = 2c;

  • the ions, charges and parity checks;

  • the order-4 rules and the Möbius inversion, against an independent recursive evaluation;

  • the force scatter, the virial and the ladders;

  • the pressures and the units.

Phase-2 notes:

  • to_seamm returns a (3, 3) stress, but SEAMM’s stress property is Voigt [6], so the step converts it.

  • The design document’s targets now carry the exact-constant values.

Next#

  • The “design” review, then the repo molssi-seamm/seamm_mbe (Paul), the PR via the release-seamm-plugin skill, the release, and GitHub Pages.

  • Phase 2 (mbe_step) with Paul’s OK.