1use std::fmt;
23use std::path::PathBuf;
24
25use anyhow::{ensure, format_err};
26use derive_builder::Builder;
27use itertools::Itertools;
28use nalgebra::Point3;
29use ndarray::{Array2, Axis};
30use num_traits::ToPrimitive;
31use rayon::prelude::*;
32use serde::{Deserialize, Serialize};
33
34use crate::auxiliary::geometry::Transform;
35use crate::auxiliary::molecule::Molecule;
36use crate::drivers::symmetry_group_detection::{
37 SymmetryGroupDetectionDriver, SymmetryGroupDetectionParams,
38};
39use crate::drivers::QSym2Driver;
40use crate::io::format::{
41 log_subtitle, log_title, nice_bool, qsym2_output, QSym2Output,
42};
43use crate::io::QSym2FileType;
44use crate::permutation::IntoPermutation;
45use crate::symmetry::symmetry_core::{PreSymmetry, Symmetry};
46
47#[cfg(test)]
48#[path = "molecule_symmetrisation_bootstrap_tests.rs"]
49mod molecule_symmetrisation_bootstrap_tests;
50
51fn default_true() -> bool {
60 true
61}
62fn default_max_iterations() -> usize {
63 50
64}
65fn default_consistent_iterations() -> usize {
66 10
67}
68fn default_loose_threshold() -> f64 {
69 1e-2
70}
71fn default_tight_threshold() -> f64 {
72 1e-7
73}
74
75#[derive(Clone, Builder, Debug, Serialize, Deserialize)]
77pub struct MoleculeSymmetrisationBootstrapParams {
78 #[builder(default = "true")]
83 #[serde(default = "default_true")]
84 pub reorientate_molecule: bool,
85
86 #[builder(default = "1e-2")]
89 #[serde(default = "default_loose_threshold")]
90 pub loose_moi_threshold: f64,
91
92 #[builder(default = "1e-2")]
95 #[serde(default = "default_loose_threshold")]
96 pub loose_distance_threshold: f64,
97
98 #[builder(default = "1e-7")]
100 #[serde(default = "default_tight_threshold")]
101 pub target_moi_threshold: f64,
102
103 #[builder(default = "1e-7")]
105 #[serde(default = "default_tight_threshold")]
106 pub target_distance_threshold: f64,
107
108 #[builder(default = "50")]
110 #[serde(default = "default_max_iterations")]
111 pub max_iterations: usize,
112
113 #[builder(default = "10")]
117 #[serde(default = "default_consistent_iterations")]
118 pub consistent_target_symmetry_iterations: usize,
119
120 #[builder(default = "None")]
123 #[serde(default)]
124 pub infinite_order_to_finite: Option<u32>,
125
126 #[builder(default = "false")]
129 #[serde(default)]
130 pub use_magnetic_group: bool,
131
132 #[builder(default = "0")]
134 #[serde(default)]
135 pub verbose: u8,
136
137 #[builder(default = "None")]
140 #[serde(default)]
141 pub symmetrised_result_xyz: Option<PathBuf>,
142
143 #[builder(default = "None")]
147 #[serde(default)]
148 pub symmetrised_result_save_name: Option<PathBuf>,
149}
150
151impl MoleculeSymmetrisationBootstrapParams {
152 pub fn builder() -> MoleculeSymmetrisationBootstrapParamsBuilder {
154 MoleculeSymmetrisationBootstrapParamsBuilder::default()
155 }
156}
157
158impl Default for MoleculeSymmetrisationBootstrapParams {
159 fn default() -> Self {
160 Self::builder()
161 .build()
162 .expect("Unable to construct a default `MoleculeSymmetrisationBootstrapParams`.")
163 }
164}
165
166impl fmt::Display for MoleculeSymmetrisationBootstrapParams {
167 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
168 writeln!(f, "Loose MoI threshold: {:.3e}", self.loose_moi_threshold)?;
169 writeln!(
170 f,
171 "Loose geo threshold: {:.3e}",
172 self.loose_distance_threshold
173 )?;
174 writeln!(f, "Target MoI threshold: {:.3e}", self.target_moi_threshold)?;
175 writeln!(
176 f,
177 "Target geo threshold: {:.3e}",
178 self.target_distance_threshold
179 )?;
180 writeln!(f)?;
181 writeln!(
182 f,
183 "Group used for symmetrisation: {}",
184 if self.use_magnetic_group {
185 "magnetic group"
186 } else {
187 "unitary group"
188 }
189 )?;
190 if let Some(finite_order) = self.infinite_order_to_finite {
191 writeln!(f, "Infinite order to finite: {finite_order}")?;
192 }
193 writeln!(
194 f,
195 "Maximum symmetrisation iterations: {}",
196 self.max_iterations
197 )?;
198 writeln!(
199 f,
200 "Target symmetry consistent iterations: {}",
201 self.consistent_target_symmetry_iterations
202 )?;
203 writeln!(f, "Output level: {}", self.verbose)?;
204 writeln!(
205 f,
206 "Save symmetrised molecule to XYZ file: {}",
207 if let Some(name) = self.symmetrised_result_xyz.as_ref() {
208 let mut path = name.clone();
209 path.set_extension("xyz");
210 path.display().to_string()
211 } else {
212 nice_bool(false)
213 }
214 )?;
215 writeln!(
216 f,
217 "Save symmetry-group detection results of symmetrised system to file: {}",
218 if let Some(name) = self.symmetrised_result_save_name.as_ref() {
219 let mut path = name.clone();
220 path.set_extension(QSym2FileType::Sym.ext());
221 path.display().to_string()
222 } else {
223 nice_bool(false)
224 }
225 )?;
226 writeln!(f)?;
227
228 Ok(())
229 }
230}
231
232#[derive(Clone, Builder, Debug)]
238pub struct MoleculeSymmetrisationBootstrapResult<'a> {
239 parameters: &'a MoleculeSymmetrisationBootstrapParams,
241
242 pub symmetrised_molecule: Molecule,
244}
245
246impl<'a> MoleculeSymmetrisationBootstrapResult<'a> {
247 fn builder() -> MoleculeSymmetrisationBootstrapResultBuilder<'a> {
248 MoleculeSymmetrisationBootstrapResultBuilder::default()
249 }
250}
251
252#[derive(Clone, Builder)]
258#[builder(build_fn(validate = "Self::validate"))]
259pub struct MoleculeSymmetrisationBootstrapDriver<'a> {
260 parameters: &'a MoleculeSymmetrisationBootstrapParams,
262
263 molecule: &'a Molecule,
265
266 #[builder(setter(skip), default = "None")]
268 result: Option<MoleculeSymmetrisationBootstrapResult<'a>>,
269}
270
271impl<'a> MoleculeSymmetrisationBootstrapDriverBuilder<'a> {
272 fn validate(&self) -> Result<(), String> {
273 let params = self
274 .parameters
275 .ok_or("No molecule symmetrisation parameters found.".to_string())?;
276 if params.consistent_target_symmetry_iterations > params.max_iterations {
277 return Err(format!(
278 "The number of consistent target-symmetry iterations, `{}`, cannot exceed the \
279 maximum number of iterations, `{}`.",
280 params.consistent_target_symmetry_iterations, params.max_iterations,
281 ));
282 }
283 if params.target_moi_threshold < 0.0
284 || params.loose_moi_threshold < 0.0
285 || params.target_distance_threshold < 0.0
286 || params.loose_distance_threshold < 0.0
287 {
288 return Err("The thresholds cannot be negative.".to_string());
289 }
290 if params.target_moi_threshold > params.loose_moi_threshold {
291 return Err(format!(
292 "The target MoI threshold, `{:.3e}`, cannot be larger than the \
293 loose MoI threshold, `{:.3e}`.",
294 params.target_moi_threshold, params.loose_moi_threshold
295 ));
296 }
297 if params.target_distance_threshold > params.loose_distance_threshold {
298 return Err(format!(
299 "The target distance threshold, `{:.3e}`, cannot be larger than the \
300 loose distance threshold, `{:.3e}`.",
301 params.target_distance_threshold, params.loose_distance_threshold
302 ));
303 }
304 Ok(())
305 }
306}
307
308impl<'a> MoleculeSymmetrisationBootstrapDriver<'a> {
309 pub fn builder() -> MoleculeSymmetrisationBootstrapDriverBuilder<'a> {
311 MoleculeSymmetrisationBootstrapDriverBuilder::default()
312 }
313
314 fn symmetrise_molecule(&mut self) -> Result<(), anyhow::Error> {
316 log_title("Molecule Symmetrisation by Bootstrapping");
317 qsym2_output!("");
318 let params = self.parameters;
319 params.log_output_display();
320
321 let mut trial_mol = self.molecule.recentre();
322
323 if params.verbose >= 1 {
324 let orig_mol = self
325 .molecule
326 .adjust_threshold(params.target_distance_threshold);
327 qsym2_output!("Unsymmetrised original molecule:");
328 orig_mol.log_output_display();
329 qsym2_output!("");
330
331 qsym2_output!("Unsymmetrised recentred molecule:");
332 trial_mol.log_output_display();
333 qsym2_output!("");
334 }
335
336 if params.reorientate_molecule {
337 trial_mol.reorientate_mut(params.target_moi_threshold);
343 qsym2_output!("Unsymmetrised recentred and reoriented molecule:");
344 trial_mol.log_output_display();
345 qsym2_output!("");
346 };
347
348 log_subtitle("Iterative molecule symmetry bootstrapping");
349 qsym2_output!("");
350 qsym2_output!("Thresholds:");
351 qsym2_output!(
352 " Loose : {:.3e} (MoI) - {:.3e} (distance)",
353 params.loose_distance_threshold,
354 params.loose_moi_threshold,
355 );
356 qsym2_output!(
357 " Target: {:.3e} (MoI) - {:.3e} (distance)",
358 params.target_moi_threshold,
359 params.target_distance_threshold
360 );
361 qsym2_output!("");
362 qsym2_output!("Convergence criteria:");
363 qsym2_output!(
364 " either: (1) when the loose-threshold symmetry agrees with the target-threshold symmetry,",
365 );
366 qsym2_output!(
367 " or : (2) when the target-threshold symmetry contains more elements than the loose-threshold symmetry and has been consistently identified for {} consecutive iteration{}.",
368 params.consistent_target_symmetry_iterations,
369 if params.consistent_target_symmetry_iterations == 1 { "" } else { "s" }
370 );
371 qsym2_output!("");
372
373 let count_length = usize::try_from(params.max_iterations.ilog10() + 1).map_err(|_| {
374 format_err!(
375 "Unable to convert `{}` to `usize`.",
376 params.max_iterations.ilog10() + 1
377 )
378 })?;
379 qsym2_output!("{}", "┈".repeat(count_length + 101));
380 qsym2_output!(
381 " {:>count_length$} {:>22} {:>19} {:>22} {:>19} {:>10}",
382 "#",
383 "Rot. sym. (loose)",
384 "Group (loose)",
385 "Rot. sym. (target)",
386 "Group (target)",
387 "Converged?",
388 );
389 qsym2_output!("{}", "┈".repeat(count_length + 101));
390
391 let mut symmetrisation_count = 0;
392 let mut consistent_target_sym_count = 0;
393 let mut loose_ops = vec![];
394 let mut prev_target_sym_group_name: Option<String> = None;
395 let mut converged = false;
396 while symmetrisation_count == 0
397 || (!converged && symmetrisation_count < params.max_iterations)
398 {
399 symmetrisation_count += 1;
400
401 let mut loose_mol =
405 trial_mol.adjust_threshold(self.parameters.loose_distance_threshold);
406 let loose_presym = PreSymmetry::builder()
407 .moi_threshold(self.parameters.loose_moi_threshold)
408 .molecule(&loose_mol)
409 .build()
410 .map_err(|_| {
411 format_err!("Cannot construct a loose-threshold pre-symmetry structure.")
412 })?;
413
414 let mut loose_sym = Symmetry::new();
415
416 let _loose_res = loose_sym.analyse(&loose_presym, self.parameters.use_magnetic_group);
418
419 loose_ops.extend_from_slice(
424 &loose_sym.generate_all_operations(self.parameters.infinite_order_to_finite),
425 );
426 let n_ops_f64 = loose_ops.len().to_f64().ok_or_else(|| {
427 format_err!("Unable to convert the number of operations to `f64`.")
428 })?;
429
430 let ts = loose_ops
432 .into_par_iter()
433 .flat_map(|op| {
434 let tmat = op
435 .get_3d_spatial_matrix()
436 .select(Axis(0), &[2, 0, 1])
437 .select(Axis(1), &[2, 0, 1])
438 .reversed_axes();
439
440 let ord_perm = op
441 .act_permute(&loose_mol.molecule_ordinary_atoms())
442 .ok_or_else(|| {
443 format_err!(
444 "Unable to determine the ordinary-atom permutation corresponding to `{op}`."
445 )
446 })?;
447 let mag_perm_opt = loose_mol
448 .molecule_magnetic_atoms()
449 .as_ref()
450 .and_then(|loose_mag_mol| op.act_permute(loose_mag_mol));
451 let elec_perm_opt = loose_mol
452 .molecule_electric_atoms()
453 .as_ref()
454 .and_then(|loose_elec_mol| op.act_permute(loose_elec_mol));
455 Ok::<_, anyhow::Error>((tmat, ord_perm, mag_perm_opt, elec_perm_opt))
456 })
457 .collect::<Vec<_>>();
458
459 let loose_ord_coords = Array2::from_shape_vec(
461 (loose_mol.atoms.len(), 3),
462 loose_mol
463 .atoms
464 .iter()
465 .flat_map(|atom| atom.coordinates.coords.iter().cloned())
466 .collect::<Vec<_>>(),
467 )?;
468 let ave_ord_coords = ts.iter().fold(
469 Array2::<f64>::zeros(loose_ord_coords.raw_dim()),
472 |acc, (tmat, ord_perm, _, _)| {
473 acc + loose_ord_coords.dot(tmat).select(Axis(0), ord_perm.image())
477 },
478 ) / n_ops_f64;
479 loose_mol
480 .atoms
481 .par_iter_mut()
482 .enumerate()
483 .for_each(|(i, atom)| {
484 atom.coordinates = Point3::<f64>::from_slice(
485 ave_ord_coords
486 .row(i)
487 .as_slice()
488 .expect("Unable to convert a row of averaged coordinates to a slice."),
489 )
490 });
491
492 if let Some(mag_atoms) = loose_mol.magnetic_atoms.as_mut() {
494 let loose_mag_coords = Array2::from_shape_vec(
495 (mag_atoms.len(), 3),
496 mag_atoms
497 .iter()
498 .flat_map(|atom| atom.coordinates.coords.iter().cloned())
499 .collect::<Vec<_>>(),
500 )?;
501 let ave_mag_coords = ts.iter().fold(
502 Ok(Array2::<f64>::zeros(loose_mag_coords.raw_dim())),
503 |acc: Result<Array2<f64>, anyhow::Error>, (tmat, _, mag_perm_opt, _)| {
504 Ok(acc?
508 + loose_mag_coords.dot(tmat).select(
509 Axis(0),
510 mag_perm_opt
511 .as_ref()
512 .ok_or_else(|| {
513 format_err!("Expected magnetic atom permutation not found.")
514 })?
515 .image(),
516 ))
517 },
518 )? / n_ops_f64;
519 mag_atoms.iter_mut().enumerate().for_each(|(i, atom)| {
520 atom.coordinates = Point3::<f64>::from_slice(
521 ave_mag_coords
522 .row(i)
523 .as_slice()
524 .expect("Unable to convert a row of averaged coordinates to a slice."),
525 )
526 });
527 }
528
529 if let Some(elec_atoms) = loose_mol.electric_atoms.as_mut() {
531 let loose_elec_coords = Array2::from_shape_vec(
532 (elec_atoms.len(), 3),
533 elec_atoms
534 .iter()
535 .flat_map(|atom| atom.coordinates.coords.iter().cloned())
536 .collect::<Vec<_>>(),
537 )?;
538 let ave_elec_coords = ts.iter().fold(
539 Ok(Array2::<f64>::zeros(loose_elec_coords.raw_dim())),
540 |acc: Result<Array2<f64>, anyhow::Error>, (tmat, _, _, elec_perm_opt)| {
541 Ok(acc?
545 + loose_elec_coords.dot(tmat).select(
546 Axis(0),
547 elec_perm_opt
548 .as_ref()
549 .ok_or_else(|| {
550 format_err!("Expected electric atom permutation not found.")
551 })?
552 .image(),
553 ))
554 },
555 )? / n_ops_f64;
556 elec_atoms.iter_mut().enumerate().for_each(|(i, atom)| {
557 atom.coordinates = Point3::<f64>::from_slice(
558 ave_elec_coords
559 .row(i)
560 .as_slice()
561 .expect("Unable to convert a row of averaged coordinates to a slice."),
562 )
563 });
564 }
565
566 trial_mol = loose_mol;
567
568 trial_mol.recentre_mut();
570 if params.reorientate_molecule {
571 trial_mol.reorientate_mut(params.target_moi_threshold);
577 };
578
579 let target_mol = trial_mol.adjust_threshold(self.parameters.target_distance_threshold);
583 let target_presym = PreSymmetry::builder()
584 .moi_threshold(self.parameters.target_moi_threshold)
585 .molecule(&target_mol)
586 .build()
587 .map_err(|_| {
588 format_err!("Cannot construct a target-threshold pre-symmetry structure.")
589 })?;
590
591 let mut target_sym = Symmetry::new();
592
593 let _ = target_sym.analyse(&target_presym, params.use_magnetic_group);
594
595 let target_loose_consistent = target_sym.n_elements() == loose_sym.n_elements()
596 && target_sym.group_name.is_some()
597 && target_sym.group_name == loose_sym.group_name;
598
599 if target_sym.group_name == prev_target_sym_group_name
600 && target_sym.n_elements() >= loose_sym.n_elements()
601 {
602 consistent_target_sym_count += 1;
603 } else {
604 consistent_target_sym_count = 0;
605 }
606 prev_target_sym_group_name = target_sym.group_name.clone();
607 let target_consistent =
608 consistent_target_sym_count >= params.consistent_target_symmetry_iterations;
609
610 converged = target_loose_consistent || target_consistent;
611 let converged_reason = [target_loose_consistent, target_consistent]
612 .iter()
613 .enumerate()
614 .filter_map(|(i, c)| {
615 if *c {
616 Some(format!("({})", i + 1))
617 } else {
618 None
619 }
620 })
621 .join("");
622
623 qsym2_output!(
624 " {:>count_length$} {:>22} {:>19} {:>22} {:>19} {:>10}",
625 symmetrisation_count,
626 loose_presym.rotational_symmetry.to_string(),
627 format!(
628 "{} ({})",
629 loose_sym.group_name.as_ref().unwrap_or(&"--".to_string()),
630 loose_sym.n_elements()
631 ),
632 target_presym.rotational_symmetry.to_string(),
633 format!(
634 "{} ({})",
635 target_sym.group_name.as_ref().unwrap_or(&"--".to_string()),
636 target_sym.n_elements()
637 ),
638 if converged {
639 "yes ".to_string() + &converged_reason
640 } else {
641 "no".to_string()
642 },
643 );
644
645 loose_ops =
646 target_sym.generate_all_operations(self.parameters.infinite_order_to_finite);
647 }
648 qsym2_output!("{}", "┈".repeat(count_length + 101));
649 qsym2_output!("");
650
651 qsym2_output!("Verifying symmetrisation results...");
655 qsym2_output!("");
656 let verifying_pd_params = SymmetryGroupDetectionParams::builder()
657 .moi_thresholds(&[params.target_moi_threshold])
658 .distance_thresholds(&[params.target_distance_threshold])
659 .time_reversal(params.use_magnetic_group)
660 .write_symmetry_elements(true)
661 .result_save_name(params.symmetrised_result_save_name.clone())
662 .build()?;
663 let mut verifying_pd_driver = SymmetryGroupDetectionDriver::builder()
664 .parameters(&verifying_pd_params)
665 .molecule(Some(&trial_mol))
666 .build()?;
667 verifying_pd_driver.run()?;
668 let verifying_pd_res = verifying_pd_driver.result()?;
669 let verifying_group_name = if params.use_magnetic_group {
670 verifying_pd_res
671 .magnetic_symmetry
672 .as_ref()
673 .and_then(|magsym| magsym.group_name.as_ref())
674 } else {
675 verifying_pd_res
676 .unitary_symmetry
677 .group_name
678 .as_ref()
679 };
680 ensure!(
681 prev_target_sym_group_name.as_ref() == verifying_group_name,
682 "Mismatched symmetry: iterative symmetry bootstrapping found {}, but verification found {}.",
683 prev_target_sym_group_name.as_ref().unwrap_or(&"--".to_string()),
684 verifying_group_name.unwrap_or(&"--".to_string()),
685
686 );
687 qsym2_output!("Verifying symmetrisation results... Done.");
688 qsym2_output!("");
689
690 self.result = Some(
694 MoleculeSymmetrisationBootstrapResult::builder()
695 .parameters(self.parameters)
696 .symmetrised_molecule(trial_mol.clone())
697 .build()?,
698 );
699
700 if let Some(xyz_name) = params.symmetrised_result_xyz.as_ref() {
701 let mut path = xyz_name.clone();
702 path.set_extension("xyz");
703 verifying_pd_res
704 .pre_symmetry
705 .recentred_molecule
706 .to_xyz(&path)?;
707 qsym2_output!("Symmetrised molecule written to: {}", path.display());
708 qsym2_output!("");
709 }
710
711 Ok(())
712 }
713}
714
715impl<'a> QSym2Driver for MoleculeSymmetrisationBootstrapDriver<'a> {
716 type Params = MoleculeSymmetrisationBootstrapParams;
717
718 type Outcome = MoleculeSymmetrisationBootstrapResult<'a>;
719
720 fn result(&self) -> Result<&Self::Outcome, anyhow::Error> {
721 self.result
722 .as_ref()
723 .ok_or_else(|| format_err!("No molecule sprucing results found."))
724 }
725
726 fn run(&mut self) -> Result<(), anyhow::Error> {
727 self.symmetrise_molecule()
728 }
729}