0.6.2 — nonpolar solvation, and three defects found while shipping it¶

This release CHANGES THE IMPLICIT-SOLVENT HAMILTONIAN, and the previous one was not what it claimed

Two separate things are true, and the second is the reason the first is not a preference.

The nonpolar term now defaults ON. implicit_solvent.nonpolar_sasa was false in 0.6.1 and is true in 0.6.2, matching OpenMM's implicit/gbn2.xml. That is worth ~16 kJ/mol (~6 kT) on ACE-ALA-NME, so an implicit run at 0.6.2 is not the same experiment as the same configuration at 0.6.1.

The setting did not work at all before this release. build_implicit_bundle_inputs accepted nonpolar_sasa and never forwarded it to build_implicit_system. The argument was declared, documented, resolved from configuration, recorded in forcefield.json — and dropped. Every implicit System md-tools has ever built has no nonpolar term in it, whatever its manifest says, and a configuration that asked for nonpolar_sasa: true got exactly the same System as one that asked for false.

Traced with git log -S: the call site has never carried the keyword, going back to the commit that introduced the parameter (1abeaa36). There is no version in which it worked.

So a manifest saying nonpolar_sasa: true from 0.6.1 or earlier describes a request, not a Hamiltonian. One saying false is accurate, because false is what was built either way. Records with the key absent are read as false, which is correct for the same reason.

A source-text test was watching this and could not see it. It asserted "nonpolar_sasa=bool(" in inspect.getsource(builders) — a string in the caller, which was there and correct the whole time, while the callee dropped the value. It has been replaced by an AST check on the call node itself, and that check fails when the keyword is removed.

What this means for existing data. Every registered implicit dataset was built without the nonpolar term. They remain internally consistent and their reference bundles still rebuild, because a bundle carries the recorded value and the standalone rebuild honours it. What must not happen is comparing an implicit result from 0.6.1 with one from 0.6.2 without checking implicit_solvent.nonpolar in both. Explicit-solvent runs are untouched.

The nonpolar term is now a parameter set, not a flag¶

implicit_solvent gained one key:

implicit_solvent:
  model: GBn2
  radii: mbondi3
  nonpolar_sasa: true         # was false in 0.6.1, and inert before that
  nonpolar_surften: null      # kcal/mol/A^2; null means ACE's own 0.0054

nonpolar_surften is the surface tension σ in E_np = σ × SASA. OpenMM offers no way to set it: customgbforces.py writes it into the ACE energy expression as the literal 28.3919551, which is 4π × 0.0054 kcal/mol/Å² converted to kJ/mol/nm². sasaMethod="ACE" and ParmEd's useSASA=True are booleans — on, off, and nothing between. md-tools locates that energy term in the built CustomGBForce and substitutes a prefactor computed from the requested tension, so surften is a real parameter here.

Measured on paracetamol (TYL, 20 atoms), GB energy against the tension:

surften (kcal/mol/Å²)	prefactor	GB energy (kJ/mol)	vs ACE default
0.0054 — ACE, OpenMM's own	28.3919551	−58.96	reproduces exactly
0.005 — Amber `pmemd` GB default	26.2888473	−60.18	−1.21
0.010	52.5776947	−45.01	+13.95
term off (`nonpolar_sasa: false`)	—	−75.34	−16.38

null costs nothing. Rewriting at 0.0054 reproduces OpenMM's literal to every printed digit, so the serialised System is byte-identical to the un-rewritten one — same digest, same checkpoint fingerprint, and a bundle built before the rewrite still verifies against one built after. A test asserts the bytes.

Zero is refused. nonpolar_surften: 0 would leave an energy term that always evaluates to zero: it costs time every step, and it changes the System's serialisation and therefore the fingerprint, so it would not be the same experiment as a run with the term off. The refusal names nonpolar_sasa: false as the way to say that. A tension given alongside nonpolar_sasa: false is refused too, before anything is loaded — that pairing is how an accepted-and-inert setting gets written into a manifest, which is the defect this release exists to remove.

What is recorded¶

forcefield.json gained implicit_solvent.nonpolar, the method by name with the tension and its units:

"nonpolar": {"method": "ACE", "applied": true, "surften": 0.0054,
             "surften_units": "kcal/mol/A^2", "prefactor_kj_mol_nm2": 28.391955111258543,
             "openmm_prefactor_kj_mol_nm2": 28.3919551, "energy_term_index": 1,
             "is_openmm_default": true}

nonpolar_sasa stays as a bool for readers that already parse it, but it stopped identifying the Hamiltonian the moment the tension became adjustable: 0.0054 and 0.005 differ by more than a kJ/mol on a 20-atom ligand, which is the size of difference that silently invalidates a comparison between two runs somebody assumed were the same.

ParmEd and OpenMM disagree about one GB fallback, in γ only¶

An element outside the GB-Neck2 fit gets a placeholder (α, β, γ). ParmEd writes (1.0, 0.8, 4.85) and OpenMM writes (1.0, 0.8, 4.851) — a difference of 0.001 in γ, and nothing else.

It reaches only atoms the fit never covered — the halogens, P, Se, B, Si — so no protein and no C/H/N/O/S ligand can see it. It is recorded rather than reconciled: neither value is derived from anything, both are a library's choice of placeholder, and picking one would build a System that is neither library's. A solute containing such an atom is already flagged by gb_parameter_coverage as carrying parameters from outside the GB-Neck2 fit, and that warning, not this constant, is the thing to act on.

A test reads both libraries' values and fails if either changes, so the note cannot outlive the state of affairs it describes.

`amber_igb8_parity_claimed` is removed¶

forcefield.json no longer carries amber_igb8_parity_claimed or amber_igb8_parity_basis.

They asserted bit-for-bit equivalence with a program this package has never run. No test, no evidence page and no release note has ever compared an md-tools implicit energy against a pmemd one. The value was computed as status == "supported", which comes from a route label plus an element-coverage check — premises about the ingredients Amber would use, which say nothing about whether the two programs produce the same energy. That is the only thing "parity" means to whoever cites it, and forcefield.json is exactly the file somebody quotes when comparing an md-tools number against a published Amber one.

Nothing branched on either field and no schema carried them, so removing them loses no capability. What the record can stand behind is already beside them and stays: implicit_model, radii, radius_assignment_method, nonpolar and parameter_coverage — all measured on the built System.

The claim was nearly entrenched instead: this release's first draft made it conditional on the nonpolar term and gave it a gating test, which would have dressed an unmeasured assertion in machinery that looks like verification. A gap invites work; a confident false entry closes the question.

The comparison has since been made, and the removed claim was wrong as worded. See 20260928-amber-igb8-cross-engine-evidence.md: one prmtop through CPU pmemd 26.0 and through build_implicit_system, nonpolar_sasa=False so gbsa=0 means the same thing on both sides.

Bonded terms agree exactly. Totals agree to 2.6e-4 relative at worst — good agreement, but not "bit-for-bit". And the disagreement is ten times larger under mbondi3 than under mbondi2 on the same molecule at the same coordinates with the same charges (3.1e-4 vs 3.2e-5 on EGB), so the old field claimed parity most confidently exactly where the agreement is weakest, mbondi3 being the default. Ruled out as causes: rgbmax, system size, formal charge, and any per-atom parameter mismatch — 0 of 124 atoms disagree in radius once GBn2's 0.195141 Å offset is applied, and the screen factors match the GBn2 table element for element.

So the field is not merely unmeasured now; it is measured and false. What md-tools may say about Amber is a tolerance with a table behind it, not the word "parity".

The second force classifier is gone¶

md_tools.openmm.system carried its own SCALED_FORCE_CLASSES, DELIBERATELY_UNSCALED_FORCE_CLASSES, ENERGY_FREE_FORCE_CLASSES, UnclassifiedForceError, audit_force_classes and build_rest2_scaled_system — a whole-solute REST2 implementation duplicating md_tools.rest2.hamiltonian. Nothing outside that file imported any of it; every caller reaches the live ones through md_tools.rest2.scaler. 232 lines removed, and the fast lane count did not change, which is what says they were covered by nothing.

Two classifiers is the defect, not the duplication. A force class added to one and not the other makes audit_force_classes refuse in one code path and scale in the other, and which one runs is decided by whichever module the caller happened to import.

Implicit AIS could not run at all¶

Under implicit solvent, AIS refused every run: "the two end states are not a parameter-only pair". Bisected to 090be0e8 (the state renumbering, unreleased). At tau = 0 the scaler returns an untouched clone with no rest2_scale_gb, while the scaled state has that global parameter and its energy terms rewritten — a difference in force LAYOUT, not parameters, so the pair check correctly refused it. Explicit solvent has no CustomGBForce, which is why every AIS tutorial passed and this hid in the implicit path alone.

Fixed by writing the GB factor as a literal in the energy expression for AIS states (gb_literal), rather than behind a global parameter:

tau = 0 stays a byte-identical untouched clone — the documented invariant is preserved, not amended;
REST2 and cMD keep the global parameter, so every saved state's bytes are unchanged;
a DEGENERATE pair — a System against its own unscaled state — still has nothing to switch and is still refused;
the real AIS pair is accepted, with no global parameter to collide.

The first attempt at this fix gave the tau = 0 state an inert rest2_scale_gb = 1.0 so the layouts matched. It made the failing tests pass, and it also made a System and its own unscaled state differ by one inert parameter — silently defeating the "nothing to switch" guard, so eight tests then passed on a null transformation. It is recorded here because the wrong fix looked exactly like the right one.

Never released, so no existing data is affected.

A bond another rule has already classified is not undecided¶

build-top --rest2-scaler refused builds with "1 item(s) could not be classified … RDKit found no ordinary-amide match for this C-N bond", over a bond it had already protected.

_amide_candidates finds every C-N bond whose carbon also bears an oxygen. A C-N bond inside an aromatic ring meets that description and can never match [CX3](=[OX1])[NX3], because SMARTS C is an ALIPHATIC carbon and RDKit aromatises the ring. Folate's pterin 4-oxo lactam is the standard case: the bond was in central_bonds as aromatic_ring, its torsions were never going to be scaled — and the refusal was unanswerable, because the md-tools-torsion-exclusions/1 file only ADDS protection, which the bond already had.

This does not weaken "a torsion nobody can classify is refused". It says a bond the aromatic or double-bond rule placed HAS been classified. A candidate no rule placed is still refused, and a whole residue with no bond orders is never dropped. The resolution is recorded in amide_candidates_resolved_by_another_rule, with the rule that placed it and the amide test's own words for why it could not.

Do not 'fix' this by widening the SMARTS

Three candidate widenings were tried and measured; all three are worse:

[#6X3](=[OX1])[#7X3] — matches the aromatic lactam, which then hits the proline-like branch (the nitrogen is in a 6-ring, inside the bound of 7) and the bond is SCALED. An aromatic ring bond that is protected today would start having its torsions scaled, silently.
[CD3](=[OX1])-!@[#7] — -!@ excludes ring bonds, and in a cyclic peptide every backbone amide IS a ring bond. Measured: 4 amide bonds found, 0 after the change.
Neither closes the case they were aimed at — a deprotonated amide (nitrogen with two connections) fails on X3, the connection count, which changing element classes does not touch.

test_relaxing_the_amide_carbon_to_any_element_would_scale_an_aromatic_ring_bond pins the first hazard. The amidate gap remains open as of 0.6.2, and is 0.6.3 work with cyclic peptides in the test set from the start.

Closed after 0.6.2 — and not by any of the three above

The widening taken was [CX3](=[OX1])[NX2,NX3]: the NITROGEN's connection count, and nothing else. CX3 stays aliphatic, which is what keeps the aromatic lactam out, so the first hazard is structurally unreachable rather than merely avoided; no !@ appears, so a cyclic peptide's backbone amides are still ring bonds and still matched. Measured on both pinned hazards — cyclo(Gly)4 4 → 4, folate's bond 26-28 unmatched — and every currently accepted molecule classifies identically. Only molecules that previously REFUSED change, so no existing run's Hamiltonian moves.

It does route a deprotonated small-ring lactam into the proline-like branch. That is the third hazard's shape but not its defect: the NEUTRAL lactam is already proline-like today, so the anion is being made to agree with its own neutral form, and the <= 7 bound still holds a macrolactam as an ordinary amide.

The pattern is now ORDINARY_AMIDE_SMARTS, defined once in openmm/peptide_map.py and imported by openmm/system.py — two literals answering one question is how they drift apart. detector_version was deliberately NOT bumped, so a record is told apart by the pattern its detection_method names. (After 0.6.2 the stamp was removed altogether: it named the algorithm rather than the decision, had been a constant since 0.6.0, and sat inside the selection digest — see the v4 identity format.)

The refusal names the atoms¶

That message printed the literal letters C and N — the ELEMENT of each end, not the atoms:

bond 29-31: FOL0 C -> FOL0 N: RDKit found no ordinary-amide match ...

Every amide candidate in a residue rendered identically, and the letters read exactly like atom names — actively misleading on folate, whose PDB component really does have atoms named C and N (its benzoylglutamate amide), so the line pointed at a different bond from the one that failed. The bond INDEX was correct throughout; only the label lied. It now names each atom:

bond 6-8: FOL0 C7 -> FOL0 N1: RDKit found no ordinary-amide match ...

Four CUDA-inventory sites were unclassified¶

The inventory guard is gpu-marked, so no CPU lane runs it, and it had four unclassified sites — it would have failed at the first gate. Two were this release's; two had landed earlier and the guard had not run since. Neither tagged release contains either function, so nothing shipped with an unclassified site. All four are matched as derive by System/Topology accessors rather than State accessors and create no Context; each entry says which accessor matched it and why that is not a device read.

Also in this release¶

The ParmEd-free construction route was measured and reproduces the package bit-identically for a parameterised ligand: charges, σ and ε to max |Δ| = 0.00000000, and 20 bonds, 31 angles and 74 torsions identical term for term. No route change is made here — the measurement records that ParmEd is not load-bearing for that path, and the γ divergence above is its only difference.

Both ParmEd-free routes were executed, not just written down:

route	result
protein: `ForceField("amber14-all.xml", "implicit/gbn2.xml")`	22 particles, 12 constraints, ACE term at index 1, −103.7465 kJ/mol
unparameterised ligand: `SMIRNOFFTemplateGenerator` + the same pair	20 particles, ACE present, −200.2069 kJ/mol

The protein figure reproduces the value already recorded in implicit.py's docstring for the ForceField route, to every printed digit. For the ligand, OpenMM assigns GB radii itself (0.1005–0.1505 nm, none zero), so changeRadii("mbondi3") is not what stands between that route and unassigned radii. That remains true only of the ParmEd route, where an OpenFF/Sage prmtop genuinely carries none — which is what the docstring says, and it still says it.

Umbrella sampling moved to 0.6.3.

Compatibility¶


explicit solvent	unchanged in every respect
implicit solvent, 0.6.2 defaults	new Hamiltonian — nonpolar term present for the first time
implicit solvent, `nonpolar_sasa: false`	byte-identical to 0.6.1
`nonpolar_surften: null` with the term on	byte-identical to OpenMM's own ACE System
reference bundles built before 0.6.2	rebuild unchanged; the recorded value is honoured
manifests from 0.6.1 or earlier saying `nonpolar_sasa: true`	describe a request, not what was built
implicit AIS	worked in no unreleased state; fixed here. Saved REST2/cMD state bytes unchanged
REST2 and cMD saved states	byte-identical to 0.6.1 — `gb_literal` is AIS-only
a build that 0.6.1 refused with "could not be classified"	may now succeed; what is SCALED is unchanged

On that last row. The classifier fix changes which builds are allowed, never which torsions are scaled. A molecule that 0.6.1 refused and 0.6.2 builds gets the same Hamiltonian 0.6.1 would have produced had it not refused — the bond in question was already protected. Nothing built on 0.6.1 is invalidated by it.