1use std::fmt;
4use std::path::PathBuf;
5
6use anyhow::{bail, format_err};
7use derive_builder::Builder;
8use itertools::Itertools;
9use log;
10use nalgebra::{Point3, Vector3};
11use serde::{Deserialize, Serialize};
12
13use crate::auxiliary::atom::{Atom, AtomKind};
14use crate::auxiliary::molecule::Molecule;
15use crate::drivers::QSym2Driver;
16use crate::io::format::{
17 log_subtitle, log_title, nice_bool, qsym2_output, write_subtitle, QSym2Output,
18};
19use crate::io::{write_qsym2_binary, QSym2FileType};
20use crate::symmetry::symmetry_core::{PreSymmetry, Symmetry};
21use crate::symmetry::symmetry_element::{AntiunitaryKind, SymmetryElementKind};
22
23#[cfg(test)]
24#[path = "symmetry_group_detection_tests.rs"]
25mod symmetry_group_detection_tests;
26
27fn default_thresholds() -> Vec<f64> {
36 vec![1.0e-4, 1e-5, 1.0e-6]
37}
38
39#[derive(Clone, Builder, Debug, Serialize, Deserialize)]
41pub struct SymmetryGroupDetectionParams {
42 #[builder(setter(custom), default = "vec![1.0e-4, 1.0e-5, 1.0e-6]")]
44 #[serde(default = "default_thresholds")]
45 pub moi_thresholds: Vec<f64>,
46
47 #[builder(setter(custom), default = "vec![1.0e-4, 1.0e-5, 1.0e-6]")]
49 #[serde(default = "default_thresholds")]
50 pub distance_thresholds: Vec<f64>,
51
52 #[builder(default = "false")]
54 #[serde(default)]
55 pub time_reversal: bool,
56
57 #[builder(default = "None")]
62 #[serde(default)]
63 pub fictitious_magnetic_fields: Option<Vec<(Point3<f64>, Vector3<f64>)>>,
64
65 #[builder(default = "None")]
69 #[serde(default)]
70 pub fictitious_electric_fields: Option<Vec<(Point3<f64>, Vector3<f64>)>>,
71
72 #[builder(default = "false")]
76 #[serde(default)]
77 pub field_origin_com: bool,
78
79 #[builder(default = "false")]
82 #[serde(default)]
83 pub write_symmetry_elements: bool,
84
85 #[builder(default = "None")]
88 #[serde(default)]
89 pub result_save_name: Option<PathBuf>,
90}
91
92impl SymmetryGroupDetectionParams {
93 pub fn builder() -> SymmetryGroupDetectionParamsBuilder {
95 SymmetryGroupDetectionParamsBuilder::default()
96 }
97}
98
99impl SymmetryGroupDetectionParamsBuilder {
100 pub fn moi_thresholds(&mut self, threshs: &[f64]) -> &mut Self {
101 self.moi_thresholds = Some(threshs.to_vec());
102 self
103 }
104
105 pub fn distance_thresholds(&mut self, threshs: &[f64]) -> &mut Self {
106 self.distance_thresholds = Some(threshs.to_vec());
107 self
108 }
109}
110
111impl Default for SymmetryGroupDetectionParams {
112 fn default() -> Self {
113 Self::builder()
114 .build()
115 .expect("Unable to construct a default `SymmetryGroupDetectionParams`.")
116 }
117}
118
119impl fmt::Display for SymmetryGroupDetectionParams {
120 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
121 let threshs = self
122 .moi_thresholds
123 .iter()
124 .cartesian_product(self.distance_thresholds.iter());
125 let nthreshs = threshs.clone().count();
126 if nthreshs == 1 {
127 writeln!(f, "Fixed thresholds:")?;
128 writeln!(f, " MoI threshold: {:.3e}", self.moi_thresholds[0])?;
129 writeln!(f, " Geo threshold: {:.3e}", self.distance_thresholds[0])?;
130 } else {
131 writeln!(f, "Variable thresholds:")?;
132 writeln!(
133 f,
134 " MoI thresholds: {}",
135 self.moi_thresholds
136 .iter()
137 .map(|v| format!("{v:.3e}"))
138 .join(", ")
139 )?;
140 writeln!(
141 f,
142 " Geo thresholds: {}",
143 self.distance_thresholds
144 .iter()
145 .map(|v| format!("{v:.3e}"))
146 .join(", ")
147 )?;
148 writeln!(f)?;
149 }
150
151 if self.fictitious_magnetic_fields.is_some() || self.fictitious_electric_fields.is_some() {
152 if self.field_origin_com {
153 writeln!(f, "Field origins relative to: molecule's centre of mass")?;
154 } else {
155 writeln!(f, "Field origins relative to: space-fixed origin")?;
156 }
157 }
158
159 if let Some(fictitious_magnetic_fields) = self.fictitious_magnetic_fields.as_ref() {
160 writeln!(f, "Fictitious magnetic fields:")?;
161 for (origin, field) in fictitious_magnetic_fields.iter() {
162 writeln!(
163 f,
164 " ({}) ± ({})",
165 origin.iter().map(|x| format!("{x:+.3}")).join(", "),
166 field.iter().map(|x| format!("{x:+.3}")).join(", "),
167 )?;
168 }
169 writeln!(f)?;
170 }
171
172 if let Some(fictitious_electric_fields) = self.fictitious_electric_fields.as_ref() {
173 writeln!(f, "Fictitious electric fields:")?;
174 for (origin, field) in fictitious_electric_fields.iter() {
175 writeln!(
176 f,
177 " ({}) + ({})",
178 origin.iter().map(|x| format!("{x:+.3}")).join(", "),
179 field.iter().map(|x| format!("{x:+.3}")).join(", "),
180 )?;
181 }
182 writeln!(f)?;
183 }
184
185 writeln!(
186 f,
187 "Consider time reversal: {}",
188 nice_bool(self.time_reversal)
189 )?;
190 writeln!(
191 f,
192 "Report symmetry elements/generators: {}",
193 nice_bool(self.write_symmetry_elements)
194 )?;
195 writeln!(
196 f,
197 "Save symmetry-group detection results to file: {}",
198 if let Some(name) = self.result_save_name.as_ref() {
199 let mut path = name.clone();
200 path.set_extension(QSym2FileType::Sym.ext());
201 path.display().to_string()
202 } else {
203 nice_bool(false)
204 }
205 )?;
206 writeln!(f)?;
207
208 Ok(())
209 }
210}
211
212#[derive(Clone, Builder, Debug, Serialize, Deserialize)]
218pub struct SymmetryGroupDetectionResult {
219 pub parameters: SymmetryGroupDetectionParams,
221
222 pub pre_symmetry: PreSymmetry,
225
226 pub unitary_symmetry: Symmetry,
228
229 #[builder(default = "None")]
232 pub magnetic_symmetry: Option<Symmetry>,
233}
234
235impl SymmetryGroupDetectionResult {
236 fn builder() -> SymmetryGroupDetectionResultBuilder {
238 SymmetryGroupDetectionResultBuilder::default()
239 }
240
241 fn write_symmetry_elements(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
244 if let Some(magnetic_symmetry) = self.magnetic_symmetry.as_ref() {
245 write_subtitle(
246 f,
247 &format!(
248 "Symmetry element report for magnetic group {}",
249 magnetic_symmetry
250 .group_name
251 .as_ref()
252 .unwrap_or(&"?".to_string()),
253 ),
254 )?;
255 writeln!(f)?;
256 write_element_table(f, magnetic_symmetry)?;
257 writeln!(f)?;
258 }
259
260 write_subtitle(
261 f,
262 &format!(
263 "Symmetry element report for unitary group {}",
264 self.unitary_symmetry
265 .group_name
266 .as_ref()
267 .unwrap_or(&"?".to_string())
268 ),
269 )?;
270 writeln!(f)?;
271 write_element_table(f, &self.unitary_symmetry)?;
272 Ok(())
273 }
274}
275
276impl fmt::Display for SymmetryGroupDetectionResult {
277 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
278 if let Some(highest_mag_sym) = self.magnetic_symmetry.as_ref() {
279 let n_mag_elements = if highest_mag_sym.is_infinite() {
280 "∞".to_string()
281 } else {
282 highest_mag_sym.n_elements().to_string()
283 };
284 writeln!(
285 f,
286 "Highest mag. group found: {} ({} {})",
287 highest_mag_sym
288 .group_name
289 .as_ref()
290 .unwrap_or(&"?".to_string()),
291 n_mag_elements,
292 if n_mag_elements != "1" {
293 "symmetry elements"
294 } else {
295 "symmetry element"
296 }
297 )?;
298 }
299
300 let n_uni_elements = if self.unitary_symmetry.is_infinite() {
301 "∞".to_string()
302 } else {
303 self.unitary_symmetry.n_elements().to_string()
304 };
305 writeln!(
306 f,
307 "Highest uni. group found: {} ({} {})",
308 self.unitary_symmetry
309 .group_name
310 .as_ref()
311 .unwrap_or(&"?".to_string()),
312 n_uni_elements,
313 if n_uni_elements != "1" {
314 "symmetry elements"
315 } else {
316 "symmetry element"
317 }
318 )?;
319 writeln!(
320 f,
321 " Associated MoI threshold: {:.3e}",
322 self.pre_symmetry.moi_threshold
323 )?;
324 writeln!(
325 f,
326 " Associated geo threshold: {:.3e}",
327 self.pre_symmetry.dist_threshold
328 )?;
329 writeln!(f)?;
330
331 if self.parameters.write_symmetry_elements {
332 self.write_symmetry_elements(f)?;
333 }
334
335 Ok(())
336 }
337}
338
339#[derive(Clone, Builder)]
345pub struct SymmetryGroupDetectionDriver<'a> {
346 parameters: &'a SymmetryGroupDetectionParams,
348
349 #[builder(default = "None")]
352 xyz: Option<PathBuf>,
353
354 #[builder(default = "None")]
356 molecule: Option<&'a Molecule>,
357
358 #[builder(setter(skip), default = "None")]
360 result: Option<SymmetryGroupDetectionResult>,
361}
362
363impl<'a> SymmetryGroupDetectionDriver<'a> {
364 pub fn builder() -> SymmetryGroupDetectionDriverBuilder<'a> {
366 SymmetryGroupDetectionDriverBuilder::default()
367 }
368
369 fn detect_symmetry_group(&mut self) -> Result<(), anyhow::Error> {
371 log_title("Symmetry-Group Detection");
372 qsym2_output!("");
373 let params = self.parameters;
374 params.log_output_display();
375
376 let smallest_dist_thresh = *params
377 .distance_thresholds
378 .iter()
379 .min_by(|x, y| {
380 x.partial_cmp(y)
381 .expect("Unable to determine the smallest distance threshold.")
382 })
383 .ok_or_else(|| format_err!("Unable to determine the smallest distance threshold."))?;
384 let target_mol = match (self.molecule, self.xyz.as_ref()) {
385 (Some(molecule), None) => Molecule::from_atoms(
386 &molecule.get_all_atoms().into_iter().cloned().collect_vec(),
387 smallest_dist_thresh,
388 ),
389 (None, Some(xyz)) => Molecule::from_xyz(xyz, smallest_dist_thresh),
390 _ => bail!("Neither or both `molecule` and `xyz` are specified."),
391 };
392 qsym2_output!("Molecule for symmetry-group detection:");
393 target_mol.log_output_display();
394 qsym2_output!("");
395
396 let threshs = params
397 .moi_thresholds
398 .iter()
399 .cartesian_product(params.distance_thresholds.iter());
400 let nthreshs = threshs.clone().count();
401
402 log_subtitle("Threshold-scanning symmetry-group detection");
403 qsym2_output!("");
404
405 let count_length = usize::try_from(nthreshs.ilog10() + 2).map_err(|_| {
406 format_err!("Unable to convert `{}` to `usize`.", nthreshs.ilog10() + 2)
407 })?;
408 qsym2_output!("{}", "┈".repeat(count_length + 75));
409 qsym2_output!(
410 "{:>width$} {:>12} {:>12} {:>14} {:>9} {:>12} {:>9}",
411 "#",
412 "MoI thresh",
413 "Geo thresh",
414 "Mag. group",
415 "Elements",
416 "Uni. group",
417 "Elements",
418 width = count_length
419 );
420 qsym2_output!("{}", "┈".repeat(count_length + 75));
421 let mut i = 0;
422 let syms = threshs.map(|(moi_thresh, dist_thresh)| {
423 let mut mol = match (self.molecule, self.xyz.as_ref()) {
425 (Some(molecule), None) => Molecule::from_atoms(
426 &molecule.get_all_atoms().into_iter().cloned().collect_vec(),
427 *dist_thresh
428 ),
429 (None, Some(xyz)) => Molecule::from_xyz(
430 xyz,
431 *dist_thresh
432 ),
433 _ => {
434 qsym2_output!("Neither or both `molecule` and `xyz` are specified.");
435 bail!("Neither or both `molecule` and `xyz` are specified.")
436 }
437 };
438
439 let global_origin = if params.field_origin_com {
441 mol.calc_com() - Point3::origin()
442 } else {
443 Vector3::zeros()
444 };
445 if let Some(fictitious_magnetic_fields) = params.fictitious_magnetic_fields.as_ref() {
446 if mol.magnetic_atoms.is_some() {
447 log::error!("Magnetic fields already present. Fictitious magnetic fields cannot be added.");
448 bail!("Magnetic fields already present. Fictitious magnetic fields cannot be added.")
449 } else {
450 let magnetic_atoms = fictitious_magnetic_fields.iter().flat_map(|(origin, vec)| {
451 Ok::<[Atom; 2], anyhow::Error>([
452 Atom::new_special(AtomKind::Magnetic(true), origin + global_origin + vec, *dist_thresh).ok_or_else(||
453 format_err!("Cannot construct a fictitious magnetic atom.")
454 )?,
455 Atom::new_special(AtomKind::Magnetic(false), origin + global_origin - vec, *dist_thresh).ok_or_else(||
456 format_err!("Cannot construct a fictitious magnetic atom.")
457 )?,
458 ])
459 }).flatten().collect_vec();
460 mol.magnetic_atoms = Some(magnetic_atoms);
461 }
462 }
463
464 if let Some(fictitious_electric_fields) = params.fictitious_electric_fields.as_ref() {
466 if mol.electric_atoms.is_some() {
467 log::error!("Electric fields already present. Fictitious electric fields cannot be added.");
468 bail!("Electric fields already present. Fictitious electric fields cannot be added.")
469 } else {
470 let electric_atoms = fictitious_electric_fields.iter().flat_map(|(origin, vec)| {
471 Atom::new_special(AtomKind::Electric(true), origin + global_origin + vec, *dist_thresh).ok_or_else(||
472 format_err!("Cannot construct a fictitious electric atom.")
473 )
474 }).collect_vec();
475 mol.electric_atoms = Some(electric_atoms);
476 }
477 }
478
479 let presym = PreSymmetry::builder()
482 .moi_threshold(*moi_thresh)
483 .molecule(&mol)
484 .build()
485 .map_err(|_| format_err!("Cannot construct a pre-symmetry structure."))?;
486 let mut uni_sym = Symmetry::new();
487 let uni_res = uni_sym.analyse(&presym, false);
488 let uni_ok = uni_res.is_ok();
489 if !uni_ok {
490 log::error!("{}", uni_res.unwrap_err())
491 }
492 let uni_group_name = uni_sym.group_name.clone().unwrap_or("?".to_string());
493 let uni_group_nele = if uni_sym.is_infinite() {
494 "∞".to_string()
495 } else {
496 uni_sym.n_elements().to_string()
497 };
498
499 let (mag_sym_opt, mag_ok) = if params.time_reversal {
500 let mut mag_sym = Symmetry::new();
501 let mag_res = mag_sym.analyse(&presym, true);
502 let mag_ok = mag_res.is_ok();
503 if !mag_ok {
504 log::error!("{}", mag_res.unwrap_err())
505 }
506 (Some(mag_sym), mag_ok)
507 } else {
508 (None, true)
509 };
510 let mag_group_name = mag_sym_opt
511 .as_ref()
512 .map(|mag_sym| {
513 mag_sym
514 .group_name
515 .clone()
516 .unwrap_or("?".to_string())
517 })
518 .unwrap_or_else(|| "--".to_string());
519 let mag_group_nele = mag_sym_opt
520 .as_ref()
521 .map(|mag_sym| if mag_sym.is_infinite() {
522 "∞".to_string()
523 } else {
524 mag_sym.n_elements().to_string()
525 })
526 .unwrap_or_else(|| "--".to_string());
527
528 i += 1;
529 if uni_ok && mag_ok {
530 qsym2_output!(
531 "{:>width$} {:>12.3e} {:>12.3e} {:>14} {:>9} {:>12} {:>9}",
532 i,
533 moi_thresh,
534 dist_thresh,
535 mag_group_name,
536 mag_group_nele,
537 uni_group_name,
538 uni_group_nele,
539 width = count_length
540 );
541 Ok((presym, uni_sym, mag_sym_opt))
542 } else {
543 if !uni_ok {
544 log::debug!(
545 "Unitary group detection with MoI threshold {:.3e} and distance threshold {:.3e} has failed.",
546 moi_thresh,
547 dist_thresh
548 );
549 }
550 if !mag_ok {
551 log::debug!(
552 "Magnetic group detection with MoI threshold {:.3e} and distance threshold {:.3e} has failed.",
553 moi_thresh,
554 dist_thresh
555 );
556 }
557 qsym2_output!(
558 "{:>width$} {:>12.3e} {:>12.3e} {:>14} {:>9} {:>12} {:>9}",
559 i,
560 moi_thresh,
561 dist_thresh,
562 "--",
563 "--",
564 "--",
565 "--",
566 width = count_length
567 );
568 log::error!(
569 "Group determination with MoI threshold {:.3e} and distance threshold {:.3e} has failed.",
570 moi_thresh,
571 dist_thresh
572 );
573 bail!(
574 "Group determination with MoI threshold {:.3e} and distance threshold {:.3e} has failed.",
575 moi_thresh,
576 dist_thresh
577 )
578 }
579 })
580 .filter_map(|res_sym| res_sym.ok())
581 .collect_vec();
582 qsym2_output!("{}", "┈".repeat(count_length + 75));
583 qsym2_output!(
584 "(The number of symmetry elements is not the same as the order of the group.)"
585 );
586 qsym2_output!("");
587
588 let (highest_presym, highest_uni_sym, highest_mag_sym_opt) = syms
589 .into_iter()
590 .max_by(
591 |(presym_a, uni_sym_a, mag_sym_opt_a), (presym_b, uni_sym_b, mag_sym_opt_b)| {
592 (
593 mag_sym_opt_a
594 .as_ref()
595 .map(|mag_sym| mag_sym.is_infinite())
596 .unwrap_or(false),
597 mag_sym_opt_a
598 .as_ref()
599 .map(|mag_sym| mag_sym.n_elements())
600 .unwrap_or(0),
601 uni_sym_a.is_infinite(),
602 uni_sym_a.n_elements(),
603 1.0 / presym_a.moi_threshold,
604 1.0 / presym_a.dist_threshold,
605 )
606 .partial_cmp(&(
607 mag_sym_opt_b
608 .as_ref()
609 .map(|mag_sym| mag_sym.is_infinite())
610 .unwrap_or(false),
611 mag_sym_opt_b
612 .as_ref()
613 .map(|mag_sym| mag_sym.n_elements())
614 .unwrap_or(0),
615 uni_sym_b.is_infinite(),
616 uni_sym_b.n_elements(),
617 1.0 / presym_b.moi_threshold,
618 1.0 / presym_b.dist_threshold,
619 ))
620 .expect("Unable to perform a comparison.")
621 },
622 )
623 .ok_or_else(|| {
624 format_err!("Unable to identify the highest-symmetry group.".to_string())
625 })?;
626
627 self.result = SymmetryGroupDetectionResult::builder()
628 .parameters(params.clone())
629 .pre_symmetry(highest_presym)
630 .unitary_symmetry(highest_uni_sym)
631 .magnetic_symmetry(highest_mag_sym_opt)
632 .build()
633 .ok();
634
635 if let Some(pd_res) = self.result.as_ref() {
637 pd_res.log_output_display();
638 if let Some(name) = params.result_save_name.as_ref() {
639 write_qsym2_binary(name, QSym2FileType::Sym, pd_res)?;
640 let mut path = name.to_path_buf();
641 path.set_extension(QSym2FileType::Sym.ext());
642 qsym2_output!(
643 "Symmetry-group detection results saved as {}.",
644 path.display().to_string()
645 );
646 qsym2_output!("");
647 }
648 }
649
650 Ok(())
651 }
652}
653
654impl QSym2Driver for SymmetryGroupDetectionDriver<'_> {
655 type Params = SymmetryGroupDetectionParams;
656
657 type Outcome = SymmetryGroupDetectionResult;
658
659 fn result(&self) -> Result<&Self::Outcome, anyhow::Error> {
660 self.result
661 .as_ref()
662 .ok_or_else(|| format_err!("No symmetry-group detection results found."))
663 }
664
665 fn run(&mut self) -> Result<(), anyhow::Error> {
666 self.detect_symmetry_group()
667 }
668}
669
670fn write_element_table(f: &mut fmt::Formatter<'_>, sym: &Symmetry) -> fmt::Result {
676 let all_elements = sym
677 .elements
678 .iter()
679 .map(|(kind, kind_elements)| (false, kind, kind_elements))
680 .chain(
681 sym.generators
682 .iter()
683 .map(|(kind, kind_generators)| (true, kind, kind_generators)),
684 );
685 all_elements
686 .sorted_by_key(|(generator, kind, _)| {
687 (
688 *generator,
689 kind.contains_antiunitary(),
690 !matches!(kind, SymmetryElementKind::Proper(_)),
691 )
692 })
693 .try_for_each(|(generator, kind, kind_elements)| {
694 if !sym.is_infinite() && generator {
695 Ok::<(), fmt::Error>(())
696 } else {
697 if generator {
698 writeln!(f, "> {kind} generators")?;
699 } else {
700 writeln!(f, "> {kind} elements")?;
701 }
702 writeln!(f, "{}", "┈".repeat(54))?;
703 writeln!(
704 f,
705 "{:>7} {:>7} {:>11} {:>11} {:>11}",
706 "", "Symbol", "x", "y", "z"
707 )?;
708 writeln!(f, "{}", "┈".repeat(54))?;
709 kind_elements.keys().sorted().try_for_each(|order| {
710 let order_elements = kind_elements.get(order).unwrap_or_else(|| {
711 panic!("Elements/generators of order `{order}` cannot be retrieved.")
712 });
713 let any_element = order_elements
714 .get_index(0)
715 .expect("Unable to retrieve an element/generator of order `{order}`.");
716 let kind_str = match any_element.kind() {
717 SymmetryElementKind::Proper(_) => "",
718 SymmetryElementKind::ImproperInversionCentre(_) => " (inversion-centre)",
719 SymmetryElementKind::ImproperMirrorPlane(_) => " (mirror-plane)",
720 };
721 let au_str = match any_element.contains_antiunitary() {
722 None => "",
723 Some(AntiunitaryKind::TimeReversal) => " (time-reversed)",
724 Some(AntiunitaryKind::ComplexConjugation) => " (complex-conjugated)",
725 };
726 writeln!(f, " Order: {order}{au_str}{kind_str}")?;
727 order_elements.iter().try_for_each(|element| {
728 let axis = element.raw_axis();
729 writeln!(
730 f,
731 "{:>7} {:>7} {:>+11.7} {:>+11.7} {:>+11.7}",
732 element.get_simplified_symbol(),
733 element.get_full_symbol(),
734 axis[0],
735 axis[1],
736 axis[2]
737 )?;
738 Ok::<(), fmt::Error>(())
739 })?;
740 Ok::<(), fmt::Error>(())
741 })?;
742 writeln!(f, "{}", "┈".repeat(54))?;
743 writeln!(f)?;
744 Ok::<(), fmt::Error>(())
745 }
746 })?;
747 Ok(())
748}