Skip to content

Add proteinbox and protein_mixedbox simulation types - #59

Merged
maxscheurer merged 40 commits into
developfrom
feature/protein-mixedbox
Nov 10, 2026
Merged

maxscheurer merged 40 commits into
developfrom
feature/protein-mixedbox

Conversation

@maxscheurer

Copy link
Copy Markdown
Collaborator

Summary

Adds two protein simulation types. Supersedes #15 (closed): this branch contains all of its commits plus the protein_mixedbox work, and has develop merged in. Implements #13 (proteinbox).

proteinbox — protein in a water box

  • Protein prep: PDBFixer cleanup, gmx pdb2gmx (CHARMM only), disulfides from an explicit list (-ss prompts answered deterministically), residue-name based protonation overrides, multi-chain proteins via an explicit chains declaration.
  • Cubic box with padding, solvation with the existing spc216 tiling, Na⁺/Cl⁻ via the existing ionize_solvated_system, restrained OpenMM relaxation, grompp validation, GROMACS run schedule (em/nvt/npt/md).
  • CHARMM36m force-field registry with pinned SHA-256 download; [gromacs] settings (GMX_PATH, FORCEFIELD_DIR) in the config and mdfactory config init.
  • CSV pipeline support (;-separated chains/disulfides) and YAML-relative pdb_path resolution.

protein_mixedbox — protein + water / ions / SMILES solutes

  • parametrization: charmm (CharmmConfig, currently charmm36m + tip3p only).
  • Solution species (SolutionSpecies: water as smiles: O, ions, arbitrary SMILES) given by count or molar concentration; concentrations are resolved on the box or protein_excluded volume (partial_specific_volume). Sizing by fixed_box or count_density.
  • Solutes parametrized with CGenFF (native CHARMM water/ions short-circuit); packmol packing around the fixed protein; OpenMM compression with a best-fit RMSD restraint; protein made whole and centered; RMSD and clash checks; ionization on the chosen volume basis; relaxation.
  • Composite protein + CGenFF topology (generate_gromacs_topology_with_protein) with atom-type de-duplication against the force field.
  • build_metadata.json (resolved/final counts, achieved concentrations), surfaced in Simulation.metadata and system_chemistry.
  • Analysis/artifact registry entries (system_chemistry, last_frame_pdb); no protein-specific analyses yet.

Notes for review

Changes to shared code (also affect existing builds)

  • relax_with_protein_restraints (used by build_proteinbox) now uses a best-fit RMSD restraint instead of absolute position restraints, and no longer wraps by residue. build_proteinbox has no re-imaging step, so its output is now unwrapped and the protein can move as a rigid body.
  • compress_box (used by the existing mixedbox build) checks the density every 25 steps (was 1000) and raises RuntimeError if the target is not reached or is overshot by more than 2% (previously silent).
  • ionize_solvated_system gained an optional solvent_volume_a3.

Merge of develop: conflicts from #43 were resolved in build.py, models/composition.py and models/input.py. build_proteinbox read the ion counts by position ([0] = Na, [1] = Cl); since #43 returns only the species actually added, a lone Cl⁻ would have been written as NA. Counts are now looked up by residue name (_ion_counts, with tests).

Tests

pytest -k "not build" mdfactory on a compute node: 2998 passed, 25 failed, 53 errors, 1 skipped. The failing/erroring test IDs are identical to the pre-merge head (0de0618), so none are caused by this branch: CGenFF parametrization tests (the local SILCSBIO license is expired) and DB push/pull tests (Column names do not match with existing table in the local environment).

Status

Draft: not yet validated end-to-end. The unit tests mock gmx, pdb2gmx, packmol and OpenMM. The test plan is in the first comment. SMILES-solute builds are currently blocked by the expired CGenFF license (renewal in progress); water/ion-only builds are not affected.

MSiggel and others added 30 commits September 3, 2026 23:00
…Config

New models for the proteinbox simulation type:
- ProteinSpecies: PDB-path-based species with disulfide/protonation annotations
- ProteinBoxComposition: protein + box_padding + ionization config
- Pdb2gmxConfig: force field and water model config for pdb2gmx
- GromacsProteinParameterSet: output paths from pdb2gmx
- BuildInput extended with proteinbox type and pdb2gmx parametrization

Closes: relates to #13

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Functions for the proteinbox build pipeline:
- check_gmx_available: verify gmx binary is on PATH
- clean_pdb: PDBFixer-based PDB standardization
- run_pdb2gmx: subprocess wrapper for gmx pdb2gmx
- extract_charge_from_topology: parse net charge from .top/.itp
- update_topology_molecules: append water/ion entries to [ molecules ]
- validate_with_grompp: dry-run topology validation

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Brief NPT equilibration with protein heavy atoms position-restrained
via CustomExternalForce. Allows water/ions to relax around the fixed
protein structure before GROMACS production runs.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Build pipeline: clean PDB -> pdb2gmx -> center in cubic box -> solvate
(reuses existing solvate()) -> ionize (reuses ionize_solvated_system())
-> update topology -> OpenMM relax with protein restraints -> validate.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
MDP files for protein-in-waterbox equilibration and production:
- em.mdp: steepest descent minimization
- nvt.mdp: NVT with -DPOSRES and Protein/Non-Protein tc-grps
- npt.mdp: NPT with -DPOSRES and Berendsen barostat
- md.mdp: production with Parrinello-Rahman, no position restraints

All use CHARMM36m-specific nonbonded settings (1.2nm cutoffs,
Force-switch VdW modifier, no dispersion correction).

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
- 13 unit tests covering ProteinSpecies, ProteinBoxComposition,
  Pdb2gmxConfig, BuildInput integration, and topology parsing
- Fix ProteinSpecies to default count=1/fraction=1.0 at field level
  (avoids parent Species validator ordering issue)
- Example YAML for lysozyme with CHARMM36m force field

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
Addresses 4 P2 review comments:
- Add missing `from pathlib import Path` to build.py module scope
- Resolve all paths to absolute in run_pdb2gmx before subprocess call
  (avoids cwd confusion with output_dir)
- Add species/total_count/charge properties to ProteinBoxComposition
  and proteinbox case in BuildInput.metadata for analysis compatibility
- Implement _apply_protonation_states: renames residues in PDB before
  pdb2gmx so protonation overrides (HIS->HIE, GLU->GLH) are applied

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
- Add check_forcefield_available() that verifies the force field directory
  exists in GROMACS share, GMXLIB, or cwd before invoking pdb2gmx. Lists
  available force fields on failure.
- Clean up mdout.mdp artifact in validate_with_grompp alongside check.tpr.
- Fix bare open() file handles in relax_with_protein_restraints (use with-blocks).
- Add CHARMM HIS alias translation (HIE→HSE, HID→HSD, HIP→HSP for charmm FFs).
- Use re.fullmatch for protonation state key parsing (reject trailing chars).
- Validate 3-character residue names for PDB column safety.
- Add disulfide bond prompt generation (_build_disulfide_prompt_input) for
  deterministic pdb2gmx -ss interaction.
- Resolve paths before working_directory context switch to avoid cwd breakage.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
- Add FORCEFIELD_REGISTRY mapping friendly names (charmm36m, charmm36m-ljpme)
  to MacKerell lab download URLs and extracted directory names.
- Auto-download missing force fields on first use from the registry.
  Stored in platformdirs user_data_dir (~/Library/Application Support/mdfactory/forcefields/).
- resolve_forcefield() translates friendly names to actual directory stems
  so users write "charmm36m" in YAML and pdb2gmx receives the correct
  directory name.
- Inject GMXLIB in subprocess env so gmx finds downloaded force fields.
- Fix extract_charge_from_topology: skip .ff/ library includes (ions.itp,
  tip3p.itp) that contain atom type templates, not system charges. Was
  causing +159 charge for lysozyme with charmm36m instead of +8.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
- Add [gromacs] section to settings.py with GMX_PATH and FORCEFIELD_DIR,
  following the same pattern as [cgenff] SILCSBIODIR.
- Settings.__init__ auto-prepends FORCEFIELD_DIR to GMXLIB env var on
  startup, matching how SILCSBIODIR is auto-set for CGenFF.
- check_gmx_available() checks configured GMX_PATH first, falls back
  to PATH lookup.
- All gmx subprocess calls (pdb2gmx, grompp) use the configured binary
  and inject GMXLIB for force field resolution.
- Add GROMACS setup to config wizard (sync_config.py): prompts for gmx
  path, forcefield dir, and offers to download CHARMM36m on setup.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
- rename Pdb2gmxConfig.ignh to ignore_hydrogens
- rename ProteinBoxComposition.box_padding to padding (matches LNP)
- accept disulfide bonds as CYS-prefixed residue references (e.g. CYS6)
  to match the protonation_states residue nomenclature

Addresses review comments on examples/proteinbox/lysozyme_charmm36m.yaml.
Address 7 issues flagged in review of the protein-box feature:

- [P1] Resolve a relative protein pdb_path against the YAML file's
  directory in run_build_from_file, not the process cwd (the CLI cd's
  into the output dir before building), so bundled paths like ./1aki.pdb
  are found.
- [P1] Declare pdbfixer in [tool.pixi.dependencies]; clean_pdb imports
  it at runtime but it was undeclared.
- [P2] Run pdb2gmx with relative output filenames so the generated
  topology has a portable relative #include "posre.itp" instead of an
  absolute path.
- [P2] Enforce simulation_type <-> parametrization <-> config-type
  consistency in BuildInput; proteinbox+cgenff and mismatched configs
  now fail validation with a clear message instead of late in the build.
- [P2] Restrict Pdb2gmxConfig to CHARMM force fields and 3-site water
  (tip3p/spc/spce), matching what solvation and the mdp files support.
- [P2] Let grompp validation failure fail the build instead of being
  swallowed as a warning.
- [P2] Center the protein on its bounding-box center (not center of
  geometry) so every box face keeps the requested padding.
Rewrite the absolute force-field #includes that pdb2gmx writes into a
relative path and copy the .ff directory next to the topology, mirroring
the CGenFF strategy in generate_gromacs_topology. Both grompp (cwd) and
OpenMM's GromacsTopFile (topology dir) then resolve the force field from
the build directory, so it is portable to another machine.

Reject CHARMM LJ-PME force-field variants (e.g. charmm36m-ljpme): the
proteinbox run schedule uses cutoff/Force-switch LJ, which is
incompatible with LJ-PME.

Move the newly added test imports to module top level.
Multi-chain PDBs now build each subunit as its own Protein_chain_<ID>
moleculetype, referenced from topol.top, rather than requiring a merge.
The chains list must be declared in the YAML and is validated against the
chains pdb2gmx produces; an undeclared multi-chain PDB is an error so
subunits are never silently dropped. GromacsProteinParameterSet carries
the per-chain include files (topol_Protein_chain_*.itp, posre*.itp), which
build_proteinbox copies alongside topology.top in both the build and
relaxation directories.

Also fixes two review findings:
- Protonation overrides match histidine tautomers as a family and raise
  when an override matches no residue; 4-character CHARMM acid states
  (GLUP/ASPP) are rejected with a clear message.
- bundle_forcefield_into_topology copies each .ff directory only once
  instead of once per include line.
Protein charge is set from the pdb2gmx topology, and every caller treats
species.charge as Optional[int] (guarding with 'is not None'). Raising
NotImplementedError broke that contract and would propagate through the
hasattr guard in _compute_charge_from_universe if a protein ever reached
it. Returning None aligns with the contract and is skipped by the charge
sums, leaving proteinbox charge results (from the topology) unchanged.
- Parse delimited protein list cells (chains 'A;B', disulfide pairs
  'CYS6-CYS127' joined by ';') and dotted protonation_states columns
  from CSV rows; drop blank optional protein cells per row and exclude
  them from the strict NaN guard so heterogeneous proteins share one CSV.
- Resolve relative protein pdb_path against the CSV directory at read time.
- Reject merge_all combined with declared protein.chains at config
  validation instead of failing deep inside pdb2gmx.
- Make relax_steps config-driven on ProteinBoxComposition.
- Remove proteinbox build intermediates on success so the output dir
  holds only final artifacts, matching small-molecule builds.
Serialize prepared per-system YAML with model_dump(mode="json") so
proteinbox Path fields (pdb_path) and tuple fields (disulfide_bonds)
become plain str/list that yaml.safe_dump can write. Small-molecule
inputs have no Path/tuple fields, so this path never surfaced before.

Correct the build_proteinbox cleanup comment: only the three
pdb2gmx-specific loose intermediates are removed; solvated.pdb and
relaxation/ are left in place, matching how the bilayer build leaves
solvated.pdb and bilayer_squeeze/.

Add examples/proteinbox/proteins.csv showing the flat dotted-key /
delimited-cell CSV encoding for single-chain (lysozyme) and
multi-chain (insulin) systems. PDBs are referenced by RCSB id
filename and fetched separately, matching the existing YAML example.
Extend compress_box with optional protein position restraints and a
target-volume stop for the protein_mixedbox pipeline. When protein_indices
is given, every listed atom is held with a harmonic CustomExternalForce
through minimization and high-pressure compression so a pdb2gmx-prepared
protein stays fixed while solvent and solutes pack around it. When
target_volume is given, compression stops on box volume instead of density
(the fixed_box sizing mode). Defaults preserve the existing mixedbox
behavior exactly.
generate_gromacs_topology_with_protein assembles one topology for a
pdb2gmx-prepared protein packed with CGenFF small molecules, water, and
ions, in strict include order: bundled CHARMM force field, merged CGenFF
parameters (atomtypes de-duplicated against the force field, [defaults]
dropped since the FF already bundles CGenFF), protein moleculetypes and
position restraints, CGenFF small-molecule moleculetypes, then the water
and ion topologies pdb2gmx already includes. The [ molecules ] section
keeps the protein chain entries and appends the solution molecules in
coordinate order so topology and coordinates agree. Native water and ions
reuse the pdb2gmx force-field moleculetypes and contribute no extra
include or parameters.
…dbox

extract_mass_from_topology sums the topology mass column (following protein
chain includes) so the build can derive the protein displacement volume.
ionize_solvated_system gains an optional solvent_volume_a3 to resolve the salt
pair count on an explicit volume basis instead of the water-count proxy.
Pack water, ions, and SMILES solutes around a pdb2gmx-prepared protein under
one CHARMM/CGenFF force-field family, compress with every protein atom
restrained, re-image the protein whole and centered across PBC, ionize on the
chosen concentration-volume basis, and relax with the protein restrained.

compress_box skips its residue wrap when protein_indices is given, since a
restrained protein ends up torn across the boundary; build_protein_mixedbox
makes it whole from the restraint reference shape and re-centers instead.
@CLAassistant

CLAassistant commented Oct 20, 2026 •

Copy link
Copy Markdown

CLA assistant check
All committers have signed the CLA.

@maxscheurer

Copy link
Copy Markdown
Collaborator Author

End-to-end run of the test plan (2026-11-02)

I ran every item of the test plan above for real, on the gpu partition (L40S, gromacs/2024.3-gpu) and the cpu partition. Passing items are ticked in the test-plan comment. Four blocking bugs turned up; they are listed below and fixes follow on this branch.

Environment notes

  • The cluster modules ship only gmx_mpi. gromacs/2024.3-cpu and -gpu have no gmx, and check_gmx_available looks only for gmx, so every protein build fails until [gromacs] gmx_path points at gmx_mpi. With that set, pdb2gmx, grompp and simulate all work.
  • mackerell.umaryland.edu is unreachable from the cluster network (TLS handshake timeout), so the mdfactory config init download fails. A manually downloaded charmm36-feb2026_cgenff-5.0.ff.tgz matches the pinned SHA-256 and works.
  • The CGenFF license is renewed: SILCSBIO 2026.1.0, CGenFF 5.0.

Blocking bugs

  1. The stricter compress_box randomly overshoots its target.
    • 13 of 27 protein_mixedbox compressions overshot the target volume or density by 3.4–7.8%, above the 2% limit. One more crashed with NaN encountered during RMSD force calculation.
    • Identical inputs give different results: the same ethanol YAML overshot by 3.8% on one build and 6.3% on the next.
    • Cause: the 1000 bar barostat. Water at 1000 bar settles about 4.5% denser than at 1 bar, so the box sweeps through the liquid-density target between two 25-step checks. In one run the box went from 385.5 to 356.5 nm³ within 150 steps, against a target of 373.2 nm³.
    • When fixed_box does pass, it still ends 0.2–0.8% under the requested volume.
  2. Both LNP examples no longer build. The core compression reaches only 0.935 and 0.922 g/mL against the 0.95 target in 200k steps. On develop, compress_box ignored whether the target was reached and carried on; it now raises RuntimeError.
  3. CGenFF solutes get the wrong atom names in system.pdb.
    • Packing keeps the OpenFF element-indexed names (C1 C2 O1 H1 … H6), while the CGenFF .itp numbers atoms by position (C1 C2 O3 H4 … H9).
    • Atom order is identical, so coordinates are correct. But grompp -maxwarn 0 reports "1498 non-matching atom names" (214 ethanols × 7) and the build validation fails on every ethanol protein_mixedbox build: 6 of 6 that got past compression.
    • The plain mixedbox CGenFF path has the same mismatch, but no build check catches it there.
  4. neutralize: false with concentration: 0 always fails. Lysozyme keeps its +8 charge, and grompp's "Ewald electrostatics in a system with net charge" warning fails the -maxwarn 0 validation. The test plan expects this build to complete.

Passed

  • proteinbox: YAML and CSV, including 4-chain insulin.
    • [ molecules ] matches system.pdb, the net charge is 0, and grompp -maxwarn 0 passes.
    • Ions are neutralization plus about 0.15 M: lysozyme +8 gives NA 23 / CL 31.
    • HIS15 → HSE is applied. The protein is whole (CA–CA gaps ≤ 4.0 Å within chains) and centered.
  • Error cases: a wrong override (HIS16), a chain mismatch ([A, B]) and a bad disulfide (CYS999) each fail with a clear message.
  • protein_mixedbox, water only:
    • fixed_box gives 55.03 M achieved against 55 M requested; count_density also builds.
    • Protein RMSD is 0.26–0.32 Å after compression and relaxation, and there are no clash errors.
  • protein_mixedbox, other cases:
    • The one-sided case adds only 8 CL.
    • Merged multi-chain insulin builds, with RMSD 0.27 Å.
  • CGenFF solutes:
    • The CGenFF atom types merge cleanly with the Feb-2026 bundle; validate_cgenff_atomtype_compatibility passes.
    • Achieved ethanol is 1.006–1.01 M against 1.0 M requested.
    • The cached ethanol .itp is reused by later builds.
  • Names in the bundle: ions.itp defines NA/CL (alongside SOD/CLA), and tip3p.itp is SOL with OW/HW1/HW2.
  • Simulation: mdfactory simulate --stages EM NVT ran for lysozyme, insulin (with per-chain posre_*.itp), and the water-only, one-sided and merged-insulin protein_mixedbox builds.
    • Every EM converged to Fmax < 1000 in 44–255 steps.
    • The 100 ps NVT runs were clean at 740–1010 ns/day.
    • tc-grps = Protein Non-Protein and -DPOSRES work for multi-chain systems.
  • Shared code: mixedbox (CGenFF and SMIRNOFF) and bilayer (CGenFF and SMIRNOFF) still build.

Minor, not blocking

  • system.pdb loses chain IDs (every atom is chain X), and auto-assigned histidines stay HIS in the coordinates while the topology has HSD/HSE.
  • Insulin built without listed disulfides has no S–S bonds and gives no warning (already listed under known limitations).
  • proteinbox and protein_mixedbox insulin differ by 8 atoms (1572 vs 1564) when merge_all is used. Not yet explained.

Amber14SB check for issue 60

  • With GROMACS 2026.3 from conda-forge, pdb2gmx -ff amber14sb -water tip3p on lysozyme loads in OpenMM 8.4 GromacsTopFile without errors.
  • Single-point energies match GROMACS 2026.3:
    • bonds and angles: < 0.02 kJ/mol
    • total torsions: < 0.001 kJ/mol
    • nonbonded, protein only: 0.2 of −20,130 kJ/mol
    • nonbonded, solvated: 33 of −519,900 kJ/mol
  • Turning LEaP reordering off changes nothing for lysozyme.

…charge grompp validation

Findings from the end-to-end run of the protein_mixedbox test plan:

- compress_box: the 1000 bar barostat sweeps past the target between checks
  (3.4-7.8% in about half of the builds). Overshoots up to 15% are now undone by
  expanding the box and molecule centers to the exact target volume
  (finish_compression / rescale_to_volume); a restrained protein moves as one
  rigid unit.
- LNP core compression warns instead of failing when the nominal density is
  not reached (require_convergence=False), restoring the develop behavior.
- protein_mixedbox writes CGenFF solutes with their ITP atom names
  (topology_atom_names), so grompp -maxwarn 0 no longer sees non-matching names.
- validate_with_grompp(allow_net_charge=...) accepts only grompp's net-charge
  warning when neutralization is disabled, in a single grompp call (a failing
  gmx_mpi call aborts through MPI and can kill an enclosing srun step).
- Reformat test_topology_utilities.py with ruff 0.9.2 (the CI version).
@maxscheurer

Copy link
Copy Markdown
Collaborator Author

Progress update: fixes for the four blocking e2e bugs

This addresses the four blocking bugs in the e2e findings comment above.

  1. compress_box overshoot
    • Fix: after compression, finish_compression expands the box and molecule centers to the exact target volume (rescale_to_volume), so the result no longer depends on how far the barostat overshot. A restrained protein is moved as one rigid unit. Overshoots above 15% still raise.
    • Rerun: 13 of 13 compression-based builds passed: 12 protein_mixedbox and mixedbox builds, plus the CGenFF LNP core. Overshoots corrected were 0.1–3.3%, and fixed_box now ends at exactly 373.248 nm³ (72³ ų).
  2. LNP regression
    • Fix: create_mixedbox_by_compression(..., require_target_density=False) for the LNP core warns instead of raising when 0.95 g/mL is not reached, which was the behavior on develop. Other system types still require the target.
    • Rerun: both LNP examples now build end to end (CGenFF: core at 0.933 g/mL, 24 min; SMIRNOFF: 15 min).
  3. CGenFF atom names
    • Fix: protein_mixedbox now writes CGenFF solutes with their .itp atom names (topology_atom_names). The ITP atoms are checked against the packed molecule by count and mass, and a mismatch raises.
    • Rerun: 5 of 5 ethanol builds (fixed box ×3, count density, CSV) pass grompp -maxwarn 0. Achieved ethanol is 0.999 M, and EM plus 100 ps NVT run cleanly.
  4. neutralize: false with concentration: 0
    • Fix: validate_with_grompp(allow_net_charge=True) accepts grompp's net-charge warning and nothing else. It makes a single grompp call, because a failing gmx_mpi aborts through MPI and can kill an enclosing srun step.
    • Rerun: the build completes, the warning is logged, and EM plus NVT run cleanly.

Also:

  • Reformatted test_topology_utilities.py with ruff 0.9.2, the version CI uses. This was the existing Formatting CI failure.
  • Added 18 regression tests in mdfactory/tests/test_e2e_regressions.py.
  • pytest -k "not build" mdfactory: 3094 passed, none failing; the unchanged branch had 3076 passing, none failing.

Still open, not addressed here:

  • Random NaN blow-up during the 1000 bar compression. Seen in 2 of about 37 builds; OpenMM reports it as NaN encountered during RMSD force calculation because the coordinates were already NaN. It looks like the same failure as issue 58. Merged insulin then built 3 of 3 times.
  • Plain mixedbox CGenFF builds still have mismatched water and solute atom names in system.pdb (pre-existing; there is no grompp check on that path).
  • gmx_mpi detection: check_gmx_available still needs [gromacs] gmx_path on this cluster, whose modules ship only gmx_mpi.
  • PR description: the line "RuntimeError if overshot by more than 2%" is now outdated.

Commit: d0c0f19


Progress tracked by mach6

@maxscheurer
maxscheurer marked this pull request as ready for review November 2, 2026 12:11
@maxscheurer

Copy link
Copy Markdown
Collaborator Author

Unverified Review Candidates — Pending Assessment

Review round: 1
Reviewed commit: d0c0f19

These findings have not been verified yet. Severity is the reviewer's confidence, not a shipping decision. None of them is a merge blocker until the assessment comment is posted.

Critical

None.

Important

Finding 1: CSV push and discovery hash the relative pdb_path, but prepare and build hash the resolved one (code-reviewer, 85)

  • Where: mdfactory/cli.py (around lines 91–106) and mdfactory/utils/push.py (load_models_from_csv, around lines 111–128 and 247–264).
  • The mismatch: df_models_from_input_csv resolves system.protein.pdb_path against the CSV's directory before it builds the models. push.load_models_from_csv reads the same CSV and calls df_to_build_input_models(df) without that step. BuildInput.hash covers pdb_path, so the two code paths get different hashes for the same row.
  • Trigger: a proteinbox or protein_mixedbox CSV with a relative pdb_path, which this PR adds support for, followed by mdfactory push --csv … or analysis/artifact discovery through discover_simulation_folders(csv=...).
  • Impact: every protein row logs "No folder found for hash", and none of those simulations is pushed.

Finding 2: single-YAML builds through --config (Parsl) skip the YAML-relative pdb_path resolution (code-reviewer, 80)

  • Where: mdfactory/cli.py _build_from_yaml (around lines 334–365), mdfactory/workflows.py (lines 54–85) and mdfactory/orchestration/apps.py (around lines 444–452).
  • The gap: _resolve_proteinbox_pdb_path runs only inside run_build_from_file. _build_from_yaml builds BuildInput(**data) from the raw YAML and hands it to build_systems. The Parsl app then changes into output/<hash>/ and calls run_build_from_dict, so a relative pdb_path is resolved against the build directory instead of the YAML's directory.
  • Impact: mdfactory build sys.yaml --config exec.yaml with pdb_path: protein.pdb fails with file-not-found. On the sequential path, build_dir is named from the hash of the unresolved path, but the model that actually gets built has a different hash.

Finding 3: neutral-lysine (LSN) residues drop out of the proteinbox relaxation restraint (error-auditor, 80)

  • Where: mdfactory/build.py, around line 856: u_ionized.select_atoms("protein and not name H*").
  • Trigger: a LYS<n>: LYN/LSN override, which the alias table accepts. pdb2gmx then writes the residue as LSN, and the MDAnalysis protein selection does not include LSN.
  • Impact: that residue's heavy atoms are not restrained during the proteinbox relaxation, and nothing warns about it. build_proteinbox has no RMSD or clash check after relaxation to catch it. protein_mixedbox selects atoms with range(n_protein_atoms) and is not affected.

Finding 4: issue 13 asks for a selectable force field including amber99sb-ildn, but only CHARMM is accepted (completeness-checker, 90)

  • What the code does: Pdb2gmxConfig._validate_charmm_forcefield_and_water rejects any force field that does not start with charmm.
  • What the PR says: the body says "CHARMM only", and the known-limitations list says "No … Amber protein path". The Amber work is tracked separately in issues 60, 61 and 62.
  • Open point: I found no maintainer comment approving the narrower scope against issue 13.

Finding 5: the LNP leniency and the stricter mixedbox default have no tests, although e2e bug 2 was exactly this regression (test-reviewer, 88)

  • Where: create_lnp_core and create_mixedbox_by_compression in mdfactory/build.py. No test calls either function.
  • What is unguarded: nothing stops a refactor from dropping require_target_density=False, which would break LNP builds again. The finish_compression interaction is also untested: require_convergence=False combined with an overshoot beyond tolerance should still raise.

Finding 6: compress_box orchestration is not tested end-to-end with mocks (test-reviewer, 85)

  • Where: mdfactory/simulation/openmm_utils.py, compress_box, which is shared with mixedbox and LNP.
  • What is unguarded:
    • the ų→nm³ conversion of target_volume;
    • the target volume derived from density;
    • forwarding of density_target, rigid_indices and require_convergence;
    • skipping the residue wrap when protein_indices is set.
  • Impact: a swapped or mis-scaled argument passes the unit suite.

Finding 7: the count and report bookkeeping in build_protein_mixedbox and build_proteinbox is untested (test-reviewer, 82)

  • What is unguarded:
    • the water count reduced by n_replaced;
    • alignment of final_input_counts when a species resolves to 0;
    • allow_net_charge=not neutralize;
    • the build_metadata.json writer together with its system_chemistry reader (the reader is only tested against hand-written JSON).
  • History: all four e2e bugs were in orchestration code like this.

Suggestions

Finding 8: out-of-date docstrings (code-reviewer, 90)

  • The protein_indices docstring of compress_box still says "harmonic position restraints".
  • The Returns section of relax_with_protein_restraints says the protein stays "at their restrained coordinates". The restraint is a best-fit RMSD restraint, so the protein can still move as a rigid body.

Finding 9: the PR body is out of date on compress_box (completeness-checker, 92). It still says overshoot beyond 2% raises. The code now rescales to the exact target and raises only above 15%, and require_convergence=False (LNP) warns instead of raising. The test plan has the same stale line.

Finding 10: test-plan checkboxes not updated (completeness-checker, 80). These are still unticked although the progress comment reports them as passing:

  • neutralize: false with concentration: 0;
  • bilayer/lnp unaffected;
  • the ethanol examples;
  • the CGenFF/Feb-2026 atom-type merge.

Finding 11: no test of the density-compression loop (test-reviewer, 80). simulate_compression_until_density_reached now checks every 25 steps and exits early. The volume variant has a FakeSimulation test; the density variant has none.

Finding 12: weak RMSD-restraint test (test-reviewer, 85)

  • test_rmsd_restraint_is_translation_invariant only shows that translation leaves the energy unchanged. A force that is identically zero would also pass.
  • No check that deforming the protein raises the energy.
  • The test uses the default platform, which can be single precision, while asserting abs=1e-8.

Finding 13: relax_with_protein_restraints has no test (test-reviewer, 80). Nothing checks that pressure=None adds no barostat, that a single RMSD restraint is added, or that the output is not wrapped.

Finding 14: merge_extra_parameter_itps(exclude_atomtypes, drop_defaults) has no direct test (test-reviewer, 80). The composite-topology test never asserts that [ defaults ] was removed.

Finding 15: four almost identical [ atoms ] column parsers (simplifier, 92). The charge and mass extractors in mdfactory/setup/protein.py (around lines 786–925) differ only in the column index. One _sum_atoms_column(path, column, follow_includes) helper could replace them.

Finding 16: the concentration-volume-basis branch appears three times (simplifier, 88). It is in build.py (around lines 1258 and 1368) and in setup_utilities.resolve_solution_counts. A shared concentration_volume_a3(...) helper would cover all three.

Finding 17: build_proteinbox and build_protein_mixedbox repeat the same prelude and epilogue (simplifier, 85). Both have the same PDB clean, run-schedule copy plus grompp validation, and intermediate cleanup.

Finding 18: topology and force-field staging copied three times, relaxation-directory staging twice, in build.py (simplifier, 82).

Finding 19: two compress_box calls that differ in one keyword argument (simplifier, 90). In build.py (around lines 1216–1235), only target_volume= vs target_density= differs.

Finding 20: the RMSD and clash verification block after compression and after relaxation is duplicated (simplifier, 80). It is in build.py, around lines 1240 and 1354.

Finding 21: one-use wrapper collect_forcefield_atomtypes (simplifier, 80). The force-field directory is parsed twice per topology generation, in topology_utilities.py and parametrize.py.

Finding 22: verbose index bookkeeping for final_input_counts (simplifier, 82). The loop in build.py (around lines 1293–1302) could be a short comprehension over an iterator.

Strengths

  • All four e2e blocking fixes are present and wired in, with regression tests in test_e2e_regressions.py:
    • finish_compression/rescale_to_volume;
    • LNP require_target_density=False;
    • CGenFF topology_atom_names, checked by count and mass;
    • validate_with_grompp(allow_net_charge): a single call that accepts only the net-charge warning.
  • Order of [ molecules ] against the coordinates is correct:
    • proteinbox: protein → SOL → NA → CL;
    • protein_mixedbox: the packed species, then the ions.
  • _ion_counts looks ions up by residue name, so one-sided neutralization comes out right.
  • Concentration resolution, conversion of protein displaced volume, and periodic-image-consistent rigid rescaling all hold up under static review.
  • Existing mixedbox, bilayer and LNP hashes and topology defaults are unchanged. ionize_solvated_system behaves exactly as before when solvent_volume_a3=None.
  • There is a thorough real end-to-end validation record on GPU and CPU nodes.

Agents run: code-reviewer (split into 4 parts after 2 stalled runs), error-auditor, test-reviewer, completeness-checker, simplifier


Reviewed by mach6

@maxscheurer

Copy link
Copy Markdown
Collaborator Author

Review Assessment

Findings comment: #59 (comment)

Scope basis: The maintainer-authored test plan, under "Known limitations / follow-ups (not in this PR)", explicitly narrows issue 13's "also supports amber99sb-ildn" bullet ("No OpenFF/SMIRNOFF or Amber protein path"). Amber is tracked separately in issues 60, 61 and 62.

Classifications

Finding Classification Reasoning
1: CSV push and discovery hash the relative pdb_path merge blocker Factual: cli.df_models_from_input_csv (cli.py:91-106) resolves system.protein.pdb_path against the CSV's directory. push.load_models_from_csv (push.py:111-128) does not. BuildInput.hash covers pdb_path, so discover_simulation_folders searches with hashes that don't match the build folders. Scope: The PR claims "CSV pipeline support", and the mismatch comes from adding the resolution to only one of the two CSV loaders. Practical: Both shipped example CSVs use relative paths (1aki.pdb, ../proteinbox/1aki.pdb). With them, mdfactory push/sync/analysis/artifact push --csv logs "No folder found for hash" for every protein row and pushes nothing.
2: single-YAML --config/--dry-run skips YAML-relative pdb_path resolution merge blocker Factual: cli._build_from_yaml (cli.py:353-361) builds BuildInput(**data) from the raw YAML. _resolve_proteinbox_pdb_path runs only in run_build_from_file (workflows.py:67). The Parsl app changes into output/<hash>/ before the build runs. Scope: The PR claims "YAML-relative pdb_path resolution", but it holds only on the local sequential path. Practical: mdfactory build examples/proteinbox/lysozyme_charmm36m.yaml --config exec.yaml, which is the HPC route, fails with file-not-found for the shipped examples. Dry-run and sequential builds also report a different hash than the model that actually gets built.
3: LSN lysines drop out of the proteinbox relaxation restraint deferred Factual: LYN→LSN is accepted (protein.py:61). MDAnalysis prot_res lacks LSN, and build.py:856 selects atoms with protein. Scope: Edge case of an opt-in override; not in the issue or test plan. Practical: Only that one bonded side chain goes unrestrained during a brief (around 100-step) relaxation, and GROMACS -DPOSRES applies during equilibration anyway. A useful follow-up is to select atoms by pdb2gmx atom index.
4: issue 13 asks for a selectable FF including Amber deferred Factual: parametrization.py:41 rejects non-CHARMM force fields. Scope: The maintainer test plan explicitly excludes Amber from this PR, and issues 60/61/62 track it. Practical: Non-CHARMM input fails with a clear error.
5: LNP leniency and mixedbox strict default untested deferred Factual: No test calls create_lnp_core or create_mixedbox_by_compression. Scope: Missing test; no acceptance criterion requires it. Practical: If the flag regressed, LNP builds would fail loudly and immediately. The finish_compression warn/raise behavior is tested. Useful follow-up.
6: compress_box orchestration not mock-tested deferred Factual: No such test exists. Scope: Missing test. Practical: The e2e run validated 13/13 compression builds with the exact target volume, so a mis-scaling would be visible immediately.
7: build count/report bookkeeping untested deferred Factual: No such test exists. Scope: Missing test. Practical: The in-build grompp -maxwarn 0 check catches [ molecules ] mismatches, and the e2e count checks passed.
8: stale docstrings (compress_box protein_indices, relax_with_protein_restraints Returns) nitpick Factual: Confirmed at openmm_utils.py:565 and :993. Scope/Practical: Documentation only.
9: PR body out of date on compress_box nitpick Factual: Confirmed; the progress comment already flags it. Scope/Practical: PR text only; update it before merge as housekeeping.
10: test-plan checkboxes not ticked nitpick Factual: The progress and e2e comments report those items as passing. Scope/Practical: Housekeeping only.
11: density-compression loop untested deferred Factual: No such test exists. Scope: Missing test. Practical: The e2e runs exercised the loop, and no undetected regression was named.
12: weak RMSD-restraint test deferred Factual: The test only checks translation invariance. Scope: Test quality. Practical: The e2e run measured an RMSD of 0.26–0.32 Å, which shows the restraint holds.
13: relax_with_protein_restraints untested deferred Factual: No such test exists. Scope: Missing test. Practical: Validated by the e2e runs and grompp.
14: merge_extra_parameter_itps drop options untested deferred Factual: No direct test exists. Scope: Missing test. Practical: A duplicated [ defaults ] would fail grompp validation inside the build, so it can't ship silently.
15: four near-identical [ atoms ] parsers nitpick Factual: The duplication exists. Scope/Practical: Refactor; no behavioral effect.
16: concentration-volume-basis branch appears 3× nitpick Factual: The duplication exists. Scope/Practical: Refactor only.
17: shared prelude/epilogue in the two builders nitpick Factual: The duplication exists. Scope/Practical: Refactor only.
18: topology/FF staging repeated nitpick Factual: The duplication exists. Scope/Practical: Refactor only.
19: two compress_box calls differing by one keyword nitpick Factual: Valid. Scope/Practical: Style.
20: RMSD/clash check block duplicated nitpick Factual: The duplication exists. Scope/Practical: Refactor only.
21: one-use collect_forcefield_atomtypes wrapper nitpick Factual: The wrapper exists, and the FF directory is parsed twice. Scope/Practical: Negligible cost.
22: verbose final_input_counts loop nitpick Factual: Valid. Scope/Practical: Style.

Action Plan

  1. Finding 1: Move the CSV-relative pdb_path resolution into one shared helper (for example in prepare.py). Call it from both cli.df_models_from_input_csv and push.load_models_from_csv, so push, sync and analysis discovery with --csv compute the same hashes as prepare and build. Add a small test that both loaders give the same BuildInput.hash for a CSV with a relative pdb_path.
  2. Finding 2: In cli._build_from_yaml, call _resolve_proteinbox_pdb_path(data, input.parent) before BuildInput(**data). The Parsl, dry-run and sequential paths will then all use the resolved path and the same hash. Add a test that the model passed to build_systems carries an absolute pdb_path.

Assessment by mach6

…try points

- Add prepare.load_build_input_csv as the single CSV loader for prepare-build,
  build and push/discovery, so relative system.protein.pdb_path cells resolve
  against the CSV's directory everywhere and hashes match the built folders.
- Move the YAML resolver to prepare.resolve_protein_pdb_path and call it in
  cli._build_from_yaml before hashing, so --config (Parsl) and --dry-run find
  the PDB and the build directory is named after the model that is built.
- Add regression tests for the shared CSV hash, the dispatched model path and
  the sequential build-directory hash.
@maxscheurer

Copy link
Copy Markdown
Collaborator Author

Progress Update: the two round 1 merge blockers are fixed

This fixes findings 1 and 2 from the round 1 assessment. Both were cases where one entry point did not resolve a relative pdb_path.

  1. Finding 1: push --csv and analysis discovery did not find protein builds.
    • Problem: push.load_models_from_csv hashed the unresolved relative pdb_path, while prepare-build hashed the resolved one. The rows never matched a build folder.
    • Fix: mdfactory/prepare.py now has one CSV loader, load_build_input_csv. It drops Unnamed columns and resolves relative system.protein.pdb_path against the CSV's own directory. cli.df_models_from_input_csv and push.load_models_from_csv both use it, so they now compute the same hashes.
  2. Finding 2: mdfactory build <yaml> --config / --dry-run skipped the YAML-relative pdb_path resolution.
    • Fix: the resolver moved from workflows.py to prepare.resolve_protein_pdb_path, and workflows.py still imports it under the old name. cli._build_from_yaml now calls it before BuildInput(**data).
    • Result: the Parsl worker finds the PDB from inside output/<hash>/. Dry-run and sequential builds now use the hash of the model that is actually built.

Tests: three new tests in TestRelativePdbPathAcrossEntryPoints (test_proteinbox.py). Each runs from an unrelated working directory and checks one of these:

  • the prepare and push CSV loaders give the same absolute path and the same hash;
  • the model passed to build_systems has an absolute pdb_path;
  • the sequential build directory is named after the hash of the built model.

All three fail without the fix and pass with it. pytest -k "not build" mdfactory on a compute node: 3093 passed, 0 failed. The pre-commit hooks (ruff 0.9.2) pass.

Not addressed here: the deferred follow-ups and nitpicks from the assessment. That includes findings 8–10, the stale docstrings, the outdated PR-body line about compress_box, and the unticked test-plan boxes.

Commit: 1248bb4


Progress tracked by mach6

…e mdp templates, consolidate tests

- protein_mixedbox reuses the identical proteinbox mdp files via a run-schedule
  template_dir key
- build.py: shared helpers for protein prep, relaxation, topology staging,
  reimage/verify and finalization; drop the redundant second force-field bundle
- setup/protein.py: one [ atoms ] column parser, shared residue resolution and
  gmx subprocess helper, drop thin settings wrappers
- openmm_utils: fixed-box compression uses the density loop; shared
  universe/reporter helpers
- models: shared CHARMM config base, Literal single-copy ProteinSpecies, drop
  unused charge properties and single-use wrappers
- tests: shared fixtures and parametrization across proteinbox/protein_mixedbox;
  rename tests that -k "not build" was silently deselecting
@maxscheurer

Copy link
Copy Markdown
Collaborator Author

Progress Update

Simplification/duplication pass over the PR (no feature changes intended). PR additions vs develop drop from ~7255 to ~5920 lines (source ~4212 → ~3540, tests 2780 → 2119).

Source

  • protein_mixedbox mdp files were byte-identical to proteinbox; removed and reused via a new template_dir key in run_schedules.yaml.
  • build.py: build_proteinbox / build_protein_mixedbox share _prepare_protein, _relax_protein_system, _finalize_protein_build, _stage_protein_topology, _reimage_and_verify; single compress_box call; concentration-volume basis computed by one helper in setup_utilities.
  • The redundant second bundle_forcefield_into_topology call in protein_mixedbox is gone (it re-resolved and re-copied the already bundled .ff via gmx search paths).
  • setup/protein.py: four copy-pasted [ atoms ] parsers → one; shared residue-reference resolution and gmx subprocess helper; settings wrappers removed.
  • openmm_utils.py: volume compression loop removed; fixed-box targets convert to an equivalent density. Shared _universe_from_simulation / reporter helpers.
  • Models: shared CHARMM config base for Pdb2gmxConfig/CharmmConfig; ProteinSpecies single-copy via Literal; unused .charge properties and single-use wrappers removed.

Small behavior changes

  • Fixed-box overshoot is now measured as density excess (target/V − 1, slightly stricter at the same 15%); one overshoot_tolerance parameter.
  • ProteinSpecies rejects the string "1" for count.
  • RMSD and concentration-volume errors use one message each; grompp scratch files are also removed on timeout.

Tests

  • Shared conftest fixtures, parametrized proteinbox/protein_mixedbox suites, literal-default tests dropped.
  • Renamed tests whose names contained "build" and were silently deselected by -k "not build" (pass count 3093 → 3099).
  • pytest -k "not build" mdfactory: 3099 passed, 1 skipped; ruff clean.

Still to do: the refactored build functions have no direct unit tests. Re-run one real proteinbox and one protein_mixedbox e2e build before merge.

Commit: 47c2298


Progress tracked by mach6

@maxscheurer

Copy link
Copy Markdown
Collaborator Author

Unverified Review Candidates — Pending Assessment

Review round: 2
Reviewed commit: 47c2298

These are unverified candidates. Severity reflects reviewer confidence; do not treat any item as a merge blocker until the assessment comment is posted.

Round 1 merge blockers, re-checked at this commit:

  • Finding 1 (CSV push/discovery hash mismatch): fixed. prepare.load_build_input_csv is the only CSV loader, used by cli.prepare_build, _build_from_csv, check_csv and push.load_models_from_csv. Both loaders give identical hashes on the example CSVs.
  • Finding 2 (--config/--dry-run YAML-relative pdb_path): fixed. cli._build_from_yaml calls resolve_protein_pdb_path before BuildInput(**data). TestRelativePdbPathAcrossEntryPoints covers both fixes.

Critical

None.

Important

Finding 1: clean_pdb quietly adds missing residues from SEQRES (error-auditor, 85)

  • Where: mdfactory/setup/protein.py:519-523, which calls findMissingResidues() and then addMissingAtoms().
  • What happens: PDBFixer's addMissingAtoms also inserts every residue in fixer.missingResidues, placed with stub geometry. The code never clears that list (fixer.missingResidues = {}) and never logs it.
  • Trigger: a PDB whose SEQRES lists residues that have no coordinates, such as disordered tails, tags or loops.
  • Impact:
    • The built protein differs from the deposited structure, and the user gets no warning.
    • It contradicts the test plan's stated limitation ("missing residues are detected but not modeled").
    • build_proteinbox has no RMSD or clash check after relaxation.
    • Neither 1aki nor 1zni exercises this.

Finding 2: no real end-to-end rerun after the 47c2298 refactor (completeness-checker, 85)

  • Scope of the refactor: 47c2298 (−2741/+1487 lines):
    • merges the builder bodies into shared helpers;
    • removes the second bundle_forcefield_into_topology call;
    • replaces the volume-compression loop with an equivalent-density target;
    • shares one mdp directory through template_dir.
  • What is missing: the last real e2e evidence comes from before it. The maintainer's own progress note lists "Re-run one real proteinbox and one protein_mixedbox e2e build before merge" as still to do.
  • What the static review found: code-reviewers A and B found no regressions. They checked:
    • the fixed-box density conversion, with a numeric OpenMM round trip;
    • hash stability for every example, against both base and d0c0f19;
    • that the merged [ atoms ] parser gives the same results as the old ones.

Finding 3: the _reimage_and_verify RMSD gate is untested, and the rotation-invariance test was deleted (test-reviewer, 85)

  • Where: mdfactory/build.py around line 996.
  • What changed: 47c2298 replaced _protein_internal_rmsd with an inline rmsd(center=True, superposition=True) call. It deleted test_internal_rmsd_ignores_rigid_rotation_and_translation without a replacement.
  • Impact: dropping superposition=True or comparing the wrong slice would still pass the suite. The gate would then reject healthy builds or accept distorted ones.

Finding 4: the fixed-box→density conversion and the argument forwarding in compress_box are untested (test-reviewer, 85)

  • What is unchecked: openmm_utils.compress_box:
    • converts ų to nm³ and derives the density target from the total mass;
    • forwards rigid_indices and require_convergence;
    • skips the residue wrap when protein_indices is set.
  • History: round 1 finding 6 raised this and it was deferred. The refactor now routes both sizing modes through this arithmetic.

Finding 5: the LNP leniency flag is still untested (test-reviewer, 84)

  • What is unguarded: create_lnp_core → require_target_density=False → require_convergence.
  • Impact: dropping it would bring back e2e blocking bug 2, and only a 15–24 minute GPU build would catch that.
  • History: this repeats round 1 finding 5.

Finding 6: the shared protein build helpers and the report bookkeeping are untested (test-reviewer, 82)

  • Helpers: _relax_protein_system and _finalize_protein_build.
  • What is unchecked:
    • the allow_net_charge=not neutralize wiring, which is the fix for e2e bug 4;
    • cleanup running only on success;
    • pressure=None for a fixed box;
    • final_input_counts alignment and the water count reduced by n_replaced;
    • the keys the build_metadata.json writer emits.
  • Also: if em.mdp is missing, grompp validation is skipped silently.
  • History: this repeats round 1 findings 7 and 13.

Finding 7: test_system_chemistry_reads_protein_mixedbox_build_metadata is still deselected by -k "not build" (test-reviewer, 92)

  • Where: mdfactory/tests/test_system_chemistry.py:154. The rename pass in 47c2298 missed it.
  • Why it matters: it is the only test asserting that the reader maps final_count and achieved_concentration_molar → achieved_concentration.
  • Also deselected (pre-existing): the whole test_build.py module, including the ionize_solvated_system zero-ion and one-sided tests.

Finding 8: the consolidated [ atoms ] parser is tested for the charge column on a single itp only (test-reviewer, 82)

  • Where: setup/protein.py _sum_atoms_column.
  • What is untested:
    • the mass column;
    • following local #include "topol_Protein_chain_*.itp", which is the real multi-chain layout;
    • skipping .ff/ includes.
  • Impact: if include-following breaks, the count_density box size and the protein_excluded concentrations come out wrong, and grompp doesn't catch it.

Finding 9: test_update_topology_molecules cannot detect swapped ion counts (test-reviewer, 80)

  • Where: mdfactory/tests/test_proteinbox.py:346.
  • The weakness: it checks "NA" in content and "10" in content separately, never as name/count pairs and never in order. So NA 18 / CL 10 would pass.
  • Why it matters: this is the area of the _ion_counts merge fix.

Finding 10: the RMSD-restraint test passes for a force that is identically zero, and relax_with_protein_restraints is never run in a test (test-reviewer, 85)

  • The weak test: it checks translation invariance only. It runs on the default (possibly single-precision) platform with abs=1e-8.
  • History: this repeats round 1 findings 12 and 13.

Suggestions

Finding 11: the PR body is out of date in several places (completeness-checker, 95)

  • compress_box: it says overshoot "beyond 2%" raises. The code rescales and raises only above 15%, and the LNP core only warns.
  • Test counts: it says "2998 passed, 25 failed"; the current count is 3099 passed, 1 skipped.
  • Status: it says "Draft: not yet validated end-to-end" and "blocked by the expired CGenFF license".
  • Test plan: it repeats the 2% line, and its "identical copies" mdp limitation is now resolved by template_dir.

Finding 12: test-plan checkboxes are still unticked for items the progress comments report as passing (completeness-checker, 90)

  • neutralize: false with concentration: 0;
  • the ethanol examples;
  • the CGenFF/Feb-2026 atom-type merge;
  • bilayer/lnp unaffected.

Finding 13: the example inputs need 1aki.pdb and 1zni.pdb, which are not in the repo (completeness-checker, 80). The test plan leaves open whether to ship them or document the download.

Finding 14: GromacsProteinParameterSet.forcefield and .water_model are written but never read (simplifier, 82). See models/parametrization.py:197-198 and setup/protein.py:724-725.

Finding 15: compress_box(overshoot_tolerance=0.15) only passes the value through (simplifier, 80). No caller overrides it, and it duplicates the finish_compression default.

Finding 16: the cutoff parameter on _check_no_protein_clashes and _check_no_solution_clashes is never overridden (simplifier, 80).

Strengths

  • Refactor (static review): neither code-reviewer found a regression from d0c0f19 to HEAD.
    • All four e2e fixes are still wired: finish_compression/rescale_to_volume, LNP require_target_density=False, topology_atom_names, and validate_with_grompp(allow_net_charge).
    • Every verification step is kept: re-image, RMSD, protein–solution clash and solution–solution clash, after both compression and relaxation.
  • Hash stability: hashes and model_dump_json are byte-identical to base for the six mixedbox, bilayer and LNP examples, and to d0c0f19 for the protein examples.
  • Fixed-box conversion: the conversion to an equivalent density round-trips exactly (373.248 nm³ → 0.0549215 g/mL → 373.248 nm³). The protein mass is part of the total mass.
  • Merged [ atoms ] parser: it gives the same results as the old parsers on a synthetic topology with comments, #ifdef POSRES, a local itp and a .ff/ include.
  • template_dir: it defaults to system_type, so existing types resolve their schedules as before.
  • [ molecules ] order: it matches the coordinates (protein, packed species, NA, CL), and ion counts are looked up by residue name.
  • grompp warning filter: with -maxwarn 1, it still rejects any single warning that is not the net-charge one.

Agents run: code-reviewer ×2 (split: build/openmm/setup_utilities; models/setup/protein/topology/CLI/push/analysis), error-auditor, test-reviewer, completeness-checker, simplifier


Reviewed by mach6

@maxscheurer

Copy link
Copy Markdown
Collaborator Author

Review Assessment

Findings comment: #59 (comment)

Review round: 2. Assessed by: independent-assessor.

Scope basis:

  • Issue 13's Phase-1 pipeline. Step 1 is "Clean PDB (PDBFixer)".
  • The PR body.
  • The maintainer test plan, including "Known limitations / follow-ups (not in this PR)". That section says "missing residues are detected but not modeled" and excludes Amber.
  • The maintainer's 47c2298 progress note: "Re-run one real proteinbox and one protein_mixedbox e2e build before merge."

Round 1 merge blockers: findings 1 and 2 are both verified fixed at 47c2298.

Classifications

Finding Classification Reasoning
1: clean_pdb quietly adds missing residues from SEQRES merge blocker Factual: setup/protein.py:519-523 runs findMissingResidues() and then addMissingAtoms(). It never clears fixer.missingResidues. In the installed PDBFixer, findMissingResidues records internal gaps and both terminal tails, and _addAtomsToTopology inserts them. PDBFixer's own CLI clears the dict with missingResidues = {} when residues are not wanted. I reproduced it: a PDB with SEQRES GLY GLY ALA ALA ALA SER SER and only ALA 3–5 resolved comes out of clean_pdb as all 7 residues, with no warning. The code has been there since the first proteinbox commit, so it is new relative to develop but not a regression from the refactor. Scope: It is issue step 1, and it does the opposite of the maintainer's documented contract ("detected but not modeled"). Practical: It affects any user with a typical RCSB PDB that has disordered tails, tags or loops. They silently simulate a different sequence with a different net charge, so the ion counts change too. Terminal tails are built as straight extensions, which inflates the proteinbox cube. In protein_mixedbox fixed_box they cause a misleading "box cannot hold the protein" error. No later check catches it, because grompp, the RMSD/clash gates and the restraint references all act on the already-modified protein.
2: no real e2e rerun after the 47c2298 refactor merge blocker Factual: 47c2298 (32 files, +1487/−2741) rewrites protein build orchestration. It also changes shared compress_box code used by mixedbox and LNP: the volume loop is removed and the fixed-box target is now converted to a density. The last real e2e evidence is for d0c0f19. Scope: The maintainer made the rerun an explicit pre-merge item. Practical: Unit tests mock gmx, pdb2gmx, packmol and OpenMM, and the maintainer notes the refactored builders have no direct tests. All four earlier e2e bugs were in exactly this orchestration code. Static review found no regression, which lowers the risk but does not show that the builds run.
3: _reimage_and_verify RMSD gate untested deferred (useful follow-up) Factual: True. The code (rmsd(center=True, superposition=True)) is correct. Scope/Practical: This is a missing test. A regression would fail builds loudly, and the e2e rerun in finding 2 exercises the gate.
4: compress_box density conversion and forwarding untested deferred (useful follow-up) Factual: True. The arithmetic is correct by inspection and by a numeric round trip. Scope/Practical: This is a missing test (round 1 finding 6). The fixed-box e2e build in finding 2 checks the exact final volume.
5: LNP leniency untested deferred (useful follow-up) Factual: True. Scope/Practical: It repeats round 1 finding 5, and the lines did not change in 47c2298. A regression fails loudly, and finish_compression's warn/raise path is tested.
6: protein helpers and bookkeeping untested; grompp skipped if em.mdp is missing deferred (useful follow-up) Factual: True (build.py:836-844 if em_mdp.is_file()). Scope/Practical: This is a missing test. test_protein_mixedbox_run_files_on_disk asserts that the mdp files resolve. Making a missing em.mdp raise would be cheap hardening.
7: system_chemistry build-metadata test deselected by -k "not build" deferred (trivial follow-up) Factual: True. Practical: The writer and reader keys currently agree (build.py:1278-1279 and analysis/utils.py:477-478). The fix is a one-line rename, and nothing is broken today.
8: _sum_atoms_column tested only for the charge column deferred Factual: True. Practical: The parser gives the same results as the old ones, include-following is covered by static review, and the e2e multi-chain insulin builds gave the correct net charge.
9: test_update_topology_molecules can't detect swapped counts deferred Factual: True. Practical: The function only writes f-strings, and _ion_counts has its own tests. A NA/CL swap gives mismatched atom names, which grompp -maxwarn 0 rejects.
10: weak RMSD-restraint test; relax_with_protein_restraints never run in a test deferred Factual: True. Practical: It repeats round 1 findings 12 and 13. The e2e RMSD of 0.26–0.32 Å shows the restraint holds.
11: PR body out of date nitpick Factual: True: the 2% line, the test counts, "Draft: not validated", the CGenFF license and the "identical copies" limitation. Update it before merge, including the missing-residue limitation from finding 1.
12: test-plan checkboxes not ticked nitpick Housekeeping. The progress comments report those items as passing.
13: examples need PDB files that are not in the repo deferred (useful follow-up) Factual: True. Practical: The user gets a clear file-not-found error. A README line or a download note is enough.
14: unused GromacsProteinParameterSet.forcefield/.water_model nitpick Factual: True. Harmless metadata.
15: compress_box(overshoot_tolerance) only passes the value through nitpick Factual: True. Style.
16: unused cutoff parameter on the clash checks nitpick Factual: True. Style.

Action Plan

  1. Finding 1: make clean_pdb's missing-residue handling deliberate, visible and documented.
    • Log: after findMissingResidues(), log a WARNING listing each insertion: chain, position, residue names, and whether it is terminal or internal.
    • Recommended policy: drop the terminal-tail entries, the keys at index 0 and index == the number of chain residues, so that tails are not built. Then either model internal gaps, with the warning, or raise on them.
    • Avoid clearing missingResidues completely. Internal gaps would then reach pdb2gmx unfilled, and pdb2gmx would bond across them.
    • Docs: update the clean_pdb docstring and the test-plan limitation text.
    • Test: add a fast unit test with real PDBFixer, using a SEQRES snippet like the one above. Keep "build" out of its name.
  2. Finding 2: run the real e2e builds on the final commit, after item 1, and post the results with the commit SHA.
    • proteinbox: lysozyme_charmm36m.yaml completes, and the in-build grompp validation runs and passes.
    • protein_mixedbox: lysozyme_fixed_box.yaml with ethanol (CGenFF). Check that:
      • the final box is 373.248 nm³;
      • RMSD is below 0.5 Å after both compression and relaxation;
      • grompp -maxwarn 0 passes;
      • build_metadata.json shows about 1.0 M ethanol.
    • mixedbox: one plain example, because compress_box is shared and was restructured. An LNP rerun is optional.

Assessment by mach6

@maxscheurer

Copy link
Copy Markdown
Collaborator Author

Progress Update: round 2 merge blockers (findings 1 and 2)

Finding 1: clean_pdb built residues that are missing from the coordinates

Fix (mdfactory/setup/protein.py): a new helper, _restrict_missing_residues_to_internal_gaps, runs right after PDBFixer's findMissingResidues().

  • Terminal tails: unresolved N- and C-terminal residues (tags, disordered termini) are dropped and not built. A warning lists the residues and the first or last resolved residue.
  • Internal gaps: these are still modeled by PDBFixer, and a warning names the residues and the two flanking residues. Leaving an internal gap open would make pdb2gmx bond straight across it.
  • Docstring: clean_pdb now describes this behavior.
  • Known-limitations line: "missing residues are detected but not modeled" in the test plan should now read "unresolved terminal residues are not modeled; internal gaps are modeled by PDBFixer; both are logged as warnings."

Test: TestCleanPdb.test_models_internal_gaps_but_not_terminal_tails runs real PDBFixer on a PDB whose SEQRES lists 9 residues, of which 4 have coordinates.

  • It expects residues 3–7, with the internal GLY 6 modeled, and warnings for the N-terminal, C-terminal and internal segments.
  • Without the fix it fails, because all 9 residues get built.

Checks: pytest -k "not build" mdfactory on a cpu node gives 3100 passed and 1 skipped. The pre-commit hooks pass.

Finding 2: real end-to-end runs on the refactored code

  • Where and what: the runs were on the gpu partition (L40S, gromacs/2024.3-gpu, gmx_mpi). The code was 47c2298 plus the fix above, identical to this commit.
  • Validation: every build ran the in-build grompp validation, and every protein build passed it with -maxwarn 0.
  • Simulation: each protein build then ran EM and 100 ps NVT with no LINCS or SETTLE warnings.
Build Result
proteinbox lysozyme_charmm36m.yaml Built in 25 s. [ molecules ]: protein, SOL 8424, NA 23, CL 31. Net charge 0. HIS15 → HSE. The protein is whole (max CA–CA 4.02 Å). EM Fmax 945.
protein_mixedbox lysozyme_fixed_box.yaml (ethanol, CGenFF) Built in 47 s. Compression overshot by 1.48% and was rescaled to exactly 373.248 nm³. RMSD was 0.28 Å after compression, and the relaxation RMSD check (< 0.5 Å) also passed. build_metadata.json: ETH 214 → 0.999 M, SOL 11716 → 54.67 M, NA 32 / CL 40. EM Fmax 858.
mixedbox water_ethanol_cgenff.yaml Built in 19 s. The density target was reached at step 15800, and a 0.17% overshoot was rescaled. The pre-existing atom-name mismatch on the plain mixedbox path is still there; this path has no grompp check.
proteinbox 1STN (finding 1, terminal tails) SEQRES covers 1–149; the coordinates cover 6–141. Both tails (5 and 8 residues) were reported and not built, giving 136 residues. Before the fix these 13 residues would have been built as extended stubs.
proteinbox 1aki with residues 101–103 deleted (finding 1, internal gap) The log reports "modeling 3 missing internal residue(s) ASP-GLY-ASN between SER100 and GLY104". The result has 129 residues, max CA–CA 3.97 Å. EM Fmax 920.

New bug, not fixed: chains that start with MET fail in pdb2gmx

The 5DFR build (DHFR, whose first residue is MET1) fails with atom C1 not found in buiding block 1MET while combining tdb and rtp.

  • Cause: the CHARMM36 Feb-2026 bundle's ethers.n.tdb defines a MET1 terminus, a methyl-ether cap. pdb2gmx matches termini on the first three characters of the residue name, so it takes MET1 as a MET-specific N-terminus.
  • C-terminal MET: ethers.c.tdb has a MET2 entry, which probably breaks chains that end with MET in the same way. I have not tested that.
  • Origin: this goes back to the original proteinbox work and is not caused by the changes here.
  • Impact: proteins whose sequence starts with Met are common, so this needs a decision: fix it in this PR, or file a follow-up issue.

Commit: b33adf3


Progress tracked by mach6

@gregorweiss gregorweiss left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

awesome new feature!

@maxscheurer

Copy link
Copy Markdown
Collaborator Author

Unverified Review Candidates — Pending Assessment

Review round: 3
Reviewed commit: eed1131

These are unverified candidates. Severity reflects reviewer confidence; do not treat any item as a merge blocker until the assessment comment is posted.

Scope of this round. This is a composite review. The target is the whole stack (PR 59, PR 63, PR 65, PR 67) as one diff against develop, taken at the stack tip eed1131 (branch feature/issue-62-generate-mdp-run-schedules), not PR 59's own head. It is posted on PR 59 because that is the PR that merges into develop. The maintainer asked for extra attention to legacy or backward-compatibility code, fallbacks, aliases and duplication. Those are tagged [COMPAT], [DEAD], [OVER] and [DUP], and they are collected in the index below.

Round 2 merge blockers on PR 59, re-checked:

  • Round 2 finding 1 (clean_pdb silently modeling SEQRES-missing residues and terminal tails): fixed. _restrict_missing_residues_to_internal_gaps drops terminal entries with a logged warning and keeps only internal gaps, before findMissingAtoms().
  • Round 2 finding 2 (no real end-to-end run after the refactor): fixed. The end-to-end runs are on the PR threads, and a further full build, simulate and analysis run through SLURM was done on this tip for the docs and examples.

No leftover legacy names in source. A grep of the final tree found no remaining reference to parametrize_smirnoff_gromacs, parametrization_config, CgenffConfig, SmirnoffConfig, CharmmConfig, run_schedules.yaml, df_models_from_input_csv or charmm_native_atom_names. Mentions remain only in a docs migration table (finding 16).

Index of requested categories

  • Legacy / compat / fallback: findings 6 (partly), 10, 16, 17, 18, 19, 32
  • Dead or unused code: finding 20
  • Over-engineering: finding 21
  • Duplication: findings 11, 17, 22 to 28

Critical

None.

Important

Finding 1: A build whose grompp validation fails is left looking simulation-ready (error-auditor, 85) [BUG][SILENT]

  • Where: mdfactory/build.py (u_relaxed.atoms.write("system.pdb") runs before _finalize_protein_build, which calls write_mdp_files then validate_with_grompp); mdfactory/cli.py discovery treats any directory with system.pdb as ready.
  • Trigger: validate_with_grompp raises (net charge, topology warning). The directory keeps system.pdb, topology.top, the mdp files and <hash>.yaml.
  • Impact: mdfactory simulate picks it up. The stage maxwarn budgets (1/1/2/2) tolerate the very warning that failed validation, so the run proceeds. With a Parsl build the failure only shows in the summary, and the CLI exits 0 (finding 2).
  • Suggested fix: write system.pdb last, or write a completion marker after validation and have discovery require it.

Finding 2: build (and simulate) exit 0 when builds or stages fail (error-auditor, 90) [SILENT]

  • Where: mdfactory/cli.py _report_build_results (logs failed, no sys.exit(1)); _report_simulation_results is the same.
  • Impact: scripts and CI see success. The stack adds new failure causes (version gate, grompp validation, chain mismatch, compression non-convergence). apply-settings already exits 1 on refused directories. The pipeline example script uses set -e and does not stop on a failed simulate.

Finding 3: Generated mdps differ from the deleted static templates for existing mixedbox, bilayer and lnp, with no token-level guard (code-reviewer, 90; test-reviewer, 82)

  • Where: mdfactory/run_schedules/gromacs.py (NONBONDED, DISPCORR, render_stage); mdfactory/tests/test_run_schedules.py:220-240, 325.
  • Diff against git show origin/develop:.../mixedbox/*.mdp and bilayer/*.mdp:
    • Identical: integrator, dt, lengths and intervals, LINCS and constraints, V-rescale tau-t, barostat per type, refcoord-scaling, gen-vel, default production lengths.
    • charmm family: old plain cutoff at 1.0 nm becomes 1.2 nm with Force-switch from 1.0 nm.
    • openff family: old plain cutoff at 1.0 nm becomes 0.9 nm with Potential-switch from 0.8 nm. Bilayer DispCorr changes from no to EnerPres, because the old file was chosen by system type and not by force field.
    • em.mdp now carries the family nonbonded block and DispCorr; ns_type = grid is gone.
    • lnp: the old code looked in a lnp/ directory that does not exist and copied nothing; lnp builds now get the mixedbox schedule.
  • Why it is a candidate: the changes are in line with the issue 62 decisions, but they alter the physics of every existing system and should be recorded in the PR description. Only the protein schedules have a full "reproduces the static files" table. For mixedbox, bilayer and lnp, nothing pins tau-p, compressibility, nstlist, the constraints block, PME order and spacing, minimizer settings or bilayer nsteps. A later edit to COUPLING or render_stage would change existing systems with no test failing.

Finding 4: CI deselects, by name, the only tests that run the real mixedbox and bilayer builds, real grompp, and the orchestration CLI (test-reviewer, 92; completeness-checker, 90)

  • Where: pytest -k "not build" matches file and test names. Of 3609 collected tests, 3524 run. Deselected: all of test_build.py (the only tests that run build_mixedbox and build_bilayer and check they call write_mdp_files and write the CHARMM schedule), test_amber_gmx_build.py (7 tests; also skip without gmx 2026), 20 tests in test_orchestration_build.py, 19 in test_orchestration_cli.py, 15 in test_orchestration_simulate.py, and test_system_chemistry_reads_protein_mixedbox_build_metadata (test_system_chemistry.py:154, round 2 finding 7, still open; the only test of the final_count and achieved_concentration_molar mapping).
  • Impact: removing write_mdp_files(inp) from a builder, or breaking family-to-mdp wiring, keeps CI green.
  • Suggested fix: rename or mark these tests so the CI selection runs the CPU-only ones.

Finding 5: sync push --csv and check-csv cannot take the settings protocol file, so hashes can mismatch (code-reviewer, 82)

  • Where: mdfactory/utils/push.py:110-124 and mdfactory/cli.py:865 call load_build_input_csv(csv) without the protocol file; build and prepare-build fold it into the hash.
  • Trigger: a CSV built with --simulation-settings protocol.yaml where the protocol sets temperature: 310; then sync push --csv.
  • Impact: the CSV-side hash is computed at 300 K and no build folder matches. This was a deferred round 1 follow-up on PR 67 and is documented in a docs callout.

Finding 6: Cache and download integrity is not enforced (error-auditor, 80; code-reviewer, 80) [BUG][SILENT] [COMPAT] for the third bullet

  • CGenFF cache: parametrize.py:208-241 reuses a cache entry when both files exist (if itp_path.is_file() and prm_path.is_file()) and writes them with plain open(..., "w"). A kill or disk-full during the second write leaves a truncated file that every later build reuses. The OpenFF cache already uses a .complete marker written last.
  • Force-field download: setup/protein.py:543-562 returns early when the target directory exists and extracts straight into the final location, so an interrupted extraction is trusted afterwards.
  • Optional hash: download_forcefield does expected_sha256 = entry.get("sha256") then if expected_sha256:. A registry entry without a hash silently skips verification. Both current entries have one.

Finding 7: Charge and mass sums silently skip bad lines and missing includes (error-auditor, 80) [SILENT]

  • Where: mdfactory/setup/protein.py:956-983 (_sum_atoms_column): if itp_path.is_file() skips a missing chain or posre include; except ValueError: pass skips a malformed charge or mass; round(sum) never checks closeness to an integer.
  • Impact: the charge drives neutralization and the mass drives the density-based box edge. A wrong value silently gives a wrongly neutralized or sized system. With neutralize: false, the tolerated net-charge warning hides it.

Finding 8: Close cysteine pairs are silently left reduced when disulfide_bonds is empty (error-auditor, 80) [SILENT]

  • Where: setup/protein.py:890 (if disulfide_bonds:) and _build_disulfide_prompt_input.
  • Trigger: a PDB with SG–SG pairs at or below 2.2 Å and no declared bonds.
  • Impact: a real disulfide-bonded protein is built with free thiols and no warning. Candidate detection only runs when bonds are requested.

Finding 9: apply-settings is not robust to bad records and writes non-atomically (error-auditor, 80) [BUG]

  • Where: run_schedules/__init__.py:84-87, 133-138, 195. yaml.safe_load of the record is outside the _Refused conversion and yaml.YAMLError and OSError are not caught. write_mdp_files then write_build_record are separate writes.
  • Impact: one corrupt record aborts the batch after earlier directories were already rewritten, with no summary. A failure between the two writes leaves mdps newer than the record. This was a deferred PR 67 round 1 follow-up.

Finding 10: Old-schema build directories silently drop out of discovery (error-auditor, 80) [SILENT][COMPAT]

  • Where: analysis/utils.py:86-89, 106-107 (except Exception: continue) and utils/push.py:207-209 (except (FileNotFoundError, ValueError): continue).
  • Why now: BuildInput is extra="forbid" with a new schema, so every pre-stack <hash>.yaml fails validation, and these unchanged blocks swallow the break. Old directories vanish from status, analysis and push with no message. A clean cut should fail loudly with "build record is from an older schema".

Finding 11: The same ionization.concentration resolves to different salt counts in the two protein types (code-reviewer A1, 80) [BUG][DUP]

  • Where: mdfactory/build.py:689-696 in ionize_solvated_system: if solvent_volume_a3 is not None: n_ions = int(round(c * V * N_A)) else np.ceil(c * n_water / 55.55).
  • Impact: proteinbox uses the water-count basis with ceil; protein_mixedbox uses the protein-excluded volume with round. One rule with an explicit volume would be consistent.

Finding 12: PR housekeeping before merge (completeness-checker, 90)

  • PR 56 is obsolete and breaks CI after the stack lands. It adds test_packaging.py, which hard-codes the eight deleted gromacs/{bilayer,mixedbox}/*.mdp wheel members and smoke-tests RunScheduleManager, and wires it into .github/workflows/ci.yml. It merges textually clean. Issue 55 stays open because only PR 56 says Closes.
  • Closing keywords: PR 67 should say Closes for issue 55 and issue 66 (the protein docs delivered in eed1131). PR 59 says "Implements" for issue 13, so it never closes. Issue 13 itself still promises amber99sb-ildn and names the deleted Pdb2gmxConfig.
  • Stale PR 59 body (requested in both earlier assessments, never done): it still says compress_box raises above 2% overshoot (now rescale, raise above 15%), "2998 passed, 25 failed", "Draft: not yet validated end-to-end", "blocked by the expired CGenFF license", and CharmmConfig. Four boxes in its test-plan comment are unticked although the progress comments report them passing. The plan-comment acceptance boxes on PR 63, PR 65 and PR 67 are unticked, and the PR 63 "identical hash" criterion is superseded by the PR 65 clean break.
  • PR 67 body does not mention apply-settings, the docs, the examples or the package decision (issue 62 acceptance criterion 4 says the issue or PR documents it; only the plan comment does).

Finding 13: The 20-system tutorial CSV is delivered in the new schema but fails validation on every row (completeness-checker, 95)

  • Where: docs/public/tutorials/mixedbox_chain_length_series_20.csv and tutorial-mixedbox-campaign.mdx:65, 278, 323 use the resname MEOH (4 characters; the limit is 3).
  • Evidence: load_build_input_csv returns 0 models and 20 errors. Flagged as pre-existing in the PR 65 progress comment. The examples-parse test covers only examples/.

Finding 14: The CHARMM MET1 pdb2gmx bug is neither fixed nor tracked (completeness-checker, 82)

  • Round 2 asked for "fix here or file a follow-up". No issue exists and no docs mention it. Any chain starting with Met fails on the CHARMM route ("atom C1 not found in building block 1MET"); the C-terminal case is probably affected too.

Finding 15: Unguarded seams in tests (test-reviewer, 80 to 85)

  • run_build_from_dict (writes <hash>.yaml before dispatch) is only ever mocked; nothing asserts the record exists before the builder runs or that its hash equals BuildInput.hash. Without the record, apply-settings, simulate and analysis refuse the directory.
  • The LNP leniency (build_lnp to create_mixedbox_by_compression(require_target_density=False) to compress_box(require_convergence=False)) is tested only at the leaf finish_compression. build_lnp is not exercised at all, including its write_mdp_files call.
  • The shipped examples/simulation_settings/*.yaml protocols, examples/pipeline/protocol.yaml and slurm_gpu.yaml are never loaded by a test (the examples tests filter on simulation_type). The docs tell users to pass them to --simulation-settings and apply-settings.

Legacy, compatibility and fallback code (maintainer policy: none wanted)

Finding 16: Docs keep an old-key to new-key migration section (code-reviewer, 80; simplifier, 90) [COMPAT]

  • docs/content/docs/user-guide/parametrization.mdx:251-272, "Migrating from the old keys": a table mapping cgenff|smirnoff|pdb2gmx and each parametrization_config.* key to the new form, plus prose on old OpenFF cache directories and old build directories. The code is a clean cut; only this docs section is migration material.

Finding 17: Settings.DEFAULT_CONFIG is a hand-kept twin of _get_defaults() "for backward compat in tests" (code-reviewer, 85; simplifier, 80; error-auditor, 80) [COMPAT][DUP]

  • mdfactory/settings.py:85-89. The stack edited it to add the gromacs block. Only tests read it (conftest.py, test_csv_backend.py, test_push.py), and the two tables already disagree on FORCEFIELD_DIR.

Finding 18: LNP core compression keeps the old non-convergence behaviour behind a flag (code-reviewer A1, 80) [COMPAT]

  • build.py:97-107, 148-152, 394-399 and openmm_utils.py:180-184, 511-514. The old compress_box ignored non-convergence; the stack makes it raise and adds require_convergence=False whose only user is create_lnp_core. The comments say "keep the pre-existing behavior" and "the behavior LNP core compression relies on", and the miss is reported with a bare print("WARNING ...") that no logger captures.

Finding 19: Optional parameters and .get defaults that no production caller omits (code-reviewer A2a, 80; simplifier, 85) [COMPAT]

  • setup/protein.py: _apply_protonation_states(..., forcefield=None), _resolve_protonation_residue_name(..., forcefield: str | None) with if forcefield and forcefield.lower().startswith("charmm"), run_pdb2gmx(..., chains=None) then chains = chains or []. If a caller ever omits forcefield, the CHARMM residue aliasing is skipped without a signal.
  • topology_utilities.py:760-771: validate_cgenff_atomtype_compatibility(forcefield_types=None) with an ff_dir re-read fallback; the only production caller always passes forcefield_types.
  • parametrize.py:744: default out_path=Path("topology.top"); both callers pass it.

Finding 20: Dead or write-only code retained (simplifier, 80 to 95; code-reviewer, 85) [DEAD]

  • SimulationSettings.identity() (simulation_settings.py:304): no caller; BuildInput.hash uses IDENTITY_FIELDS directly.
  • SimulationType Literal (models/parametrization.py:27): never used.
  • GromacsProteinParameterSet.forcefield, .water_model, .chains (models/parametrization.py:469-471): read only by tests.
  • workflows.run_build_from_file (workflows.py:56): the CLI now calls run_build_from_dict; only test_proteinbox.py and SKILLS.md reference it. The diff added a resolve_protein_pdb_path import to workflows.py solely for it.
  • compress_box(overshoot_tolerance=0.15, restraint_k=...) and relax_with_protein_restraints(restraint_k=...): forwarded or never passed; finish_compression returns a value its only production caller discards.
  • Parameters that only take their default: cutoff on _check_no_protein_clashes and _check_no_solution_clashes, cutoff_angstrom on _build_disulfide_prompt_input, water_name on update_topology_molecules, tolerance on check_defaults_compatible.
  • The charmm36m-ljpme entry in FORCEFIELD_REGISTRY: unreachable, because CharmmParametrization rejects any protein name containing ljpme, yet the config wizard still downloads it.

Suggestions

Finding 21: Over-engineering (simplifier, 80 to 85) [OVER]

  • small_molecules_model ClassVar duplicates the field annotation in four spec classes; OpenFFParametrization.protein plus _check_no_protein only improve an error message that _FamilySpec.resolve already raises; provenance.json is required by is_complete but never read outside tests; the GROMACS version is formatted to a string, passed through grompp_kwargs and re-parsed in _build_gmx_version_check.

Finding 22: Hand-rolled GROMACS section and include walkers (simplifier, 92) [DUP]

  • Nine copies of "strip, skip comments, track [ name ], split" across topology_utilities.py (collect_forcefield_atomtype_definitions, read_itp_atoms, write_itp_with_atom_names, rename_atomtypes ×2, read_topology_defaults), setup/protein.py (_parse_protein_molecule_names, _sum_atoms_column) and parametrize.py (_read_include_layout); four copies of #include parsing with the .ff/ test. Two small generators in topology_utilities.py would replace them.

Finding 23: Duplicate atom-type loop in two validators (simplifier, 90) [DUP]

  • topology_utilities.py:772-779 and :918-925 (validate_cgenff_atomtype_compatibility, check_atomtype_collisions).

Finding 24: Native water and ion tables repeated (simplifier, 80 to 85; code-reviewer, 80) [DUP]

  • The O, [Na+], [Cl-] to SOL, NA, CL mapping lives in five places (models/species.py, parametrize.py, utils/setup_utilities.py, build.py, setup/protein.py). The native GromacsSingleMoleculeParameterSet is built twice (parametrize.py:164-188 and :431-458) with only the ion file name different.

Finding 25: Repeated type sets and mappings (simplifier, 80 to 95; code-reviewer, 85) [DUP]

  • PROTEIN_PDB_SIMULATION_TYPES (prepare.py:17) and the literal tuple in models/input.py:250 duplicate PROTEIN_SIMULATION_TYPES; the simulation-type list is spelled five times; STAGE_FILES duplicates StageSpec.mdp_file and _stage_outputs bridges them with a lazy import and a linear search; the nonbonded and coupling tables are restated in parametrization.mdx, simulation-settings.mdx and the tests; 1e-27, 1e24 and 6.022e23 (models/composition.py:225) are repeated against N_AVOGADRO; relative pdb_path resolution exists twice in prepare.py; build_metadata.json is read and the protein row injected in two places (Simulation.metadata/BuildInput.metadata and system_chemistry).

Finding 26: Duplication across the builders (code-reviewer A1, 80 to 85; simplifier, 80) [DUP]

  • Parametrizer plus topology-function resolution copy-pasted in create_mixedbox_by_compression and create_bilayer_by_compression; solvate, box, zero-water guard and water species repeated in build_bilayer, build_lnp and build_proteinbox (a TODO is already there); build_proteinbox hand-rolls protein staging that _stage_protein_topology factors out for protein_mixedbox; bounding-box centering and extent computed in four places; the 0.1 dilution idiom exists as a literal and as PROTEIN_MIXEDBOX_DILUTION_FACTOR; create_sphere keeps its own copy of the PDB/mmCIF round trip that _universe_from_simulation replaced elsewhere.

Finding 27: Test duplication and global-settings leakage (simplifier, 88; test-reviewer, 80) [DUP]

  • Spec rejection and frozen tests are repeated across test_parametrization_schema.py, test_proteinbox.py and test_protein_mixedbox.py; the WATER dict, EXAMPLES_DIR, protein-type constants and three native force-field stubs are copied between files; two CGenFF tests assign config.parameter_store on the global settings object without restoring it.

Finding 28: Brittle assertions (test-reviewer, 80)

  • Tests couple to pydantic's exact message text and layout (test_parametrization_schema.py:276-312, test_proteinbox.py:280) and to functools.partial internals (test_protein_mixedbox.py:918-929). exc.errors() and observable backend calls would carry the same contract.

Finding 29: Version-gate loose ends (test-reviewer, 80; error-auditor, 80)

  • The simulate dry-run path that passes min_gmx_version and logs "Requires GROMACS >= ..." is untested, as is a gmx_mpi-only PATH case. min_gmx_version_for_topology returns None when topology.top is missing and scans only top-level includes, so the simulate-time gate can skip.

Finding 30: Docs overclaim override-key validation (completeness-checker, 85)

  • simulation-settings.mdx ("GROMACS overrides") says "grompp validates the result when the stage runs". An unknown key is only a grompp warning that stages with spare maxwarn absorb. PR 67 round 2 finding 2 asked at minimum for this sentence to be softened.

Finding 31: OpenMM nits (code-reviewer A1, 80)

  • compress_box still does app.PDBFile.writeFile(..., open("min.pdb", "w")) with an unclosed handle (fixed for create_sphere in the same diff). The RMSD restraint docs say "harmonic position restraints" that keep positions, but the force is a best-fit RMSD restraint that allows rigid motion, with energy k·Σd² rather than the conventional ½kΣd², so the effective constant is twice the documented restraint_k.

Finding 32: Remaining silent defaults (error-auditor, 80) [SILENT]

  • An empty or comment-only --simulation-settings file loads as {} (simulation_settings.py:392-395), so apply-settings reports "unchanged" and exits 0. system_chemistry reports None for resolved_count, final_count and achieved_concentration with no warning when build_metadata.json is missing or shorter than the species list, matched by position. clean_pdb removes heterogens (removeHeterogens(keepWater=False)) without logging which ligands, ions or cofactors were dropped. Pre-existing and outside this diff: analysis/utils.py:40, 59 ("Default to 'production' for backward compatibility") and utils/push.py:99 ("Fallback for folders without structure file").

Strengths

  • Clean cut in the code: no leftover old names in source (see above); BuildInput and every spec use extra="forbid"; family is required everywhere; the DB column holds only the family name.
  • Identity model: the hash covers BuildInput plus temperature and pressure only. A different temperature changes it, a different production length does not, and a <hash>.yaml round trip reloads to the same hash (checked directly), with sort_keys=False preserving dict order.
  • Caches: the OpenFF cache key includes force field, resolved charge method and format version; the completion marker is written last; locks live next to the directory so cleanup cannot remove its own lock.
  • Version gate: build-time and simulate-time checks both fail closed on an unparseable version; tar.extractall(..., filter="data") blocks traversal; the SHA-256 is verified for both registry entries.
  • Packaging: removing the hatch force-include for the mdp directory is safe now that the files are rendered in Python; no source reference to the deleted files remains.
  • Numerics: unit conversions (ų, nm³, L, g/cm³), the ion and water bookkeeping and the [ molecules ] order were verified; a positive concentration that rounds to zero molecules is rejected.
  • Verification evidence: the new test files from the CI selection pass (432 passed, 1 skipped for Packmol); a full CSV, SLURM build, simulate and analysis run completed for four protein systems on this tip, and five new example inputs built on a GPU node.

Agents run: code-reviewer ×3 (part B; part A1; part A2a), error-auditor, simplifier, test-reviewer, completeness-checker. Two further code-reviewer runs (part A, part A2b) were cancelled after stalling; their scope (setup_utilities.py, topology_utilities.py, orchestration/apps.py, orchestration/stages.py) was covered by the error-auditor and simplifier, by A2a for setup/protein.py, and by a direct read of the version-gate bash and the concentration and box-size maths, which found no defect.


Reviewed by mach6

@maxscheurer

Copy link
Copy Markdown
Collaborator Author

Review Assessment

Findings comment: #59 (comment)

Round 3, composite review of the stack tip eed1131. Three assessors ran in parallel: the independent assessor (verified every finding against the code), the developer's advocate (practical value of fixing each item before merge) and the devil's advocate (acceptance evidence). A finding is a merge blocker only if the independent assessor and the developer's advocate both find material practical impact. Neither did for any finding, so there are no merge blockers. Where the advocates added something concrete, it is noted in the reasoning.

Classifications

Finding Classification Reasoning
1: Failed grompp validation leaves a simulation-ready directory useful follow-up Factual: True for protein builds. system.pdb is written (build.py:882, 1233) before _finalize_protein_build runs mdp rendering and validate_with_grompp; discovery and simulate._missing_build_files need only system.pdb, topology.top and the mdps, and the stage maxwarn budgets (1/1/2/2) absorb a warning-level failure. Scope: Validation is new in the stack, but no criterion says a failed build must be undiscoverable, and the write-then-fail order exists on develop. Practical: Needs a Parsl batch where one system fails with a warning-level validation error, the failure is shown in the summary, and the user then runs simulate on the whole directory anyway. A hard grompp error would fail again. Delete system.pdb or write a marker on failure as a small separate change.
2: build and simulate exit 0 on failures deferred Factual: True, and identical on develop (_report_build_results existed there). Scope: Pre-existing; the stack only adds failure causes. Practical: Only the Parsl or SLURM path (sequential builds raise); failures are logged and summarised. A sys.exit(1) when failed > 0 is small but changes the CLI contract and the example scripts, so track it as its own issue.
3: Generated mdps differ from the deleted templates; no token-level guard useful follow-up Factual: Confirmed against git show origin/develop:. Old: 1.0 nm plain cutoff for both families, bilayer DispCorr = no, ns_type = grid. New: charmm 1.2 nm Force-switch from 1.0, sage 0.9 nm Potential-switch from 0.8, bilayer DispCorr follows the family. Integrator, LINCS, thermostat, barostat per type, PME and default lengths are unchanged. Scope: The change is the approved issue 62 decision. Practical: Nothing is broken; the physics of existing mixedbox and bilayer systems changed, which the PR 67 body should state. The devil's advocate notes the existing PR 67 regression report and grompp -po diff are real evidence for the change itself. The missing piece is a permanent guard: a frozen MIXEDBOX_SCHEDULE and BILAYER_SCHEDULE table mirroring test_protein_schedules_reproduce_the_static_files, with values taken from develop.
4: CI deselects build-named tests useful follow-up Factual: True (-k "not build", pyproject.toml:106). Scope: Pre-existing convention; no criterion requires CI coverage. Practical: Removing write_mdp_files(inp) from a builder would keep CI green. The devil's advocate points out that renaming would not help for test_build_mixedbox and test_build_bilayer, which need CGenFF; the useful fix is a mocked builder-wiring test for build_mixedbox, build_bilayer and build_lnp, plus renaming test_system_chemistry_reads_protein_mixedbox_build_metadata.
5: push --csv and check-csv ignore the protocol file deferred Factual: True. Scope: Round 1 deferral on PR 67, documented in a callout in simulation-settings.mdx. Practical: Only when a protocol sets temperature or pressure; fails loudly ("no folder found").
6: Cache and download integrity discarded observation Factual: The patterns exist (parametrize.py:208, early return on an existing ff directory, optional sha256). Scope: The CGenFF cache pattern and the early return predate the stack; both registry entries carry a hash. Practical: Needs a kill or full disk during a small write. If wanted later, the .complete marker already used by the OpenFF cache is the model.
7: _sum_atoms_column skips bad lines and missing includes discarded observation Factual: True (protein.py:956-983). Scope: No criterion. Practical: Input is pdb2gmx output; skipping .ff/ includes is deliberate, and grompp validation reads the same topology.
8: Close cysteine pairs left reduced when disulfide_bonds is empty discarded observation Factual: Detection runs only inside if disulfide_bonds:. Scope: Issue 13 asks for user-declared bonds; auto-detection would be a new feature. Practical: Bonds are an explicit documented input.
9: apply-settings batch not robust to bad records deferred Factual: yaml.safe_load is outside the _Refused conversion; the two writes are separate. Scope: Round 1 deferral on PR 67. Practical: Records are machine-written; re-running repairs it. Wrapping YAMLError and OSError is a trivial change.
10: Old-schema directories silently dropped from discovery deferred (maintainer decision) Factual: True, in blocks that predate the stack (analysis/utils.py:86-89, push.py:207-209). Scope: The stack's schema break is what makes every old record fail there, but it touched neither block, and the standing policy forbids special handling of old keys. Practical: The docs state old directories cannot be parsed. The devil's advocate asks for one test with a real old-schema record asserting the intended behaviour (loud failure or logged skip); the maintainer should decide which behaviour is intended.
11: Different salt-count rules in the two protein types discarded observation Factual: Both branches exist (build.py:689-696). Scope: Intentional: protein_mixedbox has an explicit concentration_volume_basis (documented); the other branch is the pre-existing behaviour. Practical: Differences are about one ion; unifying would change the ion counts that the end-to-end runs validated.
12: PR housekeeping (PR 56, closing keywords, stale bodies) useful follow-up Factual: PR 56 adds test_packaging.py (hard-codes the deleted gromacs/{bilayer,mixedbox}/*.mdp and RunScheduleManager) and wires it into CI; it merges textually clean, so merging it after the stack would fail its own CI. Issue 55 stays open and issue 66 is delivered but unreferenced. The PR 59 body still carries "Draft, not validated", the 2% overshoot rule, "2998 passed, 25 failed" and CharmmConfig. Scope: No acceptance criterion. Practical: Only a wrong merge record or an accidental PR 56 merge. Close PR 56, add Closes keywords to PR 67 for issues 55 and 66, refresh the PR 59 and PR 67 bodies, tick or supersede the plan checkboxes.
13: Tutorial CSV fails on every row useful follow-up Factual: Confirmed: 0 models, 20 errors ("Residue name must be less or equal than 3 characters long", MEOH). Scope: The 3-character limit is on develop, so the file was already broken; the stack only changed the parametrization column. Practical: Anyone following the tutorial hits 20 errors. A rename in the CSV and tutorial-mixedbox-campaign.mdx, plus adding docs/public/tutorials/*.csv to the examples-parse test, fixes it. Cheap and worth doing before merge, though it does not meet the blocker gates.
14: CHARMM MET1 pdb2gmx bug untracked useful follow-up Factual: No issue and no docs mention. Scope: A gap in the new CHARMM route, not a regression; round 2 asked for "fix or file". Practical: Fails loudly with a pdb2gmx error; the Amber route is unaffected. File the issue and add a troubleshooting line.
15: Unguarded test seams useful follow-up Factual: Plausible and matches the devil's advocate: run_build_from_dict is always mocked; build_lnp wiring is untested; the shipped example protocols are never loaded. Scope: Missing tests are not blockers on their own. Practical: Minimal evidence is three small tests: <hash>.yaml exists before dispatch and reloads to the same hash; the example protocols load through load_simulation_settings; a parametrized check that every renderer-owned key is rejected as an override.
16: Docs migration section useful follow-up Factual: True (parametrization.mdx:251-272). Scope: Contradicts the maintainer's clean-cut policy and was explicitly requested to be flagged; docs only. Practical: No harm. Delete the section and keep a one-line breaking-change note in the PR body.
17: DEFAULT_CONFIG "backward compat in tests" useful follow-up Factual: True; the comment is older, the stack added the gromacs block, and the two tables disagree on FORCEFIELD_DIR. Scope: Pre-existing attribute extended by the stack. Practical: Duplication only; swapping needs edits to three test modules.
18: LNP require_convergence=False useful follow-up (small) Factual: True; the comments say "pre-existing behavior" and the miss is a bare print. Scope: Not a regression; the flag has a real production caller (create_lnp_core), so removing it would break LNP builds. The developer's advocate rates this not compat. Practical: Reword the comments and use logger.warning instead of print.
19: Optional parameters no production caller omits nitpick Factual: True for forcefield and chains in protein.py. Scope: Not required. Practical: No harm today. Tighten when touching those functions.
20: Dead and write-only code useful follow-up Factual: Confirmed: SimulationSettings.identity(), SimulationType, run_build_from_file (only a test and SKILLS.md), write-only GromacsProteinParameterSet fields, and the unreachable charmm36m-ljpme entry that the wizard still downloads. Scope: Flagged at the maintainer's request, not required by an issue. Practical: Safe deletions in one cleanup PR.
21: Over-engineering nitpick Design preferences with no defect; re-opens validated code.
22: Hand-rolled section and include walkers useful follow-up Factual: Nine section loops and four include parsers exist. Scope: Duplication was requested to be flagged. Practical: The parsers work and are validated end to end; consolidate in a dedicated refactor, not in this stack.
23: Duplicate atom-type loop in two validators useful follow-up Factual: True. Scope/Practical: Small, low-risk cleanup; the two sites apply different rules to the same loop.
24: Native water and ion tables repeated useful follow-up Factual: Five copies of the SMILES to resname mapping and two native parameter-set constructors. Practical: Constants only; consolidate in the cleanup PR.
25: Repeated type sets, stage map, constants useful follow-up Factual: PROTEIN_PDB_SIMULATION_TYPES and the tuple at models/input.py:250 duplicate PROTEIN_SIMULATION_TYPES; the others are style. Practical: The type-set de-duplication is cheap and safe.
26: Duplication across the builders useful follow-up Factual: True; build.py carries an existing TODO. Practical: All five builders passed end-to-end runs; refactoring now would invalidate that evidence.
27: Test duplication and settings leakage useful follow-up Factual: True. Practical: Cosmetic apart from the two CGenFF tests that leak parameter_store; fix those with finding 17.
28: Brittle assertions nitpick Style point consistent with the maintainer's preference for contract tests; low priority.
29: Version-gate loose ends discarded observation Factual: min_gmx_version_for_topology returns None for a missing file. Scope/Practical: _missing_build_files already refuses a build without topology.top, and builds write the force-field include at the top level, so this is not reachable in supported use.
30: Docs overclaim override validation nitpick Factual: simulation-settings.mdx:124 says grompp validates the result; unknown keys are only a warning that spare maxwarn absorbs. Practical: Minor misreading risk; soften as requested in PR 67 round 2.
31: OpenMM nits nitpick The unclosed handle is pre-existing; the restraint constant wording is a docs point.
32: Remaining silent defaults discarded observation An empty settings file means "no change" and the user sees no change; heterogen removal is documented; the other items are outside this diff.

Action Plan

No merge blockers; there is nothing to implement before merge.

The assessors suggest these non-blocking items, cheapest first, if you want to tidy the stack before merging (none of them gates the merge):

  1. Merge hygiene (finding 12): close PR 56 and add Closes keywords for issues 55 and 66 on PR 67, refresh the PR 59 and PR 67 bodies, record the mdp physics changes from finding 3, and tick or supersede the stale plan checkboxes.
  2. Docs only: rename MEOH in the tutorial CSV and tutorial (finding 13), delete the migration section (finding 16), soften the override-validation sentence (finding 30).
  3. File issues for the follow-ups: the Met-start bug (14), exit codes (2), failed-validation directories (1), CI wiring and the mixedbox/bilayer schedule table (3, 4, 15), plus one cleanup issue for findings 17 to 20 and 22 to 27.

Assessment by mach6

@maxscheurer
maxscheurer merged commit ebfeb58 into develop Nov 10, 2026
5 checks passed
@maxscheurer
maxscheurer deleted the feature/protein-mixedbox branch November 10, 2026 16:40
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.

4 participants