From bd77d4ed19097fa57cfd5eee4cf9bf244c63bf0b Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 7 Sep 2023 09:09:22 +0100 Subject: [PATCH 01/11] Remove note regarding reduction rounding bias. --- doc/source/api/index_Box.rst | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/doc/source/api/index_Box.rst b/doc/source/api/index_Box.rst index 32e7e7677..892320348 100644 --- a/doc/source/api/index_Box.rst +++ b/doc/source/api/index_Box.rst @@ -11,11 +11,7 @@ To support triclinic boxes that work across the range of molecular simulation engines that support, we represent the triclinic space in reduced form, using the approach documented in Appendix A of Chapter 3 from `"Molecular dynamics of sense and sensibility in processing and analysis of data" `__ -by Tsjerk A. Wassenaar. Due to the fixed-width format used to represent box -dimensions and angles in the various molecular input files, repeated reading -and writing can lead to oscillation of the box angles on reduction due to -rounding precision errors. To account for this, we add a small bias so that -we always round in a consistent direction. +by Tsjerk A. Wassenaar. .. automodule:: BioSimSpace.Box From 09305736b6a0e402cb3d08e7bcd3c9628593baf5 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 7 Sep 2023 11:58:37 +0100 Subject: [PATCH 02/11] Update expected box angles in test. --- tests/Sandpit/Exscientia/_SireWrappers/test_system.py | 2 +- tests/_SireWrappers/test_system.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py index 86d4e0f85..f715ee05a 100644 --- a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py +++ b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py @@ -274,7 +274,7 @@ def test_set_box(system): # Store the expected box dimensions and angles. expected_box = [30, 30, 30] * BSS.Units.Length.angstrom - expected_angles = [70.5288, 109.4712, 70.5288] * BSS.Units.Angle.degree + expected_angles = [109.4712, 109.4712, 109.4712] * BSS.Units.Angle.degree # Check that the box dimensions match. for b0, b1 in zip(box, expected_box): diff --git a/tests/_SireWrappers/test_system.py b/tests/_SireWrappers/test_system.py index acd45eb39..acfe504bf 100644 --- a/tests/_SireWrappers/test_system.py +++ b/tests/_SireWrappers/test_system.py @@ -274,7 +274,7 @@ def test_set_box(system): # Store the expected box dimensions and angles. expected_box = [30, 30, 30] * BSS.Units.Length.angstrom - expected_angles = [70.5288, 109.4712, 70.5288] * BSS.Units.Angle.degree + expected_angles = [109.4712, 109.4712, 109.4712] * BSS.Units.Angle.degree # Check that the box dimensions match. for b0, b1 in zip(box, expected_box): From 7e79b8b46a98972c0d622d34d4aa2ae28f54d8ed Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 7 Sep 2023 12:00:11 +0100 Subject: [PATCH 03/11] Add methods to rotate and reduce box vectors. --- .../Exscientia/_SireWrappers/_system.py | 134 ++++++++++++++++++ python/BioSimSpace/_SireWrappers/_system.py | 134 ++++++++++++++++++ 2 files changed, 268 insertions(+) diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py index 8793472ef..680c612d6 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py @@ -1115,6 +1115,140 @@ def nMLMolecules(self): """ return len(self.getMLMolecules()) + def rotateBoxVectors(self, precision=0.0, property_map={}): + """ + Rotate the box vectors of the system so that the first vector is + aligned with the x-axis, the second vector lies in the x-y plane, + and the third vector has a positive z-component. This is a requirement + of certain molecular dynamics engines and is used for simulation + efficiency. All vector properties of the system, e.g. positions, + velocities, etc., will be rotated accordingly. + + Parameters + ---------- + + precision : float + The precision to use when sorting box vectors by their magnitude. + This can be used to prevent unecessary box rotation when the + vectors were read from a text file with limited precision. + + property_map : dict + A dictionary that maps system "properties" to their user defined + values. This allows the user to refer to properties with their + own naming scheme, e.g. { "charge" : "my-charge" } + """ + + # First check that the system has a space. + space_prop = property_map.get("space", "space") + try: + space = self._sire_object.property(space_prop) + except: + return + + # A rotation is only applicable to a triclinic space. + if not isinstance(space, _SireVol.TriclinicBox): + return + + # Rotate the space according to the desired precision. + space.rotate(precision) + + # Get the rotation matrix. + rotation_matrix = space.rotationMatrix() + + from sire.system import System + + # Create a cursor. + cursor = System(self._sire_object).cursor() + + # Rotate all vector properties. + + # Coordinates. + try: + prop_name = property_map.get("coordinates", "coordinates") + cursor = cursor.rotate( + matrix=rotation_matrix, map={"coordinates": prop_name} + ) + except: + pass + + # Velocities. + try: + prop_name = property_map.get("velocity", "velocity") + cursor = cursor.rotate( + matrix=rotation_matrix, map={"coordinates": prop_name} + ) + except: + pass + + # Now deal with any perturbable molecules. + if self.nPerturbableMolecules() > 0: + # Coordinates. + try: + prop_name = property_map.get("coordinates", "coordinates") + "0" + cursor = cursor.rotate( + matrix=rotation_matrix, map={"coordinates": prop_name} + ) + prop_name = property_map.get("coordinates", "coordinates") + "1" + cursor = cursor.rotate( + matrix=rotation_matrix, map={"coordinates": prop_name} + ) + except: + pass + + # Velocities. + try: + prop_name = property_map.get("velocity", "velocity") + "0" + cursor = cursor.rotate( + matrix=rotation_matrix, map={"coordinates": prop_name} + ) + prop_name = property_map.get("velocity", "velocity") + "1" + cursor = cursor.rotate( + matrix=rotation_matrix, map={"coordinates": prop_name} + ) + except: + pass + + # Commit the changes. + self._sire_object = cursor.commit()._system + + def reduceBoxVectors(self, bias=0, property_map={}): + """ + Reduce the box vectors. + + Parameters + ---------- + + bias : float + The bias to use when rounding during the lattice reduction. Negative + values biases towards left-tilting boxes, whereas positive values + biases towards right-tilting boxes. This can be used to ensure that + rounding is performed in a consistent direction, avoiding oscillation + when the box is instantiated from box vectors, or dimensions and + angles, that have been read from fixed-precision input files. + + property_map : dict + A dictionary that maps system "properties" to their user defined + values. This allows the user to refer to properties with their + own naming scheme, e.g. { "charge" : "my-charge" } + """ + + # First check that the system has a space. + space_prop = property_map.get("space", "space") + try: + space = self._sire_object.property(space_prop) + except: + return + + # A rotation is only applicable to a triclinic space. + if not isinstance(space, _SireVol.TriclinicBox): + return + + # Reduce the space using the desired bias. + space.reduce(bias) + + # Update the space property in the sire object. + self._sire_object.setProperty(space_prop, space) + def repartitionHydrogenMass( self, factor=4, water="no", use_coordinates=False, property_map={} ): diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py index 14ef3ba65..7d42a5c32 100644 --- a/python/BioSimSpace/_SireWrappers/_system.py +++ b/python/BioSimSpace/_SireWrappers/_system.py @@ -1063,6 +1063,140 @@ def nPerturbableMolecules(self): """ return len(self.getPerturbableMolecules()) + def rotateBoxVectors(self, precision=0.0, property_map={}): + """ + Rotate the box vectors of the system so that the first vector is + aligned with the x-axis, the second vector lies in the x-y plane, + and the third vector has a positive z-component. This is a requirement + of certain molecular dynamics engines and is used for simulation + efficiency. All vector properties of the system, e.g. positions, + velocities, etc., will be rotated accordingly. + + Parameters + ---------- + + precision : float + The precision to use when sorting box vectors by their magnitude. + This can be used to prevent unecessary box rotation when the + vectors were read from a text file with limited precision. + + property_map : dict + A dictionary that maps system "properties" to their user defined + values. This allows the user to refer to properties with their + own naming scheme, e.g. { "charge" : "my-charge" } + """ + + # First check that the system has a space. + space_prop = property_map.get("space", "space") + try: + space = self._sire_object.property(space_prop) + except: + return + + # A rotation is only applicable to a triclinic space. + if not isinstance(space, _SireVol.TriclinicBox): + return + + # Rotate the space according to the desired precision. + space.rotate(precision) + + # Get the rotation matrix. + rotation_matrix = space.rotationMatrix() + + from sire.system import System + + # Create a cursor. + cursor = System(self._sire_object).cursor() + + # Rotate all vector properties. + + # Coordinates. + try: + prop_name = property_map.get("coordinates", "coordinates") + cursor = cursor.rotate( + matrix=rotation_matrix, map={"coordinates": prop_name} + ) + except: + pass + + # Velocities. + try: + prop_name = property_map.get("velocity", "velocity") + cursor = cursor.rotate( + matrix=rotation_matrix, map={"coordinates": prop_name} + ) + except: + pass + + # Now deal with any perturbable molecules. + if self.nPerturbableMolecules() > 0: + # Coordinates. + try: + prop_name = property_map.get("coordinates", "coordinates") + "0" + cursor = cursor.rotate( + matrix=rotation_matrix, map={"coordinates": prop_name} + ) + prop_name = property_map.get("coordinates", "coordinates") + "1" + cursor = cursor.rotate( + matrix=rotation_matrix, map={"coordinates": prop_name} + ) + except: + pass + + # Velocities. + try: + prop_name = property_map.get("velocity", "velocity") + "0" + cursor = cursor.rotate( + matrix=rotation_matrix, map={"coordinates": prop_name} + ) + prop_name = property_map.get("velocity", "velocity") + "1" + cursor = cursor.rotate( + matrix=rotation_matrix, map={"coordinates": prop_name} + ) + except: + pass + + # Commit the changes. + self._sire_object = cursor.commit()._system + + def reduceBoxVectors(self, bias=0, property_map={}): + """ + Reduce the box vectors. + + Parameters + ---------- + + bias : float + The bias to use when rounding during the lattice reduction. Negative + values biases towards left-tilting boxes, whereas positive values + biases towards right-tilting boxes. This can be used to ensure that + rounding is performed in a consistent direction, avoiding oscillation + when the box is instantiated from box vectors, or dimensions and + angles, that have been read from fixed-precision input files. + + property_map : dict + A dictionary that maps system "properties" to their user defined + values. This allows the user to refer to properties with their + own naming scheme, e.g. { "charge" : "my-charge" } + """ + + # First check that the system has a space. + space_prop = property_map.get("space", "space") + try: + space = self._sire_object.property(space_prop) + except: + return + + # A rotation is only applicable to a triclinic space. + if not isinstance(space, _SireVol.TriclinicBox): + return + + # Reduce the space using the desired bias. + space.reduce(bias) + + # Update the space property in the sire object. + self._sire_object.setProperty(space_prop, space) + def repartitionHydrogenMass( self, factor=4, water="no", use_coordinates=False, property_map={} ): From d01d82ff491661d5837956fe29e6cef995634624 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 7 Sep 2023 12:13:29 +0100 Subject: [PATCH 04/11] Allow user to rotate and/or reduce box on read. --- python/BioSimSpace/IO/_io.py | 40 ++++++++++++++++++- .../BioSimSpace/Sandpit/Exscientia/IO/_io.py | 40 ++++++++++++++++++- 2 files changed, 76 insertions(+), 4 deletions(-) diff --git a/python/BioSimSpace/IO/_io.py b/python/BioSimSpace/IO/_io.py index 922426d11..0ca90bce8 100644 --- a/python/BioSimSpace/IO/_io.py +++ b/python/BioSimSpace/IO/_io.py @@ -342,7 +342,13 @@ def readPDB(id, pdb4amber=False, work_dir=None, show_warnings=False, property_ma def readMolecules( - files, make_whole=False, show_warnings=False, download_dir=None, property_map={} + files, + make_whole=False, + rotate_box=False, + reduce_box=False, + show_warnings=False, + download_dir=None, + property_map={}, ): """ Read a molecular system from file. @@ -359,6 +365,17 @@ def readMolecules( Whether to make molecules whole, i.e. unwrap those that are split across the periodic boundary. + rotate_box : bool + Rotate the box vectors of the system so that the first vector is + aligned with the x-axis, the second vector lies in the x-y plane, + and the third vector has a positive z-component. This is a requirement + of certain molecular dynamics engines and is used for simulation + efficiency. All vector properties of the system, e.g. positions, + velocities, etc., will be rotated accordingly. + + reduce_box : bool + Reduce the box vectors. + show_warnings : bool Whether to show any warnings raised during parsing of the input files. @@ -439,6 +456,14 @@ def readMolecules( if make_whole: property_map["make_whole"] = _SireBase.wrap(make_whole) + # Validate the box rotation flag. + if not isinstance(rotate_box, bool): + raise TypeError("'rotate_box' must be of type 'bool'.") + + # Validate the box reduction flag. + if not isinstance(reduce_box, bool): + raise TypeError("'reduce_box' must be of type 'bool'.") + # Validate the warning message flag. if not isinstance(show_warnings, bool): raise TypeError("'show_warnings' must be of type 'bool'.") @@ -545,7 +570,18 @@ def readMolecules( except: pass - return _System(system) + # Wrap the Sire system. + system = _System(system) + + # Rotate the box vectors. + if rotate_box: + system.rotateBoxVectors() + + # Reduce the box vectors. + if reduce_box: + system.reduceBoxVectors() + + return system def saveMolecules(filebase, system, fileformat, property_map={}, **kwargs): diff --git a/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py b/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py index 922426d11..0ca90bce8 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py +++ b/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py @@ -342,7 +342,13 @@ def readPDB(id, pdb4amber=False, work_dir=None, show_warnings=False, property_ma def readMolecules( - files, make_whole=False, show_warnings=False, download_dir=None, property_map={} + files, + make_whole=False, + rotate_box=False, + reduce_box=False, + show_warnings=False, + download_dir=None, + property_map={}, ): """ Read a molecular system from file. @@ -359,6 +365,17 @@ def readMolecules( Whether to make molecules whole, i.e. unwrap those that are split across the periodic boundary. + rotate_box : bool + Rotate the box vectors of the system so that the first vector is + aligned with the x-axis, the second vector lies in the x-y plane, + and the third vector has a positive z-component. This is a requirement + of certain molecular dynamics engines and is used for simulation + efficiency. All vector properties of the system, e.g. positions, + velocities, etc., will be rotated accordingly. + + reduce_box : bool + Reduce the box vectors. + show_warnings : bool Whether to show any warnings raised during parsing of the input files. @@ -439,6 +456,14 @@ def readMolecules( if make_whole: property_map["make_whole"] = _SireBase.wrap(make_whole) + # Validate the box rotation flag. + if not isinstance(rotate_box, bool): + raise TypeError("'rotate_box' must be of type 'bool'.") + + # Validate the box reduction flag. + if not isinstance(reduce_box, bool): + raise TypeError("'reduce_box' must be of type 'bool'.") + # Validate the warning message flag. if not isinstance(show_warnings, bool): raise TypeError("'show_warnings' must be of type 'bool'.") @@ -545,7 +570,18 @@ def readMolecules( except: pass - return _System(system) + # Wrap the Sire system. + system = _System(system) + + # Rotate the box vectors. + if rotate_box: + system.rotateBoxVectors() + + # Reduce the box vectors. + if reduce_box: + system.reduceBoxVectors() + + return system def saveMolecules(filebase, system, fileformat, property_map={}, **kwargs): From a04ec64425acb6582f11e35ac4987536635f6f96 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 7 Sep 2023 12:15:28 +0100 Subject: [PATCH 05/11] Use correct property name to avoid confusion. --- python/BioSimSpace/IO/_io.py | 2 +- python/BioSimSpace/Sandpit/Exscientia/IO/_io.py | 2 +- python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py | 2 +- python/BioSimSpace/_SireWrappers/_system.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/python/BioSimSpace/IO/_io.py b/python/BioSimSpace/IO/_io.py index 0ca90bce8..193747290 100644 --- a/python/BioSimSpace/IO/_io.py +++ b/python/BioSimSpace/IO/_io.py @@ -370,7 +370,7 @@ def readMolecules( aligned with the x-axis, the second vector lies in the x-y plane, and the third vector has a positive z-component. This is a requirement of certain molecular dynamics engines and is used for simulation - efficiency. All vector properties of the system, e.g. positions, + efficiency. All vector properties of the system, e.g. coordinates, velocities, etc., will be rotated accordingly. reduce_box : bool diff --git a/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py b/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py index 0ca90bce8..193747290 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py +++ b/python/BioSimSpace/Sandpit/Exscientia/IO/_io.py @@ -370,7 +370,7 @@ def readMolecules( aligned with the x-axis, the second vector lies in the x-y plane, and the third vector has a positive z-component. This is a requirement of certain molecular dynamics engines and is used for simulation - efficiency. All vector properties of the system, e.g. positions, + efficiency. All vector properties of the system, e.g. coordinates, velocities, etc., will be rotated accordingly. reduce_box : bool diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py index 680c612d6..97dac39bb 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py @@ -1121,7 +1121,7 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): aligned with the x-axis, the second vector lies in the x-y plane, and the third vector has a positive z-component. This is a requirement of certain molecular dynamics engines and is used for simulation - efficiency. All vector properties of the system, e.g. positions, + efficiency. All vector properties of the system, e.g. coordinates, velocities, etc., will be rotated accordingly. Parameters diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py index 7d42a5c32..d350b6907 100644 --- a/python/BioSimSpace/_SireWrappers/_system.py +++ b/python/BioSimSpace/_SireWrappers/_system.py @@ -1069,7 +1069,7 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): aligned with the x-axis, the second vector lies in the x-y plane, and the third vector has a positive z-component. This is a requirement of certain molecular dynamics engines and is used for simulation - efficiency. All vector properties of the system, e.g. positions, + efficiency. All vector properties of the system, e.g. coordinates, velocities, etc., will be rotated accordingly. Parameters From e3dced0e7e553eb7b79d9276a680e7c8ac94688c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 7 Sep 2023 14:47:55 +0100 Subject: [PATCH 06/11] Update space property in system following rotation. --- python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py | 3 +++ python/BioSimSpace/_SireWrappers/_system.py | 3 +++ 2 files changed, 6 insertions(+) diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py index 97dac39bb..8adb2010a 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py @@ -1152,6 +1152,9 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): # Rotate the space according to the desired precision. space.rotate(precision) + # Update the space property in the sire object. + self._sire_object.setProperty(space_prop, space) + # Get the rotation matrix. rotation_matrix = space.rotationMatrix() diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py index d350b6907..02975301e 100644 --- a/python/BioSimSpace/_SireWrappers/_system.py +++ b/python/BioSimSpace/_SireWrappers/_system.py @@ -1100,6 +1100,9 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): # Rotate the space according to the desired precision. space.rotate(precision) + # Update the space property in the sire object. + self._sire_object.setProperty(space_prop, space) + # Get the rotation matrix. rotation_matrix = space.rotationMatrix() From 1abdc4e057c42490f566afbc4f25d47ec44a1c9c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 7 Sep 2023 14:50:56 +0100 Subject: [PATCH 07/11] Grammar tweak. --- .../BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py | 4 ++-- python/BioSimSpace/_SireWrappers/_system.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py index 8adb2010a..a7ddccefb 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py @@ -1145,7 +1145,7 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): except: return - # A rotation is only applicable to a triclinic space. + # A rotation is only applicable to triclinic spaces. if not isinstance(space, _SireVol.TriclinicBox): return @@ -1242,7 +1242,7 @@ def reduceBoxVectors(self, bias=0, property_map={}): except: return - # A rotation is only applicable to a triclinic space. + # A rotation is only applicable to triclinic spaces. if not isinstance(space, _SireVol.TriclinicBox): return diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py index 02975301e..490414475 100644 --- a/python/BioSimSpace/_SireWrappers/_system.py +++ b/python/BioSimSpace/_SireWrappers/_system.py @@ -1093,7 +1093,7 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): except: return - # A rotation is only applicable to a triclinic space. + # A rotation is only applicable to triclinic spaces. if not isinstance(space, _SireVol.TriclinicBox): return @@ -1190,7 +1190,7 @@ def reduceBoxVectors(self, bias=0, property_map={}): except: return - # A rotation is only applicable to a triclinic space. + # A rotation is only applicable to triclinic spaces. if not isinstance(space, _SireVol.TriclinicBox): return From f70ee551787184d9ea07d555c6ace41e504e8865 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 12 Sep 2023 11:58:23 +0100 Subject: [PATCH 08/11] Allow user to specify the box origin for rotation of vector properties. --- .../Exscientia/_SireWrappers/_system.py | 47 ++++++++++++++++--- python/BioSimSpace/_SireWrappers/_system.py | 47 ++++++++++++++++--- 2 files changed, 80 insertions(+), 14 deletions(-) diff --git a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py index a7ddccefb..29a8d44c4 100644 --- a/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py +++ b/python/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_system.py @@ -41,6 +41,7 @@ from .. import _isVerbose from .._Exceptions import IncompatibleError as _IncompatibleError from ..Types import Angle as _Angle +from ..Types import Coordinate as _Coordinate from ..Types import Length as _Length from .. import Units as _Units @@ -1115,7 +1116,16 @@ def nMLMolecules(self): """ return len(self.getMLMolecules()) - def rotateBoxVectors(self, precision=0.0, property_map={}): + def rotateBoxVectors( + self, + origin=_Coordinate( + 0 * _Units.Length.angstrom, + 0 * _Units.Length.angstrom, + 0 * _Units.Length.angstrom, + ), + precision=0.0, + property_map={}, + ): """ Rotate the box vectors of the system so that the first vector is aligned with the x-axis, the second vector lies in the x-y plane, @@ -1127,6 +1137,10 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): Parameters ---------- + origin : :class:`Coordinate ` + The origin of the triclinic space. This is the point about which + the lattice vectors are rotated. + precision : float The precision to use when sorting box vectors by their magnitude. This can be used to prevent unecessary box rotation when the @@ -1149,6 +1163,14 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): if not isinstance(space, _SireVol.TriclinicBox): return + # Validate input. + + if not isinstance(origin, _Coordinate): + raise TypeError("'origin' must be of type 'BioSimSpace.Types.Coordinate'") + + if not isinstance(precision, float): + raise TypeError("'precision' must be of type 'float'") + # Rotate the space according to the desired precision. space.rotate(precision) @@ -1158,6 +1180,9 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): # Get the rotation matrix. rotation_matrix = space.rotationMatrix() + # Get the center of rotation, as a Sire vector. + center = origin.toVector()._sire_object + from sire.system import System # Create a cursor. @@ -1169,7 +1194,7 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): try: prop_name = property_map.get("coordinates", "coordinates") cursor = cursor.rotate( - matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, matrix=rotation_matrix, map={"coordinates": prop_name} ) except: pass @@ -1178,7 +1203,7 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): try: prop_name = property_map.get("velocity", "velocity") cursor = cursor.rotate( - matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, matrix=rotation_matrix, map={"coordinates": prop_name} ) except: pass @@ -1189,11 +1214,15 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): try: prop_name = property_map.get("coordinates", "coordinates") + "0" cursor = cursor.rotate( - matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, + matrix=rotation_matrix, + map={"coordinates": prop_name}, ) prop_name = property_map.get("coordinates", "coordinates") + "1" cursor = cursor.rotate( - matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, + matrix=rotation_matrix, + map={"coordinates": prop_name}, ) except: pass @@ -1202,11 +1231,15 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): try: prop_name = property_map.get("velocity", "velocity") + "0" cursor = cursor.rotate( - matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, + matrix=rotation_matrix, + map={"coordinates": prop_name}, ) prop_name = property_map.get("velocity", "velocity") + "1" cursor = cursor.rotate( - matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, + matrix=rotation_matrix, + map={"coordinates": prop_name}, ) except: pass diff --git a/python/BioSimSpace/_SireWrappers/_system.py b/python/BioSimSpace/_SireWrappers/_system.py index 490414475..2ff728270 100644 --- a/python/BioSimSpace/_SireWrappers/_system.py +++ b/python/BioSimSpace/_SireWrappers/_system.py @@ -41,6 +41,7 @@ from .. import _isVerbose from .._Exceptions import IncompatibleError as _IncompatibleError from ..Types import Angle as _Angle +from ..Types import Coordinate as _Coordinate from ..Types import Length as _Length from .. import Units as _Units @@ -1063,7 +1064,16 @@ def nPerturbableMolecules(self): """ return len(self.getPerturbableMolecules()) - def rotateBoxVectors(self, precision=0.0, property_map={}): + def rotateBoxVectors( + self, + origin=_Coordinate( + 0 * _Units.Length.angstrom, + 0 * _Units.Length.angstrom, + 0 * _Units.Length.angstrom, + ), + precision=0.0, + property_map={}, + ): """ Rotate the box vectors of the system so that the first vector is aligned with the x-axis, the second vector lies in the x-y plane, @@ -1075,6 +1085,10 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): Parameters ---------- + origin : :class:`Coordinate ` + The origin of the triclinic space. This is the point about which + the lattice vectors are rotated. + precision : float The precision to use when sorting box vectors by their magnitude. This can be used to prevent unecessary box rotation when the @@ -1097,6 +1111,14 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): if not isinstance(space, _SireVol.TriclinicBox): return + # Validate input. + + if not isinstance(origin, _Coordinate): + raise TypeError("'origin' must be of type 'BioSimSpace.Types.Coordinate'") + + if not isinstance(precision, float): + raise TypeError("'precision' must be of type 'float'") + # Rotate the space according to the desired precision. space.rotate(precision) @@ -1106,6 +1128,9 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): # Get the rotation matrix. rotation_matrix = space.rotationMatrix() + # Get the center of rotation, as a Sire vector. + center = origin.toVector()._sire_object + from sire.system import System # Create a cursor. @@ -1117,7 +1142,7 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): try: prop_name = property_map.get("coordinates", "coordinates") cursor = cursor.rotate( - matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, matrix=rotation_matrix, map={"coordinates": prop_name} ) except: pass @@ -1126,7 +1151,7 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): try: prop_name = property_map.get("velocity", "velocity") cursor = cursor.rotate( - matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, matrix=rotation_matrix, map={"coordinates": prop_name} ) except: pass @@ -1137,11 +1162,15 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): try: prop_name = property_map.get("coordinates", "coordinates") + "0" cursor = cursor.rotate( - matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, + matrix=rotation_matrix, + map={"coordinates": prop_name}, ) prop_name = property_map.get("coordinates", "coordinates") + "1" cursor = cursor.rotate( - matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, + matrix=rotation_matrix, + map={"coordinates": prop_name}, ) except: pass @@ -1150,11 +1179,15 @@ def rotateBoxVectors(self, precision=0.0, property_map={}): try: prop_name = property_map.get("velocity", "velocity") + "0" cursor = cursor.rotate( - matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, + matrix=rotation_matrix, + map={"coordinates": prop_name}, ) prop_name = property_map.get("velocity", "velocity") + "1" cursor = cursor.rotate( - matrix=rotation_matrix, map={"coordinates": prop_name} + center=center, + matrix=rotation_matrix, + map={"coordinates": prop_name}, ) except: pass From 38a88e0a6193456edac28c9bfdc8c375693fb0c7 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 12 Sep 2023 12:27:34 +0100 Subject: [PATCH 09/11] Add unit test for triclinic box vector rotation. --- .../Exscientia/_SireWrappers/test_system.py | 51 +++++++++++++++++++ tests/_SireWrappers/test_system.py | 51 +++++++++++++++++++ 2 files changed, 102 insertions(+) diff --git a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py index f715ee05a..d5a5a0a50 100644 --- a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py +++ b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py @@ -1,6 +1,8 @@ import math import pytest +from sire.legacy.Vol import TriclinicBox + import BioSimSpace.Sandpit.Exscientia as BSS from tests.Sandpit.Exscientia.conftest import url, has_amber, has_openff @@ -392,3 +394,52 @@ def test_velocity_removal(): # Check that no molecules have a velocity property. assert len(new_system.search("not mol with property velocity").molecules()) == 2 + + +def test_rotate_box_vectors(system): + # Make sure that the triclinic space is rotated correctly and that vector + # properties (here coordinates) are updated accordingly. + + # Store the initial coordinates for the first molecule. + coords0 = system[0].coordinates() + + # Extract the space property. + space = system._sire_object.property("space") + + # Store the box matrix. + box_matrix = space.boxMatrix() + + # Create a triclinic box using the box matrix. This will create: + # TriclinicBox( ( 31.3979, 0, 0 ), ( 0, 34.1, 0 ), ( 0, 0, 29.273 ) ) + box = TriclinicBox(box_matrix.column0(), box_matrix.column1(), box_matrix.column2()) + + # Update the space property. + system._sire_object.setProperty("space", box) + + # Store the current box. + box0, angles0 = system.getBox() + + # Rotate the box. This will create an "optimal" triclinic space, where the + # shortest vector is aligned with the x axis, the second vector is in the + # x-y plane, and the largest has positive z component. + system.rotateBoxVectors() + + # Store the updated box. + box1, angles1 = system.getBox() + + # Store the updated coordinates for the first molecule. + coords1 = system[0].coordinates() + + # Make sure the box dimensions have been rotated correctly. + assert math.isclose(box1[0].value(), box0[2].value()) + assert math.isclose(box1[1].value(), box0[0].value()) + assert math.isclose(box1[2].value(), box0[1].value()) + assert math.isclose(angles1[0].value(), angles0[0].value()) + assert math.isclose(angles1[1].value(), angles0[1].value()) + assert math.isclose(angles1[2].value(), angles0[2].value()) + + # Make sure the coordinates have been rotated correctly. + for c0, c1 in zip(coords0, coords1): + 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()) diff --git a/tests/_SireWrappers/test_system.py b/tests/_SireWrappers/test_system.py index acfe504bf..3a351711c 100644 --- a/tests/_SireWrappers/test_system.py +++ b/tests/_SireWrappers/test_system.py @@ -1,6 +1,8 @@ import math import pytest +from sire.legacy.Vol import TriclinicBox + import BioSimSpace as BSS from tests.conftest import url, has_amber, has_openff @@ -391,3 +393,52 @@ def test_velocity_removal(): # Check that no molecules have a velocity property. assert len(new_system.search("not mol with property velocity").molecules()) == 2 + + +def test_rotate_box_vectors(system): + # Make sure that the triclinic space is rotated correctly and that vector + # properties (here coordinates) are updated accordingly. + + # Store the initial coordinates for the first molecule. + coords0 = system[0].coordinates() + + # Extract the space property. + space = system._sire_object.property("space") + + # Store the box matrix. + box_matrix = space.boxMatrix() + + # Create a triclinic box using the box matrix. This will create: + # TriclinicBox( ( 31.3979, 0, 0 ), ( 0, 34.1, 0 ), ( 0, 0, 29.273 ) ) + box = TriclinicBox(box_matrix.column0(), box_matrix.column1(), box_matrix.column2()) + + # Update the space property. + system._sire_object.setProperty("space", box) + + # Store the current box. + box0, angles0 = system.getBox() + + # Rotate the box. This will create an "optimal" triclinic space, where the + # shortest vector is aligned with the x axis, the second vector is in the + # x-y plane, and the largest has positive z component. + system.rotateBoxVectors() + + # Store the updated box. + box1, angles1 = system.getBox() + + # Store the updated coordinates for the first molecule. + coords1 = system[0].coordinates() + + # Make sure the box dimensions have been rotated correctly. + assert math.isclose(box1[0].value(), box0[2].value()) + assert math.isclose(box1[1].value(), box0[0].value()) + assert math.isclose(box1[2].value(), box0[1].value()) + assert math.isclose(angles1[0].value(), angles0[0].value()) + assert math.isclose(angles1[1].value(), angles0[1].value()) + assert math.isclose(angles1[2].value(), angles0[2].value()) + + # Make sure the coordinates have been rotated correctly. + for c0, c1 in zip(coords0, coords1): + 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()) From b392af957a7ec025923653fee58c1fb4d6850403 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 12 Sep 2023 14:23:40 +0100 Subject: [PATCH 10/11] Add comment regarding unchanged angles. --- tests/Sandpit/Exscientia/_SireWrappers/test_system.py | 2 ++ tests/_SireWrappers/test_system.py | 2 ++ 2 files changed, 4 insertions(+) diff --git a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py index d5a5a0a50..79054ff9d 100644 --- a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py +++ b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py @@ -434,6 +434,8 @@ def test_rotate_box_vectors(system): assert math.isclose(box1[0].value(), box0[2].value()) assert math.isclose(box1[1].value(), box0[0].value()) assert math.isclose(box1[2].value(), box0[1].value()) + # All angles should remain the same since this is a cubic system with + # 90 degree angles. assert math.isclose(angles1[0].value(), angles0[0].value()) assert math.isclose(angles1[1].value(), angles0[1].value()) assert math.isclose(angles1[2].value(), angles0[2].value()) diff --git a/tests/_SireWrappers/test_system.py b/tests/_SireWrappers/test_system.py index 3a351711c..358c4dec6 100644 --- a/tests/_SireWrappers/test_system.py +++ b/tests/_SireWrappers/test_system.py @@ -433,6 +433,8 @@ def test_rotate_box_vectors(system): assert math.isclose(box1[0].value(), box0[2].value()) assert math.isclose(box1[1].value(), box0[0].value()) assert math.isclose(box1[2].value(), box0[1].value()) + # All angles should remain the same since this is a cubic system with + # 90 degree angles. assert math.isclose(angles1[0].value(), angles0[0].value()) assert math.isclose(angles1[1].value(), angles0[1].value()) assert math.isclose(angles1[2].value(), angles0[2].value()) From c74693a4afb989df5268a5003dd83232739412f9 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 22 Sep 2023 12:51:12 +0100 Subject: [PATCH 11/11] Operate on a copy of the system to avoid breaking later test. --- tests/Sandpit/Exscientia/_SireWrappers/test_system.py | 3 +++ tests/_SireWrappers/test_system.py | 3 +++ 2 files changed, 6 insertions(+) diff --git a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py index 79054ff9d..454f15b03 100644 --- a/tests/Sandpit/Exscientia/_SireWrappers/test_system.py +++ b/tests/Sandpit/Exscientia/_SireWrappers/test_system.py @@ -265,6 +265,9 @@ def test_get_box(system): def test_set_box(system): + # Make a copy of the system. + system = system.copy() + # Generate box dimensions and angles for a truncated octahedron. box, angles = BSS.Box.truncatedOctahedron(30 * BSS.Units.Length.angstrom) diff --git a/tests/_SireWrappers/test_system.py b/tests/_SireWrappers/test_system.py index 358c4dec6..3f8cdb126 100644 --- a/tests/_SireWrappers/test_system.py +++ b/tests/_SireWrappers/test_system.py @@ -268,6 +268,9 @@ def test_set_box(system): # Generate box dimensions and angles for a truncated octahedron. box, angles = BSS.Box.truncatedOctahedron(30 * BSS.Units.Length.angstrom) + # Make a copy of the system. + system = system.copy() + # Set the box dimensions in the system. system.setBox(box, angles)