From 9238a0b6b87106fd8534a2a9529148536e2d9cd9 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 27 Aug 2026 19:32:36 +0100 Subject: [PATCH 1/5] Fix issues with multi-molecule REST2 selections. --- doc/source/changelog.rst | 10 ++++ src/sire/mol/_dynamics.py | 42 +++++++-------- tests/convert/test_openmm_rest2.py | 63 ++++++++++++++++++++++ wrapper/Convert/SireOpenMM/_sommcontext.py | 8 ++- 4 files changed, 99 insertions(+), 24 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 9fac4d157..28499b2aa 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -62,6 +62,16 @@ organisation on `GitHub `__. but can be set to ``False`` to fall back on the internal bond inference heuristic, which is much faster for large molecules, e.g. proteins. +* Fixed a bug where a ``rest2_selection`` spanning more than one molecule raised + ``list.remove(x): x not in list``, since ``selection_to_atoms`` returns a + ``SelectorM`` whose ``to_list()`` gives one view per molecule, rather than a flat + list of atoms. + +* Fixed the atom index offset used when preparing the REST2 data structures, which + was applied to every atom in the selection rather than only those belonging to the + molecule being processed. This gave incorrect indices when a ``rest2_selection`` + spanned more than one non-perturbable molecule. + `2026.1.0 `__ - June 2026 ----------------------------------------------------------------------------------------- diff --git a/src/sire/mol/_dynamics.py b/src/sire/mol/_dynamics.py index 13ea7e2fd..c34ae4492 100644 --- a/src/sire/mol/_dynamics.py +++ b/src/sire/mol/_dynamics.py @@ -248,40 +248,38 @@ def __init__(self, mols=None, map=None, **kwargs): ) # Store all the perturbable molecules associated with the selection - # and remove perturbable atoms from the selection. Remove alchemical ions - # from the selection. + # and exclude perturbable atoms and alchemical ions from the selection. pert_mols = {} - non_pert_atoms = atoms.to_list() + non_pert_atoms = [] for atom in atoms: mol = atom.molecule() if mol.has_property("is_alchemical_ion"): - non_pert_atoms.remove(atom) + continue elif mol.has_property("is_perturbable"): - non_pert_atoms.remove(atom) if mol.number() not in pert_mols: pert_mols[mol.number()] = [atom] else: pert_mols[mol.number()].append(atom) + else: + non_pert_atoms.append(atom) # Now create a boolean is_rest2 mask for the atoms in the perturbable molecules. - # Only do this if there are perturbable atoms in the selection. - if len(non_pert_atoms) != len(atoms): - for num in pert_mols: - mol = self._sire_mols[num] - is_rest2 = [False] * mol.num_atoms() - for atom in pert_mols[num]: - is_rest2[atom.index().value()] = True - - # Set the is_rest2 property for each perturbable molecule. - mol = ( - mol.edit() - .set_property("is_rest2", is_rest2) - .molecule() - .commit() - ) + for num in pert_mols: + mol = self._sire_mols[num] + is_rest2 = [False] * mol.num_atoms() + for atom in pert_mols[num]: + is_rest2[atom.index().value()] = True - # Update the system. - self._sire_mols.update(mol) + # Set the is_rest2 property for each perturbable molecule. + mol = ( + mol.edit() + .set_property("is_rest2", is_rest2) + .molecule() + .commit() + ) + + # Update the system. + self._sire_mols.update(mol) # Search for alchemical ions and exclude them via a REST2 mask. try: diff --git a/tests/convert/test_openmm_rest2.py b/tests/convert/test_openmm_rest2.py index 7ac498a52..4c81745c8 100644 --- a/tests/convert/test_openmm_rest2.py +++ b/tests/convert/test_openmm_rest2.py @@ -13,6 +13,69 @@ def toluene_methane(): return sr.load_test_files("toluene_methane.s3") +@pytest.mark.parametrize("mols", ["ala_mols", "merged_ethane_methanol"]) +def test_rest2_selection_multiple_molecules(mols, request): + """ + Test that a REST2 selection spanning multiple molecules is applied to the + atoms in each of the selected molecules. + """ + + mols = request.getfixturevalue(mols) + + # Link to the reference state. + try: + mols = sr.morph.link_to_reference(mols) + except: + pass + + # The REST2 region is the union of the first two molecules. Work out the + # system indices of their atoms. Perturbable molecules are scaled via the + # lambda lever rather than the NonbondedForce, so are excluded here. + scaled_atoms = set() + num_selected_atoms = 0 + for mol in [mols[0], mols[1]]: + if not mol.has_property("is_perturbable"): + scaled_atoms.update( + range(num_selected_atoms, num_selected_atoms + mol.num_atoms()) + ) + num_selected_atoms += mol.num_atoms() + + # Create a dynamics object, selecting the first two molecules. + d = mols.dynamics(platform="Reference", rest2_selection="molidx 0 or molidx 1") + + # Find the NonbondedForce. + for force in d.context().getSystem().getForces(): + if force.getName() == "NonbondedForce": + break + + # Store the initial parameters. + nonbonded_params_initial = [ + force.getParticleParameters(i) for i in range(force.getNumParticles()) + ] + + # Update the REST2 scaling factor. + d.set_lambda(0.0, rest2_scale=2.0) + + # Find the NonbondedForce. + for force in d.context().getSystem().getForces(): + if force.getName() == "NonbondedForce": + break + + # Store the scaling factor. + scale = 0.5 + + # Only the atoms in the two selected molecules should be scaled. + for i in range(force.getNumParticles()): + charge, _, epsilon = nonbonded_params_initial[i] + charge_modified, _, epsilon_modified = force.getParticleParameters(i) + if i in scaled_atoms: + assert isclose(charge_modified._value, charge._value * scale**0.5) + assert isclose(epsilon_modified._value, epsilon._value * scale) + elif i >= num_selected_atoms: + assert isclose(charge_modified._value, charge._value) + assert isclose(epsilon_modified._value, epsilon._value) + + @pytest.mark.parametrize( ["mols", "rest2_selection", "excluded_atoms"], [ diff --git a/wrapper/Convert/SireOpenMM/_sommcontext.py b/wrapper/Convert/SireOpenMM/_sommcontext.py index b5d8233cf..f2c8e2e78 100644 --- a/wrapper/Convert/SireOpenMM/_sommcontext.py +++ b/wrapper/Convert/SireOpenMM/_sommcontext.py @@ -510,8 +510,12 @@ def _prepare_rest2(self, system, atoms): for i in range(mol_idx): num_atoms += system_mols[i].num_atoms() - # Create a list of atom indices. - atom_idxs = [atom.index().value() + num_atoms for atom in atoms] + # Create a list of system indices for the selected atoms in this molecule. + atom_idxs = [ + atom.index().value() + num_atoms + for atom in atoms + if atom.molecule().number() == mol.number() + ] # Gather the nonbonded parameters for the atoms in the selection. for idx in atom_idxs: From 2bd50fe8b1f46f0c8db00b36f033674c4cdc5bc4 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 27 Aug 2026 19:38:35 +0100 Subject: [PATCH 2/5] Remove kwarg from generic .to() method. --- src/sire/convert/__init__.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/src/sire/convert/__init__.py b/src/sire/convert/__init__.py index e036d9d05..7cf31d335 100644 --- a/src/sire/convert/__init__.py +++ b/src/sire/convert/__init__.py @@ -48,7 +48,7 @@ def supported_formats(): return _supported_formats() -def to(obj, format: str = "sire", map=None, determine_bond_orders: bool = True): +def to(obj, format: str = "sire", map=None): """ Convert the passed object from its current object format to the specified object format (default "sire"). Typically this will be converting @@ -62,19 +62,13 @@ def to(obj, format: str = "sire", map=None, determine_bond_orders: bool = True): The format to convert to map: The property map to use for the conversion - determine_bond_orders: bool (default True) - Whether to use RDKit's ``determineBondOrders`` function when bond - orders need to be inferred during conversion to rdkit format. This - is more robust than the internal heuristic, but can be slow for - large molecules, e.g. proteins. (Only used when converting to - rdkit format.) """ format = format.lower() if format == "sire": return to_sire(obj, map=map) elif format == "rdkit": - return to_rdkit(obj, map=map, determine_bond_orders=determine_bond_orders) + return to_rdkit(obj, map=map) elif format == "gemmi": return to_gemmi(obj, map=map) elif format == "biosimspace": From 0b8dd333afffacac000d922e7e7a13b67dd490a3 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 27 Aug 2026 19:40:29 +0100 Subject: [PATCH 3/5] Use a set for REST2 atom indices to avoid quadratic membership tests. --- doc/source/changelog.rst | 4 ++++ wrapper/Convert/SireOpenMM/_sommcontext.py | 8 ++++---- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 28499b2aa..f2d0be394 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -72,6 +72,10 @@ organisation on `GitHub `__. molecule being processed. This gave incorrect indices when a ``rest2_selection`` spanned more than one non-perturbable molecule. +* Used a set rather than a list for the atom indices when preparing the REST2 data + structures. The indices are membership tested against every exception and torsion + in the system, which was quadratic for large REST2 regions, e.g. proteins. + `2026.1.0 `__ - June 2026 ----------------------------------------------------------------------------------------- diff --git a/wrapper/Convert/SireOpenMM/_sommcontext.py b/wrapper/Convert/SireOpenMM/_sommcontext.py index f2c8e2e78..54289232c 100644 --- a/wrapper/Convert/SireOpenMM/_sommcontext.py +++ b/wrapper/Convert/SireOpenMM/_sommcontext.py @@ -510,15 +510,15 @@ def _prepare_rest2(self, system, atoms): for i in range(mol_idx): num_atoms += system_mols[i].num_atoms() - # Create a list of system indices for the selected atoms in this molecule. - atom_idxs = [ + # Create a set of system indices for the selected atoms in this molecule. + atom_idxs = { atom.index().value() + num_atoms for atom in atoms if atom.molecule().number() == mol.number() - ] + } # Gather the nonbonded parameters for the atoms in the selection. - for idx in atom_idxs: + for idx in sorted(atom_idxs): self._nonbonded_params[idx] = nonbonded_force.getParticleParameters(idx) # Store the exception parameters. From cadbadc4b156b88260b8bbf5a5e60d253a947a1f Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 27 Aug 2026 19:45:26 +0100 Subject: [PATCH 4/5] Add tests for additive and narrowing REST2 selection semantics. --- tests/convert/test_openmm_rest2.py | 104 +++++++++++++++++++++++------ 1 file changed, 82 insertions(+), 22 deletions(-) diff --git a/tests/convert/test_openmm_rest2.py b/tests/convert/test_openmm_rest2.py index 4c81745c8..67d767a38 100644 --- a/tests/convert/test_openmm_rest2.py +++ b/tests/convert/test_openmm_rest2.py @@ -13,32 +13,16 @@ def toluene_methane(): return sr.load_test_files("toluene_methane.s3") -@pytest.mark.parametrize("mols", ["ala_mols", "merged_ethane_methanol"]) -def test_rest2_selection_multiple_molecules(mols, request): +def test_rest2_selection_multiple_molecules(ala_mols): """ Test that a REST2 selection spanning multiple molecules is applied to the atoms in each of the selected molecules. """ - mols = request.getfixturevalue(mols) + mols = ala_mols - # Link to the reference state. - try: - mols = sr.morph.link_to_reference(mols) - except: - pass - - # The REST2 region is the union of the first two molecules. Work out the - # system indices of their atoms. Perturbable molecules are scaled via the - # lambda lever rather than the NonbondedForce, so are excluded here. - scaled_atoms = set() - num_selected_atoms = 0 - for mol in [mols[0], mols[1]]: - if not mol.has_property("is_perturbable"): - scaled_atoms.update( - range(num_selected_atoms, num_selected_atoms + mol.num_atoms()) - ) - num_selected_atoms += mol.num_atoms() + # The REST2 region is the union of the first two molecules. + num_rest2_atoms = mols[0].num_atoms() + mols[1].num_atoms() # Create a dynamics object, selecting the first two molecules. d = mols.dynamics(platform="Reference", rest2_selection="molidx 0 or molidx 1") @@ -68,10 +52,86 @@ def test_rest2_selection_multiple_molecules(mols, request): for i in range(force.getNumParticles()): charge, _, epsilon = nonbonded_params_initial[i] charge_modified, _, epsilon_modified = force.getParticleParameters(i) - if i in scaled_atoms: + if i < num_rest2_atoms: assert isclose(charge_modified._value, charge._value * scale**0.5) assert isclose(epsilon_modified._value, epsilon._value * scale) - elif i >= num_selected_atoms: + else: + assert isclose(charge_modified._value, charge._value) + assert isclose(epsilon_modified._value, epsilon._value) + + +@pytest.mark.parametrize( + ["rest2_selection", "pert_atoms", "extra_mols"], + [ + # No selection, so the region is the entire perturbable molecule. + (None, None, []), + # A whole non-perturbable molecule, which is added to the entire + # perturbable molecule. + ("molidx 1", None, [1]), + # Part of the perturbable molecule, which narrows the region to those + # atoms alone. + ("molidx 0 and atomidx 0,1", [0, 1], []), + # Part of the perturbable molecule plus a whole non-perturbable + # molecule, which are combined. + ("(molidx 0 and atomidx 0,1) or molidx 1", [0, 1], [1]), + ], +) +def test_rest2_selection_semantics( + merged_ethane_methanol, rest2_selection, pert_atoms, extra_mols +): + """ + Test that a REST2 selection adds to the default region of the whole + perturbable molecule, and that selecting atoms within the perturbable + molecule narrows the region to those atoms. + """ + + mols = sr.morph.link_to_reference(merged_ethane_methanol) + + # Work out the system index of the first atom of each molecule. + offsets = [] + offset = 0 + for mol in mols: + offsets.append(offset) + offset += mol.num_atoms() + + # Work out the system indices of the atoms in the REST2 region. The + # perturbable molecule is molecule zero. + if pert_atoms is None: + pert_atoms = range(mols[0].num_atoms()) + rest2_atoms = {offsets[0] + i for i in pert_atoms} + for i in extra_mols: + rest2_atoms.update(range(offsets[i], offsets[i] + mols[i].num_atoms())) + + # Create a dynamics object. + d = mols.dynamics(platform="Reference", rest2_selection=rest2_selection) + + # Find the NonbondedForce. + for force in d.context().getSystem().getForces(): + if force.getName() == "NonbondedForce": + break + + # Store the unscaled parameters at the same lambda value, so that the + # comparison isolates the REST2 scaling from the lambda lever. + d.set_lambda(0.0, rest2_scale=1.0) + nonbonded_params_initial = [ + force.getParticleParameters(i) for i in range(force.getNumParticles()) + ] + + # Update the REST2 scaling factor. + d.set_lambda(0.0, rest2_scale=2.0) + + # Store the scaling factor. + scale = 0.5 + + # Only the atoms in the REST2 region should be scaled. + for i in range(force.getNumParticles()): + charge, _, epsilon = nonbonded_params_initial[i] + charge_modified, _, epsilon_modified = force.getParticleParameters(i) + if i in rest2_atoms: + assert isclose(charge_modified._value, charge._value * scale**0.5) + if epsilon._value > 1e-6: + assert isclose(epsilon_modified._value, epsilon._value * scale) + else: assert isclose(charge_modified._value, charge._value) assert isclose(epsilon_modified._value, epsilon._value) From 6c09d79b1d3d3f453170f3c4758400b83b92e457 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 27 Aug 2026 20:15:20 +0100 Subject: [PATCH 5/5] Add missing includes for std::getenv and std::exit. --- corelib/src/apps/sire/main.cpp | 2 ++ corelib/src/libs/SireBase/tempdir.cpp | 2 ++ corelib/src/libs/SireCluster/mpi/mpicluster.cpp | 2 ++ corelib/src/libs/SireIO/trajectorymonitor.cpp | 2 ++ corelib/src/libs/SireMove/simstore.cpp | 2 ++ 5 files changed, 10 insertions(+) diff --git a/corelib/src/apps/sire/main.cpp b/corelib/src/apps/sire/main.cpp index bb206a109..ce084d5a2 100644 --- a/corelib/src/apps/sire/main.cpp +++ b/corelib/src/apps/sire/main.cpp @@ -33,6 +33,8 @@ #include "sire_version.h" +#include + using namespace SireCluster; using namespace SireMove; using namespace SireSystem; diff --git a/corelib/src/libs/SireBase/tempdir.cpp b/corelib/src/libs/SireBase/tempdir.cpp index d3269c406..58f1597ec 100644 --- a/corelib/src/libs/SireBase/tempdir.cpp +++ b/corelib/src/libs/SireBase/tempdir.cpp @@ -36,6 +36,8 @@ #include +#include + using namespace SireBase; static QString getUserName() diff --git a/corelib/src/libs/SireCluster/mpi/mpicluster.cpp b/corelib/src/libs/SireCluster/mpi/mpicluster.cpp index 8ef624596..a302d8849 100644 --- a/corelib/src/libs/SireCluster/mpi/mpicluster.cpp +++ b/corelib/src/libs/SireCluster/mpi/mpicluster.cpp @@ -53,6 +53,8 @@ #include +#include + using namespace SireCluster; using namespace SireCluster::MPI; diff --git a/corelib/src/libs/SireIO/trajectorymonitor.cpp b/corelib/src/libs/SireIO/trajectorymonitor.cpp index 70ff1b906..7f67009b8 100644 --- a/corelib/src/libs/SireIO/trajectorymonitor.cpp +++ b/corelib/src/libs/SireIO/trajectorymonitor.cpp @@ -45,6 +45,8 @@ #include +#include + using std::shared_ptr; using namespace SireIO; diff --git a/corelib/src/libs/SireMove/simstore.cpp b/corelib/src/libs/SireMove/simstore.cpp index 56d76cae0..0ccff7949 100644 --- a/corelib/src/libs/SireMove/simstore.cpp +++ b/corelib/src/libs/SireMove/simstore.cpp @@ -38,6 +38,8 @@ #include +#include + using namespace SireMove; using namespace SireSystem; using namespace SireStream;