Search [-1, 1] transformation matrix first

bench/mp/mp.py

@@ -2,18 +2,26 @@
 
 from time import perf_counter
 
+import click
 import pandas as pd
-import spglib
 from matbench_discovery.data import DataFiles
 from pymatgen.entries.computed_entries import ComputedStructureEntry
 from pymatviz.enums import Key
 from spglib import get_symmetry_dataset
 from tqdm.auto import tqdm
 
 import moyopy
+from moyopy.interface import MoyoAdapter
 
+SYMPREC_LIST = [1e-4, 3e-4, 1e-3, 3e-3, 1e-2, 3e-2, 1e-1]
 
-def main():
+
+@click.group()
+def cli(): ...
+
+
+@cli.command()
+def perf_spglib():
     data_path = DataFiles.mp_computed_structure_entries.path
     df = pd.read_json(data_path).set_index(Key.mat_id)
 
@@ -25,31 +33,58 @@ def main():
             basis = structure.lattice.matrix
             positions = structure.frac_coords
             numbers = [site.specie.Z for site in structure]
-
-            moyopy_cell = moyopy.Cell(basis.tolist(), positions.tolist(), numbers)
             spglib_cell = (basis, positions, numbers)
 
-            for symprec in [1e-4, 3e-4, 1e-3, 3e-3, 1e-2, 3e-2, 1e-1]:
+            for i, symprec in enumerate(SYMPREC_LIST):
                 try:
-                    start = perf_counter()
-                    moyopy_dataset = moyopy.MoyoDataset(moyopy_cell, symprec=symprec)
-                    time_moyopy = perf_counter() - start
-
                     start = perf_counter()
                     spglib_dataset = get_symmetry_dataset(spglib_cell, symprec=symprec)
                     time_spglib = perf_counter() - start
 
                     all_stats.append(
                         {
-                            "time_moyopy": time_moyopy,
-                            "time_spglib": time_spglib,
+                            "id": f"{material_id}_{i}",
+                            "material_id": material_id,
                             "symprec": symprec,
+                            "time_spglib": time_spglib,
                             "num_atoms": len(numbers),
+                            "number_spglib": spglib_dataset["number"],
+                        }
+                    )
+                except:  # noqa: E722
+                    print(f"Abort: {material_id=} {symprec=}")
+
+    df_stats = pd.DataFrame(all_stats)
+    df_stats.to_json("stats_spglib.json")
+
+
+@cli.command()
+def perf_moyopy():
+    data_path = DataFiles.mp_computed_structure_entries.path
+    df = pd.read_json(data_path).set_index(Key.mat_id)
+
+    all_stats = []
+
+    with tqdm(df.iterrows(), total=len(df)) as pbar:
+        for material_id, row in pbar:
+            structure = ComputedStructureEntry.from_dict(row["entry"]).structure
+            numbers = [site.specie.Z for site in structure]
+            moyopy_cell = MoyoAdapter.from_structure(structure)
+
+            for i, symprec in enumerate(SYMPREC_LIST):
+                try:
+                    start = perf_counter()
+                    moyopy_dataset = moyopy.MoyoDataset(moyopy_cell, symprec=symprec)
+                    time_moyopy = perf_counter() - start
+
+                    all_stats.append(
+                        {
+                            "id": f"{material_id}_{i}",
                             "material_id": material_id,
+                            "symprec": symprec,
+                            "time_moyopy": time_moyopy,
+                            "num_atoms": len(numbers),
                             "number_moyopy": moyopy_dataset.number,
-                            "number_spglib": spglib_dataset["number"],
-                            "moyopy": moyopy.__version__,
-                            "spglib": spglib.__version__,
                         }
                     )
                 except:  # noqa: E722
@@ -59,14 +94,9 @@ def main():
                         f.write(moyopy_cell.serialize_json())
                     return
 
-                # if moyopy_dataset.number != spglib_dataset["number"]:
-                #     print(
-                #         f"Inconsistent: {material_id=} {moyopy_dataset.number=} {spglib_dataset['number']=}"  # noqa: E501
-                #     )
-
     df_stats = pd.DataFrame(all_stats)
-    df_stats.to_json("stats.json")
+    df_stats.to_json("stats_moyopy.json")
 
 
 if __name__ == "__main__":
-    main()
+    cli()

moyo/src/base/lattice.rs

@@ -91,6 +91,18 @@ impl Lattice {
         self.basis.determinant().abs()
     }
 
+    #[allow(dead_code)]
+    pub(crate) fn lattice_constant(&self) -> [f64; 6] {
+        let g = self.metric_tensor();
+        let a = g[(0, 0)].sqrt();
+        let b = g[(1, 1)].sqrt();
+        let c = g[(2, 2)].sqrt();
+        let alpha = (g[(1, 2)] / (b * c)).acos().to_degrees();
+        let beta = (g[(0, 2)] / (a * c)).acos().to_degrees();
+        let gamma = (g[(0, 1)] / (a * b)).acos().to_degrees();
+        [a, b, c, alpha, beta, gamma]
+    }
+
     /// Rotate the lattice by the given rotation matrix
     pub fn rotate(&self, rotation_matrix: &Matrix3<f64>) -> Self {
         Self {

moyo/src/identify/point_group.rs

@@ -258,22 +258,28 @@ pub fn iter_trans_mat_basis(
 pub fn iter_unimodular_trans_mat(
     trans_mat_basis: Vec<Matrix3<i32>>,
 ) -> impl Iterator<Item = UnimodularLinear> {
-    (0..trans_mat_basis.len())
+    // First try with coefficients in [-1, 1]
+    let iter_multi_1 = (0..trans_mat_basis.len())
+        .map(|_| -1..=1)
+        .multi_cartesian_product();
+    let iter_multi_2 = (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
-            }
-        })
+        .filter(|comb| comb.iter().any(|&e| (e as i32).abs() == 2));
+
+    iter_multi_1.chain(iter_multi_2).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)]

moyo/tests/assets/wyckoff_edge_case.json

@@ -0,0 +1 @@
+{"lattice":{"basis":[6.57275068253895,0.0,1.1376274714501016,-6.610122181503525e-16,10.79514875,6.610122181503525e-16,0.0,0.0,9.91142379]},"positions":[[0.0,0.59536646,0.25],[0.0,0.40463354,0.75],[0.5,0.09536646,0.25],[0.5,0.90463354,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

@@ -609,3 +609,20 @@ fn test_with_high_symprec_and_angle_tolerance() {
 
     let _ = MoyoDataset::new(&cell, symprec, angle_tolerance, setting).unwrap();
 }
+
+#[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']
+    );
+}