Fix #208

doc/source/changelog.rst

@@ -20,6 +20,8 @@ Devel
 * Minor internal updates due to Sire API fixes (`#203 <https://github.com/OpenBioSim/biosimspace/pull/203>`__).
 * Fix bug in the BSS Boresch restraint search code (`@fjclark <https://github.com/fjclark>`_) (`#204 <https://github.com/OpenBioSim/biosimspace/pull/204>`__).
 * Added SOMD and GROMACS support for multiple distance restraints for ABFE calculations (`#178 <https://github.com/OpenBioSim/biosimspace/pull/178>`__).
+* Fix ``renumber`` option in :meth:`extract <BioSimSpace._SireWrappers.Molecule.extract>` method (`#210 <https://github.com/OpenBioSim/biosimspace/pull/210>`__).
+* Add workaround for fixing reconstruction of intrascale matrix in :func:`readPerturbableSystem <BioSimSpace.IO.readPerturbableSystem>` function (`#210 <https://github.com/OpenBioSim/biosimspace/pull/210>`__).
 
 `2023.4.0 <https://github.com/openbiosim/biosimspace/compare/2023.3.1...2023.4.0>`_ - Oct 13 2023
 -------------------------------------------------------------------------------------------------

python/BioSimSpace/IO/_io.py

@@ -1100,9 +1100,13 @@ def readPerturbableSystem(top0, coords0, top1, coords1, property_map={}):
     # Flag that the molecule is perturbable.
     mol.setProperty("is_perturbable", _SireBase.wrap(True))
 
+    # Get the two molecules.
+    mol0 = system0[idx]._sire_object
+    mol1 = system1[idx]._sire_object
+
     # Add the molecule0 and molecule1 properties.
-    mol.setProperty("molecule0", system0[idx]._sire_object)
-    mol.setProperty("molecule1", system1[idx]._sire_object)
+    mol.setProperty("molecule0", mol0)
+    mol.setProperty("molecule1", mol1)
 
     # Get the connectivity property name.
     conn_prop = property_map.get("connectivity", "connectivity")
@@ -1121,6 +1125,23 @@ def readPerturbableSystem(top0, coords0, top1, coords1, property_map={}):
         mol = mol.removeProperty(conn_prop + "0").molecule()
         mol = mol.removeProperty(conn_prop + "1").molecule()
 
+    # Reconstruct the intrascale matrices using the GroTop parser.
+    intra0 = (
+        _SireIO.GroTop(_Molecule(mol0).toSystem()._sire_object)
+        .toSystem()[0]
+        .property("intrascale")
+    )
+    intra1 = (
+        _SireIO.GroTop(_Molecule(mol1).toSystem()._sire_object)
+        .toSystem()[0]
+        .property("intrascale")
+    )
+
+    # Set the "intrascale" properties.
+    intrascale_prop = property_map.get("intrascale", "intrascale")
+    mol.setProperty(intrascale_prop + "0", intra0)
+    mol.setProperty(intrascale_prop + "1", intra0)
+
     # Commit the changes.
     mol = _Molecule(mol.commit())
 

python/BioSimSpace/Sandpit/Exscientia/IO/_io.py

@@ -1100,9 +1100,13 @@ def readPerturbableSystem(top0, coords0, top1, coords1, property_map={}):
     # Flag that the molecule is perturbable.
     mol.setProperty("is_perturbable", _SireBase.wrap(True))
 
+    # Get the two molecules.
+    mol0 = system0[idx]._sire_object
+    mol1 = system1[idx]._sire_object
+
     # Add the molecule0 and molecule1 properties.
-    mol.setProperty("molecule0", system0[idx]._sire_object)
-    mol.setProperty("molecule1", system1[idx]._sire_object)
+    mol.setProperty("molecule0", mol0)
+    mol.setProperty("molecule1", mol1)
 
     # Get the connectivity property name.
     conn_prop = property_map.get("connectivity", "connectivity")
@@ -1121,6 +1125,23 @@ def readPerturbableSystem(top0, coords0, top1, coords1, property_map={}):
         mol = mol.removeProperty(conn_prop + "0").molecule()
         mol = mol.removeProperty(conn_prop + "1").molecule()
 
+    # Reconstruct the intrascale matrices using the GroTop parser.
+    intra0 = (
+        _SireIO.GroTop(_Molecule(mol0).toSystem()._sire_object)
+        .toSystem()[0]
+        .property("intrascale")
+    )
+    intra1 = (
+        _SireIO.GroTop(_Molecule(mol1).toSystem()._sire_object)
+        .toSystem()[0]
+        .property("intrascale")
+    )
+
+    # Set the "intrascale" properties.
+    intrascale_prop = property_map.get("intrascale", "intrascale")
+    mol.setProperty(intrascale_prop + "0", intra0)
+    mol.setProperty(intrascale_prop + "1", intra0)
+
     # Commit the changes.
     mol = _Molecule(mol.commit())
 

python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py

@@ -340,23 +340,18 @@ def extract(self, indices, renumber=False, property_map={}):
             # Append to the updated atom index.
             indices_.append(_SireMol.AtomIdx(x))
 
-        if renumber:
-            mol = self.copy()
-        else:
-            mol = self
-
         # Extract a partial molecule.
         try:
             # Create an empty atom selection for this molecule.
-            selection = mol._sire_object.selection()
+            selection = self._sire_object.selection()
             selection.selectNone()
 
             # Add the atom indices to the selection.
             for idx in indices_:
                 selection.select(idx)
 
             partial_mol = (
-                _SireMol.PartialMolecule(mol._sire_object, selection)
+                _SireMol.PartialMolecule(self._sire_object, selection)
                 .extract()
                 .molecule()
             )
@@ -371,7 +366,7 @@ def extract(self, indices, renumber=False, property_map={}):
         intrascale = property_map.get("intrascale", "intrascale")
 
         # Flag whether the molecule has an intrascale property.
-        has_intrascale = mol._sire_object.hasProperty(intrascale)
+        has_intrascale = self._sire_object.hasProperty(intrascale)
 
         # Remove the "intrascale" property, since this doesn't correspond to the
         # extracted molecule.
@@ -409,6 +404,15 @@ def extract(self, indices, renumber=False, property_map={}):
         else:
             mol = Molecule(partial_mol)
 
+        # Keep the same MolNum.
+        if not renumber:
+            mol._sire_object = (
+                mol._sire_object.edit()
+                .renumber(self._sire_object.number())
+                .commit()
+                .molecule()
+            )
+
         return mol
 
     def molecule0(self):

python/BioSimSpace/_SireWrappers/_molecule.py

@@ -340,23 +340,18 @@ def extract(self, indices, renumber=False, property_map={}):
             # Append to the updated atom index.
             indices_.append(_SireMol.AtomIdx(x))
 
-        if renumber:
-            mol = self.copy()
-        else:
-            mol = self
-
         # Extract a partial molecule.
         try:
             # Create an empty atom selection for this molecule.
-            selection = mol._sire_object.selection()
+            selection = self._sire_object.selection()
             selection.selectNone()
 
             # Add the atom indices to the selection.
             for idx in indices_:
                 selection.select(idx)
 
             partial_mol = (
-                _SireMol.PartialMolecule(mol._sire_object, selection)
+                _SireMol.PartialMolecule(self._sire_object, selection)
                 .extract()
                 .molecule()
             )
@@ -371,7 +366,7 @@ def extract(self, indices, renumber=False, property_map={}):
         intrascale = property_map.get("intrascale", "intrascale")
 
         # Flag whether the molecule has an intrascale property.
-        has_intrascale = mol._sire_object.hasProperty(intrascale)
+        has_intrascale = self._sire_object.hasProperty(intrascale)
 
         # Remove the "intrascale" property, since this doesn't correspond to the
         # extracted molecule.
@@ -409,6 +404,15 @@ def extract(self, indices, renumber=False, property_map={}):
         else:
             mol = Molecule(partial_mol)
 
+        # Keep the same MolNum.
+        if not renumber:
+            mol._sire_object = (
+                mol._sire_object.edit()
+                .renumber(self._sire_object.number())
+                .commit()
+                .molecule()
+            )
+
         return mol
 
     def molecule0(self):

tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py

@@ -102,3 +102,27 @@ def test_hydrogen_mass_repartitioning(system, ignore_waters):
 
     # Assert the the masses are approximately the same.
     assert final_mass == pytest.approx(initial_mass)
+
+
+def test_extract(system):
+    """Test the extract method."""
+
+    # A list of atom indices to extract.
+    idxs = [0, 1, 2, 3]
+
+    # Extract the first molecule from the system.
+    mol = system[0]
+
+    # Extract the atoms.
+    partial_mol = mol.extract(idxs)
+
+    assert partial_mol.nAtoms() == 4
+
+    # Make sure the numbers match.
+    assert partial_mol.number() == mol.number()
+
+    # Extract and renumber.
+    partial_mol = mol.extract(idxs, renumber=True)
+
+    # Make sure the numbers are different.
+    assert partial_mol.number() != mol.number()

tests/_SireWrappers/test_molecule.py

@@ -96,3 +96,27 @@ def test_hydrogen_mass_repartitioning(system, ignore_waters):
 
     # Assert the the masses are approximately the same.
     assert final_mass == pytest.approx(initial_mass)
+
+
+def test_extract(system):
+    """Test the extract method."""
+
+    # A list of atom indices to extract.
+    idxs = [0, 1, 2, 3]
+
+    # Extract the first molecule from the system.
+    mol = system[0]
+
+    # Extract the atoms.
+    partial_mol = mol.extract(idxs)
+
+    assert partial_mol.nAtoms() == 4
+
+    # Make sure the numbers match.
+    assert partial_mol.number() == mol.number()
+
+    # Extract and renumber.
+    partial_mol = mol.extract(idxs, renumber=True)
+
+    # Make sure the numbers are different.
+    assert partial_mol.number() != mol.number()