WIP: conjugator for type 4

moyo/src/identify/magnetic_space_group.rs

@@ -1,11 +1,15 @@
 use std::collections::HashMap;
 
+use itertools::Itertools;
 use log::debug;
 
 use super::normalizer::integral_normalizer;
-use super::space_group::SpaceGroup;
+use super::point_group::iter_trans_mat_basis;
+use super::rotation_type::identify_rotation_type;
+use super::space_group::{match_origin_shift, SpaceGroup};
 use crate::base::{
-    MagneticOperations, MoyoError, Operations, Rotation, Translation, UnimodularTransformation,
+    project_rotations, MagneticOperations, MoyoError, Operations, Rotation, Translation,
+    UnimodularLinear, UnimodularTransformation,
 };
 use crate::data::{
     get_magnetic_space_group_type, hall_symbol_entry, magnetic_hall_symbol_entry, uni_number_range,
@@ -67,6 +71,7 @@ impl MagneticSpaceGroup {
                     new_transformation.transform_magnetic_operations(prim_mag_operations);
 
                 // TODO:
+
                 // debug!("Matched with UNI number {}", uni_number);
                 // return Ok(Self {
                 //     uni_number,
@@ -184,6 +189,61 @@ fn db_reference_space_group_primitive(entry: &MagneticHallSymbolEntry) -> (Opera
     (ref_prim_operations, ref_prim_generators)
 }
 
+/// Find a unimodular transformation that transforms (E, src_translation) to (E, dst_translation) while keeping `stabilized_prim_operations`, which are generated by `stabilized_prim_generators`.
+fn find_conjugator_type4(
+    stabilized_prim_generators: &Operations,
+    stabilized_prim_operations: &Operations,
+    src_translation: &Translation,
+    dst_translation: &Translation,
+    epsilon: f64,
+) -> Option<UnimodularTransformation> {
+    let stabilized_prim_rotations = project_rotations(stabilized_prim_operations);
+    let stabilized_prim_rotation_generators = project_rotations(stabilized_prim_generators);
+
+    let rotation_types = stabilized_prim_rotations
+        .iter()
+        .map(identify_rotation_type)
+        .collect::<Vec<_>>();
+    for trans_mat_basis in iter_trans_mat_basis(
+        stabilized_prim_rotations,
+        rotation_types,
+        stabilized_prim_rotation_generators,
+    ) {
+        // Search integer linear combination such that the transformation matrix is unimodular
+        // Consider coefficients in [-2, 2], which will be sufficient for Delaunay reduced basis
+        for comb in (0..trans_mat_basis.len())
+            .map(|_| -2..=2)
+            .multi_cartesian_product()
+        {
+            let mut prim_trans_mat = UnimodularLinear::zeros();
+            for (i, matrix) in trans_mat_basis.iter().enumerate() {
+                prim_trans_mat += comb[i] * matrix;
+            }
+            let det = prim_trans_mat.map(|e| e as f64).determinant().round() as i32;
+            if det < 0 {
+                prim_trans_mat *= -1;
+            }
+            if det == 1 {
+                // (P, p)^-1 (E, c_src) (P, p) = (P^-1, -P^-1 p) (P, p + c_src) = (E, P^-1 c_src) == (E, c_dst)
+                let diff = prim_trans_mat.map(|e| e as f64) * dst_translation - src_translation;
+                if !diff.iter().all(|e| (e - e.round()).abs() < epsilon) {
+                    continue;
+                }
+
+                if let Some(origin_shift) = match_origin_shift(
+                    stabilized_prim_operations,
+                    &prim_trans_mat,
+                    stabilized_prim_generators,
+                    epsilon,
+                ) {
+                    return Some(UnimodularTransformation::new(prim_trans_mat, origin_shift));
+                }
+            }
+        }
+    }
+    None
+}
+
 #[cfg(test)]
 mod tests {
     use rstest::rstest;