qsym2/drivers/representation_analysis/
vibrational_coordinate.rs

1//! Driver for symmetry analysis of vibrational coordinates.
2
3use std::fmt;
4use std::ops::Mul;
5
6use anyhow::{self, bail, format_err};
7use derive_builder::Builder;
8use duplicate::duplicate_item;
9use ndarray_linalg::types::Lapack;
10use num_complex::{Complex, ComplexFloat};
11use num_traits::Float;
12use rayon::prelude::*;
13use serde::{Deserialize, Serialize};
14
15use crate::analysis::{EigenvalueComparisonMode, RepAnalysis};
16use crate::chartab::chartab_group::CharacterProperties;
17use crate::chartab::SubspaceDecomposable;
18use crate::drivers::representation_analysis::angular_function::{
19    find_angular_function_representation, find_spinor_function_representation,
20    AngularFunctionRepAnalysisParams,
21};
22use crate::drivers::representation_analysis::{
23    fn_construct_magnetic_group, fn_construct_unitary_group, log_cc_transversal,
24    CharacterTableDisplay, MagneticSymmetryAnalysisKind,
25};
26use crate::drivers::symmetry_group_detection::SymmetryGroupDetectionResult;
27use crate::drivers::QSym2Driver;
28use crate::group::{GroupProperties, MagneticRepresentedGroup, UnitaryRepresentedGroup};
29use crate::io::format::{
30    log_subtitle, nice_bool, qsym2_output, write_subtitle, write_title, QSym2Output,
31};
32use crate::symmetry::symmetry_group::{
33    MagneticRepresentedSymmetryGroup, SymmetryGroupProperties, UnitaryRepresentedSymmetryGroup,
34};
35use crate::symmetry::symmetry_transformation::SymmetryTransformationKind;
36use crate::target::vibration::vibration_analysis::VibrationalCoordinateSymmetryOrbit;
37use crate::target::vibration::VibrationalCoordinateCollection;
38
39#[cfg(test)]
40#[path = "vibrational_coordinate_tests.rs"]
41mod vibrational_coordinate_tests;
42
43// ==================
44// Struct definitions
45// ==================
46
47// ----------
48// Parameters
49// ----------
50
51const fn default_true() -> bool {
52    true
53}
54const fn default_symbolic() -> Option<CharacterTableDisplay> {
55    Some(CharacterTableDisplay::Symbolic)
56}
57
58/// Structure containing control parameters for vibrational coordinate representation analysis.
59#[derive(Clone, Builder, Debug, Serialize, Deserialize)]
60pub struct VibrationalCoordinateRepAnalysisParams<T: From<f64>> {
61    /// Threshold for checking if subspace multiplicities are integral.
62    pub integrality_threshold: T,
63
64    /// Threshold for determining zero eigenvalues in the orbit overlap matrix.
65    pub linear_independence_threshold: T,
66
67    /// Option indicating if the magnetic group is to be used for symmetry analysis, and if so,
68    /// whether unitary representations or unitary-antiunitary corepresentations should be used.
69    #[builder(default = "None")]
70    #[serde(default)]
71    pub use_magnetic_group: Option<MagneticSymmetryAnalysisKind>,
72
73    /// Boolean indicating if the double group is to be used for symmetry analysis.
74    #[builder(default = "false")]
75    #[serde(default)]
76    pub use_double_group: bool,
77
78    /// Boolean indicating if the Cayley table of the group, if available, should be used to speed
79    /// up the computation of orbit overlap matrices.
80    #[builder(default = "true")]
81    #[serde(default = "default_true")]
82    pub use_cayley_table: bool,
83
84    /// The kind of symmetry transformation to be applied on the reference vibrational coordinate to
85    /// generate the orbit for symmetry analysis.
86    #[builder(default = "SymmetryTransformationKind::Spatial")]
87    #[serde(default)]
88    pub symmetry_transformation_kind: SymmetryTransformationKind,
89
90    /// Option indicating if the character table of the group used for symmetry analysis is to be
91    /// printed out.
92    #[builder(default = "Some(CharacterTableDisplay::Symbolic)")]
93    #[serde(default = "default_symbolic")]
94    pub write_character_table: Option<CharacterTableDisplay>,
95
96    /// The comparison mode for filtering out orbit overlap eigenvalues.
97    #[builder(default = "EigenvalueComparisonMode::Modulus")]
98    #[serde(default)]
99    pub eigenvalue_comparison_mode: EigenvalueComparisonMode,
100
101    /// The finite order to which any infinite-order symmetry element is reduced, so that a finite
102    /// subgroup of an infinite group can be used for the symmetry analysis.
103    #[builder(default = "None")]
104    #[serde(default)]
105    pub infinite_order_to_finite: Option<u32>,
106}
107
108impl<T> VibrationalCoordinateRepAnalysisParams<T>
109where
110    T: Float + From<f64>,
111{
112    /// Returns a builder to construct a [`VibrationalCoordinateRepAnalysisParams`] structure.
113    pub fn builder() -> VibrationalCoordinateRepAnalysisParamsBuilder<T> {
114        VibrationalCoordinateRepAnalysisParamsBuilder::default()
115    }
116}
117
118impl Default for VibrationalCoordinateRepAnalysisParams<f64> {
119    fn default() -> Self {
120        Self::builder()
121            .integrality_threshold(1e-7)
122            .linear_independence_threshold(1e-7)
123            .build()
124            .expect("Unable to construct a default `VibrationalCoordinateRepAnalysisParams<f64>`.")
125    }
126}
127
128impl<T> fmt::Display for VibrationalCoordinateRepAnalysisParams<T>
129where
130    T: From<f64> + fmt::LowerExp + fmt::Debug,
131{
132    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
133        writeln!(
134            f,
135            "Integrality threshold: {:.3e}",
136            self.integrality_threshold
137        )?;
138        writeln!(
139            f,
140            "Linear independence threshold: {:.3e}",
141            self.linear_independence_threshold
142        )?;
143        writeln!(
144            f,
145            "Orbit eigenvalue comparison mode: {}",
146            self.eigenvalue_comparison_mode
147        )?;
148        writeln!(f)?;
149        writeln!(
150            f,
151            "Use magnetic group for analysis: {}",
152            match self.use_magnetic_group {
153                None => "no",
154                Some(MagneticSymmetryAnalysisKind::Representation) =>
155                    "yes, using unitary representations",
156                Some(MagneticSymmetryAnalysisKind::Corepresentation) =>
157                    "yes, using magnetic corepresentations",
158            }
159        )?;
160        writeln!(
161            f,
162            "Use double group for analysis: {}",
163            nice_bool(self.use_double_group)
164        )?;
165        writeln!(
166            f,
167            "Use Cayley table for orbit overlap matrices: {}",
168            nice_bool(self.use_cayley_table)
169        )?;
170        if let Some(finite_order) = self.infinite_order_to_finite {
171            writeln!(f, "Infinite order to finite: {finite_order}")?;
172        }
173        writeln!(
174            f,
175            "Symmetry transformation kind: {}",
176            self.symmetry_transformation_kind
177        )?;
178        writeln!(f)?;
179        writeln!(
180            f,
181            "Write character table: {}",
182            if let Some(chartab_display) = self.write_character_table.as_ref() {
183                format!("yes, {}", chartab_display.to_string().to_lowercase())
184            } else {
185                "no".to_string()
186            }
187        )?;
188
189        Ok(())
190    }
191}
192
193// ------
194// Result
195// ------
196
197/// Structure to contain vibrational coordinate representation analysis results.
198#[derive(Clone, Builder)]
199pub struct VibrationalCoordinateRepAnalysisResult<'a, G, T>
200where
201    G: SymmetryGroupProperties + Clone,
202    G::CharTab: SubspaceDecomposable<T>,
203    T: ComplexFloat + Lapack,
204    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
205{
206    /// The control parameters used to obtain this set of vibrational coordinate representation
207    /// analysis results.
208    parameters: &'a VibrationalCoordinateRepAnalysisParams<<T as ComplexFloat>::Real>,
209
210    /// The collection of vibrational coordinates being analysed.
211    vibrational_coordinate_collection: &'a VibrationalCoordinateCollection<'a, T>,
212
213    /// The group used for the representation analysis.
214    group: G,
215
216    /// The deduced symmetries of the vibrational coordinates.
217    vib_symmetries: Vec<Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String>>,
218
219    /// The overlap eigenvalues above and below the linear independence threshold for each
220    /// vibrational coordinate symmetry deduction.
221    vib_symmetries_thresholds: Option<Vec<(Option<T>, Option<T>)>>,
222}
223
224impl<'a, G, T> VibrationalCoordinateRepAnalysisResult<'a, G, T>
225where
226    G: SymmetryGroupProperties + Clone,
227    G::CharTab: SubspaceDecomposable<T>,
228    T: ComplexFloat + Lapack,
229    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
230{
231    /// Returns a builder to construct a new [`VibrationalCoordinateRepAnalysisResultBuilder`]
232    /// structure.
233    fn builder() -> VibrationalCoordinateRepAnalysisResultBuilder<'a, G, T> {
234        VibrationalCoordinateRepAnalysisResultBuilder::default()
235    }
236}
237
238impl<'a, G, T> fmt::Display for VibrationalCoordinateRepAnalysisResult<'a, G, T>
239where
240    G: SymmetryGroupProperties + Clone,
241    G::CharTab: SubspaceDecomposable<T>,
242    T: ComplexFloat + Lapack,
243    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + fmt::Display,
244{
245    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
246        write_subtitle(f, "Orbit-based symmetry analysis results")?;
247        writeln!(f)?;
248        writeln!(
249            f,
250            "> Group: {} ({})",
251            self.group
252                .finite_subgroup_name()
253                .map(|subgroup_name| format!("{} > {}", self.group.name(), subgroup_name))
254                .unwrap_or(self.group.name()),
255            self.group.group_type().to_string().to_lowercase()
256        )?;
257        writeln!(f)?;
258
259        let vib_index_length =
260            (usize::try_from(self.vibrational_coordinate_collection.n_modes().ilog10())
261                .unwrap_or(4)
262                + 2)
263            .max(4);
264        let vib_frequency_length = self
265            .vibrational_coordinate_collection
266            .frequencies()
267            .iter()
268            .map(|freq| format!("{freq:+.3}").chars().count())
269            .max()
270            .unwrap_or(9)
271            .max(9);
272        let vib_symmetry_length = self
273            .vib_symmetries
274            .iter()
275            .map(|vib_sym_res| {
276                vib_sym_res
277                    .as_ref()
278                    .map(|sym| sym.to_string())
279                    .unwrap_or("--".to_string())
280                    .chars()
281                    .count()
282            })
283            .max()
284            .unwrap_or(8)
285            .max(8);
286
287        let vib_eig_above_length: usize = self
288            .vib_symmetries_thresholds
289            .as_ref()
290            .map(|vib_symmetries_thresholds| {
291                vib_symmetries_thresholds
292                    .iter()
293                    .map(|(above, _)| {
294                        above
295                            .as_ref()
296                            .map(|eig| format!("{eig:+.3e}"))
297                            .unwrap_or("--".to_string())
298                            .chars()
299                            .count()
300                    })
301                    .max()
302                    .unwrap_or(10)
303                    .max(10)
304            })
305            .unwrap_or(10);
306        let vib_eig_below_length: usize = self
307            .vib_symmetries_thresholds
308            .as_ref()
309            .map(|vib_symmetries_thresholds| {
310                vib_symmetries_thresholds
311                    .iter()
312                    .map(|(_, below)| {
313                        below
314                            .as_ref()
315                            .map(|eig| format!("{eig:+.3e}"))
316                            .unwrap_or("--".to_string())
317                            .chars()
318                            .count()
319                    })
320                    .max()
321                    .unwrap_or(10)
322                    .max(10)
323            })
324            .unwrap_or(10);
325
326        writeln!(f, "> Vibrational coordinate results")?;
327
328        writeln!(
329            f,
330            "{}",
331            "┈".repeat(
332                10 + vib_index_length
333                    + vib_frequency_length
334                    + vib_symmetry_length
335                    + vib_eig_above_length
336                    + vib_eig_below_length
337            )
338        )?;
339        writeln!(
340            f,
341            " {:>vib_index_length$}  {:<vib_frequency_length$}  {:<vib_symmetry_length$}  {:<vib_eig_above_length$}  Eig. below",
342            "Mode", "Frequency", "Symmetry", "Eig. above"
343        )?;
344        writeln!(
345            f,
346            "{}",
347            "┈".repeat(
348                10 + vib_index_length
349                    + vib_frequency_length
350                    + vib_symmetry_length
351                    + vib_eig_above_length
352                    + vib_eig_below_length
353            )
354        )?;
355
356        for (vibi, vib_sym) in self.vib_symmetries.iter().enumerate() {
357            let vib_frequency_str = self
358                .vibrational_coordinate_collection
359                .frequencies()
360                .get(vibi)
361                .map(|freq| format!("{freq:>+vib_frequency_length$.3}"))
362                .unwrap_or("--".to_string());
363            let vib_sym_str = vib_sym
364                .as_ref()
365                .map(|sym| sym.to_string())
366                .unwrap_or("--".to_string());
367            let (eig_above_str, eig_below_str) = self
368                .vib_symmetries_thresholds
369                .as_ref()
370                .map(|vib_symmetries_thresholds| {
371                    vib_symmetries_thresholds
372                        .get(vibi)
373                        .map(|(eig_above_opt, eig_below_opt)| {
374                            (
375                                eig_above_opt
376                                    .map(|eig_above| format!("{eig_above:>+.3e}"))
377                                    .unwrap_or("--".to_string()),
378                                eig_below_opt
379                                    .map(|eig_below| format!("{eig_below:>+.3e}"))
380                                    .unwrap_or("--".to_string()),
381                            )
382                        })
383                        .unwrap_or(("--".to_string(), "--".to_string()))
384                })
385                .unwrap_or(("--".to_string(), "--".to_string()));
386            writeln!(
387                f,
388                " {vibi:>vib_index_length$}  \
389                {vib_frequency_str:<vib_frequency_length$}  \
390                {vib_sym_str:<vib_symmetry_length$}  \
391                {eig_above_str:<vib_eig_above_length$}  \
392                {eig_below_str}"
393            )?;
394        }
395        writeln!(
396            f,
397            "{}",
398            "┈".repeat(
399                10 + vib_index_length
400                    + vib_frequency_length
401                    + vib_symmetry_length
402                    + vib_eig_above_length
403                    + vib_eig_below_length
404            )
405        )?;
406        writeln!(f)?;
407
408        Ok(())
409    }
410}
411
412impl<'a, G, T> fmt::Debug for VibrationalCoordinateRepAnalysisResult<'a, G, T>
413where
414    G: SymmetryGroupProperties + Clone,
415    G::CharTab: SubspaceDecomposable<T>,
416    T: ComplexFloat + Lapack,
417    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + fmt::Display,
418{
419    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
420        writeln!(f, "{self}")
421    }
422}
423
424impl<'a, G, T> VibrationalCoordinateRepAnalysisResult<'a, G, T>
425where
426    G: SymmetryGroupProperties + Clone,
427    G::CharTab: SubspaceDecomposable<T>,
428    T: ComplexFloat + Lapack,
429    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + fmt::Display,
430{
431    /// Returns the vibrational coordinate symmetries obtained from the analysis result.
432    pub fn vibrational_coordinate_symmetries(
433        &self,
434    ) -> &Vec<Result<<G::CharTab as SubspaceDecomposable<T>>::Decomposition, String>> {
435        &self.vib_symmetries
436    }
437}
438
439// ------
440// Driver
441// ------
442
443// ~~~~~~~~~~~~~~~~~
444// Struct definition
445// ~~~~~~~~~~~~~~~~~
446
447/// Driver structure for performing representation analysis on vibrational coordinates.
448#[derive(Clone, Builder)]
449#[builder(build_fn(validate = "Self::validate"))]
450pub struct VibrationalCoordinateRepAnalysisDriver<'a, G, T>
451where
452    G: SymmetryGroupProperties + Clone,
453    G::CharTab: SubspaceDecomposable<T>,
454    T: ComplexFloat + Lapack,
455    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
456{
457    /// The control parameters for vibrational coordinate representation analysis.
458    parameters: &'a VibrationalCoordinateRepAnalysisParams<<T as ComplexFloat>::Real>,
459
460    /// The collection of vibrational coordinates to be analysed.
461    vibrational_coordinate_collection: &'a VibrationalCoordinateCollection<'a, T>,
462
463    /// The result from symmetry-group detection on the underlying molecular structure of the
464    /// vibrational coordinates.
465    symmetry_group: &'a SymmetryGroupDetectionResult,
466
467    /// The control parameters for symmetry analysis of angular functions.
468    angular_function_parameters: &'a AngularFunctionRepAnalysisParams,
469
470    /// The result of the vibrational coordinate representation analysis.
471    #[builder(setter(skip), default = "None")]
472    result: Option<VibrationalCoordinateRepAnalysisResult<'a, G, T>>,
473}
474
475impl<'a, G, T> VibrationalCoordinateRepAnalysisDriverBuilder<'a, G, T>
476where
477    G: SymmetryGroupProperties + Clone,
478    G::CharTab: SubspaceDecomposable<T>,
479    T: ComplexFloat + Lapack,
480    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
481{
482    fn validate(&self) -> Result<(), String> {
483        let params = self.parameters.ok_or(
484            "No vibrational coordinate representation analysis parameters found.".to_string(),
485        )?;
486
487        let sym_res = self
488            .symmetry_group
489            .ok_or("No symmetry group information found.".to_string())?;
490
491        let sym = if params.use_magnetic_group.is_some() {
492            sym_res
493                .magnetic_symmetry
494                .as_ref()
495                .ok_or("Magnetic symmetry requested for representation analysis, but no magnetic symmetry found.")?
496        } else {
497            &sym_res.unitary_symmetry
498        };
499
500        if sym.is_infinite() && params.infinite_order_to_finite.is_none() {
501            Err(
502                format!(
503                    "Representation analysis cannot be performed using the entirety of the infinite group `{}`. \
504                    Consider setting the parameter `infinite_order_to_finite` to restrict to a finite subgroup instead.",
505                    sym.group_name.as_ref().expect("No symmetry group name found.")
506                )
507            )
508        } else {
509            Ok(())
510        }
511    }
512}
513
514// ~~~~~~~~~~~~~~~~~~~~~~
515// Struct implementations
516// ~~~~~~~~~~~~~~~~~~~~~~
517
518// Generic for all symmetry groups G and determinant numeric type T
519// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
520
521impl<'a, G, T> VibrationalCoordinateRepAnalysisDriver<'a, G, T>
522where
523    G: SymmetryGroupProperties + Clone,
524    G::CharTab: SubspaceDecomposable<T>,
525    T: ComplexFloat + Lapack,
526    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
527{
528    /// Returns a builder to construct a [`VibrationalCoordinateRepAnalysisDriver`] structure.
529    pub fn builder() -> VibrationalCoordinateRepAnalysisDriverBuilder<'a, G, T> {
530        VibrationalCoordinateRepAnalysisDriverBuilder::default()
531    }
532}
533
534// Specific for unitary-represented symmetry groups, but generic for determinant numeric type T
535// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
536
537impl<'a, T> VibrationalCoordinateRepAnalysisDriver<'a, UnitaryRepresentedSymmetryGroup, T>
538where
539    T: ComplexFloat + Lapack + Sync + Send,
540    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug + Sync + Send,
541    for<'b> Complex<f64>: Mul<&'b T, Output = Complex<f64>>,
542{
543    fn_construct_unitary_group!(
544        /// Constructs the unitary-represented group (which itself can be unitary or magnetic) ready
545        /// for vibrational coordinate representation analysis.
546        construct_unitary_group
547    );
548}
549
550// Specific for magnetic-represented symmetry groups, but generic for determinant numeric type T
551// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
552
553impl<'a, T> VibrationalCoordinateRepAnalysisDriver<'a, MagneticRepresentedSymmetryGroup, T>
554where
555    T: ComplexFloat + Lapack + Sync + Send,
556    <T as ComplexFloat>::Real: From<f64> + Sync + Send + fmt::LowerExp + fmt::Debug,
557    for<'b> Complex<f64>: Mul<&'b T, Output = Complex<f64>>,
558{
559    fn_construct_magnetic_group!(
560        /// Constructs the magnetic-represented group (which itself can only be magnetic) ready for
561        /// vibrational coordinate corepresentation analysis.
562        construct_magnetic_group
563    );
564}
565
566// Specific for unitary-represented and magnetic-represented symmetry groups and determinant numeric types f64 and C128
567// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
568
569#[duplicate_item(
570    duplicate!{
571        [ dtype_nested; [f64]; [Complex<f64>] ]
572        [
573            gtype_ [ UnitaryRepresentedSymmetryGroup ]
574            dtype_ [ dtype_nested ]
575            doc_sub_ [ "Performs representation analysis using a unitary-represented group and stores the result." ]
576            analyse_fn_ [ analyse_representation ]
577            construct_group_ [ self.construct_unitary_group()? ]
578        ]
579    }
580    duplicate!{
581        [ dtype_nested; [f64]; [Complex<f64>] ]
582        [
583            gtype_ [ MagneticRepresentedSymmetryGroup ]
584            dtype_ [ dtype_nested ]
585            doc_sub_ [ "Performs corepresentation analysis using a magnetic-represented group and stores the result." ]
586            analyse_fn_ [ analyse_corepresentation ]
587            construct_group_ [ self.construct_magnetic_group()? ]
588        ]
589    }
590)]
591impl<'a> VibrationalCoordinateRepAnalysisDriver<'a, gtype_, dtype_> {
592    #[doc = doc_sub_]
593    fn analyse_fn_(&mut self) -> Result<(), anyhow::Error> {
594        let params = self.parameters;
595        let group = construct_group_;
596        log_cc_transversal(&group);
597        let _ = find_angular_function_representation(&group, self.angular_function_parameters);
598        if group.is_double_group() {
599            let _ = find_spinor_function_representation(&group, self.angular_function_parameters);
600        }
601
602        let (vib_symmetries, vib_symmetries_thresholds) = {
603            let vibs = self
604                .vibrational_coordinate_collection
605                .to_vibrational_coordinates();
606            let vib_orbits = vibs.par_iter().map(|vib| {
607                let mut vib_orbit = VibrationalCoordinateSymmetryOrbit::builder()
608                    .group(&group)
609                    .origin(vib)
610                    .integrality_threshold(params.integrality_threshold)
611                    .linear_independence_threshold(params.linear_independence_threshold)
612                    .symmetry_transformation_kind(params.symmetry_transformation_kind.clone())
613                    .eigenvalue_comparison_mode(params.eigenvalue_comparison_mode.clone())
614                    .build()?;
615                vib_orbit
616                    .calc_smat(None, None, params.use_cayley_table)
617                    .and_then(|orbit| orbit.calc_xmat(false))?;
618                Ok::<_, anyhow::Error>(vib_orbit)
619            });
620            let (vib_symmetries, vib_symmetries_thresholds): (Vec<_>, Vec<_>) = vib_orbits
621                .map(|vib_orbit_res| {
622                    vib_orbit_res
623                        .as_ref()
624                        .map(|vib_orbit| {
625                            let sym_res = vib_orbit.analyse_rep().map_err(|err| err.to_string());
626                            let eigs = vib_orbit
627                                .smat_eigvals
628                                .as_ref()
629                                .map(|eigvals| {
630                                    let mut eigvals_vec = eigvals.iter().collect::<Vec<_>>();
631                                    match vib_orbit.eigenvalue_comparison_mode {
632                                        EigenvalueComparisonMode::Modulus => {
633                                            eigvals_vec.sort_by(|a, b| {
634                                                a.abs().partial_cmp(&b.abs()).expect("Unable to compare two eigenvalues based on their moduli.")
635                                            });
636                                        }
637                                        EigenvalueComparisonMode::Real => {
638                                            eigvals_vec.sort_by(|a, b| {
639                                                a.re().partial_cmp(&b.re()).expect("Unable to compare two eigenvalues based on their real parts.")
640                                            });
641                                        }
642                                    }
643                                    let eigval_above = match vib_orbit.eigenvalue_comparison_mode {
644                                        EigenvalueComparisonMode::Modulus => eigvals_vec
645                                            .iter()
646                                            .find(|val| {
647                                                val.abs() >= vib_orbit.linear_independence_threshold
648                                            })
649                                            .copied()
650                                            .copied(),
651                                        EigenvalueComparisonMode::Real => eigvals_vec
652                                            .iter()
653                                            .find(|val| {
654                                                val.re() >= vib_orbit.linear_independence_threshold
655                                            })
656                                            .copied()
657                                            .copied(),
658                                    };
659                                    eigvals_vec.reverse();
660                                    let eigval_below = match vib_orbit.eigenvalue_comparison_mode {
661                                        EigenvalueComparisonMode::Modulus => eigvals_vec
662                                            .iter()
663                                            .find(|val| {
664                                                val.abs() < vib_orbit.linear_independence_threshold
665                                            })
666                                            .copied()
667                                            .copied(),
668                                        EigenvalueComparisonMode::Real => eigvals_vec
669                                            .iter()
670                                            .find(|val| {
671                                                val.re() < vib_orbit.linear_independence_threshold
672                                            })
673                                            .copied()
674                                            .copied(),
675                                    };
676                                    (eigval_above, eigval_below)
677                                })
678                                .unwrap_or_else(|| (None, None));
679                            (sym_res, eigs)
680                        })
681                        .unwrap_or((Err("--".to_string()), (None, None)))
682                })
683                .unzip();
684            (vib_symmetries, vib_symmetries_thresholds)
685        };
686
687        let result = VibrationalCoordinateRepAnalysisResult::builder()
688            .parameters(params)
689            .vibrational_coordinate_collection(self.vibrational_coordinate_collection)
690            .group(group)
691            .vib_symmetries(vib_symmetries)
692            .vib_symmetries_thresholds(Some(vib_symmetries_thresholds))
693            .build()?;
694        self.result = Some(result);
695
696        Ok(())
697    }
698}
699
700// ~~~~~~~~~~~~~~~~~~~~~
701// Trait implementations
702// ~~~~~~~~~~~~~~~~~~~~~
703
704// Generic for all symmetry groups G and determinant numeric type T
705// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
706
707impl<'a, G, T> fmt::Display for VibrationalCoordinateRepAnalysisDriver<'a, G, T>
708where
709    G: SymmetryGroupProperties + Clone,
710    G::CharTab: SubspaceDecomposable<T>,
711    T: ComplexFloat + Lapack,
712    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
713{
714    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
715        write_title(f, "Vibrational Coordinate Symmetry Analysis")?;
716        writeln!(f)?;
717        writeln!(f, "{}", self.parameters)?;
718        Ok(())
719    }
720}
721
722impl<'a, G, T> fmt::Debug for VibrationalCoordinateRepAnalysisDriver<'a, G, T>
723where
724    G: SymmetryGroupProperties + Clone,
725    G::CharTab: SubspaceDecomposable<T>,
726    T: ComplexFloat + Lapack,
727    <T as ComplexFloat>::Real: From<f64> + fmt::LowerExp + fmt::Debug,
728{
729    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
730        writeln!(f, "{self}")
731    }
732}
733
734// Specific for unitary/magnetic-represented groups and vibration numeric type f64/Complex<f64>
735// ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
736
737#[duplicate_item(
738    duplicate!{
739        [ dtype_nested; [f64]; [Complex<f64>] ]
740        [
741            gtype_ [ UnitaryRepresentedSymmetryGroup ]
742            dtype_ [ dtype_nested ]
743            analyse_fn_ [ analyse_representation ]
744        ]
745    }
746    duplicate!{
747        [ dtype_nested; [f64]; [Complex<f64>] ]
748        [
749            gtype_ [ MagneticRepresentedSymmetryGroup ]
750            dtype_ [ dtype_nested ]
751            analyse_fn_ [ analyse_corepresentation ]
752        ]
753    }
754)]
755impl<'a> QSym2Driver for VibrationalCoordinateRepAnalysisDriver<'a, gtype_, dtype_> {
756    type Params = VibrationalCoordinateRepAnalysisParams<f64>;
757
758    type Outcome = VibrationalCoordinateRepAnalysisResult<'a, gtype_, dtype_>;
759
760    fn result(&self) -> Result<&Self::Outcome, anyhow::Error> {
761        self.result
762            .as_ref()
763            .ok_or_else(|| format_err!("No vibrational representation analysis results found."))
764    }
765
766    fn run(&mut self) -> Result<(), anyhow::Error> {
767        self.log_output_display();
768        self.analyse_fn_()?;
769        self.result()?.log_output_display();
770        Ok(())
771    }
772}