Merge pull request #1234 from glotzerlab/feat_from_to_lattice_vectors

Add class method for Box class which create a box from Lx,Ly,Lz and angles

ChangeLog.md

@@ -8,6 +8,7 @@ and this project adheres to
 
 ### Added
 * New continuous coordination number compute `freud.order.ContinuousCoordination`.
+* New methods for conversion of box lengths and angles to/from `freud.box.Box`.
 
 ### Removed
 * `freud.order.Translational`.

doc/source/reference/credits.rst

@@ -369,6 +369,10 @@ Domagoj Fijan
 * Contributed code, design, documentation, and testing for ``freud.locality.FilterSANN`` class.
 * Contributed code, design, documentation, and testing for ``freud.locality.FilterRAD`` class.
 * Added support for ``gsd.hoomd.Frame`` in ``NeighborQuery.from_system`` calls.
+* Added support for ``freud.box.Box`` class methods for construction of boxes from cell
+  lengths and angles (``freud.box.Box.from_box_lengths_and_angles``), as well as a
+  method for returning box vector lengths and angles
+  (``freud.box.Box.to_box_lengths_and_angles``).
 
 Andrew Kerr
 

freud/box.pyx

@@ -689,6 +689,27 @@ cdef class Box:
                            [0, self.Ly, self.yz * self.Lz],
                            [0, 0, self.Lz]])
 
+    def to_box_lengths_and_angles(self):
+        r"""Return the box lengths and angles.
+
+        Returns:
+            tuple:
+                The box vector lengths and angles in radians
+                :math:`(L_1, L_2, L_3, \alpha, \beta, \gamma)`.
+        """
+        alpha = np.arccos(
+            (self.xy * self.xz + self.yz)
+            / (np.sqrt(1 + self.xy**2) * np.sqrt(1 + self.xz**2 + self.yz**2))
+        )
+        beta = np.arccos(self.xz/np.sqrt(1+self.xz**2+self.yz**2))
+        gamma = np.arccos(self.xy/np.sqrt(1+self.xy**2))
+        L1 = self.Lx
+        a2 = [self.Ly*self.xy, self.Ly, 0]
+        a3 = [self.Lz*self.xz, self.Lz*self.yz, self.Lz]
+        L2 = np.linalg.norm(a2)
+        L3 = np.linalg.norm(a3)
+        return (L1, L2, L3, alpha, beta, gamma)
+
     def __repr__(self):
         return ("freud.box.{cls}(Lx={Lx}, Ly={Ly}, Lz={Lz}, "
                 "xy={xy}, xz={xz}, yz={yz}, "
@@ -921,6 +942,55 @@ cdef class Box:
                             "positional argument: L")
         return cls(Lx=L, Ly=L, Lz=0, xy=0, xz=0, yz=0, is2D=True)
 
+    @classmethod
+    def from_box_lengths_and_angles(
+        cls, L1, L2, L3, alpha, beta, gamma, dimensions=None,
+    ):
+        r"""Construct a box from lengths and angles (in radians).
+
+        All the angles provided must be between 0 and :math:`\pi`.
+
+        Args:
+            L1 (float):
+                The length of the first lattice vector.
+            L2 (float):
+                The length of the second lattice vector.
+            L3 (float):
+                The length of the third lattice vector.
+            alpha (float):
+                The angle between second and third lattice vector in radians (must be
+                between 0 and :math:`\pi`).
+            beta (float):
+                The angle between first and third lattice vector in radians (must be
+                between 0 and :math:`\pi`).
+            gamma (float):
+                The angle between the first and second lattice vector in radians (must
+                be between 0 and :math:`\pi`).
+            dimensions (int):
+                The number of dimensions (Default value = :code:`None`).
+
+        Returns:
+            :class:`freud.box.Box`: The resulting box object.
+        """
+        if not 0 < alpha < np.pi:
+            raise ValueError("alpha must be between 0 and pi.")
+        if not 0 < beta < np.pi:
+            raise ValueError("beta must be between 0 and pi.")
+        if not 0 < gamma < np.pi:
+            raise ValueError("gamma must be between 0 and pi.")
+        a1 = np.array([L1, 0, 0])
+        a2 = np.array([L2 * np.cos(gamma), L2 * np.sin(gamma), 0])
+        a3x = np.cos(beta)
+        a3y = (np.cos(alpha) - np.cos(beta) * np.cos(gamma)) / np.sin(gamma)
+        under_sqrt = 1 - a3x**2 - a3y**2
+        if under_sqrt < 0:
+            raise ValueError("The provided angles can not form a valid box.")
+        a3z = np.sqrt(under_sqrt)
+        a3 = np.array([L3 * a3x, L3 * a3y, L3 * a3z])
+        if dimensions is None:
+            dimensions = 2 if L3 == 0 else 3
+        return cls.from_matrix(np.array([a1, a2, a3]).T, dimensions=dimensions)
+
 
 cdef BoxFromCPP(const freud._box.Box & cppbox):
     b = Box(cppbox.getLx(), cppbox.getLy(), cppbox.getLz(),

tests/test_box_Box.py

@@ -485,6 +485,52 @@ def test_from_box(self):
         box7 = freud.box.Box.from_matrix(box.to_matrix())
         assert np.isclose(box.to_matrix(), box7.to_matrix()).all()
 
+    def test_standard_orthogonal_box(self):
+        box = freud.box.Box.from_box((1, 2, 3, 0, 0, 0))
+        Lx, Ly, Lz, alpha, beta, gamma = box.to_box_lengths_and_angles()
+        npt.assert_allclose(
+            (Lx, Ly, Lz, alpha, beta, gamma), (1, 2, 3, np.pi / 2, np.pi / 2, np.pi / 2)
+        )
+
+    def test_to_and_from_box_lengths_and_angles(self):
+        original_box_lengths_and_angles = (
+            np.random.uniform(0, 100000),
+            np.random.uniform(0, 100000),
+            np.random.uniform(0, 100000),
+            np.random.uniform(0, np.pi),
+            np.random.uniform(0, np.pi),
+            np.random.uniform(0, np.pi),
+        )
+        if (
+            1
+            - np.cos(original_box_lengths_and_angles[4]) ** 2
+            - (
+                (
+                    np.cos(original_box_lengths_and_angles[3])
+                    - np.cos(original_box_lengths_and_angles[4])
+                    * np.cos(original_box_lengths_and_angles[5])
+                )
+                / np.sin(original_box_lengths_and_angles[5])
+            )
+            ** 2
+            < 0
+        ):
+            with pytest.raises(ValueError):
+                freud.box.Box.from_box_lengths_and_angles(
+                    *original_box_lengths_and_angles
+                )
+        else:
+            box = freud.box.Box.from_box_lengths_and_angles(
+                *original_box_lengths_and_angles
+            )
+            lengths_and_angles_computed = box.to_box_lengths_and_angles()
+            np.testing.assert_allclose(
+                lengths_and_angles_computed,
+                original_box_lengths_and_angles,
+                rtol=1e-6,
+                atol=1e-14,
+            )
+
     def test_matrix(self):
         box = freud.box.Box(2, 2, 2, 1, 0.5, 0.1)
         box2 = freud.box.Box.from_matrix(box.to_matrix())