qsym2/drivers/molecule_symmetrisation_bootstrap/
mod.rs

1//! Driver for molecule symmetrisation by bootstrapping in QSym².
2//!
3//! This algorithm symmetrises a molecule iteratively by defining two threshold levels: a `loose`
4//! level and a `target` level.
5//!
6//! In every iteration, the following steps are performed:
7//!
8//! 1. The molecule is symmetry-analysed at the `target` level; any symmetry elements found are stashed and the symmetry group name, if any, is registered.
9//! 2. The molecule is symmetry-analysed at the `loose` level; any symmetry elements found are added to the stash and the symmetry group name, if any, is registered.
10//! 3. The convergence criteria (see below) are checked.
11//!     - If convergence has been reached, the symmetrisation procedure is terminated.
12//!     - If convergence has not been reached, the following steps are carried out.
13//! 4. All symmetry elements found in the stash are used to generate all possible symmetry operations which are then used to symmetrise the molecule: each symmetry operation is applied on the original molecule to produce a symmetry-equivalent copy, then all symmetry-equivalent copies are averaged to give the symmetrised molecule.
14//! 5. Repeat steps 1 to 4 above until convergence is reached.
15//!
16//! There are two convergence criteria for the symmetrisation procedure:
17//! - **either** when the loose-threshold symmetry agrees with the target-threshold symmetry,
18//! - **or** when the target-threshold symmetry contains more elements than the loose-threshold symmetry and has been consistently identified for a pre-specified number of consecutive iterations.
19//!
20//! At least one criterion must be satisfied in order for convergence to be reached.
21
22use 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::{log_subtitle, log_title, nice_bool, qsym2_output, QSym2Output};
41use crate::io::QSym2FileType;
42use crate::permutation::IntoPermutation;
43use crate::symmetry::symmetry_core::{PreSymmetry, Symmetry};
44
45#[cfg(test)]
46#[path = "molecule_symmetrisation_bootstrap_tests.rs"]
47mod molecule_symmetrisation_bootstrap_tests;
48
49// ==================
50// Struct definitions
51// ==================
52
53// ----------
54// Parameters
55// ----------
56
57fn default_true() -> bool {
58    true
59}
60fn default_max_iterations() -> usize {
61    50
62}
63fn default_consistent_iterations() -> usize {
64    10
65}
66fn default_loose_threshold() -> f64 {
67    1e-2
68}
69fn default_tight_threshold() -> f64 {
70    1e-7
71}
72
73/// Structure containing control parameters for molecule symmetrisation by bootstrapping.
74#[derive(Clone, Builder, Debug, Serialize, Deserialize)]
75pub struct MoleculeSymmetrisationBootstrapParams {
76    /// Boolean indicating if the molecule is also reoriented to align its principal axes with the
77    /// space-fixed Cartesian axes at every iteration.
78    ///
79    /// See [`Molecule::reorientate`] for more information.
80    #[builder(default = "true")]
81    #[serde(default = "default_true")]
82    pub reorientate_molecule: bool,
83
84    /// The `loose` moment-of-inertia threshold for the symmetrisation. The symmetry elements found
85    /// at this threshold level will be used to bootstrap the symmetry of the molecule.
86    #[builder(default = "1e-2")]
87    #[serde(default = "default_loose_threshold")]
88    pub loose_moi_threshold: f64,
89
90    /// The `loose` distance threshold for the symmetrisation. The symmetry elements found at this
91    /// threshold level will be used to bootstrap the symmetry of the molecule.
92    #[builder(default = "1e-2")]
93    #[serde(default = "default_loose_threshold")]
94    pub loose_distance_threshold: f64,
95
96    /// The `target` moment-of-inertia threshold for the symmetrisation.
97    #[builder(default = "1e-7")]
98    #[serde(default = "default_tight_threshold")]
99    pub target_moi_threshold: f64,
100
101    /// The `target` distance threshold for the symmetrisation.
102    #[builder(default = "1e-7")]
103    #[serde(default = "default_tight_threshold")]
104    pub target_distance_threshold: f64,
105
106    /// The maximum number of symmetrisation iterations.
107    #[builder(default = "50")]
108    #[serde(default = "default_max_iterations")]
109    pub max_iterations: usize,
110
111    /// The number of consecutive iterations during which the symmetry group at the `target` level
112    /// of threshold must be consistently found for convergence to be reached, *if this group
113    /// cannot become identical to the symmetry group at the `loose` level of threshold*.
114    #[builder(default = "10")]
115    #[serde(default = "default_consistent_iterations")]
116    pub consistent_target_symmetry_iterations: usize,
117
118    /// The finite order to which any infinite-order symmetry element is reduced, so that a finite
119    /// number of symmetry operations can be used for the symmetrisation.
120    #[builder(default = "None")]
121    #[serde(default)]
122    pub infinite_order_to_finite: Option<u32>,
123
124    /// Boolean indicating if any available magnetic group should be used for symmetrisation instead
125    /// of the unitary group.
126    #[builder(default = "false")]
127    #[serde(default)]
128    pub use_magnetic_group: bool,
129
130    /// The output verbosity level.
131    #[builder(default = "0")]
132    #[serde(default)]
133    pub verbose: u8,
134
135    /// Optional name (without the `.xyz` extension) for writing the symmetrised molecule to an XYZ
136    /// file. If `None`, no XYZ files will be written.
137    #[builder(default = "None")]
138    #[serde(default)]
139    pub symmetrised_result_xyz: Option<PathBuf>,
140
141    /// Optional name for saving the symmetry-group detection verification result of the symmetrised
142    /// system as a binary file of type [`QSym2FileType::Sym`]. If `None`, the result will not be
143    /// saved.
144    #[builder(default = "None")]
145    #[serde(default)]
146    pub symmetrised_result_save_name: Option<PathBuf>,
147}
148
149impl MoleculeSymmetrisationBootstrapParams {
150    /// Returns a builder to construct a [`MoleculeSymmetrisationBootstrapParams`] structure.
151    pub fn builder() -> MoleculeSymmetrisationBootstrapParamsBuilder {
152        MoleculeSymmetrisationBootstrapParamsBuilder::default()
153    }
154}
155
156impl Default for MoleculeSymmetrisationBootstrapParams {
157    fn default() -> Self {
158        Self::builder()
159            .build()
160            .expect("Unable to construct a default `MoleculeSymmetrisationBootstrapParams`.")
161    }
162}
163
164impl fmt::Display for MoleculeSymmetrisationBootstrapParams {
165    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
166        writeln!(f, "Loose MoI threshold: {:.3e}", self.loose_moi_threshold)?;
167        writeln!(
168            f,
169            "Loose geo threshold: {:.3e}",
170            self.loose_distance_threshold
171        )?;
172        writeln!(f, "Target MoI threshold: {:.3e}", self.target_moi_threshold)?;
173        writeln!(
174            f,
175            "Target geo threshold: {:.3e}",
176            self.target_distance_threshold
177        )?;
178        writeln!(f)?;
179        writeln!(
180            f,
181            "Group used for symmetrisation: {}",
182            if self.use_magnetic_group {
183                "magnetic group"
184            } else {
185                "unitary group"
186            }
187        )?;
188        if let Some(finite_order) = self.infinite_order_to_finite {
189            writeln!(f, "Infinite order to finite: {finite_order}")?;
190        }
191        writeln!(
192            f,
193            "Maximum symmetrisation iterations: {}",
194            self.max_iterations
195        )?;
196        writeln!(
197            f,
198            "Target symmetry consistent iterations: {}",
199            self.consistent_target_symmetry_iterations
200        )?;
201        writeln!(f, "Output level: {}", self.verbose)?;
202        writeln!(
203            f,
204            "Save symmetrised molecule to XYZ file: {}",
205            if let Some(name) = self.symmetrised_result_xyz.as_ref() {
206                let mut path = name.clone();
207                path.set_extension("xyz");
208                path.display().to_string()
209            } else {
210                nice_bool(false)
211            }
212        )?;
213        writeln!(
214            f,
215            "Save symmetry-group detection results of symmetrised system to file: {}",
216            if let Some(name) = self.symmetrised_result_save_name.as_ref() {
217                let mut path = name.clone();
218                path.set_extension(QSym2FileType::Sym.ext());
219                path.display().to_string()
220            } else {
221                nice_bool(false)
222            }
223        )?;
224        writeln!(f)?;
225
226        Ok(())
227    }
228}
229
230// ------
231// Result
232// ------
233
234/// Structure to contain molecule symmetrisation by bootstrapping results.
235#[derive(Clone, Builder, Debug)]
236pub struct MoleculeSymmetrisationBootstrapResult<'a> {
237    /// The control parameters used to obtain this set of molecule symmetrisation results.
238    parameters: &'a MoleculeSymmetrisationBootstrapParams,
239
240    /// The symmetrised molecule.
241    pub symmetrised_molecule: Molecule,
242}
243
244impl<'a> MoleculeSymmetrisationBootstrapResult<'a> {
245    fn builder() -> MoleculeSymmetrisationBootstrapResultBuilder<'a> {
246        MoleculeSymmetrisationBootstrapResultBuilder::default()
247    }
248}
249
250// ------
251// Driver
252// ------
253
254/// Driver for molecule symmetrisation by bootstrapping in QSym².
255#[derive(Clone, Builder)]
256#[builder(build_fn(validate = "Self::validate"))]
257pub struct MoleculeSymmetrisationBootstrapDriver<'a> {
258    /// The control parameters for molecule symmetrisation by bootstrapping.
259    parameters: &'a MoleculeSymmetrisationBootstrapParams,
260
261    /// The molecule to be symmetrised.
262    molecule: &'a Molecule,
263
264    /// The result of the symmetrisation.
265    #[builder(setter(skip), default = "None")]
266    result: Option<MoleculeSymmetrisationBootstrapResult<'a>>,
267}
268
269impl<'a> MoleculeSymmetrisationBootstrapDriverBuilder<'a> {
270    fn validate(&self) -> Result<(), String> {
271        let params = self
272            .parameters
273            .ok_or("No molecule symmetrisation parameters found.".to_string())?;
274        if params.consistent_target_symmetry_iterations > params.max_iterations {
275            return Err(format!(
276                "The number of consistent target-symmetry iterations, `{}`, cannot exceed the \
277                    maximum number of iterations, `{}`.",
278                params.consistent_target_symmetry_iterations, params.max_iterations,
279            ));
280        }
281        if params.target_moi_threshold < 0.0
282            || params.loose_moi_threshold < 0.0
283            || params.target_distance_threshold < 0.0
284            || params.loose_distance_threshold < 0.0
285        {
286            return Err("The thresholds cannot be negative.".to_string());
287        }
288        if params.target_moi_threshold > params.loose_moi_threshold {
289            return Err(format!(
290                "The target MoI threshold, `{:.3e}`, cannot be larger than the \
291                    loose MoI threshold, `{:.3e}`.",
292                params.target_moi_threshold, params.loose_moi_threshold
293            ));
294        }
295        if params.target_distance_threshold > params.loose_distance_threshold {
296            return Err(format!(
297                "The target distance threshold, `{:.3e}`, cannot be larger than the \
298                    loose distance threshold, `{:.3e}`.",
299                params.target_distance_threshold, params.loose_distance_threshold
300            ));
301        }
302        Ok(())
303    }
304}
305
306impl<'a> MoleculeSymmetrisationBootstrapDriver<'a> {
307    /// Returns a builder to construct a [`MoleculeSymmetrisationBootstrapDriver`] structure.
308    pub fn builder() -> MoleculeSymmetrisationBootstrapDriverBuilder<'a> {
309        MoleculeSymmetrisationBootstrapDriverBuilder::default()
310    }
311
312    /// Executes molecule symmetrisation by bootstrapping.
313    fn symmetrise_molecule(&mut self) -> Result<(), anyhow::Error> {
314        log_title("Molecule Symmetrisation by Bootstrapping");
315        qsym2_output!("");
316        let params = self.parameters;
317        params.log_output_display();
318
319        let mut trial_mol = self.molecule.recentre();
320
321        if params.verbose >= 1 {
322            let orig_mol = self
323                .molecule
324                .adjust_threshold(params.target_distance_threshold);
325            qsym2_output!("Unsymmetrised original molecule:");
326            orig_mol.log_output_display();
327            qsym2_output!("");
328
329            qsym2_output!("Unsymmetrised recentred molecule:");
330            trial_mol.log_output_display();
331            qsym2_output!("");
332        }
333
334        if params.reorientate_molecule {
335            // If reorientation is requested, the trial molecule is reoriented prior to
336            // symmetrisation, so that the symmetrisation procedure acts on the reoriented molecule
337            // itself. The molecule might become disoriented during the symmetrisation process, but
338            // any such disorientation is likely to be fairly small, and post-symmetrisation
339            // corrections on small disorientation are better than on large disorientation.
340            trial_mol.reorientate_mut(params.target_moi_threshold);
341            qsym2_output!("Unsymmetrised recentred and reoriented molecule:");
342            trial_mol.log_output_display();
343            qsym2_output!("");
344        };
345
346        log_subtitle("Iterative molecule symmetry bootstrapping");
347        qsym2_output!("");
348        qsym2_output!("Thresholds:");
349        qsym2_output!(
350            "  Loose : {:.3e} (MoI) - {:.3e} (distance)",
351            params.loose_distance_threshold,
352            params.loose_moi_threshold,
353        );
354        qsym2_output!(
355            "  Target: {:.3e} (MoI) - {:.3e} (distance)",
356            params.target_moi_threshold,
357            params.target_distance_threshold
358        );
359        qsym2_output!("");
360        qsym2_output!("Convergence criteria:");
361        qsym2_output!(
362            "  either: (1) when the loose-threshold symmetry agrees with the target-threshold symmetry,",
363        );
364        qsym2_output!(
365            "  or    : (2) when the target-threshold symmetry contains more elements than the loose-threshold symmetry and has been consistently identified for {} consecutive iteration{}.",
366            params.consistent_target_symmetry_iterations,
367            if params.consistent_target_symmetry_iterations == 1 { "" } else { "s" }
368        );
369        qsym2_output!("");
370
371        let count_length = usize::try_from(params.max_iterations.ilog10() + 1).map_err(|_| {
372            format_err!(
373                "Unable to convert `{}` to `usize`.",
374                params.max_iterations.ilog10() + 1
375            )
376        })?;
377        qsym2_output!("{}", "┈".repeat(count_length + 101));
378        qsym2_output!(
379            " {:>count_length$} {:>22} {:>19}  {:>22} {:>19}  {:>10}",
380            "#",
381            "Rot. sym. (loose)",
382            "Group (loose)",
383            "Rot. sym. (target)",
384            "Group (target)",
385            "Converged?",
386        );
387        qsym2_output!("{}", "┈".repeat(count_length + 101));
388
389        let mut symmetrisation_count = 0;
390        let mut consistent_target_sym_count = 0;
391        let mut loose_ops = vec![];
392        let mut prev_target_sym_group_name: Option<String> = None;
393        let mut converged = false;
394        while symmetrisation_count == 0
395            || (!converged && symmetrisation_count < params.max_iterations)
396        {
397            symmetrisation_count += 1;
398
399            // -------------------------------
400            // Loose threshold symmetry search
401            // -------------------------------
402            let mut loose_mol =
403                trial_mol.adjust_threshold(self.parameters.loose_distance_threshold);
404            let loose_presym = PreSymmetry::builder()
405                .moi_threshold(self.parameters.loose_moi_threshold)
406                .molecule(&loose_mol)
407                .build()
408                .map_err(|_| {
409                    format_err!("Cannot construct a loose-threshold pre-symmetry structure.")
410                })?;
411
412            let mut loose_sym = Symmetry::new();
413
414            // This might fail, but that's fine. We are bootstrapping.
415            let _loose_res = loose_sym.analyse(&loose_presym, self.parameters.use_magnetic_group);
416
417            // Only the operations are needed for the symmetrisation. We avoid constructing the
418            // full abstract group here, as group closure might not be fulfilled due to the low
419            // thresholds.
420
421            loose_ops.extend_from_slice(
422                &loose_sym.generate_all_operations(self.parameters.infinite_order_to_finite),
423            );
424            let n_ops_f64 = loose_ops.len().to_f64().ok_or_else(|| {
425                format_err!("Unable to convert the number of operations to `f64`.")
426            })?;
427
428            // Generate transformation matrix and atom permutations for each operation
429            let ts = loose_ops
430                .into_par_iter()
431                .flat_map(|op| {
432                    let tmat = op
433                        .get_3d_spatial_matrix()
434                        .select(Axis(0), &[2, 0, 1])
435                        .select(Axis(1), &[2, 0, 1])
436                        .reversed_axes();
437
438                    let ord_perm = op
439                        .act_permute(&loose_mol.molecule_ordinary_atoms())
440                        .ok_or_else(|| {
441                            format_err!(
442                                "Unable to determine the ordinary-atom permutation corresponding to `{op}`."
443                            )
444                        })?;
445                    let mag_perm_opt = loose_mol
446                        .molecule_magnetic_atoms()
447                        .as_ref()
448                        .and_then(|loose_mag_mol| op.act_permute(loose_mag_mol));
449                    let elec_perm_opt = loose_mol
450                        .molecule_electric_atoms()
451                        .as_ref()
452                        .and_then(|loose_elec_mol| op.act_permute(loose_elec_mol));
453                    Ok::<_, anyhow::Error>((tmat, ord_perm, mag_perm_opt, elec_perm_opt))
454                })
455                .collect::<Vec<_>>();
456
457            // Apply symmetry operations to the ordinary atoms
458            let loose_ord_coords = Array2::from_shape_vec(
459                (loose_mol.atoms.len(), 3),
460                loose_mol
461                    .atoms
462                    .iter()
463                    .flat_map(|atom| atom.coordinates.coords.iter().cloned())
464                    .collect::<Vec<_>>(),
465            )?;
466            let ave_ord_coords = ts.iter().fold(
467                // Parallelisation here does not improve performance, and even causes more
468                // numerical instability.
469                Array2::<f64>::zeros(loose_ord_coords.raw_dim()),
470                |acc, (tmat, ord_perm, _, _)| {
471                    // coords.dot(tmat) gives the atom positions transformed in R^3 by tmat.
472                    // .select(Axis(0), perm.image()) then permutes the rows so that the atom positions
473                    // go back to approximately where they were originally.
474                    acc + loose_ord_coords.dot(tmat).select(Axis(0), ord_perm.image())
475                },
476            ) / n_ops_f64;
477            loose_mol
478                .atoms
479                .par_iter_mut()
480                .enumerate()
481                .for_each(|(i, atom)| {
482                    atom.coordinates = Point3::<f64>::from_slice(
483                        ave_ord_coords
484                            .row(i)
485                            .as_slice()
486                            .expect("Unable to convert a row of averaged coordinates to a slice."),
487                    )
488                });
489
490            // Apply symmetry operations to the magnetic atoms, if any
491            if let Some(mag_atoms) = loose_mol.magnetic_atoms.as_mut() {
492                let loose_mag_coords = Array2::from_shape_vec(
493                    (mag_atoms.len(), 3),
494                    mag_atoms
495                        .iter()
496                        .flat_map(|atom| atom.coordinates.coords.iter().cloned())
497                        .collect::<Vec<_>>(),
498                )?;
499                let ave_mag_coords = ts.iter().fold(
500                    Ok(Array2::<f64>::zeros(loose_mag_coords.raw_dim())),
501                    |acc: Result<Array2<f64>, anyhow::Error>, (tmat, _, mag_perm_opt, _)| {
502                        // coords.dot(tmat) gives the atom positions transformed in R^3 by tmat.
503                        // .select(Axis(0), perm.image()) then permutes the rows so that the atom positions
504                        // go back to approximately where they were originally.
505                        Ok(acc?
506                            + loose_mag_coords.dot(tmat).select(
507                                Axis(0),
508                                mag_perm_opt
509                                    .as_ref()
510                                    .ok_or_else(|| {
511                                        format_err!("Expected magnetic atom permutation not found.")
512                                    })?
513                                    .image(),
514                            ))
515                    },
516                )? / n_ops_f64;
517                mag_atoms.iter_mut().enumerate().for_each(|(i, atom)| {
518                    atom.coordinates = Point3::<f64>::from_slice(
519                        ave_mag_coords
520                            .row(i)
521                            .as_slice()
522                            .expect("Unable to convert a row of averaged coordinates to a slice."),
523                    )
524                });
525            }
526
527            // Apply symmetry operations to the electric atoms, if any
528            if let Some(elec_atoms) = loose_mol.electric_atoms.as_mut() {
529                let loose_elec_coords = Array2::from_shape_vec(
530                    (elec_atoms.len(), 3),
531                    elec_atoms
532                        .iter()
533                        .flat_map(|atom| atom.coordinates.coords.iter().cloned())
534                        .collect::<Vec<_>>(),
535                )?;
536                let ave_elec_coords = ts.iter().fold(
537                    Ok(Array2::<f64>::zeros(loose_elec_coords.raw_dim())),
538                    |acc: Result<Array2<f64>, anyhow::Error>, (tmat, _, _, elec_perm_opt)| {
539                        // coords.dot(tmat) gives the atom positions transformed in R^3 by tmat.
540                        // .select(Axis(0), perm.image()) then permutes the rows so that the atom positions
541                        // go back to approximately where they were originally.
542                        Ok(acc?
543                            + loose_elec_coords.dot(tmat).select(
544                                Axis(0),
545                                elec_perm_opt
546                                    .as_ref()
547                                    .ok_or_else(|| {
548                                        format_err!("Expected electric atom permutation not found.")
549                                    })?
550                                    .image(),
551                            ))
552                    },
553                )? / n_ops_f64;
554                elec_atoms.iter_mut().enumerate().for_each(|(i, atom)| {
555                    atom.coordinates = Point3::<f64>::from_slice(
556                        ave_elec_coords
557                            .row(i)
558                            .as_slice()
559                            .expect("Unable to convert a row of averaged coordinates to a slice."),
560                    )
561                });
562            }
563
564            trial_mol = loose_mol;
565
566            // Recentre and reorientate after symmetrisation
567            trial_mol.recentre_mut();
568            if params.reorientate_molecule {
569                // If reorientation is requested, the trial molecule is reoriented prior to
570                // symmetrisation, so that the symmetrisation procedure acts on the reoriented molecule
571                // itself. The molecule might become disoriented during the symmetrisation process, but
572                // any such disorientation is likely to be fairly small, and post-symmetrisation
573                // corrections on small disorientation are better than on large disorientation.
574                trial_mol.reorientate_mut(params.target_moi_threshold);
575            };
576
577            // -------------------------------
578            // Target threshold symmetry check
579            // -------------------------------
580            let target_mol = trial_mol.adjust_threshold(self.parameters.target_distance_threshold);
581            let target_presym = PreSymmetry::builder()
582                .moi_threshold(self.parameters.target_moi_threshold)
583                .molecule(&target_mol)
584                .build()
585                .map_err(|_| {
586                    format_err!("Cannot construct a target-threshold pre-symmetry structure.")
587                })?;
588
589            let mut target_sym = Symmetry::new();
590
591            let _ = target_sym.analyse(&target_presym, params.use_magnetic_group);
592
593            let target_loose_consistent = target_sym.n_elements() == loose_sym.n_elements()
594                && target_sym.group_name.is_some()
595                && target_sym.group_name == loose_sym.group_name;
596
597            if target_sym.group_name == prev_target_sym_group_name
598                && target_sym.n_elements() >= loose_sym.n_elements()
599            {
600                consistent_target_sym_count += 1;
601            } else {
602                consistent_target_sym_count = 0;
603            }
604            prev_target_sym_group_name = target_sym.group_name.clone();
605            let target_consistent =
606                consistent_target_sym_count >= params.consistent_target_symmetry_iterations;
607
608            converged = target_loose_consistent || target_consistent;
609            let converged_reason = [target_loose_consistent, target_consistent]
610                .iter()
611                .enumerate()
612                .filter_map(|(i, c)| {
613                    if *c {
614                        Some(format!("({})", i + 1))
615                    } else {
616                        None
617                    }
618                })
619                .join("");
620
621            qsym2_output!(
622                " {:>count_length$} {:>22} {:>19}  {:>22} {:>19}  {:>10}",
623                symmetrisation_count,
624                loose_presym.rotational_symmetry.to_string(),
625                format!(
626                    "{} ({})",
627                    loose_sym.group_name.as_ref().unwrap_or(&"--".to_string()),
628                    loose_sym.n_elements()
629                ),
630                target_presym.rotational_symmetry.to_string(),
631                format!(
632                    "{} ({})",
633                    target_sym.group_name.as_ref().unwrap_or(&"--".to_string()),
634                    target_sym.n_elements()
635                ),
636                if converged {
637                    "yes ".to_string() + &converged_reason
638                } else {
639                    "no".to_string()
640                },
641            );
642
643            loose_ops =
644                target_sym.generate_all_operations(self.parameters.infinite_order_to_finite);
645        }
646        qsym2_output!("{}", "┈".repeat(count_length + 101));
647        qsym2_output!("");
648
649        // --------------------------------
650        // Verify the symmetrisation result
651        // --------------------------------
652        qsym2_output!("Verifying symmetrisation results...");
653        qsym2_output!("");
654        let verifying_pd_params = SymmetryGroupDetectionParams::builder()
655            .moi_thresholds(&[params.target_moi_threshold])
656            .distance_thresholds(&[params.target_distance_threshold])
657            .time_reversal(params.use_magnetic_group)
658            .write_symmetry_elements(true)
659            .result_save_name(params.symmetrised_result_save_name.clone())
660            .build()?;
661        let mut verifying_pd_driver = SymmetryGroupDetectionDriver::builder()
662            .parameters(&verifying_pd_params)
663            .molecule(Some(&trial_mol))
664            .build()?;
665        verifying_pd_driver.run()?;
666        let verifying_pd_res = verifying_pd_driver.result()?;
667        let verifying_group_name = if params.use_magnetic_group {
668            verifying_pd_res
669                .magnetic_symmetry
670                .as_ref()
671                .and_then(|magsym| magsym.group_name.as_ref())
672        } else {
673            verifying_pd_res.unitary_symmetry.group_name.as_ref()
674        };
675        ensure!(
676            prev_target_sym_group_name.as_ref() == verifying_group_name,
677            "Mismatched symmetry: iterative symmetry bootstrapping found {}, but verification found {}.",
678            prev_target_sym_group_name.as_ref().unwrap_or(&"--".to_string()),
679            verifying_group_name.unwrap_or(&"--".to_string()),
680
681        );
682        qsym2_output!("Verifying symmetrisation results... Done.");
683        qsym2_output!("");
684
685        // --------------
686        // Saving results
687        // --------------
688        self.result = Some(
689            MoleculeSymmetrisationBootstrapResult::builder()
690                .parameters(self.parameters)
691                .symmetrised_molecule(trial_mol.clone())
692                .build()?,
693        );
694
695        if let Some(xyz_name) = params.symmetrised_result_xyz.as_ref() {
696            let mut path = xyz_name.clone();
697            path.set_extension("xyz");
698            verifying_pd_res
699                .pre_symmetry
700                .recentred_molecule
701                .to_xyz(&path)?;
702            qsym2_output!("Symmetrised molecule written to: {}", path.display());
703            qsym2_output!("");
704        }
705
706        Ok(())
707    }
708}
709
710impl<'a> QSym2Driver for MoleculeSymmetrisationBootstrapDriver<'a> {
711    type Params = MoleculeSymmetrisationBootstrapParams;
712
713    type Outcome = MoleculeSymmetrisationBootstrapResult<'a>;
714
715    fn result(&self) -> Result<&Self::Outcome, anyhow::Error> {
716        self.result
717            .as_ref()
718            .ok_or_else(|| format_err!("No molecule sprucing results found."))
719    }
720
721    fn run(&mut self) -> Result<(), anyhow::Error> {
722        self.symmetrise_molecule()
723    }
724}