lammps: rotate cell to lower_triangle before writing to lammps (#842)

<!-- This is an auto-generated comment: release notes by coderabbit.ai
-->
## Summary by CodeRabbit

* **New Features**
* Cell matrices and atomic coordinates are now automatically transformed
to a lower triangular form, ensuring consistent output formatting.

* **Improvements**
* All cell parameter and atomic position outputs now reflect the
rotated, standardized cell representation.

* **Tests**
* Added comprehensive tests verifying the lower triangular
transformation and preservation of atomic distances across various cell
configurations.
<!-- end of auto-generated comment: release notes by coderabbit.ai -->

---------

Co-authored-by: root <pxlxingliang>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: pxlxingliang

dpdata/lammps/lmp.py

@@ -193,31 +193,64 @@ def to_system_data(lines, type_map=None, type_idx_zero=True):
     return system_data(lines, type_map=type_map, type_idx_zero=type_idx_zero)
 
 
+def rotate_to_lower_triangle(
+    cell: np.ndarray, coord: np.ndarray
+) -> tuple[np.ndarray, np.ndarray]:
+    """Rotate the cell to lower triangular and ensure the diagonal elements are non-negative.
+
+    Args:
+        cell (np.ndarray): The original cell matrix.
+        coord (np.ndarray): The coordinates of the atoms.
+
+    Returns
+    -------
+    tuple[np.ndarray, np.ndarray]: The rotated cell and adjusted coordinates.
+    """
+    q, _ = np.linalg.qr(cell.T)
+    cell = np.matmul(cell, q)
+    coord = np.matmul(coord, q)
+
+    # Ensure the diagonal elements of the cell are non-negative
+    rot = np.eye(3)
+    if cell[0][0] < 0:
+        rot[0][0] = -1
+    if cell[1][1] < 0:
+        rot[1][1] = -1
+    if cell[2][2] < 0:
+        rot[2][2] = -1
+    cell = np.matmul(cell, rot)
+    coord = np.matmul(coord, rot)
+    return cell, coord
+
+
 def from_system_data(system, f_idx=0):
     ret = ""
     ret += "\n"
     natoms = sum(system["atom_numbs"])
     ntypes = len(system["atom_numbs"])
+    cell, coord = rotate_to_lower_triangle(
+        system["cells"][f_idx], system["coords"][f_idx]
+    )
     ret += "%d atoms\n" % natoms  # noqa: UP031
     ret += "%d atom types\n" % ntypes  # noqa: UP031
     ret += (ptr_float_fmt + " " + ptr_float_fmt + " xlo xhi\n") % (
         0,
-        system["cells"][f_idx][0][0],
+        cell[0][0],
     )  # noqa: UP031
     ret += (ptr_float_fmt + " " + ptr_float_fmt + " ylo yhi\n") % (
         0,
-        system["cells"][f_idx][1][1],
+        cell[1][1],
     )  # noqa: UP031
     ret += (ptr_float_fmt + " " + ptr_float_fmt + " zlo zhi\n") % (
         0,
-        system["cells"][f_idx][2][2],
+        cell[2][2],
     )  # noqa: UP031
     ret += (
         ptr_float_fmt + " " + ptr_float_fmt + " " + ptr_float_fmt + " xy xz yz\n"
     ) % (
-        system["cells"][f_idx][1][0],
-        system["cells"][f_idx][2][0],
-        system["cells"][f_idx][2][1],
+        cell[1][0],
+        cell[2][0],
+        cell[2][1],
     )  # noqa: UP031
     ret += "\n"
     ret += "Atoms # atomic\n"
@@ -255,9 +288,9 @@ def from_system_data(system, f_idx=0):
                 ret += coord_fmt % (
                     ii + 1,
                     system["atom_types"][ii] + 1,
-                    system["coords"][f_idx][ii][0] - system["orig"][0],
-                    system["coords"][f_idx][ii][1] - system["orig"][1],
-                    system["coords"][f_idx][ii][2] - system["orig"][2],
+                    coord[ii][0] - system["orig"][0],
+                    coord[ii][1] - system["orig"][1],
+                    coord[ii][2] - system["orig"][2],
                     system["spins"][f_idx][ii][0] / spins_norm[ii],
                     system["spins"][f_idx][ii][1] / spins_norm[ii],
                     system["spins"][f_idx][ii][2] / spins_norm[ii],
@@ -267,9 +300,9 @@ def from_system_data(system, f_idx=0):
                 ret += coord_fmt % (
                     ii + 1,
                     system["atom_types"][ii] + 1,
-                    system["coords"][f_idx][ii][0] - system["orig"][0],
-                    system["coords"][f_idx][ii][1] - system["orig"][1],
-                    system["coords"][f_idx][ii][2] - system["orig"][2],
+                    coord[ii][0] - system["orig"][0],
+                    coord[ii][1] - system["orig"][1],
+                    coord[ii][2] - system["orig"][2],
                     system["spins"][f_idx][ii][0],
                     system["spins"][f_idx][ii][1],
                     system["spins"][f_idx][ii][2] + 1,
@@ -279,9 +312,9 @@ def from_system_data(system, f_idx=0):
             ret += coord_fmt % (
                 ii + 1,
                 system["atom_types"][ii] + 1,
-                system["coords"][f_idx][ii][0] - system["orig"][0],
-                system["coords"][f_idx][ii][1] - system["orig"][1],
-                system["coords"][f_idx][ii][2] - system["orig"][2],
+                coord[ii][0] - system["orig"][0],
+                coord[ii][1] - system["orig"][1],
+                coord[ii][2] - system["orig"][2],
             )  # noqa: UP031
     return ret
 

tests/test_lammps_lmp_dump.py

@@ -3,9 +3,12 @@
 import os
 import unittest
 
+import numpy as np
 from context import dpdata
 from poscars.poscar_ref_oh import TestPOSCARoh
 
+from dpdata.lammps.lmp import rotate_to_lower_triangle
+
 
 class TestLmpDump(unittest.TestCase, TestPOSCARoh):
     def setUp(self):
@@ -27,5 +30,65 @@ def setUp(self):
         self.system.from_fmt("tmp.lmp", fmt="lammps/lmp", type_map=["O", "H"])
 
 
+class TestLmpRotateTriAngle(unittest.TestCase):
+    def test_simple_cubic(self):
+        cubic_cell = np.diag([5, 5, 5])
+        cubic_coords = np.array([[0, 0, 0], [2.5, 2.5, 2.5]])
+        rotated_cell, rotated_coords = rotate_to_lower_triangle(
+            cubic_cell, cubic_coords
+        )
+
+        self.assertTrue(np.allclose(rotated_cell, np.diag([5, 5, 5])))
+        self.assertTrue(np.allclose(rotated_coords, cubic_coords))
+
+    def test_triclinic_cell(self):
+        triclinic_cell = np.array([[4.0, 0.0, 0.0], [1.0, 3.0, 0.0], [0.5, 0.5, 2.5]])
+
+        triclinic_coords = np.array([[0.5, 0.5, 0.5], [3.0, 2.0, 1.0]])
+        rotated_cell, rotated_coords = rotate_to_lower_triangle(
+            triclinic_cell.copy(), triclinic_coords.copy()
+        )
+
+        self.assertTrue(np.isclose(rotated_cell[0, 1], 0, atol=1e-10))
+        self.assertTrue(np.isclose(rotated_cell[0, 2], 0, atol=1e-10))
+        self.assertTrue(np.isclose(rotated_cell[1, 2], 0, atol=1e-10))
+
+        self.assertTrue(rotated_cell[0, 0] > 0)
+        self.assertTrue(rotated_cell[1, 1] > 0)
+        self.assertTrue(rotated_cell[2, 2] > 0)
+
+        # check the distance between two atoms
+        self.assertTrue(
+            np.isclose(
+                np.linalg.norm(rotated_coords[0] - rotated_coords[1]),
+                np.linalg.norm(triclinic_coords[0] - triclinic_coords[1]),
+                atol=1e-10,
+            )
+        )
+
+    def test_negative_diagonal(self):
+        negative_cell = np.array([[-3.0, 1.0, 0.5], [0.0, -2.0, 0.3], [0.0, 0.0, -4.0]])
+        negative_coords = np.array([[1.0, 0.5, -1.0], [0.5, 1.0, -2.0]])
+        rotated_cell, rotated_coords = rotate_to_lower_triangle(
+            negative_cell.copy(), negative_coords.copy()
+        )
+
+        self.assertTrue(np.isclose(rotated_cell[0, 1], 0, atol=1e-10))
+        self.assertTrue(np.isclose(rotated_cell[0, 2], 0, atol=1e-10))
+        self.assertTrue(np.isclose(rotated_cell[1, 2], 0, atol=1e-10))
+
+        self.assertTrue(rotated_cell[0, 0] > 0)
+        self.assertTrue(rotated_cell[1, 1] > 0)
+        self.assertTrue(rotated_cell[2, 2] > 0)
+
+        self.assertTrue(
+            np.isclose(
+                np.linalg.norm(negative_coords[0] - negative_coords[1]),
+                np.linalg.norm(rotated_coords[0] - rotated_coords[1]),
+                atol=1e-6,
+            )
+        )
+
+
 if __name__ == "__main__":
     unittest.main()