Skip to content

Stale coordinates1 when perturbed_system is set #188

Description

@mark-mackey-cresset

From a Claude review that ended up examining the somd2 code:

When replica_exchange is enabled and perturbed_system is set, somd2 seeds the λ > 0.5 replicas at somd2/runner/_base.py:604-617 by
flattening the perturbed system's λ=1 coordinates with sire.io.get_coords_array(sire.morph.link_to_perturbed(...)) and stamping them
onto a clone of the reference system via sire.legacy.IO.setCoordinates(self._system._system, pert_coords.tolist()). That call omits
the third argument; the signature is setCoordinates(system, coordinates, is_lambda1=False, map={}), so for a perturbable (merged)
molecule the incoming coordinates are written to coordinates0 only, and coordinates1 silently retains whatever the reference system
had. Ordinary molecules have a single coordinates property and update correctly, so the inconsistency is confined to the merged
molecule(s) — but it is not cosmetic whenever the perturbed system is not in the same coordinate frame as the reference system, which
is common, since the two are equilibrated independently and can differ by an arbitrary rotation, translation, or periodic image. In
our case (a 59.9° rotation plus a 51 Å translation between the two inputs) the seeded merged ligand ends up with coordinates0 1.78 Å
from the nearest protein atom (correct) and coordinates1 22.47 Å away, still in the reference system's frame — the two end states of
a single 35-atom molecule sitting ~53 Å apart. Nothing repairs this afterwards: sire's dynamics only ever writes back the
currently-linked coordinate set (running 0.2 ps at λ=0.75 linked to reference moved coordinates0 by 2.28 Å and coordinates1 by
xactly 0.0000 Å; linked to perturbed, exactly the reverse), and the divergence survives a sire.stream round-trip, so checkpoints
carry it forward. The OpenMM starting positions themselves are correct — they come from coordinates0 through link_to_reference — so
we have not identified a path where this alters a trajectory; the risk is that any consumer reading the λ=1 end state of a seeded or
restarted replica gets coordinates in the wrong frame. The most likely candidate we found is _base.py:1430/1443, where checkpoint
systems are passed through link_to_perturbed and then get_coords_array when rebuilding GCMC restart positions. The minimal fix would
be to populate coordinates1 as well when seeding.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions