Merge pull request #64 from spglib/monoclinic-standardization

more offset search for wyckoff position assignment

moyo/src/math/cycle_checker.rs

@@ -1,26 +1,33 @@
 use std::collections::HashSet;
 
-use nalgebra::Matrix3;
+use nalgebra::base::allocator::Allocator;
+use nalgebra::{DefaultAllocator, Dim, OMatrix};
 
 /// Record transformation matrices during lattice reduction
 #[derive(Debug)]
-pub struct CycleChecker {
-    visited: HashSet<Matrix3<i32>>,
+pub struct CycleChecker<N: Dim>
+where
+    DefaultAllocator: Allocator<i32, N, N>,
+{
+    visited: HashSet<OMatrix<i32, N, N>>,
 }
 
-impl CycleChecker {
+impl<N: Dim> CycleChecker<N>
+where
+    DefaultAllocator: Allocator<i32, N, N>,
+{
     pub fn new() -> Self {
         Self {
             visited: HashSet::new(),
         }
     }
 
     /// If `matrix` is not visited, insert it and return true.
-    pub fn insert(&mut self, matrix: &Matrix3<i32>) -> bool {
+    pub fn insert(&mut self, matrix: &OMatrix<i32, N, N>) -> bool {
         if self.visited.contains(matrix) {
             return false;
         }
-        self.visited.insert(*matrix);
+        self.visited.insert(matrix.clone());
         true
     }
 }

moyo/src/math/minkowski.rs

@@ -1,5 +1,8 @@
 use itertools::Itertools;
-use nalgebra::{DMatrix, DVector, Matrix3, Vector3, U3};
+use nalgebra::allocator::Allocator;
+use nalgebra::{
+    DMatrix, DVector, DefaultAllocator, Dim, DimName, Matrix3, OMatrix, OVector, Vector3, U3,
+};
 
 use super::cycle_checker::CycleChecker;
 use super::elementary::swapping_column_matrix;
@@ -10,7 +13,7 @@ const EPS: f64 = 1e-8;
 pub fn minkowski_reduce(basis: &Matrix3<f64>) -> (Matrix3<f64>, Matrix3<i32>) {
     let mut reduced_basis = *basis;
     let mut trans_mat = Matrix3::<i32>::identity();
-    minkowski_reduce_greedy(&mut reduced_basis, &mut trans_mat, 3);
+    minkowski_reduce_greedy(U3, &mut reduced_basis, &mut trans_mat, 3);
 
     // Preserve parity
     if trans_mat.map(|e| e as f64).determinant() < 0. {
@@ -25,68 +28,75 @@ pub fn minkowski_reduce(basis: &Matrix3<f64>) -> (Matrix3<f64>, Matrix3<i32>) {
 /// basis * trans_mat is always preserved
 ///
 /// [1] https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1992.0004
-fn minkowski_reduce_greedy(basis: &mut Matrix3<f64>, trans_mat: &mut Matrix3<i32>, dim: usize) {
+fn minkowski_reduce_greedy<N: Dim + DimName>(
+    dim: N,
+    basis: &mut OMatrix<f64, N, N>,
+    trans_mat: &mut OMatrix<i32, N, N>,
+    rank: usize,
+) where
+    DefaultAllocator: Allocator<f64, N, N> + Allocator<i32, N, N> + Allocator<f64, N>,
+{
     // Line 1: exit condition
-    if dim == 1 {
+    if rank == 1 {
         return;
     }
 
     let mut cc = CycleChecker::new();
     loop {
         // Line 3: sort basis vectors by their lengths
         let lengths: Vec<f64> = basis.column_iter().map(|column| column.norm()).collect();
-        for i in 0..dim {
-            for j in 0..(dim - 1 - i) {
+        for i in 0..rank {
+            for j in 0..(rank - 1 - i) {
                 if lengths[j] > lengths[j + 1] + EPS {
                     basis.swap_columns(j, j + 1);
-                    *trans_mat *= swapping_column_matrix(U3, j, j + 1);
+                    *trans_mat *= swapping_column_matrix(dim, j, j + 1);
                 }
             }
         }
 
         // Line 4: Recursive call
-        minkowski_reduce_greedy(basis, trans_mat, dim - 1);
+        minkowski_reduce_greedy(dim, basis, trans_mat, rank - 1);
 
         // Line 5: Solve the closest vector problem (CVP) for basis.column(d - 1)
         // linear projection (Gram-Schmidt)
-        let h = DMatrix::from_fn(dim - 1, dim - 1, |i, j| {
+        let h = DMatrix::from_fn(rank - 1, rank - 1, |i, j| {
             basis.column(i).dot(&basis.column(j)) / basis.column(i).norm_squared()
         });
-        let u = DVector::from_fn(dim - 1, |i, _| {
-            basis.column(i).dot(&basis.column(dim - 1)) / basis.column(i).norm_squared()
+        let u = DVector::from_fn(rank - 1, |i, _| {
+            basis.column(i).dot(&basis.column(rank - 1)) / basis.column(i).norm_squared()
         });
         let gs_coeffs = h.try_inverse().unwrap() * u;
-        let gs_coeffs_rint = DVector::from_fn(dim - 1, |i, _| gs_coeffs[(i, 0)].round() as i32);
+        let gs_coeffs_rint = DVector::from_fn(rank - 1, |i, _| gs_coeffs[(i, 0)].round() as i32);
 
-        // Since basis[0..(dim-1)] are already Minkowski reduced, we only need to check Voronoi relevant vectors from y_int
+        // Since basis[0..(rank-1)] are already Minkowski reduced, we only need to check Voronoi relevant vectors from y_int
         let mut cvp_min = f64::INFINITY;
-        let mut coeffs_argmin = DVector::<i32>::zeros(dim - 1);
-        let mut c_argmin = Vector3::<f64>::zeros();
+        let mut coeffs_argmin = DVector::<i32>::zeros(rank - 1);
+        let mut c_argmin = OVector::<f64, N>::zeros();
 
         // [-1, 0, 1] is sufficient for up to three dimensional Minkowski-reduced basis
-        for voronoi_vector in (0..dim - 1).map(|_| -1..=1).multi_cartesian_product() {
-            let coeffs = gs_coeffs_rint.clone() + DVector::from_iterator(dim - 1, voronoi_vector);
-            let c = basis.columns(0, dim - 1) * coeffs.map(|e| e as f64);
-            let cvp = (c - basis.column(dim - 1)).norm();
+        for voronoi_vector in (0..rank - 1).map(|_| -1..=1).multi_cartesian_product() {
+            let coeffs = gs_coeffs_rint.clone() + DVector::from_iterator(rank - 1, voronoi_vector);
+            let c = basis.columns(0, rank - 1) * coeffs.map(|e| e as f64);
+            let cvp = (c.clone() - basis.column(rank - 1)).norm();
             if cvp < cvp_min {
                 cvp_min = cvp;
                 coeffs_argmin = coeffs.clone();
                 c_argmin = c;
             }
         }
 
-        // Line 6: update basis.column(dim - 1)
+        // Line 6: update basis.column(rank - 1)
         for j in 0..3 {
-            basis[(j, dim - 1)] -= c_argmin[j];
+            basis[(j, rank - 1)] -= c_argmin[j];
         }
-        let mut add_mat = Matrix3::<i32>::identity();
-        for i in 0..(dim - 1) {
-            add_mat[(i, dim - 1)] = -coeffs_argmin[i];
+        let mut add_mat = OMatrix::<i32, N, N>::identity();
+        for i in 0..(rank - 1) {
+            add_mat[(i, rank - 1)] = -coeffs_argmin[i];
         }
         *trans_mat *= add_mat;
 
         // Line 7: loop until length ordering is changed
-        if basis.column(dim - 1).norm() + EPS > basis.column(dim - 2).norm() {
+        if basis.column(rank - 1).norm() + EPS > basis.column(rank - 2).norm() {
             break;
         }
 

moyo/src/symmetrize/standardize.rs

@@ -294,7 +294,13 @@ fn assign_wyckoff_position(
         //    = lattice * (D * R^-1 * y - L * (offset + position - space.origin)
         let space = WyckoffPositionSpace::new(wyckoff.coordinates);
         let snf = SNF::new(&space.linear);
-        for offset in iproduct!(-1..=1, -1..=1, -1..=1) {
+
+        let iter_multi_1 = iproduct!(-1..=1, -1..=1, -1..=1);
+        let iter_multi_2 = iproduct!(-2..=2, -2..=2, -2..=2).filter(|&(n1, n2, n3)| {
+            (n1 as i32).abs() == 2 || (n2 as i32).abs() == 2 || (n3 as i32).abs() == 2
+        });
+
+        for offset in iter_multi_1.chain(iter_multi_2) {
             let offset = Vector3::new(offset.0 as f64, offset.1 as f64, offset.2 as f64);
             let b = snf.l.map(|e| e as f64) * (offset + position - space.origin);
             let mut rinvy = Vector3::zeros();

moyo/tests/assets/AB_mC8_15_e_a.json

@@ -0,0 +1 @@
+{"lattice":{"basis":[4.807835971596717,0.0,5.95538517065568,-8.07000216883193e-16,13.1793137,8.07000216883193e-16,0.0,0.0,11.74174177]},"positions":[[0.0,0.57413983,0.25],[0.0,0.42586017,0.75],[0.5,0.07413983,0.25],[0.5,0.92586017,0.75],[0.0,0.0,0.0],[0.0,0.0,0.5],[0.5,0.5,0.0],[0.5,0.5,0.5]],"numbers":[20,20,20,20,26,26,26,26]}

moyo/tests/test_moyo_dataset.rs

@@ -612,17 +612,30 @@ fn test_with_high_symprec_and_angle_tolerance() {
 
 #[test]
 fn test_wyckoff_position_assignment() {
-    // https://github.com/spglib/moyo/issues/54
-    let path = Path::new("tests/assets/wyckoff_edge_case.json");
-    let cell: Cell = serde_json::from_str(&fs::read_to_string(&path).unwrap()).unwrap();
-
     let symprec = 1e-4;
     let angle_tolerance = AngleTolerance::Default;
     let setting = Setting::Standard;
 
-    let dataset = MoyoDataset::new(&cell, symprec, angle_tolerance, setting).unwrap();
-    assert_eq!(
-        dataset.wyckoffs,
-        vec!['e', 'e', 'e', 'e', 'a', 'a', 'a', 'a']
-    );
+    {
+        // https://github.com/CompRhys/aviary/pull/96#issuecomment-2628703353
+        let path = Path::new("tests/assets/AB_mC8_15_e_a.json");
+        let cell: Cell = serde_json::from_str(&fs::read_to_string(&path).unwrap()).unwrap();
+        let dataset = MoyoDataset::new(&cell, symprec, angle_tolerance, setting).unwrap();
+        assert_eq!(
+            dataset.wyckoffs,
+            vec!['e', 'e', 'e', 'e', 'a', 'a', 'a', 'a']
+        );
+    }
+
+    {
+        // https://github.com/spglib/moyo/issues/54
+        let path = Path::new("tests/assets/wyckoff_edge_case.json");
+        let cell: Cell = serde_json::from_str(&fs::read_to_string(&path).unwrap()).unwrap();
+
+        let dataset = MoyoDataset::new(&cell, symprec, angle_tolerance, setting).unwrap();
+        assert_eq!(
+            dataset.wyckoffs,
+            vec!['e', 'e', 'e', 'e', 'a', 'a', 'a', 'a']
+        );
+    }
 }