From 7b80ccf3e5946c93087ff59ad306cb782bedb636 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 24 Oct 2023 14:02:39 +0100 Subject: [PATCH 01/20] Backport fix from PR #190. [ci skip] --- .../BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py | 1 + python/BioSimSpace/Solvent/_solvent.py | 1 + tests/Sandpit/Exscientia/Solvent/test_solvent.py | 10 ++++++++-- tests/Solvent/test_solvent.py | 10 ++++++++-- 4 files changed, 18 insertions(+), 4 deletions(-) diff --git a/python/BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py b/python/BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py index 7e3d4e944..0bcb4a682 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Solvent/_solvent.py @@ -142,6 +142,7 @@ def solvate( ion_conc, is_neutral, is_aligned, + match_water, work_dir, property_map, ) diff --git a/python/BioSimSpace/Solvent/_solvent.py b/python/BioSimSpace/Solvent/_solvent.py index 7e3d4e944..0bcb4a682 100644 --- a/python/BioSimSpace/Solvent/_solvent.py +++ b/python/BioSimSpace/Solvent/_solvent.py @@ -142,6 +142,7 @@ def solvate( ion_conc, is_neutral, is_aligned, + match_water, work_dir, property_map, ) diff --git a/tests/Sandpit/Exscientia/Solvent/test_solvent.py b/tests/Sandpit/Exscientia/Solvent/test_solvent.py index f3a0b8978..7e6f7de25 100644 --- a/tests/Sandpit/Exscientia/Solvent/test_solvent.py +++ b/tests/Sandpit/Exscientia/Solvent/test_solvent.py @@ -1,6 +1,8 @@ import pytest import tempfile +from functools import partial + from sire.legacy.IO import GroTop import BioSimSpace.Sandpit.Exscientia as BSS @@ -18,8 +20,12 @@ def system(): @pytest.mark.parametrize("match_water", [True, False]) +@pytest.mark.parametrize( + "function", + [partial(BSS.Solvent.solvate, "tip3p"), BSS.Solvent.tip3p], +) @pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") -def test_crystal_water(system, match_water): +def test_crystal_water(system, match_water, function): """ Test that user defined crystal waters can be preserved during solvation and on write to GroTop format. @@ -35,7 +41,7 @@ def test_crystal_water(system, match_water): box, angles = BSS.Box.cubic(5.5 * BSS.Units.Length.nanometer) # Create the solvated system. - solvated = BSS.Solvent.tip3p(system, box, angles, match_water=match_water) + solvated = function(system, box, angles, match_water=match_water) # Search for the crystal waters in the solvated system. try: diff --git a/tests/Solvent/test_solvent.py b/tests/Solvent/test_solvent.py index 7edf20cbf..185bb6373 100644 --- a/tests/Solvent/test_solvent.py +++ b/tests/Solvent/test_solvent.py @@ -1,6 +1,8 @@ import pytest import tempfile +from functools import partial + from sire.legacy.IO import GroTop import BioSimSpace as BSS @@ -18,8 +20,12 @@ def system(): @pytest.mark.parametrize("match_water", [True, False]) +@pytest.mark.parametrize( + "function", + [partial(BSS.Solvent.solvate, "tip3p"), BSS.Solvent.tip3p], +) @pytest.mark.skipif(not has_gromacs, reason="Requires GROMACS to be installed") -def test_crystal_water(system, match_water): +def test_crystal_water(system, match_water, function): """ Test that user defined crystal waters can be preserved during solvation and on write to GroTop format. @@ -35,7 +41,7 @@ def test_crystal_water(system, match_water): box, angles = BSS.Box.cubic(5.5 * BSS.Units.Length.nanometer) # Create the solvated system. - solvated = BSS.Solvent.tip3p(system, box, angles, match_water=match_water) + solvated = function(system, box, angles, match_water=match_water) # Search for the crystal waters in the solvated system. try: From 45e58975086c083b8c9042cecbd4735f4f57a899 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 31 Oct 2023 08:54:32 +0000 Subject: [PATCH 02/20] Backport fix from PR #193. [ci skip] --- .../Exscientia/_SireWrappers/_system.py | 27 ++++++++++++++++++- python/BioSimSpace/_SireWrappers/_system.py | 27 ++++++++++++++++++- 2 files changed, 52 insertions(+), 2 deletions(-) diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py index 0cb99ef12..b22d18dc5 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py @@ -628,7 +628,7 @@ def addMolecules(self, molecules): # Remove velocities if any molecules are missing them. if self.nMolecules() > 1: - # Search for water molecules in the system. + # Search for molecules with a velocity property. try: mols_with_velocities = self.search( f"mols with property velocity" @@ -637,6 +637,19 @@ def addMolecules(self, molecules): except: num_vels = 0 + # 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 + # Not all molecules have velocities. if num_vels > 0 and num_vels != self.nMolecules(): _warnings.warn( @@ -650,6 +663,18 @@ def addMolecules(self, molecules): _warnings.warn( "Failed to remove 'velocity' property from all molecules!" ) + if num_pert_vels > 0: + try: + self._sire_object = _SireIO.removeProperty( + self._sire_object, "velocity0" + ) + self._sire_object = _SireIO.removeProperty( + self._sire_object, "velocity1" + ) + except: + _warnings.warn( + "Failed to remove 'velocity0' and 'velocity1' property from molecules!" + ) def removeMolecules(self, molecules): """ diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py index d8200e55e..8137a1478 100644 --- a/python/BioSimSpace/_SireWrappers/_system.py +++ b/python/BioSimSpace/_SireWrappers/_system.py @@ -628,7 +628,7 @@ def addMolecules(self, molecules): # Remove velocities if any molecules are missing them. if self.nMolecules() > 1: - # Search for water molecules in the system. + # Search for molecules with a velocity property. try: mols_with_velocities = self.search( f"mols with property velocity" @@ -637,6 +637,19 @@ def addMolecules(self, molecules): except: num_vels = 0 + # 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 + # Not all molecules have velocities. if num_vels > 0 and num_vels != self.nMolecules(): _warnings.warn( @@ -650,6 +663,18 @@ def addMolecules(self, molecules): _warnings.warn( "Failed to remove 'velocity' property from all molecules!" ) + if num_pert_vels > 0: + try: + self._sire_object = _SireIO.removeProperty( + self._sire_object, "velocity0" + ) + self._sire_object = _SireIO.removeProperty( + self._sire_object, "velocity1" + ) + except: + _warnings.warn( + "Failed to remove 'velocity0' and 'velocity1' property from molecules!" + ) def removeMolecules(self, molecules): """ From 9226af3834a7d165c1a48a42266101ce2205b115 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 8 Nov 2023 10:20:50 +0000 Subject: [PATCH 03/20] Backport fixes from PR #197. [ci skip] --- doc/source/tutorials/crystal_water.rst | 2 +- .../Exscientia/_SireWrappers/_system.py | 24 +++++++++---------- python/BioSimSpace/_SireWrappers/_system.py | 24 +++++++++---------- .../Exscientia/_SireWrappers/test_system.py | 21 ++++++++++++++++ tests/_SireWrappers/test_system.py | 21 ++++++++++++++++ 5 files changed, 67 insertions(+), 25 deletions(-) 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") From 8f3309c677155d78a50467dd36baa310ccd5d8a3 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 13 Nov 2023 09:36:33 +0000 Subject: [PATCH 04/20] Backport fix from PR #201. [ci skip] --- .../Sandpit/Exscientia/Types/_type.py | 20 ------------------- python/BioSimSpace/Types/_type.py | 20 ------------------- tests/Sandpit/Exscientia/Types/test_types.py | 18 +++++++++++++++++ tests/Types/test_types.py | 18 +++++++++++++++++ 4 files changed, 36 insertions(+), 40 deletions(-) diff --git a/python/BioSimSpace/Sandpit/Exscientia/Types/_type.py b/python/BioSimSpace/Sandpit/Exscientia/Types/_type.py index 25912a900..f01635631 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Types/_type.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Types/_type.py @@ -388,26 +388,6 @@ def __eq__(self, other): else: return False - def __ne__(self, other): - """Not equals to operator.""" - - # Compare to another object of the same type. - if type(other) is type(self): - return _math.isclose( - self._to_default_unit().value(), - other._to_default_unit().value(), - ) - - # Compare with a string. - elif isinstance(other, str): - return not _math.isclose( - self._to_default_unit().value(), - self._from_string(other)._to_default_unit().value(), - ) - - else: - return True - def __ge__(self, other): """Greater than or equal to operator.""" diff --git a/python/BioSimSpace/Types/_type.py b/python/BioSimSpace/Types/_type.py index 25912a900..f01635631 100644 --- a/python/BioSimSpace/Types/_type.py +++ b/python/BioSimSpace/Types/_type.py @@ -388,26 +388,6 @@ def __eq__(self, other): else: return False - def __ne__(self, other): - """Not equals to operator.""" - - # Compare to another object of the same type. - if type(other) is type(self): - return _math.isclose( - self._to_default_unit().value(), - other._to_default_unit().value(), - ) - - # Compare with a string. - elif isinstance(other, str): - return not _math.isclose( - self._to_default_unit().value(), - self._from_string(other)._to_default_unit().value(), - ) - - else: - return True - def __ge__(self, other): """Greater than or equal to operator.""" diff --git a/tests/Sandpit/Exscientia/Types/test_types.py b/tests/Sandpit/Exscientia/Types/test_types.py index f7549286b..b3cc909f5 100644 --- a/tests/Sandpit/Exscientia/Types/test_types.py +++ b/tests/Sandpit/Exscientia/Types/test_types.py @@ -134,3 +134,21 @@ def test_container_mul(Type): # Assert we have a container of the original type. assert all(isinstance(x, Type) for x in container) + + +@pytest.mark.parametrize("Type", types) +def test_operators(Type): + """Test equality/inequality operators.""" + + t0 = Type(1.0, Type._default_unit) + t1 = Type(2.0, Type._default_unit) + t2 = Type(3.0, Type._default_unit) + + assert t0 == t0 + assert t0 <= t0 + assert t0 >= t0 + assert t0 != t1 + assert t0 < t1 + assert t1 > t0 + assert t1 < t2 + assert t2 > t1 diff --git a/tests/Types/test_types.py b/tests/Types/test_types.py index 32098b77f..5aa9a7b78 100644 --- a/tests/Types/test_types.py +++ b/tests/Types/test_types.py @@ -134,3 +134,21 @@ def test_container_mul(Type): # Assert we have a container of the original type. assert all(isinstance(x, Type) for x in container) + + +@pytest.mark.parametrize("Type", types) +def test_operators(Type): + """Test equality/inequality operators.""" + + t0 = Type(1.0, Type._default_unit) + t1 = Type(2.0, Type._default_unit) + t2 = Type(3.0, Type._default_unit) + + assert t0 == t0 + assert t0 <= t0 + assert t0 >= t0 + assert t0 != t1 + assert t0 < t1 + assert t1 > t0 + assert t1 < t2 + assert t2 > t1 From 954c5459c5808e192343fabf0850f77354a901e2 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 14 Nov 2023 09:52:18 +0000 Subject: [PATCH 05/20] Backport fixes from PR #203. [ci skip] --- .../Sandpit/Exscientia/Trajectory/_trajectory.py | 5 ++--- .../Sandpit/Exscientia/_SireWrappers/_system.py | 14 ++++++++++++-- python/BioSimSpace/Trajectory/_trajectory.py | 5 ++--- python/BioSimSpace/_SireWrappers/_system.py | 14 ++++++++++++-- 4 files changed, 28 insertions(+), 10 deletions(-) diff --git a/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py b/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py index e3d747ca0..b49a05d7d 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py @@ -1120,11 +1120,10 @@ def _split_molecules(frame, pdb, reference, work_dir, property_map={}): # Create a triclinic space from the information in the frame file. if isinstance(frame, _SireIO.AmberRst7): - # Get the box dimensions and angles. Take the values, since the - # units are wrong. + # Get the box dimensions and angles. degree = _SireUnits.degree dimensions = [x.value() for x in frame.box_dimensions()] - angles = [x.value() * degree for x in frame.box_angles()] + angles = [x.to(degree) * degree for x in frame.box_angles()] box = _SireVol.TriclinicBox(*dimensions, *angles) else: box = _SireVol.TriclinicBox(frame.box_v1(), frame.box_v2(), frame.box_v3()) diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py index 72c3cd4be..cafb3dde0 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py @@ -1218,7 +1218,10 @@ def rotateBoxVectors( try: prop_name = property_map.get("coordinates", "coordinates") cursor = cursor.rotate( - center=center, matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, + matrix=rotation_matrix, + rotate_velocities=False, + map={"coordinates": prop_name}, ) except: pass @@ -1227,7 +1230,10 @@ def rotateBoxVectors( try: prop_name = property_map.get("velocity", "velocity") cursor = cursor.rotate( - center=center, matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, + matrix=rotation_matrix, + rotate_velocities=False, + map={"coordinates": prop_name}, ) except: pass @@ -1240,12 +1246,14 @@ def rotateBoxVectors( cursor = cursor.rotate( center=center, matrix=rotation_matrix, + rotate_velocities=False, map={"coordinates": prop_name}, ) prop_name = property_map.get("coordinates", "coordinates") + "1" cursor = cursor.rotate( center=center, matrix=rotation_matrix, + rotate_velocities=False, map={"coordinates": prop_name}, ) except: @@ -1257,12 +1265,14 @@ def rotateBoxVectors( cursor = cursor.rotate( center=center, matrix=rotation_matrix, + rotate_velocities=False, map={"coordinates": prop_name}, ) prop_name = property_map.get("velocity", "velocity") + "1" cursor = cursor.rotate( center=center, matrix=rotation_matrix, + rotate_velocities=False, map={"coordinates": prop_name}, ) except: diff --git a/python/BioSimSpace/Trajectory/_trajectory.py b/python/BioSimSpace/Trajectory/_trajectory.py index e3d747ca0..b49a05d7d 100644 --- a/python/BioSimSpace/Trajectory/_trajectory.py +++ b/python/BioSimSpace/Trajectory/_trajectory.py @@ -1120,11 +1120,10 @@ def _split_molecules(frame, pdb, reference, work_dir, property_map={}): # Create a triclinic space from the information in the frame file. if isinstance(frame, _SireIO.AmberRst7): - # Get the box dimensions and angles. Take the values, since the - # units are wrong. + # Get the box dimensions and angles. degree = _SireUnits.degree dimensions = [x.value() for x in frame.box_dimensions()] - angles = [x.value() * degree for x in frame.box_angles()] + angles = [x.to(degree) * degree for x in frame.box_angles()] box = _SireVol.TriclinicBox(*dimensions, *angles) else: box = _SireVol.TriclinicBox(frame.box_v1(), frame.box_v2(), frame.box_v3()) diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py index c6f0e52ed..0a5bf5157 100644 --- a/python/BioSimSpace/_SireWrappers/_system.py +++ b/python/BioSimSpace/_SireWrappers/_system.py @@ -1166,7 +1166,10 @@ def rotateBoxVectors( try: prop_name = property_map.get("coordinates", "coordinates") cursor = cursor.rotate( - center=center, matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, + matrix=rotation_matrix, + rotate_velocities=False, + map={"coordinates": prop_name}, ) except: pass @@ -1175,7 +1178,10 @@ def rotateBoxVectors( try: prop_name = property_map.get("velocity", "velocity") cursor = cursor.rotate( - center=center, matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, + matrix=rotation_matrix, + rotate_velocities=False, + map={"coordinates": prop_name}, ) except: pass @@ -1188,12 +1194,14 @@ def rotateBoxVectors( cursor = cursor.rotate( center=center, matrix=rotation_matrix, + rotate_velocities=False, map={"coordinates": prop_name}, ) prop_name = property_map.get("coordinates", "coordinates") + "1" cursor = cursor.rotate( center=center, matrix=rotation_matrix, + rotate_velocities=False, map={"coordinates": prop_name}, ) except: @@ -1205,12 +1213,14 @@ def rotateBoxVectors( cursor = cursor.rotate( center=center, matrix=rotation_matrix, + rotate_velocities=False, map={"coordinates": prop_name}, ) prop_name = property_map.get("velocity", "velocity") + "1" cursor = cursor.rotate( center=center, matrix=rotation_matrix, + rotate_velocities=False, map={"coordinates": prop_name}, ) except: From e821df21826985a66c43ad9e7498a4ab39be96a4 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 14 Nov 2023 10:21:43 +0000 Subject: [PATCH 06/20] Backport fix from PR #204. [ci skip] --- .../Sandpit/Exscientia/FreeEnergy/_restraint_search.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py b/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py index 1eb666d1a..b4da8e8d3 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py +++ b/python/BioSimSpace/Sandpit/Exscientia/FreeEnergy/_restraint_search.py @@ -1374,7 +1374,6 @@ def _findOrderedBoresch( pair_list[:no_pairs], desc="Scoring candidate Boresch anchor points. Anchor set no: ", ): - boresch_dof_data[pair] = {} l1_idx, r1_idx = pair try: _, l2_idx, l3_idx = _getAnchorAts(l1_idx, ligand_selection_str, u) @@ -1383,6 +1382,8 @@ def _findOrderedBoresch( _AnalysisError ): # Failed to find full set of anchor points for this pair continue + + boresch_dof_data[pair] = {} boresch_dof_data[pair]["anchor_ats"] = [ l1_idx, l2_idx, From 2258775804411be603fef05df4573c953bf81586 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 14 Nov 2023 10:22:19 +0000 Subject: [PATCH 07/20] Fix typos in CHANGELOG. [ci skip] --- doc/source/changelog.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 782fab3e9..a78abcd71 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -17,7 +17,7 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Add unified free-energy perturbation analysis using ``alchemlyb`` (`@annamherz `_) (`#155 `__). * Fix handling of connectivity changes during molecular perturbations (`#157 `__). * Fix issues related to new shared properties in Sire (`#160 `__). -* Fix issues in SOMD perturbation files for absolute binding free-energy simulations (`@fkclark `_) (`#164 `__). +* Fix issues in SOMD perturbation files for absolute binding free-energy simulations (`@fjclark `_) (`#164 `__). * Don't generate velocities when performing a continuation with GROMACS (`#169 `__). * Decouple custom parameters and additional commands in ``LEaP`` input (`#170 `__). * Check for periodic space when updating box information from restart file or trajectory (`#173 `__). @@ -49,9 +49,9 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Fix order of imports in ``prepareFEP`` node (`#90 `__). * Recenter molecules following vacuum simulation with GROMACS to avoid precision overflow with molecular coordinates on write (`#95 `__). * Fix expected angles used in unit test following updates to triclinic box code in Sire (`#99 `__). -* Add absolute binding free-energy support for SOMD (`@fkclark `_) (`#104 `__). +* Add absolute binding free-energy support for SOMD (`@fjclark `_) (`#104 `__). * Avoid streaming issues when reading binary AMBER restart files for a single frame (`#105 `__). -* Improve overlap matrix plotting functionality (`@fkclark `_) (`#107 `__). +* Improve overlap matrix plotting functionality (`@fjclark `_) (`#107 `__). * Handle updates to Sire parser format naming (`#108 `__). * Wrap new Sire units grammar to improve parsing of units from strings (`#109 `__). * Expose ``make_whole`` option in Sire to allow un-wrapping of molecular coordinates on read (`#110 `__). @@ -60,7 +60,7 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Fix bug in ``plumed`` version requirement check (`#113 `__). * Reinstate temperature control for all GROMACS simulation protocols (`#115 `__). * Fix pre-processing selector in test section of ``conda`` recipe (`#117 `__). -* Fixed bug in SOMD free-energy perturbation analysis (`@fkclark `_) (`#119 `__). +* Fixed bug in SOMD free-energy perturbation analysis (`@fjclark `_) (`#119 `__). * Catch exception when vacuum system has a cartesian space (`#120 `__). * Add support for Sire as a trajectory backend (`#121 `__). From 89be7f7db0015186e634de2ce5f1313c7e663954 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 21 Nov 2023 11:35:53 +0000 Subject: [PATCH 08/20] Fix a typo introduced in PR #198. [ci skip] --- .../BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py index cafb3dde0..67e0c3d92 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py @@ -639,11 +639,11 @@ def addMolecules(self, molecules): # Search for perturbable molecules with a velocity property. # Only consider the lambda = 0 end state. - has_pertrubable = False + has_perturbable = False for mol in self.getPerturbableMolecules(): # Add perturbable velocities. if mol._sire_object.hasProperty("velocity0"): - has_pertrubable = True + has_perturbable = True num_vels += 1 # Remove non-perturbable velocities to avoid double counting. elif mol._sire_object.hasProperty("velocity"): From 08d9f85feff5d632ebd6cce650ebf6a63dc40185 Mon Sep 17 00:00:00 2001 From: "William (Zhiyi) Wu" Date: Wed, 18 Oct 2023 19:07:57 +0800 Subject: [PATCH 09/20] Sire 2023.4 --- .github/workflows/Sandpit_exs.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/Sandpit_exs.yml b/.github/workflows/Sandpit_exs.yml index 4f7f9fad2..76b9a3242 100644 --- a/.github/workflows/Sandpit_exs.yml +++ b/.github/workflows/Sandpit_exs.yml @@ -35,7 +35,7 @@ jobs: - name: Install dependency run: | - mamba install -c conda-forge -c openbiosim/label/main biosimspace python=3.10 ambertools gromacs "sire=2023.3" "alchemlyb>=2.1" pytest openff-interchange pint=0.21 rdkit "jaxlib>0.3.7" tqdm + mamba install -c conda-forge -c openbiosim/label/main biosimspace python=3.10 ambertools gromacs "sire=2023.4" "alchemlyb>=2.1" pytest openff-interchange pint=0.21 rdkit "jaxlib>0.3.7" tqdm python -m pip install git+https://github.com/Exscientia/MDRestraintsGenerator.git # For the testing of BSS.FreeEnergy.AlchemicalFreeEnergy.analysis python -m pip install https://github.com/alchemistry/alchemtest/archive/master.zip From a6296e432e1c776abe41bc900e7febd6e6c04fdf Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 4 Dec 2023 10:15:31 +0000 Subject: [PATCH 10/20] Backport fixes from PR #210. --- python/BioSimSpace/IO/_io.py | 25 +++++++++++++++++-- .../BioSimSpace/Sandpit/Exscientia/IO/_io.py | 25 +++++++++++++++++-- .../Exscientia/_SireWrappers/_molecule.py | 20 +++++++++------ python/BioSimSpace/_SireWrappers/_molecule.py | 20 +++++++++------ .../Exscientia/_SireWrappers/test_molecule.py | 24 ++++++++++++++++++ tests/_SireWrappers/test_molecule.py | 24 ++++++++++++++++++ 6 files changed, 118 insertions(+), 20 deletions(-) diff --git a/python/BioSimSpace/IO/_io.py b/python/BioSimSpace/IO/_io.py index fce3193e0..705e8f9d6 100644 --- a/python/BioSimSpace/IO/_io.py +++ b/python/BioSimSpace/IO/_io.py @@ -1100,9 +1100,13 @@ def readPerturbableSystem(top0, coords0, top1, coords1, property_map={}): # Flag that the molecule is perturbable. mol.setProperty("is_perturbable", _SireBase.wrap(True)) + # Get the two molecules. + mol0 = system0[idx]._sire_object + mol1 = system1[idx]._sire_object + # Add the molecule0 and molecule1 properties. - mol.setProperty("molecule0", system0[idx]._sire_object) - mol.setProperty("molecule1", system1[idx]._sire_object) + mol.setProperty("molecule0", mol0) + mol.setProperty("molecule1", mol1) # Get the connectivity property name. conn_prop = property_map.get("connectivity", "connectivity") @@ -1121,6 +1125,23 @@ def readPerturbableSystem(top0, coords0, top1, coords1, property_map={}): mol = mol.removeProperty(conn_prop + "0").molecule() mol = mol.removeProperty(conn_prop + "1").molecule() + # Reconstruct the intrascale matrices using the GroTop parser. + intra0 = ( + _SireIO.GroTop(_Molecule(mol0).toSystem()._sire_object) + .toSystem()[0] + .property("intrascale") + ) + intra1 = ( + _SireIO.GroTop(_Molecule(mol1).toSystem()._sire_object) + .toSystem()[0] + .property("intrascale") + ) + + # Set the "intrascale" properties. + intrascale_prop = property_map.get("intrascale", "intrascale") + mol.setProperty(intrascale_prop + "0", intra0) + mol.setProperty(intrascale_prop + "1", intra0) + # Commit the changes. mol = _Molecule(mol.commit()) diff --git a/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py b/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py index fce3193e0..705e8f9d6 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py +++ b/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py @@ -1100,9 +1100,13 @@ def readPerturbableSystem(top0, coords0, top1, coords1, property_map={}): # Flag that the molecule is perturbable. mol.setProperty("is_perturbable", _SireBase.wrap(True)) + # Get the two molecules. + mol0 = system0[idx]._sire_object + mol1 = system1[idx]._sire_object + # Add the molecule0 and molecule1 properties. - mol.setProperty("molecule0", system0[idx]._sire_object) - mol.setProperty("molecule1", system1[idx]._sire_object) + mol.setProperty("molecule0", mol0) + mol.setProperty("molecule1", mol1) # Get the connectivity property name. conn_prop = property_map.get("connectivity", "connectivity") @@ -1121,6 +1125,23 @@ def readPerturbableSystem(top0, coords0, top1, coords1, property_map={}): mol = mol.removeProperty(conn_prop + "0").molecule() mol = mol.removeProperty(conn_prop + "1").molecule() + # Reconstruct the intrascale matrices using the GroTop parser. + intra0 = ( + _SireIO.GroTop(_Molecule(mol0).toSystem()._sire_object) + .toSystem()[0] + .property("intrascale") + ) + intra1 = ( + _SireIO.GroTop(_Molecule(mol1).toSystem()._sire_object) + .toSystem()[0] + .property("intrascale") + ) + + # Set the "intrascale" properties. + intrascale_prop = property_map.get("intrascale", "intrascale") + mol.setProperty(intrascale_prop + "0", intra0) + mol.setProperty(intrascale_prop + "1", intra0) + # Commit the changes. mol = _Molecule(mol.commit()) diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py index 05ca21681..60f43224d 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py @@ -340,15 +340,10 @@ def extract(self, indices, renumber=False, property_map={}): # Append to the updated atom index. indices_.append(_SireMol.AtomIdx(x)) - if renumber: - mol = self.copy() - else: - mol = self - # Extract a partial molecule. try: # Create an empty atom selection for this molecule. - selection = mol._sire_object.selection() + selection = self._sire_object.selection() selection.selectNone() # Add the atom indices to the selection. @@ -356,7 +351,7 @@ def extract(self, indices, renumber=False, property_map={}): selection.select(idx) partial_mol = ( - _SireMol.PartialMolecule(mol._sire_object, selection) + _SireMol.PartialMolecule(self._sire_object, selection) .extract() .molecule() ) @@ -371,7 +366,7 @@ def extract(self, indices, renumber=False, property_map={}): intrascale = property_map.get("intrascale", "intrascale") # Flag whether the molecule has an intrascale property. - has_intrascale = mol._sire_object.hasProperty(intrascale) + has_intrascale = self._sire_object.hasProperty(intrascale) # Remove the "intrascale" property, since this doesn't correspond to the # extracted molecule. @@ -409,6 +404,15 @@ def extract(self, indices, renumber=False, property_map={}): else: mol = Molecule(partial_mol) + # Keep the same MolNum. + if not renumber: + mol._sire_object = ( + mol._sire_object.edit() + .renumber(self._sire_object.number()) + .commit() + .molecule() + ) + return mol def molecule0(self): diff --git a/python/BioSimSpace/_SireWrappers/_molecule.py b/python/BioSimSpace/_SireWrappers/_molecule.py index d619f1d48..424d6196d 100644 --- a/python/BioSimSpace/_SireWrappers/_molecule.py +++ b/python/BioSimSpace/_SireWrappers/_molecule.py @@ -340,15 +340,10 @@ def extract(self, indices, renumber=False, property_map={}): # Append to the updated atom index. indices_.append(_SireMol.AtomIdx(x)) - if renumber: - mol = self.copy() - else: - mol = self - # Extract a partial molecule. try: # Create an empty atom selection for this molecule. - selection = mol._sire_object.selection() + selection = self._sire_object.selection() selection.selectNone() # Add the atom indices to the selection. @@ -356,7 +351,7 @@ def extract(self, indices, renumber=False, property_map={}): selection.select(idx) partial_mol = ( - _SireMol.PartialMolecule(mol._sire_object, selection) + _SireMol.PartialMolecule(self._sire_object, selection) .extract() .molecule() ) @@ -371,7 +366,7 @@ def extract(self, indices, renumber=False, property_map={}): intrascale = property_map.get("intrascale", "intrascale") # Flag whether the molecule has an intrascale property. - has_intrascale = mol._sire_object.hasProperty(intrascale) + has_intrascale = self._sire_object.hasProperty(intrascale) # Remove the "intrascale" property, since this doesn't correspond to the # extracted molecule. @@ -409,6 +404,15 @@ def extract(self, indices, renumber=False, property_map={}): else: mol = Molecule(partial_mol) + # Keep the same MolNum. + if not renumber: + mol._sire_object = ( + mol._sire_object.edit() + .renumber(self._sire_object.number()) + .commit() + .molecule() + ) + return mol def molecule0(self): diff --git a/tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py b/tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py index d19b385f8..7f6bc0aa3 100644 --- a/tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py +++ b/tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py @@ -102,3 +102,27 @@ def test_hydrogen_mass_repartitioning(system, ignore_waters): # Assert the the masses are approximately the same. assert final_mass == pytest.approx(initial_mass) + + +def test_extract(system): + """Test the extract method.""" + + # A list of atom indices to extract. + idxs = [0, 1, 2, 3] + + # Extract the first molecule from the system. + mol = system[0] + + # Extract the atoms. + partial_mol = mol.extract(idxs) + + assert partial_mol.nAtoms() == 4 + + # Make sure the numbers match. + assert partial_mol.number() == mol.number() + + # Extract and renumber. + partial_mol = mol.extract(idxs, renumber=True) + + # Make sure the numbers are different. + assert partial_mol.number() != mol.number() diff --git a/tests/_SireWrappers/test_molecule.py b/tests/_SireWrappers/test_molecule.py index 867bac2c6..28b6c28df 100644 --- a/tests/_SireWrappers/test_molecule.py +++ b/tests/_SireWrappers/test_molecule.py @@ -96,3 +96,27 @@ def test_hydrogen_mass_repartitioning(system, ignore_waters): # Assert the the masses are approximately the same. assert final_mass == pytest.approx(initial_mass) + + +def test_extract(system): + """Test the extract method.""" + + # A list of atom indices to extract. + idxs = [0, 1, 2, 3] + + # Extract the first molecule from the system. + mol = system[0] + + # Extract the atoms. + partial_mol = mol.extract(idxs) + + assert partial_mol.nAtoms() == 4 + + # Make sure the numbers match. + assert partial_mol.number() == mol.number() + + # Extract and renumber. + partial_mol = mol.extract(idxs, renumber=True) + + # Make sure the numbers are different. + assert partial_mol.number() != mol.number() From 8bd52c91fe9f60d09d0c1329d7402075209e2366 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 4 Dec 2023 11:31:04 +0000 Subject: [PATCH 11/20] Fix typo. [ci skip] --- python/BioSimSpace/_SireWrappers/_system.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py index 0a5bf5157..cc41b3020 100644 --- a/python/BioSimSpace/_SireWrappers/_system.py +++ b/python/BioSimSpace/_SireWrappers/_system.py @@ -662,7 +662,7 @@ def addMolecules(self, molecules): _warnings.warn( "Failed to remove 'velocity' property from all molecules!" ) - if has_perturnable: + if has_perturbable: try: self._sire_object = _SireIO.removeProperty( self._sire_object, "velocity0" From 838077f3b3db400dc0198c2bf6e43b1adc66ef5b Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 11 Dec 2023 09:47:46 +0000 Subject: [PATCH 12/20] Backport fix from PR #217. [ci skip] --- python/BioSimSpace/Metadynamics/_aux/metadynamics.py | 6 ++---- python/BioSimSpace/Process/_openmm.py | 12 ++++++------ .../Exscientia/Metadynamics/_aux/metadynamics.py | 6 ++---- .../Sandpit/Exscientia/Process/_openmm.py | 12 ++++++------ 4 files changed, 16 insertions(+), 20 deletions(-) diff --git a/python/BioSimSpace/Metadynamics/_aux/metadynamics.py b/python/BioSimSpace/Metadynamics/_aux/metadynamics.py index 83fe58c91..c6c37be1a 100644 --- a/python/BioSimSpace/Metadynamics/_aux/metadynamics.py +++ b/python/BioSimSpace/Metadynamics/_aux/metadynamics.py @@ -28,10 +28,8 @@ USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from ..._Utils import _try_import - -mm = _try_import("openmm") -unit = _try_import("openmm.unit", "conda install -c conda-forge openmm") +import openmm as mm +from openmm import unit from collections import namedtuple from functools import reduce diff --git a/python/BioSimSpace/Process/_openmm.py b/python/BioSimSpace/Process/_openmm.py index 8a913b667..92a40be04 100644 --- a/python/BioSimSpace/Process/_openmm.py +++ b/python/BioSimSpace/Process/_openmm.py @@ -877,13 +877,13 @@ def _generate_config(self): self.addToConfig("\n# Upper wall.") self.addToConfig( - "upper_wall_rest = CustomCentroidBondForce(3, '(k/2)*max(distance(g1,g2)*cos(angle(g1,g2,g3)) - upper_wall, 0)^2')" + "upper_wall_rest = CustomCentroidBondForce(3, '(k1/2)*max(distance(g1,g2)*cos(angle(g1,g2,g3)) - upper_wall, 0)^2')" ) self.addToConfig("upper_wall_rest.addGroup(lig)") self.addToConfig("upper_wall_rest.addGroup(p1)") self.addToConfig("upper_wall_rest.addGroup(p2)") self.addToConfig("upper_wall_rest.addBond([0,1,2])") - self.addToConfig("upper_wall_rest.addGlobalParameter('k', k1)") + self.addToConfig("upper_wall_rest.addGlobalParameter('k1', k1)") self.addToConfig( "upper_wall_rest.addGlobalParameter('upper_wall', upper_wall)" ) @@ -903,13 +903,13 @@ def _generate_config(self): ) self.addToConfig( - "dist_restraint = CustomCentroidBondForce(3, '(k/2)*max(distance(g1,g2)*sin(angle(g1,g2,g3)) - (a/(1+exp(b*(distance(g1,g2)*cos(angle(g1,g2,g3))-c)))+d), 0)^2')" + "dist_restraint = CustomCentroidBondForce(3, '(k2/2)*max(distance(g1,g2)*sin(angle(g1,g2,g3)) - (a/(1+exp(b*(distance(g1,g2)*cos(angle(g1,g2,g3))-c)))+d), 0)^2')" ) self.addToConfig("dist_restraint.addGroup(lig)") self.addToConfig("dist_restraint.addGroup(p1)") self.addToConfig("dist_restraint.addGroup(p2)") self.addToConfig("dist_restraint.addBond([0,1,2])") - self.addToConfig("dist_restraint.addGlobalParameter('k', k2)") + self.addToConfig("dist_restraint.addGlobalParameter('k2', k2)") self.addToConfig("dist_restraint.addGlobalParameter('a', wall_width)") self.addToConfig("dist_restraint.addGlobalParameter('b', beta_cent)") self.addToConfig("dist_restraint.addGlobalParameter('c', s_cent)") @@ -919,13 +919,13 @@ def _generate_config(self): self.addToConfig("\n# Lower wall.") self.addToConfig( - "lower_wall_rest = CustomCentroidBondForce(3, '(k/2)*min(distance(g1,g2)*cos(angle(g1,g2,g3)) - lower_wall, 0)^2')" + "lower_wall_rest = CustomCentroidBondForce(3, '(k1/2)*min(distance(g1,g2)*cos(angle(g1,g2,g3)) - lower_wall, 0)^2')" ) self.addToConfig("lower_wall_rest.addGroup(lig)") self.addToConfig("lower_wall_rest.addGroup(p1)") self.addToConfig("lower_wall_rest.addGroup(p2)") self.addToConfig("lower_wall_rest.addBond([0,1,2])") - self.addToConfig("lower_wall_rest.addGlobalParameter('k', k1)") + self.addToConfig("lower_wall_rest.addGlobalParameter('k1', k1)") self.addToConfig( "lower_wall_rest.addGlobalParameter('lower_wall', lower_wall)" ) diff --git a/python/BioSimSpace/Sandpit/Exscientia/Metadynamics/_aux/metadynamics.py b/python/BioSimSpace/Sandpit/Exscientia/Metadynamics/_aux/metadynamics.py index 83fe58c91..c6c37be1a 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Metadynamics/_aux/metadynamics.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Metadynamics/_aux/metadynamics.py @@ -28,10 +28,8 @@ USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from ..._Utils import _try_import - -mm = _try_import("openmm") -unit = _try_import("openmm.unit", "conda install -c conda-forge openmm") +import openmm as mm +from openmm import unit from collections import namedtuple from functools import reduce diff --git a/python/BioSimSpace/Sandpit/Exscientia/Process/_openmm.py b/python/BioSimSpace/Sandpit/Exscientia/Process/_openmm.py index 374e6700e..07dc8c7a2 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Process/_openmm.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Process/_openmm.py @@ -878,13 +878,13 @@ def _generate_config(self): self.addToConfig("\n# Upper wall.") self.addToConfig( - "upper_wall_rest = CustomCentroidBondForce(3, '(k/2)*max(distance(g1,g2)*cos(angle(g1,g2,g3)) - upper_wall, 0)^2')" + "upper_wall_rest = CustomCentroidBondForce(3, '(k1/2)*max(distance(g1,g2)*cos(angle(g1,g2,g3)) - upper_wall, 0)^2')" ) self.addToConfig("upper_wall_rest.addGroup(lig)") self.addToConfig("upper_wall_rest.addGroup(p1)") self.addToConfig("upper_wall_rest.addGroup(p2)") self.addToConfig("upper_wall_rest.addBond([0,1,2])") - self.addToConfig("upper_wall_rest.addGlobalParameter('k', k1)") + self.addToConfig("upper_wall_rest.addGlobalParameter('k1', k1)") self.addToConfig( "upper_wall_rest.addGlobalParameter('upper_wall', upper_wall)" ) @@ -904,13 +904,13 @@ def _generate_config(self): ) self.addToConfig( - "dist_restraint = CustomCentroidBondForce(3, '(k/2)*max(distance(g1,g2)*sin(angle(g1,g2,g3)) - (a/(1+exp(b*(distance(g1,g2)*cos(angle(g1,g2,g3))-c)))+d), 0)^2')" + "dist_restraint = CustomCentroidBondForce(3, '(k2/2)*max(distance(g1,g2)*sin(angle(g1,g2,g3)) - (a/(1+exp(b*(distance(g1,g2)*cos(angle(g1,g2,g3))-c)))+d), 0)^2')" ) self.addToConfig("dist_restraint.addGroup(lig)") self.addToConfig("dist_restraint.addGroup(p1)") self.addToConfig("dist_restraint.addGroup(p2)") self.addToConfig("dist_restraint.addBond([0,1,2])") - self.addToConfig("dist_restraint.addGlobalParameter('k', k2)") + self.addToConfig("dist_restraint.addGlobalParameter('k2', k2)") self.addToConfig("dist_restraint.addGlobalParameter('a', wall_width)") self.addToConfig("dist_restraint.addGlobalParameter('b', beta_cent)") self.addToConfig("dist_restraint.addGlobalParameter('c', s_cent)") @@ -920,13 +920,13 @@ def _generate_config(self): self.addToConfig("\n# Lower wall.") self.addToConfig( - "lower_wall_rest = CustomCentroidBondForce(3, '(k/2)*min(distance(g1,g2)*cos(angle(g1,g2,g3)) - lower_wall, 0)^2')" + "lower_wall_rest = CustomCentroidBondForce(3, '(k1/2)*min(distance(g1,g2)*cos(angle(g1,g2,g3)) - lower_wall, 0)^2')" ) self.addToConfig("lower_wall_rest.addGroup(lig)") self.addToConfig("lower_wall_rest.addGroup(p1)") self.addToConfig("lower_wall_rest.addGroup(p2)") self.addToConfig("lower_wall_rest.addBond([0,1,2])") - self.addToConfig("lower_wall_rest.addGlobalParameter('k', k1)") + self.addToConfig("lower_wall_rest.addGlobalParameter('k1', k1)") self.addToConfig( "lower_wall_rest.addGlobalParameter('lower_wall', lower_wall)" ) From 34d3d9b726499a9b52c690e206d51a3609472af4 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 11 Dec 2023 09:49:04 +0000 Subject: [PATCH 13/20] Backport fix from PR #219. [ci skip] --- .../BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py | 4 ++-- python/BioSimSpace/Trajectory/_trajectory.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py b/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py index b49a05d7d..1dc7e7d24 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Trajectory/_trajectory.py @@ -921,9 +921,9 @@ def nFrames(self): The number of trajectory frames. """ - # First get the current MDTraj object. + # First get the current trajectory object using the existing backend. if self._process is not None and self._process.isRunning(): - self._trajectory = self.getTrajectory() + self._trajectory = self.getTrajectory(format=self._backend) # There is no trajectory. if self._trajectory is None: diff --git a/python/BioSimSpace/Trajectory/_trajectory.py b/python/BioSimSpace/Trajectory/_trajectory.py index b49a05d7d..1dc7e7d24 100644 --- a/python/BioSimSpace/Trajectory/_trajectory.py +++ b/python/BioSimSpace/Trajectory/_trajectory.py @@ -921,9 +921,9 @@ def nFrames(self): The number of trajectory frames. """ - # First get the current MDTraj object. + # First get the current trajectory object using the existing backend. if self._process is not None and self._process.isRunning(): - self._trajectory = self.getTrajectory() + self._trajectory = self.getTrajectory(format=self._backend) # There is no trajectory. if self._trajectory is None: From 565e0d913915ae7a1f0c15bb5e959d05e9ad77ea Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 11 Dec 2023 09:51:46 +0000 Subject: [PATCH 14/20] Backport fix from PR #220. [ci skip] --- python/BioSimSpace/_Config/_amber.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/BioSimSpace/_Config/_amber.py b/python/BioSimSpace/_Config/_amber.py index 744de459f..55537109c 100644 --- a/python/BioSimSpace/_Config/_amber.py +++ b/python/BioSimSpace/_Config/_amber.py @@ -178,7 +178,7 @@ def createConfig( protocol_dict["cut"] = "999." if is_pmemd: # Use vacuum generalised Born model. - self.addToConfig(" igb=6,") + protocol_dict["igb"] = "6" else: # Non-bonded cut-off. protocol_dict["cut"] = "8.0" From 322062ad9418a0e376df2aac477ae4bd6e1a196d Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 11 Dec 2023 10:06:29 +0000 Subject: [PATCH 15/20] Backport fix from PR #222. [ci skip] --- python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py | 2 +- python/BioSimSpace/_Config/_amber.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py b/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py index e71e1369d..20ac4e3fe 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Protocol/_config.py @@ -278,7 +278,7 @@ def generateAmberConfig(self, extra_options=None, extra_lines=None): if restraint == "backbone": restraint_mask = "@CA,C,O,N" elif restraint == "heavy": - restraint_mask = "!:WAT & !@H" + restraint_mask = "!:WAT & !@H=" elif restraint == "all": restraint_mask = "!:WAT" diff --git a/python/BioSimSpace/_Config/_amber.py b/python/BioSimSpace/_Config/_amber.py index 55537109c..5941bc51a 100644 --- a/python/BioSimSpace/_Config/_amber.py +++ b/python/BioSimSpace/_Config/_amber.py @@ -212,7 +212,7 @@ def createConfig( if restraint == "backbone": restraint_mask = "@CA,C,O,N" elif restraint == "heavy": - restraint_mask = "!:WAT & !@H" + restraint_mask = "!:WAT & !@H=" elif restraint == "all": restraint_mask = "!:WAT" From 8311ce247fb2ae423114f0b115b954ef3677bb28 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 12 Dec 2023 11:13:46 +0000 Subject: [PATCH 16/20] Backport fix from PR #224. [ci skip] --- python/BioSimSpace/Parameters/_Protocol/_amber.py | 2 +- .../Sandpit/Exscientia/Parameters/_Protocol/_amber.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/python/BioSimSpace/Parameters/_Protocol/_amber.py b/python/BioSimSpace/Parameters/_Protocol/_amber.py index d83d016aa..d5ef645d3 100644 --- a/python/BioSimSpace/Parameters/_Protocol/_amber.py +++ b/python/BioSimSpace/Parameters/_Protocol/_amber.py @@ -725,7 +725,7 @@ def _get_disulphide_bonds( # Try searching for disulphide bonds. try: - disulphides = query(mol, property_map) + disulphides = query(mol, property_map).bonds() except: disulphides = [] diff --git a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py index d83d016aa..d5ef645d3 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py @@ -725,7 +725,7 @@ def _get_disulphide_bonds( # Try searching for disulphide bonds. try: - disulphides = query(mol, property_map) + disulphides = query(mol, property_map).bonds() except: disulphides = [] From 38e48607eddf9806bbfe034c51d68e591cccffbe Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 14 Dec 2023 09:06:13 +0000 Subject: [PATCH 17/20] Update Sire version. --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 39fffe478..b501e7840 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,7 +1,7 @@ # BioSimSpace runtime requirements. # main -sire~=2023.4.0 +sire~=2023.4.2 # devel #sire==2023.5.0.dev From 89c48571d4f95fcfaaf8215090968e8cfe88d35c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 14 Dec 2023 12:11:20 +0000 Subject: [PATCH 18/20] Backport fix from PR #226. --- .../Parameters/_Protocol/_amber.py | 76 +++++++-------- python/BioSimSpace/Parameters/_parameters.py | 97 +++++++++++-------- .../Exscientia/Parameters/_Protocol/_amber.py | 76 +++++++-------- .../Exscientia/Parameters/_parameters.py | 97 +++++++++++-------- tests/Parameters/test_parameters.py | 60 ++++++++++++ .../Exscientia/Parameters/test_parameters.py | 60 ++++++++++++ 6 files changed, 298 insertions(+), 168 deletions(-) diff --git a/python/BioSimSpace/Parameters/_Protocol/_amber.py b/python/BioSimSpace/Parameters/_Protocol/_amber.py index d5ef645d3..5b1601686 100644 --- a/python/BioSimSpace/Parameters/_Protocol/_amber.py +++ b/python/BioSimSpace/Parameters/_Protocol/_amber.py @@ -122,8 +122,8 @@ def __init__( tolerance=1.2, max_distance=_Length(6, "A"), water_model=None, - custom_parameters=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, property_map={}, @@ -155,14 +155,17 @@ def __init__( Run 'BioSimSpace.Solvent.waterModels()' to see the supported water models. - custom_parameters: [str] - A list of paths to custom parameter files. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -218,34 +221,25 @@ def __init__( else: self._water_model = water_model - # Validate the custom parameter file list. - if custom_parameters is not None: - if not isinstance(custom_parameters, (list, tuple)): - raise TypeError("'custom_parameters' must be a 'list' of 'str' types.") + # Validate the additional leap commands. + if pre_mol_commands is not None: + if not isinstance(pre_mol_commands, (list, tuple)): + raise TypeError("'pre_mol_commands must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in custom_parameters): + if not all(isinstance(x, str) for x in pre_mol_commands): raise TypeError( - "'custom_parameters' must be a 'list' of 'str' types." + "'pre_mol_commands' must be a 'list' of 'str' types." ) - for x in custom_parameters: - if not os.path.isfile(x): - raise ValueError(f"Custom parameter file does not exist: '{x}'") - - # Convert to absolute paths. - self._custom_parameters = [] - for x in enumerate(custom_parameters): - self._custom_parameters.append(_os.path.abspath(x)) - else: - self._custom_parameters = None - - # Validate the additional leap commands. - if leap_commands is not None: - if not isinstance(leap_commands, (list, tuple)): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") + self._pre_mol_commands = pre_mol_commands + if post_mol_commands is not None: + if not isinstance(post_mol_commands, (list, tuple)): + raise TypeError("'post_mol_commands must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in leap_commands): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") - self._leap_commands = leap_commands + if not all(isinstance(x, str) for x in post_mol_commands): + raise TypeError( + "'post_mol_commands' must be a 'list' of 'str' types." + ) + self._post_mol_commands = post_mol_commands # Validate the bond records. if bonds is not None: @@ -273,7 +267,7 @@ def __init__( # Set the compatibility flags. self._tleap = True - if self._custom_parameters is not None or self._leap_commands is not None: + if self._pre_mol_commands is not None or self._post_mol_commands is not None: self._pdb2gmx = False def run(self, molecule, work_dir=None, queue=None): @@ -512,10 +506,10 @@ def _run_tleap(self, molecule, work_dir): file.write("source leaprc.water.tip4pew\n") else: file.write("source leaprc.water.%s\n" % self._water_model) - # Write custom parameters. - if self._custom_parameters is not None: - for param in self._custom_parameters: - file.write("%s\n" % param) + # Write pre-mol user leap commands. + if self._pre_mol_commands is not None: + for command in self._pre_mol_commands: + file.write("%s\n" % command) file.write("mol = loadPdb leap.pdb\n") # Add any disulphide bond records. for bond in disulphide_bonds: @@ -523,9 +517,9 @@ def _run_tleap(self, molecule, work_dir): # Add any additional bond records. for bond in pruned_bond_records: file.write("%s\n" % bond) - # Write user leap commands. - if self._leap_commands is not None: - for command in self._leap_commands: + # Write post-mol user leap commands. + if self._post_mol_commands is not None: + for command in self._post_mol_commands: file.write("%s\n" % command) file.write("saveAmberParm mol leap.top leap.crd\n") file.write("quit") diff --git a/python/BioSimSpace/Parameters/_parameters.py b/python/BioSimSpace/Parameters/_parameters.py index f62d5630e..945319d99 100644 --- a/python/BioSimSpace/Parameters/_parameters.py +++ b/python/BioSimSpace/Parameters/_parameters.py @@ -135,8 +135,8 @@ def _parameterise_amber_protein( tolerance=1.2, max_distance=_Length(6, "A"), water_model=None, - custom_parameters=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, work_dir=None, @@ -172,14 +172,17 @@ def _parameterise_amber_protein( or lone oxygen atoms corresponding to structural (crystal) water molecules. - custom_parameters: [str] - A list of paths to custom parameter files. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -241,8 +244,8 @@ def _parameterise_amber_protein( max_distance=max_distance, water_model=water_model, check_ions=True, - custom_parameters=custom_parameters, - leap_commands=leap_commands, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, bonds=bonds, ensure_compatible=ensure_compatible, property_map=property_map, @@ -255,8 +258,8 @@ def _parameterise_amber_protein( tolerance=tolerance, water_model=water_model, max_distance=max_distance, - custom_parameters=custom_parameters, - leap_commands=leap_commands, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, bonds=bonds, ensure_compatible=ensure_compatible, property_map=property_map, @@ -802,8 +805,8 @@ def _validate( max_distance=None, water_model=None, check_ions=False, - custom_parameters=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, work_dir=None, @@ -839,14 +842,17 @@ def _validate( Whether to check for the presence of structural ions. This is only required when parameterising with protein force fields. - custom_parameters: [str] - A list of paths to custom parameter files. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -905,22 +911,19 @@ def _validate( "Please choose a 'water_model' for the ion parameters." ) - if custom_parameters is not None: - if not isinstance(custom_parameters, (list, tuple)): - raise TypeError("'custom_parameters' must be a 'list' of 'str' types.") + if pre_mol_commands is not None: + if not isinstance(pre_mol_commands, (list, tuple)): + raise TypeError("'pre_mol_commands' must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in custom_parameters): - raise TypeError("'custom_parameters' must be a 'list' of 'str' types.") - for x in custom_parameters: - if not os.path.isfile(x): - raise ValueError(f"Custom parameter file does not exist: '{x}'") - - if leap_commands is not None: - if not isinstance(leap_commands, (list, tuple)): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") + if not all(isinstance(x, str) for x in pre_mol_commands): + raise TypeError("'pre_mol_commands' must be a 'list' of 'str' types.") + + if post_mol_commands is not None: + if not isinstance(post_mol_commands, (list, tuple)): + raise TypeError("'post_mol_commands' must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in leap_commands): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") + if not all(isinstance(x, str) for x in post_mol_commands): + raise TypeError("'post_mol_commands' must be a 'list' of 'str' types.") if bonds is not None: if molecule is None: @@ -981,7 +984,8 @@ def _function( tolerance=1.2, max_distance=_Length(6, "A"), water_model=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, work_dir=None, @@ -1013,11 +1017,17 @@ def _function( or lone oxygen atoms corresponding to structural (crystal) water molecules. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands and can be used to, e.g., - load additional parameter files. When this option is set, we can no - longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. + + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -1052,7 +1062,8 @@ def _function( tolerance=tolerance, max_distance=max_distance, water_model=water_model, - leap_commands=leap_commands, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, bonds=bonds, ensure_compatible=ensure_compatible, work_dir=work_dir, diff --git a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py index d5ef645d3..5b1601686 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_Protocol/_amber.py @@ -122,8 +122,8 @@ def __init__( tolerance=1.2, max_distance=_Length(6, "A"), water_model=None, - custom_parameters=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, property_map={}, @@ -155,14 +155,17 @@ def __init__( Run 'BioSimSpace.Solvent.waterModels()' to see the supported water models. - custom_parameters: [str] - A list of paths to custom parameter files. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -218,34 +221,25 @@ def __init__( else: self._water_model = water_model - # Validate the custom parameter file list. - if custom_parameters is not None: - if not isinstance(custom_parameters, (list, tuple)): - raise TypeError("'custom_parameters' must be a 'list' of 'str' types.") + # Validate the additional leap commands. + if pre_mol_commands is not None: + if not isinstance(pre_mol_commands, (list, tuple)): + raise TypeError("'pre_mol_commands must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in custom_parameters): + if not all(isinstance(x, str) for x in pre_mol_commands): raise TypeError( - "'custom_parameters' must be a 'list' of 'str' types." + "'pre_mol_commands' must be a 'list' of 'str' types." ) - for x in custom_parameters: - if not os.path.isfile(x): - raise ValueError(f"Custom parameter file does not exist: '{x}'") - - # Convert to absolute paths. - self._custom_parameters = [] - for x in enumerate(custom_parameters): - self._custom_parameters.append(_os.path.abspath(x)) - else: - self._custom_parameters = None - - # Validate the additional leap commands. - if leap_commands is not None: - if not isinstance(leap_commands, (list, tuple)): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") + self._pre_mol_commands = pre_mol_commands + if post_mol_commands is not None: + if not isinstance(post_mol_commands, (list, tuple)): + raise TypeError("'post_mol_commands must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in leap_commands): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") - self._leap_commands = leap_commands + if not all(isinstance(x, str) for x in post_mol_commands): + raise TypeError( + "'post_mol_commands' must be a 'list' of 'str' types." + ) + self._post_mol_commands = post_mol_commands # Validate the bond records. if bonds is not None: @@ -273,7 +267,7 @@ def __init__( # Set the compatibility flags. self._tleap = True - if self._custom_parameters is not None or self._leap_commands is not None: + if self._pre_mol_commands is not None or self._post_mol_commands is not None: self._pdb2gmx = False def run(self, molecule, work_dir=None, queue=None): @@ -512,10 +506,10 @@ def _run_tleap(self, molecule, work_dir): file.write("source leaprc.water.tip4pew\n") else: file.write("source leaprc.water.%s\n" % self._water_model) - # Write custom parameters. - if self._custom_parameters is not None: - for param in self._custom_parameters: - file.write("%s\n" % param) + # Write pre-mol user leap commands. + if self._pre_mol_commands is not None: + for command in self._pre_mol_commands: + file.write("%s\n" % command) file.write("mol = loadPdb leap.pdb\n") # Add any disulphide bond records. for bond in disulphide_bonds: @@ -523,9 +517,9 @@ def _run_tleap(self, molecule, work_dir): # Add any additional bond records. for bond in pruned_bond_records: file.write("%s\n" % bond) - # Write user leap commands. - if self._leap_commands is not None: - for command in self._leap_commands: + # Write post-mol user leap commands. + if self._post_mol_commands is not None: + for command in self._post_mol_commands: file.write("%s\n" % command) file.write("saveAmberParm mol leap.top leap.crd\n") file.write("quit") diff --git a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_parameters.py b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_parameters.py index f62d5630e..945319d99 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/Parameters/_parameters.py +++ b/python/BioSimSpace/Sandpit/Exscientia/Parameters/_parameters.py @@ -135,8 +135,8 @@ def _parameterise_amber_protein( tolerance=1.2, max_distance=_Length(6, "A"), water_model=None, - custom_parameters=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, work_dir=None, @@ -172,14 +172,17 @@ def _parameterise_amber_protein( or lone oxygen atoms corresponding to structural (crystal) water molecules. - custom_parameters: [str] - A list of paths to custom parameter files. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -241,8 +244,8 @@ def _parameterise_amber_protein( max_distance=max_distance, water_model=water_model, check_ions=True, - custom_parameters=custom_parameters, - leap_commands=leap_commands, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, bonds=bonds, ensure_compatible=ensure_compatible, property_map=property_map, @@ -255,8 +258,8 @@ def _parameterise_amber_protein( tolerance=tolerance, water_model=water_model, max_distance=max_distance, - custom_parameters=custom_parameters, - leap_commands=leap_commands, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, bonds=bonds, ensure_compatible=ensure_compatible, property_map=property_map, @@ -802,8 +805,8 @@ def _validate( max_distance=None, water_model=None, check_ions=False, - custom_parameters=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, work_dir=None, @@ -839,14 +842,17 @@ def _validate( Whether to check for the presence of structural ions. This is only required when parameterising with protein force fields. - custom_parameters: [str] - A list of paths to custom parameter files. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands. When this option is set, - we can no longer fall back on GROMACS's pdb2gmx. + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -905,22 +911,19 @@ def _validate( "Please choose a 'water_model' for the ion parameters." ) - if custom_parameters is not None: - if not isinstance(custom_parameters, (list, tuple)): - raise TypeError("'custom_parameters' must be a 'list' of 'str' types.") + if pre_mol_commands is not None: + if not isinstance(pre_mol_commands, (list, tuple)): + raise TypeError("'pre_mol_commands' must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in custom_parameters): - raise TypeError("'custom_parameters' must be a 'list' of 'str' types.") - for x in custom_parameters: - if not os.path.isfile(x): - raise ValueError(f"Custom parameter file does not exist: '{x}'") - - if leap_commands is not None: - if not isinstance(leap_commands, (list, tuple)): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") + if not all(isinstance(x, str) for x in pre_mol_commands): + raise TypeError("'pre_mol_commands' must be a 'list' of 'str' types.") + + if post_mol_commands is not None: + if not isinstance(post_mol_commands, (list, tuple)): + raise TypeError("'post_mol_commands' must be a 'list' of 'str' types.") else: - if not all(isinstance(x, str) for x in leap_commands): - raise TypeError("'leap_commands' must be a 'list' of 'str' types.") + if not all(isinstance(x, str) for x in post_mol_commands): + raise TypeError("'post_mol_commands' must be a 'list' of 'str' types.") if bonds is not None: if molecule is None: @@ -981,7 +984,8 @@ def _function( tolerance=1.2, max_distance=_Length(6, "A"), water_model=None, - leap_commands=None, + pre_mol_commands=None, + post_mol_commands=None, bonds=None, ensure_compatible=True, work_dir=None, @@ -1013,11 +1017,17 @@ def _function( or lone oxygen atoms corresponding to structural (crystal) water molecules. - leap_commands : [str] - An optional list of extra commands for the LEaP program. These - will be added after any default commands and can be used to, e.g., - load additional parameter files. When this option is set, we can no - longer fall back on GROMACS's pdb2gmx. + pre_mol_commands : [str] + A list of custom LEaP commands to be executed before loading the molecule. + This can be used for loading custom parameter files, or sourcing additional + scripts. Make sure to use absolute paths when specifying any external files. + When this option is set, we can no longer fall back on GROMACS's pdb2gmx. + + post_mol_commands : [str] + A list of custom LEaP commands to be executed after loading the molecule. + This allows the use of additional commands that operate on the molecule, + which is named "mol". When this option is set, we can no longer fall + back on GROMACS's pdb2gmx. bonds : ((class:`Atom `, class:`Atom `)) An optional tuple of atom pairs to specify additional atoms that @@ -1052,7 +1062,8 @@ def _function( tolerance=tolerance, max_distance=max_distance, water_model=water_model, - leap_commands=leap_commands, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, bonds=bonds, ensure_compatible=ensure_compatible, work_dir=work_dir, diff --git a/tests/Parameters/test_parameters.py b/tests/Parameters/test_parameters.py index 55ab15afb..8acd41ff7 100644 --- a/tests/Parameters/test_parameters.py +++ b/tests/Parameters/test_parameters.py @@ -1,4 +1,5 @@ import pytest +import tempfile import BioSimSpace as BSS @@ -77,3 +78,62 @@ def test_molecule_rename(): # Make sure the name is correct. assert mol._sire_object.name().value() == "smiles:[C@@H](C(F)(F)F)(OC(F)F)Cl" + + +@pytest.mark.skipif( + has_antechamber is False or has_tleap is False, + reason="Requires AmberTools/antechamber and tLEaP to be installed.", +) +def test_leap_commands(molecule0): + """ + Test that custom commands are correctly inserted into the LEaP input + script. + """ + + # Create lists of pre- and post-commands. + pre_mol_commands = ["command1", "command2"] + post_mol_commands = ["command3", "command4"] + + with tempfile.TemporaryDirectory() as tmp: + # This will fail, but we only want to check the LEaP script. + try: + mol = BSS.Parameters.ff14SB( + molecule0, + work_dir=tmp, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, + ).getMolecule() + except: + pass + + # Load the script and check that the custom parameters are present and + # in the correct order. + with open(f"{tmp}/leap.txt", "r") as f: + script = f.readlines() + + # Create lists to store the indices of the custom commands. + line_pre = [-1 for _ in range(len(pre_mol_commands))] + line_post = [-1 for _ in range(len(post_mol_commands))] + + # Loop over the lines in the script and store the line numbers + # where the custom commands are found. + for x, line in enumerate(script): + for y, command in enumerate(pre_mol_commands): + if command in line: + line_pre[y] = x + for y, command in enumerate(post_mol_commands): + if command in line: + line_post[y] = x + + # Make sure the lines are found. + for line in line_pre: + assert line != -1 + for line in line_post: + assert line != -1 + + # Make sure the lines are in the correct order. + for x in range(len(line_pre) - 1): + assert line_pre[x] < line_pre[x + 1] + assert line_pre[-1] < line_post[0] + for x in range(len(line_post) - 1): + assert line_post[x] < line_post[x + 1] diff --git a/tests/Sandpit/Exscientia/Parameters/test_parameters.py b/tests/Sandpit/Exscientia/Parameters/test_parameters.py index db89a8025..f42f2d539 100644 --- a/tests/Sandpit/Exscientia/Parameters/test_parameters.py +++ b/tests/Sandpit/Exscientia/Parameters/test_parameters.py @@ -1,4 +1,5 @@ import pytest +import tempfile import BioSimSpace.Sandpit.Exscientia as BSS @@ -82,3 +83,62 @@ def test_molecule_rename(): # Make sure the name is correct. assert mol._sire_object.name().value() == "smiles:[C@@H](C(F)(F)F)(OC(F)F)Cl" + + +@pytest.mark.skipif( + has_antechamber is False or has_tleap is False, + reason="Requires AmberTools/antechamber and tLEaP to be installed.", +) +def test_leap_commands(molecule0): + """ + Test that custom commands are correctly inserted into the LEaP input + script. + """ + + # Create lists of pre- and post-commands. + pre_mol_commands = ["command1", "command2"] + post_mol_commands = ["command3", "command4"] + + with tempfile.TemporaryDirectory() as tmp: + # This will fail, but we only want to check the LEaP script. + try: + mol = BSS.Parameters.ff14SB( + molecule0, + work_dir=tmp, + pre_mol_commands=pre_mol_commands, + post_mol_commands=post_mol_commands, + ).getMolecule() + except: + pass + + # Load the script and check that the custom parameters are present and + # in the correct order. + with open(f"{tmp}/leap.txt", "r") as f: + script = f.readlines() + + # Create lists to store the indices of the custom commands. + line_pre = [-1 for _ in range(len(pre_mol_commands))] + line_post = [-1 for _ in range(len(post_mol_commands))] + + # Loop over the lines in the script and store the line numbers + # where the custom commands are found. + for x, line in enumerate(script): + for y, command in enumerate(pre_mol_commands): + if command in line: + line_pre[y] = x + for y, command in enumerate(post_mol_commands): + if command in line: + line_post[y] = x + + # Make sure the lines are found. + for line in line_pre: + assert line != -1 + for line in line_post: + assert line != -1 + + # Make sure the lines are in the correct order. + for x in range(len(line_pre) - 1): + assert line_pre[x] < line_pre[x + 1] + assert line_pre[-1] < line_post[0] + for x in range(len(line_post) - 1): + assert line_post[x] < line_post[x + 1] From 3f97e006f5636b7e2fcd4da62422ecea5e2530cd Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 14 Dec 2023 12:27:07 +0000 Subject: [PATCH 19/20] Update CHANGELOG for the 2023.4.1 release. [ci skip] --- doc/source/changelog.rst | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index a78abcd71..77dde7a7b 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -9,6 +9,25 @@ company supporting open-source development of fostering academic/industrial coll within the biomolecular simulation community. Our software is hosted via the `OpenBioSim` `GitHub `__ organisation. +`2023.4.1 `_ - Dec 14 2023 +------------------------------------------------------------------------------------------------- + +* Make sure ``match_water`` keyword argument is passed to specialised solvation functions (`#190 `__). +* Check perturbable molecules for velocities when combining molecules (`#192 `__). +* Make sure velocities are double counted when searching for velocity properties when combining molecules (`#197 `__). +* Remove redundant ``BioSimSpace.Types.Type.__ne__`` operator (`#201 `__). +* Minor internal updates due to Sire API fixes (`#203 `__). +* Fixed bug in the BSS Boresch restraint search code (`@fjclark `_) (`#204 `__). +* Fixed ``renumber`` option in :meth:`extract ` method (`#210 `__). +* Add workaround for fixing reconstruction of intrascale matrix in :func:`readPerturbableSystem ` function (`#210 `__). +* Remove incorrect ``try_import`` statement in metadynamics driver script and make sure that global parameters in OpenMM script are unique (`#217 `__). +* Ensure the existing trajectory backend is used when getting the number of trajectory frames from a running process (`#219 `__). +* Fixed setting of ``igb`` config parameter for PMEMD simulations (`@annamherz `_) (`#220 `__). +* Make sure AMBER restraint mask matches all hydrogen atoms (`#222 `__). +* Ensure all searches for disulphide bonds are convert to a ``SelectorBond`` object (`#224 `__). +* Fix injection of custom commands into ``LEaP`` script (`#226 `__). + + `2023.4.0 `_ - Oct 13 2023 ------------------------------------------------------------------------------------------------- From 2cb653526d14e83a9a35d7c9ce1d7d1b3a3a1a35 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 15 Dec 2023 16:29:04 +0000 Subject: [PATCH 20/20] Fix typo in test name. --- tests/IO/test_file_cache.py | 2 +- tests/Sandpit/Exscientia/IO/test_file_cache.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/IO/test_file_cache.py b/tests/IO/test_file_cache.py index 25d831b2f..9d3866c99 100644 --- a/tests/IO/test_file_cache.py +++ b/tests/IO/test_file_cache.py @@ -87,7 +87,7 @@ def test_file_cache(): has_amber is False or has_openff is False, reason="Requires AMBER and OpenFF to be installed", ) -def test_file_cache_mol_nuums(): +def test_file_cache_mol_nums(): """ Make sure that systems can be cached if they have the same UID, but contain different MolNUms. diff --git a/tests/Sandpit/Exscientia/IO/test_file_cache.py b/tests/Sandpit/Exscientia/IO/test_file_cache.py index 92e0a8d91..20ed2f9f3 100644 --- a/tests/Sandpit/Exscientia/IO/test_file_cache.py +++ b/tests/Sandpit/Exscientia/IO/test_file_cache.py @@ -88,7 +88,7 @@ def test_file_cache(): has_amber is False or has_openff is False, reason="Requires AMBER and OpenFF to be installed", ) -def test_file_cache_mol_nuums(): +def test_file_cache_mol_nums(): """ Make sure that systems can be cached if they have the same UID, but contain different MolNUms.