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::{
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
51// ==================
52// Struct definitions
53// ==================
54
55// ----------
56// Parameters
57// ----------
58
59fn 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/// Structure containing control parameters for molecule symmetrisation by bootstrapping.
76#[derive(Clone, Builder, Debug, Serialize, Deserialize)]
77pub struct MoleculeSymmetrisationBootstrapParams {
78    /// Boolean indicating if the molecule is also reoriented to align its principal axes with the
79    /// space-fixed Cartesian axes at every iteration.
80    ///
81    /// See [`Molecule::reorientate`] for more information.
82    #[builder(default = "true")]
83    #[serde(default = "default_true")]
84    pub reorientate_molecule: bool,
85
86    /// The `loose` moment-of-inertia threshold for the symmetrisation. The symmetry elements found
87    /// at this threshold level will be used to bootstrap the symmetry of the molecule.
88    #[builder(default = "1e-2")]
89    #[serde(default = "default_loose_threshold")]
90    pub loose_moi_threshold: f64,
91
92    /// The `loose` distance threshold for the symmetrisation. The symmetry elements found at this
93    /// threshold level will be used to bootstrap the symmetry of the molecule.
94    #[builder(default = "1e-2")]
95    #[serde(default = "default_loose_threshold")]
96    pub loose_distance_threshold: f64,
97
98    /// The `target` moment-of-inertia threshold for the symmetrisation.
99    #[builder(default = "1e-7")]
100    #[serde(default = "default_tight_threshold")]
101    pub target_moi_threshold: f64,
102
103    /// The `target` distance threshold for the symmetrisation.
104    #[builder(default = "1e-7")]
105    #[serde(default = "default_tight_threshold")]
106    pub target_distance_threshold: f64,
107
108    /// The maximum number of symmetrisation iterations.
109    #[builder(default = "50")]
110    #[serde(default = "default_max_iterations")]
111    pub max_iterations: usize,
112
113    /// The number of consecutive iterations during which the symmetry group at the `target` level
114    /// of threshold must be consistently found for convergence to be reached, *if this group
115    /// cannot become identical to the symmetry group at the `loose` level of threshold*.
116    #[builder(default = "10")]
117    #[serde(default = "default_consistent_iterations")]
118    pub consistent_target_symmetry_iterations: usize,
119
120    /// The finite order to which any infinite-order symmetry element is reduced, so that a finite
121    /// number of symmetry operations can be used for the symmetrisation.
122    #[builder(default = "None")]
123    #[serde(default)]
124    pub infinite_order_to_finite: Option<u32>,
125
126    /// Boolean indicating if any available magnetic group should be used for symmetrisation instead
127    /// of the unitary group.
128    #[builder(default = "false")]
129    #[serde(default)]
130    pub use_magnetic_group: bool,
131
132    /// The output verbosity level.
133    #[builder(default = "0")]
134    #[serde(default)]
135    pub verbose: u8,
136
137    /// Optional name (without the `.xyz` extension) for writing the symmetrised molecule to an XYZ
138    /// file. If `None`, no XYZ files will be written.
139    #[builder(default = "None")]
140    #[serde(default)]
141    pub symmetrised_result_xyz: Option<PathBuf>,
142
143    /// Optional name for saving the symmetry-group detection verification result of the symmetrised
144    /// system as a binary file of type [`QSym2FileType::Sym`]. If `None`, the result will not be
145    /// saved.
146    #[builder(default = "None")]
147    #[serde(default)]
148    pub symmetrised_result_save_name: Option<PathBuf>,
149}
150
151impl MoleculeSymmetrisationBootstrapParams {
152    /// Returns a builder to construct a [`MoleculeSymmetrisationBootstrapParams`] structure.
153    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// ------
233// Result
234// ------
235
236/// Structure to contain molecule symmetrisation by bootstrapping results.
237#[derive(Clone, Builder, Debug)]
238pub struct MoleculeSymmetrisationBootstrapResult<'a> {
239    /// The control parameters used to obtain this set of molecule symmetrisation results.
240    parameters: &'a MoleculeSymmetrisationBootstrapParams,
241
242    /// The symmetrised molecule.
243    pub symmetrised_molecule: Molecule,
244}
245
246impl<'a> MoleculeSymmetrisationBootstrapResult<'a> {
247    fn builder() -> MoleculeSymmetrisationBootstrapResultBuilder<'a> {
248        MoleculeSymmetrisationBootstrapResultBuilder::default()
249    }
250}
251
252// ------
253// Driver
254// ------
255
256/// Driver for molecule symmetrisation by bootstrapping in QSym².
257#[derive(Clone, Builder)]
258#[builder(build_fn(validate = "Self::validate"))]
259pub struct MoleculeSymmetrisationBootstrapDriver<'a> {
260    /// The control parameters for molecule symmetrisation by bootstrapping.
261    parameters: &'a MoleculeSymmetrisationBootstrapParams,
262
263    /// The molecule to be symmetrised.
264    molecule: &'a Molecule,
265
266    /// The result of the symmetrisation.
267    #[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    /// Returns a builder to construct a [`MoleculeSymmetrisationBootstrapDriver`] structure.
310    pub fn builder() -> MoleculeSymmetrisationBootstrapDriverBuilder<'a> {
311        MoleculeSymmetrisationBootstrapDriverBuilder::default()
312    }
313
314    /// Executes molecule symmetrisation by bootstrapping.
315    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            // If reorientation is requested, the trial molecule is reoriented prior to
338            // symmetrisation, so that the symmetrisation procedure acts on the reoriented molecule
339            // itself. The molecule might become disoriented during the symmetrisation process, but
340            // any such disorientation is likely to be fairly small, and post-symmetrisation
341            // corrections on small disorientation are better than on large disorientation.
342            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            // -------------------------------
402            // Loose threshold symmetry search
403            // -------------------------------
404            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            // This might fail, but that's fine. We are bootstrapping.
417            let _loose_res = loose_sym.analyse(&loose_presym, self.parameters.use_magnetic_group);
418
419            // Only the operations are needed for the symmetrisation. We avoid constructing the
420            // full abstract group here, as group closure might not be fulfilled due to the low
421            // thresholds.
422
423            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            // Generate transformation matrix and atom permutations for each operation
431            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            // Apply symmetry operations to the ordinary atoms
460            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                // Parallelisation here does not improve performance, and even causes more
470                // numerical instability.
471                Array2::<f64>::zeros(loose_ord_coords.raw_dim()),
472                |acc, (tmat, ord_perm, _, _)| {
473                    // coords.dot(tmat) gives the atom positions transformed in R^3 by tmat.
474                    // .select(Axis(0), perm.image()) then permutes the rows so that the atom positions
475                    // go back to approximately where they were originally.
476                    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            // Apply symmetry operations to the magnetic atoms, if any
493            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                        // coords.dot(tmat) gives the atom positions transformed in R^3 by tmat.
505                        // .select(Axis(0), perm.image()) then permutes the rows so that the atom positions
506                        // go back to approximately where they were originally.
507                        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            // Apply symmetry operations to the electric atoms, if any
530            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                        // coords.dot(tmat) gives the atom positions transformed in R^3 by tmat.
542                        // .select(Axis(0), perm.image()) then permutes the rows so that the atom positions
543                        // go back to approximately where they were originally.
544                        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            // Recentre and reorientate after symmetrisation
569            trial_mol.recentre_mut();
570            if params.reorientate_molecule {
571                // If reorientation is requested, the trial molecule is reoriented prior to
572                // symmetrisation, so that the symmetrisation procedure acts on the reoriented molecule
573                // itself. The molecule might become disoriented during the symmetrisation process, but
574                // any such disorientation is likely to be fairly small, and post-symmetrisation
575                // corrections on small disorientation are better than on large disorientation.
576                trial_mol.reorientate_mut(params.target_moi_threshold);
577            };
578
579            // -------------------------------
580            // Target threshold symmetry check
581            // -------------------------------
582            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        // --------------------------------
652        // Verify the symmetrisation result
653        // --------------------------------
654        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        // --------------
691        // Saving results
692        // --------------
693        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}