diff --git a/doc/source/tutorials/crystal_water.rst b/doc/source/tutorials/crystal_water.rst
index 32a148815..2e957b023 100644
--- a/doc/source/tutorials/crystal_water.rst
+++ b/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
diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py
index 30f7915e1..c1fed655f 100644
--- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py
+++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py
@@ -640,16 +640,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():
@@ -664,7 +663,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"
@@ -2384,7 +2383,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
diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py
index 1820bcd76..8208eb0da 100644
--- a/python/BioSimSpace/_SireWrappers/_system.py
+++ b/python/BioSimSpace/_SireWrappers/_system.py
@@ -640,16 +640,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():
@@ -664,7 +663,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"
@@ -2304,7 +2303,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
diff --git a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py
index 8e5d6b596..ee98bfac0 100644
--- a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py
+++ b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py
@@ -473,3 +473,24 @@ def test_residue_searches_rna(rna_system):
     assert isinstance(nucleotides[0], BSS._SireWrappers._residue.Residue)
     assert len(rna_system.getAminoAcids()) == rna_system.nAminoAcids() == 0
     assert len(nucleotides) == rna_system.nNucleotides() == 29
+
+
+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")
diff --git a/tests/_SireWrappers/test_system.py b/tests/_SireWrappers/test_system.py
index a7234cb00..81ff55a5f 100644
--- a/tests/_SireWrappers/test_system.py
+++ b/tests/_SireWrappers/test_system.py
@@ -469,3 +469,24 @@ def test_residue_searches_rna(rna_system):
     assert isinstance(nucleotides[0], BSS._SireWrappers._residue.Residue)
     assert len(rna_system.getAminoAcids()) == rna_system.nAminoAcids() == 0
     assert len(nucleotides) == rna_system.nNucleotides() == 29
+
+
+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")