Iteratively change tolerance

moyo/src/base/error.rs

@@ -1,21 +1,23 @@
 use thiserror::Error;
 
-#[derive(Error, Debug)]
+#[derive(Error, Debug, PartialEq, Eq, Clone, Copy)]
 pub enum MoyoError {
     #[error("Minkowski reduction failed")]
     MinkowskiReductionError,
     #[error("Niggli reduction failed")]
     NiggliReductionError,
     #[error("Delaunay reduction failed")]
     DelaunayReductionError,
-    #[error("Too small symprec")]
-    TooSmallSymprecError,
+    #[error("Too small tolerance")]
+    TooSmallToleranceError,
+    #[error("Too large tolerance")]
+    TooLargeToleranceError,
     #[error("Primitive cell search failed")]
     PrimitiveCellError,
-    #[error("Bravais group search failed")]
-    BravaisGroupSearchError,
     #[error("Primitive symmetry search failed")]
     PrimitiveSymmetrySearchError,
+    #[error("Bravais group search failed")]
+    BravaisGroupSearchError,
     #[error("Geometric crystal class identification failed")]
     GeometricCrystalClassIdentificationError,
     #[error("Arithmetic crystal class identification failed")]

moyo/src/base/transformation.rs

@@ -43,6 +43,10 @@ impl UnimodularTransformation {
         Self::new(linear, OriginShift::zeros())
     }
 
+    pub fn from_origin_shift(origin_shift: OriginShift) -> Self {
+        Self::new(UnimodularLinear::identity(), origin_shift)
+    }
+
     pub fn linear_as_f64(&self) -> Matrix3<f64> {
         self.linear.map(|e| e as f64)
     }
@@ -111,6 +115,10 @@ impl Transformation {
         Self::new(linear, OriginShift::zeros())
     }
 
+    pub fn from_origin_shift(origin_shift: OriginShift) -> Self {
+        Self::new(Linear::identity(), origin_shift)
+    }
+
     pub fn linear_as_f64(&self) -> Matrix3<f64> {
         self.linear.map(|e| e as f64)
     }

moyo/src/identify/point_group.rs

@@ -25,6 +25,7 @@ impl PointGroup {
     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)?;
+        debug!("Geometric crystal class: {:?}", geometric_crystal_class);
 
         let crystal_system = CrystalSystem::from_geometric_crystal_class(geometric_crystal_class);
         match crystal_system {

moyo/src/identify/space_group.rs

@@ -1,5 +1,6 @@
 use std::collections::HashMap;
 
+use log::debug;
 use nalgebra::{Dyn, Matrix3, OMatrix, OVector, Vector3, U3};
 
 use super::point_group::PointGroup;
@@ -28,6 +29,10 @@ impl SpaceGroup {
     ) -> Result<Self, MoyoError> {
         // point_group.trans_mat: self -> primitive
         let point_group = PointGroup::new(&prim_operations.rotations)?;
+        debug!(
+            "Arithmetic crystal class: No. {}",
+            point_group.arithmetic_number
+        );
 
         for hall_number in setting.hall_numbers() {
             let entry = hall_symbol_entry(hall_number);
@@ -44,6 +49,10 @@ impl SpaceGroup {
                 if let Some(origin_shift) =
                     match_origin_shift(prim_operations, &trans_mat, &db_prim_generators, epsilon)
                 {
+                    debug!(
+                        "Matched with Hall number {} (No. {})",
+                        hall_number, entry.number
+                    );
                     return Ok(Self {
                         number: entry.number,
                         hall_number,
@@ -198,6 +207,7 @@ fn match_origin_shift(
     match solve_mod1(&a, &b, epsilon) {
         Some(s) => {
             let origin_shift = (trans_mat.map(|e| e as f64) * s).map(|e| e % 1.);
+
             Some(origin_shift)
         }
         None => None,

moyo/src/lib.rs

@@ -9,6 +9,7 @@ pub mod math;
 pub mod search;
 pub mod symmetrize;
 
+use log::warn;
 use nalgebra::Matrix3;
 
 use crate::base::{AngleTolerance, Cell, MoyoError, Operations, OriginShift, Transformation};
@@ -64,8 +65,8 @@ pub struct MoyoDataset {
     // ------------------------------------------------------------------------
     // Final parameters
     // ------------------------------------------------------------------------
-    // pub symprec: f64,
-    // pub angle_tolerance: AngleTolerance,
+    pub symprec: f64,
+    pub angle_tolerance: AngleTolerance,
 }
 
 impl MoyoDataset {
@@ -75,10 +76,8 @@ impl MoyoDataset {
         angle_tolerance: AngleTolerance,
         setting: Setting,
     ) -> Result<Self, MoyoError> {
-        // Symmetry search
-        let prim_cell = PrimitiveCell::new(cell, symprec)?;
-        let symmetry_search =
-            PrimitiveSymmetrySearch::new(&prim_cell.cell, symprec, angle_tolerance)?;
+        let (prim_cell, symmetry_search, symprec, angle_tolerance) =
+            iterative_symmetry_search(cell, symprec, angle_tolerance)?;
         let operations = operations_in_cell(&prim_cell, &symmetry_search.operations);
 
         // Space-group type identification
@@ -143,6 +142,9 @@ impl MoyoDataset {
                 .iter()
                 .map(|w| w.site_symmetry.to_string())
                 .collect(),
+            // Final parameters
+            symprec,
+            angle_tolerance,
         })
     }
 
@@ -151,6 +153,119 @@ impl MoyoDataset {
     }
 }
 
+const MAX_SYMMETRY_SEARCH_TRIALS: usize = 16;
+const INITIAL_SYMMETRY_SEARCH_STRIDE: f64 = 2.0;
+
+fn iterative_symmetry_search(
+    cell: &Cell,
+    symprec: f64,
+    angle_tolerance: AngleTolerance,
+) -> Result<(PrimitiveCell, PrimitiveSymmetrySearch, f64, AngleTolerance), MoyoError> {
+    let mut symprec = symprec;
+    let mut angle_tolerance = angle_tolerance;
+    let mut stride: f64 = INITIAL_SYMMETRY_SEARCH_STRIDE;
+    let mut prev_error: Option<MoyoError> = None;
+
+    fn _increase_tolerance(
+        symprec: f64,
+        angle_tolerance: AngleTolerance,
+        stride: f64,
+    ) -> (f64, AngleTolerance) {
+        let symprec = symprec * stride;
+        let angle_tolerance = if let AngleTolerance::Radian(angle) = angle_tolerance {
+            AngleTolerance::Radian(angle * stride)
+        } else {
+            AngleTolerance::Default
+        };
+        warn!(
+            "Increase tolerance to symprec={}, angle_tolerance={:?}",
+            symprec, angle_tolerance
+        );
+        (symprec, angle_tolerance)
+    }
+
+    fn _reduce_tolerance(
+        symprec: f64,
+        angle_tolerance: AngleTolerance,
+        stride: f64,
+    ) -> (f64, AngleTolerance) {
+        let symprec = symprec / stride;
+        let angle_tolerance = if let AngleTolerance::Radian(angle) = angle_tolerance {
+            AngleTolerance::Radian(angle / stride)
+        } else {
+            AngleTolerance::Default
+        };
+        warn!(
+            "Reduce tolerance to symprec={}, angle_tolerance={:?}",
+            symprec, angle_tolerance
+        );
+        (symprec, angle_tolerance)
+    }
+
+    fn _update_stride(
+        err: MoyoError,
+        prev_error: Option<MoyoError>,
+        stride: f64,
+    ) -> (Option<MoyoError>, f64) {
+        let stride = if !prev_error.is_none() && prev_error != Some(err) {
+            stride.sqrt()
+        } else {
+            stride
+        };
+        (Some(err), stride)
+    }
+
+    fn handle_tolerance_error(
+        err: MoyoError,
+        prev_error: Option<MoyoError>,
+        stride: f64,
+        symprec: f64,
+        angle_tolerance: AngleTolerance,
+    ) -> (Option<MoyoError>, f64, f64, AngleTolerance) {
+        let (new_error, new_stride) = _update_stride(err, prev_error, stride);
+        let (new_symprec, new_angle_tolerance) = match err {
+            MoyoError::TooSmallToleranceError => {
+                _increase_tolerance(symprec, angle_tolerance, new_stride)
+            }
+            MoyoError::TooLargeToleranceError => {
+                _reduce_tolerance(symprec, angle_tolerance, new_stride)
+            }
+            _ => (symprec, angle_tolerance),
+        };
+        (new_error, new_stride, new_symprec, new_angle_tolerance)
+    }
+
+    for _ in 0..MAX_SYMMETRY_SEARCH_TRIALS {
+        match PrimitiveCell::new(cell, symprec) {
+            Ok(prim_cell) => {
+                match PrimitiveSymmetrySearch::new(&prim_cell.cell, symprec, angle_tolerance) {
+                    Ok(symmetry_search) => {
+                        return Ok((prim_cell, symmetry_search, symprec, angle_tolerance));
+                    }
+                    Err(err @ MoyoError::TooSmallToleranceError)
+                    | Err(err @ MoyoError::TooLargeToleranceError) => {
+                        (prev_error, stride, symprec, angle_tolerance) = handle_tolerance_error(
+                            err,
+                            prev_error,
+                            stride,
+                            symprec,
+                            angle_tolerance,
+                        );
+                    }
+                    Err(err) => return Err(err),
+                }
+            }
+            Err(err @ MoyoError::TooSmallToleranceError)
+            | Err(err @ MoyoError::TooLargeToleranceError) => {
+                (prev_error, stride, symprec, angle_tolerance) =
+                    handle_tolerance_error(err, prev_error, stride, symprec, angle_tolerance);
+            }
+            Err(err) => return Err(err),
+        }
+    }
+    Err(MoyoError::PrimitiveSymmetrySearchError)
+}
+
 fn operations_in_cell(prim_cell: &PrimitiveCell, prim_operations: &Operations) -> Operations {
     let mut rotations = vec![];
     let mut translations = vec![];

moyo/src/search/primitive_cell.rs

@@ -1,5 +1,6 @@
 use std::collections::BTreeMap;
 
+use log::error;
 use nalgebra::{Dyn, Matrix3, OMatrix, Vector3, U3};
 
 use super::solve::{
@@ -44,7 +45,10 @@ impl PrimitiveCell {
             .fold(f64::INFINITY, f64::min);
         let rough_symprec = 2.0 * symprec;
         if rough_symprec > minimum_basis_norm / 2.0 {
-            return Err(MoyoError::TooSmallSymprecError);
+            error!(
+                "symprec is too large compared to the basis vectors. Consider reducing symprec."
+            );
+            return Err(MoyoError::TooLargeToleranceError);
         }
 
         // Try possible translations: overlap the `src`the site to the `dst`th site

moyo/src/search/symmetry_search.rs

@@ -1,7 +1,7 @@
 use itertools::iproduct;
-use std::collections::{HashSet, VecDeque};
+use std::collections::{HashMap, HashSet, VecDeque};
 
-use log::{debug, warn};
+use log::{debug, error, warn};
 use nalgebra::{Matrix3, Vector3};
 
 use super::solve::{
@@ -24,6 +24,8 @@ pub struct PrimitiveSymmetrySearch {
 impl PrimitiveSymmetrySearch {
     /// Return coset representatives of the space group w.r.t. its translation subgroup.
     /// Assume `primitive_cell` is a primitive cell and its basis vectors are Minkowski reduced.
+    /// Returned operations are guaranteed to form a group.
+    /// If the group closure and tolerance (symprec and angle_tolerance) are incompatible, the former is prioritized.
     /// Possible replacements for spglib/src/spacegroup.h::spa_search_spacegroup
     pub fn new(
         primitive_cell: &Cell,
@@ -34,7 +36,10 @@ impl PrimitiveSymmetrySearch {
         let minimum_basis_norm = primitive_cell.lattice.basis.column(0).norm();
         let rough_symprec = 2.0 * symprec;
         if rough_symprec > minimum_basis_norm / 2.0 {
-            return Err(MoyoError::TooSmallSymprecError);
+            error!(
+                "symprec is too large compared to the basis vectors. Consider reducing symprec."
+            );
+            return Err(MoyoError::TooLargeToleranceError);
         }
 
         // Search symmetry operations
@@ -80,7 +85,10 @@ impl PrimitiveSymmetrySearch {
             }
         }
         if operations_and_permutations.is_empty() {
-            return Err(MoyoError::PrimitiveSymmetrySearchError);
+            error!(
+                "No symmetry operations are found. Consider increasing symprec and angle_tolerance."
+            );
+            return Err(MoyoError::TooSmallToleranceError);
         }
 
         // Recover operations by group multiplication
@@ -94,7 +102,6 @@ impl PrimitiveSymmetrySearch {
             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) {
@@ -111,13 +118,39 @@ impl PrimitiveSymmetrySearch {
                 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());
+                    .map(|e| e - e.round());
                 let new_permutation = permutation_lhs.clone() * permutation_rhs.clone();
                 queue.push_back((new_rotation, new_translation, new_permutation));
             }
         }
         if rotations.len() != operations_and_permutations.len() {
-            warn!("Found operations do not form a group. symprec and angle_tolerance may be too large.");
+            warn!("Found operations do not form a group. Consider reducing symprec and angle_tolerance.");
+        }
+
+        // Check closure
+        let mut translations_map = HashMap::new();
+        for (rotation, translation) in rotations.iter().zip(translations.iter()) {
+            translations_map.insert(rotation.clone(), translation.clone());
+        }
+        let mut closed = true;
+        for (r1, t1) in rotations.iter().zip(translations.iter()) {
+            if !closed {
+                break;
+            }
+            for (r2, t2) in rotations.iter().zip(translations.iter()) {
+                // (r1, t1) * (r2, t2) = (r1 * r2, r1 * t2 + t1)
+                let r = r1 * r2;
+                let t = r1.map(|e| e as f64) * t2 + t1;
+                let diff = (translations_map[&r] - t).map(|e| e - e.round());
+                if primitive_cell.lattice.cartesian_coords(&diff).norm() > rough_symprec {
+                    closed = false;
+                    break;
+                }
+            }
+        }
+        if !closed {
+            error!("Some centering translations are missing. Consider increasing symprec and angle_tolerance.");
+            return Err(MoyoError::TooSmallToleranceError);
         }
 
         debug!("Order of point group: {}", rotations.len());
@@ -229,7 +262,7 @@ fn search_bravais_group(
     // 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.");
+        warn!("Found automorphisms for the lattice do not form a group. Consider reducing symprec and angle_tolerance.");
     }
     debug!("Order of Bravais group: {}", complemented_rotations.len());
     Ok(complemented_rotations)