1use 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#[derive(Clone, Builder, Debug, Serialize, Deserialize)]
49pub struct PreSymmetry {
50 #[builder(setter(custom))]
52 pub(crate) original_molecule: Molecule,
53
54 #[builder(setter(custom))]
56 pub(crate) recentred_molecule: Molecule,
57
58 #[builder(setter(custom))]
60 pub(crate) moi_threshold: f64,
61
62 #[builder(setter(skip), default = "self.calc_rotational_symmetry()")]
65 pub(crate) rotational_symmetry: RotationalSymmetry,
66
67 #[builder(setter(skip), default = "self.calc_sea_groups()")]
69 pub(crate) sea_groups: Vec<Vec<Atom>>,
70
71 #[builder(setter(skip), default = "self.get_dist_threshold()")]
73 pub(crate) dist_threshold: f64,
74}
75
76impl PreSymmetryBuilder {
77 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 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 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 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 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 #[must_use]
171 pub(crate) fn builder() -> PreSymmetryBuilder {
172 PreSymmetryBuilder::default()
173 }
174
175 #[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 #[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#[derive(Builder, Clone, Debug, Serialize, Deserialize)]
279pub struct Symmetry {
280 #[builder(setter(skip, strip_option), default = "None")]
282 pub group_name: Option<String>,
283
284 #[builder(setter(skip), default = "HashMap::new()")]
292 pub(crate) elements:
293 HashMap<SymmetryElementKind, HashMap<ElementOrder, IndexSet<SymmetryElement>>>,
294
295 #[builder(setter(skip), default = "HashMap::new()")]
303 pub(crate) generators:
304 HashMap<SymmetryElementKind, HashMap<ElementOrder, IndexSet<SymmetryElement>>>,
305}
306
307impl Symmetry {
308 #[must_use]
314 fn builder() -> SymmetryBuilder {
315 SymmetryBuilder::default()
316 }
317
318 #[must_use]
320 pub fn new() -> Self {
321 Symmetry::builder()
322 .build()
323 .expect("Unable to construct a `Symmetry` structure.")
324 }
325
326 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 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 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 if presym.recentred_molecule == presym.recentred_molecule.reverse_time() {
378 log::debug!("Time reversal is a symmetry element. This is a grey group.");
379 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 })
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 self.add_proper(ORDER_1, c1.raw_axis(), true, presym.dist_threshold, true);
403
404 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 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 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 #[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 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 #[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 #[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 #[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 #[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 #[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 #[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 #[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 #[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 #[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 #[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 #[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 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 #[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 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 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 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 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 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 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 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 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 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#[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 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 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 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 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 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 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 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 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 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 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 } 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 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 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 } log::debug!("++++++++++++++++++++++++");
1927 log::debug!("SEA group analysis ends.");
1928 log::debug!("++++++++++++++++++++++++");
1929
1930 if asymmetric && count_c2 == 3 {
1932 } else {
1933 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;