Merge pull request #175 from OpenBioSim/feature_reduce

Feature reduce

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" <https://research.rug.nl/files/2839528/01_c1.pdf>`__
-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
 

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. coordinates,
+        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(

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. coordinates,
+        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(

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,6 +1116,175 @@ def nMLMolecules(self):
         """
         return len(self.getMLMolecules())
 
+    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,
+        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. coordinates,
+        velocities, etc., will be rotated accordingly.
+
+        Parameters
+        ----------
+
+        origin : :class:`Coordinate <BioSimSpace.Types.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
+            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 triclinic spaces.
+        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)
+
+        # Update the space property in the sire object.
+        self._sire_object.setProperty(space_prop, space)
+
+        # 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.
+        cursor = System(self._sire_object).cursor()
+
+        # Rotate all vector properties.
+
+        # Coordinates.
+        try:
+            prop_name = property_map.get("coordinates", "coordinates")
+            cursor = cursor.rotate(
+                center=center, matrix=rotation_matrix, map={"coordinates": prop_name}
+            )
+        except:
+            pass
+
+        # Velocities.
+        try:
+            prop_name = property_map.get("velocity", "velocity")
+            cursor = cursor.rotate(
+                center=center, 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(
+                    center=center,
+                    matrix=rotation_matrix,
+                    map={"coordinates": prop_name},
+                )
+                prop_name = property_map.get("coordinates", "coordinates") + "1"
+                cursor = cursor.rotate(
+                    center=center,
+                    matrix=rotation_matrix,
+                    map={"coordinates": prop_name},
+                )
+            except:
+                pass
+
+            # Velocities.
+            try:
+                prop_name = property_map.get("velocity", "velocity") + "0"
+                cursor = cursor.rotate(
+                    center=center,
+                    matrix=rotation_matrix,
+                    map={"coordinates": prop_name},
+                )
+                prop_name = property_map.get("velocity", "velocity") + "1"
+                cursor = cursor.rotate(
+                    center=center,
+                    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 triclinic spaces.
+        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={}
     ):