qsym2/sandbox/target/real_space_function/
real_space_function_analysis.rs

1//! Implementation of symmetry analysis for real-space functions.
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 log;
11use nalgebra::Point3;
12use ndarray::{Array1, Array2, Axis, Ix1};
13use ndarray_linalg::{
14    eig::Eig,
15    eigh::Eigh,
16    types::{Lapack, Scalar},
17    UPLO,
18};
19use num_complex::{Complex, ComplexFloat};
20use num_traits::{Float, ToPrimitive, Zero};
21use rayon::prelude::*;
22
23use crate::analysis::{
24    fn_calc_xmat_complex, fn_calc_xmat_real, EigenvalueComparisonMode, Orbit, OrbitIterator,
25    Overlap, RepAnalysis,
26};
27use crate::auxiliary::misc::complex_modified_gram_schmidt;
28use crate::chartab::chartab_group::CharacterProperties;
29use crate::chartab::{DecompositionError, SubspaceDecomposable};
30use crate::io::format::{log_subtitle, qsym2_output, QSym2Output};
31use crate::sandbox::target::real_space_function::RealSpaceFunction;
32use crate::symmetry::symmetry_element::symmetry_operation::SpecialSymmetryTransformation;
33use crate::symmetry::symmetry_group::SymmetryGroupProperties;
34use crate::symmetry::symmetry_transformation::{SymmetryTransformable, SymmetryTransformationKind};
35
36// =======
37// Overlap
38// =======
39
40impl<T, F> Overlap<T, Ix1> for RealSpaceFunction<T, F>
41where
42    T: ComplexFloat + Lapack + Send + Sync,
43    F: Clone + Sync + Send + Fn(&Point3<f64>) -> T,
44{
45    fn complex_symmetric(&self) -> bool {
46        false
47    }
48
49    /// Computes the overlap between two real-space functions.
50    ///
51    /// No attempts are made to check that the grid points between the two real-space functions are
52    /// 'compatible'. The overlap is simply evaluated as
53    ///
54    /// ```math
55    ///     \sum_i [\hat{\iota} f_1(\mathbf{r}_{1i})]^* f_2(\mathbf{r}_{2i}) w_i,
56    /// ```
57    ///
58    /// where  $`w_i`$ are weight values.
59    ///
60    ///
61    /// # Errors
62    ///
63    /// Errors if `self` and `other` have mismatched numbers of grid points or if the number of
64    /// weight values does not match the number of grid points.
65    fn overlap(
66        &self,
67        other: &Self,
68        metric: Option<&Array1<T>>,
69        _: Option<&Array1<T>>,
70    ) -> Result<T, anyhow::Error> {
71        let weight = metric.ok_or_else(|| {
72            format_err!("No weights found for `RealSpaceFunction` overlap calculation.")
73        })?;
74
75        ensure!(
76            self.grid_points.len() == other.grid_points.len(),
77            "Inconsistent number of grid points between `self` and `other`."
78        );
79        ensure!(
80            self.grid_points.len() == weight.len(),
81            "Inconsistent number of weight values and grid points."
82        );
83
84        let overlap = (0..weight.len())
85            .into_par_iter()
86            .map(|i| {
87                let s_pt = self.grid_points[i];
88                let o_pt = other.grid_points[i];
89                let w = weight[i];
90                match (self.complex_conjugated, other.complex_conjugated) {
91                    (false, false) => self.function()(&s_pt).conj() * other.function()(&o_pt) * w,
92                    (false, true) => {
93                        self.function()(&s_pt).conj() * other.function()(&o_pt).conj() * w
94                    }
95                    (true, false) => self.function()(&s_pt) * other.function()(&o_pt) * w,
96                    (true, true) => self.function()(&s_pt) * other.function()(&o_pt).conj() * w,
97                }
98            })
99            .sum();
100        Ok(overlap)
101    }
102
103    /// Returns the mathematical definition of the overlap between two real-space functions.
104    fn overlap_definition(&self) -> String {
105        let k = if self.complex_symmetric() { "κ " } else { "" };
106        format!("⟨{k}f_1|f_2⟩ = ∫ [{k}f_1(r)]* w(r) f_2(r) dr    where w(r) is a required weight")
107    }
108}
109
110// ==============================
111// RealSpaceFunctionSymmetryOrbit
112// ==============================
113
114// -----------------
115// Struct definition
116// -----------------
117
118/// Structure to manage symmetry orbits (*i.e.* orbits generated by symmetry groups) of real-space
119/// functions.
120#[derive(Builder, Clone)]
121pub struct RealSpaceFunctionSymmetryOrbit<'a, G, T, F>
122where
123    G: SymmetryGroupProperties,
124    T: ComplexFloat + fmt::Debug + Lapack,
125    F: Fn(&Point3<f64>) -> T,
126    RealSpaceFunction<T, F>: SymmetryTransformable,
127{
128    /// The generating symmetry group.
129    group: &'a G,
130
131    /// The origin [`RealSpaceFunction`] of the orbit.
132    origin: &'a RealSpaceFunction<T, F>,
133
134    /// The threshold for determining zero eigenvalues in the orbit overlap matrix.
135    pub(crate) linear_independence_threshold: <T as ComplexFloat>::Real,
136
137    /// The threshold for determining if calculated multiplicities in representation analysis are
138    /// integral.
139    integrality_threshold: <T as ComplexFloat>::Real,
140
141    /// The kind of transformation determining the way the symmetry operations in the generating
142    /// group act on [`Self::origin`].
143    symmetry_transformation_kind: SymmetryTransformationKind,
144
145    /// The overlap matrix between the symmetry-equivalent real-space functions in the orbit.
146    #[builder(setter(skip), default = "None")]
147    smat: Option<Array2<T>>,
148
149    /// The eigenvalues of the overlap matrix between the symmetry-equivalent real-space functions
150    /// in the orbit.
151    #[builder(setter(skip), default = "None")]
152    pub(crate) smat_eigvals: Option<Array1<T>>,
153
154    /// The $`\mathbf{X}`$ matrix for the overlap matrix between the symmetry-equivalent real-space
155    /// functions in the orbit.
156    ///
157    /// See [`RepAnalysis::xmat`] for further information.
158    #[builder(setter(skip), default = "None")]
159    xmat: Option<Array2<T>>,
160
161    /// An enumerated type specifying the comparison mode for filtering out orbit overlap
162    /// eigenvalues.
163    eigenvalue_comparison_mode: EigenvalueComparisonMode,
164}
165
166// ----------------------------
167// Struct method implementation
168// ----------------------------
169
170impl<'a, G, T, F> RealSpaceFunctionSymmetryOrbit<'a, G, T, F>
171where
172    G: SymmetryGroupProperties + Clone,
173    T: ComplexFloat + fmt::Debug + Lapack,
174    F: Clone + Fn(&Point3<f64>) -> T,
175    RealSpaceFunction<T, F>: SymmetryTransformable,
176{
177    /// Returns a builder for constructing a new real-space function symmetry orbit.
178    pub fn builder() -> RealSpaceFunctionSymmetryOrbitBuilder<'a, G, T, F> {
179        RealSpaceFunctionSymmetryOrbitBuilder::default()
180    }
181}
182
183impl<'a, G, F> RealSpaceFunctionSymmetryOrbit<'a, G, f64, F>
184where
185    G: SymmetryGroupProperties + Clone,
186    F: Clone + Fn(&Point3<f64>) -> f64,
187{
188    fn_calc_xmat_real!(
189        /// Calculates the $`\mathbf{X}`$ matrix for real and symmetric overlap matrix
190        /// $`\mathbf{S}`$ between the symmetry-equivalent real-space functions in the orbit.
191        ///
192        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
193        ///
194        /// # Arguments
195        ///
196        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
197        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit real-space
198        /// functions. If `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is
199        /// already of full rank.
200        pub calc_xmat
201    );
202}
203
204impl<'a, G, T, F> RealSpaceFunctionSymmetryOrbit<'a, G, Complex<T>, F>
205where
206    G: SymmetryGroupProperties + Clone,
207    T: Float + Scalar<Complex = Complex<T>>,
208    Complex<T>: ComplexFloat<Real = T> + Scalar<Real = T, Complex = Complex<T>> + Lapack,
209    F: Clone + Fn(&Point3<f64>) -> Complex<T>,
210    RealSpaceFunction<Complex<T>, F>: SymmetryTransformable + Overlap<Complex<T>, Ix1>,
211{
212    fn_calc_xmat_complex!(
213        /// Calculates the $`\mathbf{X}`$ matrix for complex and symmetric or Hermitian overlap
214        /// matrix $`\mathbf{S}`$ between the symmetry-equivalent real-space functions in the orbit.
215        ///
216        /// The resulting $`\mathbf{X}`$ is stored in the orbit.
217        ///
218        /// # Arguments
219        ///
220        /// * `preserves_full_rank` - If `true`, when $`\mathbf{S}`$ is already of full rank, then
221        /// $`\mathbf{X}`$ is set to be the identity matrix to avoid mixing the orbit real-space
222        /// functions. If `false`, $`\mathbf{X}`$ also orthogonalises $`\mathbf{S}`$ even when it is
223        /// already of full rank.
224        pub calc_xmat
225    );
226}
227
228// ---------------------
229// Trait implementations
230// ---------------------
231
232// ~~~~~
233// Orbit
234// ~~~~~
235
236impl<'a, G, T, F> Orbit<G, RealSpaceFunction<T, F>> for RealSpaceFunctionSymmetryOrbit<'a, G, T, F>
237where
238    G: SymmetryGroupProperties + Clone,
239    T: ComplexFloat + fmt::Debug + Lapack,
240    F: Fn(&Point3<f64>) -> T,
241    RealSpaceFunction<T, F>: SymmetryTransformable,
242{
243    type OrbitIter = OrbitIterator<'a, G, RealSpaceFunction<T, F>>;
244
245    fn group(&self) -> &'a G {
246        self.group
247    }
248
249    fn origin(&self) -> &RealSpaceFunction<T, F> {
250        self.origin
251    }
252
253    fn iter(&self) -> Self::OrbitIter {
254        OrbitIterator::new(
255            self.group,
256            self.origin,
257            match self.symmetry_transformation_kind {
258                SymmetryTransformationKind::Spatial => |op, real_space_function| {
259                    real_space_function.sym_transform_spatial(op).with_context(|| {
260                        format!("Unable to apply `{op}` spatially on the origin real-space function")
261                    })
262                },
263                SymmetryTransformationKind::SpatialWithSpinTimeReversal => {
264                    |op, real_space_function| {
265                        real_space_function.sym_transform_spatial_with_spintimerev(op).with_context(|| {
266                        format!("Unable to apply `{op}` spatially (with spin-including time-reversal) on the origin real-space function")
267                    })
268                    }
269                }
270                SymmetryTransformationKind::Spin => {
271                    |op, real_space_function| {
272                        real_space_function.sym_transform_spin(op).with_context(|| {
273                        format!("Unable to apply `{op}` spin-wise on the origin real-space function")
274                    })
275                    }
276                }
277                SymmetryTransformationKind::SpinSpatial => |op, real_space_function| {
278                    real_space_function.sym_transform_spin_spatial(op).with_context(|| {
279                        format!("Unable to apply `{op}` spin-spatially on the origin real-space function")
280                    })
281                },
282            },
283        )
284    }
285}
286
287// ~~~~~~~~~~~
288// RepAnalysis
289// ~~~~~~~~~~~
290
291impl<'a, G, T, F> RepAnalysis<G, RealSpaceFunction<T, F>, T, Ix1>
292    for RealSpaceFunctionSymmetryOrbit<'a, G, T, F>
293where
294    G: SymmetryGroupProperties + Clone,
295    G::CharTab: SubspaceDecomposable<T>,
296    T: Lapack
297        + ComplexFloat<Real = <T as Scalar>::Real>
298        + fmt::Debug
299        + Send
300        + Sync
301        + Mul<<T as ComplexFloat>::Real, Output = T>,
302    <T as ComplexFloat>::Real: fmt::Debug
303        + Zero
304        + From<u16>
305        + ToPrimitive
306        + approx::RelativeEq<<T as ComplexFloat>::Real>
307        + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
308    F: Clone + Sync + Send + Fn(&Point3<f64>) -> T,
309    RealSpaceFunction<T, F>: SymmetryTransformable,
310{
311    fn set_smat(&mut self, smat: Array2<T>) {
312        self.smat = Some(smat)
313    }
314
315    fn smat(&self) -> Option<&Array2<T>> {
316        self.smat.as_ref()
317    }
318
319    fn xmat(&self) -> &Array2<T> {
320        self.xmat
321            .as_ref()
322            .expect("Orbit overlap orthogonalisation matrix not found.")
323    }
324
325    fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error> {
326        if self.origin.complex_symmetric() {
327            Err(format_err!("`norm_preserving_scalar_map` is currently not implemented for complex symmetric overlaps."))
328        } else {
329            if self
330                .group
331                .get_index(i)
332                .unwrap_or_else(|| panic!("Group operation index `{i}` not found."))
333                .contains_time_reversal()
334            {
335                Ok(ComplexFloat::conj)
336            } else {
337                Ok(|x| x)
338            }
339        }
340    }
341
342    fn integrality_threshold(&self) -> <T as ComplexFloat>::Real {
343        self.integrality_threshold
344    }
345
346    fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode {
347        &self.eigenvalue_comparison_mode
348    }
349
350    /// Reduces the representation or corepresentation spanned by the real-space functions in the
351    /// orbit to a direct sum of the irreducible representations or corepresentations of the
352    /// generating symmetry group.
353    ///
354    /// # Returns
355    ///
356    /// The decomposed result.
357    ///
358    /// # Errors
359    ///
360    /// Errors if the decomposition fails, *e.g.* because one or more calculated multiplicities
361    /// are non-integral.
362    fn analyse_rep(
363        &self,
364    ) -> Result<
365        <<G as CharacterProperties>::CharTab as SubspaceDecomposable<T>>::Decomposition,
366        DecompositionError,
367    > {
368        log::debug!("Analysing representation symmetry for a real-space function...");
369        let chis = self
370            .calc_characters()
371            .map_err(|err| DecompositionError(err.to_string()))?;
372        log::debug!("Characters calculated.");
373        log_subtitle("Real-space function orbit characters");
374        qsym2_output!("");
375        self.characters_to_string(&chis, self.integrality_threshold)
376            .log_output_display();
377        qsym2_output!("");
378
379        let res = self.group().character_table().reduce_characters(
380            &chis.iter().map(|(cc, chi)| (cc, *chi)).collect::<Vec<_>>(),
381            self.integrality_threshold(),
382        );
383        log::debug!("Characters reduced.");
384        log::debug!("Analysing representation symmetry for a real-space function... Done.");
385        res
386    }
387}
qsym2/sandbox/target/real_space_function/real_space_function_analysis.rs

qsym2/sandbox/target/real_space_function/
real_space_function_analysis.rs
