Backport fix from PR #467

corelib/src/apps/sire/main.cpp

@@ -33,6 +33,8 @@
 
 #include "sire_version.h"
 
+#include <cstdlib>
+
 using namespace SireCluster;
 using namespace SireMove;
 using namespace SireSystem;

corelib/src/libs/SireBase/tempdir.cpp

@@ -36,6 +36,8 @@
 
 #include <QDebug>
 
+#include <cstdlib>
+
 using namespace SireBase;
 
 static QString getUserName()

corelib/src/libs/SireCluster/mpi/mpicluster.cpp

@@ -53,6 +53,8 @@
 
 #include <QDebug>
 
+#include <cstdlib>
+
 using namespace SireCluster;
 using namespace SireCluster::MPI;
 

corelib/src/libs/SireIO/trajectorymonitor.cpp

@@ -45,6 +45,8 @@
 
 #include <QDebug>
 
+#include <cstdlib>
+
 using std::shared_ptr;
 
 using namespace SireIO;

corelib/src/libs/SireMove/simstore.cpp

@@ -38,6 +38,8 @@
 
 #include <QDebug>
 
+#include <cstdlib>
+
 using namespace SireMove;
 using namespace SireSystem;
 using namespace SireStream;

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":

src/sire/mol/_dynamics.py

@@ -218,40 +218,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:

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"],
     [

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.