Merge pull request #281 from OpenBioSim/fix_280

Fix issue #280

python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py

@@ -153,6 +153,14 @@ def getFrame(trajectory, topology, index, system=None, property_map={}):
         # Update the water topology to match topology/trajectory.
         system = _update_water_topology(system, topology, trajectory, property_map)
 
+        # Copy the system.
+        renumbered_system = system.copy()
+
+        # Make sure the constituents of the system are numbered in ascending order.
+        renumbered_system._sire_object = _SireIO.renumberConstituents(
+            system._sire_object
+        )
+
     # Try to load the frame with Sire.
     errors = []
     is_sire = False
@@ -219,9 +227,7 @@ def getFrame(trajectory, topology, index, system=None, property_map={}):
                 if "space" in new_system.propertyKeys():
                     box = new_system.property("space")
                     if box.isPeriodic():
-                        sire_system.setProperty(
-                            self._property_map.get("space", "space"), box
-                        )
+                        sire_system.setProperty(property_map.get("space", "space"), box)
 
                 new_system = _System(sire_system)
 
@@ -265,7 +271,7 @@ def getFrame(trajectory, topology, index, system=None, property_map={}):
             # coordinates of all of the atoms in the reference. As such, we
             # will need to split the system into molecules.
             new_system = _split_molecules(
-                frame, pdb, system, str(work_dir), property_map
+                frame, pdb, renumbered_system, str(work_dir), property_map
             )
             try:
                 sire_system, _ = _SireIO.updateCoordinatesAndVelocities(
@@ -276,9 +282,7 @@ def getFrame(trajectory, topology, index, system=None, property_map={}):
                 if "space" in new_system.propertyKeys():
                     box = new_system.property("space")
                     if box.isPeriodic():
-                        sire_system.setProperty(
-                            self._property_map.get("space", "space"), box
-                        )
+                        sire_system.setProperty(property_map.get("space", "space"), box)
 
                 new_system = _System(sire_system)
             except Exception as e:
@@ -288,11 +292,6 @@ def getFrame(trajectory, topology, index, system=None, property_map={}):
                 else:
                     raise IOError(msg) from None
 
-            else:
-                raise IOError(
-                    "The trajectory frame is incompatible with the passed system!"
-                )
-
     # Load the frame directly to create a new System object.
     else:
         try:
@@ -474,6 +473,14 @@ def __init__(
                     self._system, self._top_file, self._traj_file, self._property_map
                 )
 
+            # Copy the system.
+            self._renumbered_system = self._system.copy()
+
+            # Make sure the constituents of the system are numbered in ascending order.
+            self._renumbered_system._sire_object = _SireIO.renumberConstituents(
+                self._system._sire_object
+            )
+
         if not isinstance(backend, str):
             raise TypeError("'backend' must be of type 'str'")
 
@@ -849,7 +856,7 @@ def getFrames(self, indices=None):
                     new_system = _split_molecules(
                         frame,
                         pdb,
-                        self._system,
+                        self._renumbered_system,
                         str(self._work_dir),
                         self._property_map,
                     )
@@ -1115,6 +1122,10 @@ def _split_molecules(frame, pdb, reference, work_dir, property_map={}):
     # Store the formats associated with the reference system.
     formats = reference.fileFormat()
 
+    # Convert NoneType to empty list.
+    if formats is None:
+        formats = []
+
     # Write the frame coordinates/velocities to file.
     coord_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.coords")
     top_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.top")

python/BioSimSpace/Trajectory/_trajectory.py

@@ -153,6 +153,14 @@ def getFrame(trajectory, topology, index, system=None, property_map={}):
         # Update the water topology to match topology/trajectory.
         system = _update_water_topology(system, topology, trajectory, property_map)
 
+        # Copy the system.
+        renumbered_system = system.copy()
+
+        # Make sure the constituents of the system are numbered in ascending order.
+        renumbered_system._sire_object = _SireIO.renumberConstituents(
+            system._sire_object
+        )
+
     # Try to load the frame with Sire.
     errors = []
     is_sire = False
@@ -219,9 +227,7 @@ def getFrame(trajectory, topology, index, system=None, property_map={}):
                 if "space" in new_system.propertyKeys():
                     box = new_system.property("space")
                     if box.isPeriodic():
-                        sire_system.setProperty(
-                            self._property_map.get("space", "space"), box
-                        )
+                        sire_system.setProperty(property_map.get("space", "space"), box)
 
                 new_system = _System(sire_system)
 
@@ -265,7 +271,7 @@ def getFrame(trajectory, topology, index, system=None, property_map={}):
             # coordinates of all of the atoms in the reference. As such, we
             # will need to split the system into molecules.
             new_system = _split_molecules(
-                frame, pdb, system, str(work_dir), property_map
+                frame, pdb, renumbered_system, str(work_dir), property_map
             )
             try:
                 sire_system, _ = _SireIO.updateCoordinatesAndVelocities(
@@ -276,9 +282,7 @@ def getFrame(trajectory, topology, index, system=None, property_map={}):
                 if "space" in new_system.propertyKeys():
                     box = new_system.property("space")
                     if box.isPeriodic():
-                        sire_system.setProperty(
-                            self._property_map.get("space", "space"), box
-                        )
+                        sire_system.setProperty(property_map.get("space", "space"), box)
 
                 new_system = _System(sire_system)
             except Exception as e:
@@ -288,11 +292,6 @@ def getFrame(trajectory, topology, index, system=None, property_map={}):
                 else:
                     raise IOError(msg) from None
 
-            else:
-                raise IOError(
-                    "The trajectory frame is incompatible with the passed system!"
-                )
-
     # Load the frame directly to create a new System object.
     else:
         try:
@@ -474,6 +473,14 @@ def __init__(
                     self._system, self._top_file, self._traj_file, self._property_map
                 )
 
+            # Copy the system.
+            self._renumbered_system = self._system.copy()
+
+            # Make sure the constituents of the system are numbered in ascending order.
+            self._renumbered_system._sire_object = _SireIO.renumberConstituents(
+                self._system._sire_object
+            )
+
         if not isinstance(backend, str):
             raise TypeError("'backend' must be of type 'str'")
 
@@ -849,7 +856,7 @@ def getFrames(self, indices=None):
                     new_system = _split_molecules(
                         frame,
                         pdb,
-                        self._system,
+                        self._renumbered_system,
                         str(self._work_dir),
                         self._property_map,
                     )
@@ -1115,6 +1122,10 @@ def _split_molecules(frame, pdb, reference, work_dir, property_map={}):
     # Store the formats associated with the reference system.
     formats = reference.fileFormat()
 
+    # Convert NoneType to empty list.
+    if formats is None:
+        formats = []
+
     # Write the frame coordinates/velocities to file.
     coord_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.coords")
     top_file = _os.path.join(str(work_dir), f"{str(_uuid.uuid4())}.top")

tests/Sandpit/Exscientia/Trajectory/test_trajectory.py

@@ -170,3 +170,24 @@ def test_rmsd(traj_sire, traj_mdtraj, traj_mdanalysis):
     for v0, v1, v2 in zip(rmsd0, rmsd1, rmsd2):
         assert v0.value() == pytest.approx(v1.value(), abs=1e-2)
         assert v0.value() == pytest.approx(v2.value(), abs=1e-2)
+
+
+@pytest.mark.skipif(has_mdtraj is False, reason="Requires mdtraj to be installed.")
+def test_getFrame(system):
+    """Regression test to make sure the getFrame function works."""
+
+    # Just load the trajectory directly.
+    frame = BSS.Trajectory.getFrame(
+        "tests/input/ala.trr", topology="tests/input/ala.gro", index=5
+    )
+
+    assert frame.nAtoms() == system.nAtoms()
+
+    # Load using the system as a reference.
+    frame = BSS.Trajectory.getFrame(
+        "tests/input/ala.trr", "tests/input/ala.gro", system=system, index=5
+    )
+
+    assert frame.nMolecules() == system.nMolecules()
+    assert frame.nResidues() == system.nResidues()
+    assert frame.nAtoms() == system.nAtoms()

tests/Trajectory/test_trajectory.py

@@ -161,3 +161,24 @@ def test_rmsd(traj_sire, traj_mdtraj, traj_mdanalysis):
     for v0, v1, v2 in zip(rmsd0, rmsd1, rmsd2):
         assert v0.value() == pytest.approx(v1.value(), abs=1e-2)
         assert v0.value() == pytest.approx(v2.value(), abs=1e-2)
+
+
+@pytest.mark.skipif(has_mdtraj is False, reason="Requires mdtraj to be installed.")
+def test_getFrame(system):
+    """Regression test to make sure the getFrame function works."""
+
+    # Just load the trajectory directly.
+    frame = BSS.Trajectory.getFrame(
+        "tests/input/ala.trr", topology="tests/input/ala.gro", index=5
+    )
+
+    assert frame.nAtoms() == system.nAtoms()
+
+    # Load using the system as a reference.
+    frame = BSS.Trajectory.getFrame(
+        "tests/input/ala.trr", "tests/input/ala.gro", system=system, index=5
+    )
+
+    assert frame.nMolecules() == system.nMolecules()
+    assert frame.nResidues() == system.nResidues()
+    assert frame.nAtoms() == system.nAtoms()