Skip to content

Commit

Permalink
Point group normalizer is not enough...
Browse files Browse the repository at this point in the history
  • Loading branch information
@lan496
lan496 committed Jan 4, 2025
1 parent 53933f5 commit 55e17bd
Show file tree
Hide file tree
Showing 6 changed files with 146 additions and 153 deletions.
1 change: 1 addition & 0 deletions moyo/src/identify.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
132 changes: 61 additions & 71 deletions moyo/src/identify/magnetic_space_group.rs
View file
Original file line number Diff line number Diff line change
@@ -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,
Expand Down Expand Up @@ -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)
}
Expand Down Expand Up @@ -110,27 +117,6 @@ impl MagneticSpaceGroup {
}
}

fn correction_transformation_matrices(
db_ref_prim_generators: &Operations,
construct_type: &ConstructType,
) -> Option<Vec<UnimodularLinear>> {
// 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 <db_prim_mag_generators>
/// TODO: unify with identify/space_group.rs::match_origin_shift
fn match_origin_shift(
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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;
Expand All @@ -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)]
Expand All @@ -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(
Expand All @@ -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();
Expand Down
99 changes: 22 additions & 77 deletions moyo/src/identify/normalizer.rs
View file
Original file line number Diff line number Diff line change
@@ -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<Vec<UnimodularLinear>> {
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<Vec<Vec<UnimodularLinear>>> = 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<UnimodularLinear> {
prim_operations: &Operations,
prim_generators: &Operations,
epsilon: f64,
) -> Vec<UnimodularTransformation> {
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::<Vec<_>>();
let prim_generator_rotation_types = prim_generators
let prim_generator_rotation_types = prim_rotation_generators
.iter()
.map(identify_rotation_type)
.collect::<Vec<_>>();
Expand All @@ -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::<Vec<_>>(),
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
Expand All @@ -90,53 +69,19 @@ 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;
}
}
}
}
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));
}
}
}
}
}
2 changes: 1 addition & 1 deletion moyo/src/identify/space_group.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,7 @@ fn correction_transformation_matrices(
}

/// Search for origin_shift such that (trans_mat, origin_shift) transforms <prim_operations> into <db_prim_generators>
fn match_origin_shift(
pub fn match_origin_shift(
prim_operations: &Operations,
trans_mat: &UnimodularLinear,
db_prim_generators: &Operations,
Expand Down
Loading

0 comments on commit 55e17bd

Please sign in to comment.