qsym2/target/vibration/
vibration_analysis.rs

1//! Implementation of symmetry analysis for vibrational coordinates.
2
3use std::fmt;
4use std::ops::Mul;
5
6use anyhow::{self, ensure, format_err, Context};
7use approx;
8use derive_builder::Builder;
9use itertools::Itertools;
10use ndarray::{Array1, Array2, Axis, Ix2};
11use ndarray_linalg::{
12    eig::Eig,
13    eigh::Eigh,
14    types::{Lapack, Scalar},
15    UPLO,
16};
17use num_complex::{Complex, ComplexFloat};
18use num_traits::{Float, Zero};
19
20use crate::analysis::{
21    fn_calc_xmat_complex, fn_calc_xmat_real, EigenvalueComparisonMode, Orbit, OrbitIterator,
22    Overlap, RepAnalysis,
23};
24use crate::auxiliary::misc::complex_modified_gram_schmidt;
25use crate::chartab::SubspaceDecomposable;
26use crate::symmetry::symmetry_element::symmetry_operation::SpecialSymmetryTransformation;
27use crate::symmetry::symmetry_group::SymmetryGroupProperties;
28use crate::symmetry::symmetry_transformation::{SymmetryTransformable, SymmetryTransformationKind};
29use crate::target::vibration::VibrationalCoordinate;
30
31// -------
32// Overlap
33// -------
34
35impl<'a, T> Overlap<T, Ix2> for VibrationalCoordinate<'a, T>
36where
37    T: Lapack
38        + ComplexFloat<Real = <T as Scalar>::Real>
39        + fmt::Debug
40        + Mul<<T as ComplexFloat>::Real, Output = T>,
41    <T as ComplexFloat>::Real: fmt::Debug
42        + approx::RelativeEq<<T as ComplexFloat>::Real>
43        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
44{
45    fn complex_symmetric(&self) -> bool {
46        false
47    }
48
49    /// Computes the overlap between two vibrational coordinates.
50    ///
51    /// Typically, `metric` is `None`, which means that the local orthogonal Cartesian coordinates
52    /// are non-overlapping.
53    ///
54    /// # Panics
55    ///
56    /// Panics if `self` and `other` have mismatched coefficient array lengths.
57    fn overlap(
58        &self,
59        other: &Self,
60        metric: Option<&Array2<T>>,
61        _: Option<&Array2<T>>,
62    ) -> Result<T, anyhow::Error> {
63        ensure!(
64            self.coefficients.len() == other.coefficients.len(),
65            "Inconsistent numbers of coefficient matrices between `self` and `other`."
66        );
67
68        let ov = if let Some(s) = metric {
69            self.coefficients
70                .t()
71                .mapv(|x| x.conj())
72                .dot(s)
73                .dot(&other.coefficients)
74        } else {
75            self.coefficients
76                .t()
77                .mapv(|x| x.conj())
78                .dot(&other.coefficients)
79        };
80        Ok(ov)
81    }
82
83    /// Returns the mathematical definition of the overlap between two vibrational coordinates.
84    fn overlap_definition(&self) -> String {
85        let k = if self.complex_symmetric() { "κ " } else { "" };
86        format!("⟨{k}v_1|v_2⟩ = [{k}v_1]† g v_2    where g is an optional metric")
87    }
88}
89
90// ==================================
91// VibrationalCoordinateSymmetryOrbit
92// ==================================
93
94// -----------------
95// Struct definition
96// -----------------
97
98/// Structure to manage symmetry orbits (*i.e.* orbits generated by symmetry groups) of
99/// vibrational coordinates.
100#[derive(Builder, Clone)]
101pub struct VibrationalCoordinateSymmetryOrbit<'a, G, T>
102where
103    G: SymmetryGroupProperties,
104    T: ComplexFloat + fmt::Debug + Lapack,
105    VibrationalCoordinate<'a, T>: SymmetryTransformable,
106{
107    /// The generating symmetry group.
108    group: &'a G,
109
110    /// The origin vibrational coordinate of the orbit.
111    origin: &'a VibrationalCoordinate<'a, T>,
112
113    /// The threshold for determining if calculated multiplicities in representation analysis are
114    /// integral.
115    integrality_threshold: <T as ComplexFloat>::Real,
116
117    /// The threshold for determining zero eigenvalues in the orbit overlap matrix.
118    pub(crate) linear_independence_threshold: <T as ComplexFloat>::Real,
119
120    /// The kind of transformation determining the way the symmetry operations in `group` act on
121    /// [`Self::origin`].
122    symmetry_transformation_kind: SymmetryTransformationKind,
123
124    /// The overlap matrix between the symmetry-equivalent vibrational coordinates in the orbit.
125    #[builder(setter(skip), default = "None")]
126    smat: Option<Array2<T>>,
127
128    /// The eigenvalues of the overlap matrix between the symmetry-equivalent vibrational
129    /// coordinates in the orbit.
130    #[builder(setter(skip), default = "None")]
131    pub(crate) smat_eigvals: Option<Array1<T>>,
132
133    /// The $`\mathbf{X}`$ matrix for the overlap matrix between the symmetry-equivalent
134    /// vibrational coordinates in the orbit.
135    ///
136    /// See [`RepAnalysis::xmat`] for further information.
137    #[builder(setter(skip), default = "None")]
138    xmat: Option<Array2<T>>,
139
140    /// An enumerated type specifying the comparison mode for filtering out orbit overlap
141    /// eigenvalues.
142    pub(crate) eigenvalue_comparison_mode: EigenvalueComparisonMode,
143}
144
145// ----------------------------
146// Struct method implementation
147// ----------------------------
148
149impl<'a, G, T> VibrationalCoordinateSymmetryOrbit<'a, G, T>
150where
151    G: SymmetryGroupProperties + Clone,
152    T: ComplexFloat + fmt::Debug + Lapack,
153    VibrationalCoordinate<'a, T>: SymmetryTransformable,
154{
155    /// Returns a builder to construct a new [`VibrationalCoordinateSymmetryOrbit`].
156    pub fn builder() -> VibrationalCoordinateSymmetryOrbitBuilder<'a, G, T> {
157        VibrationalCoordinateSymmetryOrbitBuilder::default()
158    }
159}
160
161impl<'a, G> VibrationalCoordinateSymmetryOrbit<'a, G, f64>
162where
163    G: SymmetryGroupProperties,
164{
165    fn_calc_xmat_real!(
166        /// Calculates the $`\mathbf{X}`$ matrix for real and symmetric overlap matrix
167        /// $`\mathbf{S}`$ between the symmetry-equivalent vibrational coordinates in the orbit.
168        ///
169        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
170        ///
171        /// # Arguments
172        ///
173        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
174        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit molecular
175        /// orbitals. If `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is
176        /// already of full rank.
177        pub calc_xmat
178    );
179}
180
181impl<'a, G, T> VibrationalCoordinateSymmetryOrbit<'a, G, Complex<T>>
182where
183    G: SymmetryGroupProperties,
184    T: Float + Scalar<Complex = Complex<T>>,
185    Complex<T>: ComplexFloat<Real = T> + Scalar<Real = T, Complex = Complex<T>> + Lapack,
186    VibrationalCoordinate<'a, Complex<T>>: SymmetryTransformable + Overlap<Complex<T>, Ix2>,
187{
188    fn_calc_xmat_complex!(
189        /// Calculates the $`\mathbf{X}`$ matrix for complex and symmetric or Hermitian overlap
190        /// matrix $`\mathbf{S}`$ between the symmetry-equivalent vibrational coordinates in the
191        /// orbit.
192        ///
193        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
194        ///
195        /// # Arguments
196        ///
197        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
198        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit molecular
199        /// orbitals. If `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is
200        /// already of full rank.
201        pub calc_xmat
202    );
203}
204
205// ---------------------
206// Trait implementations
207// ---------------------
208
209// ~~~~~
210// Orbit
211// ~~~~~
212
213impl<'a, G, T> Orbit<G, VibrationalCoordinate<'a, T>>
214    for VibrationalCoordinateSymmetryOrbit<'a, G, T>
215where
216    G: SymmetryGroupProperties,
217    T: ComplexFloat + fmt::Debug + Lapack,
218    VibrationalCoordinate<'a, T>: SymmetryTransformable,
219{
220    type OrbitIter = OrbitIterator<'a, G, VibrationalCoordinate<'a, T>>;
221
222    fn group(&self) -> &G {
223        self.group
224    }
225
226    fn origin(&self) -> &VibrationalCoordinate<'a, T> {
227        self.origin
228    }
229
230    fn iter(&self) -> Self::OrbitIter {
231        OrbitIterator::new(
232            self.group,
233            self.origin,
234            match self.symmetry_transformation_kind {
235                SymmetryTransformationKind::Spatial
236                | SymmetryTransformationKind::SpatialWithSpinTimeReversal => |op, vib| {
237                    // Vibrational coordinates are time-even, so both `sym_transform_spatial` and
238                    // `sym_transform_spatial_with_spintimerev` would give the same thing.
239                    vib.sym_transform_spatial(op).with_context(|| {
240                        format_err!(
241                            "Unable to apply `{op}` spatially on the origin vibrational coordinate"
242                        )
243                    })
244                },
245                SymmetryTransformationKind::Spin => |op, vib| {
246                    vib.sym_transform_spin(op).with_context(|| {
247                        format_err!(
248                            "Unable to apply `{op}` spin-wise on the origin vibrational coordinate"
249                        )
250                    })
251                },
252                SymmetryTransformationKind::SpinSpatial => |op, vib| {
253                    vib.sym_transform_spin_spatial(op).with_context(|| {
254                        format_err!(
255                            "Unable to apply `{op}` spin-spatially on the origin vibrational coordinate"
256                        )
257                    })
258                },
259            },
260        )
261    }
262}
263
264// ~~~~~~~~~~~
265// RepAnalysis
266// ~~~~~~~~~~~
267
268impl<'a, G, T> RepAnalysis<G, VibrationalCoordinate<'a, T>, T, Ix2>
269    for VibrationalCoordinateSymmetryOrbit<'a, G, T>
270where
271    G: SymmetryGroupProperties,
272    G::CharTab: SubspaceDecomposable<T>,
273    T: Lapack
274        + ComplexFloat<Real = <T as Scalar>::Real>
275        + fmt::Debug
276        + Mul<<T as ComplexFloat>::Real, Output = T>,
277    <T as ComplexFloat>::Real: fmt::Debug
278        + Zero
279        + approx::RelativeEq<<T as ComplexFloat>::Real>
280        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
281    VibrationalCoordinate<'a, T>: SymmetryTransformable,
282{
283    fn set_smat(&mut self, smat: Array2<T>) {
284        self.smat = Some(smat)
285    }
286
287    fn smat(&self) -> Option<&Array2<T>> {
288        self.smat.as_ref()
289    }
290
291    fn xmat(&self) -> &Array2<T> {
292        self.xmat
293            .as_ref()
294            .expect("Orbit overlap orthogonalisation matrix not found.")
295    }
296
297    fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error> {
298        if self.origin.complex_symmetric() {
299            Err(format_err!("`norm_preserving_scalar_map` is currently not implemented for complex symmetric overlaps."))
300        } else {
301            if self
302                .group
303                .get_index(i)
304                .unwrap_or_else(|| panic!("Group operation index `{i}` not found."))
305                .contains_time_reversal()
306            {
307                Ok(ComplexFloat::conj)
308            } else {
309                Ok(|x| x)
310            }
311        }
312    }
313
314    fn integrality_threshold(&self) -> <T as ComplexFloat>::Real {
315        self.integrality_threshold
316    }
317
318    fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode {
319        &self.eigenvalue_comparison_mode
320    }
321}
qsym2/target/vibration/vibration_analysis.rs

qsym2/target/vibration/
vibration_analysis.rs
