diff --git a/moyo/src/identify.rs b/moyo/src/identify.rs index 2e92d48..af9ffd0 100644 --- a/moyo/src/identify.rs +++ b/moyo/src/identify.rs @@ -4,4 +4,5 @@ mod point_group; mod rotation_type; mod space_group; +pub(super) use magnetic_space_group::MagneticSpaceGroup; pub(super) use space_group::SpaceGroup; diff --git a/moyo/src/identify/magnetic_space_group.rs b/moyo/src/identify/magnetic_space_group.rs index 0b3e492..142266f 100644 --- a/moyo/src/identify/magnetic_space_group.rs +++ b/moyo/src/identify/magnetic_space_group.rs @@ -1,13 +1,12 @@ -use std::collections::{HashMap, HashSet}; +use std::collections::HashMap; use log::debug; use nalgebra::{Dyn, OMatrix, OVector, U3}; -use super::normalizer::integral_normalizer; use super::space_group::{solve_mod1, SpaceGroup}; use crate::base::{ - project_rotations, traverse, MagneticOperations, MoyoError, Operations, Rotation, Translation, - UnimodularLinear, UnimodularTransformation, + MagneticOperations, MoyoError, Operations, Rotation, Translation, UnimodularLinear, + UnimodularTransformation, }; use crate::data::{ get_magnetic_space_group_type, uni_number_range, ConstructType, MagneticHallSymbol, Setting, @@ -45,33 +44,41 @@ impl MagneticSpaceGroup { continue; } + if construct_type == ConstructType::Type1 || construct_type == ConstructType::Type2 { + // No need to further check the magnetic operations + return Ok(Self { + uni_number, + transformation: std_ref_spg.transformation, + }); + } + + // TODO: + let mhs = MagneticHallSymbol::from_uni_number(uni_number) .ok_or(MoyoError::MagneticSpaceGroupTypeIdentificationError)?; let db_prim_mag_generators = mhs.primitive_generators(); let db_ref_prim_generators = Self::get_db_reference_space_group_primitive_generators(&mhs, &construct_type); - // The correction transformation matrices keep the point group of `tmp_prim_mag_operations` + // The correction transformations keep the reference space group of `tmp_prim_mag_operations` // TODO: precompute the correction transformation matrices - let correction_transformation_matrices = - Self::correction_transformation_matrices(&db_ref_prim_generators, &construct_type) - .ok_or(MoyoError::MagneticSpaceGroupTypeIdentificationError)?; - for corr_trans_mat in correction_transformation_matrices { - // (trans_mat, origin_shift): primitive input -> primitive DB - let trans_mat = std_ref_spg.transformation.linear * corr_trans_mat; - if let Some(origin_shift) = Self::match_origin_shift( - prim_mag_operations, - &trans_mat, - &db_prim_mag_generators, - epsilon, - ) { - debug!("Matched with UNI number {}", uni_number); - return Ok(Self { - uni_number, - transformation: UnimodularTransformation::new(trans_mat, origin_shift), - }); - } - } + let correction_transformations = todo!(); + // for corr_trans_mat in correction_transformation_matrices { + // // (trans_mat, origin_shift): primitive input -> primitive DB + // let trans_mat = std_ref_spg.transformation.linear * corr_trans_mat; + // if let Some(origin_shift) = Self::match_origin_shift( + // prim_mag_operations, + // &trans_mat, + // &db_prim_mag_generators, + // epsilon, + // ) { + // debug!("Matched with UNI number {}", uni_number); + // return Ok(Self { + // uni_number, + // transformation: UnimodularTransformation::new(trans_mat, origin_shift), + // }); + // } + // } } Err(MoyoError::MagneticSpaceGroupTypeIdentificationError) } @@ -110,27 +117,6 @@ impl MagneticSpaceGroup { } } - fn correction_transformation_matrices( - db_ref_prim_generators: &Operations, - construct_type: &ConstructType, - ) -> Option> { - // Family space group (FSG): F(M) - // Maximal space subgroup (XSG): D(M) - match construct_type { - ConstructType::Type1 | ConstructType::Type2 => Some(vec![UnimodularLinear::identity()]), - ConstructType::Type3 | ConstructType::Type4 => { - // For type-III, try to map D(M) into D(M_std) while keeping F(M) = F(M_std) - // For type-IV, try to map F(M) into F(M_std) while keeping D(M) = D(M_std) - let db_ref_prim_rotation_generators = project_rotations(&db_ref_prim_generators); - let db_ref_prim_rotations = traverse(&db_ref_prim_rotation_generators); - Some(integral_normalizer( - &db_ref_prim_rotations, - &db_ref_prim_rotation_generators, - )) - } - } - } - /// Search for origin_shift such that (trans_mat, origin_shift) transforms `prim_mag_operations` into /// TODO: unify with identify/space_group.rs::match_origin_shift fn match_origin_shift( @@ -179,8 +165,9 @@ fn identify_reference_space_group( prim_mag_operations: &MagneticOperations, epsilon: f64, ) -> Option<(Operations, ConstructType)> { - let prim_xsg = maximal_space_subgroup_from_magnetic_space_group(prim_mag_operations); - let fsg = family_space_group_from_magnetic_space_group(prim_mag_operations); + let prim_xsg = primitive_maximal_space_subgroup_from_magnetic_space_group(prim_mag_operations); + let (fsg, is_type2) = + family_space_group_from_magnetic_space_group(prim_mag_operations, epsilon); if (prim_mag_operations.len() % prim_xsg.len() != 0) || (prim_mag_operations.len() % fsg.len() != 0) @@ -189,25 +176,12 @@ fn identify_reference_space_group( return None; } - let identity = Rotation::identity(); - let has_pure_time_reversal = prim_mag_operations.iter().any(|mops| { - mops.time_reversal - && mops.operation.rotation == identity - && mops - .operation - .translation - .iter() - .all(|e| (e - e.round()).abs() < epsilon) - }); - - let construct_type = match ( - prim_mag_operations.len() / prim_xsg.len(), - has_pure_time_reversal, - ) { + let construct_type = match (prim_mag_operations.len() / prim_xsg.len(), is_type2) { (1, false) => ConstructType::Type1, (2, true) => ConstructType::Type2, (2, false) => { // Find coset representatives of MSG/XSG + let identity = Rotation::identity(); if prim_mag_operations .iter() .any(|mops| mops.time_reversal && mops.operation.rotation == identity) @@ -239,10 +213,11 @@ fn identify_reference_space_group( } /// XSG: take only operations without time-reversal -fn maximal_space_subgroup_from_magnetic_space_group( +fn primitive_maximal_space_subgroup_from_magnetic_space_group( prim_mag_operations: &MagneticOperations, ) -> Operations { let mut xsg = vec![]; + for mops in prim_mag_operations { if mops.time_reversal { continue; @@ -256,12 +231,25 @@ fn maximal_space_subgroup_from_magnetic_space_group( /// Returned operations may contain duplicated rotation parts (for type-IV). fn family_space_group_from_magnetic_space_group( prim_mag_operations: &MagneticOperations, -) -> Operations { + epsilon: f64, +) -> (Operations, bool) { let mut fsg = vec![]; + let mut hm_translation = HashMap::new(); + let mut is_type2 = false; + for mops in prim_mag_operations { + if let Some(&other_translation) = hm_translation.get(&mops.operation.rotation) { + let diff: Translation = mops.operation.translation - other_translation; + if diff.iter().all(|e| (e - e.round()).abs() < epsilon) { + is_type2 = true; + continue; + } + } + fsg.push(mops.operation.clone()); + hm_translation.insert(mops.operation.rotation.clone(), mops.operation.translation); } - fsg + (fsg, is_type2) } #[cfg(test)] @@ -284,9 +272,9 @@ mod tests { } #[rstest] - #[case(2, ConstructType::Type2, 2, 1, 2)] + #[case(2, ConstructType::Type2, 2, 1, 1)] #[case(1594, ConstructType::Type1, 48, 48, 48)] - #[case(1595, ConstructType::Type2, 96, 48, 96)] + #[case(1595, ConstructType::Type2, 96, 48, 48)] #[case(1596, ConstructType::Type3, 48, 24, 48)] #[case(1599, ConstructType::Type4, 96, 48, 96)] // -P 4 2 3 1abc' (UNI No. 1599) fn test_xsg_and_fsg( @@ -299,20 +287,22 @@ mod tests { let prim_mag_operations = get_prim_mag_operations(uni_number); assert_eq!(prim_mag_operations.len(), order_msg); - let xsg = maximal_space_subgroup_from_magnetic_space_group(&prim_mag_operations); + let xsg = primitive_maximal_space_subgroup_from_magnetic_space_group(&prim_mag_operations); assert_eq!(xsg.len(), order_xsg); - let fsg = family_space_group_from_magnetic_space_group(&prim_mag_operations); + let epsilon = 1e-8; + let (fsg, _) = family_space_group_from_magnetic_space_group(&prim_mag_operations, epsilon); assert_eq!(fsg.len(), order_fsg); let (_, construct_type_actual) = - identify_reference_space_group(&prim_mag_operations, 1e-8).unwrap(); + identify_reference_space_group(&prim_mag_operations, epsilon).unwrap(); assert_eq!(construct_type_actual, construct_type); } #[test_with_log] fn test_identify_magnetic_space_group() { - for uni_number in 1..=NUM_MAGNETIC_SPACE_GROUP_TYPES { + // for uni_number in 1..=NUM_MAGNETIC_SPACE_GROUP_TYPES { + for uni_number in 3..=NUM_MAGNETIC_SPACE_GROUP_TYPES { dbg!(uni_number); let prim_mag_operations = get_prim_mag_operations(uni_number as UNINumber); let magnetic_space_group = MagneticSpaceGroup::new(&prim_mag_operations, 1e-8).unwrap(); diff --git a/moyo/src/identify/normalizer.rs b/moyo/src/identify/normalizer.rs index c3c5932..17c6e3f 100644 --- a/moyo/src/identify/normalizer.rs +++ b/moyo/src/identify/normalizer.rs @@ -1,50 +1,29 @@ use itertools::Itertools; -use once_cell::sync::Lazy; -use std::collections::HashSet; use super::rotation_type::identify_rotation_type; -use crate::base::{traverse, Rotations, UnimodularLinear}; -use crate::data::{ArithmeticNumber, PointGroupRepresentative}; +use super::space_group::match_origin_shift; +use crate::base::{project_rotations, Operations, UnimodularLinear, UnimodularTransformation}; use crate::math::sylvester3; /// Return integral normalizer of the point group representative up to its centralizer. -pub fn get_point_group_normalizer( - arithmetic_number: ArithmeticNumber, -) -> Option> { - POINT_GROUP_NORMALIZERS - .get((arithmetic_number - 1) as usize) - .cloned() -} /// Integral normalizers for all arithmetic point groups up to their centralizers. -static POINT_GROUP_NORMALIZERS: Lazy>> = Lazy::new(|| { - (1..=73) - .map(|arithmetic_number| { - let point_group_db = - PointGroupRepresentative::from_arithmetic_crystal_class(arithmetic_number); - let prim_generators = point_group_db.primitive_generators(); - let prim_rotations = traverse(&prim_generators); - let mut prim_rotations_set = HashSet::new(); - for rotation in prim_rotations.iter() { - prim_rotations_set.insert(rotation.clone()); - } - integral_normalizer(&prim_rotations, &prim_generators) - }) - .collect() -}); - -/// Generate integral normalizer of the given point group up to its centralizer. +/// Generate integral normalizer of the given space group up to its centralizer. /// Because the factor group of the integral normalizer by the centralizer is isomorphic to a finite permutation group, the output is guaranteed to be finite. pub fn integral_normalizer( - prim_rotations: &Rotations, - prim_generators: &Rotations, -) -> Vec { + prim_operations: &Operations, + prim_generators: &Operations, + epsilon: f64, +) -> Vec { + let prim_rotations = project_rotations(prim_operations); + let prim_rotation_generators = project_rotations(prim_generators); + let rotation_types = prim_rotations .iter() .map(identify_rotation_type) .collect::>(); - let prim_generator_rotation_types = prim_generators + let prim_generator_rotation_types = prim_rotation_generators .iter() .map(identify_rotation_type) .collect::>(); @@ -67,13 +46,13 @@ pub fn integral_normalizer( .map(|v| v.iter()) .multi_cartesian_product() { - // Solve P^-1 * prim_rotations[pivot[i]] * P = prim_generators[i] (for all i) + // Solve P^-1 * prim_rotations[pivot[i]] * P = prim_rotation_generators[i] (for all i) if let Some(trans_mat_basis) = sylvester3( &pivot .iter() .map(|&i| prim_rotations[*i]) .collect::>(), - prim_generators, + &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 @@ -90,7 +69,15 @@ pub fn integral_normalizer( prim_trans_mat *= -1; } if det == 1 { - conjugators.push(prim_trans_mat); + if let Some(origin_shift) = match_origin_shift( + prim_operations, + &prim_trans_mat, + prim_generators, + epsilon, + ) { + conjugators + .push(UnimodularTransformation::new(prim_trans_mat, origin_shift)); + } break; } } @@ -98,45 +85,3 @@ pub fn integral_normalizer( } conjugators } - -#[cfg(test)] -mod tests { - use std::collections::HashSet; - - use super::*; - use crate::base::traverse; - use crate::data::PointGroupRepresentative; - - #[test] - fn test_integral_normalizer() { - for arithmetic_number in 1..=73 { - let point_group_db = - PointGroupRepresentative::from_arithmetic_crystal_class(arithmetic_number); - let prim_generators = point_group_db.primitive_generators(); - let prim_rotations = traverse(&prim_generators); - - let mut prim_rotations_set = HashSet::new(); - for rotation in prim_rotations.iter() { - prim_rotations_set.insert(rotation.clone()); - } - - let normalizer = get_point_group_normalizer(arithmetic_number).unwrap(); - assert!(normalizer.len() > 0); - - for linear in normalizer.iter() { - let linear_inv = linear - .map(|e| e as f64) - .try_inverse() - .unwrap() - .map(|e| e.round() as i32); - let prim_rotations_actual: Vec<_> = prim_rotations - .iter() - .map(|r| linear_inv * r * linear) - .collect(); - for rotation in prim_rotations_actual { - assert!(prim_rotations_set.contains(&rotation)); - } - } - } - } -} diff --git a/moyo/src/identify/space_group.rs b/moyo/src/identify/space_group.rs index 1608556..a664ef6 100644 --- a/moyo/src/identify/space_group.rs +++ b/moyo/src/identify/space_group.rs @@ -163,7 +163,7 @@ fn correction_transformation_matrices( } /// Search for origin_shift such that (trans_mat, origin_shift) transforms into -fn match_origin_shift( +pub fn match_origin_shift( prim_operations: &Operations, trans_mat: &UnimodularLinear, db_prim_generators: &Operations, diff --git a/moyo/src/lib.rs b/moyo/src/lib.rs index a85d6a9..c108fa3 100644 --- a/moyo/src/lib.rs +++ b/moyo/src/lib.rs @@ -71,13 +71,18 @@ pub mod search; // Public for benchmarking mod identify; mod symmetrize; -use base::{AngleTolerance, Cell, MoyoError, Operations, OriginShift}; +use base::{ + AngleTolerance, Cell, MagneticCell, MagneticMoment, MoyoError, Operations, OriginShift, + RotationMagneticMomentAction, +}; use data::{HallNumber, Number, Setting}; use nalgebra::Matrix3; -use crate::identify::SpaceGroup; -use crate::search::{iterative_symmetry_search, operations_in_cell}; +use crate::identify::{MagneticSpaceGroup, SpaceGroup}; +use crate::search::{ + iterative_magnetic_symmetry_search, iterative_symmetry_search, operations_in_cell, +}; use crate::symmetrize::{orbits_in_cell, StandardizedCell}; #[derive(Debug)] @@ -227,3 +232,55 @@ impl MoyoDataset { self.operations.len() } } + +#[derive(Debug)] +pub struct MoyoMagneticDataset { + // ------------------------------------------------------------------------ + // Final parameters + // ------------------------------------------------------------------------ + /// Actually used `symprec` in iterative symmetry search. + pub symprec: f64, + /// Actually used `angle_tolerance` in iterative symmetry search. + pub angle_tolerance: AngleTolerance, + /// Actually used `mag_symprec` in iterative symmetry search. + pub mag_symprec: f64, +} + +impl MoyoMagneticDataset { + pub fn new( + magnetic_cell: &MagneticCell, + symprec: f64, + angle_tolerance: AngleTolerance, + mag_symprec: Option, + action: RotationMagneticMomentAction, + ) -> Result + where + M: MagneticMoment + Clone, + { + let (prim_mag_cell, magnetic_symmetry_search, symprec, angle_tolerance, mag_symprec) = + iterative_magnetic_symmetry_search( + magnetic_cell, + symprec, + angle_tolerance, + mag_symprec, + action, + )?; + + // Magnetic space-group type identification + let epsilon = symprec + / prim_mag_cell + .magnetic_cell + .cell + .lattice + .volume() + .powf(1.0 / 3.0); + let space_group = + MagneticSpaceGroup::new(&magnetic_symmetry_search.magnetic_operations, epsilon)?; + + Ok(Self { + symprec, + angle_tolerance, + mag_symprec, + }) + } +} diff --git a/moyo/src/search.rs b/moyo/src/search.rs index 2629e26..ac6061c 100644 --- a/moyo/src/search.rs +++ b/moyo/src/search.rs @@ -9,4 +9,4 @@ pub use solve::{ pub(super) use primitive_cell::PrimitiveCell; pub(super) use primitive_symmetry_search::{operations_in_cell, PrimitiveSymmetrySearch}; -pub(super) use symmetry_search::iterative_symmetry_search; +pub(super) use symmetry_search::{iterative_magnetic_symmetry_search, iterative_symmetry_search};