Recover point group by group multiplication

.github/workflows/check.yaml

@@ -20,4 +20,5 @@ jobs:
     runs-on: ubuntu-latest
     steps:
     - uses: actions/checkout@v4
+    - uses: Swatinem/rust-cache@v2
     - uses: EmbarkStudios/cargo-deny-action@v1

.github/workflows/ci-python.yaml

@@ -51,7 +51,7 @@ jobs:
           pytest -v
       - name: pytest
         if: ${{ !startsWith(matrix.target, 'x86') && matrix.target != 'ppc64' }}
-        uses: uraimo/run-on-arch-action@v2.5.0
+        uses: uraimo/run-on-arch-action@v2
         with:
           arch: ${{ matrix.target }}
           distro: ubuntu22.04

moyo/src/base/error.rs

@@ -16,6 +16,8 @@ pub enum MoyoError {
     BravaisGroupSearchError,
     #[error("Primitive symmetry search failed")]
     PrimitiveSymmetrySearchError,
+    #[error("Geometric crystal class identification failed")]
+    GeometricCrystalClassIdentificationError,
     #[error("Arithmetic crystal class identification failed")]
     ArithmeticCrystalClassIdentificationError,
     #[error("Space group type identification failed")]

moyo/src/identify/point_group.rs

@@ -1,6 +1,7 @@
 use std::cmp::Ordering;
 
 use itertools::Itertools;
+use log::debug;
 use nalgebra::{Dyn, Matrix3, OMatrix, OVector, U9};
 
 use crate::base::{MoyoError, Rotation, UnimodularLinear};
@@ -23,7 +24,7 @@ impl PointGroup {
     /// Assume the rotations are given in the (reduced) primitive basis
     pub fn new(prim_rotations: &Vec<Rotation>) -> Result<Self, MoyoError> {
         let rotation_types = prim_rotations.iter().map(identify_rotation_type).collect();
-        let geometric_crystal_class = identify_geometric_crystal_class(&rotation_types);
+        let geometric_crystal_class = identify_geometric_crystal_class(&rotation_types)?;
 
         let crystal_system = CrystalSystem::from_geometric_crystal_class(geometric_crystal_class);
         match crystal_system {
@@ -346,7 +347,9 @@ fn identify_rotation_type(rotation: &Rotation) -> RotationType {
 }
 
 /// Use look up table in Table 6 of https://arxiv.org/pdf/1808.01590.pdf
-fn identify_geometric_crystal_class(rotation_types: &Vec<RotationType>) -> GeometricCrystalClass {
+fn identify_geometric_crystal_class(
+    rotation_types: &Vec<RotationType>,
+) -> Result<GeometricCrystalClass, MoyoError> {
     // count RotationTypes in point_group
     let mut rotation_types_count = [0; 10];
     for rotation_type in rotation_types {
@@ -365,46 +368,52 @@ fn identify_geometric_crystal_class(rotation_types: &Vec<RotationType>) -> Geome
     }
     match rotation_types_count {
         // Triclinic
-        [0, 0, 0, 0, 0, 1, 0, 0, 0, 0] => GeometricCrystalClass::C1,
-        [0, 0, 0, 0, 1, 1, 0, 0, 0, 0] => GeometricCrystalClass::Ci,
+        [0, 0, 0, 0, 0, 1, 0, 0, 0, 0] => Ok(GeometricCrystalClass::C1),
+        [0, 0, 0, 0, 1, 1, 0, 0, 0, 0] => Ok(GeometricCrystalClass::Ci),
         // Monoclinic
-        [0, 0, 0, 0, 0, 1, 1, 0, 0, 0] => GeometricCrystalClass::C2,
-        [0, 0, 0, 1, 0, 1, 0, 0, 0, 0] => GeometricCrystalClass::C1h,
-        [0, 0, 0, 1, 1, 1, 1, 0, 0, 0] => GeometricCrystalClass::C2h,
+        [0, 0, 0, 0, 0, 1, 1, 0, 0, 0] => Ok(GeometricCrystalClass::C2),
+        [0, 0, 0, 1, 0, 1, 0, 0, 0, 0] => Ok(GeometricCrystalClass::C1h),
+        [0, 0, 0, 1, 1, 1, 1, 0, 0, 0] => Ok(GeometricCrystalClass::C2h),
         // Orthorhombic
-        [0, 0, 0, 0, 0, 1, 3, 0, 0, 0] => GeometricCrystalClass::D2,
-        [0, 0, 0, 2, 0, 1, 1, 0, 0, 0] => GeometricCrystalClass::C2v,
-        [0, 0, 0, 3, 1, 1, 3, 0, 0, 0] => GeometricCrystalClass::D2h,
+        [0, 0, 0, 0, 0, 1, 3, 0, 0, 0] => Ok(GeometricCrystalClass::D2),
+        [0, 0, 0, 2, 0, 1, 1, 0, 0, 0] => Ok(GeometricCrystalClass::C2v),
+        [0, 0, 0, 3, 1, 1, 3, 0, 0, 0] => Ok(GeometricCrystalClass::D2h),
         // Tetragonal
-        [0, 0, 0, 0, 0, 1, 1, 0, 2, 0] => GeometricCrystalClass::C4,
-        [0, 2, 0, 0, 0, 1, 1, 0, 0, 0] => GeometricCrystalClass::S4,
-        [0, 2, 0, 1, 1, 1, 1, 0, 2, 0] => GeometricCrystalClass::C4h,
-        [0, 0, 0, 0, 0, 1, 5, 0, 2, 0] => GeometricCrystalClass::D4,
-        [0, 0, 0, 4, 0, 1, 1, 0, 2, 0] => GeometricCrystalClass::C4v,
-        [0, 2, 0, 2, 0, 1, 3, 0, 0, 0] => GeometricCrystalClass::D2d,
-        [0, 2, 0, 5, 1, 1, 5, 0, 2, 0] => GeometricCrystalClass::D4h,
+        [0, 0, 0, 0, 0, 1, 1, 0, 2, 0] => Ok(GeometricCrystalClass::C4),
+        [0, 2, 0, 0, 0, 1, 1, 0, 0, 0] => Ok(GeometricCrystalClass::S4),
+        [0, 2, 0, 1, 1, 1, 1, 0, 2, 0] => Ok(GeometricCrystalClass::C4h),
+        [0, 0, 0, 0, 0, 1, 5, 0, 2, 0] => Ok(GeometricCrystalClass::D4),
+        [0, 0, 0, 4, 0, 1, 1, 0, 2, 0] => Ok(GeometricCrystalClass::C4v),
+        [0, 2, 0, 2, 0, 1, 3, 0, 0, 0] => Ok(GeometricCrystalClass::D2d),
+        [0, 2, 0, 5, 1, 1, 5, 0, 2, 0] => Ok(GeometricCrystalClass::D4h),
         // Trigonal
-        [0, 0, 0, 0, 0, 1, 0, 2, 0, 0] => GeometricCrystalClass::C3,
-        [0, 0, 2, 0, 1, 1, 0, 2, 0, 0] => GeometricCrystalClass::C3i,
-        [0, 0, 0, 0, 0, 1, 3, 2, 0, 0] => GeometricCrystalClass::D3,
-        [0, 0, 0, 3, 0, 1, 0, 2, 0, 0] => GeometricCrystalClass::C3v,
-        [0, 0, 2, 3, 1, 1, 3, 2, 0, 0] => GeometricCrystalClass::D3d,
+        [0, 0, 0, 0, 0, 1, 0, 2, 0, 0] => Ok(GeometricCrystalClass::C3),
+        [0, 0, 2, 0, 1, 1, 0, 2, 0, 0] => Ok(GeometricCrystalClass::C3i),
+        [0, 0, 0, 0, 0, 1, 3, 2, 0, 0] => Ok(GeometricCrystalClass::D3),
+        [0, 0, 0, 3, 0, 1, 0, 2, 0, 0] => Ok(GeometricCrystalClass::C3v),
+        [0, 0, 2, 3, 1, 1, 3, 2, 0, 0] => Ok(GeometricCrystalClass::D3d),
         // Hexagonal
-        [0, 0, 0, 0, 0, 1, 1, 2, 0, 2] => GeometricCrystalClass::C6,
-        [2, 0, 0, 1, 0, 1, 0, 2, 0, 0] => GeometricCrystalClass::C3h,
-        [2, 0, 2, 1, 1, 1, 1, 2, 0, 2] => GeometricCrystalClass::C6h,
-        [0, 0, 0, 0, 0, 1, 7, 2, 0, 2] => GeometricCrystalClass::D6,
-        [0, 0, 0, 6, 0, 1, 1, 2, 0, 2] => GeometricCrystalClass::C6v,
-        [2, 0, 0, 4, 0, 1, 3, 2, 0, 0] => GeometricCrystalClass::D3h,
-        [2, 0, 2, 7, 1, 1, 7, 2, 0, 2] => GeometricCrystalClass::D6h,
+        [0, 0, 0, 0, 0, 1, 1, 2, 0, 2] => Ok(GeometricCrystalClass::C6),
+        [2, 0, 0, 1, 0, 1, 0, 2, 0, 0] => Ok(GeometricCrystalClass::C3h),
+        [2, 0, 2, 1, 1, 1, 1, 2, 0, 2] => Ok(GeometricCrystalClass::C6h),
+        [0, 0, 0, 0, 0, 1, 7, 2, 0, 2] => Ok(GeometricCrystalClass::D6),
+        [0, 0, 0, 6, 0, 1, 1, 2, 0, 2] => Ok(GeometricCrystalClass::C6v),
+        [2, 0, 0, 4, 0, 1, 3, 2, 0, 0] => Ok(GeometricCrystalClass::D3h),
+        [2, 0, 2, 7, 1, 1, 7, 2, 0, 2] => Ok(GeometricCrystalClass::D6h),
         // Cubic
-        [0, 0, 0, 0, 0, 1, 3, 8, 0, 0] => GeometricCrystalClass::T,
-        [0, 0, 8, 3, 1, 1, 3, 8, 0, 0] => GeometricCrystalClass::Th,
-        [0, 0, 0, 0, 0, 1, 9, 8, 6, 0] => GeometricCrystalClass::O,
-        [0, 6, 0, 6, 0, 1, 3, 8, 0, 0] => GeometricCrystalClass::Td,
-        [0, 6, 8, 9, 1, 1, 9, 8, 6, 0] => GeometricCrystalClass::Oh,
+        [0, 0, 0, 0, 0, 1, 3, 8, 0, 0] => Ok(GeometricCrystalClass::T),
+        [0, 0, 8, 3, 1, 1, 3, 8, 0, 0] => Ok(GeometricCrystalClass::Th),
+        [0, 0, 0, 0, 0, 1, 9, 8, 6, 0] => Ok(GeometricCrystalClass::O),
+        [0, 6, 0, 6, 0, 1, 3, 8, 0, 0] => Ok(GeometricCrystalClass::Td),
+        [0, 6, 8, 9, 1, 1, 9, 8, 6, 0] => Ok(GeometricCrystalClass::Oh),
         // Unknown
-        _ => unreachable!("Unknown point group"),
+        _ => {
+            debug!(
+                "Unknown geometric crystal class: {:?}",
+                rotation_types_count
+            );
+            Err(MoyoError::GeometricCrystalClassIdentificationError)
+        }
     }
 }
 

moyo/src/search/symmetry_search.rs

@@ -1,4 +1,5 @@
 use itertools::iproduct;
+use std::collections::{HashMap, HashSet, VecDeque};
 
 use log::{debug, warn};
 use nalgebra::{Matrix3, Vector3};
@@ -9,7 +10,7 @@ use super::solve::{
 };
 use crate::base::{
     traverse, AngleTolerance, Cell, Lattice, MoyoError, Operations, Permutation, Position,
-    Rotation, EPS,
+    Rotation, Translation, EPS,
 };
 
 #[derive(Debug)]
@@ -66,9 +67,7 @@ impl PrimitiveSymmetrySearch {
         assert!(!symmetries_tmp.is_empty());
 
         // Purify symmetry operations by permutations
-        let mut rotations = vec![];
-        let mut translations = vec![];
-        let mut permutations = vec![];
+        let mut translations_and_permutations = HashMap::new();
         for (rotation, rough_translation, permutation) in symmetries_tmp.iter() {
             let (translation, distance) = symmetrize_translation_from_permutation(
                 primitive_cell,
@@ -77,15 +76,51 @@ impl PrimitiveSymmetrySearch {
                 rough_translation,
             );
             if distance < symprec {
-                rotations.push(*rotation);
-                translations.push(translation);
-                permutations.push(permutation.clone());
+                translations_and_permutations.insert(*rotation, (translation, permutation.clone()));
             }
         }
-        if rotations.is_empty() {
+        if translations_and_permutations.is_empty() {
             return Err(MoyoError::PrimitiveSymmetrySearchError);
         }
 
+        // Recover operations by group multiplication
+        let mut queue = VecDeque::new();
+        let mut visited = HashSet::new();
+        let mut rotations = vec![];
+        let mut translations = vec![];
+        let mut permutations = vec![];
+        queue.push_back((
+            Rotation::identity(),
+            Translation::zeros(),
+            Permutation::identity(primitive_cell.num_atoms()),
+        ));
+
+        while !queue.is_empty() {
+            let (rotation_lhs, translation_lhs, permutation_lhs) = queue.pop_front().unwrap();
+            if visited.contains(&rotation_lhs) {
+                continue;
+            }
+            visited.insert(rotation_lhs);
+            rotations.push(rotation_lhs);
+            translations.push(translation_lhs);
+            permutations.push(permutation_lhs.clone());
+
+            for (&rotation_rhs, (translation_rhs, permutation_rhs)) in
+                translations_and_permutations.iter()
+            {
+                let new_rotation = rotation_lhs * rotation_rhs;
+                let new_translation = (rotation_lhs.map(|e| e as f64) * translation_rhs
+                    + translation_lhs)
+                    .map(|e| e - e.floor());
+                let new_permutation = permutation_lhs.clone() * permutation_rhs.clone();
+                queue.push_back((new_rotation, new_translation, new_permutation));
+            }
+        }
+        if rotations.len() != translations_and_permutations.len() {
+            warn!("Found operations do not form a group. symprec and angle_tolerance may be too large.");
+        }
+
+        debug!("Order of point group: {}", rotations.len());
         Ok(Self {
             operations: Operations::new(rotations, translations),
             permutations,
@@ -191,7 +226,7 @@ fn search_bravais_group(
         return Err(MoyoError::BravaisGroupSearchError);
     }
 
-    // Complement rotations by group multiplication
+    // Recover rotations by group multiplication
     let complemented_rotations = traverse(&rotations);
     if complemented_rotations.len() != rotations.len() {
         warn!("Found automorphisms for the lattice do not form a group. symprec and angle_tolerance may be too large.");

moyo/tests/test_moyo_dataset.rs

@@ -461,3 +461,22 @@ fn test_with_mp_1197586() {
     assert_eq!(dataset.hall_number, 488);
     assert_eq!(dataset.num_operations(), 24);
 }
+
+#[test]
+fn test_with_mp_1185639() {
+    // https://next-gen.materialsproject.org/materials/mp-1185639
+    let path = Path::new("tests/assets/mp-1185639.json");
+    let cell: Cell = serde_json::from_str(&fs::read_to_string(&path).unwrap()).unwrap();
+
+    let symprec = 1e-4;
+    let angle_tolerance = AngleTolerance::Default;
+    let setting = Setting::Standard;
+
+    let dataset = assert_dataset(&cell, symprec, angle_tolerance, setting);
+    assert_dataset(&dataset.std_cell, symprec, angle_tolerance, setting);
+    assert_dataset(&dataset.prim_std_cell, symprec, angle_tolerance, setting);
+
+    assert_eq!(dataset.number, 187); // P-6m2
+    assert_eq!(dataset.hall_number, 481);
+    assert_eq!(dataset.num_operations(), 12);
+}