Add support for creation of box object from non-triangular matrices

hoomd/box.py

@@ -160,8 +160,8 @@ class Box:
     .. rubric:: Factory Methods
 
     `Box` has factory methods to enable easier creation of boxes: `cube`,
-    `square`, `from_matrix`, and `from_box`. See each method's documentation for
-    more details.
+    `square`, `from_matrix`, `from_basis_vectors`, and `from_box`. See each
+    method's documentation for more details.
 
     .. rubric:: Example:
 
@@ -216,6 +216,86 @@ def square(cls, L):
         """
         return cls(L, L, 0, 0, 0, 0)
 
+    @classmethod
+    def from_basis_vectors(cls, box_matrix):
+        r"""Initialize a Box instance from a box matrix.
+
+        Args:
+            box_matrix ((3, 3) `numpy.ndarray` of `float`): A 3x3 matrix
+                or list of lists representing a set of lattice basis vectors.
+
+        Note:
+           The created box will be rotated with respect to the lattice basis. As
+           a consequence the output of `to_matrix` will not be the same as the
+           input provided to this function. This function also returns a
+           rotation matrix comensurate with this transformation. Using this
+           rotation matrix users can rotate the original points into the new box
+           by applying the rotation to the points.
+
+        Note:
+           When passing a 2D basis vectors, the third vector should be set to
+           all zeros, while first two vectors should have the last element set
+           to zero.
+
+        Returns:
+            tuple[hoomd.Box, numpy.ndarray]: A tuple containing:
+                - hoomd.Box: The created box configured according to the given
+                  basis vectors.
+                - numpy.ndarray: A 3x3 floating-point rotation matrix that can
+                  be used to transform the original basis vectors to align with
+                  the new box basis vectors.
+
+        .. rubric:: Example:
+
+        .. code-block:: python
+
+            points = np.array([[0, 0, 0], [0.5, 0, 0], [0.25, 0.25, 0]])
+            box, rotation = hoomd.Box.from_basis_vectors(
+                box_matrix = [[ 1,  1,  0],
+                              [ 1, -1,  0],
+                              [ 0,  0,  1]])
+            rotated_points = rotation @ points
+        """
+        box_matrix = np.asarray(box_matrix, dtype=np.float64)
+        if box_matrix.shape != (3, 3):
+            raise ValueError("Box matrix must be a 3x3 matrix.")
+        v0 = box_matrix[:, 0]
+        v1 = box_matrix[:, 1]
+        v2 = box_matrix[:, 2]
+        Lx = np.sqrt(np.dot(v0, v0))
+        a2x = np.dot(v0, v1) / Lx
+        Ly = np.sqrt(np.dot(v1, v1) - a2x * a2x)
+        xy = a2x / Ly
+        v0xv1 = np.cross(v0, v1)
+        v0xv1mag = np.sqrt(np.dot(v0xv1, v0xv1))
+        Lz = np.dot(v2, v0xv1) / v0xv1mag
+        if Lz != 0:
+            a3x = np.dot(v0, v2) / Lx
+            xz = a3x / Lz
+            yz = (np.dot(v1, v2) - a2x * a3x) / (Ly * Lz)
+            upper_triangular_box_matrix = np.array([[Lx, Ly * xy, Lz * xz],
+                                                    [0, Ly, Lz * yz],
+                                                    [0, 0, Lz]])
+        else:
+            xz = yz = 0
+            if not (np.allclose(v2, [0, 0, 0]) and np.allclose(v0[2], 0)
+                    and np.allclose(v1[2], 0)):
+                error_string = ("A 2D box matrix must have a third vector and"
+                                "third component of first two vectors set to"
+                                "zero.")
+                raise ValueError(error_string)
+            upper_triangular_box_matrix = np.array([[Lx, Ly * xy], [0, Ly]])
+            box_matrix = box_matrix[:2, :2]
+
+        rotation = np.linalg.solve(upper_triangular_box_matrix, box_matrix)
+
+        if Lz == 0:
+            rotation = np.zeros((3, 3))
+            rotation[:2, :2] = box_matrix
+            rotation[2, 2] = 1
+
+        return cls(Lx=Lx, Ly=Ly, Lz=Lz, xy=xy, xz=xz, yz=yz), rotation
+
     @classmethod
     def from_matrix(cls, box_matrix):
         r"""Create a box from an upper triangular matrix.
@@ -247,7 +327,7 @@ def from_matrix(cls, box_matrix):
                               [0, 8, 16],
                               [0, 0, 18]])
         """
-        box_matrix = np.asarray(box_matrix)
+        box_matrix = np.asarray(box_matrix, dtype=np.float64)
         if box_matrix.shape != (3, 3):
             raise ValueError("Box matrix must be a 3x3 matrix.")
         if not np.allclose(box_matrix, np.triu(box_matrix)):

hoomd/pytest/test_box.py

@@ -2,7 +2,9 @@
 # Part of HOOMD-blue, released under the BSD 3-Clause License.
 
 from math import isclose
+
 import numpy as np
+import pytest
 from pytest import fixture
 
 from hoomd.box import Box
@@ -172,6 +174,45 @@ def test_from_matrix(new_box_matrix_dict):
     ])
 
 
+@pytest.mark.parametrize("theta",
+                         [np.pi, np.pi / 2, np.pi / 3, np.pi / 4, np.pi * 1.23])
+def test_from_basis_vectors_non_triangular(theta):
+    box_matrix = np.array([[np.cos(theta), -np.sin(theta), 0],
+                           [np.sin(theta), np.cos(theta), 0], [0, 0, 1]])
+    box, rotation = Box.from_basis_vectors(box_matrix.T)
+    assert np.allclose([box.Lx, box.Ly, box.Lz, box.xy, box.xz, box.yz],
+                       [1, 1, 1, 0, 0, 0],
+                       atol=1e-6)
+    rotated_matrix = box.to_matrix()
+    rotated_points = rotation @ box_matrix
+    assert np.allclose(rotated_matrix, rotated_points)
+
+
+def test_from_matrix_two_dimensional():
+    box_matrix = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 0]])
+    box, _ = Box.from_basis_vectors(box_matrix)
+    assert box.is2D and box.dimensions == 2
+
+
+@pytest.mark.parametrize("theta",
+                         [np.pi, np.pi / 2, np.pi / 3, np.pi / 4, np.pi * 1.23])
+def test_rotation_matrix_from_basis_vectors_two_dimensional(theta):
+    box_matrix = np.array([[np.cos(theta), -np.sin(theta), 0],
+                           [np.sin(theta), np.cos(theta), 0], [0, 0, 0]])
+    box, rotation = Box.from_basis_vectors(box_matrix.T)
+    rotated_matrix = box.to_matrix()
+    rotated_points = rotation @ box_matrix
+    assert np.allclose(rotated_matrix, rotated_points)
+    assert box.is2D and box.dimensions == 2
+
+
+def test_invalid_from_basis_vectors_two_dimensional():
+    box_matrix = np.array([[1, 0, 0], [0, 1, 1], [0, 0, 0]])
+    import pytest
+    with pytest.raises(ValueError):
+        Box.from_basis_vectors(box_matrix)
+
+
 def test_eq(base_box, box_dict):
     box2 = Box(**box_dict)
     assert base_box == box2