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; diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 9fac4d157..f2d0be394 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -62,6 +62,20 @@ 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. + +* 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/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": 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..67d767a38 100644 --- a/tests/convert/test_openmm_rest2.py +++ b/tests/convert/test_openmm_rest2.py @@ -13,6 +13,129 @@ def toluene_methane(): return sr.load_test_files("toluene_methane.s3") +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 = ala_mols + + # 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") + + # 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 < num_rest2_atoms: + assert isclose(charge_modified._value, charge._value * scale**0.5) + assert isclose(epsilon_modified._value, epsilon._value * scale) + 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) + + @pytest.mark.parametrize( ["mols", "rest2_selection", "excluded_atoms"], [ diff --git a/wrapper/Convert/SireOpenMM/_sommcontext.py b/wrapper/Convert/SireOpenMM/_sommcontext.py index b5d8233cf..54289232c 100644 --- a/wrapper/Convert/SireOpenMM/_sommcontext.py +++ b/wrapper/Convert/SireOpenMM/_sommcontext.py @@ -510,11 +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 atom indices. - atom_idxs = [atom.index().value() + num_atoms for atom in atoms] + # 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.