Merge pull request #341 from KenHino/fix/hexagonal-site-symmetry

support Z'=1/3, 1/6

Author: qzhu2017

pyxtal/molecule.py

@@ -2137,7 +2137,7 @@ def get_orientations_in_wp(self, wp, rtol=1e-2):
         elif wp.index > 1 and self.pga.sch_symbol == "C1":
             return []
 
-        symm_w = wp.get_site_symm_wo_translation()  # symmetry without translation
+        symm_w = wp.get_site_symm_wo_translation(space="euclidean")  # symmetry without translation
         # molecule has fewer symops
         if len(self.pg[0]) < len(symm_w):
             return []
@@ -2637,7 +2637,7 @@ def is_compatible_symmetry(mol, wp):
     if len(mol) == 1 or wp.index == 0:
         return True
     pga = PointGroupAnalyzer(mol)
-    return all(pga.is_valid_op(op) for op in wp.get_site_symm_wo_translation())
+    return all(pga.is_valid_op(op) for op in wp.get_site_symm_wo_translation(space="euclidean"))
 
 
 def make_graph(mol, tol=0.2, ignore_HH=False, debug=False):

pyxtal/symmetry.py

@@ -15,6 +15,7 @@
 import re
 from copy import deepcopy
 from ast import literal_eval
+from typing import Literal
 
 import numpy as np
 from monty.serialization import loadfn
@@ -2465,8 +2466,27 @@ def equivalent_set(self, index):
         return self
 
     # =============================Get functions===========================
-    def get_site_symm_wo_translation(self):
-        return [SymmOp.from_rotation_and_translation(op.rotation_matrix, [0, 0, 0]) for op in self.symmetry[0]]
+    def get_site_symm_wo_translation(self, space: Literal["fractional", "euclidean"] = "fractional"):
+        """
+        Site symmetry operations of the first point, with translations removed.
+
+        Args:
+            space: "fractional" returns the rotations as stored (fractional
+                coordinates). "euclidean" returns them in Cartesian space, which
+                differs only for trigonal/hexagonal groups (``self.euclidean``),
+                where a 3-fold rotation such as (-y, x-y, z) is not orthogonal.
+                Use "euclidean" when comparing against a molecule's point group.
+
+        Returns:
+            list of pymatgen SymmOp objects
+        """
+        if space == "fractional":
+            ops = self.symmetry[0]
+        elif space == "euclidean":
+            ops = self.get_site_symm_ops()
+        else:
+            raise ValueError(f"space must be 'fractional' or 'euclidean', got {space!r}")
+        return [SymmOp.from_rotation_and_translation(op.rotation_matrix, [0, 0, 0]) for op in ops]
 
     def get_site_symmetry_object(self, idx=0):
         ops = self.get_site_symm_ops(idx)#; print(self.number, self.index, self.letter)

tests/test_molecule.py

@@ -307,5 +307,25 @@ def test_special_sites(self):
         struc.from_random(3, 61, ["Benzene"], [4])
         assert struc.valid
 
+    def test_special_sites_in_hexagonal_setting(self):
+        """Benzene (D6h) on the special positions of R-3 (148, hexagonal axes)."""
+        cases = [
+            (9, [["9d"]]),  # Z' = 1/2: inversion centre (-1)
+            (9, None),
+            (6, [["6c"]]),  # Z' = 1/3: three-fold axis (3)
+            (6, None),
+            (3, [["3a"]]),  # Z' = 1/6: rotoinversion axis (-3)
+            (3, None),
+        ]
+        for num_mols, site in cases:
+            with self.subTest(site=site):
+                struc = pyxtal(molecular=True)
+                struc.from_random(3, 148, ["Benzene"], [num_mols], sites=site)
+                assert struc.valid
+                if site is not None:
+                    assert [[s.wp.get_label() for s in struc.mol_sites]] == site
+                sga = SpacegroupAnalyzer(struc.to_pymatgen(), symprec=1e-2)
+                assert sga.get_space_group_number() == 148
+
 if __name__ == "__main__":
     unittest.main()