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 b22d18dc5..72c3cd4be 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py +++ b/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 diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py index 8137a1478..c6f0e52ed 100644 --- a/python/BioSimSpace/_SireWrappers/_system.py +++ b/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 diff --git a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py index a0f78b116..ed0391534 100644 --- a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py +++ b/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") diff --git a/tests/_SireWrappers/test_system.py b/tests/_SireWrappers/test_system.py index 3f8cdb126..e715b637f 100644 --- a/tests/_SireWrappers/test_system.py +++ b/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")