Merge pull request #138 from OpenBioSim/fix-137

Fix 137

[ci skip]

doc/source/changelog.rst

@@ -66,6 +66,9 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
   means that we can catch out-of-memory errors and raise a more
   informative exception.
 
+* Fixed the bug where the wrong return type from ``.minimisation()`` and
+  ``.dynamics()`` was returned. This fixes issue #137.
+
 * Please add an item to this changelog when you create your PR
 
 `2023.4.1 <https://github.com/openbiosim/sire/compare/2023.4.0...2023.4.1>`__ - October 2023

doc/source/cheatsheet/openmm.rst

@@ -161,6 +161,9 @@ Available keys and allowable values are listed below.
 | timestep                     | Time between integration steps, e.g.                     |
 |                              | ``2 * sr.units.femtosecond``                             |
 +------------------------------+----------------------------------------------------------+
+| tolerance                    | The tolerance to use for the PME calculation, e.g.       |
+|                              | ``0.0001``                                               |
++------------------------------+----------------------------------------------------------+
 | use_dispersion_correction    | Whether or not to use the dispersion correction to       |
 |                              | deal with cutoff issues. This is very expensive.         |
 +------------------------------+----------------------------------------------------------+

src/sire/mol/_dynamics.py

@@ -13,6 +13,10 @@ def __init__(self, mols=None, map=None, **kwargs):
 
         map = create_map(map, kwargs)
 
+        # Save the original view, so that we can return an object
+        # of the same type in .commit()
+        self._orig_mols = mols
+
         if mols is not None:
             self._map = map  # this is already a PropertyMap
 
@@ -45,12 +49,8 @@ def __init__(self, mols=None, map=None, **kwargs):
             else:
                 # create a system to work on
                 self._sire_mols = System()
-                self._sire_mols._system.add(
-                    mols.molecules().to_molecule_group()
-                )
-                self._sire_mols._system.set_property(
-                    "space", self._ffinfo.space()
-                )
+                self._sire_mols._system.add(mols.molecules().to_molecule_group())
+                self._sire_mols._system.set_property("space", self._ffinfo.space())
 
             # find the existing energy trajectory - we will build on this
             self._energy_trajectory = self._sire_mols.energy_trajectory(
@@ -154,9 +154,7 @@ def _update_from(self, state, state_has_cv, nsteps_completed):
 
     def _enter_dynamics_block(self):
         if self._is_running:
-            raise SystemError(
-                "Cannot start dynamics while it is already running!"
-            )
+            raise SystemError("Cannot start dynamics while it is already running!")
 
         self._is_running = True
         self._omm_state = None
@@ -189,8 +187,7 @@ def _exit_dynamics_block(
             self._omm_state_has_cv = (False, False)
 
         current_time = (
-            self._omm_state.getTime().value_in_unit(openmm.unit.nanosecond)
-            * nanosecond
+            self._omm_state.getTime().value_in_unit(openmm.unit.nanosecond) * nanosecond
         )
 
         delta = current_time - self._elapsed_time
@@ -358,9 +355,7 @@ def set_temperature(self, temperature, rescale_velocities: bool = True):
 
         ensemble.set_temperature(temperature)
 
-        self.set_ensemble(
-            ensemble=ensemble, rescale_velocities=rescale_velocities
-        )
+        self.set_ensemble(ensemble=ensemble, rescale_velocities=rescale_velocities)
 
     def set_pressure(self, pressure):
         """
@@ -600,9 +595,7 @@ def run(
                 lambda_windows = [lambda_windows]
 
         try:
-            steps_to_run = int(
-                time.to(picosecond) / self.timestep().to(picosecond)
-            )
+            steps_to_run = int(time.to(picosecond) / self.timestep().to(picosecond))
         except Exception:
             # passed in the number of steps instead
             steps_to_run = int(time)
@@ -617,9 +610,7 @@ def run(
         if save_frequency != 0:
             if save_frequency is None:
                 if self._map.specified("save_frequency"):
-                    save_frequency = (
-                        self._map["save_frequency"].value().to(picosecond)
-                    )
+                    save_frequency = self._map["save_frequency"].value().to(picosecond)
                 else:
                     save_frequency = 25
             else:
@@ -672,13 +663,9 @@ def process_block(state, state_has_cv, nsteps_completed):
 
         completed = 0
 
-        frame_frequency_steps = int(
-            frame_frequency / self.timestep().to(picosecond)
-        )
+        frame_frequency_steps = int(frame_frequency / self.timestep().to(picosecond))
 
-        energy_frequency_steps = int(
-            energy_frequency / self.timestep().to(picosecond)
-        )
+        energy_frequency_steps = int(energy_frequency / self.timestep().to(picosecond))
 
         def get_steps_till_save(completed: int, total: int):
             """Internal function to calculate the number of steps
@@ -703,17 +690,15 @@ def get_steps_till_save(completed: int, total: int):
 
             if frame_frequency_steps > 0:
                 n_to_frame = min(
-                    frame_frequency_steps
-                    - (completed % frame_frequency_steps),
+                    frame_frequency_steps - (completed % frame_frequency_steps),
                     n_to_end,
                 )
             else:
                 n_to_frame = total - completed
 
             if energy_frequency_steps > 0:
                 n_to_energy = min(
-                    energy_frequency_steps
-                    - (completed % energy_frequency_steps),
+                    energy_frequency_steps - (completed % energy_frequency_steps),
                     n_to_end,
                 )
             else:
@@ -829,10 +814,8 @@ class NeedsMinimiseError(Exception):
 
                         saved_last_frame = False
 
-                        kinetic_energy = (
-                            state.getKineticEnergy().value_in_unit(
-                                openmm.unit.kilocalorie_per_mole
-                            )
+                        kinetic_energy = state.getKineticEnergy().value_in_unit(
+                            openmm.unit.kilocalorie_per_mole
                         )
 
                         ke_per_atom = kinetic_energy / self._num_atoms
@@ -855,9 +838,7 @@ class NeedsMinimiseError(Exception):
                                 "and run again."
                             )
 
-                self._walltime += (
-                    datetime.now() - start_time
-                ).total_seconds() * second
+                self._walltime += (datetime.now() - start_time).total_seconds() * second
 
             if state is not None and not saved_last_frame:
                 # we can process the last block in the main thread
@@ -883,19 +864,27 @@ def commit(self):
             return
 
         self._update_from(
-            state=self._get_current_state(
-                include_coords=True, include_velocities=True
-            ),
+            state=self._get_current_state(include_coords=True, include_velocities=True),
             state_has_cv=(True, True),
             nsteps_completed=self._current_step,
         )
 
-        self._sire_mols.set_energy_trajectory(
-            self._energy_trajectory, map=self._map
-        )
+        self._sire_mols.set_energy_trajectory(self._energy_trajectory, map=self._map)
 
         self._sire_mols.set_ensemble(self.ensemble())
 
+        if self._orig_mols is not None:
+            from ..system import System
+
+            if System.is_system(self._orig_mols):
+                return self._sire_mols
+            else:
+                r = self._orig_mols.clone()
+                r.update(self._sire_mols.molecules())
+                return r
+        else:
+            return self._sire_mols.clone()
+
 
 def _add_extra(extras, key, value):
     if value is not None:
@@ -1080,9 +1069,7 @@ def run(
         if not self._d.is_null():
             if save_velocities is None:
                 if self._d._map.specified("save_velocities"):
-                    save_velocities = (
-                        self._d._map["save_velocities"].value().as_bool()
-                    )
+                    save_velocities = self._d._map["save_velocities"].value().as_bool()
                 else:
                     save_velocities = False
 
@@ -1180,9 +1167,7 @@ def randomise_velocities(self, temperature=None, random_seed: int = None):
         - random_seed (int): The random seed to use. If None, then
           a random seed will be generated
         """
-        self._d.randomise_velocities(
-            temperature=temperature, random_seed=random_seed
-        )
+        self._d.randomise_velocities(temperature=temperature, random_seed=random_seed)
 
     def constraint(self):
         """
@@ -1341,9 +1326,7 @@ def current_kinetic_energy(self):
         """
         return self._d.current_kinetic_energy()
 
-    def energy_trajectory(
-        self, to_pandas: bool = False, to_alchemlyb: bool = False
-    ):
+    def energy_trajectory(self, to_pandas: bool = False, to_alchemlyb: bool = False):
         """
         Return the energy trajectory. This is the trajectory of
         energy values that have been captured during dynamics.
@@ -1363,9 +1346,9 @@ def energy_trajectory(
 
     def commit(self):
         if not self._d.is_null():
-            self._d.commit()
-
-        return self._d._sire_mols
+            return self._d.commit()
+        else:
+            return None
 
     def __call__(
         self,

src/sire/mol/_minimisation.py

@@ -70,9 +70,9 @@ def commit(self):
         minimised molecules.
         """
         if not self._d.is_null():
-            self._d.commit()
-
-        return self._d._sire_mols
+            return self._d.commit()
+        else:
+            return None
 
     def __call__(self, max_iterations: int = 10000):
         """

tests/mol/test_dynamics.py

@@ -0,0 +1,50 @@
+import sire as sr
+import pytest
+
+
+@pytest.mark.skipif(
+    "openmm" not in sr.convert.supported_formats(),
+    reason="openmm support is not available",
+)
+def test_dynamics_return_type(ala_mols):
+    mols = ala_mols
+
+    m = mols.dynamics(platform="Reference").run("1 fs").commit()
+
+    assert isinstance(m, type(mols))
+
+    mol = mols[0]
+
+    m = mol.dynamics(platform="Reference").run("1 fs").commit()
+
+    assert isinstance(m, type(mol))
+
+    atom = mol[0]
+
+    a = atom.dynamics(platform="Reference").run("1 fs").commit()
+
+    assert isinstance(a, type(atom))
+
+
+@pytest.mark.skipif(
+    "openmm" not in sr.convert.supported_formats(),
+    reason="openmm support is not available",
+)
+def test_minimisation_return_type(ala_mols):
+    mols = ala_mols
+
+    m = mols.minimisation(platform="Reference").run(1).commit()
+
+    assert isinstance(m, type(mols))
+
+    mol = mols[0]
+
+    m = mol.minimisation(platform="Reference").run(1).commit()
+
+    assert isinstance(m, type(mol))
+
+    atom = mol[0]
+
+    a = atom.minimisation(platform="Reference").run(1).commit()
+
+    assert isinstance(a, type(atom))