qsym2/target/tensor/
axialvector_analysis.rs

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

qsym2/target/tensor/
axialvector_analysis.rs
