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, Context, format_err};
7use approx;
8use derive_builder::Builder;
9use itertools::Itertools;
10use log;
11use ndarray::{Array1, Array2, Axis, Ix2};
12use ndarray_linalg::{
13    UPLO,
14    eig::Eig,
15    eigh::Eigh,
16    types::{Lapack, Scalar},
17};
18use num_complex::{Complex, ComplexFloat};
19use num_traits::{Float, Zero};
20
21use crate::analysis::{
22    EigenvalueComparisonMode, Orbit, OrbitIterator, Overlap, RepAnalysis, fn_calc_xmat_complex,
23    fn_calc_xmat_real,
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 threshold for checking symmetry of the orbit overlap matrix.
110    orbit_overlap_symmetry_threshold: <T as ComplexFloat>::Real,
111
112    /// The kind of transformation determining the way the symmetry operations in `group` act on
113    /// [`Self::origin`].
114    symmetry_transformation_kind: SymmetryTransformationKind,
115
116    /// The overlap matrix between the symmetry-equivalent axial vectors in the orbit.
117    #[builder(setter(skip), default = "None")]
118    smat: Option<Array2<T>>,
119
120    /// The eigenvalues of the overlap matrix between the symmetry-equivalent axial vectors in the
121    /// orbit.
122    #[builder(setter(skip), default = "None")]
123    pub(crate) smat_eigvals: Option<Array1<T>>,
124
125    /// The $`\mathbf{X}`$ matrix for the overlap matrix between the symmetry-equivalent axial
126    /// vectors in the orbit.
127    ///
128    /// See [`RepAnalysis::xmat`] for further information.
129    #[builder(setter(skip), default = "None")]
130    xmat: Option<Array2<T>>,
131
132    /// An enumerated type specifying the comparison mode for filtering out orbit overlap
133    /// eigenvalues.
134    pub(crate) eigenvalue_comparison_mode: EigenvalueComparisonMode,
135}
136
137// ----------------------------
138// Struct method implementation
139// ----------------------------
140
141impl<'a, G, T> AxialVector3SymmetryOrbit<'a, G, T>
142where
143    G: SymmetryGroupProperties + Clone,
144    T: ComplexFloat + fmt::Debug + Lapack,
145    AxialVector3<T>: SymmetryTransformable,
146{
147    /// Returns a builder to construct a new [`AxialVector3SymmetryOrbit`] structure.
148    pub fn builder() -> AxialVector3SymmetryOrbitBuilder<'a, G, T> {
149        AxialVector3SymmetryOrbitBuilder::default()
150    }
151}
152
153impl<'a, G> AxialVector3SymmetryOrbit<'a, G, f64>
154where
155    G: SymmetryGroupProperties,
156{
157    fn_calc_xmat_real!(
158        /// Calculates the $`\mathbf{X}`$ matrix for real and symmetric overlap matrix
159        /// $`\mathbf{S}`$ between the symmetry-equivalent axial vectors in the orbit.
160        ///
161        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
162        ///
163        /// # Arguments
164        ///
165        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
166        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit axial vectors.
167        /// If `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is already of
168        /// full rank.
169        pub calc_xmat
170    );
171}
172
173impl<'a, G, T> AxialVector3SymmetryOrbit<'a, G, Complex<T>>
174where
175    G: SymmetryGroupProperties,
176    T: Float + Scalar<Complex = Complex<T>>,
177    Complex<T>: ComplexFloat<Real = T> + Scalar<Real = T, Complex = Complex<T>> + Lapack,
178    AxialVector3<Complex<T>>: SymmetryTransformable + Overlap<Complex<T>, Ix2>,
179{
180    fn_calc_xmat_complex!(
181        /// Calculates the $`\mathbf{X}`$ matrix for complex and symmetric or Hermitian overlap
182        /// matrix $`\mathbf{S}`$ between the symmetry-equivalent axial vectors in the orbit.
183        ///
184        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
185        ///
186        /// # Arguments
187        ///
188        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
189        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit axial vectors.
190        /// If `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is already of
191        /// full rank.
192        pub calc_xmat
193    );
194}
195
196// ---------------------
197// Trait implementations
198// ---------------------
199
200// ~~~~~
201// Orbit
202// ~~~~~
203
204impl<'a, G, T> Orbit<G, AxialVector3<T>> for AxialVector3SymmetryOrbit<'a, G, T>
205where
206    G: SymmetryGroupProperties,
207    T: ComplexFloat + fmt::Debug + Lapack,
208    AxialVector3<T>: SymmetryTransformable,
209{
210    type OrbitIter = OrbitIterator<'a, G, AxialVector3<T>>;
211
212    fn group(&self) -> &G {
213        self.group
214    }
215
216    fn origin(&self) -> &AxialVector3<T> {
217        self.origin
218    }
219
220    fn iter(&self) -> Self::OrbitIter {
221        OrbitIterator::new(
222            self.group,
223            self.origin,
224            match self.symmetry_transformation_kind {
225                SymmetryTransformationKind::Spatial => |op, axvec| {
226                    axvec.sym_transform_spatial(op).with_context(|| {
227                        format!("Unable to apply `{op}` spatially on the origin axial vector")
228                    })
229                },
230                SymmetryTransformationKind::SpatialWithSpinTimeReversal => |op, axvec| {
231                    axvec.sym_transform_spatial_with_spintimerev(op).with_context(|| {
232                        format!("Unable to apply `{op}` spatially (with potentially direction-reversing time reversal) on the origin axial vector")
233                    })
234                },
235                SymmetryTransformationKind::Spin => |op, axvec| {
236                    axvec.sym_transform_spin(op).with_context(|| {
237                        format!("Unable to apply `{op}` spin-wise on the origin axial vector")
238                    })
239                },
240                SymmetryTransformationKind::SpinSpatial => |op, axvec| {
241                    axvec.sym_transform_spin_spatial(op).with_context(|| {
242                        format!("Unable to apply `{op}` spin-spatially on the origin axial vector")
243                    })
244                },
245            },
246        )
247    }
248}
249
250// ~~~~~~~~~~~
251// RepAnalysis
252// ~~~~~~~~~~~
253
254impl<'a, G, T> RepAnalysis<G, AxialVector3<T>, T, Ix2> for AxialVector3SymmetryOrbit<'a, G, T>
255where
256    G: SymmetryGroupProperties,
257    G::CharTab: SubspaceDecomposable<T>,
258    T: Lapack
259        + ComplexFloat<Real = <T as Scalar>::Real>
260        + fmt::Debug
261        + Mul<<T as ComplexFloat>::Real, Output = T>,
262    <T as ComplexFloat>::Real: fmt::Debug
263        + Zero
264        + approx::RelativeEq<<T as ComplexFloat>::Real>
265        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
266    AxialVector3<T>: SymmetryTransformable,
267{
268    fn set_smat(&mut self, smat: Array2<T>) {
269        self.smat = Some(smat)
270    }
271
272    fn smat(&self) -> Option<&Array2<T>> {
273        self.smat.as_ref()
274    }
275
276    fn xmat(&self) -> &Array2<T> {
277        self.xmat
278            .as_ref()
279            .expect("Orbit overlap orthogonalisation matrix not found.")
280    }
281
282    fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error> {
283        if self.origin.complex_symmetric() {
284            Err(format_err!(
285                "`norm_preserving_scalar_map` is currently not implemented for complex-symmetric overlaps. This thus precludes the use of the Cayley table to speed up the computation of the orbit overlap matrix."
286            ))
287        } else if self
288            .group
289            .get_index(i)
290            .unwrap_or_else(|| panic!("Group operation index `{i}` not found."))
291            .contains_time_reversal()
292        {
293            Ok(ComplexFloat::conj)
294        } else {
295            Ok(|x| x)
296        }
297    }
298
299    fn orbit_overlap_symmetry_threshold(&self) -> <T as ComplexFloat>::Real {
300        self.orbit_overlap_symmetry_threshold
301    }
302
303    fn integrality_threshold(&self) -> <T as ComplexFloat>::Real {
304        self.integrality_threshold
305    }
306
307    fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode {
308        &self.eigenvalue_comparison_mode
309    }
310
311    /// Reduces the representation or corepresentation spanned by the axial vectors in the
312    /// orbit to a direct sum of the irreducible representations or corepresentations of the
313    /// generating symmetry group.
314    ///
315    /// # Returns
316    ///
317    /// The decomposed result.
318    ///
319    /// # Errors
320    ///
321    /// Errors if the decomposition fails, *e.g.* because one or more calculated multiplicities
322    /// are non-integral.
323    fn analyse_rep(
324        &self,
325    ) -> Result<
326        <<G as CharacterProperties>::CharTab as SubspaceDecomposable<T>>::Decomposition,
327        DecompositionError,
328    > {
329        let chis = self
330            .calc_characters()
331            .map_err(|err| DecompositionError(err.to_string()))?;
332        self.group().character_table().reduce_characters(
333            &chis.iter().map(|(cc, chi)| (cc, *chi)).collect::<Vec<_>>(),
334            self.integrality_threshold(),
335        )
336    }
337}
qsym2/target/tensor/axialvector_analysis.rs

qsym2/target/tensor/
axialvector_analysis.rs
