diff --git a/Cargo.lock b/Cargo.lock index ebf09d6..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", @@ -626,6 +670,7 @@ dependencies = [ "kiddo", "log", "nalgebra", + "once_cell", "rand", "rstest", "serde", diff --git a/moyo/Cargo.toml b/moyo/Cargo.toml index b42233c..6b385e5 100644 --- a/moyo/Cargo.toml +++ b/moyo/Cargo.toml @@ -24,7 +24,8 @@ 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] rand = "0.8" diff --git a/moyo/src/base/error.rs b/moyo/src/base/error.rs index 880ea11..a4621b8 100644 --- a/moyo/src/base/error.rs +++ b/moyo/src/base/error.rs @@ -27,6 +27,10 @@ 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")] StandardizationError, #[error("Wyckoff position assignment failed")] diff --git a/moyo/src/base/transformation.rs b/moyo/src/base/transformation.rs index 89c2f99..ca1eff2 100644 --- a/moyo/src/base/transformation.rs +++ b/moyo/src/base/transformation.rs @@ -1,10 +1,12 @@ +use std::ops::Mul; + use itertools::iproduct; 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 +75,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) @@ -98,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 { @@ -188,6 +219,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..61400c0 100644 --- a/moyo/src/data.rs +++ b/moyo/src/data.rs @@ -11,13 +11,18 @@ 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_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, +}; 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..338e870 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, @@ -102,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) @@ -206,6 +214,25 @@ 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) + } else { + None + } + } + + pub fn primitive_generators(&self) -> MagneticOperations { + Transformation::from_linear(self.centering.linear()) + .inverse_transform_magnetic_operations(&self.generators) + } } /// Tokenize string by whitespaces @@ -562,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::*; @@ -587,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/data/magnetic_hall_symbol_database.rs b/moyo/src/data/magnetic_hall_symbol_database.rs index 32b8dd3..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; +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,15 +17,20 @@ impl MagneticHallSymbolEntry { uni_number, } } -} -// TODO: -// smallest and largest UNI numbers for each Hall number -// let hall_number_to_uni_numbers: HashMap; + pub fn construct_type(&self) -> ConstructType { + get_magnetic_space_group_type(self.uni_number) + .unwrap() + .construct_type + } -// TODO: -// smallest and largest Hall numbers for XSG of each UNI number -// let uni_number_to_xsg_hall_numbers: HashMap; + 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 { MAGNETIC_HALL_SYMBOL_DATABASE @@ -29,7 +38,8 @@ pub fn magnetic_hall_symbol_entry(uni_number: UNINumber) -> Option impl Iterator { MAGNETIC_HALL_SYMBOL_DATABASE.iter() diff --git a/moyo/src/data/magnetic_space_group.rs b/moyo/src/data/magnetic_space_group.rs index 1ecb0d9..b241152 100644 --- a/moyo/src/data/magnetic_space_group.rs +++ b/moyo/src/data/magnetic_space_group.rs @@ -1,9 +1,15 @@ +use std::ops::RangeInclusive; + +use once_cell::sync::Lazy; + use super::hall_symbol_database::Number; +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/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.rs b/moyo/src/identify.rs index 1e5ca66..af9ffd0 100644 --- a/moyo/src/identify.rs +++ b/moyo/src/identify.rs @@ -1,4 +1,8 @@ +mod magnetic_space_group; +mod normalizer; 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 new file mode 100644 index 0000000..c658b8f --- /dev/null +++ b/moyo/src/identify/magnetic_space_group.rs @@ -0,0 +1,446 @@ +use std::collections::HashMap; + +use itertools::Itertools; +use log::debug; + +use super::normalizer::integral_normalizer; +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::{ + 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, + ConstructType, HallSymbol, MagneticHallSymbol, MagneticHallSymbolEntry, 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 (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 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 { + if get_magnetic_space_group_type(uni_number) + .unwrap() + .construct_type + != construct_type + { + 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, + }); + } + + 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_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); + + 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( + prim_mag_operations: &MagneticOperations, + epsilon: f64, +) -> Option<(Operations, ConstructType)> { + 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) + { + debug!("Input magnetic operations are incomplete."); + return None; + } + + 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) + { + // 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)) +} + +/// XSG: take only operations without time-reversal +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; + } + xsg.push(mops.operation.clone()); + } + xsg +} + +/// 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, + 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, is_type2) +} + +/// 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(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(); + + 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::>(); + 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, + ) { + 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 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; + use test_log::test as test_with_log; + + use super::*; + 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(); + mhs.primitive_traverse() + } + + #[rstest] + #[case(2, ConstructType::Type2, 2, 1, 1)] + #[case(1594, ConstructType::Type1, 48, 48, 48)] + #[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( + #[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 = primitive_maximal_space_subgroup_from_magnetic_space_group(&prim_mag_operations); + assert_eq!(xsg.len(), order_xsg); + + 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, epsilon).unwrap(); + 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(&entry); + + let mhs = MagneticHallSymbol::new(&entry.magnetic_hall_symbol).unwrap(); + let identity = Rotation::identity(); + let mut 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(), + }; + 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() { + 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 { + // 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 new file mode 100644 index 0000000..e9842b6 --- /dev/null +++ b/moyo/src/identify/normalizer.rs @@ -0,0 +1,41 @@ +use itertools::Itertools; + +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}; + +/// Return integral normalizer of the point group representative up to its centralizer. + +/// Integral normalizers for all arithmetic point groups up to their centralizers. + +/// 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_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 mut conjugators = vec![]; + for trans_mat_basis in + iter_trans_mat_basis(prim_rotations, rotation_types, prim_rotation_generators) + { + 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; + } + } + } + conjugators +} diff --git a/moyo/src/identify/point_group.rs b/moyo/src/identify/point_group.rs index d8ec6ea..96949dc 100644 --- a/moyo/src/identify/point_group.rs +++ b/moyo/src/identify/point_group.rs @@ -2,14 +2,15 @@ use std::cmp::Ordering; use itertools::Itertools; use log::debug; -use nalgebra::{Dyn, Matrix3, OMatrix, OVector, U9}; +use nalgebra::Matrix3; -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, }; -use crate::math::IntegerLinearSystem; +use crate::math::sylvester3; /// Crystallographic point group with group-type information #[derive(Debug)] @@ -51,20 +52,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 -} - /// Faster matching algorithm for cubic point groups fn match_with_cubic_point_group( prim_rotations: &Rotations, @@ -85,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 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) = sylvester( - &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 => {} + 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); + + // 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, + }); } } @@ -162,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; @@ -173,178 +131,21 @@ 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) = sylvester( - &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, - }); - } - } - } - } - } - - 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, + 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() - { - 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)]; + 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, + }); } } } - 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; - - 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"), - } + Err(MoyoError::ArithmeticCrystalClassIdentificationError) } /// Use look up table in Table 6 of https://arxiv.org/pdf/1808.01590.pdf @@ -418,11 +219,68 @@ 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, + ) + }) +} + +/// 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; - use super::{integral_normalizer, PointGroup, PointGroupRepresentative}; + use super::*; use crate::base::traverse; #[test] @@ -464,37 +322,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"), + } +} diff --git a/moyo/src/identify/space_group.rs b/moyo/src/identify/space_group.rs index 0ad35e2..7ff20ae 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, @@ -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, @@ -249,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}; @@ -273,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(); @@ -310,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(); @@ -334,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() { 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/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}; 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}; 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, - }) } }