qsym2/symmetry/symmetry_core/
mod.rs

1//! Molecular symmetry elements and their detections.
2
3use std::collections::hash_map::Entry::Vacant;
4use std::collections::{HashMap, HashSet};
5use std::error::Error;
6use std::fmt;
7
8use anyhow::{self, ensure, format_err};
9use derive_builder::Builder;
10use indexmap::IndexSet;
11use itertools::Itertools;
12use log;
13use nalgebra::{Point3, Vector3};
14use rayon::prelude::*;
15use serde::{Deserialize, Serialize};
16
17use crate::auxiliary::atom::Atom;
18use crate::auxiliary::geometry::{self, Transform};
19use crate::auxiliary::molecule::Molecule;
20use crate::rotsym::{self, RotationalSymmetry};
21use crate::symmetry::symmetry_element::symmetry_operation::{
22    sort_operations, SpecialSymmetryTransformation, SymmetryOperation,
23};
24use crate::symmetry::symmetry_element::{
25    AntiunitaryKind, SymmetryElement, SymmetryElementKind, ROT, SIG, SO3, TRROT, TRSIG,
26};
27use crate::symmetry::symmetry_element_order::{ElementOrder, ORDER_1, ORDER_2};
28use crate::symmetry::symmetry_symbols::deduce_sigma_symbol;
29
30#[cfg(test)]
31mod symmetry_core_tests;
32
33#[cfg(test)]
34mod symmetry_group_detection_tests;
35
36#[derive(Debug)]
37pub struct PointGroupDetectionError(pub String);
38
39impl fmt::Display for PointGroupDetectionError {
40    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
41        write!(f, "Point-group detection error: {}.", self.0)
42    }
43}
44
45impl Error for PointGroupDetectionError {}
46
47/// Structure for storing and managing information required for symmetry analysis.
48#[derive(Clone, Builder, Debug, Serialize, Deserialize)]
49pub struct PreSymmetry {
50    /// The original molecule.
51    #[builder(setter(custom))]
52    pub(crate) original_molecule: Molecule,
53
54    /// The recentred molecule to be symmetry-analysed.
55    #[builder(setter(custom))]
56    pub(crate) recentred_molecule: Molecule,
57
58    /// Threshold for relative comparisons of moments of inertia.
59    #[builder(setter(custom))]
60    pub(crate) moi_threshold: f64,
61
62    /// The rotational symmetry of [`Self::recentred_molecule`] based on its moments of
63    /// inertia.
64    #[builder(setter(skip), default = "self.calc_rotational_symmetry()")]
65    pub(crate) rotational_symmetry: RotationalSymmetry,
66
67    /// The groups of symmetry-equivalent atoms in [`Self::recentred_molecule`].
68    #[builder(setter(skip), default = "self.calc_sea_groups()")]
69    pub(crate) sea_groups: Vec<Vec<Atom>>,
70
71    /// Threshold for relative distance comparisons.
72    #[builder(setter(skip), default = "self.get_dist_threshold()")]
73    pub(crate) dist_threshold: f64,
74}
75
76impl PreSymmetryBuilder {
77    /// Initialises the molecule to be symmetry-analysed.
78    ///
79    /// The original molecule and a recentred copy of the molecule will be stored. The recentred
80    /// copy will be used for subsequent symmetry analyses.
81    ///
82    /// # Arguments
83    ///
84    /// * molecule - The molecule to be symmetry-analysed.
85    ///
86    /// # Returns
87    /// A mutable borrow of `[Self]`.
88    pub fn molecule(&mut self, molecule: &Molecule) -> &mut Self {
89        self.recentred_molecule = Some(molecule.recentre());
90        self.original_molecule = Some(molecule.clone());
91        self
92    }
93
94    /// Initialises the threshold for moment-of-inertia comparisons.
95    ///
96    /// # Arguments
97    ///
98    /// * thresh - The threshold for moment-of-inertia comparisons.
99    ///
100    /// # Returns
101    /// A mutable borrow of `[Self]`.
102    pub fn moi_threshold(&mut self, thresh: f64) -> &mut Self {
103        if thresh >= f64::EPSILON {
104            self.moi_threshold = Some(thresh);
105        } else {
106            log::error!(
107                "Threshold value {} is invalid. Threshold must be at least the machine epsilon.",
108                thresh
109            );
110            self.moi_threshold = None;
111        }
112        self
113    }
114
115    /// Calculates the rotational symmetry of [`Self::recentred_molecule`].
116    fn calc_rotational_symmetry(&self) -> RotationalSymmetry {
117        let com = self
118            .recentred_molecule
119            .as_ref()
120            .expect("A molecule has not been set.")
121            .calc_com();
122        let inertia = self
123            .recentred_molecule
124            .as_ref()
125            .expect("A molecule has not been set.")
126            .calc_inertia_tensor(&com);
127        approx::assert_relative_eq!(
128            com,
129            Point3::origin(),
130            epsilon = self
131                .recentred_molecule
132                .as_ref()
133                .expect("A molecule has not been set.")
134                .threshold,
135            max_relative = self
136                .recentred_molecule
137                .as_ref()
138                .expect("A molecule has not been set.")
139                .threshold
140        );
141        rotsym::calc_rotational_symmetry(
142            &inertia,
143            self.moi_threshold.expect("MoI threshold has not been set."),
144        )
145    }
146
147    /// Calculates the symmetry-equivalent groups of the atoms in [`Self::recentred_molecule`].
148    fn calc_sea_groups(&self) -> Vec<Vec<Atom>> {
149        self.recentred_molecule
150            .as_ref()
151            .expect("A molecule has not been set.")
152            .calc_sea_groups()
153    }
154
155    /// Returns the distance threshold of [`Self::recentred_molecule`].
156    fn get_dist_threshold(&self) -> f64 {
157        self.recentred_molecule
158            .as_ref()
159            .expect("A molecule has not been set.")
160            .threshold
161    }
162}
163
164impl PreSymmetry {
165    /// Returns a builder to construct a new pre-symmetry struct.
166    ///
167    /// # Returns
168    ///
169    /// A builder to construct a new pre-symmetry struct.
170    #[must_use]
171    pub(crate) fn builder() -> PreSymmetryBuilder {
172        PreSymmetryBuilder::default()
173    }
174
175    /// Checks for the existence of the proper symmetry element $`C_n`$  or $`\theta C_n`$ along
176    /// `axis` in [`Self::recentred_molecule`].
177    ///
178    /// Non-time-reversed elements are always preferred.
179    ///
180    /// # Arguments
181    ///
182    /// * `order` - The geometrical order $`n`$ of the rotation axis. Only finite
183    /// orders are supported.
184    /// * `axis` - The rotation axis.
185    /// * `tr` - A flag indicating if time reversal should also be considered in case the
186    /// non-time-reversed symmetry element does not exist.
187    ///
188    /// # Returns
189    ///
190    /// An [`Option`] containing the proper kind if the $`C_n`$ or $`\theta C_n`$ element exists in
191    /// [`Self::recentred_molecule`]. If not, [`None`] is returned.
192    #[allow(clippy::trivially_copy_pass_by_ref)]
193    fn check_proper(
194        &self,
195        order: &ElementOrder,
196        axis: &Vector3<f64>,
197        tr: bool,
198    ) -> Option<SymmetryElementKind> {
199        assert_ne!(
200            *order,
201            ElementOrder::Inf,
202            "This method does not work for infinite-order elements."
203        );
204        let angle = 2.0 * std::f64::consts::PI / order.to_float();
205        let rotated_mol = self.recentred_molecule.rotate(angle, axis);
206        if rotated_mol == self.recentred_molecule {
207            Some(ROT)
208        } else if tr {
209            let tr_rotated_mol = rotated_mol.reverse_time();
210            if tr_rotated_mol == self.recentred_molecule {
211                Some(TRROT)
212            } else {
213                None
214            }
215        } else {
216            None
217        }
218    }
219
220    /// Checks for the existence of the improper symmetry element $`S_n`$, $`\dot{S}_n`$,
221    /// $`\theta S_n`$, or $`\theta \dot{S}_n`$ along `axis` in [`Self::recentred_molecule`].
222    ///
223    /// Non-time-reversed elements are always preferred.
224    ///
225    /// # Arguments
226    ///
227    /// * `order` - The geometrical order $`n`$ of the improper rotation axis. Only
228    /// finite orders are supported.
229    /// * `axis` - The rotation axis.
230    /// * `kind` - The convention in which the improper element is defined. The time reversal
231    /// property of this does not matter.
232    /// * `tr` - A flag indicating if time reversal should also be considered in case the
233    /// non-time-reversed symmetry element does not exist.
234    ///
235    /// # Returns
236    ///
237    /// An [`Option`] containing the improper kind if the $`S_n`$, $`\theta S_n`$, $`\theta S_n`$,
238    /// or $`\theta \dot{S}_n`$ element exists in [`Self::recentred_molecule`]. If not, [`None`] is
239    /// returned.
240    #[allow(clippy::trivially_copy_pass_by_ref)]
241    fn check_improper(
242        &self,
243        order: &ElementOrder,
244        axis: &Vector3<f64>,
245        kind: &SymmetryElementKind,
246        tr: bool,
247    ) -> Option<SymmetryElementKind> {
248        assert_ne!(
249            *order,
250            ElementOrder::Inf,
251            "This method does not work for infinite-order elements."
252        );
253        let angle = 2.0 * std::f64::consts::PI / order.to_float();
254        let transformed_mol = self.recentred_molecule.improper_rotate(
255            angle,
256            axis,
257            &kind.to_tr(false).try_into().unwrap_or_else(|err| {
258                log::error!("Error detected: {err}.");
259                panic!("Error detected: {err}.")
260            }),
261        );
262        if transformed_mol == self.recentred_molecule {
263            Some(kind.to_tr(false))
264        } else if tr {
265            let tr_transformed_mol = transformed_mol.reverse_time();
266            if tr_transformed_mol == self.recentred_molecule {
267                Some(kind.to_tr(true))
268            } else {
269                None
270            }
271        } else {
272            None
273        }
274    }
275}
276
277/// Structure for storing and managing symmetry analysis results.
278#[derive(Builder, Clone, Debug, Serialize, Deserialize)]
279pub struct Symmetry {
280    /// The determined point group in Schönflies notation.
281    #[builder(setter(skip, strip_option), default = "None")]
282    pub group_name: Option<String>,
283
284    /// The symmetry elements found.
285    ///
286    /// Each entry in the hash map is for one kind of symmetry elements: the key gives the kind,
287    /// and the value is a hash map where each key gives the order and the corresponding value
288    /// gives the [`HashSet`] of the elements with that order.
289    ///
290    /// Note that for improper elements, the mirror-plane convention is preferred.
291    #[builder(setter(skip), default = "HashMap::new()")]
292    pub(crate) elements:
293        HashMap<SymmetryElementKind, HashMap<ElementOrder, IndexSet<SymmetryElement>>>,
294
295    /// The symmetry generators found.
296    ///
297    /// Each entry in the hash map is for one kind of symmetry generators: the key gives the kind,
298    /// and the value is a hash map where each key gives the order and the corresponding value
299    /// gives the [`HashSet`] of the generators with that order.
300    ///
301    /// Note that for improper generatrors, the mirror-plane convention is preferred.
302    #[builder(setter(skip), default = "HashMap::new()")]
303    pub(crate) generators:
304        HashMap<SymmetryElementKind, HashMap<ElementOrder, IndexSet<SymmetryElement>>>,
305}
306
307impl Symmetry {
308    /// Returns a builder to construct a new symmetry struct.
309    ///
310    /// # Returns
311    ///
312    /// A builder to construct a new symmetry struct.
313    #[must_use]
314    fn builder() -> SymmetryBuilder {
315        SymmetryBuilder::default()
316    }
317
318    /// Construct a new and empty symmetry struct.
319    #[must_use]
320    pub fn new() -> Self {
321        Symmetry::builder()
322            .build()
323            .expect("Unable to construct a `Symmetry` structure.")
324    }
325
326    /// Performs point-group detection analysis and populates the fields in this structure with the
327    /// results.
328    ///
329    /// # Arguments
330    ///
331    /// * `presym` - A pre-symmetry-analysis structure containing the molecule and its rotational
332    /// symmetry required for point-group detection.
333    /// * `tr` - A boolean indicating if time reversal should also be considered. A time-reversed
334    /// symmetry element will only be considered if its non-time-reversed version turns out to be
335    /// not a symmetry element.
336    pub fn analyse(&mut self, presym: &PreSymmetry, tr: bool) -> Result<&mut Self, anyhow::Error> {
337        log::debug!("Rotational symmetry found: {}", presym.rotational_symmetry);
338
339        if tr {
340            log::debug!("Antiunitary symmetry generated by time reversal will be considered.");
341        };
342
343        // Add the identity, which must always exist.
344        let c1 = SymmetryElement::builder()
345            .threshold(presym.dist_threshold)
346            .proper_order(ORDER_1)
347            .proper_power(1)
348            .raw_axis(Vector3::new(0.0, 0.0, 1.0))
349            .kind(ROT)
350            .rotation_group(SO3)
351            .build()
352            .map_err(|err| format_err!(err))?;
353        self.add_proper(ORDER_1, c1.raw_axis(), false, presym.dist_threshold, false);
354
355        // Identify all symmetry elements and generators
356        match &presym.rotational_symmetry {
357            RotationalSymmetry::Spherical => self.analyse_spherical(presym, tr)?,
358            RotationalSymmetry::ProlateLinear => self.analyse_linear(presym, tr)?,
359            RotationalSymmetry::OblatePlanar
360            | RotationalSymmetry::OblateNonPlanar
361            | RotationalSymmetry::ProlateNonLinear => self.analyse_symmetric(presym, tr)?,
362            RotationalSymmetry::AsymmetricPlanar | RotationalSymmetry::AsymmetricNonPlanar => {
363                self.analyse_asymmetric(presym, tr)?
364            }
365        }
366
367        if tr {
368            if self.get_elements(&TRROT).is_none()
369                && self.get_elements(&TRSIG).is_none()
370                && self.get_generators(&TRROT).is_none()
371                && self.get_generators(&TRSIG).is_none()
372            {
373                log::debug!("Antiunitary symmetry requested, but so far only non-time-reversed elements found.");
374                // Time-reversal requested, but the above analysis gives only non-time-reversed
375                // elements, which means the system must also contain time reversal as a symmetry
376                // operation. This implies that the group is a grey group.
377                if presym.recentred_molecule == presym.recentred_molecule.reverse_time() {
378                    log::debug!("Time reversal is a symmetry element. This is a grey group.");
379                    // Add time-reversed copies of proper elements
380                    self.elements.insert(
381                        TRROT,
382                        self.get_elements(&ROT)
383                            .ok_or_else(|| format_err!("No proper elements found."))?
384                            .iter()
385                            .map(|(order, proper_elements)| {
386                                let tr_proper_elements = proper_elements
387                                    .iter()
388                                    .map(|proper_element| {
389                                        proper_element.to_tr(true)
390                                        // let mut tr_proper_element = proper_element.clone();
391                                        // tr_proper_element.kind = proper_element.kind.to_tr(true);
392                                        // tr_proper_element
393                                    })
394                                    .collect::<IndexSet<_>>();
395                                (*order, tr_proper_elements)
396                            })
397                            .collect::<HashMap<_, _>>(),
398                    );
399                    log::debug!("Time-reversed copies of all proper elements added.");
400
401                    // Add the time-reversal element as a generator
402                    self.add_proper(ORDER_1, c1.raw_axis(), true, presym.dist_threshold, true);
403
404                    // Add time-reversed copies of improper elements, if any
405                    if self.get_elements(&SIG).is_some() {
406                        self.elements.insert(
407                            TRSIG,
408                            self.get_elements(&SIG)
409                                .ok_or_else(|| format_err!("No improper elements found."))?
410                                .iter()
411                                .map(|(order, improper_elements)| {
412                                    let tr_improper_elements = improper_elements
413                                        .iter()
414                                        .map(|improper_element| improper_element.to_tr(true))
415                                        .collect::<IndexSet<_>>();
416                                    (*order, tr_improper_elements)
417                                })
418                                .collect::<HashMap<_, _>>(),
419                        );
420                        log::debug!("Time-reversed copies of all improper elements added.");
421                    }
422
423                    // Rename the group to include the antiunitary coset generated by time reversal
424                    let unitary_group = self
425                        .group_name
426                        .as_ref()
427                        .expect("No point groups found.")
428                        .clone();
429                    self.set_group_name(format!("{unitary_group} + θ·{unitary_group}"));
430                } else {
431                    log::debug!(
432                        "Time reversal is not a symmetry element. This is an ordinary group."
433                    );
434                }
435            } else {
436                log::debug!(
437                    "Antiunitary symmetry requested and some non-time-reversed elements found."
438                );
439                log::debug!(
440                    "Time reversal is not a symmetry element. This is a black-and-white group."
441                );
442            }
443        }
444        Ok(self)
445    }
446
447    /// Adds a proper symmetry element to this struct.
448    ///
449    /// # Arguments
450    ///
451    /// * `order` - The order of the proper symmetry element.
452    /// * `axis` - The axis of rotation of the proper symmetry element.
453    /// * `generator` - A boolean indicating if this element should be added as a generator.
454    /// * `threshold` - A threshold for element comparisons.
455    /// * `tr` - A boolean indicating if time reversal is to be considered.
456    ///
457    /// # Returns
458    ///
459    /// `true` if the specified element is not present and has just been added,
460    /// `false` otherwise.
461    fn add_proper(
462        &mut self,
463        order: ElementOrder,
464        axis: &Vector3<f64>,
465        generator: bool,
466        threshold: f64,
467        tr: bool,
468    ) -> bool {
469        let positive_axis = geometry::get_standard_positive_pole(axis, threshold).normalize();
470        let proper_kind = if tr { TRROT } else { ROT };
471        let element = SymmetryElement::builder()
472            .threshold(threshold)
473            .proper_order(order)
474            .proper_power(1)
475            .raw_axis(positive_axis)
476            .kind(proper_kind)
477            .rotation_group(SO3)
478            .generator(generator)
479            .build()
480            .expect("Unable to construct a proper element.");
481        let simplified_symbol = element.get_simplified_symbol();
482        let full_symbol = element.get_full_symbol();
483        let result = if generator {
484            if let Vacant(proper_generators) = self.generators.entry(proper_kind) {
485                proper_generators.insert(HashMap::from([(order, IndexSet::from([element]))]));
486                true
487            } else {
488                let proper_generators =
489                    self.generators.get_mut(&proper_kind).unwrap_or_else(|| {
490                        panic!(
491                            "{} generators not found.",
492                            if tr { "Time-reversed proper" } else { "Proper" }
493                        )
494                    });
495
496                if let Vacant(proper_generators_order) = proper_generators.entry(order) {
497                    proper_generators_order.insert(IndexSet::from([element]));
498                    true
499                } else {
500                    proper_generators
501                        .get_mut(&order)
502                        .unwrap_or_else(|| {
503                            panic!(
504                                "Proper generators {}C{order} not found.",
505                                if tr { "θ" } else { "" }
506                            )
507                        })
508                        .insert(element)
509                }
510            }
511        } else if let Vacant(proper_elements) = self.elements.entry(proper_kind) {
512            proper_elements.insert(HashMap::from([(order, IndexSet::from([element]))]));
513            true
514        } else {
515            let proper_elements = self.elements.get_mut(&proper_kind).unwrap_or_else(|| {
516                panic!(
517                    "{} elements not found.",
518                    if tr { "Time-reversed proper" } else { "Proper" }
519                )
520            });
521
522            if let Vacant(e) = proper_elements.entry(order) {
523                e.insert(IndexSet::from([element]));
524                true
525            } else {
526                proper_elements
527                    .get_mut(&order)
528                    .unwrap_or_else(|| {
529                        panic!(
530                            "Proper elements {}C{order} not found.",
531                            if tr { "θ" } else { "" }
532                        )
533                    })
534                    .insert(element)
535            }
536        };
537
538        let dest_str = if generator {
539            "generator".to_owned()
540        } else {
541            "element".to_owned()
542        };
543        if result {
544            log::debug!(
545                "Proper rotation {} ({}): {} axis along ({:+.3}, {:+.3}, {:+.3}) added.",
546                dest_str,
547                simplified_symbol,
548                full_symbol,
549                positive_axis[0] + 0.0,
550                positive_axis[1] + 0.0,
551                positive_axis[2] + 0.0,
552            );
553        }
554        result
555    }
556
557    /// Adds an improper symmetry element to this struct.
558    ///
559    /// The improper symmetry element can be defined in whichever convention
560    /// specified by `kind`, but it will always be added in the mirror-plane
561    /// convention.
562    ///
563    /// # Arguments
564    ///
565    /// * `order` - The order of the improper symmetry element in the convention
566    ///     specified by `kind`.
567    /// * `axis` - The axis of the improper symmetry element.
568    /// * `generator` - A flag indicating if this element should be added as a generator.
569    /// * `kind` - The convention in which the improper symmetry element is defined.
570    /// * `sigma` - An optional additional string indicating the type of mirror
571    ///     plane in the case the improper element is a mirror plane.
572    /// * `threshold` - A threshold for element comparisons.
573    /// * `tr` - A boolean indicating if time reversal is to be considered.
574    ///
575    /// # Returns
576    ///
577    /// `true` if the specified element is not present and has just been added,
578    /// `false` otherwise.
579    #[allow(clippy::too_many_lines)]
580    fn add_improper(
581        &mut self,
582        order: ElementOrder,
583        axis: &Vector3<f64>,
584        generator: bool,
585        kind: SymmetryElementKind,
586        sigma: Option<String>,
587        threshold: f64,
588        tr: bool,
589    ) -> bool {
590        let positive_axis = geometry::get_standard_positive_pole(axis, threshold).normalize();
591        let mirror_kind = if tr { TRSIG } else { SIG };
592        let element = if let Some(sigma_str) = sigma {
593            assert!(sigma_str == "d" || sigma_str == "v" || sigma_str == "h");
594            let mut sym_ele = SymmetryElement::builder()
595                .threshold(threshold)
596                .proper_order(order)
597                .proper_power(1)
598                .raw_axis(positive_axis)
599                .kind(kind.to_tr(tr))
600                .rotation_group(SO3)
601                .generator(generator)
602                .build()
603                .expect("Unable to construct an improper symmetry element.")
604                .convert_to_improper_kind(&mirror_kind, false);
605            if *sym_ele.raw_proper_order() == ElementOrder::Int(1) {
606                sym_ele.additional_subscript = sigma_str;
607            }
608            sym_ele
609        } else {
610            SymmetryElement::builder()
611                .threshold(threshold)
612                .proper_order(order)
613                .proper_power(1)
614                .raw_axis(positive_axis)
615                .kind(kind.to_tr(tr))
616                .rotation_group(SO3)
617                .generator(generator)
618                .build()
619                .expect("Unable to construct an improper symmetry element.")
620                .convert_to_improper_kind(&mirror_kind, false)
621        };
622        let order = *element.raw_proper_order();
623        let simplified_symbol = element.get_simplified_symbol();
624        let full_symbol = element.get_full_symbol();
625        let au = if tr {
626            Some(AntiunitaryKind::TimeReversal)
627        } else {
628            None
629        };
630        let is_o3_mirror_plane = element.is_o3_mirror_plane(au);
631        let is_o3_inversion_centre = element.is_o3_inversion_centre(au);
632        let improper_kind = if tr { TRSIG } else { SIG };
633        let result = if generator {
634            if let Vacant(improper_generators) = self.generators.entry(improper_kind) {
635                improper_generators.insert(HashMap::from([(order, IndexSet::from([element]))]));
636                true
637            } else {
638                let improper_generators =
639                    self.generators.get_mut(&improper_kind).unwrap_or_else(|| {
640                        panic!(
641                            "{} generators not found.",
642                            if tr {
643                                "Time-reversed improper"
644                            } else {
645                                "Improper"
646                            }
647                        )
648                    });
649
650                if let Vacant(e) = improper_generators.entry(order) {
651                    e.insert(IndexSet::from([element]));
652                    true
653                } else {
654                    improper_generators
655                        .get_mut(&order)
656                        .unwrap_or_else(|| {
657                            panic!(
658                                "Improper generators {}S{order} not found.",
659                                if tr { "θ" } else { "" }
660                            )
661                        })
662                        .insert(element)
663                }
664            }
665        } else if let Vacant(improper_elements) = self.elements.entry(improper_kind) {
666            improper_elements.insert(HashMap::from([(order, IndexSet::from([element]))]));
667            true
668        } else {
669            let improper_elements = self.elements.get_mut(&improper_kind).unwrap_or_else(|| {
670                panic!(
671                    "{} elements not found.",
672                    if tr {
673                        "Time-reversed improper"
674                    } else {
675                        "Improper"
676                    }
677                )
678            });
679
680            if let Vacant(e) = improper_elements.entry(order) {
681                e.insert(IndexSet::from([element]));
682                true
683            } else {
684                improper_elements
685                    .get_mut(&order)
686                    .unwrap_or_else(|| {
687                        panic!(
688                            "Improper elements {}S{order} not found.",
689                            if tr { "θ" } else { "" }
690                        )
691                    })
692                    .insert(element)
693            }
694        };
695        let dest_str = if generator {
696            "generator".to_owned()
697        } else {
698            "element".to_owned()
699        };
700        if result {
701            if is_o3_mirror_plane {
702                log::debug!(
703                    "Mirror plane {} ({}): {} axis along ({:+.3}, {:+.3}, {:+.3}) added.",
704                    dest_str,
705                    simplified_symbol,
706                    full_symbol,
707                    positive_axis[0] + 0.0,
708                    positive_axis[1] + 0.0,
709                    positive_axis[2] + 0.0,
710                );
711            } else if is_o3_inversion_centre {
712                log::debug!(
713                    "Inversion centre {} ({}): {} axis along ({:+.3}, {:+.3}, {:+.3}) added.",
714                    dest_str,
715                    simplified_symbol,
716                    full_symbol,
717                    positive_axis[0] + 0.0,
718                    positive_axis[1] + 0.0,
719                    positive_axis[2] + 0.0,
720                );
721            } else {
722                log::debug!(
723                    "Improper rotation {} ({}): {} axis along ({:+.3}, {:+.3}, {:+.3}) added.",
724                    dest_str,
725                    simplified_symbol,
726                    full_symbol,
727                    positive_axis[0] + 0.0,
728                    positive_axis[1] + 0.0,
729                    positive_axis[2] + 0.0,
730                );
731            }
732        }
733        result
734    }
735
736    /// Sets the name of the symmetry group.
737    ///
738    /// # Arguments
739    ///
740    /// * `name` - The name to be given to the symmetry group.
741    fn set_group_name(&mut self, name: String) {
742        self.group_name = Some(name);
743        log::debug!(
744            "Symmetry group determined: {}",
745            self.group_name.as_ref().expect("No symmetry groups found.")
746        );
747    }
748
749    /// Obtains elements of a particular kind.
750    ///
751    /// # Arguments
752    ///
753    /// * `kind` - An element kind to be obtained.
754    ///
755    /// # Returns
756    ///
757    /// An optional shared reference to the hash map of the required element kind.
758    #[must_use]
759    pub fn get_elements(
760        &self,
761        kind: &SymmetryElementKind,
762    ) -> Option<&HashMap<ElementOrder, IndexSet<SymmetryElement>>> {
763        self.elements.get(kind)
764    }
765
766    /// Obtains elements of a particular kind (mutable).
767    ///
768    /// # Arguments
769    ///
770    /// * `kind` - An element kind to be obtained.
771    ///
772    /// # Returns
773    ///
774    /// An optional exclusive reference to the hash map of the required element kind.
775    #[must_use]
776    pub fn get_elements_mut(
777        &mut self,
778        kind: &SymmetryElementKind,
779    ) -> Option<&mut HashMap<ElementOrder, IndexSet<SymmetryElement>>> {
780        self.elements.get_mut(kind)
781    }
782
783    /// Obtains generators of a particular kind.
784    ///
785    /// # Arguments
786    ///
787    /// * `kind` - A generator kind to be obtained.
788    ///
789    /// # Returns
790    ///
791    /// An optional shared reference to the hash map of the required generator kind.
792    #[must_use]
793    pub fn get_generators(
794        &self,
795        kind: &SymmetryElementKind,
796    ) -> Option<&HashMap<ElementOrder, IndexSet<SymmetryElement>>> {
797        self.generators.get(kind)
798    }
799
800    /// Obtains generators of a particular kind (mutable).
801    ///
802    /// # Arguments
803    ///
804    /// * `kind` - A generator kind to be obtained.
805    ///
806    /// # Returns
807    ///
808    /// An optional exclusive reference to the hash map of the required generator kind.
809    #[must_use]
810    pub fn get_generators_mut(
811        &mut self,
812        kind: &SymmetryElementKind,
813    ) -> Option<&mut HashMap<ElementOrder, IndexSet<SymmetryElement>>> {
814        self.generators.get_mut(kind)
815    }
816
817    /// Obtains mirror-plane elements by their type (`"h"`, `"v"`, `"d"`, or `""`), including both
818    /// time-reversed and non-time-reversed variants.
819    ///
820    /// # Arguments
821    ///
822    /// * `sigma` - The mirror-plane type which is one of `"h"`, `"v"`, `"d"`, or `""`.
823    ///
824    /// # Returns
825    ///
826    /// An option containing the set of the required mirror-plane element type, if exists. If not,
827    /// then `None` is returned.
828    #[must_use]
829    pub fn get_sigma_elements(&self, sigma: &str) -> Option<HashSet<&SymmetryElement>> {
830        self.get_improper(&ORDER_1)
831            .map(|sigma_elements| {
832                sigma_elements
833                    .iter()
834                    .filter_map(|ele| {
835                        if ele.additional_subscript == sigma {
836                            Some(*ele)
837                        } else {
838                            None
839                        }
840                    })
841                    .collect::<HashSet<_>>()
842            })
843            .filter(|sigma_elements| !sigma_elements.is_empty())
844    }
845
846    /// Obtains mirror-plane generators by their type (`"h"`, `"v"`, `"d"`, or `""`), including both
847    /// time-reversed and non-time-reversed variants.
848    ///
849    /// # Arguments
850    ///
851    /// * `sigma` - The mirror-plane type which is one of `"h"`, `"v"`, `"d"`, or `""`.
852    ///
853    /// # Returns
854    ///
855    /// A set of the required mirror-plane generator type, if exists.
856    #[must_use]
857    pub fn get_sigma_generators(&self, sigma: &str) -> Option<HashSet<&SymmetryElement>> {
858        let mut sigma_generators: HashSet<&SymmetryElement> = HashSet::new();
859        if let Some(improper_generators) = self.get_generators(&SIG) {
860            if let Some(sigmas) = improper_generators.get(&ORDER_1) {
861                sigma_generators.extend(
862                    sigmas
863                        .iter()
864                        .filter(|ele| ele.additional_subscript == sigma),
865                );
866            }
867        }
868        if let Some(tr_improper_generators) = self.get_generators(&TRSIG) {
869            if let Some(sigmas) = tr_improper_generators.get(&ORDER_1) {
870                sigma_generators.extend(
871                    sigmas
872                        .iter()
873                        .filter(|ele| ele.additional_subscript == sigma),
874                );
875            }
876        }
877        if sigma_generators.is_empty() {
878            None
879        } else {
880            Some(sigma_generators)
881        }
882    }
883
884    /// Obtains the highest proper rotation order.
885    ///
886    /// # Returns
887    ///
888    /// The highest proper rotation order.
889    #[must_use]
890    pub fn get_max_proper_order(&self) -> ElementOrder {
891        *self
892            .get_generators(&ROT)
893            .unwrap_or(&HashMap::new())
894            .keys()
895            .chain(self.get_elements(&ROT).unwrap_or(&HashMap::new()).keys())
896            .chain(
897                self.get_generators(&TRROT)
898                    .unwrap_or(&HashMap::new())
899                    .keys(),
900            )
901            .chain(self.get_elements(&TRROT).unwrap_or(&HashMap::new()).keys())
902            .max()
903            .expect("No highest proper rotation order could be obtained.")
904    }
905
906    /// Obtains all proper elements of a certain order (both time-reversed and non-time-reversed).
907    ///
908    /// # Arguments
909    ///
910    /// * `order` - The required order of elements.
911    ///
912    /// # Returns
913    ///
914    /// An optional hash set of proper elements of the required order. If no such elements exist,
915    /// `None` will be returned.
916    #[must_use]
917    pub fn get_proper(&self, order: &ElementOrder) -> Option<HashSet<&SymmetryElement>> {
918        let opt_proper_elements = self
919            .get_elements(&ROT)
920            .map(|proper_elements| proper_elements.get(order))
921            .unwrap_or_default();
922        let opt_tr_proper_elements = self
923            .get_elements(&TRROT)
924            .map(|tr_proper_elements| tr_proper_elements.get(order))
925            .unwrap_or_default();
926
927        match (opt_proper_elements, opt_tr_proper_elements) {
928            (None, None) => None,
929            (Some(proper_elements), None) => Some(proper_elements.iter().collect::<HashSet<_>>()),
930            (None, Some(tr_proper_elements)) => {
931                Some(tr_proper_elements.iter().collect::<HashSet<_>>())
932            }
933            (Some(proper_elements), Some(tr_proper_elements)) => Some(
934                proper_elements
935                    .iter()
936                    .chain(tr_proper_elements.iter())
937                    .collect::<HashSet<_>>(),
938            ),
939        }
940    }
941
942    /// Obtains all improper elements of a certain order (both time-reversed and non-time-reversed).
943    ///
944    /// # Arguments
945    ///
946    /// * `order` - The required order of elements.
947    ///
948    /// # Returns
949    ///
950    /// An optional hash set of improper elements of the required order. If no such elements exist,
951    /// `None` will be returned.
952    #[must_use]
953    pub fn get_improper(&self, order: &ElementOrder) -> Option<HashSet<&SymmetryElement>> {
954        let opt_improper_elements = self
955            .get_elements(&SIG)
956            .map(|improper_elements| improper_elements.get(order))
957            .unwrap_or_default();
958        let opt_tr_improper_elements = self
959            .get_elements(&TRSIG)
960            .map(|tr_improper_elements| tr_improper_elements.get(order))
961            .unwrap_or_default();
962
963        match (opt_improper_elements, opt_tr_improper_elements) {
964            (None, None) => None,
965            (Some(improper_elements), None) => {
966                Some(improper_elements.iter().collect::<HashSet<_>>())
967            }
968            (None, Some(tr_improper_elements)) => {
969                Some(tr_improper_elements.iter().collect::<HashSet<_>>())
970            }
971            (Some(improper_elements), Some(tr_improper_elements)) => Some(
972                improper_elements
973                    .iter()
974                    .chain(tr_improper_elements.iter())
975                    .collect::<HashSet<_>>(),
976            ),
977        }
978    }
979
980    /// Obtains a proper principal element, *i.e.* a time-reversed or non-time-reversed proper
981    /// element with the highest order.
982    ///
983    /// If there are several such elements, the element to be returned will be randomly chosen but
984    /// with any non-time-reversed ones prioritised.
985    ///
986    /// # Returns
987    ///
988    /// A proper principal element.
989    ///
990    /// # Panics
991    ///
992    /// Panics if no proper elements or generators can be found.
993    #[must_use]
994    pub fn get_proper_principal_element(&self) -> &SymmetryElement {
995        let max_ord = self.get_max_proper_order();
996        let principal_elements = self.get_proper(&max_ord).unwrap_or_else(|| {
997            let opt_proper_generators = self
998                .get_generators(&ROT)
999                .map(|proper_generators| proper_generators.get(&max_ord))
1000                .unwrap_or_default();
1001            let opt_tr_proper_generators = self
1002                .get_elements(&TRROT)
1003                .map(|tr_proper_generators| tr_proper_generators.get(&max_ord))
1004                .unwrap_or_default();
1005
1006            match (opt_proper_generators, opt_tr_proper_generators) {
1007                (None, None) => panic!("No proper elements found."),
1008                (Some(proper_generators), None) => proper_generators.iter().collect::<HashSet<_>>(),
1009                (None, Some(tr_proper_generators)) => {
1010                    tr_proper_generators.iter().collect::<HashSet<_>>()
1011                }
1012                (Some(proper_generators), Some(tr_proper_generators)) => proper_generators
1013                    .iter()
1014                    .chain(tr_proper_generators.iter())
1015                    .collect::<HashSet<_>>(),
1016            }
1017        });
1018        principal_elements
1019            .iter()
1020            .find(|ele| !ele.contains_time_reversal())
1021            .unwrap_or_else(|| {
1022                principal_elements
1023                    .iter()
1024                    .next()
1025                    .expect("No proper principal elements found.")
1026            })
1027    }
1028
1029    /// Determines if this group is an infinite group.
1030    ///
1031    /// # Returns
1032    ///
1033    /// A boolean indicating if this group is an infinite group.
1034    #[must_use]
1035    pub fn is_infinite(&self) -> bool {
1036        self.get_max_proper_order() == ElementOrder::Inf
1037            || *self
1038                .get_generators(&ROT)
1039                .unwrap_or(&HashMap::new())
1040                .keys()
1041                .chain(self.get_elements(&ROT).unwrap_or(&HashMap::new()).keys())
1042                .chain(
1043                    self.get_generators(&TRROT)
1044                        .unwrap_or(&HashMap::new())
1045                        .keys(),
1046                )
1047                .chain(self.get_elements(&TRROT).unwrap_or(&HashMap::new()).keys())
1048                .max()
1049                .unwrap_or(&ElementOrder::Int(0))
1050                == ElementOrder::Inf
1051    }
1052
1053    /// Returns the total number of symmetry elements (*NOT* symmetry operations). In
1054    /// infinite-order groups, this is the sum of the number of discrete symmetry elements and the
1055    /// number of discrete symmetry generators.
1056    pub fn n_elements(&self) -> usize {
1057        let n_elements = self
1058            .elements
1059            .values()
1060            .flat_map(|kind_elements| kind_elements.values())
1061            .flatten()
1062            .count();
1063        if self.is_infinite() {
1064            n_elements
1065                + self
1066                    .generators
1067                    .values()
1068                    .flat_map(|kind_elements| kind_elements.values())
1069                    .flatten()
1070                    .count()
1071        } else {
1072            n_elements
1073        }
1074    }
1075
1076    /// Generates all possible symmetry operations from the available symmetry elements.
1077    ///
1078    /// # Arguments
1079    ///
1080    /// * `infinite_order_to_finite` - A finite order to interpret infinite-order generators of
1081    /// infinite groups.
1082    ///
1083    /// # Returns
1084    ///
1085    /// A vector of generated symmetry operations.
1086    ///
1087    /// # Panics
1088    ///
1089    /// Panics if the group is infinite but `infinite_order_to_finite` is `None`, or if the finite
1090    /// order specified in `infinite_order_to_finite` is incompatible with the infinite group.
1091    #[allow(clippy::too_many_lines)]
1092    pub fn generate_all_operations(
1093        &self,
1094        infinite_order_to_finite: Option<u32>,
1095    ) -> Vec<SymmetryOperation> {
1096        let handles_infinite_group = if self.is_infinite() {
1097            assert!(infinite_order_to_finite.is_some());
1098            infinite_order_to_finite
1099        } else {
1100            None
1101        };
1102
1103        if let Some(finite_order) = handles_infinite_group {
1104            let group_name = self.group_name.as_ref().expect("Group name not found.");
1105            if group_name.contains("O(3)") {
1106                if !matches!(finite_order, 2 | 4) {
1107                    log::error!(
1108                        "Finite order of `{finite_order}` is not yet supported for `{group_name}`."
1109                    );
1110                }
1111                assert!(
1112                    matches!(finite_order, 2 | 4),
1113                    "Finite order of `{finite_order}` is not yet supported for `{group_name}`."
1114                );
1115            }
1116        }
1117
1118        let id_element = self
1119            .get_elements(&ROT)
1120            .unwrap_or(&HashMap::new())
1121            .get(&ORDER_1)
1122            .expect("No identity elements found.")
1123            .iter()
1124            .next()
1125            .expect("No identity elements found.")
1126            .clone();
1127
1128        let id_operation = SymmetryOperation::builder()
1129            .generating_element(id_element)
1130            .power(1)
1131            .build()
1132            .expect("Unable to construct an identity operation.");
1133
1134        let empty_elements: HashMap<ElementOrder, IndexSet<SymmetryElement>> = HashMap::new();
1135
1136        // Finite proper operations
1137        let mut proper_orders = self
1138            .get_elements(&ROT)
1139            .unwrap_or(&empty_elements)
1140            .keys()
1141            .collect::<Vec<_>>();
1142        proper_orders.sort_by(|a, b| {
1143            a.partial_cmp(b)
1144                .unwrap_or_else(|| panic!("`{a}` and `{b}` cannot be compared."))
1145        });
1146        let proper_operations =
1147            proper_orders
1148                .iter()
1149                .fold(vec![id_operation], |mut acc, proper_order| {
1150                    self.get_elements(&ROT)
1151                        .unwrap_or(&empty_elements)
1152                        .get(proper_order)
1153                        .unwrap_or_else(|| panic!("Proper elements C{proper_order} not found."))
1154                        .iter()
1155                        .for_each(|proper_element| {
1156                            if let ElementOrder::Int(io) = proper_order {
1157                                acc.extend((1..*io).map(|power| {
1158                                    SymmetryOperation::builder()
1159                                        .generating_element(proper_element.clone())
1160                                        .power(power.try_into().unwrap_or_else(|_| {
1161                                            panic!("Unable to convert `{power}` to `i32`.")
1162                                        }))
1163                                        .build()
1164                                        .expect("Unable to construct a symmetry operation.")
1165                                }));
1166                            }
1167                        });
1168                    acc
1169                });
1170
1171        // Finite proper operations from generators
1172        let proper_operations_from_generators = if let Some(fin_ord) = handles_infinite_group {
1173            self.get_generators(&ROT)
1174                .unwrap_or(&empty_elements)
1175                .par_iter()
1176                .fold(std::vec::Vec::new, |mut acc, (order, proper_generators)| {
1177                    for proper_generator in proper_generators.iter() {
1178                        let finite_order = match order {
1179                            ElementOrder::Int(io) => *io,
1180                            ElementOrder::Inf => fin_ord,
1181                        };
1182                        let finite_proper_element = SymmetryElement::builder()
1183                            .threshold(proper_generator.threshold())
1184                            .proper_order(ElementOrder::Int(finite_order))
1185                            .proper_power(1)
1186                            .raw_axis(*proper_generator.raw_axis())
1187                            .kind(*proper_generator.kind())
1188                            .rotation_group(proper_generator.rotation_group().clone())
1189                            .additional_superscript(proper_generator.additional_superscript.clone())
1190                            .additional_subscript(proper_generator.additional_subscript.clone())
1191                            .build()
1192                            .expect("Unable to construct a symmetry element.");
1193                        acc.extend((1..finite_order).map(|power| {
1194                            SymmetryOperation::builder()
1195                                .generating_element(finite_proper_element.clone())
1196                                .power(power.try_into().unwrap_or_else(|_| {
1197                                    panic!("Unable to convert `{power}` to `i32`.")
1198                                }))
1199                                .build()
1200                                .expect("Unable to construct a symmetry operation.")
1201                        }));
1202                    }
1203                    acc
1204                })
1205                .reduce(std::vec::Vec::new, |mut acc, vec| {
1206                    acc.extend(vec);
1207                    acc
1208                })
1209        } else {
1210            vec![]
1211        };
1212
1213        // Finite time-reversed proper operations
1214        let mut tr_proper_orders = self
1215            .get_elements(&TRROT)
1216            .unwrap_or(&empty_elements)
1217            .keys()
1218            .collect::<Vec<_>>();
1219        tr_proper_orders.sort_by(|a, b| {
1220            a.partial_cmp(b)
1221                .unwrap_or_else(|| panic!("`{a}` and `{b}` cannot be compared."))
1222        });
1223        let tr_proper_operations =
1224            tr_proper_orders
1225                .iter()
1226                .fold(vec![], |mut acc, tr_proper_order| {
1227                    self.get_elements(&TRROT)
1228                        .unwrap_or(&empty_elements)
1229                        .get(tr_proper_order)
1230                        .unwrap_or_else(|| {
1231                            panic!("Proper elements θ·C{tr_proper_order} not found.")
1232                        })
1233                        .iter()
1234                        .for_each(|tr_proper_element| {
1235                            if let ElementOrder::Int(io) = tr_proper_order {
1236                                acc.extend((1..(2 * *io)).step_by(2).map(|power| {
1237                                    SymmetryOperation::builder()
1238                                        .generating_element(tr_proper_element.clone())
1239                                        .power(power.try_into().unwrap_or_else(|_| {
1240                                            panic!("Unable to convert `{power}` to `i32`.")
1241                                        }))
1242                                        .build()
1243                                        .expect("Unable to construct a symmetry operation.")
1244                                }));
1245                            }
1246                        });
1247                    acc
1248                });
1249
1250        // Finite time-reversed proper operations from generators
1251        let tr_proper_operations_from_generators = if let Some(fin_ord) = handles_infinite_group {
1252            self.get_generators(&TRROT)
1253                .unwrap_or(&empty_elements)
1254                .par_iter()
1255                .fold(
1256                    std::vec::Vec::new,
1257                    |mut acc, (order, tr_proper_generators)| {
1258                        for tr_proper_generator in tr_proper_generators.iter() {
1259                            let finite_order = match order {
1260                                ElementOrder::Int(io) => *io,
1261                                ElementOrder::Inf => fin_ord,
1262                            };
1263                            let finite_tr_proper_element = SymmetryElement::builder()
1264                                .threshold(tr_proper_generator.threshold())
1265                                .proper_order(ElementOrder::Int(finite_order))
1266                                .proper_power(1)
1267                                .raw_axis(*tr_proper_generator.raw_axis())
1268                                .kind(*tr_proper_generator.kind())
1269                                .rotation_group(tr_proper_generator.rotation_group().clone())
1270                                .additional_superscript(
1271                                    tr_proper_generator.additional_superscript.clone(),
1272                                )
1273                                .additional_subscript(
1274                                    tr_proper_generator.additional_subscript.clone(),
1275                                )
1276                                .build()
1277                                .expect("Unable to construct a symmetry element.");
1278                            acc.extend((1..finite_order).map(|power| {
1279                                SymmetryOperation::builder()
1280                                    .generating_element(finite_tr_proper_element.clone())
1281                                    .power(power.try_into().unwrap_or_else(|_| {
1282                                        panic!("Unable to convert `{power}` to `i32`.")
1283                                    }))
1284                                    .build()
1285                                    .expect("Unable to construct a symmetry operation.")
1286                            }));
1287                        }
1288                        acc
1289                    },
1290                )
1291                .reduce(std::vec::Vec::new, |mut acc, vec| {
1292                    acc.extend(vec);
1293                    acc
1294                })
1295        } else {
1296            vec![]
1297        };
1298
1299        // Finite improper operations
1300        let mut improper_orders = self
1301            .get_elements(&SIG)
1302            .unwrap_or(&empty_elements)
1303            .keys()
1304            .collect::<Vec<_>>();
1305        improper_orders.sort_by(|a, b| {
1306            a.partial_cmp(b)
1307                .unwrap_or_else(|| panic!("`{a}` and `{b}` cannot be compared."))
1308        });
1309        let improper_operations = improper_orders
1310            .iter()
1311            .fold(vec![], |mut acc, improper_order| {
1312                self.get_elements(&SIG)
1313                    .unwrap_or(&empty_elements)
1314                    .get(improper_order)
1315                    .unwrap_or_else(|| panic!("Improper elements S{improper_order} not found."))
1316                    .iter()
1317                    .for_each(|improper_element| {
1318                        if let ElementOrder::Int(io) = improper_order {
1319                            acc.extend((1..(2 * *io)).step_by(2).map(|power| {
1320                                SymmetryOperation::builder()
1321                                    .generating_element(improper_element.clone())
1322                                    .power(power.try_into().unwrap_or_else(|_| {
1323                                        panic!("Unable to convert `{power}` to `i32`.")
1324                                    }))
1325                                    .build()
1326                                    .expect("Unable to construct a symmetry operation.")
1327                            }));
1328                        }
1329                    });
1330                acc
1331            });
1332
1333        // Finite improper operations from generators
1334        let improper_operations_from_generators = if let Some(fin_ord) = handles_infinite_group {
1335            self.get_generators(&SIG)
1336                .unwrap_or(&empty_elements)
1337                .par_iter()
1338                .fold(
1339                    std::vec::Vec::new,
1340                    |mut acc, (order, improper_generators)| {
1341                        for improper_generator in improper_generators.iter() {
1342                            let finite_order = match order {
1343                                ElementOrder::Int(io) => *io,
1344                                ElementOrder::Inf => fin_ord,
1345                            };
1346                            let finite_improper_element = SymmetryElement::builder()
1347                                .threshold(improper_generator.threshold())
1348                                .proper_order(ElementOrder::Int(finite_order))
1349                                .proper_power(1)
1350                                .raw_axis(*improper_generator.raw_axis())
1351                                .kind(*improper_generator.kind())
1352                                .rotation_group(improper_generator.rotation_group().clone())
1353                                .additional_superscript(
1354                                    improper_generator.additional_superscript.clone(),
1355                                )
1356                                .additional_subscript(
1357                                    improper_generator.additional_subscript.clone(),
1358                                )
1359                                .build()
1360                                .expect("Unable to construct a symmetry element.");
1361                            acc.extend((1..(2 * finite_order)).step_by(2).map(|power| {
1362                                SymmetryOperation::builder()
1363                                    .generating_element(finite_improper_element.clone())
1364                                    .power(power.try_into().unwrap_or_else(|_| {
1365                                        panic!("Unable to convert `{power}` to `i32`.")
1366                                    }))
1367                                    .build()
1368                                    .expect("Unable to construct a symmetry operation.")
1369                            }));
1370                        }
1371                        acc
1372                    },
1373                )
1374                .reduce(std::vec::Vec::new, |mut acc, vec| {
1375                    acc.extend(vec);
1376                    acc
1377                })
1378        } else {
1379            vec![]
1380        };
1381
1382        // Finite time-reversed improper operations
1383        let mut tr_improper_orders = self
1384            .get_elements(&TRSIG)
1385            .unwrap_or(&empty_elements)
1386            .keys()
1387            .collect::<Vec<_>>();
1388        tr_improper_orders.sort_by(|a, b| {
1389            a.partial_cmp(b)
1390                .unwrap_or_else(|| panic!("`{a}` and `{b}` cannot be compared."))
1391        });
1392        let tr_improper_operations =
1393            tr_improper_orders
1394                .iter()
1395                .fold(vec![], |mut acc, tr_improper_order| {
1396                    self.get_elements(&TRSIG)
1397                        .unwrap_or(&empty_elements)
1398                        .get(tr_improper_order)
1399                        .unwrap_or_else(|| {
1400                            panic!("Improper elements θ·S{tr_improper_order} not found.")
1401                        })
1402                        .iter()
1403                        .for_each(|tr_improper_element| {
1404                            if let ElementOrder::Int(io) = tr_improper_order {
1405                                acc.extend((1..(2 * *io)).step_by(2).map(|power| {
1406                                    SymmetryOperation::builder()
1407                                        .generating_element(tr_improper_element.clone())
1408                                        .power(power.try_into().unwrap_or_else(|_| {
1409                                            panic!("Unable to convert `{power}` to `i32`.")
1410                                        }))
1411                                        .build()
1412                                        .expect("Unable to construct a symmetry operation.")
1413                                }));
1414                            }
1415                        });
1416                    acc
1417                });
1418
1419        // Finite time-reversed improper operations from generators
1420        let tr_improper_operations_from_generators = if let Some(fin_ord) = handles_infinite_group {
1421            self.get_generators(&TRSIG)
1422                .unwrap_or(&empty_elements)
1423                .par_iter()
1424                .fold(
1425                    std::vec::Vec::new,
1426                    |mut acc, (order, tr_improper_generators)| {
1427                        for tr_improper_generator in tr_improper_generators.iter() {
1428                            let finite_order = match order {
1429                                ElementOrder::Int(io) => *io,
1430                                ElementOrder::Inf => fin_ord,
1431                            };
1432                            let finite_tr_improper_element = SymmetryElement::builder()
1433                                .threshold(tr_improper_generator.threshold())
1434                                .proper_order(ElementOrder::Int(finite_order))
1435                                .proper_power(1)
1436                                .raw_axis(*tr_improper_generator.raw_axis())
1437                                .kind(*tr_improper_generator.kind())
1438                                .rotation_group(tr_improper_generator.rotation_group().clone())
1439                                .additional_superscript(
1440                                    tr_improper_generator.additional_superscript.clone(),
1441                                )
1442                                .additional_subscript(
1443                                    tr_improper_generator.additional_subscript.clone(),
1444                                )
1445                                .build()
1446                                .expect("Unable to construct a symmetry element.");
1447                            acc.extend((1..(2 * finite_order)).step_by(2).map(|power| {
1448                                SymmetryOperation::builder()
1449                                    .generating_element(finite_tr_improper_element.clone())
1450                                    .power(power.try_into().unwrap_or_else(|_| {
1451                                        panic!("Unable to convert `{power}` to `i32`.")
1452                                    }))
1453                                    .build()
1454                                    .expect("Unable to construct a symmetry operation.")
1455                            }));
1456                        }
1457                        acc
1458                    },
1459                )
1460                .reduce(std::vec::Vec::new, |mut acc, vec| {
1461                    acc.extend(vec);
1462                    acc
1463                })
1464        } else {
1465            vec![]
1466        };
1467
1468        let operations: IndexSet<_> = if handles_infinite_group.is_none() {
1469            proper_operations
1470                .into_iter()
1471                .chain(proper_operations_from_generators)
1472                .chain(improper_operations)
1473                .chain(improper_operations_from_generators)
1474                .chain(tr_proper_operations)
1475                .chain(tr_proper_operations_from_generators)
1476                .chain(tr_improper_operations)
1477                .chain(tr_improper_operations_from_generators)
1478                .collect()
1479        } else {
1480            // Fulfil group closure
1481            log::debug!("Fulfilling closure for a finite subgroup of an infinite group...");
1482            let mut existing_operations: IndexSet<_> = proper_operations
1483                .into_iter()
1484                .chain(proper_operations_from_generators)
1485                .chain(improper_operations)
1486                .chain(improper_operations_from_generators)
1487                .chain(tr_proper_operations)
1488                .chain(tr_proper_operations_from_generators)
1489                .chain(tr_improper_operations)
1490                .chain(tr_improper_operations_from_generators)
1491                .collect();
1492            let mut extra_operations = HashSet::<SymmetryOperation>::new();
1493            let mut npasses = 0;
1494            let mut nstable = 0;
1495
1496            let principal_element = self.get_proper_principal_element();
1497            while nstable < 2 || npasses == 0 {
1498                let n_extra_operations = extra_operations.len();
1499                existing_operations.extend(extra_operations);
1500
1501                npasses += 1;
1502                log::debug!(
1503                    "Generating all group elements: {} pass{}, {} element{} (of which {} {} new)",
1504                    npasses,
1505                    {
1506                        if npasses > 1 {
1507                            "es"
1508                        } else {
1509                            ""
1510                        }
1511                    }
1512                    .to_string(),
1513                    existing_operations.len(),
1514                    {
1515                        if existing_operations.len() > 1 {
1516                            "s"
1517                        } else {
1518                            ""
1519                        }
1520                    }
1521                    .to_string(),
1522                    n_extra_operations,
1523                    {
1524                        if n_extra_operations == 1 {
1525                            "is"
1526                        } else {
1527                            "are"
1528                        }
1529                    }
1530                    .to_string(),
1531                );
1532
1533                extra_operations = existing_operations
1534                    .iter()
1535                    .combinations_with_replacement(2)
1536                    .par_bridge()
1537                    .filter_map(|op_pairs| {
1538                        let op_i_ref = op_pairs[0];
1539                        let op_j_ref = op_pairs[1];
1540                        let op_k = op_i_ref * op_j_ref;
1541                        if existing_operations.contains(&op_k) {
1542                            None
1543                        } else if op_k.is_proper() {
1544                            Some(op_k)
1545                        } else if op_k.is_spatial_reflection()
1546                            && op_k.generating_element.additional_subscript.is_empty()
1547                        {
1548                            if let Some(sigma_symbol) = deduce_sigma_symbol(
1549                                op_k.generating_element.raw_axis(),
1550                                principal_element,
1551                                op_k.generating_element.threshold(),
1552                                false,
1553                            ) {
1554                                let mut op_k_sym = op_k.convert_to_improper_kind(&SIG);
1555                                op_k_sym.generating_element.additional_subscript = sigma_symbol;
1556                                Some(op_k_sym)
1557                            } else {
1558                                Some(op_k.convert_to_improper_kind(&SIG))
1559                            }
1560                        } else {
1561                            Some(op_k.convert_to_improper_kind(&SIG))
1562                        }
1563                    })
1564                    .collect();
1565                if extra_operations.is_empty() {
1566                    nstable += 1;
1567                } else {
1568                    nstable = 0;
1569                }
1570            }
1571
1572            assert_eq!(extra_operations.len(), 0);
1573            log::debug!(
1574                "Group closure reached with {} elements.",
1575                existing_operations.len()
1576            );
1577            existing_operations
1578        };
1579
1580        let mut sorted_operations: Vec<SymmetryOperation> = operations.into_iter().collect();
1581        sort_operations(&mut sorted_operations);
1582        sorted_operations
1583    }
1584}
1585
1586impl Default for Symmetry {
1587    fn default() -> Self {
1588        Self::new()
1589    }
1590}
1591
1592/// Locates all proper rotation elements present in [`PreSymmetry::recentred_molecule`]
1593///
1594/// # Arguments
1595///
1596/// * `presym` - A pre-symmetry-analysis struct containing information about
1597/// the molecular system.
1598/// * `sym` - A symmetry struct to store the proper rotation elements found.
1599/// * `asymmetric` - If `true`, the search assumes that the group is one of the
1600/// Abelian point groups for which the highest possible rotation order is $`2`$
1601/// and there can be at most three $`C_2`$ axes.
1602/// * `tr` - A flag indicating if time reversal should also be considered.
1603#[allow(clippy::too_many_lines)]
1604fn _search_proper_rotations(
1605    presym: &PreSymmetry,
1606    sym: &mut Symmetry,
1607    asymmetric: bool,
1608    tr: bool,
1609) -> Result<(), anyhow::Error> {
1610    log::debug!("==============================");
1611    log::debug!("Proper rotation search begins.");
1612    log::debug!("==============================");
1613    let mut linear_sea_groups: Vec<&Vec<Atom>> = vec![];
1614    let mut count_c2: usize = 0;
1615    log::debug!("++++++++++++++++++++++++++");
1616    log::debug!("SEA group analysis begins.");
1617    log::debug!("++++++++++++++++++++++++++");
1618    for sea_group in &presym.sea_groups {
1619        if asymmetric && count_c2 == 3 {
1620            break;
1621        }
1622        let k_sea = sea_group.len();
1623        match k_sea {
1624            1 => {
1625                continue;
1626            }
1627            2 => {
1628                log::debug!("A linear SEA set detected: {:?}.", sea_group);
1629                linear_sea_groups.push(sea_group);
1630            }
1631            _ => {
1632                let sea_mol = Molecule::from_atoms(sea_group, presym.dist_threshold);
1633                let (sea_mois, sea_axes) = sea_mol.calc_moi();
1634                // Search for high-order rotation axes
1635                if approx::relative_eq!(
1636                    sea_mois[0] + sea_mois[1],
1637                    sea_mois[2],
1638                    epsilon = presym.moi_threshold,
1639                    max_relative = presym.moi_threshold,
1640                ) {
1641                    // Planar SEA
1642                    let k_fac_range: Vec<_> = if approx::relative_eq!(
1643                        sea_mois[0],
1644                        sea_mois[1],
1645                        epsilon = presym.moi_threshold,
1646                        max_relative = presym.moi_threshold,
1647                    ) {
1648                        // Regular k-sided polygon
1649                        log::debug!(
1650                            "A regular {}-sided polygon SEA set detected: {:?}.",
1651                            k_sea,
1652                            sea_group
1653                        );
1654                        let mut divisors = divisors::get_divisors(k_sea);
1655                        divisors.push(k_sea);
1656                        divisors
1657                    } else {
1658                        // Irregular k-sided polygon
1659                        log::debug!(
1660                            "An irregular {}-sided polygon SEA set detected: {:?}.",
1661                            k_sea,
1662                            sea_group
1663                        );
1664                        divisors::get_divisors(k_sea)
1665                    };
1666                    for k_fac in &k_fac_range {
1667                        if let Some(proper_kind) =
1668                            presym.check_proper(
1669                                &ElementOrder::Int((*k_fac).try_into().unwrap_or_else(|_| {
1670                                    panic!("Unable to convert {k_fac} to `u32`.")
1671                                })),
1672                                &sea_axes[2],
1673                                tr,
1674                            )
1675                        {
1676                            match *k_fac {
1677                                2 => {
1678                                    count_c2 += usize::from(sym.add_proper(
1679                                        ElementOrder::Int((*k_fac).try_into().unwrap_or_else(
1680                                            |_| panic!("Unable to convert {k_fac} to `u32`."),
1681                                        )),
1682                                        &sea_axes[2],
1683                                        false,
1684                                        presym.dist_threshold,
1685                                        proper_kind.contains_time_reversal(),
1686                                    ));
1687                                }
1688                                _ => {
1689                                    sym.add_proper(
1690                                        ElementOrder::Int((*k_fac).try_into().unwrap_or_else(
1691                                            |_| panic!("Unable to convert {k_fac} to `u32`."),
1692                                        )),
1693                                        &sea_axes[2],
1694                                        false,
1695                                        presym.dist_threshold,
1696                                        proper_kind.contains_time_reversal(),
1697                                    );
1698                                }
1699                            }
1700                        }
1701                    }
1702                } else {
1703                    // Polyhedral SEA
1704                    if approx::relative_eq!(
1705                        sea_mois[1],
1706                        sea_mois[2],
1707                        epsilon = presym.moi_threshold,
1708                        max_relative = presym.moi_threshold,
1709                    ) {
1710                        // The number of atoms in this SEA group must be even.
1711                        ensure!(
1712                            k_sea % 2 == 0,
1713                            "Unexpected odd number of atoms in this SEA group."
1714                        );
1715                        if approx::relative_eq!(
1716                            sea_mois[0],
1717                            sea_mois[1],
1718                            epsilon = presym.moi_threshold,
1719                            max_relative = presym.moi_threshold,
1720                        ) {
1721                            // Spherical top SEA
1722                            log::debug!("A spherical top SEA set detected: {:?}", sea_group);
1723                            let sea_presym = PreSymmetry::builder()
1724                                .moi_threshold(presym.moi_threshold)
1725                                .molecule(&sea_mol)
1726                                .build()
1727                                .expect("Unable to construct a `PreSymmetry` structure.");
1728                            let mut sea_sym = Symmetry::builder()
1729                                .build()
1730                                .expect("Unable to construct a default `Symmetry` structure.");
1731                            log::debug!("-----------------------------------------------");
1732                            log::debug!("Symmetry analysis for spherical top SEA begins.");
1733                            log::debug!("-----------------------------------------------");
1734                            sea_sym.analyse(&sea_presym, tr)?;
1735                            log::debug!("---------------------------------------------");
1736                            log::debug!("Symmetry analysis for spherical top SEA ends.");
1737                            log::debug!("---------------------------------------------");
1738                            for (order, proper_elements) in sea_sym
1739                                .get_elements(&ROT)
1740                                .unwrap_or(&HashMap::new())
1741                                .iter()
1742                                .chain(
1743                                    sea_sym
1744                                        .get_elements(&TRROT)
1745                                        .unwrap_or(&HashMap::new())
1746                                        .iter(),
1747                                )
1748                            {
1749                                for proper_element in proper_elements {
1750                                    if let Some(proper_kind) =
1751                                        presym.check_proper(order, proper_element.raw_axis(), tr)
1752                                    {
1753                                        sym.add_proper(
1754                                            *order,
1755                                            proper_element.raw_axis(),
1756                                            false,
1757                                            presym.dist_threshold,
1758                                            proper_kind.contains_time_reversal(),
1759                                        );
1760                                    }
1761                                }
1762                            }
1763                        } else {
1764                            // Prolate symmetric top
1765                            log::debug!("A prolate symmetric top SEA set detected.");
1766                            for k_fac in divisors::get_divisors(k_sea / 2)
1767                                .iter()
1768                                .chain(vec![k_sea / 2].iter())
1769                            {
1770                                let k_fac_order =
1771                                    ElementOrder::Int((*k_fac).try_into().unwrap_or_else(|_| {
1772                                        panic!("Unable to convert {k_fac} to u32.")
1773                                    }));
1774                                if let Some(proper_kind) =
1775                                    presym.check_proper(&k_fac_order, &sea_axes[0], tr)
1776                                {
1777                                    if *k_fac == 2 {
1778                                        count_c2 += usize::from(sym.add_proper(
1779                                            k_fac_order,
1780                                            &sea_axes[0],
1781                                            false,
1782                                            presym.dist_threshold,
1783                                            proper_kind.contains_time_reversal(),
1784                                        ));
1785                                    } else {
1786                                        sym.add_proper(
1787                                            k_fac_order,
1788                                            &sea_axes[0],
1789                                            false,
1790                                            presym.dist_threshold,
1791                                            proper_kind.contains_time_reversal(),
1792                                        );
1793                                    }
1794                                }
1795                            }
1796                        }
1797                    } else if approx::relative_eq!(
1798                        sea_mois[0],
1799                        sea_mois[1],
1800                        epsilon = presym.moi_threshold,
1801                        max_relative = presym.moi_threshold,
1802                    ) {
1803                        // Oblate symmetry top
1804                        log::debug!("An oblate symmetric top SEA set detected.");
1805                        ensure!(
1806                            k_sea % 2 == 0,
1807                            "Unexpected odd number of atoms in this SEA group."
1808                        );
1809                        for k_fac in divisors::get_divisors(k_sea / 2)
1810                            .iter()
1811                            .chain(vec![k_sea / 2].iter())
1812                        {
1813                            let k_fac_order =
1814                                ElementOrder::Int((*k_fac).try_into().map_err(|_| {
1815                                    format_err!("Unable to convert `{k_fac}` to `u32`.")
1816                                })?);
1817                            if let Some(proper_kind) =
1818                                presym.check_proper(&k_fac_order, &sea_axes[2], tr)
1819                            {
1820                                if *k_fac == 2 {
1821                                    count_c2 += usize::from(sym.add_proper(
1822                                        k_fac_order,
1823                                        &sea_axes[2],
1824                                        false,
1825                                        presym.dist_threshold,
1826                                        proper_kind.contains_time_reversal(),
1827                                    ));
1828                                } else {
1829                                    sym.add_proper(
1830                                        k_fac_order,
1831                                        &sea_axes[2],
1832                                        false,
1833                                        presym.dist_threshold,
1834                                        proper_kind.contains_time_reversal(),
1835                                    );
1836                                }
1837                            }
1838                        }
1839                    } else {
1840                        // Asymmetric top
1841                        log::debug!("An asymmetric top SEA set detected.");
1842                        for sea_axis in &sea_axes {
1843                            if let Some(proper_kind) = presym.check_proper(&ORDER_2, sea_axis, tr) {
1844                                count_c2 += usize::from(sym.add_proper(
1845                                    ORDER_2,
1846                                    sea_axis,
1847                                    false,
1848                                    presym.dist_threshold,
1849                                    proper_kind.contains_time_reversal(),
1850                                ));
1851                            }
1852                        }
1853                    }
1854                }
1855            }
1856        } // end match k_sea
1857
1858        // Search for any remaining C2 axes
1859        log::debug!("Searching for any remaining C2 axes in this SEA group...");
1860        for atom2s in sea_group.iter().combinations(2) {
1861            if asymmetric && count_c2 == 3 {
1862                log::debug!("Three C2 axes have been found for an asymmetric top. No more C2 axes will be found. The C2 search has completed.");
1863                break;
1864            }
1865            let atom_i_pos = atom2s[0].coordinates;
1866            let atom_j_pos = atom2s[1].coordinates;
1867
1868            // Case B: C2 might cross through any two atoms
1869            if let Some(proper_kind) = presym.check_proper(&ORDER_2, &atom_i_pos.coords, tr) {
1870                count_c2 += usize::from(sym.add_proper(
1871                    ORDER_2,
1872                    &atom_i_pos.coords,
1873                    false,
1874                    presym.dist_threshold,
1875                    proper_kind.contains_time_reversal(),
1876                ));
1877            }
1878
1879            // Case A: C2 might cross through the midpoint of two atoms
1880            // Care needs to be taken if the midpoint is at the centre of mass of the molecule,
1881            // because then `midvec` would be a zero vector and `c2_check` would be `None`,
1882            // regardless of whether there is actually a C2 axis passing through this midpoint or
1883            // not.
1884            let midvec = 0.5 * (atom_i_pos.coords + atom_j_pos.coords);
1885            let c2_check = presym.check_proper(&ORDER_2, &midvec, tr);
1886            if midvec.norm() > presym.dist_threshold && c2_check.is_some() {
1887                count_c2 += usize::from(
1888                    sym.add_proper(
1889                        ORDER_2,
1890                        &midvec,
1891                        false,
1892                        presym.dist_threshold,
1893                        c2_check
1894                            .expect("Expected C2 not found.")
1895                            .contains_time_reversal(),
1896                    ),
1897                );
1898            } else if let Some(electric_atoms) = &presym.recentred_molecule.electric_atoms {
1899                let com = presym.recentred_molecule.calc_com();
1900                let e_vector = electric_atoms[0].coordinates - com;
1901                if let Some(proper_kind) = presym.check_proper(&ORDER_2, &e_vector, tr) {
1902                    count_c2 += usize::from(sym.add_proper(
1903                        ORDER_2,
1904                        &e_vector,
1905                        false,
1906                        presym.dist_threshold,
1907                        proper_kind.contains_time_reversal(),
1908                    ));
1909                }
1910            } else if midvec.norm() <= presym.dist_threshold {
1911                for principal_axis in &presym.recentred_molecule.calc_moi().1 {
1912                    if let Some(proper_kind) = presym.check_proper(&ORDER_2, principal_axis, tr) {
1913                        count_c2 += usize::from(sym.add_proper(
1914                            ORDER_2,
1915                            principal_axis,
1916                            false,
1917                            presym.dist_threshold,
1918                            proper_kind.contains_time_reversal(),
1919                        ));
1920                    }
1921                }
1922            }
1923        }
1924        log::debug!("Searching for any remaining C2 axes in this SEA group... Done.");
1925    } // end for sea_group in presym.sea_groups.iter()
1926    log::debug!("++++++++++++++++++++++++");
1927    log::debug!("SEA group analysis ends.");
1928    log::debug!("++++++++++++++++++++++++");
1929
1930    // Search for any remaining C2 axes.
1931    if asymmetric && count_c2 == 3 {
1932    } else {
1933        // Case C: Molecules with two or more sets of non-parallel linear diatomic SEA groups
1934        if linear_sea_groups.len() >= 2 {
1935            let normal_option = linear_sea_groups.iter().combinations(2).find_map(|pair| {
1936                let vec_0 = pair[0][1].coordinates - pair[0][0].coordinates;
1937                let vec_1 = pair[1][1].coordinates - pair[1][0].coordinates;
1938                let trial_normal = vec_0.cross(&vec_1);
1939                if trial_normal.norm() > presym.dist_threshold {
1940                    Some(trial_normal)
1941                } else {
1942                    None
1943                }
1944            });
1945            if let Some(normal) = normal_option {
1946                if let Some(proper_kind) = presym.check_proper(&ORDER_2, &normal, tr) {
1947                    sym.add_proper(
1948                        ORDER_2,
1949                        &normal,
1950                        false,
1951                        presym.dist_threshold,
1952                        proper_kind.contains_time_reversal(),
1953                    );
1954                }
1955            }
1956        }
1957    }
1958    log::debug!("============================");
1959    log::debug!("Proper rotation search ends.");
1960    log::debug!("============================");
1961
1962    Ok(())
1963}
1964
1965mod symmetry_core_asymmetric;
1966mod symmetry_core_linear;
1967mod symmetry_core_spherical;
1968mod symmetry_core_symmetric;