updated

Author: mushroomfire

mdapy/cell_opt.py

@@ -111,14 +111,33 @@ def compute(self):
             lmp.commands_string(f"atom_style {self.atomic_style}")
             num_type = self.type_list.max()
             create_box = f"""
-            region 1 block 0 1 0 1 0 1
+            region 1 prism 0 1 0 1 0 1 0 0 0
             create_box {num_type} 1
             """
             lmp.commands_string(create_box)
-            lmp.reset_box(*self.to_lammps_box(self.box))
+            pos = self.pos
+            box = self.box
+            if box[0, 1] != 0 or box[0, 2] != 0 or box[1, 2] != 0:
+                old_box = box.copy()
+                ax = np.linalg.norm(box[0])
+                bx = box[1] @ (box[0] / ax)
+                by = np.sqrt(np.linalg.norm(box[1]) ** 2 - bx**2)
+                cx = box[2] @ (box[0] / ax)
+                cy = (box[1] @ box[2] - bx * cx) / by
+                cz = np.sqrt(np.linalg.norm(box[2]) ** 2 - cx**2 - cy**2)
+                box = np.array([[ax, bx, cx], [0, by, cy], [0, 0, cz]]).T
+                rotation = np.linalg.solve(old_box[:-1], box)
+                pos = self.pos @ rotation
+                box = np.r_[box, box[-1].reshape(1, -1)]
+            # lmp.reset_box(*self.to_lammps_box(box))
+            lo, hi, xy, xz, yz = self.to_lammps_box(box)
+            lmp.commands_string(
+                f"change_box all x final {lo[0]} {hi[0]} y final {lo[1]} {hi[1]} z final {lo[2]} {hi[2]} xy final {xy} xz final {xz} yz final {yz}"
+            )
+
             N = self.pos.shape[0]
             N_lmp = lmp.create_atoms(
-                N, np.arange(1, N + 1), self.type_list, self.pos.flatten()
+                N, np.arange(1, N + 1), self.type_list, pos.flatten()
             )
             assert N == N_lmp, "Wrong atom numbers."
 

mdapy/system.py

@@ -728,6 +728,7 @@ def cal_phono_dispersion(
         labels,
         potential,
         elements_list,
+        symprec=1e-5,
         replicate=None,
     ):
         """This function can be used to calculate the phono dispersion based on Phonopy (https://phonopy.github.io/phonopy/). We support NEP and
@@ -739,6 +740,7 @@ def cal_phono_dispersion(
             potential (BasePotential): base potential class defined in mdapy, which must including a compute method to calculate the energy, force, virial.
             elements_list (list[str]): element list, such as ['Al']
             pair_style (str, optional): pair style, selected in ['nep', 'eam/alloy']. Defaults to "eam/alloy".
+            symprec (float): this is used to set geometric tolerance to find symmetry of crystal structure. Defaults to 1e-5.
             replicate (list, optional): replication to pos, such as [3, 3, 3]. If not given, we will replicate it exceeding 15 A per directions. Defaults to None.
 
         Outputs:
@@ -761,6 +763,7 @@ def cal_phono_dispersion(
             self.box,
             elements_list,
             self.__data["type"].to_numpy(),
+            symprec,
             replicate,
         )
         self.Phon.compute()