[lohedges]

Closes #519.

Reverts the rule added in feba446, which constrained unperturbed heavy-atom bonds at the alchemical junction to match SOMD1. For molecules with large ghost regions, it constrained almost every bond, and the constraint solver diverged within a few steps. SOMD1 only did this because a mass check is missing in its hbonds-notperturbed code. The SOMD1 paper describes constraining only bonds to hydrogen that aren't perturbed.

• Removed the junction rule from the bond constraint condition.
• Inverted test_cyclopropyl_anchor_bond_constraint to expect the anchor bonds to be unconstrained.
• Added test_macrocycle_ghost_constraints (marked slow), which runs dynamics at λ = 0 and 1 for a partially mapped macrocycle that blew up before this fix.

Also closes #521 by removing constraints for bonds that are broken during a perturbation.

Closes #522

Interactions between ghost atoms in the same molecule were calculated by GhostGhostNonbondedForce, a CustomNonbondedForce with a periodic cutoff and no switching. Its Coulomb term is q1*q2*(1/sqrt(δα + r²) - κ/r), so when a schedule sets κ below 1 to keep a decoupled molecule's intramolecular electrostatics, the force holds a bare 1/r that is truncated at the cutoff. For molecules larger than the cutoff this made the energy discontinuous: ABFE decoupling of a ~16 Å ligand failed to minimise at most λ windows and blew up in the first cycle of dynamics. These pairs are now excluded from the cutoff forces. Their hard Coulomb is a NonbondedForce exception, and their soft-core terms are in a new GhostGhostBondForce with no cutoff, which is in the ghost/ghost force group and controlled by the ghost/ghost levers. Pairs in different bonded fragments stay in the cutoff force. Energies for other schedules change slightly, since intramolecular ghost LJ beyond the cutoff is no longer truncated (-0.26 kJ/mol for the 58-atom test ligand). For that system, MD throughput is 3-5% higher and set_lambda is about 0.2-0.3 ms slower per call, mostly in OpenMM's parameter uploads.

Closes #523

After a line search failure, L-BFGS reverted its coordinates but the context kept the rejected trial point, so dynamics could start from it when minimisation errors were ignored. The minimiser now re-evaluates the accepted point after each L-BFGS call. On failure it keeps the last accepted positions, unless their energy is no lower than the starting energy. Timeouts are reported as such, and line search retries are capped rather than repeating until the timeout.

Closes #524

For ring-breaking and ring-making perturbations, the pair whose bond is broken, and the 1-3 and 1-4 pairs across that bond, are NonbondedForce exceptions in the merged molecule, since they are excluded or scaled while the bond exists. OpenMM calculates exceptions as plain Coulomb, without reaction field terms, so the end state where the bond doesn't exist was wrong when using reaction field (0.53 kcal/mol in vacuum and 0.38 kcal/mol with a 10 Å periodic cutoff for a chk1 scaffold hop). PME and no cutoff were unaffected. When using reaction field, correction forces in the same force group as the NonbondedForce now add the difference between each exception and a normal reaction field pair. For the broken pair this follows the exception's coul_kappa scaled charge. For the other pairs it applies to exceptions that are normal pairs, with both the Coulomb and LJ scale equal to 1, in either end state, and is interpolated with the charge_scale lever. The end states now match the plain end-state molecules for all cutoff types, for ring-breaking and ring-making.

Closes #525

The LJ interaction of the pair whose bond is broken uses the epsilon from the end state where it isn't bonded, with a soft-core alpha that is 1 while the bond exists. The Beutler form scales this term by (1 - alpha), but the Zacharias and Taylor forms, which are used for RBFE, only soften it, so the bonded pair still had a softened LJ interaction at the bonded end state (0.058 kcal/mol for the chk1 scaffold hop). The Zacharias and Taylor forms are now also scaled by (1 - alpha), so the interaction is zero while the pair is bonded and the bonded end state matches the plain molecule. This changes the soft-core path at intermediate λ for ring-breaking perturbations using these forms.

Other minor fixes:

Lever values. get_lever_values() reported levers that only changed by floating point round-off, such as bonded terms that are the same in both end states. It now compares values with a tolerance.

Scale kappa in standard_decouple() again. This reverts the change in 2025.1.0. With κ fixed, the decoupled molecule's intramolecular electrostatic interactions are turned off, so a single-leg decoupling calculation, e.g. of a hydration free energy, misses the vacuum contribution (31.9 kcal/mol for alanine dipeptide). Water is unaffected, since all of its intramolecular pairs are excluded. The cutoff problem that made scaled κ unstable for large molecules is fixed above.

• I confirm that I have merged the latest version of devel into this branch before issuing this pull request (e.g. by running git pull origin devel): [y]
• I confirm that I have added a changelog entry to the changelog (we will add a link to this PR as part of the review): [y]
• I confirm that I have permission to release this code under the MIT license (the Python API) or the LGPLv3+ license (everything else), as described in the LICENSE file: [y]