Merge pull request #198 from OpenBioSim/backport_197

Backport fixes from PR #197

doc/source/tutorials/crystal_water.rst

@@ -52,7 +52,7 @@ protein:
 
 And then the water:
 
->>> xtal_water = BSS.Parameters.tip3p(
+>>> xtal_water = BSS.Parameters.ff14SB(
 ...     xtal_water,
 ...     water_model="tip3p",
 ...     ensure_compatible=False

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

@@ -639,16 +639,15 @@ def addMolecules(self, molecules):
 
             # Search for perturbable molecules with a velocity property.
             # Only consider the lambda = 0 end state.
-            try:
-                pert_mols_with_velocities = self.search(
-                    f"mols with property velocity0"
-                ).molecules()
-                num_pert_vels = len(pert_mols_with_velocities)
-            except:
-                num_pert_vels = 0
-
-            # Compute the total number of molecules with velocities.
-            num_vels = num_vels + num_pert_vels
+            has_pertrubable = False
+            for mol in self.getPerturbableMolecules():
+                # Add perturbable velocities.
+                if mol._sire_object.hasProperty("velocity0"):
+                    has_pertrubable = True
+                    num_vels += 1
+                # Remove non-perturbable velocities to avoid double counting.
+                elif mol._sire_object.hasProperty("velocity"):
+                    num_vels -= 1
 
             # Not all molecules have velocities.
             if num_vels > 0 and num_vels != self.nMolecules():
@@ -663,7 +662,7 @@ def addMolecules(self, molecules):
                     _warnings.warn(
                         "Failed to remove 'velocity' property from all molecules!"
                     )
-                if num_pert_vels > 0:
+                if has_perturbable:
                     try:
                         self._sire_object = _SireIO.removeProperty(
                             self._sire_object, "velocity0"
@@ -2484,7 +2483,8 @@ def _set_water_topology(self, format, is_crystal=False, property_map={}):
             The format to convert to: either "AMBER" or "GROMACS".
 
         is_crystal : bool
-            Whether to label as rystal waters.
+            Whether to label as crystal waters. This is only used when solvating
+            so that crystal waters aren't removed when adding ions.
 
         property_map : dict
            A dictionary that maps system "properties" to their user defined

python/BioSimSpace/_SireWrappers/_system.py

@@ -639,16 +639,15 @@ def addMolecules(self, molecules):
 
             # Search for perturbable molecules with a velocity property.
             # Only consider the lambda = 0 end state.
-            try:
-                pert_mols_with_velocities = self.search(
-                    f"mols with property velocity0"
-                ).molecules()
-                num_pert_vels = len(pert_mols_with_velocities)
-            except:
-                num_pert_vels = 0
-
-            # Compute the total number of molecules with velocities.
-            num_vels = num_vels + num_pert_vels
+            has_perturbable = False
+            for mol in self.getPerturbableMolecules():
+                # Add perturbable velocities.
+                if mol._sire_object.hasProperty("velocity0"):
+                    has_perturbable = True
+                    num_vels += 1
+                # Remove non-perturbable velocities to avoid double counting.
+                elif mol._sire_object.hasProperty("velocity"):
+                    num_vels -= 1
 
             # Not all molecules have velocities.
             if num_vels > 0 and num_vels != self.nMolecules():
@@ -663,7 +662,7 @@ def addMolecules(self, molecules):
                     _warnings.warn(
                         "Failed to remove 'velocity' property from all molecules!"
                     )
-                if num_pert_vels > 0:
+                if has_perturnable:
                     try:
                         self._sire_object = _SireIO.removeProperty(
                             self._sire_object, "velocity0"
@@ -2405,7 +2404,8 @@ def _set_water_topology(self, format, is_crystal=False, property_map={}):
             The format to convert to: either "AMBER" or "GROMACS".
 
         is_crystal : bool
-            Whether to label as rystal waters.
+            Whether to label as crystal waters. This is only used when solvating
+            so that crystal waters aren't removed when adding ions.
 
         property_map : dict
            A dictionary that maps system "properties" to their user defined

tests/Sandpit/Exscientia/_SireWrappers/test_system.py

@@ -451,3 +451,24 @@ def test_rotate_box_vectors(system):
         assert math.isclose(c1.x().value(), c0.z().value())
         assert math.isclose(c1.y().value(), c0.x().value())
         assert math.isclose(c1.z().value(), c0.y().value())
+
+
+def test_set_water_property_preserve(system):
+    # Make sure that unique molecular properties are preserved when swapping
+    # water topology.
+
+    # Make a copy of the system.
+    system = system.copy()
+
+    # Flag one water molecule with a unique property.
+    mol = system[-1]
+    mol._sire_object = (
+        mol._sire_object.edit().setProperty("test", True).molecule().commit()
+    )
+    system.updateMolecules(mol)
+
+    # Swap the water topology to GROMACS format.
+    system._set_water_topology("GROMACS")
+
+    # Make sure the property is preserved.
+    assert system[-1]._sire_object.hasProperty("test")

tests/_SireWrappers/test_system.py

@@ -447,3 +447,24 @@ def test_rotate_box_vectors(system):
         assert math.isclose(c1.x().value(), c0.z().value())
         assert math.isclose(c1.y().value(), c0.x().value())
         assert math.isclose(c1.z().value(), c0.y().value())
+
+
+def test_set_water_property_preserve(system):
+    # Make sure that unique molecular properties are preserved when swapping
+    # water topology.
+
+    # Make a copy of the system.
+    system = system.copy()
+
+    # Flag one water molecule with a unique property.
+    mol = system[-1]
+    mol._sire_object = (
+        mol._sire_object.edit().setProperty("test", True).molecule().commit()
+    )
+    system.updateMolecules(mol)
+
+    # Swap the water topology to GROMACS format.
+    system._set_water_topology("GROMACS")
+
+    # Make sure the property is preserved.
+    assert system[-1]._sire_object.hasProperty("test")