Skip to content

0.6.5 — ligand parameterisation

This release is about preparing a small molecule: how it is protonated, how its torsions are parameterised, and which of its torsions REST2 heats. The umbrella release is postponed and is not in it, and MD-tools is, for now, a conformational-sampling toolkit: the alchemical pages and tutorials are on hold (see below). Three changes, each a decision about what the Hamiltonian is, so each changes results for some inputs; this note says which.

1. Open Babel is in openmm-env

openbabel 3.2.1 was installed from conda-forge into the shared openmm-env (the solver also moved openssl 3.6.4 → 3.6.5, and nothing else). It is listed in environment.yml. Open Babel was not in any environment before; the only copy on the machine was a 3.0.0 binary at /usr/local/bin/obabel that the environments cannot see.

2. A small molecule is protonated at a stated pH, by default 7.0

protonation:
  ph: 7.0                   # default; the same key the protein half uses
  ligand_method: openbabel  # default; `none` keeps the input's protomer

OBMol.AddHydrogens(False, True, pH) is what obabel -p calls, so the CLI and the library cannot disagree; the record carries the library version. It runs on the molecule as supplied — before the conformer is embedded, before the charges are calculated, before the package identity is computed — for kind: ligand and kind: peptide-like, for --parameterize, and for .smi, .sdf and .mol2 alike.

What changes: CC(=O)O is now built as acetate, an aliphatic amine as the ammonium, a guanidine as the guanidinium, at pH 7. What does not: paracetamol, benzene, tertiary amides — nothing in the table applies, and the molecule you supplied is used as the same object, so a docked pose keeps its hydrogens.

Limits, stated because a default is an action. It is a table of functional-group rules and not a pKa predictor. Measured with the installed version: a phosphate (OP(=O)(O)O) and an imidazole (c1c[nH]cn1) are left as written; a carboxylic acid is ionised and an amine protonated at 7.0, and reversed at pH 2 and 12. It never changes the heavy-atom graph or the number of specified stereo elements — a result that does is refused, which the tests exercise with a molecule that does — and a 3D input's heavy atoms stay within 1e-3 Å (measured 5e-5 Å, the molfile precision). Open Babel missing is an error naming ligand_method: none; it never falls back.

Refused by name: protonation.ligand_method on kind: peptide (PROPKA decides there) and on kind: complex (a mapped ligand's hydrogens are its package's). protonation.method, overrides and the histidine keys are still refused for a molecule.

What this does to existing work. A package is one exact chemical state, so a .smi that Open Babel now changes no longer matches a package built from the unchanged molecule. That package is not edited and is still reusable by parameters: <id> from a molecule that really is that protomer; a search simply finds nothing and parameterises the new species. Anything built from an input the model leaves alone is unchanged.

3. REST2 convention v4: the ω protection is for peptides

rest2-unscaled-torsions goes from v3 to v4. An ordinary amide C–N inside a non-protein residue is now scaled, and recorded under ligand_amide_scaled_bonds (as ureas already are, so that a decision never looks like a miss). A protein residue keeps the protection, and so does a residue named in the new scaler key peptide_like_residues.

To the scaler a solute.kind: peptide-like molecule is one non-protein residue read from an SDF, which is also what a ligand looks like. So peptide_like_residues is unset by default and read from the build record: build/built.log carries resolved_config.solute.kind, and for peptide-like the solute's non-protein residues keep the protection; any other kind names none. A stated list wins, an empty one included. With no readable record the answer is "none" and the scaler's ligand amides log line says so and where the list came from (peptide_like_residues_source in scaler.yaml). The build record for a peptide-like build is classified the same way.

The rule is "inside one residue". A bond between two residues is a peptide bond whatever either is called: a hydroxyproline in a protein chain still reads bond orders or refuses, as before. The first version of this change got that wrong — it scaled HYP's backbone — and test_proline_like_is_configurable_and_hydroxyproline_can_be_declared caught it.

To protect one ligand amide anyway: its central bond in unscaled_list.

The guanidinium is classified as a group (decided 2026-10-13). Its three C–N bonds are partly double whichever one a SMILES writes as the double bond, so that notation must not decide what REST2 heats. In a protein or peptide-like residue all three are protected, as a protein's ARG always was; in a ligand all three are scaled and recorded under guanidinium_scaled_bonds. Before this, the one bond the SDF wrote as double was protected and the other two were scaled.

v3 ladders are continued, not refused. The identity a ladder stores embeds the convention and is compared whole on resume and extension, so a plain bump would have made every in-flight v3 ladder unresumable once this is installed. rest2.scaler.rest2_implementations_agree is the one named compatibility branch: recorded v3, current v4, every other field equal. The states are files, so the Hamiltonian a v3 ladder integrates is unchanged, and ReplicaRun.adopt_recorded_convention keeps the continuation recording v3 rather than relabelling states it did not build.

Tutorials. The paracetamol REST2 and AIS pages and the TYK2 REST2 page record v3 runs and say so in a notice at the top; their numbers are unchanged and were not regenerated, because that needs the GPU runs.

4. A bespoke-torsion .offxml enters a package

solute:
  torsion_refinement: mace-omol-bespokefit      # none | bespokefit | mace-omol-bespokefit | mace-off23-bespokefit
  refinement_offxml: ligand_sage_mace.offxml    # the fitted force field

BespokeFit refits a molecule's torsions against a reference and writes one .offxml. build-top does not run the fit: it takes hours and GPUs, and BespokeFit, QCEngine, Psi4 and the MACE models are not in this package's environment. The user installs them separately and names the file.

The file is verified: it must hold every proper torsion of the base force field unchanged and add at least one, so a stock file and a different base are both refused. The package records the method and the file's sha256 and a search treats them as identity, so a refined package never matches a stock request or another refinement. A method with no file is an error, never a quiet fall back to Sage's torsions. The package's AM1-BCC charges are still generated once and still win over the file's own library charges. A package written before this reads as stock, byte for byte.

Shown on the ligand-parameterization tutorial: paracetamol with Sage 2.2.1's torsions and with torsions refitted to MACE-OMOL-0. Bespoke torsion fit shows how to install BespokeFit and MACE-OMOL and run the fit, with the scripts it ships; the driver ran for paracetamol on two GPUs (success, 23 minutes) and reproduces the earlier fit's six torsions to within 0.003 kcal/mol.

5. Tutorials, and what is on hold

  • New: Ligand parameterization, with paracetamol as a ligand and c(RGDfV) as a peptide-like molecule.
  • On hold: the alchemical method pages and the two alchemical tutorials (ethanol hydration, ethane → chloroethane) moved to On hold with the plans that wait on them. They were last run against 0.6.4 and are not part of this release's gate.
  • On hold too: the barnase–barstar tutorial, because MD-tools is now dedicated to sampling a ligand's conformations in solvent or in a protein; a protein–protein complex is outside that. Its code is not removed.
  • The paracetamol REST2 and AIS pages and the TYK2 REST2 page record v3 runs and say so in a notice.

What was checked, and what was not

Fast lane (not slow and not gpu, CPU, cores 35–38) on 613b7ed6, which is this release merged with origin/main's symmetry-first analysis work: 3156 passed, 3 skipped, 5 xfailed, 0 failed. The two CUDA-lane tables below were run before that merge, on 5c039d18; the merge touches analysis/, its tests and docs, and no code the CUDA lanes exercise.

gpu or slow lane on CUDA on 5c039d18, whose only later commit is documentation (the barnase–barstar tutorial moved under On hold). The alchemy test files (tests/test_alchemy_*.py, 42 tests of the lane's 638) were excluded by decision, 2026-10-13: alchemical transformation is on hold and is not part of this release's gate. They were not run and nothing here says they pass; three of them (CUDA tolerance tests in test_alchemy_hamiltonian_cuda.py) were already failing on main before this release.

run cards result
main lane cards 1–4 leased, but the lane saw two 569 passed, 11 skipped; its multi-card tests were refused for want of cards and are covered by the two rows below
multi-card files, no MPS 1–4 63 passed, 1 failed (six ranks on four cards needs MPS), 3 skipped
multi-card files, private MPS daemon 1–4 64 passed, 0 failed, 3 skipped
MPS placement tests, evidence file written 1–4 30 passed, 1 failed

The one failure in the last row is test_twelve_workers_on_four_cards_run_together_under_mps. Its placement assertions pass (12 workers, MPS verified, three to a card, a trajectory per state); its last two assume all 48 CPUs of the machine are available (four per worker, 48 distinct), which a shared host cannot give. It is a defect of the test, not of the code, and is to be made to scale to the CPUs it is allowed.

Not run, and why. The reasons for the main lane's 11 skips were not captured (the report flag that lists them was overridden), so they are not claimed. The bespokefit (Psi4) and mace-off23-bespokefit refinement methods are accepted by the same verifier but no fit was run for them; only a MACE-OMOL fit was. build-top does not run any fit.