From 7f63e6a2d84405c1923f1f2c22d671674f262a36 Mon Sep 17 00:00:00 2001 From: lan496 Date: Tue, 31 Dec 2024 19:18:14 +0900 Subject: [PATCH 01/11] move sylvester3 --- moyo/src/identify/point_group.rs | 175 ++++++++++++------------------- moyo/src/math.rs | 2 +- moyo/src/math/integer_system.rs | 40 ++++++- 3 files changed, 109 insertions(+), 108 deletions(-) diff --git a/moyo/src/identify/point_group.rs b/moyo/src/identify/point_group.rs index d8ec6ea..e3b34c0 100644 --- a/moyo/src/identify/point_group.rs +++ b/moyo/src/identify/point_group.rs @@ -2,14 +2,13 @@ use std::cmp::Ordering; use itertools::Itertools; use log::debug; -use nalgebra::{Dyn, Matrix3, OMatrix, OVector, U9}; use crate::base::{MoyoError, Rotation, Rotations, UnimodularLinear}; use crate::data::{ iter_arithmetic_crystal_entry, ArithmeticNumber, Centering, CrystalSystem, GeometricCrystalClass, PointGroupRepresentative, }; -use crate::math::IntegerLinearSystem; +use crate::math::sylvester3; /// Crystallographic point group with group-type information #[derive(Debug)] @@ -65,6 +64,72 @@ enum RotationType { RotoInversion6, // -6 = S3^-1 } +#[allow(dead_code)] +/// Generate integral normalizer of the given point 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. +fn integral_normalizer( + prim_rotations: &Rotations, + prim_generators: &Rotations, +) -> Vec { + let rotation_types = prim_rotations + .iter() + .map(identify_rotation_type) + .collect::>(); + let prim_generator_rotation_types = prim_generators + .iter() + .map(identify_rotation_type) + .collect::>(); + + // Try to map generators + let order = prim_rotations.len(); + let candidates: Vec> = prim_generator_rotation_types + .iter() + .map(|&rotation_type| { + (0..order) + .filter(|&i| rotation_types[i] == rotation_type) + .collect() + }) + .collect(); + + // TODO: unify with match_with_point_group + let mut conjugators = vec![]; + for pivot in candidates + .iter() + .map(|v| v.iter()) + .multi_cartesian_product() + { + // Solve P^-1 * prim_rotations[pivot[i]] * P = prim_generators[i] (for all i) + if let Some(trans_mat_basis) = sylvester3( + &pivot + .iter() + .map(|&i| prim_rotations[*i]) + .collect::>(), + prim_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 { + conjugators.push(prim_trans_mat); + break; + } + } + } + } + conjugators +} + /// Faster matching algorithm for cubic point groups fn match_with_cubic_point_group( prim_rotations: &Rotations, @@ -109,7 +174,7 @@ fn match_with_cubic_point_group( .multi_cartesian_product() { // Solve P^-1 * self.rotations[self.rotations[pivot[i]]] * P = other.generators[i] (for all i) - if let Some(trans_mat_basis) = sylvester( + if let Some(trans_mat_basis) = sylvester3( &pivot .iter() .map(|&i| prim_rotations[*i]) @@ -193,7 +258,7 @@ fn match_with_point_group( .multi_cartesian_product() { // Solve P^-1 * self.rotations[self.rotations[pivot[i]]] * P = other.generators[i] (for all i) - if let Some(trans_mat_basis) = sylvester( + if let Some(trans_mat_basis) = sylvester3( &pivot .iter() .map(|&i| prim_rotations[*i]) @@ -226,108 +291,6 @@ fn match_with_point_group( Err(MoyoError::ArithmeticCrystalClassIdentificationError) } -#[allow(dead_code)] -pub fn integral_normalizer( - prim_rotations: &Rotations, - prim_generators: &Rotations, -) -> Vec { - let rotation_types = prim_rotations - .iter() - .map(identify_rotation_type) - .collect::>(); - let prim_generator_rotation_types = prim_generators - .iter() - .map(identify_rotation_type) - .collect::>(); - - // Try to map generators - let order = prim_rotations.len(); - let candidates: Vec> = prim_generator_rotation_types - .iter() - .map(|&rotation_type| { - (0..order) - .filter(|&i| rotation_types[i] == rotation_type) - .collect() - }) - .collect(); - - // TODO: unify with match_with_point_group - let mut conjugators = vec![]; - for pivot in candidates - .iter() - .map(|v| v.iter()) - .multi_cartesian_product() - { - // Solve P^-1 * prim_rotations[prim_rotations[pivot[i]]] * P = prim_generators[i] (for all i) - if let Some(trans_mat_basis) = sylvester( - &pivot - .iter() - .map(|&i| prim_rotations[*i]) - .collect::>(), - prim_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 { - conjugators.push(prim_trans_mat); - break; - } - } - } - } - conjugators -} - -/// Solve P^-1 * A[i] * P = B[i] (for all i) -/// vec(A * P - P * B) = (I_3 \otimes A - B^T \otimes I_3) * vec(P) -fn sylvester(a: &[Matrix3], b: &[Matrix3]) -> Option>> { - let size = a.len(); - assert_eq!(size, b.len()); - - let mut coeffs = OMatrix::::zeros(9 * size); - let identity = Rotation::identity(); - for k in 0..size { - let adj = identity.kronecker(&a[k]) - b[k].transpose().kronecker(&identity); - for i in 0..9 { - for j in 0..9 { - coeffs[(9 * k + i, j)] = adj[(i, j)]; - } - } - } - let solution = IntegerLinearSystem::new(&coeffs, &OVector::::zeros(coeffs.nrows())); - - if let Some(solution) = solution { - let basis: Vec<_> = solution - .nullspace - .row_iter() - .map(|e| { - // Vectorization operator is column-major - UnimodularLinear::new( - e[0], e[1], e[2], // - e[3], e[4], e[5], // - e[6], e[7], e[8], // - ) - .transpose() - }) - .collect(); - Some(basis) - } else { - None - } -} - fn identify_rotation_type(rotation: &Rotation) -> RotationType { let tr = rotation.trace(); let det = rotation.map(|e| e as f64).determinant().round() as i32; diff --git a/moyo/src/math.rs b/moyo/src/math.rs index 76252ed..44f028e 100644 --- a/moyo/src/math.rs +++ b/moyo/src/math.rs @@ -11,6 +11,6 @@ pub use hnf::HNF; pub use snf::SNF; pub(super) use delaunay::delaunay_reduce; -pub(super) use integer_system::IntegerLinearSystem; +pub(super) use integer_system::sylvester3; pub(super) use minkowski::{is_minkowski_reduced, minkowski_reduce}; pub(super) use niggli::{is_niggli_reduced, niggli_reduce}; diff --git a/moyo/src/math/integer_system.rs b/moyo/src/math/integer_system.rs index 5aff764..f041836 100644 --- a/moyo/src/math/integer_system.rs +++ b/moyo/src/math/integer_system.rs @@ -1,5 +1,5 @@ use nalgebra::base::allocator::Allocator; -use nalgebra::{DefaultAllocator, Dim, DimMin, Dyn, OMatrix, OVector, U1}; +use nalgebra::{DefaultAllocator, Dim, DimMin, Dyn, Matrix3, OMatrix, OVector, U1, U9}; use super::snf::SNF; @@ -53,6 +53,44 @@ where } } +/// Solve P^-1 * A[i] * P = B[i] (for all i) +/// vec(A * P - P * B) = (I_3 \otimes A - B^T \otimes I_3) * vec(P) +pub fn sylvester3(a: &[Matrix3], b: &[Matrix3]) -> Option>> { + let size = a.len(); + assert_eq!(size, b.len()); + + let mut coeffs = OMatrix::::zeros(9 * size); + let identity = Matrix3::::identity(); + for k in 0..size { + let adj = identity.kronecker(&a[k]) - b[k].transpose().kronecker(&identity); + for i in 0..9 { + for j in 0..9 { + coeffs[(9 * k + i, j)] = adj[(i, j)]; + } + } + } + let solution = IntegerLinearSystem::new(&coeffs, &OVector::::zeros(coeffs.nrows())); + + if let Some(solution) = solution { + let basis: Vec<_> = solution + .nullspace + .row_iter() + .map(|e| { + // Vectorization operator is column-major + Matrix3::::new( + e[0], e[1], e[2], // + e[3], e[4], e[5], // + e[6], e[7], e[8], // + ) + .transpose() + }) + .collect(); + Some(basis) + } else { + None + } +} + #[cfg(test)] mod tests { use nalgebra::{matrix, vector}; From 19dbef41d0cd3cabb7fdd16d69dc1564d6e3a6d1 Mon Sep 17 00:00:00 2001 From: lan496 Date: Tue, 31 Dec 2024 20:01:59 +0900 Subject: [PATCH 02/11] Lazy-init point group normalizer --- Cargo.lock | 1 + moyo/Cargo.toml | 1 + moyo/src/identify.rs | 2 + moyo/src/identify/normalizer.rs | 143 +++++++++++++++++++++++++++++ moyo/src/identify/point_group.rs | 137 +-------------------------- moyo/src/identify/rotation_type.rs | 34 +++++++ 6 files changed, 184 insertions(+), 134 deletions(-) create mode 100644 moyo/src/identify/normalizer.rs create mode 100644 moyo/src/identify/rotation_type.rs diff --git a/Cargo.lock b/Cargo.lock index ebf09d6..7d6ad3b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -626,6 +626,7 @@ dependencies = [ "kiddo", "log", "nalgebra", + "once_cell", "rand", "rstest", "serde", diff --git a/moyo/Cargo.toml b/moyo/Cargo.toml index b42233c..de9a80b 100644 --- a/moyo/Cargo.toml +++ b/moyo/Cargo.toml @@ -25,6 +25,7 @@ union-find = "0.4" strum = "0.25" strum_macros = "0.25" kiddo = "4.2" +once_cell = "1.20.2" [dev-dependencies] rand = "0.8" diff --git a/moyo/src/identify.rs b/moyo/src/identify.rs index 1e5ca66..e9b1e89 100644 --- a/moyo/src/identify.rs +++ b/moyo/src/identify.rs @@ -1,4 +1,6 @@ +mod normalizer; mod point_group; +mod rotation_type; mod space_group; pub(super) use space_group::SpaceGroup; diff --git a/moyo/src/identify/normalizer.rs b/moyo/src/identify/normalizer.rs new file mode 100644 index 0000000..5303447 --- /dev/null +++ b/moyo/src/identify/normalizer.rs @@ -0,0 +1,143 @@ +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 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() +}); + +#[allow(dead_code)] +/// Generate integral normalizer of the given point 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. +fn integral_normalizer( + prim_rotations: &Rotations, + prim_generators: &Rotations, +) -> Vec { + let rotation_types = prim_rotations + .iter() + .map(identify_rotation_type) + .collect::>(); + let prim_generator_rotation_types = prim_generators + .iter() + .map(identify_rotation_type) + .collect::>(); + + // Try to map generators + let order = prim_rotations.len(); + let candidates: Vec> = prim_generator_rotation_types + .iter() + .map(|&rotation_type| { + (0..order) + .filter(|&i| rotation_types[i] == rotation_type) + .collect() + }) + .collect(); + + // TODO: unify with match_with_point_group + let mut conjugators = vec![]; + for pivot in candidates + .iter() + .map(|v| v.iter()) + .multi_cartesian_product() + { + // Solve P^-1 * prim_rotations[pivot[i]] * P = prim_generators[i] (for all i) + if let Some(trans_mat_basis) = sylvester3( + &pivot + .iter() + .map(|&i| prim_rotations[*i]) + .collect::>(), + prim_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 { + conjugators.push(prim_trans_mat); + break; + } + } + } + } + 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/point_group.rs b/moyo/src/identify/point_group.rs index e3b34c0..c148fe6 100644 --- a/moyo/src/identify/point_group.rs +++ b/moyo/src/identify/point_group.rs @@ -3,7 +3,8 @@ use std::cmp::Ordering; use itertools::Itertools; use log::debug; -use crate::base::{MoyoError, Rotation, Rotations, UnimodularLinear}; +use super::rotation_type::{identify_rotation_type, RotationType}; +use crate::base::{MoyoError, Rotations, UnimodularLinear}; use crate::data::{ iter_arithmetic_crystal_entry, ArithmeticNumber, Centering, CrystalSystem, GeometricCrystalClass, PointGroupRepresentative, @@ -50,86 +51,6 @@ impl PointGroup { } } -#[derive(Debug, Copy, Clone, PartialEq)] -enum RotationType { - Rotation1, // 1 - Rotation2, // 2 - Rotation3, // 3 - Rotation4, // 4 - Rotation6, // 6 - RotoInversion1, // -1 = S2 - RotoInversion2, // -2 = m = S1 - RotoInversion3, // -3 = S6^-1 - RotoInversion4, // -4 = S4^-1 - RotoInversion6, // -6 = S3^-1 -} - -#[allow(dead_code)] -/// Generate integral normalizer of the given point 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. -fn integral_normalizer( - prim_rotations: &Rotations, - prim_generators: &Rotations, -) -> Vec { - let rotation_types = prim_rotations - .iter() - .map(identify_rotation_type) - .collect::>(); - let prim_generator_rotation_types = prim_generators - .iter() - .map(identify_rotation_type) - .collect::>(); - - // Try to map generators - let order = prim_rotations.len(); - let candidates: Vec> = prim_generator_rotation_types - .iter() - .map(|&rotation_type| { - (0..order) - .filter(|&i| rotation_types[i] == rotation_type) - .collect() - }) - .collect(); - - // TODO: unify with match_with_point_group - let mut conjugators = vec![]; - for pivot in candidates - .iter() - .map(|v| v.iter()) - .multi_cartesian_product() - { - // Solve P^-1 * prim_rotations[pivot[i]] * P = prim_generators[i] (for all i) - if let Some(trans_mat_basis) = sylvester3( - &pivot - .iter() - .map(|&i| prim_rotations[*i]) - .collect::>(), - prim_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 { - conjugators.push(prim_trans_mat); - break; - } - } - } - } - conjugators -} - /// Faster matching algorithm for cubic point groups fn match_with_cubic_point_group( prim_rotations: &Rotations, @@ -291,25 +212,6 @@ fn match_with_point_group( Err(MoyoError::ArithmeticCrystalClassIdentificationError) } -fn identify_rotation_type(rotation: &Rotation) -> RotationType { - let tr = rotation.trace(); - let det = rotation.map(|e| e as f64).determinant().round() as i32; - - match (tr, det) { - (3, 1) => RotationType::Rotation1, - (-1, 1) => RotationType::Rotation2, - (0, 1) => RotationType::Rotation3, - (1, 1) => RotationType::Rotation4, - (2, 1) => RotationType::Rotation6, - (-3, -1) => RotationType::RotoInversion1, - (1, -1) => RotationType::RotoInversion2, - (0, -1) => RotationType::RotoInversion3, - (-1, -1) => RotationType::RotoInversion4, - (-2, -1) => RotationType::RotoInversion6, - _ => unreachable!("Unknown rotation type"), - } -} - /// Use look up table in Table 6 of https://arxiv.org/pdf/1808.01590.pdf fn identify_geometric_crystal_class( rotation_types: &Vec, @@ -385,7 +287,7 @@ fn identify_geometric_crystal_class( mod tests { use std::collections::HashSet; - use super::{integral_normalizer, PointGroup, PointGroupRepresentative}; + use super::*; use crate::base::traverse; #[test] @@ -427,37 +329,4 @@ mod tests { } } } - - #[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 conjugators = integral_normalizer(&prim_rotations, &prim_generators); - assert!(conjugators.len() > 0); - - for linear in conjugators.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/rotation_type.rs b/moyo/src/identify/rotation_type.rs new file mode 100644 index 0000000..fdc0eac --- /dev/null +++ b/moyo/src/identify/rotation_type.rs @@ -0,0 +1,34 @@ +use crate::base::Rotation; + +#[derive(Debug, Copy, Clone, PartialEq)] +pub enum RotationType { + Rotation1, // 1 + Rotation2, // 2 + Rotation3, // 3 + Rotation4, // 4 + Rotation6, // 6 + RotoInversion1, // -1 = S2 + RotoInversion2, // -2 = m = S1 + RotoInversion3, // -3 = S6^-1 + RotoInversion4, // -4 = S4^-1 + RotoInversion6, // -6 = S3^-1 +} + +pub fn identify_rotation_type(rotation: &Rotation) -> RotationType { + let tr = rotation.trace(); + let det = rotation.map(|e| e as f64).determinant().round() as i32; + + match (tr, det) { + (3, 1) => RotationType::Rotation1, + (-1, 1) => RotationType::Rotation2, + (0, 1) => RotationType::Rotation3, + (1, 1) => RotationType::Rotation4, + (2, 1) => RotationType::Rotation6, + (-3, -1) => RotationType::RotoInversion1, + (1, -1) => RotationType::RotoInversion2, + (0, -1) => RotationType::RotoInversion3, + (-1, -1) => RotationType::RotoInversion4, + (-2, -1) => RotationType::RotoInversion6, + _ => unreachable!("Unknown rotation type"), + } +} From e53d26b13c1c6acbfd573e1cac3d09a38c56af54 Mon Sep 17 00:00:00 2001 From: lan496 Date: Tue, 31 Dec 2024 22:26:30 +0900 Subject: [PATCH 03/11] Bump up kiddo --- Cargo.lock | 48 +++++++++++++++++++- moyo/Cargo.toml | 2 +- moyo/src/search/primitive_symmetry_search.rs | 2 +- moyo/src/search/solve.rs | 37 ++++++++------- 4 files changed, 66 insertions(+), 23 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 7d6ad3b..ba206af 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -11,6 +11,15 @@ dependencies = [ "memchr", ] +[[package]] +name = "aligned-vec" +version = "0.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7e0966165eaf052580bd70eb1b32cb3d6245774c0104d1b2793e9650bf83b52a" +dependencies = [ + "equator", +] + [[package]] name = "anes" version = "0.1.6" @@ -81,6 +90,12 @@ version = "1.7.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "69f7f8c3906b62b754cd5326047894316021dcfe5a194c8ea52bdd94934a3457" +[[package]] +name = "array-init" +version = "2.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3d62b7694a562cdf5a74227903507c56ab2cc8bdd1f781ed5cb4cf9c9f810bfc" + [[package]] name = "autocfg" version = "1.4.0" @@ -181,6 +196,12 @@ version = "0.7.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f46ad14479a25103f283c0f10005961cf086d8dc42205bb44c46ac563475dca6" +[[package]] +name = "cmov" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1b1dc960ba75d543267db9254da8ec1cb318a037beb3f8d2497520e410096fab" + [[package]] name = "colorchoice" version = "1.0.3" @@ -301,6 +322,26 @@ dependencies = [ "log", ] +[[package]] +name = "equator" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c35da53b5a021d2484a7cc49b2ac7f2d840f8236a286f84202369bd338d761ea" +dependencies = [ + "equator-macro", +] + +[[package]] +name = "equator-macro" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3bf679796c0322556351f287a51b49e48f7c4986e727b5dd78c972d30e2e16cc" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + [[package]] name = "fixed" version = "1.28.0" @@ -532,11 +573,14 @@ dependencies = [ [[package]] name = "kiddo" -version = "4.2.1" +version = "5.0.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "60c5fcd3044b774e2c80a502b2387b75d1baa95e99b2bceeb5db00f2e2d27fe9" +checksum = "c6e580022eee0bcd05b8c1c43d3512e79bfa001ccf88ef3dcfd1a3ac74b4e390" dependencies = [ + "aligned-vec", + "array-init", "az", + "cmov", "divrem", "doc-comment", "elapsed", diff --git a/moyo/Cargo.toml b/moyo/Cargo.toml index de9a80b..6b385e5 100644 --- a/moyo/Cargo.toml +++ b/moyo/Cargo.toml @@ -24,7 +24,7 @@ thiserror = "1.0" union-find = "0.4" strum = "0.25" strum_macros = "0.25" -kiddo = "4.2" +kiddo = "5.0.3" once_cell = "1.20.2" [dev-dependencies] diff --git a/moyo/src/search/primitive_symmetry_search.rs b/moyo/src/search/primitive_symmetry_search.rs index 8ca2676..14f65cf 100644 --- a/moyo/src/search/primitive_symmetry_search.rs +++ b/moyo/src/search/primitive_symmetry_search.rs @@ -13,7 +13,7 @@ use super::{ }; use crate::base::{ traverse, AngleTolerance, Cell, Lattice, MagneticCell, MagneticMoment, MagneticOperation, - MagneticOperations, MoyoError, Operation, Operations, Permutation, Position, Rotation, + MagneticOperations, MoyoError, Operation, Operations, Permutation, Rotation, RotationMagneticMomentAction, Rotations, Transformation, EPS, }; diff --git a/moyo/src/search/solve.rs b/moyo/src/search/solve.rs index 8809087..6a7b110 100644 --- a/moyo/src/search/solve.rs +++ b/moyo/src/search/solve.rs @@ -1,7 +1,7 @@ -use std::collections::BTreeMap; +use std::{collections::BTreeMap, num::NonZero}; use itertools::iproduct; -use kiddo::{KdTree, SquaredEuclidean}; +use kiddo::{ImmutableKdTree, SquaredEuclidean}; use nalgebra::{Rotation3, Vector3}; use crate::base::{AtomicSpecie, Cell, Lattice, Permutation, Position, Rotation, Translation}; @@ -10,7 +10,7 @@ use crate::base::{AtomicSpecie, Cell, Lattice, Permutation, Position, Rotation, pub struct PeriodicKdTree { num_sites: usize, lattice: Lattice, - kdtree: KdTree, + kdtree: ImmutableKdTree, indices: Vec, symprec: f64, } @@ -59,7 +59,7 @@ impl PeriodicKdTree { Self { num_sites: reduced_cell.num_atoms(), lattice: new_lattice, - kdtree: (&entries).into(), + kdtree: ImmutableKdTree::new_from_slice(&entries), indices, symprec, } @@ -70,26 +70,25 @@ impl PeriodicKdTree { let mut wrapped_position = *position; wrapped_position -= wrapped_position.map(|e| e.floor()); // [0, 1) let cart_coords = self.lattice.cartesian_coords(&wrapped_position); - let within = self.kdtree.nearest_n_within::( + let mut within = self.kdtree.best_n_within::( &[cart_coords.x, cart_coords.y, cart_coords.z], self.symprec.powi(2), // squared distance for KdTree - 1, - false, + NonZero::new(1).unwrap(), ); - if within.is_empty() { - return None; - } + if let Some(entry) = within.next() { + let item = entry.item as usize; + let distance = entry.distance.sqrt(); + if distance > self.symprec { + return None; + } - let item = within[0].item as usize; - let distance = within[0].distance.sqrt(); - if distance > self.symprec { - return None; + Some(PeriodicNeighbor { + index: self.indices[item], + distance, + }) + } else { + None } - - Some(PeriodicNeighbor { - index: self.indices[item], - distance, - }) } } From 6dc2139483179934527b73af8f0b4c017a4817e8 Mon Sep 17 00:00:00 2001 From: lan496 Date: Fri, 3 Jan 2025 10:16:44 +0900 Subject: [PATCH 04/11] WIP: magnetic space group type identification --- moyo/src/base/error.rs | 2 + moyo/src/base/transformation.rs | 62 +++++++- moyo/src/data.rs | 4 +- moyo/src/data/hall_symbol.rs | 15 ++ .../src/data/magnetic_hall_symbol_database.rs | 12 +- moyo/src/data/magnetic_space_group.rs | 43 +++++- moyo/src/identify.rs | 1 + moyo/src/identify/magnetic_space_group.rs | 137 ++++++++++++++++++ 8 files changed, 257 insertions(+), 19 deletions(-) create mode 100644 moyo/src/identify/magnetic_space_group.rs diff --git a/moyo/src/base/error.rs b/moyo/src/base/error.rs index 880ea11..a721f30 100644 --- a/moyo/src/base/error.rs +++ b/moyo/src/base/error.rs @@ -27,6 +27,8 @@ pub enum MoyoError { ArithmeticCrystalClassIdentificationError, #[error("Space group type identification failed")] SpaceGroupTypeIdentificationError, + #[error("Magnetic space group type identification failed")] + MagneticSpaceGroupTypeIdentificationError, #[error("Standardization failed")] StandardizationError, #[error("Wyckoff position assignment failed")] diff --git a/moyo/src/base/transformation.rs b/moyo/src/base/transformation.rs index 89c2f99..aef1e92 100644 --- a/moyo/src/base/transformation.rs +++ b/moyo/src/base/transformation.rs @@ -4,7 +4,7 @@ use nalgebra::base::{Matrix3, Vector3}; use super::cell::Cell; use super::lattice::Lattice; use super::magnetic_cell::{MagneticCell, MagneticMoment}; -use super::operation::Operation; +use super::operation::{MagneticOperation, Operation}; use crate::math::SNF; pub type UnimodularLinear = Matrix3; @@ -73,6 +73,24 @@ impl UnimodularTransformation { .collect() } + pub fn transform_magnetic_operation( + &self, + magnetic_operation: &MagneticOperation, + ) -> MagneticOperation { + let new_operation = self.transform_operation(&magnetic_operation.operation); + MagneticOperation::from_operation(new_operation, magnetic_operation.time_reversal) + } + + pub fn transform_magnetic_operations( + &self, + magnetic_operations: &[MagneticOperation], + ) -> Vec { + magnetic_operations + .iter() + .map(|mops| self.transform_magnetic_operation(mops)) + .collect() + } + pub fn transform_cell(&self, cell: &Cell) -> Cell { let new_lattice = self.transform_lattice(&cell.lattice); // (P, p)^-1 x = P^-1 (x - p) @@ -188,6 +206,48 @@ impl Transformation { .collect() } + pub fn transform_magnetic_operation( + &self, + magnetic_operation: &MagneticOperation, + ) -> Option { + let new_operation = self.transform_operation(&magnetic_operation.operation)?; + Some(MagneticOperation::from_operation( + new_operation, + magnetic_operation.time_reversal, + )) + } + + pub fn transform_magnetic_operations( + &self, + magnetic_operations: &[MagneticOperation], + ) -> Vec { + magnetic_operations + .iter() + .filter_map(|mops| self.transform_magnetic_operation(mops)) + .collect() + } + + pub fn inverse_transform_magnetic_operation( + &self, + magnetic_operation: &MagneticOperation, + ) -> Option { + let new_operation = self.inverse_transform_operation(&magnetic_operation.operation)?; + Some(MagneticOperation::from_operation( + new_operation, + magnetic_operation.time_reversal, + )) + } + + pub fn inverse_transform_magnetic_operations( + &self, + magnetic_operations: &[MagneticOperation], + ) -> Vec { + magnetic_operations + .iter() + .filter_map(|mops| self.inverse_transform_magnetic_operation(mops)) + .collect() + } + // The transformation may increase the number of atoms in the cell. // Return the transformed cell and mapping from sites in the transformed cell to sites in the original cell. pub fn transform_cell(&self, cell: &Cell) -> (Cell, Vec) { diff --git a/moyo/src/data.rs b/moyo/src/data.rs index 6985a92..85d86ad 100644 --- a/moyo/src/data.rs +++ b/moyo/src/data.rs @@ -11,13 +11,15 @@ mod wyckoff; pub use arithmetic_crystal_class::ArithmeticNumber; pub use centering::Centering; -pub use hall_symbol::HallSymbol; +pub use hall_symbol::{HallSymbol, MagneticHallSymbol}; pub use hall_symbol_database::{hall_symbol_entry, HallNumber, HallSymbolEntry, Number}; +pub use magnetic_space_group::{ConstructType, UNINumber, NUM_MAGNETIC_SPACE_GROUP_TYPES}; pub use setting::Setting; pub(super) use arithmetic_crystal_class::{ arithmetic_crystal_class_entry, iter_arithmetic_crystal_entry, }; pub(super) use classification::{CrystalSystem, GeometricCrystalClass, LatticeSystem}; +pub(super) use magnetic_space_group::uni_number_range; pub(super) use point_group::PointGroupRepresentative; pub(super) use wyckoff::{iter_wyckoff_positions, WyckoffPosition, WyckoffPositionSpace}; diff --git a/moyo/src/data/hall_symbol.rs b/moyo/src/data/hall_symbol.rs index c610fae..708ed7b 100644 --- a/moyo/src/data/hall_symbol.rs +++ b/moyo/src/data/hall_symbol.rs @@ -6,6 +6,8 @@ use nalgebra::{matrix, Vector3}; use super::centering::Centering; use super::hall_symbol_database::{hall_symbol_entry, HallNumber}; +use super::magnetic_hall_symbol_database::magnetic_hall_symbol_entry; +use super::magnetic_space_group::UNINumber; use crate::base::{ MagneticOperation, MagneticOperations, Operation, Operations, OriginShift, Rotation, TimeReversal, Transformation, Translation, EPS, @@ -206,6 +208,19 @@ impl MagneticHallSymbol { operations } + + pub fn from_uni_number(uni_number: UNINumber) -> Option { + if let Some(entry) = magnetic_hall_symbol_entry(uni_number) { + Self::new(entry.magnetic_hall_symbol) + } else { + None + } + } + + pub fn primitive_generators(&self) -> MagneticOperations { + Transformation::from_linear(self.centering.linear()) + .inverse_transform_magnetic_operations(&self.generators) + } } /// Tokenize string by whitespaces diff --git a/moyo/src/data/magnetic_hall_symbol_database.rs b/moyo/src/data/magnetic_hall_symbol_database.rs index 32b8dd3..e790ef1 100644 --- a/moyo/src/data/magnetic_hall_symbol_database.rs +++ b/moyo/src/data/magnetic_hall_symbol_database.rs @@ -1,4 +1,4 @@ -use super::magnetic_space_group::UNINumber; +use super::magnetic_space_group::{UNINumber, NUM_MAGNETIC_SPACE_GROUP_TYPES}; #[derive(Debug, Clone)] pub struct MagneticHallSymbolEntry { @@ -15,21 +15,13 @@ impl MagneticHallSymbolEntry { } } -// TODO: -// smallest and largest UNI numbers for each Hall number -// let hall_number_to_uni_numbers: HashMap; - -// TODO: -// smallest and largest Hall numbers for XSG of each UNI number -// let uni_number_to_xsg_hall_numbers: HashMap; - pub fn magnetic_hall_symbol_entry(uni_number: UNINumber) -> Option { MAGNETIC_HALL_SYMBOL_DATABASE .get((uni_number - 1) as usize) .cloned() } -const MAGNETIC_HALL_SYMBOL_DATABASE: [MagneticHallSymbolEntry; 1651] = [ +const MAGNETIC_HALL_SYMBOL_DATABASE: [MagneticHallSymbolEntry; NUM_MAGNETIC_SPACE_GROUP_TYPES] = [ MagneticHallSymbolEntry::new("P 1", 1), MagneticHallSymbolEntry::new("P 1 1'", 2), MagneticHallSymbolEntry::new("P 1 1c'", 3), diff --git a/moyo/src/data/magnetic_space_group.rs b/moyo/src/data/magnetic_space_group.rs index 1ecb0d9..6eecaca 100644 --- a/moyo/src/data/magnetic_space_group.rs +++ b/moyo/src/data/magnetic_space_group.rs @@ -1,9 +1,15 @@ -use super::hall_symbol_database::Number; +use std::ops::RangeInclusive; + +use once_cell::sync::Lazy; + +use super::{hall_symbol_database::Number, HallNumber}; + +pub const NUM_MAGNETIC_SPACE_GROUP_TYPES: usize = 1651; /// UNI Number for magnetic space groups (1 - 1651) pub type UNINumber = i32; -#[derive(Debug, Clone)] +#[derive(Debug, Clone, PartialEq, Eq)] pub enum ConstructType { Type1, Type2, @@ -17,8 +23,8 @@ pub struct MagneticSpaceGroupType { pub litvin_number: i32, pub bns_number: &'static str, pub og_number: &'static str, - /// ITA number for maximal space subgroup (XSG) - pub xsg_number: Number, + /// ITA number for reference space group in BNS setting + pub number: Number, pub construct_type: ConstructType, } @@ -28,7 +34,7 @@ impl MagneticSpaceGroupType { litvin_number: i32, bns_number: &'static str, og_number: &'static str, - xsg_number: Number, + number: Number, construct_type: ConstructType, ) -> Self { Self { @@ -36,19 +42,42 @@ impl MagneticSpaceGroupType { litvin_number, bns_number, og_number, - xsg_number, + number, construct_type, } } } +/// Return inclusive range for UNI numbers associated with the given ITA number. +pub fn uni_number_range(number: Number) -> Option> { + ITA_NUMBER_TO_UNI_NUMBERS + .get((number - 1) as usize) + .cloned() +} + +static ITA_NUMBER_TO_UNI_NUMBERS: Lazy>> = Lazy::new(|| { + let mut ret = vec![]; + let mut start = 1; + for uni_number in 1..=NUM_MAGNETIC_SPACE_GROUP_TYPES { + if (uni_number == NUM_MAGNETIC_SPACE_GROUP_TYPES) + || (MAGNETIC_SPACE_GROUP_TYPES[uni_number - 1].number + != MAGNETIC_SPACE_GROUP_TYPES[uni_number].number) + { + ret.push(start..=uni_number as UNINumber); + start = uni_number as UNINumber + 1; + } + } + assert_eq!(ret.len(), 230); + ret +}); + pub fn get_magnetic_space_group_type(uni_number: UNINumber) -> Option { MAGNETIC_SPACE_GROUP_TYPES .get((uni_number - 1) as usize) .cloned() } -const MAGNETIC_SPACE_GROUP_TYPES: [MagneticSpaceGroupType; 1651] = [ +const MAGNETIC_SPACE_GROUP_TYPES: [MagneticSpaceGroupType; NUM_MAGNETIC_SPACE_GROUP_TYPES] = [ MagneticSpaceGroupType::new(1, 1, "1.1", "1.1.1", 1, ConstructType::Type1), MagneticSpaceGroupType::new(2, 2, "1.2", "1.2.2", 1, ConstructType::Type2), MagneticSpaceGroupType::new(3, 3, "1.3", "1.3.3", 1, ConstructType::Type4), diff --git a/moyo/src/identify.rs b/moyo/src/identify.rs index e9b1e89..2e92d48 100644 --- a/moyo/src/identify.rs +++ b/moyo/src/identify.rs @@ -1,3 +1,4 @@ +mod magnetic_space_group; mod normalizer; mod point_group; mod rotation_type; diff --git a/moyo/src/identify/magnetic_space_group.rs b/moyo/src/identify/magnetic_space_group.rs new file mode 100644 index 0000000..0691ac7 --- /dev/null +++ b/moyo/src/identify/magnetic_space_group.rs @@ -0,0 +1,137 @@ +use std::collections::HashSet; + +use log::debug; + +use super::space_group::SpaceGroup; +use crate::base::{MagneticOperations, MoyoError, Operations, Rotation, UnimodularTransformation}; +use crate::data::{uni_number_range, ConstructType, Setting, UNINumber}; + +#[derive(Debug)] +pub struct MagneticSpaceGroup { + pub uni_number: UNINumber, + /// Transformation to the representative for `uni_number` in primitive + pub transformation: UnimodularTransformation, +} + +impl MagneticSpaceGroup { + /// Identify the magnetic space group type from the primitive magnetic operations. + /// epsilon: tolerance for comparing translation parts + pub fn new(prim_mag_operations: &MagneticOperations, epsilon: f64) -> Result { + let prim_xsg = maximal_space_subgroup_from_magnetic_space_group(prim_mag_operations) + .ok_or(MoyoError::MagneticSpaceGroupTypeIdentificationError)?; + let prim_fsg = family_space_group_from_magnetic_space_group(prim_mag_operations) + .ok_or(MoyoError::MagneticSpaceGroupTypeIdentificationError)?; + + let construct_type = match ( + prim_mag_operations.len() % prim_xsg.len(), + prim_fsg.len() % prim_xsg.len(), + ) { + (1, 1) => ConstructType::Type1, + (2, 1) => ConstructType::Type2, + (2, 2) => { + // Find coset representatives of MSG/XSG + let identify = Rotation::identity(); + if prim_mag_operations + .iter() + .any(|mops| mops.time_reversal && mops.operation.rotation == identify) + { + // Anti-translation + ConstructType::Type4 + } else { + ConstructType::Type3 + } + } + _ => { + return Err(MoyoError::MagneticSpaceGroupTypeIdentificationError); + } + }; + + // BNS setting + // std_ref_spg.transformation: primitive input -> primitive BNS setting + let ref_spg = if construct_type == ConstructType::Type4 { + prim_xsg + } else { + prim_fsg + }; + let std_ref_spg = SpaceGroup::new(&ref_spg, Setting::Standard, epsilon)?; + let tmp_prim_mag_operations = std_ref_spg + .transformation + .transform_magnetic_operations(&prim_mag_operations); + + let uni_number_range = uni_number_range(std_ref_spg.number) + .ok_or(MoyoError::SpaceGroupTypeIdentificationError)?; + for uni_number in uni_number_range { + todo!() + } + + todo!() + } +} + +fn maximal_space_subgroup_from_magnetic_space_group( + prim_mag_operations: &MagneticOperations, +) -> Option { + let mut xsg = vec![]; + let mut visited = HashSet::new(); + + for mops in prim_mag_operations { + if mops.time_reversal { + continue; + } + + if visited.contains(&mops.operation.rotation) { + debug!("Input magnetic operations are not in primitive."); + return None; + } + xsg.push(mops.operation.clone()); + visited.insert(mops.operation.rotation.clone()); + } + + if (xsg.len() != prim_mag_operations.len()) && (xsg.len() * 2 != prim_mag_operations.len()) { + debug!("Input magnetic operations are incomplete."); + return None; + } + Some(xsg) +} + +fn family_space_group_from_magnetic_space_group( + prim_mag_operations: &MagneticOperations, +) -> Option { + let mut fsg = vec![]; + let mut visited = HashSet::new(); + + for mops in prim_mag_operations { + if visited.contains(&mops.operation.rotation) { + continue; + } + fsg.push(mops.operation.clone()); + visited.insert(mops.operation.rotation.clone()); + } + + if (fsg.len() != prim_mag_operations.len()) && (fsg.len() * 2 != prim_mag_operations.len()) { + debug!("Input magnetic operations are incomplete."); + return None; + } + Some(fsg) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::base::Transformation; + use crate::data::MagneticHallSymbol; + + #[test] + fn test_type4() { + // -P 4 2 3 1abc' (UNI No. 1599) + // P_I m-3m (BNS No. 221.97) + let mhs = MagneticHallSymbol::from_uni_number(1599).unwrap(); + let magnetic_operations = mhs.traverse(); + + // conventional -> primitive + let prim_mag_operations = Transformation::from_linear(mhs.centering.linear()) + .inverse_transform_magnetic_operations(&magnetic_operations); + + let magnetic_space_group = MagneticSpaceGroup::new(&prim_mag_operations, 1e-8).unwrap(); + } +} From 91e46aca596d5697b221f917850c1af12f0b6af2 Mon Sep 17 00:00:00 2001 From: lan496 Date: Sat, 4 Jan 2025 17:02:44 +0900 Subject: [PATCH 05/11] Fix construct type identification --- moyo/src/base/error.rs | 2 + moyo/src/data.rs | 4 +- .../src/data/magnetic_hall_symbol_database.rs | 1 + moyo/src/data/magnetic_space_group.rs | 2 +- moyo/src/identify/magnetic_space_group.rs | 345 ++++++++++++++---- moyo/src/identify/normalizer.rs | 3 +- moyo/src/identify/space_group.rs | 3 +- 7 files changed, 274 insertions(+), 86 deletions(-) diff --git a/moyo/src/base/error.rs b/moyo/src/base/error.rs index a721f30..a4621b8 100644 --- a/moyo/src/base/error.rs +++ b/moyo/src/base/error.rs @@ -27,6 +27,8 @@ pub enum MoyoError { ArithmeticCrystalClassIdentificationError, #[error("Space group type identification failed")] SpaceGroupTypeIdentificationError, + #[error("Construct type identification failed")] + ConstructTypeIdentificationError, #[error("Magnetic space group type identification failed")] MagneticSpaceGroupTypeIdentificationError, #[error("Standardization failed")] diff --git a/moyo/src/data.rs b/moyo/src/data.rs index 85d86ad..8f2a796 100644 --- a/moyo/src/data.rs +++ b/moyo/src/data.rs @@ -13,7 +13,9 @@ pub use arithmetic_crystal_class::ArithmeticNumber; pub use centering::Centering; pub use hall_symbol::{HallSymbol, MagneticHallSymbol}; pub use hall_symbol_database::{hall_symbol_entry, HallNumber, HallSymbolEntry, Number}; -pub use magnetic_space_group::{ConstructType, UNINumber, NUM_MAGNETIC_SPACE_GROUP_TYPES}; +pub use magnetic_space_group::{ + get_magnetic_space_group_type, ConstructType, UNINumber, NUM_MAGNETIC_SPACE_GROUP_TYPES, +}; pub use setting::Setting; pub(super) use arithmetic_crystal_class::{ diff --git a/moyo/src/data/magnetic_hall_symbol_database.rs b/moyo/src/data/magnetic_hall_symbol_database.rs index e790ef1..5bfbb79 100644 --- a/moyo/src/data/magnetic_hall_symbol_database.rs +++ b/moyo/src/data/magnetic_hall_symbol_database.rs @@ -21,6 +21,7 @@ pub fn magnetic_hall_symbol_entry(uni_number: UNINumber) -> Option Result { - let prim_xsg = maximal_space_subgroup_from_magnetic_space_group(prim_mag_operations) - .ok_or(MoyoError::MagneticSpaceGroupTypeIdentificationError)?; - let prim_fsg = family_space_group_from_magnetic_space_group(prim_mag_operations) - .ok_or(MoyoError::MagneticSpaceGroupTypeIdentificationError)?; - - let construct_type = match ( - prim_mag_operations.len() % prim_xsg.len(), - prim_fsg.len() % prim_xsg.len(), - ) { - (1, 1) => ConstructType::Type1, - (2, 1) => ConstructType::Type2, - (2, 2) => { - // Find coset representatives of MSG/XSG - let identify = Rotation::identity(); - if prim_mag_operations - .iter() - .any(|mops| mops.time_reversal && mops.operation.rotation == identify) - { - // Anti-translation - ConstructType::Type4 - } else { - ConstructType::Type3 - } - } - _ => { - return Err(MoyoError::MagneticSpaceGroupTypeIdentificationError); - } - }; - - // BNS setting + let (ref_spg, construct_type) = + identify_reference_space_group(prim_mag_operations, epsilon) + .ok_or(MoyoError::ConstructTypeIdentificationError)?; + debug!("Construct type: {:?}", construct_type); + let setting = Setting::Standard; // std_ref_spg.transformation: primitive input -> primitive BNS setting - let ref_spg = if construct_type == ConstructType::Type4 { - prim_xsg - } else { - prim_fsg - }; - let std_ref_spg = SpaceGroup::new(&ref_spg, Setting::Standard, epsilon)?; - let tmp_prim_mag_operations = std_ref_spg - .transformation - .transform_magnetic_operations(&prim_mag_operations); + let std_ref_spg = SpaceGroup::new(&ref_spg, setting, epsilon)?; + debug!("Reference space group: {:?}", std_ref_spg.number); let uni_number_range = uni_number_range(std_ref_spg.number) .ok_or(MoyoError::SpaceGroupTypeIdentificationError)?; for uni_number in uni_number_range { - todo!() + if get_magnetic_space_group_type(uni_number) + .unwrap() + .construct_type + != construct_type + { + continue; + } + + 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` + // 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), + }); + } + } } + Err(MoyoError::MagneticSpaceGroupTypeIdentificationError) + } - todo!() + fn get_db_reference_space_group_primitive_generators( + mhs: &MagneticHallSymbol, + construct_type: &ConstructType, + ) -> Operations { + let db_prim_mag_operations = mhs.primitive_generators(); + match construct_type { + ConstructType::Type1 | ConstructType::Type2 | ConstructType::Type3 => { + // Reference space group: FSG + // -> Remove time-reversal parts + db_prim_mag_operations + .iter() + .map(|mops| mops.operation.clone()) + .collect() + } + ConstructType::Type4 => { + // Reference space group: XSG + // -> Remove anti-translation operation + let identity = Rotation::identity(); + db_prim_mag_operations + .iter() + .filter_map(|mops| { + // Here we assume the magnetic Hall symbol composes of XSG generators and one anti-translation operation + if mops.time_reversal { + assert_eq!(mops.operation.rotation, identity); + None + } else { + Some(mops.operation.clone()) + } + }) + .collect() + } + } } -} -fn maximal_space_subgroup_from_magnetic_space_group( - prim_mag_operations: &MagneticOperations, -) -> Option { - let mut xsg = vec![]; - let mut visited = HashSet::new(); + 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, + )) + } + } + } - for mops in prim_mag_operations { - if mops.time_reversal { - continue; + /// 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( + prim_mag_operations: &MagneticOperations, + trans_mat: &UnimodularLinear, + db_prim_mag_generators: &MagneticOperations, + epsilon: f64, + ) -> Option { + let new_prim_mag_operations = UnimodularTransformation::from_linear(*trans_mat) + .transform_magnetic_operations(prim_mag_operations); + let mut hm_translations = HashMap::new(); + for mops in new_prim_mag_operations.iter() { + hm_translations.insert( + (mops.operation.rotation, mops.time_reversal), + mops.operation.translation, + ); } - if visited.contains(&mops.operation.rotation) { - debug!("Input magnetic operations are not in primitive."); - return None; + let mut a = OMatrix::::zeros(3 * db_prim_mag_generators.len()); + let mut b = OVector::::zeros(3 * db_prim_mag_generators.len()); + for (k, mops) in db_prim_mag_generators.iter().enumerate() { + let target_translation = + hm_translations.get(&(mops.operation.rotation, mops.time_reversal))?; + + let ak = mops.operation.rotation - Rotation::identity(); + let bk = mops.operation.translation - target_translation; + for i in 0..3 { + for j in 0..3 { + a[(3 * k + i, j)] = ak[(i, j)]; + } + b[3 * k + i] = bk[i]; + } + } + + 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, } - xsg.push(mops.operation.clone()); - visited.insert(mops.operation.rotation.clone()); } +} + +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); - if (xsg.len() != prim_mag_operations.len()) && (xsg.len() * 2 != prim_mag_operations.len()) { + if (prim_mag_operations.len() % prim_xsg.len() != 0) + || (prim_mag_operations.len() % fsg.len() != 0) + { debug!("Input magnetic operations are incomplete."); return None; } - Some(xsg) + + 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, + ) { + (1, false) => ConstructType::Type1, + (2, true) => ConstructType::Type2, + (2, false) => { + // Find coset representatives of MSG/XSG + if prim_mag_operations + .iter() + .any(|mops| mops.time_reversal && mops.operation.rotation == identity) + { + // Anti-translation + ConstructType::Type4 + } else { + ConstructType::Type3 + } + } + _ => { + debug!( + "Unreachable combination: |MSG/XSG|={}, |FSG/XSG|={}", + prim_mag_operations.len() / prim_xsg.len(), + fsg.len() / prim_xsg.len(), + ); + return None; + } + }; + + // BNS setting + let ref_spg = if construct_type == ConstructType::Type4 { + prim_xsg + } else { + // For type I, II, III, `fsg` is in primitive + fsg + }; + Some((ref_spg, construct_type)) } -fn family_space_group_from_magnetic_space_group( +/// XSG: take only operations without time-reversal +fn maximal_space_subgroup_from_magnetic_space_group( prim_mag_operations: &MagneticOperations, -) -> Option { - let mut fsg = vec![]; - let mut visited = HashSet::new(); - +) -> Operations { + let mut xsg = vec![]; for mops in prim_mag_operations { - if visited.contains(&mops.operation.rotation) { + if mops.time_reversal { continue; } - fsg.push(mops.operation.clone()); - visited.insert(mops.operation.rotation.clone()); + xsg.push(mops.operation.clone()); } + xsg +} - if (fsg.len() != prim_mag_operations.len()) && (fsg.len() * 2 != prim_mag_operations.len()) { - debug!("Input magnetic operations are incomplete."); - return None; +/// FSG: take all operations ignoring time-reversal parts +/// Returned operations may contain duplicated rotation parts (for type-IV). +fn family_space_group_from_magnetic_space_group( + prim_mag_operations: &MagneticOperations, +) -> Operations { + let mut fsg = vec![]; + for mops in prim_mag_operations { + fsg.push(mops.operation.clone()); } - Some(fsg) + fsg } #[cfg(test)] mod tests { + use rstest::rstest; + use test_log::test as test_with_log; + use super::*; use crate::base::Transformation; - use crate::data::MagneticHallSymbol; + use crate::data::{MagneticHallSymbol, NUM_MAGNETIC_SPACE_GROUP_TYPES}; - #[test] - fn test_type4() { - // -P 4 2 3 1abc' (UNI No. 1599) - // P_I m-3m (BNS No. 221.97) - let mhs = MagneticHallSymbol::from_uni_number(1599).unwrap(); + fn get_prim_mag_operations(uni_number: UNINumber) -> MagneticOperations { + let mhs = MagneticHallSymbol::from_uni_number(uni_number).unwrap(); let magnetic_operations = mhs.traverse(); // conventional -> primitive let prim_mag_operations = Transformation::from_linear(mhs.centering.linear()) .inverse_transform_magnetic_operations(&magnetic_operations); + prim_mag_operations + } + + #[rstest] + #[case(2, ConstructType::Type2, 2, 1, 2)] + #[case(1594, ConstructType::Type1, 48, 48, 48)] + #[case(1595, ConstructType::Type2, 96, 48, 96)] + #[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( + #[case] uni_number: UNINumber, + #[case] construct_type: ConstructType, + #[case] order_msg: usize, + #[case] order_xsg: usize, + #[case] order_fsg: usize, + ) { + 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); + assert_eq!(xsg.len(), order_xsg); - let magnetic_space_group = MagneticSpaceGroup::new(&prim_mag_operations, 1e-8).unwrap(); + let fsg = family_space_group_from_magnetic_space_group(&prim_mag_operations); + assert_eq!(fsg.len(), order_fsg); + + let (_, construct_type_actual) = + identify_reference_space_group(&prim_mag_operations, 1e-8).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 { + 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(); + assert_eq!(magnetic_space_group.uni_number, uni_number as UNINumber); + } } } diff --git a/moyo/src/identify/normalizer.rs b/moyo/src/identify/normalizer.rs index 5303447..c3c5932 100644 --- a/moyo/src/identify/normalizer.rs +++ b/moyo/src/identify/normalizer.rs @@ -34,10 +34,9 @@ static POINT_GROUP_NORMALIZERS: Lazy>> = Lazy::new(|| .collect() }); -#[allow(dead_code)] /// Generate integral normalizer of the given point 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. -fn integral_normalizer( +pub fn integral_normalizer( prim_rotations: &Rotations, prim_generators: &Rotations, ) -> Vec { diff --git a/moyo/src/identify/space_group.rs b/moyo/src/identify/space_group.rs index 0ad35e2..1608556 100644 --- a/moyo/src/identify/space_group.rs +++ b/moyo/src/identify/space_group.rs @@ -203,7 +203,6 @@ 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, @@ -211,7 +210,7 @@ fn match_origin_shift( } /// Solve a * x = b (mod 1) -fn solve_mod1( +pub fn solve_mod1( a: &OMatrix, b: &OVector, epsilon: f64, From 750e629c2b1ae0b3ef385b0d205a2e44d1e5eac9 Mon Sep 17 00:00:00 2001 From: lan496 Date: Sat, 4 Jan 2025 18:55:16 +0900 Subject: [PATCH 06/11] Point group normalizer is not enough... --- moyo/src/identify.rs | 1 + moyo/src/identify/magnetic_space_group.rs | 132 ++++++++++------------ moyo/src/identify/normalizer.rs | 99 ++++------------ moyo/src/identify/space_group.rs | 2 +- moyo/src/lib.rs | 63 ++++++++++- moyo/src/search.rs | 2 +- 6 files changed, 146 insertions(+), 153 deletions(-) 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}; From bcd2e78d9592ba30c0a55da16cd11f440744b080 Mon Sep 17 00:00:00 2001 From: lan496 Date: Sat, 4 Jan 2025 19:35:53 +0900 Subject: [PATCH 07/11] iterator for trans_mat --- moyo/src/identify/magnetic_space_group.rs | 11 +- moyo/src/identify/normalizer.rs | 77 +++------ moyo/src/identify/point_group.rs | 197 ++++++++++------------ 3 files changed, 118 insertions(+), 167 deletions(-) diff --git a/moyo/src/identify/magnetic_space_group.rs b/moyo/src/identify/magnetic_space_group.rs index 142266f..536a336 100644 --- a/moyo/src/identify/magnetic_space_group.rs +++ b/moyo/src/identify/magnetic_space_group.rs @@ -302,11 +302,10 @@ mod tests { #[test_with_log] fn test_identify_magnetic_space_group() { // 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(); - assert_eq!(magnetic_space_group.uni_number, uni_number as UNINumber); - } + // 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(); + // assert_eq!(magnetic_space_group.uni_number, uni_number as UNINumber); + // } } } diff --git a/moyo/src/identify/normalizer.rs b/moyo/src/identify/normalizer.rs index 17c6e3f..beb15ae 100644 --- a/moyo/src/identify/normalizer.rs +++ b/moyo/src/identify/normalizer.rs @@ -1,9 +1,9 @@ use itertools::Itertools; +use super::point_group::iter_trans_mat_basis; use super::rotation_type::identify_rotation_type; 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. @@ -23,63 +23,32 @@ pub fn integral_normalizer( .iter() .map(identify_rotation_type) .collect::>(); - let prim_generator_rotation_types = prim_rotation_generators - .iter() - .map(identify_rotation_type) - .collect::>(); - - // Try to map generators - let order = prim_rotations.len(); - let candidates: Vec> = prim_generator_rotation_types - .iter() - .map(|&rotation_type| { - (0..order) - .filter(|&i| rotation_types[i] == rotation_type) - .collect() - }) - .collect(); - // TODO: unify with match_with_point_group let mut conjugators = vec![]; - for pivot in candidates - .iter() - .map(|v| v.iter()) - .multi_cartesian_product() + for trans_mat_basis in + iter_trans_mat_basis(prim_rotations, rotation_types, prim_rotation_generators) { - // 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_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 { - 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; + // 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 { + 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; } } } diff --git a/moyo/src/identify/point_group.rs b/moyo/src/identify/point_group.rs index c148fe6..4f6150a 100644 --- a/moyo/src/identify/point_group.rs +++ b/moyo/src/identify/point_group.rs @@ -2,6 +2,7 @@ use std::cmp::Ordering; use itertools::Itertools; use log::debug; +use nalgebra::Matrix3; use super::rotation_type::{identify_rotation_type, RotationType}; use crate::base::{MoyoError, Rotations, UnimodularLinear}; @@ -71,72 +72,45 @@ fn match_with_cubic_point_group( .iter() .find(|(_, point_group_db)| point_group_db.centering == Centering::P) .unwrap(); - - let order = prim_rotations.len(); let other_prim_generators = primitive_arithmetic_crystal_class.1.primitive_generators(); - // Try to map generators - let other_generator_rotation_types = other_prim_generators - .iter() - .map(identify_rotation_type) - .collect::>(); - let candidates: Vec> = other_generator_rotation_types - .iter() - .map(|&rotation_type| { - (0..order) - .filter(|&i| rotation_types[i] == rotation_type) - .collect() - }) - .collect(); + for trans_mat_basis in iter_trans_mat_basis( + prim_rotations.clone(), + rotation_types.to_vec(), + other_prim_generators, + ) { + // conv_trans_mat: self -> P-centering (primitive) + // The dimension of linear integer system should be one for cubic. + assert_eq!(trans_mat_basis.len(), 1); + let mut conv_trans_mat = trans_mat_basis[0].map(|e| e as f64); - for pivot in candidates - .iter() - .map(|v| v.iter()) - .multi_cartesian_product() - { - // Solve P^-1 * self.rotations[self.rotations[pivot[i]]] * P = other.generators[i] (for all i) - if let Some(trans_mat_basis) = sylvester3( - &pivot - .iter() - .map(|&i| prim_rotations[*i]) - .collect::>(), - &other_prim_generators, - ) { - // conv_trans_mat: self -> P-centering (primitive) - // The dimension of linear integer system should be one for cubic. - assert_eq!(trans_mat_basis.len(), 1); - let mut conv_trans_mat = trans_mat_basis[0].map(|e| e as f64); - - // Guarantee det > 0 - let mut det = conv_trans_mat.determinant().round() as i32; - match det.cmp(&0) { - Ordering::Less => { - conv_trans_mat *= -1.0; - det *= -1; - } - Ordering::Equal => continue, - Ordering::Greater => {} + // Guarantee det > 0 + let mut det = conv_trans_mat.determinant().round() as i32; + match det.cmp(&0) { + Ordering::Less => { + conv_trans_mat *= -1.0; + det *= -1; } + Ordering::Equal => continue, + Ordering::Greater => {} + } - for (arithmetic_crystal_class, point_group_db) in - arithmetic_crystal_class_candidates.iter() - { - let centering = point_group_db.centering; - if centering.order() as i32 != det { - continue; - } - // conv_trans_mat: self -> conventional - let prim_trans_mat = - (conv_trans_mat * centering.inverse()).map(|e| e.round() as i32); - if prim_trans_mat.map(|e| e as f64).determinant().round() as i32 != 1 { - return Err(MoyoError::ArithmeticCrystalClassIdentificationError); - } - - return Ok(PointGroup { - arithmetic_number: *arithmetic_crystal_class, - prim_trans_mat, - }); + for (arithmetic_crystal_class, point_group_db) in arithmetic_crystal_class_candidates.iter() + { + let centering = point_group_db.centering; + if centering.order() as i32 != det { + continue; + } + // conv_trans_mat: self -> conventional + let prim_trans_mat = (conv_trans_mat * centering.inverse()).map(|e| e.round() as i32); + if prim_trans_mat.map(|e| e as f64).determinant().round() as i32 != 1 { + return Err(MoyoError::ArithmeticCrystalClassIdentificationError); } + + return Ok(PointGroup { + arithmetic_number: *arithmetic_crystal_class, + prim_trans_mat, + }); } } @@ -148,8 +122,6 @@ fn match_with_point_group( rotation_types: &[RotationType], geometric_crystal_class: GeometricCrystalClass, ) -> Result { - let order = prim_rotations.len(); - for entry in iter_arithmetic_crystal_entry() { if entry.geometric_crystal_class != geometric_crystal_class { continue; @@ -159,51 +131,28 @@ fn match_with_point_group( PointGroupRepresentative::from_arithmetic_crystal_class(entry.arithmetic_number); let other_prim_generators = point_group_db.primitive_generators(); - // Try to map generators - let other_generator_rotation_types = other_prim_generators - .iter() - .map(identify_rotation_type) - .collect::>(); - let candidates: Vec> = other_generator_rotation_types - .iter() - .map(|&rotation_type| { - (0..order) - .filter(|&i| rotation_types[i] == rotation_type) - .collect() - }) - .collect(); - - for pivot in candidates - .iter() - .map(|v| v.iter()) - .multi_cartesian_product() - { - // Solve P^-1 * self.rotations[self.rotations[pivot[i]]] * P = other.generators[i] (for all i) - if let Some(trans_mat_basis) = sylvester3( - &pivot - .iter() - .map(|&i| prim_rotations[*i]) - .collect::>(), - &other_prim_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() - { - // prim_trans_mat: self -> DB(primitive) - 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 == 1 { - return Ok(PointGroup { - arithmetic_number: entry.arithmetic_number, - prim_trans_mat, - }); - } + for trans_mat_basis in iter_trans_mat_basis( + prim_rotations.clone(), + rotation_types.to_vec(), + other_prim_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() + { + // prim_trans_mat: self -> DB(primitive) + 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 == 1 { + return Ok(PointGroup { + arithmetic_number: entry.arithmetic_number, + prim_trans_mat, + }); } } } @@ -283,6 +232,40 @@ fn identify_geometric_crystal_class( } } +/// Return iterator of basis for transformation matrices that map the point group `prim_rotations` to a point group formed by `other_prim_rotation_generators`. +pub fn iter_trans_mat_basis( + prim_rotations: Rotations, + rotation_types: Vec, + other_prim_rotation_generators: Rotations, +) -> impl Iterator>> { + let other_prim_generator_rotation_types = other_prim_rotation_generators + .iter() + .map(identify_rotation_type) + .collect::>(); + + let order = prim_rotations.len(); + let candidates: Vec> = other_prim_generator_rotation_types + .iter() + .map(|&rotation_type| { + (0..order) + .filter(|&i| rotation_types[i] == rotation_type) + .collect() + }) + .collect(); + + candidates + .into_iter() + .map(|v| v.into_iter()) + .multi_cartesian_product() + .filter_map(move |pivot| { + // Solve P^-1 * prim_rotations[pivot[i]] * P = prim_rotation_generators[i] (for all i) + sylvester3( + &pivot.iter().map(|&i| prim_rotations[i]).collect::>(), + &other_prim_rotation_generators, + ) + }) +} + #[cfg(test)] mod tests { use std::collections::HashSet; From fd887a12a4903cae80bce8a9953bf5649068ea1a Mon Sep 17 00:00:00 2001 From: lan496 Date: Sun, 5 Jan 2025 10:46:12 +0900 Subject: [PATCH 08/11] Get generators of reference space group --- moyo/src/data.rs | 1 + .../src/data/magnetic_hall_symbol_database.rs | 21 +++- moyo/src/data/setting.rs | 75 +++++------ moyo/src/identify/magnetic_space_group.rs | 117 ++++++++++++------ 4 files changed, 137 insertions(+), 77 deletions(-) diff --git a/moyo/src/data.rs b/moyo/src/data.rs index 8f2a796..61400c0 100644 --- a/moyo/src/data.rs +++ b/moyo/src/data.rs @@ -13,6 +13,7 @@ pub use arithmetic_crystal_class::ArithmeticNumber; pub use centering::Centering; pub use hall_symbol::{HallSymbol, MagneticHallSymbol}; pub use hall_symbol_database::{hall_symbol_entry, HallNumber, HallSymbolEntry, Number}; +pub use magnetic_hall_symbol_database::{magnetic_hall_symbol_entry, MagneticHallSymbolEntry}; pub use magnetic_space_group::{ get_magnetic_space_group_type, ConstructType, UNINumber, NUM_MAGNETIC_SPACE_GROUP_TYPES, }; diff --git a/moyo/src/data/magnetic_hall_symbol_database.rs b/moyo/src/data/magnetic_hall_symbol_database.rs index 5bfbb79..d4a251c 100644 --- a/moyo/src/data/magnetic_hall_symbol_database.rs +++ b/moyo/src/data/magnetic_hall_symbol_database.rs @@ -1,4 +1,8 @@ -use super::magnetic_space_group::{UNINumber, NUM_MAGNETIC_SPACE_GROUP_TYPES}; +use super::hall_symbol_database::HallNumber; +use super::magnetic_space_group::{ + get_magnetic_space_group_type, ConstructType, UNINumber, NUM_MAGNETIC_SPACE_GROUP_TYPES, +}; +use super::setting::Setting; #[derive(Debug, Clone)] pub struct MagneticHallSymbolEntry { @@ -13,6 +17,19 @@ impl MagneticHallSymbolEntry { uni_number, } } + + pub fn construct_type(&self) -> ConstructType { + get_magnetic_space_group_type(self.uni_number) + .unwrap() + .construct_type + } + + pub fn reference_hall_number(&self) -> HallNumber { + let number = get_magnetic_space_group_type(self.uni_number) + .unwrap() + .number; + Setting::Standard.hall_number(number).unwrap() + } } pub fn magnetic_hall_symbol_entry(uni_number: UNINumber) -> Option { @@ -1679,7 +1696,7 @@ const MAGNETIC_HALL_SYMBOL_DATABASE: [MagneticHallSymbolEntry; NUM_MAGNETIC_SPAC #[cfg(test)] mod tests { use super::*; - use crate::data::hall_symbol::MagneticHallSymbol; + use crate::data::MagneticHallSymbol; fn iter_magnetic_hall_symbol_entry() -> impl Iterator { MAGNETIC_HALL_SYMBOL_DATABASE.iter() diff --git a/moyo/src/data/setting.rs b/moyo/src/data/setting.rs index 84c67f1..23dc4fa 100644 --- a/moyo/src/data/setting.rs +++ b/moyo/src/data/setting.rs @@ -1,4 +1,4 @@ -use super::hall_symbol_database::HallNumber; +use super::hall_symbol_database::{HallNumber, Number}; #[derive(Debug, Copy, Clone, PartialEq)] /// Preference for the setting of the space group. @@ -11,44 +11,49 @@ pub enum Setting { Standard, } +const SPGLIB_HALL_NUMBERS: [HallNumber; 230] = [ + 1, 2, 3, 6, 9, 18, 21, 30, 39, 57, 60, 63, 72, 81, 90, 108, 109, 112, 115, 116, 119, 122, 123, + 124, 125, 128, 134, 137, 143, 149, 155, 161, 164, 170, 173, 176, 182, 185, 191, 197, 203, 209, + 212, 215, 218, 221, 227, 228, 230, 233, 239, 245, 251, 257, 263, 266, 269, 275, 278, 284, 290, + 292, 298, 304, 310, 313, 316, 322, 334, 335, 337, 338, 341, 343, 349, 350, 351, 352, 353, 354, + 355, 356, 357, 358, 359, 361, 363, 364, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, + 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, + 396, 397, 398, 399, 400, 401, 402, 404, 406, 407, 408, 410, 412, 413, 414, 416, 418, 419, 420, + 422, 424, 425, 426, 428, 430, 431, 432, 433, 435, 436, 438, 439, 440, 441, 442, 443, 444, 446, + 447, 448, 449, 450, 452, 454, 455, 456, 457, 458, 460, 462, 463, 464, 465, 466, 467, 468, 469, + 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, + 489, 490, 491, 492, 493, 494, 495, 497, 498, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, + 510, 511, 512, 513, 514, 515, 516, 517, 518, 520, 521, 523, 524, 525, 527, 529, 530, +]; +const STANDARD_HALL_NUMBERS: [HallNumber; 230] = [ + 1, 2, 3, 6, 9, 18, 21, 30, 39, 57, 60, 63, 72, 81, 90, 108, 109, 112, 115, 116, 119, 122, 123, + 124, 125, 128, 134, 137, 143, 149, 155, 161, 164, 170, 173, 176, 182, 185, 191, 197, 203, 209, + 212, 215, 218, 221, 227, 229, 230, 234, 239, 245, 251, 257, 263, 266, 269, 275, 279, 284, 290, + 292, 298, 304, 310, 313, 316, 323, 334, 336, 337, 338, 341, 343, 349, 350, 351, 352, 353, 354, + 355, 356, 357, 358, 360, 362, 363, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, + 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, + 396, 397, 398, 399, 400, 401, 403, 405, 406, 407, 409, 411, 412, 413, 415, 417, 418, 419, 421, + 423, 424, 425, 427, 429, 430, 431, 432, 433, 435, 436, 438, 439, 440, 441, 442, 443, 444, 446, + 447, 448, 449, 450, 452, 454, 455, 456, 457, 458, 460, 462, 463, 464, 465, 466, 467, 468, 469, + 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, + 489, 490, 491, 492, 493, 494, 496, 497, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, + 510, 511, 512, 513, 514, 515, 516, 517, 519, 520, 522, 523, 524, 526, 528, 529, 530, +]; + impl Setting { pub fn hall_numbers(&self) -> Vec { match self { Setting::HallNumber(hall_number) => vec![*hall_number], - Setting::Spglib => vec![ - 1, 2, 3, 6, 9, 18, 21, 30, 39, 57, 60, 63, 72, 81, 90, 108, 109, 112, 115, 116, - 119, 122, 123, 124, 125, 128, 134, 137, 143, 149, 155, 161, 164, 170, 173, 176, - 182, 185, 191, 197, 203, 209, 212, 215, 218, 221, 227, 228, 230, 233, 239, 245, - 251, 257, 263, 266, 269, 275, 278, 284, 290, 292, 298, 304, 310, 313, 316, 322, - 334, 335, 337, 338, 341, 343, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, - 359, 361, 363, 364, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, - 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, - 394, 395, 396, 397, 398, 399, 400, 401, 402, 404, 406, 407, 408, 410, 412, 413, - 414, 416, 418, 419, 420, 422, 424, 425, 426, 428, 430, 431, 432, 433, 435, 436, - 438, 439, 440, 441, 442, 443, 444, 446, 447, 448, 449, 450, 452, 454, 455, 456, - 457, 458, 460, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, - 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, - 491, 492, 493, 494, 495, 497, 498, 500, 501, 502, 503, 504, 505, 506, 507, 508, - 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 520, 521, 523, 524, 525, 527, - 529, 530, - ], - Setting::Standard => vec![ - 1, 2, 3, 6, 9, 18, 21, 30, 39, 57, 60, 63, 72, 81, 90, 108, 109, 112, 115, 116, - 119, 122, 123, 124, 125, 128, 134, 137, 143, 149, 155, 161, 164, 170, 173, 176, - 182, 185, 191, 197, 203, 209, 212, 215, 218, 221, 227, 229, 230, 234, 239, 245, - 251, 257, 263, 266, 269, 275, 279, 284, 290, 292, 298, 304, 310, 313, 316, 323, - 334, 336, 337, 338, 341, 343, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, - 360, 362, 363, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, - 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, - 394, 395, 396, 397, 398, 399, 400, 401, 403, 405, 406, 407, 409, 411, 412, 413, - 415, 417, 418, 419, 421, 423, 424, 425, 427, 429, 430, 431, 432, 433, 435, 436, - 438, 439, 440, 441, 442, 443, 444, 446, 447, 448, 449, 450, 452, 454, 455, 456, - 457, 458, 460, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, - 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, - 491, 492, 493, 494, 496, 497, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, - 509, 510, 511, 512, 513, 514, 515, 516, 517, 519, 520, 522, 523, 524, 526, 528, - 529, 530, - ], + Setting::Spglib => SPGLIB_HALL_NUMBERS.to_vec(), + Setting::Standard => STANDARD_HALL_NUMBERS.to_vec(), + } + } + + pub fn hall_number(&self, number: Number) -> Option { + match self { + Setting::HallNumber(_) => None, + Setting::Spglib => SPGLIB_HALL_NUMBERS.get(number as usize - 1).cloned(), + Setting::Standard => STANDARD_HALL_NUMBERS.get(number as usize - 1).cloned(), } } } diff --git a/moyo/src/identify/magnetic_space_group.rs b/moyo/src/identify/magnetic_space_group.rs index 536a336..5357cba 100644 --- a/moyo/src/identify/magnetic_space_group.rs +++ b/moyo/src/identify/magnetic_space_group.rs @@ -9,8 +9,8 @@ use crate::base::{ UnimodularTransformation, }; use crate::data::{ - get_magnetic_space_group_type, uni_number_range, ConstructType, MagneticHallSymbol, Setting, - UNINumber, + get_magnetic_space_group_type, hall_symbol_entry, magnetic_hall_symbol_entry, uni_number_range, + ConstructType, HallSymbol, MagneticHallSymbol, MagneticHallSymbolEntry, Setting, UNINumber, }; #[derive(Debug)] @@ -54,11 +54,11 @@ impl MagneticSpaceGroup { // TODO: - let mhs = MagneticHallSymbol::from_uni_number(uni_number) + let entry = magnetic_hall_symbol_entry(uni_number).unwrap(); + let mhs = MagneticHallSymbol::new(&entry.magnetic_hall_symbol) .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); + let db_ref_prim_generators = db_reference_space_group_primitive_generators(&entry); // The correction transformations keep the reference space group of `tmp_prim_mag_operations` // TODO: precompute the correction transformation matrices @@ -83,40 +83,6 @@ impl MagneticSpaceGroup { Err(MoyoError::MagneticSpaceGroupTypeIdentificationError) } - fn get_db_reference_space_group_primitive_generators( - mhs: &MagneticHallSymbol, - construct_type: &ConstructType, - ) -> Operations { - let db_prim_mag_operations = mhs.primitive_generators(); - match construct_type { - ConstructType::Type1 | ConstructType::Type2 | ConstructType::Type3 => { - // Reference space group: FSG - // -> Remove time-reversal parts - db_prim_mag_operations - .iter() - .map(|mops| mops.operation.clone()) - .collect() - } - ConstructType::Type4 => { - // Reference space group: XSG - // -> Remove anti-translation operation - let identity = Rotation::identity(); - db_prim_mag_operations - .iter() - .filter_map(|mops| { - // Here we assume the magnetic Hall symbol composes of XSG generators and one anti-translation operation - if mops.time_reversal { - assert_eq!(mops.operation.rotation, identity); - None - } else { - Some(mops.operation.clone()) - } - }) - .collect() - } - } - } - /// 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( @@ -252,6 +218,20 @@ fn family_space_group_from_magnetic_space_group( (fsg, is_type2) } +/// Return generators of the reference space group of magnetic space group `entry`. +/// This function assumes the magnetic Hall symbol is extended from the Hall symbol in the standard setting. +fn db_reference_space_group_primitive_generators(entry: &MagneticHallSymbolEntry) -> Operations { + let ref_hall_number = entry.reference_hall_number(); + let ref_hall_entry = hall_symbol_entry(ref_hall_number).unwrap(); + let identity = Rotation::identity(); + HallSymbol::new(&ref_hall_entry.hall_symbol) + .unwrap() + .primitive_generators() + .into_iter() + .filter(|ops| ops.rotation != identity) // In primitive, if rotation part is identity, it is a pure translation + .collect() +} + #[cfg(test)] mod tests { use rstest::rstest; @@ -259,7 +239,9 @@ mod tests { use super::*; use crate::base::Transformation; - use crate::data::{MagneticHallSymbol, NUM_MAGNETIC_SPACE_GROUP_TYPES}; + use crate::data::{ + magnetic_hall_symbol_entry, MagneticHallSymbol, NUM_MAGNETIC_SPACE_GROUP_TYPES, + }; fn get_prim_mag_operations(uni_number: UNINumber) -> MagneticOperations { let mhs = MagneticHallSymbol::from_uni_number(uni_number).unwrap(); @@ -299,6 +281,61 @@ mod tests { assert_eq!(construct_type_actual, construct_type); } + // Check generators of reference space group by two methods: + // 1. From the magnetic Hall symbol + // 2. From the Hall symbol with the corresponding Hall number + #[test_with_log] + fn test_db_reference_space_group_primitive_generators() { + for uni_number in 1..=NUM_MAGNETIC_SPACE_GROUP_TYPES { + let entry = magnetic_hall_symbol_entry(uni_number as UNINumber).unwrap(); + let actual = db_reference_space_group_primitive_generators(&entry); + + let mhs = MagneticHallSymbol::new(&entry.magnetic_hall_symbol).unwrap(); + let identity = Rotation::identity(); + let expect: Operations = match entry.construct_type() { + ConstructType::Type1 | ConstructType::Type2 => mhs + .primitive_generators() + .iter() + .filter_map(|mops| { + if mops.time_reversal || mops.operation.rotation == identity { + // Ignore 1' for Type2 + None + } else { + Some(mops.operation.clone()) + } + }) + .collect(), + ConstructType::Type3 => mhs + .primitive_generators() + .iter() + .map(|mops| mops.operation.clone()) // Ignore time-reversal parts + .collect(), + ConstructType::Type4 => mhs + .primitive_generators() + .iter() + .filter_map(|mops| { + if mops.operation.rotation == identity { + // Ignore anti-translation + None + } else { + Some(mops.operation.clone()) + } + }) + .collect(), + }; + assert_eq!(actual.len(), expect.len()); + let mut hm_translation = HashMap::new(); + for ops1 in actual.iter() { + hm_translation.insert(ops1.rotation.clone(), ops1.translation); + } + for ops2 in expect.iter() { + let translation1 = hm_translation.get(&ops2.rotation).unwrap(); + let diff = ops2.translation - translation1; + assert_relative_eq!(diff.map(|e| (e - e.round().abs())).max(), 0.0); + } + } + } + #[test_with_log] fn test_identify_magnetic_space_group() { // for uni_number in 1..=NUM_MAGNETIC_SPACE_GROUP_TYPES { From 1c15e79681507a5a2bb8d1b3e6956fa505d67b3a Mon Sep 17 00:00:00 2001 From: lan496 Date: Sun, 5 Jan 2025 11:25:57 +0900 Subject: [PATCH 09/11] traverse in primitive cell --- moyo/src/base/transformation.rs | 13 +++ moyo/src/data/hall_symbol.rs | 12 +++ moyo/src/identify/magnetic_space_group.rs | 118 ++++++---------------- moyo/src/identify/space_group.rs | 19 +--- 4 files changed, 61 insertions(+), 101 deletions(-) diff --git a/moyo/src/base/transformation.rs b/moyo/src/base/transformation.rs index aef1e92..ca1eff2 100644 --- a/moyo/src/base/transformation.rs +++ b/moyo/src/base/transformation.rs @@ -1,3 +1,5 @@ +use std::ops::Mul; + use itertools::iproduct; use nalgebra::base::{Matrix3, Vector3}; @@ -116,6 +118,17 @@ impl UnimodularTransformation { } } +impl Mul for UnimodularTransformation { + type Output = Self; + + // (P_lhs, p_lhs) * (P_rhs, p_rhs) = (P_lhs * P_rhs, P_lhs * p_rhs + p_lhs) + fn mul(self, rhs: Self) -> Self::Output { + let new_linear = self.linear * rhs.linear; + let new_origin_shift = self.linear.map(|e| e as f64) * rhs.origin_shift + self.origin_shift; + Self::new(new_linear, new_origin_shift) + } +} + /// Represent change of origin and basis for an affine space #[derive(Debug, Clone)] pub struct Transformation { diff --git a/moyo/src/data/hall_symbol.rs b/moyo/src/data/hall_symbol.rs index 708ed7b..d1805ec 100644 --- a/moyo/src/data/hall_symbol.rs +++ b/moyo/src/data/hall_symbol.rs @@ -104,6 +104,12 @@ impl HallSymbol { operations } + pub fn primitive_traverse(&self) -> Operations { + let operations = self.traverse(); + Transformation::from_linear(self.centering.linear()) + .inverse_transform_operations(&operations) + } + pub fn from_hall_number(hall_number: HallNumber) -> Option { if let Some(entry) = hall_symbol_entry(hall_number) { Self::new(entry.hall_symbol) @@ -209,6 +215,12 @@ impl MagneticHallSymbol { operations } + pub fn primitive_traverse(&self) -> MagneticOperations { + let operations = self.traverse(); + Transformation::from_linear(self.centering.linear()) + .inverse_transform_magnetic_operations(&operations) + } + pub fn from_uni_number(uni_number: UNINumber) -> Option { if let Some(entry) = magnetic_hall_symbol_entry(uni_number) { Self::new(entry.magnetic_hall_symbol) diff --git a/moyo/src/identify/magnetic_space_group.rs b/moyo/src/identify/magnetic_space_group.rs index 5357cba..5c17285 100644 --- a/moyo/src/identify/magnetic_space_group.rs +++ b/moyo/src/identify/magnetic_space_group.rs @@ -1,12 +1,11 @@ use std::collections::HashMap; use log::debug; -use nalgebra::{Dyn, OMatrix, OVector, U3}; -use super::space_group::{solve_mod1, SpaceGroup}; +use super::normalizer::integral_normalizer; +use super::space_group::SpaceGroup; use crate::base::{ - MagneticOperations, MoyoError, Operations, Rotation, Translation, UnimodularLinear, - UnimodularTransformation, + MagneticOperations, MoyoError, Operations, Rotation, Translation, UnimodularTransformation, }; use crate::data::{ get_magnetic_space_group_type, hall_symbol_entry, magnetic_hall_symbol_entry, uni_number_range, @@ -52,78 +51,30 @@ impl MagneticSpaceGroup { }); } - // TODO: - let entry = magnetic_hall_symbol_entry(uni_number).unwrap(); let mhs = MagneticHallSymbol::new(&entry.magnetic_hall_symbol) .ok_or(MoyoError::MagneticSpaceGroupTypeIdentificationError)?; - let db_prim_mag_generators = mhs.primitive_generators(); - let db_ref_prim_generators = db_reference_space_group_primitive_generators(&entry); - - // The correction transformations keep the reference space group of `tmp_prim_mag_operations` - // TODO: precompute the correction transformation matrices - 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) - } - - /// 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( - prim_mag_operations: &MagneticOperations, - trans_mat: &UnimodularLinear, - db_prim_mag_generators: &MagneticOperations, - epsilon: f64, - ) -> Option { - let new_prim_mag_operations = UnimodularTransformation::from_linear(*trans_mat) - .transform_magnetic_operations(prim_mag_operations); - let mut hm_translations = HashMap::new(); - for mops in new_prim_mag_operations.iter() { - hm_translations.insert( - (mops.operation.rotation, mops.time_reversal), - mops.operation.translation, - ); - } - - let mut a = OMatrix::::zeros(3 * db_prim_mag_generators.len()); - let mut b = OVector::::zeros(3 * db_prim_mag_generators.len()); - for (k, mops) in db_prim_mag_generators.iter().enumerate() { - let target_translation = - hm_translations.get(&(mops.operation.rotation, mops.time_reversal))?; - - let ak = mops.operation.rotation - Rotation::identity(); - let bk = mops.operation.translation - target_translation; - for i in 0..3 { - for j in 0..3 { - a[(3 * k + i, j)] = ak[(i, j)]; - } - b[3 * k + i] = bk[i]; + let db_prim_mag_operations = mhs.primitive_traverse(); + let (db_ref_prim_operations, db_ref_prim_generators) = + db_reference_space_group_primitive(&entry); + let db_ref_prim_normalizer = + integral_normalizer(&db_ref_prim_operations, &db_ref_prim_generators, epsilon); // TODO: precompute the normalizer + + for corr_trans in db_ref_prim_normalizer { + // new_transformation * corr_trans: primitive input -> primitive DB + let new_transformation = std_ref_spg.transformation.clone() * corr_trans; + let new_prim_mag_operations = + new_transformation.transform_magnetic_operations(prim_mag_operations); + + // TODO: + // debug!("Matched with UNI number {}", uni_number); + // return Ok(Self { + // uni_number, + // transformation: new_transformation, + // }); } } - - 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, - } + Err(MoyoError::MagneticSpaceGroupTypeIdentificationError) } } @@ -218,18 +169,19 @@ fn family_space_group_from_magnetic_space_group( (fsg, is_type2) } -/// Return generators of the reference space group of magnetic space group `entry`. +/// Return operations and generators of the reference space group of magnetic space group `entry`. /// This function assumes the magnetic Hall symbol is extended from the Hall symbol in the standard setting. -fn db_reference_space_group_primitive_generators(entry: &MagneticHallSymbolEntry) -> Operations { - let ref_hall_number = entry.reference_hall_number(); - let ref_hall_entry = hall_symbol_entry(ref_hall_number).unwrap(); +fn db_reference_space_group_primitive(entry: &MagneticHallSymbolEntry) -> (Operations, Operations) { + let ref_hall_entry = hall_symbol_entry(entry.reference_hall_number()).unwrap(); + let ref_hall_symbol = HallSymbol::new(&ref_hall_entry.hall_symbol).unwrap(); + let ref_prim_operations = ref_hall_symbol.primitive_traverse(); let identity = Rotation::identity(); - HallSymbol::new(&ref_hall_entry.hall_symbol) - .unwrap() + let ref_prim_generators = ref_hall_symbol .primitive_generators() .into_iter() .filter(|ops| ops.rotation != identity) // In primitive, if rotation part is identity, it is a pure translation - .collect() + .collect(); + (ref_prim_operations, ref_prim_generators) } #[cfg(test)] @@ -238,19 +190,13 @@ mod tests { use test_log::test as test_with_log; use super::*; - use crate::base::Transformation; use crate::data::{ magnetic_hall_symbol_entry, MagneticHallSymbol, NUM_MAGNETIC_SPACE_GROUP_TYPES, }; fn get_prim_mag_operations(uni_number: UNINumber) -> MagneticOperations { let mhs = MagneticHallSymbol::from_uni_number(uni_number).unwrap(); - let magnetic_operations = mhs.traverse(); - - // conventional -> primitive - let prim_mag_operations = Transformation::from_linear(mhs.centering.linear()) - .inverse_transform_magnetic_operations(&magnetic_operations); - prim_mag_operations + mhs.primitive_traverse() } #[rstest] @@ -288,7 +234,7 @@ mod tests { fn test_db_reference_space_group_primitive_generators() { for uni_number in 1..=NUM_MAGNETIC_SPACE_GROUP_TYPES { let entry = magnetic_hall_symbol_entry(uni_number as UNINumber).unwrap(); - let actual = db_reference_space_group_primitive_generators(&entry); + let (_, actual) = db_reference_space_group_primitive(&entry); let mhs = MagneticHallSymbol::new(&entry.magnetic_hall_symbol).unwrap(); let identity = Rotation::identity(); diff --git a/moyo/src/identify/space_group.rs b/moyo/src/identify/space_group.rs index a664ef6..7ff20ae 100644 --- a/moyo/src/identify/space_group.rs +++ b/moyo/src/identify/space_group.rs @@ -248,7 +248,7 @@ mod tests { use rstest::rstest; use std::collections::HashMap; - use crate::base::{Transformation, UnimodularTransformation, EPS}; + use crate::base::{UnimodularTransformation, EPS}; use crate::data::{hall_symbol_entry, HallSymbol, Setting}; use super::{correction_transformation_matrices, solve_mod1, SpaceGroup}; @@ -272,11 +272,7 @@ mod tests { fn test_correction_transformation_matrices() { let hall_number = 21; // P 1 c 1 let hall_symbol = HallSymbol::from_hall_number(hall_number).unwrap(); - let operations = hall_symbol.traverse(); - - // conventional -> primitive - let prim_operations = Transformation::from_linear(hall_symbol.centering.linear()) - .inverse_transform_operations(&operations); + let prim_operations = hall_symbol.primitive_traverse(); // The correction transformation matrices should change the group into P1c1, P1a1, and P1n1 let entry = hall_symbol_entry(hall_number).unwrap(); @@ -309,11 +305,7 @@ mod tests { fn test_identify_space_group(#[case] setting: Setting) { for hall_number in 1..=530 { let hall_symbol = HallSymbol::from_hall_number(hall_number).unwrap(); - let operations = hall_symbol.traverse(); - - // conventional -> primitive - let prim_operations = Transformation::from_linear(hall_symbol.centering.linear()) - .inverse_transform_operations(&operations); + let prim_operations = hall_symbol.primitive_traverse(); let space_group = SpaceGroup::new(&prim_operations, setting, 1e-8).unwrap(); @@ -333,10 +325,7 @@ mod tests { let matched_hall_symbol = HallSymbol::from_hall_number(space_group.hall_number).unwrap(); - let matched_operations = matched_hall_symbol.traverse(); - let matched_prim_operations = - Transformation::from_linear(matched_hall_symbol.centering.linear()) - .inverse_transform_operations(&matched_operations); + let matched_prim_operations = matched_hall_symbol.primitive_traverse(); let mut hm_translations = HashMap::new(); for operation in matched_prim_operations.iter() { From 01f4a706aa2bca0bc728a86d906e8b6f8331e57b Mon Sep 17 00:00:00 2001 From: lan496 Date: Sat, 11 Jan 2025 09:10:05 +0900 Subject: [PATCH 10/11] conjugator for type 4 --- moyo/src/identify/magnetic_space_group.rs | 208 +++++++++++++++++++--- moyo/src/identify/normalizer.rs | 2 +- 2 files changed, 185 insertions(+), 25 deletions(-) diff --git a/moyo/src/identify/magnetic_space_group.rs b/moyo/src/identify/magnetic_space_group.rs index 5c17285..99a95bb 100644 --- a/moyo/src/identify/magnetic_space_group.rs +++ b/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, Operation, Operations, Rotation, Translation, + UnimodularLinear, UnimodularTransformation, }; use crate::data::{ get_magnetic_space_group_type, hall_symbol_entry, magnetic_hall_symbol_entry, uni_number_range, @@ -34,6 +38,7 @@ impl MagneticSpaceGroup { let uni_number_range = uni_number_range(std_ref_spg.number) .ok_or(MoyoError::SpaceGroupTypeIdentificationError)?; + for uni_number in uni_number_range { if get_magnetic_space_group_type(uni_number) .unwrap() @@ -57,25 +62,119 @@ impl MagneticSpaceGroup { let db_prim_mag_operations = mhs.primitive_traverse(); let (db_ref_prim_operations, db_ref_prim_generators) = db_reference_space_group_primitive(&entry); - let db_ref_prim_normalizer = - integral_normalizer(&db_ref_prim_operations, &db_ref_prim_generators, epsilon); // TODO: precompute the normalizer - - for corr_trans in db_ref_prim_normalizer { - // new_transformation * corr_trans: primitive input -> primitive DB - let new_transformation = std_ref_spg.transformation.clone() * corr_trans; - let new_prim_mag_operations = - new_transformation.transform_magnetic_operations(prim_mag_operations); - - // TODO: - // debug!("Matched with UNI number {}", uni_number); - // return Ok(Self { - // uni_number, - // transformation: new_transformation, - // }); + + match construct_type { + ConstructType::Type3 => { + // Centralizer of FSG also keep XSG invariant. Thus, we only need to consider the normalizer of FSG up to the centralizer. + // TODO: precompute the normalizer + let db_fsg_prim_normalizer = integral_normalizer( + &db_ref_prim_operations, + &db_ref_prim_generators, + epsilon, + ); + for corr_trans in db_fsg_prim_normalizer { + // new_transformation * corr_trans: primitive input -> primitive DB + let new_transformation = std_ref_spg.transformation.clone() * corr_trans; + let new_prim_mag_operations = + new_transformation.transform_magnetic_operations(prim_mag_operations); + + if Self::match_prim_mag_operations( + &new_prim_mag_operations, + &db_prim_mag_operations, + epsilon, + ) { + debug!("Matched with UNI number {}", uni_number); + return Ok(Self { + uni_number, + transformation: new_transformation, + }); + } + } + } + ConstructType::Type4 => { + // Find conjugator which transforms the anti-translation while keeping XSGs + let identity = Rotation::identity(); + let original_anti_translation = prim_mag_operations + .iter() + .filter_map(|mops| { + if (mops.operation.rotation == identity) && mops.time_reversal { + Some(mops) + } else { + None + } + }) + .nth(0) + .unwrap(); + let src_translation = std_ref_spg + .transformation + .transform_magnetic_operation(original_anti_translation) + .operation + .translation; + // TODO: refactor filtering anti-translation + let dst_translation = db_prim_mag_operations + .iter() + .filter_map(|mops| { + if (mops.operation.rotation == identity) && mops.time_reversal { + Some(mops) + } else { + None + } + }) + .nth(0) + .unwrap() + .operation + .translation; + if let Some(corr_trans) = find_conjugator_type4( + &db_ref_prim_generators, + &db_ref_prim_operations, + &src_translation, + &dst_translation, + epsilon, + ) { + let new_transformation = std_ref_spg.transformation.clone() * corr_trans; + return Ok(Self { + uni_number, + transformation: new_transformation, + }); + } + } + _ => unreachable!(), } } Err(MoyoError::MagneticSpaceGroupTypeIdentificationError) } + + fn match_prim_mag_operations( + prim_mag_operations1: &MagneticOperations, + prim_mag_operations2: &MagneticOperations, + epsilon: f64, + ) -> bool { + if prim_mag_operations1.len() != prim_mag_operations2.len() { + return false; + } + + let mut hm_translation = HashMap::new(); + for mops1 in prim_mag_operations1 { + hm_translation.insert( + (mops1.operation.rotation.clone(), mops1.time_reversal), + mops1.operation.translation, + ); + } + + for mops2 in prim_mag_operations2 { + if let Some(translation1) = + hm_translation.get(&(mops2.operation.rotation.clone(), mops2.time_reversal)) + { + let diff = mops2.operation.translation - translation1; + if !diff.iter().all(|e| (e - e.round()).abs() < epsilon) { + return false; + } + } else { + return false; + } + } + true + } } fn identify_reference_space_group( @@ -176,14 +275,74 @@ fn db_reference_space_group_primitive(entry: &MagneticHallSymbolEntry) -> (Opera let ref_hall_symbol = HallSymbol::new(&ref_hall_entry.hall_symbol).unwrap(); let ref_prim_operations = ref_hall_symbol.primitive_traverse(); let identity = Rotation::identity(); - let ref_prim_generators = ref_hall_symbol + + let mut ref_prim_generators = ref_hall_symbol .primitive_generators() .into_iter() .filter(|ops| ops.rotation != identity) // In primitive, if rotation part is identity, it is a pure translation - .collect(); + .collect::>(); + if ref_prim_generators.is_empty() { + ref_prim_generators.push(Operation::identity()); + } + (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 { + 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::>(); + 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; @@ -285,10 +444,11 @@ mod tests { #[test_with_log] fn test_identify_magnetic_space_group() { // for uni_number in 1..=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(); - // assert_eq!(magnetic_space_group.uni_number, uni_number as UNINumber); - // } + for uni_number in 282..=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(); + assert_eq!(magnetic_space_group.uni_number, uni_number as UNINumber); + } } } diff --git a/moyo/src/identify/normalizer.rs b/moyo/src/identify/normalizer.rs index beb15ae..5cbddfc 100644 --- a/moyo/src/identify/normalizer.rs +++ b/moyo/src/identify/normalizer.rs @@ -47,8 +47,8 @@ pub fn integral_normalizer( match_origin_shift(prim_operations, &prim_trans_mat, prim_generators, epsilon) { conjugators.push(UnimodularTransformation::new(prim_trans_mat, origin_shift)); + break; } - break; } } } From fbc502b69e354cd86de948440160eaf03cddf7ea Mon Sep 17 00:00:00 2001 From: lan496 Date: Fri, 17 Jan 2025 21:25:15 +0900 Subject: [PATCH 11/11] iter for unimodular transformation matrix --- moyo/src/data/hall_symbol.rs | 3 +- moyo/src/identify/magnetic_space_group.rs | 54 ++++++++++------------- moyo/src/identify/normalizer.rs | 29 +++--------- moyo/src/identify/point_group.rs | 46 +++++++++++-------- 4 files changed, 60 insertions(+), 72 deletions(-) diff --git a/moyo/src/data/hall_symbol.rs b/moyo/src/data/hall_symbol.rs index d1805ec..338e870 100644 --- a/moyo/src/data/hall_symbol.rs +++ b/moyo/src/data/hall_symbol.rs @@ -589,6 +589,7 @@ fn purify_translation_mod1(translation: &Translation) -> Translation { mod tests { use nalgebra::{matrix, vector}; use rstest::rstest; + use test_log::test as test_with_log; use super::*; @@ -614,7 +615,7 @@ mod tests { assert_eq!(operations.len(), num_operations); } - #[test] + #[test_with_log] fn test_hall_symbol_generators() { // No. 178 let hs = HallSymbol::new("P 61 2 (0 0 5)").unwrap(); diff --git a/moyo/src/identify/magnetic_space_group.rs b/moyo/src/identify/magnetic_space_group.rs index 99a95bb..c658b8f 100644 --- a/moyo/src/identify/magnetic_space_group.rs +++ b/moyo/src/identify/magnetic_space_group.rs @@ -4,7 +4,7 @@ use itertools::Itertools; use log::debug; use super::normalizer::integral_normalizer; -use super::point_group::iter_trans_mat_basis; +use super::point_group::{iter_trans_mat_basis, iter_unimodular_trans_mat}; use super::rotation_type::identify_rotation_type; use super::space_group::{match_origin_shift, SpaceGroup}; use crate::base::{ @@ -60,6 +60,9 @@ impl MagneticSpaceGroup { let mhs = MagneticHallSymbol::new(&entry.magnetic_hall_symbol) .ok_or(MoyoError::MagneticSpaceGroupTypeIdentificationError)?; let db_prim_mag_operations = mhs.primitive_traverse(); + dbg!(&uni_number); + dbg!(&entry.magnetic_hall_symbol); + dbg!(&mhs.primitive_generators()); let (db_ref_prim_operations, db_ref_prim_generators) = db_reference_space_group_primitive(&entry); @@ -308,35 +311,20 @@ fn find_conjugator_type4( 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; + for prim_trans_mat in iter_unimodular_trans_mat(trans_mat_basis) { + // (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 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)); - } + 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)); } } } @@ -397,7 +385,7 @@ mod tests { let mhs = MagneticHallSymbol::new(&entry.magnetic_hall_symbol).unwrap(); let identity = Rotation::identity(); - let expect: Operations = match entry.construct_type() { + let mut expect: Operations = match entry.construct_type() { ConstructType::Type1 | ConstructType::Type2 => mhs .primitive_generators() .iter() @@ -428,6 +416,10 @@ mod tests { }) .collect(), }; + if expect.is_empty() { + expect.push(Operation::identity()); + } + assert_eq!(actual.len(), expect.len()); let mut hm_translation = HashMap::new(); for ops1 in actual.iter() { @@ -443,8 +435,8 @@ mod tests { #[test_with_log] fn test_identify_magnetic_space_group() { - // for uni_number in 1..=NUM_MAGNETIC_SPACE_GROUP_TYPES { - for uni_number in 282..=NUM_MAGNETIC_SPACE_GROUP_TYPES { + for uni_number in 1..=NUM_MAGNETIC_SPACE_GROUP_TYPES { + // for uni_number in 282..=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 5cbddfc..e9842b6 100644 --- a/moyo/src/identify/normalizer.rs +++ b/moyo/src/identify/normalizer.rs @@ -1,6 +1,6 @@ use itertools::Itertools; -use super::point_group::iter_trans_mat_basis; +use super::point_group::{iter_trans_mat_basis, iter_unimodular_trans_mat}; use super::rotation_type::identify_rotation_type; use super::space_group::match_origin_shift; use crate::base::{project_rotations, Operations, UnimodularLinear, UnimodularTransformation}; @@ -28,27 +28,12 @@ pub fn integral_normalizer( for trans_mat_basis in iter_trans_mat_basis(prim_rotations, rotation_types, 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 { - 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; - } + for prim_trans_mat in iter_unimodular_trans_mat(trans_mat_basis) { + 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; } } } diff --git a/moyo/src/identify/point_group.rs b/moyo/src/identify/point_group.rs index 4f6150a..96949dc 100644 --- a/moyo/src/identify/point_group.rs +++ b/moyo/src/identify/point_group.rs @@ -136,24 +136,11 @@ fn match_with_point_group( rotation_types.to_vec(), other_prim_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() - { - // prim_trans_mat: self -> DB(primitive) - 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 == 1 { - return Ok(PointGroup { - arithmetic_number: entry.arithmetic_number, - prim_trans_mat, - }); - } + if let Some(prim_trans_mat) = iter_unimodular_trans_mat(trans_mat_basis).nth(0) { + return Ok(PointGroup { + arithmetic_number: entry.arithmetic_number, + prim_trans_mat, + }); } } } @@ -266,6 +253,29 @@ pub fn iter_trans_mat_basis( }) } +/// Search integer linear combination such that the transformation matrix is unimodular +/// Consider coefficients in [-2, 2], which will be sufficient for Delaunay reduced basis +pub fn iter_unimodular_trans_mat( + trans_mat_basis: Vec>, +) -> impl Iterator { + (0..trans_mat_basis.len()) + .map(|_| -2..=2) + .multi_cartesian_product() + .filter_map(move |comb| { + // prim_trans_mat: self -> DB(primitive) + 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 == 1 { + Some(prim_trans_mat) + } else { + None + } + }) +} + #[cfg(test)] mod tests { use std::collections::HashSet;