qsym2/symmetry/symmetry_transformation/
mod.rs

1//! Transformations under symmetry operations.
2
3use std::error::Error;
4use std::fmt;
5
6use anyhow::format_err;
7use itertools::Itertools;
8use nalgebra::Vector3;
9use ndarray::{Array, Array2, Axis, RemoveAxis};
10use num_complex::Complex;
11#[cfg(feature = "python")]
12use pyo3::prelude::*;
13use serde::{Deserialize, Serialize};
14
15use crate::angmom::sh_conversion::{sh_cart2r, sh_r2cart};
16use crate::angmom::sh_rotation_3d::rlmat;
17use crate::angmom::spinor_rotation_3d::dmat_angleaxis;
18use crate::basis::ao::{BasisAngularOrder, CartOrder, PureOrder, ShellOrder, SpinorOrder};
19use crate::permutation::{PermutableCollection, Permutation};
20use crate::symmetry::symmetry_element::symmetry_operation::{
21    SpecialSymmetryTransformation, SymmetryOperation,
22};
23
24#[cfg(test)]
25#[path = "symmetry_transformation_tests.rs"]
26mod symmetry_transformation_tests;
27
28// ================
29// Enum definitions
30// ================
31
32/// Enumerated type for managing the kind of symmetry transformation on an object.
33#[derive(Debug, Clone, Serialize, Deserialize)]
34#[cfg_attr(feature = "python", pyclass)]
35pub enum SymmetryTransformationKind {
36    /// Spatial-only transformation.
37    Spatial,
38
39    /// Spatial-only transformation but with spin-including time reversal.
40    SpatialWithSpinTimeReversal,
41
42    /// Spin-only transformation.
43    Spin,
44
45    /// Spin-spatial coupled transformation.
46    SpinSpatial,
47}
48
49impl Default for SymmetryTransformationKind {
50    fn default() -> Self {
51        SymmetryTransformationKind::Spatial
52    }
53}
54
55impl fmt::Display for SymmetryTransformationKind {
56    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
57        match self {
58            Self::Spatial => write!(f, "Spatial-only transformation"),
59            Self::SpatialWithSpinTimeReversal => write!(
60                f,
61                "Spatial-only transformation but with spin-including time reversal"
62            ),
63            Self::Spin => write!(f, "Spin-only transformation"),
64            Self::SpinSpatial => write!(f, "Spin-spatial coupled transformation"),
65        }
66    }
67}
68
69// =================
70// Trait definitions
71// =================
72
73#[derive(Debug, Clone)]
74pub struct TransformationError(pub String);
75
76impl fmt::Display for TransformationError {
77    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
78        write!(f, "Transformation error: {}", self.0)
79    }
80}
81
82impl Error for TransformationError {}
83
84/// Trait for spatial unitary transformation. A spatial unitary transformation also permutes
85/// off-origin sites.
86pub trait SpatialUnitaryTransformable: Clone {
87    // ----------------
88    // Required methods
89    // ----------------
90    /// Performs a spatial transformation in-place.
91    ///
92    /// # Arguments
93    ///
94    /// * `rmat` - The three-dimensional representation matrix of the transformation in the basis
95    /// of coordinate *functions* $`(y, z, x)`$.
96    /// * `perm` - An optional permutation describing how any off-origin sites are permuted amongst
97    /// each other under the transformation.
98    fn transform_spatial_mut(
99        &mut self,
100        rmat: &Array2<f64>,
101        perm: Option<&Permutation<usize>>,
102    ) -> Result<&mut Self, TransformationError>;
103
104    // ----------------
105    // Provided methods
106    // ----------------
107    /// Performs a spatial transformation and returns the transformed result.
108    ///
109    /// # Arguments
110    ///
111    /// * `rmat` - The three-dimensional representation matrix of the transformation in the basis
112    /// of coordinate *functions* $`(y, z, x)`$.
113    ///
114    /// # Returns
115    ///
116    /// The transformed result.
117    fn transform_spatial(
118        &self,
119        rmat: &Array2<f64>,
120        perm: Option<&Permutation<usize>>,
121    ) -> Result<Self, TransformationError> {
122        let mut tself = self.clone();
123        tself.transform_spatial_mut(rmat, perm)?;
124        Ok(tself)
125    }
126}
127
128/// Trait for spin unitary transformations. A spin unitary transformation has no spatial effects.
129pub trait SpinUnitaryTransformable: Clone {
130    // ----------------
131    // Required methods
132    // ----------------
133    /// Performs a spin transformation in-place.
134    ///
135    /// # Arguments
136    ///
137    /// * `dmat` - The two-dimensional representation matrix of the transformation in the basis of
138    /// the $`\{ \alpha, \beta \}`$ spinors (*i.e.* decreasing $`m`$ order).
139    fn transform_spin_mut(
140        &mut self,
141        dmat: &Array2<Complex<f64>>,
142    ) -> Result<&mut Self, TransformationError>;
143
144    // ----------------
145    // Provided methods
146    // ----------------
147    /// Performs a spin transformation and returns the transformed result.
148    ///
149    /// # Arguments
150    ///
151    /// * `dmat` - The two-dimensional representation matrix of the transformation in the basis of
152    /// the $`\{ \alpha, \beta \}`$ spinors (*i.e.* decreasing $`m`$ order).
153    ///
154    /// # Returns
155    ///
156    /// The transformed result.
157    fn transform_spin(&self, dmat: &Array2<Complex<f64>>) -> Result<Self, TransformationError> {
158        let mut tself = self.clone();
159        tself.transform_spin_mut(dmat)?;
160        Ok(tself)
161    }
162}
163
164/// Trait for complex-conjugation transformations.
165pub trait ComplexConjugationTransformable: Clone {
166    // ----------------
167    // Required methods
168    // ----------------
169    /// Performs a complex conjugation in-place.
170    fn transform_cc_mut(&mut self) -> Result<&mut Self, TransformationError>;
171
172    // ----------------
173    // Provided methods
174    // ----------------
175    /// Performs a complex conjugation and returns the complex-conjugated result.
176    ///
177    /// # Returns
178    ///
179    /// The complex-conjugated result.
180    fn transform_cc(&self) -> Result<Self, TransformationError> {
181        let mut tself = self.clone();
182        tself.transform_cc_mut()?;
183        Ok(tself)
184    }
185}
186
187/// Trait for time-reversal transformations.
188///
189/// This trait has a blanket implementation for any implementor of the [`SpinUnitaryTransformable`]
190/// trait and the [`ComplexConjugationTransformable`] trait together with the
191/// [`DefaultTimeReversalTransformable`] marker trait.
192pub trait TimeReversalTransformable: ComplexConjugationTransformable {
193    // ----------------
194    // Required methods
195    // ----------------
196    /// Performs a time-reversal transformation in-place.
197    fn transform_timerev_mut(&mut self) -> Result<&mut Self, TransformationError>;
198
199    // ----------------
200    // Provided methods
201    // ----------------
202    /// Performs a time-reversal transformation and returns the time-reversed result.
203    ///
204    /// # Returns
205    ///
206    /// The time-reversed result.
207    fn transform_timerev(&self) -> Result<Self, TransformationError> {
208        let mut tself = self.clone();
209        tself.transform_timerev_mut()?;
210        Ok(tself)
211    }
212}
213
214/// Trait for transformations using [`SymmetryOperation`].
215pub trait SymmetryTransformable:
216    SpatialUnitaryTransformable + SpinUnitaryTransformable + TimeReversalTransformable
217{
218    // ----------------
219    // Required methods
220    // ----------------
221    /// Determines the permutation of sites (*e.g.* atoms in molecules) due to the action of a
222    /// symmetry operation.
223    ///
224    /// # Arguments
225    ///
226    /// * `symop` - A symmetry operation.
227    ///
228    /// # Returns
229    ///
230    /// The resultant site permutation under the action of `symop`, or an error if no such
231    /// permutation can be found.
232    fn sym_permute_sites_spatial(
233        &self,
234        symop: &SymmetryOperation,
235    ) -> Result<Permutation<usize>, TransformationError>;
236
237    // ----------------
238    // Provided methods
239    // ----------------
240    /// Performs a spatial transformation according to a specified symmetry operation in-place.
241    ///
242    /// Note that both $`\mathsf{SO}(3)`$ and $`\mathsf{SU}(2)`$ rotations effect the same spatial
243    /// transformation. Also note that, if the transformation contains time reversal, it will be
244    /// accompanied by a complex conjugation.
245    ///
246    /// # Arguments
247    ///
248    /// * `symop` - A symmetry operation.
249    fn sym_transform_spatial_mut(
250        &mut self,
251        symop: &SymmetryOperation,
252    ) -> Result<&mut Self, TransformationError> {
253        let rmat = symop.get_3d_spatial_matrix();
254        let perm = self.sym_permute_sites_spatial(symop)?;
255        self.transform_spatial_mut(&rmat, Some(&perm))
256            .map_err(|err| TransformationError(err.to_string()))?;
257        if symop.contains_time_reversal() {
258            self.transform_cc_mut()
259        } else {
260            Ok(self)
261        }
262    }
263
264    /// Performs a spatial transformation according to a specified symmetry operation and returns
265    /// the transformed result.
266    ///
267    /// Note that both $`\mathsf{SO}(3)`$ and $`\mathsf{SU}(2)`$ rotations effect the same spatial
268    /// transformation. Also note that, if the transformation contains time reversal, it will be
269    /// accompanied by a complex conjugation.
270    ///
271    /// # Arguments
272    ///
273    /// * `symop` - A symmetry operation.
274    ///
275    /// # Returns
276    ///
277    /// The transformed result.
278    fn sym_transform_spatial(
279        &self,
280        symop: &SymmetryOperation,
281    ) -> Result<Self, TransformationError> {
282        let mut tself = self.clone();
283        tself.sym_transform_spatial_mut(symop)?;
284        Ok(tself)
285    }
286
287    /// Performs a spatial transformation according to a specified symmetry operation in-place, but
288    /// with spin-including time reversal.
289    ///
290    /// Note that both $`\mathsf{SO}(3)`$ and $`\mathsf{SU}(2)`$ rotations effect the same spatial
291    /// transformation. Also note that, if the transformation contains time reversal, it will be
292    /// accompanied by a rotation by $`\pi`$ about the space-fixed $`y`$-axis followed by a complex
293    /// conjugation.
294    ///
295    /// # Arguments
296    ///
297    /// * `symop` - A symmetry operation.
298    fn sym_transform_spatial_with_spintimerev_mut(
299        &mut self,
300        symop: &SymmetryOperation,
301    ) -> Result<&mut Self, TransformationError> {
302        let rmat = symop.get_3d_spatial_matrix();
303        let perm = self.sym_permute_sites_spatial(symop)?;
304        self.transform_spatial_mut(&rmat, Some(&perm))
305            .map_err(|err| TransformationError(err.to_string()))?;
306        if symop.contains_time_reversal() {
307            self.transform_timerev_mut()?;
308        }
309        Ok(self)
310    }
311
312    /// Performs a spatial transformation according to a specified symmetry operation but with
313    /// spin-including time reversal and returns the transformed result.
314    ///
315    /// Note that both $`\mathsf{SO}(3)`$ and $`\mathsf{SU}(2)`$ rotations effect the same spatial
316    /// transformation. Also note that, if the transformation contains time reversal, it will be
317    /// accompanied by a rotation by $`\pi`$ about the space-fixed $`y`$-axis followed by a complex
318    /// conjugation.
319    ///
320    /// # Arguments
321    ///
322    /// * `symop` - A symmetry operation.
323    ///
324    /// # Returns
325    ///
326    /// The transformed result.
327    fn sym_transform_spatial_with_spintimerev(
328        &self,
329        symop: &SymmetryOperation,
330    ) -> Result<Self, TransformationError> {
331        let mut tself = self.clone();
332        tself.sym_transform_spatial_with_spintimerev_mut(symop)?;
333        Ok(tself)
334    }
335
336    /// Performs a spin transformation according to a specified symmetry operation in-place.
337    ///
338    /// Note that only $`\mathsf{SU}(2)`$ rotations can effect spin transformations. Also note
339    /// that, if the transformation contains a time reversal, the corresponding explicit time
340    /// reveral action will also be carried out.
341    ///
342    /// # Arguments
343    ///
344    /// * `symop` - A symmetry operation.
345    fn sym_transform_spin_mut(
346        &mut self,
347        symop: &SymmetryOperation,
348    ) -> Result<&mut Self, TransformationError> {
349        if symop.is_su2() {
350            let angle = symop.calc_pole_angle();
351            let axis = symop.calc_pole().coords;
352            let dmat = if symop.is_su2_class_1() {
353                -dmat_angleaxis(angle, axis, false)
354            } else {
355                dmat_angleaxis(angle, axis, false)
356            };
357            self.transform_spin_mut(&dmat)?;
358        }
359        if symop.contains_time_reversal() {
360            self.transform_timerev_mut()?;
361        }
362        Ok(self)
363    }
364
365    /// Performs a spin transformation according to a specified symmetry operation and returns the
366    /// transformed result.
367    ///
368    /// Note that only $`\mathsf{SU}(2)`$ rotations can effect spin transformations. Also note
369    /// that, if the transformation is antiunitary, it will be accompanied by a time reversal.
370    ///
371    /// # Arguments
372    ///
373    /// * `symop` - A symmetry operation.
374    ///
375    /// # Returns
376    ///
377    /// The transformed result.
378    fn sym_transform_spin(&self, symop: &SymmetryOperation) -> Result<Self, TransformationError> {
379        let mut tself = self.clone();
380        tself.sym_transform_spin_mut(symop)?;
381        Ok(tself)
382    }
383
384    /// Performs a coupled spin-spatial transformation according to a specified symmetry operation
385    /// in-place.
386    ///
387    /// Note that only $`\mathsf{SU}(2)`$ rotations can effect spin transformations.
388    ///
389    /// # Arguments
390    ///
391    /// * `symop` - A symmetry operation.
392    fn sym_transform_spin_spatial_mut(
393        &mut self,
394        symop: &SymmetryOperation,
395    ) -> Result<&mut Self, TransformationError> {
396        // We cannot do the following, because each of the two methods carries out its own
397        // antiunitary action, so we'd be double-acting the antiunitary action.
398        // self.sym_transform_spatial_mut(symop)?
399        //     .sym_transform_spin_mut(symop)
400
401        // Spatial
402        let rmat = symop.get_3d_spatial_matrix();
403        let perm = self.sym_permute_sites_spatial(symop)?;
404        self.transform_spatial_mut(&rmat, Some(&perm))
405            .map_err(|err| TransformationError(err.to_string()))?;
406
407        // Spin -- only SU(2) rotations can effect spin transformations.
408        if symop.is_su2() {
409            let angle = symop.calc_pole_angle();
410            let axis = symop.calc_pole().coords;
411            let dmat = if symop.is_su2_class_1() {
412                -dmat_angleaxis(angle, axis, false)
413            } else {
414                dmat_angleaxis(angle, axis, false)
415            };
416            self.transform_spin_mut(&dmat)?;
417        }
418
419        // Time reversal, if any.
420        if symop.contains_time_reversal() {
421            self.transform_timerev_mut()?;
422        }
423        Ok(self)
424    }
425
426    /// Performs a coupled spin-spatial transformation according to a specified symmetry operation
427    /// and returns the transformed result.
428    ///
429    /// Note that only $`\mathsf{SU}(2)`$ rotations can effect spin transformations.
430    ///
431    /// # Arguments
432    ///
433    /// * `symop` - A symmetry operation.
434    ///
435    /// # Returns
436    ///
437    /// The transformed result.
438    fn sym_transform_spin_spatial(
439        &self,
440        symop: &SymmetryOperation,
441    ) -> Result<Self, TransformationError> {
442        let mut tself = self.clone();
443        tself.sym_transform_spin_spatial_mut(symop)?;
444        Ok(tself)
445    }
446}
447
448// ----------------------
449// Blanket implementation
450// ----------------------
451
452/// Marker trait indicating that the implementing type should get the blanket implementation for
453/// [`TimeReversalTransformable`].
454pub trait DefaultTimeReversalTransformable {}
455
456impl<T> TimeReversalTransformable for T
457where
458    T: DefaultTimeReversalTransformable
459        + SpinUnitaryTransformable
460        + ComplexConjugationTransformable,
461{
462    /// Performs a time-reversal transformation in-place.
463    ///
464    /// The default implementation of the time-reversal transformation for any type that implements
465    /// [`SpinUnitaryTransformable`] and [`ComplexConjugationTransformable`] is a spin rotation by
466    /// $`\pi`$ about the space-fixed $`y`$-axis followed by a complex conjugation.
467    fn transform_timerev_mut(&mut self) -> Result<&mut Self, TransformationError> {
468        let dmat_y = dmat_angleaxis(std::f64::consts::PI, Vector3::y(), false);
469        self.transform_spin_mut(&dmat_y)?.transform_cc_mut()
470    }
471}
472
473// =========
474// Functions
475// =========
476
477/// Permutes the generalised rows of an array along one or more dimensions.
478///
479/// Each generalised row corresponds to a basis function, and consecutive generalised rows
480/// corresponding to basis functions localised on a single atom are grouped together and then
481/// permuted according to the permutation of the atoms.
482///
483/// # Arguments
484///
485/// * `arr` - A coefficient array of any dimensions.
486/// * `atom_perm` - A permutation for the atoms.
487/// * `axes` - The dimensions along which the generalised rows are to be permuted. The number of
488/// generalised rows along each of these dimensions *must* be equal to the number of functions in
489/// the basis.
490/// * `bao` - A structure specifying the angular order of the underlying basis.
491///
492/// # Returns
493///
494/// The permuted array.
495///
496/// # Panics
497///
498/// Panics if the number of generalised rows along any of the dimensions in `axes` does not match
499/// the number of functions in the basis, or if the permutation rank does not match the number of
500/// atoms in the basis.
501pub(crate) fn permute_array_by_atoms<T, D>(
502    arr: &Array<T, D>,
503    atom_perm: &Permutation<usize>,
504    axes: &[Axis],
505    bao: &BasisAngularOrder,
506) -> Array<T, D>
507where
508    D: RemoveAxis,
509    T: Clone,
510{
511    assert_eq!(
512        atom_perm.rank(),
513        bao.n_atoms(),
514        "The rank of permutation does not match the number of atoms in the basis."
515    );
516    let atom_boundary_indices = bao.atom_boundary_indices();
517    let permuted_shell_indices: Vec<usize> = atom_perm
518        .image()
519        .iter()
520        .flat_map(|&i| {
521            let (shell_min, shell_max) = atom_boundary_indices[i];
522            shell_min..shell_max
523        })
524        .collect();
525
526    let mut r = arr.clone();
527    for axis in axes {
528        assert_eq!(
529            arr.shape()[axis.0],
530            bao.n_funcs(),
531            "The number of generalised rows along {axis:?} in the given array does not match the number of basis functions, {}.", bao.n_funcs()
532        );
533        r = r.select(*axis, &permuted_shell_indices);
534    }
535    r
536}
537
538/// Assembles spherical-harmonic rotation matrices for all shells.
539///
540/// # Arguments
541///
542/// * `bao` - A structure specifying the angular order of the underlying basis.
543/// * `rmat` - The three-dimensional representation matrix of the transformation in the basis
544/// of coordinate *functions* $`(y, z, x)`$.
545/// * `perm` - An optional permutation describing how any off-origin sites are permuted amongst
546/// each other under the transformation.
547///
548/// # Returns
549///
550/// A vector of spherical-harmonic rotation matrices, one for each shells in `bao`. Non-standard
551/// orderings of functions in shells are taken into account.
552pub fn assemble_sh_rotation_3d_matrices(
553    bao: &BasisAngularOrder,
554    rmat: &Array2<f64>,
555    perm: Option<&Permutation<usize>>,
556) -> Result<Vec<Array2<f64>>, anyhow::Error> {
557    let pbao = if let Some(p) = perm {
558        bao.permute(p)?
559    } else {
560        bao.clone()
561    };
562    let mut rls = vec![Array2::<f64>::eye(1), rmat.clone()];
563    let lmax = pbao
564        .basis_shells()
565        .map(|shl| shl.l)
566        .max()
567        .expect("The maximum angular momentum cannot be found.");
568    for l in 2..=lmax {
569        let rl = rlmat(
570            l,
571            rmat,
572            rls.last()
573                .expect("The representation matrix for the last angular momentum cannot be found."),
574        );
575        rls.push(rl);
576    }
577
578    // All matrices in `rls` are in increasing-m order by default. See the function `rlmat` for
579    // the origin of this order. Hence, conversion matrices must also honour this.
580    let cart2rss_lex: Vec<Vec<Array2<f64>>> = (0..=lmax)
581        .map(|lcart| sh_cart2r(lcart, &CartOrder::lex(lcart), true, PureOrder::increasingm))
582        .collect();
583    let r2cartss_lex: Vec<Vec<Array2<f64>>> = (0..=lmax)
584        .map(|lcart| sh_r2cart(lcart, &CartOrder::lex(lcart), true, PureOrder::increasingm))
585        .collect();
586
587    let rmats_res = pbao.basis_shells()
588        .map(|shl| {
589            let l = usize::try_from(shl.l).unwrap_or_else(|_| {
590                panic!(
591                    "Unable to convert the angular momentum order `{}` to `usize`.",
592                    shl.l
593                );
594            });
595            let po_il = PureOrder::increasingm(shl.l);
596            match &shl.shell_order {
597                ShellOrder::Pure(pureorder) => {
598                    // Spherical functions.
599                    let rl = rls[l].clone();
600                    if *pureorder != po_il {
601                        // `rl` is in increasing-m order by default. See the function `rlmat` for
602                        // the origin of this order.
603                        let perm = pureorder
604                            .get_perm_of(&po_il)
605                            .expect("Unable to obtain the permutation that maps `pureorder` to the increasing order.");
606                        Ok(rl.select(Axis(0), &perm.image()).select(Axis(1), &perm.image()))
607                    } else {
608                        Ok(rl)
609                    }
610                }
611                ShellOrder::Cart(cart_order) => {
612                    // Cartesian functions. Convert them to real solid harmonics first, then
613                    // applying the transformation, then convert back.
614                    // The actual Cartesian order will be taken into account.
615
616                    // Perform the conversion using lexicographic order first. This allows for the
617                    // conversion matrices to be computed only once in the lexicographic order.
618                    let cart2rs = &cart2rss_lex[l];
619                    let r2carts = &r2cartss_lex[l];
620                    let rl = cart2rs.iter().zip(r2carts.iter()).enumerate().fold(
621                        Array2::zeros((cart_order.ncomps(), cart_order.ncomps())),
622                        |acc, (i, (xmat, wmat))| {
623                            let lpure = l - 2 * i;
624                            acc + wmat.dot(&rls[lpure]).dot(xmat)
625                        },
626                    );
627                    let lex_cart_order = CartOrder::lex(shl.l);
628
629                    // Now deal with the actual Cartesian order by permutations.
630                    if *cart_order != lex_cart_order {
631                        // `rl` is in lexicographic order (because of `wmat` and `xmat`) by default.
632                        // Consider a transformation R and its representation matrix D in a
633                        // lexicographically-ordered Cartesian basis b collected in a row vector.
634                        // Then,
635                        //      R b = b D.
636                        // If we now permute the basis functions in b by a permutation π, then the
637                        // representation matrix for R changes:
638                        //      R πb = πb D(π).
639                        // To relate D(π) to D, we first note the representation matrix for π, P:
640                        //      πb = π b = b P,
641                        // which, when acts on a left row vector, permutes its entries normally, but
642                        // when acts on a right column vector, permutes its entries inversely.
643                        // Then,
644                        //      R πb = R b P = b P D(π) => R b = b PD(π)P^(-1).
645                        // Thus,
646                        //      D(π) = P^(-1)DP,
647                        // i.e., to obtain D(π), we permute the rows and columns of D normally
648                        // according to π.
649                        let perm = lex_cart_order
650                            .get_perm_of(cart_order)
651                            .unwrap_or_else(
652                                || panic!("Unable to find a permutation to map `{lex_cart_order}` to `{cart_order}`.")
653                            );
654                        Ok(rl.select(Axis(0), perm.image())
655                            .select(Axis(1), perm.image()))
656                    } else {
657                        Ok(rl)
658                    }
659                }
660                ShellOrder::Spinor(_) => Err(format_err!("Spinor shells cannot have transformation matrices with spherical-harmonic bases."))
661            }
662        })
663        .collect::<Result<Vec<Array2<f64>>, _>>();
664    rmats_res
665}
666
667/// Assembles spinor rotation matrices for all shells, which also include the unitary part of time
668/// reversal, if any. This also handles whether a spinor shell is even or odd with respect to
669/// spatial inversion.
670///
671/// # Arguments
672///
673/// * `bao` - A structure specifying the angular order of the underlying basis.
674/// * `symop` - A symmetry operation representing the transformation.
675/// * `perm` - An optional permutation describing how any off-origin sites are permuted amongst
676/// each other under the transformation.
677///
678/// # Returns
679///
680/// A vector of spinor rotation matrices, one for each shells in `bao`. Non-standard orderings of
681/// functions in shells are taken into account.
682pub(crate) fn assemble_spinor_rotation_matrices(
683    bao: &BasisAngularOrder,
684    symop: &SymmetryOperation,
685    perm: Option<&Permutation<usize>>,
686) -> Result<Vec<Array2<Complex<f64>>>, anyhow::Error> {
687    let pbao = if let Some(p) = perm {
688        bao.permute(p)?
689    } else {
690        bao.clone()
691    };
692    let two_j_max = pbao
693        .basis_shells()
694        .map(|shl| match shl.shell_order {
695            ShellOrder::Pure(_) => 2 * shl.l,
696            ShellOrder::Cart(_) => 2 * shl.l,
697            ShellOrder::Spinor(_) => shl.l,
698        })
699        .max()
700        .expect("The maximum rank cannot be found.");
701    let r2js = (0..=two_j_max)
702        .map(|two_j| symop.get_wigner_matrix(two_j, true))
703        .collect::<Result<Vec<_>, _>>()?;
704
705    // All matrices in `rls` are in increasing-m order by default. See the function `get_wigner_matrix` for
706    // the origin of this order. Hence, conversion matrices must also honour this.
707    let cart2rss_lex: Vec<Vec<Array2<Complex<f64>>>> = (0..=two_j_max)
708        .step_by(2)
709        .map(|two_lcart| {
710            sh_cart2r(
711                two_lcart.div_euclid(2),
712                &CartOrder::lex(two_lcart.div_euclid(2)),
713                true,
714                PureOrder::increasingm,
715            )
716            .into_iter()
717            .map(|mat| mat.mapv(Complex::<f64>::from))
718            .collect_vec()
719        })
720        .collect();
721    let r2cartss_lex: Vec<Vec<Array2<Complex<f64>>>> = (0..=two_j_max)
722        .map(|two_lcart| {
723            sh_r2cart(
724                two_lcart.div_euclid(2),
725                &CartOrder::lex(two_lcart.div_euclid(2)),
726                true,
727                PureOrder::increasingm,
728            )
729            .into_iter()
730            .map(|mat| mat.mapv(Complex::<f64>::from))
731            .collect_vec()
732        })
733        .collect();
734
735    let rmats_res = pbao.basis_shells()
736        .map(|shl| {
737            let l =
738                usize::try_from(shl.l).unwrap_or_else(|_| {
739                panic!(
740                    "Unable to convert the rank `{}` to `usize`.",
741                    shl.l
742                );
743            });
744            match &shl.shell_order {
745                ShellOrder::Pure(pure_order) => {
746                    // Spherical functions.
747                    let po_il = PureOrder::increasingm(shl.l);
748                    let rl = r2js[2*l].clone();
749                    if *pure_order != po_il {
750                        // `rl` is in increasing-m order by default. See the function `rlmat` for
751                        // the origin of this order.
752                        let perm = pure_order
753                            .get_perm_of(&po_il)
754                            .expect("Unable to obtain the permutation that maps `pureorder` to the increasing order.");
755                        Ok(rl.select(Axis(0), &perm.image()).select(Axis(1), &perm.image()))
756                    } else {
757                        Ok(rl)
758                    }
759                }
760                ShellOrder::Cart(cart_order) => {
761                    // Cartesian functions. Convert them to real solid harmonics first, then
762                    // applying the transformation, then convert back.
763                    // The actual Cartesian order will be taken into account.
764
765                    // Perform the conversion using lexicographic order first. This allows for the
766                    // conversion matrices to be computed only once in the lexicographic order.
767                    let cart2rs = &cart2rss_lex[l];
768                    let r2carts = &r2cartss_lex[l];
769                    let rl = cart2rs.iter().zip(r2carts.iter()).enumerate().fold(
770                        Array2::zeros((cart_order.ncomps(), cart_order.ncomps())),
771                        |acc, (i, (xmat, wmat))| {
772                            let lpure = l - 2 * i;
773                            acc + wmat.dot(&r2js[2*lpure]).dot(xmat)
774                        },
775                    );
776                    let lex_cart_order = CartOrder::lex(shl.l);
777
778                    // Now deal with the actual Cartesian order by permutations.
779                    if *cart_order != lex_cart_order {
780                        // `rl` is in lexicographic order (because of `wmat` and `xmat`) by default.
781                        // Consider a transformation R and its representation matrix D in a
782                        // lexicographically-ordered Cartesian basis b collected in a row vector.
783                        // Then,
784                        //      R b = b D.
785                        // If we now permute the basis functions in b by a permutation π, then the
786                        // representation matrix for R changes:
787                        //      R πb = πb D(π).
788                        // To relate D(π) to D, we first note the representation matrix for π, P:
789                        //      πb = π b = b P,
790                        // which, when acts on a left row vector, permutes its entries normally, but
791                        // when acts on a right column vector, permutes its entries inversely.
792                        // Then,
793                        //      R πb = R b P = b P D(π) => R b = b PD(π)P^(-1).
794                        // Thus,
795                        //      D(π) = P^(-1)DP,
796                        // i.e., to obtain D(π), we permute the rows and columns of D normally
797                        // according to π.
798                        let perm = lex_cart_order
799                            .get_perm_of(cart_order)
800                            .unwrap_or_else(
801                                || panic!("Unable to find a permutation to map `{lex_cart_order}` to `{cart_order}`.")
802                            );
803                        Ok(rl.select(Axis(0), perm.image())
804                            .select(Axis(1), perm.image()))
805                    } else {
806                        Ok(rl)
807                    }
808                }
809                ShellOrder::Spinor(spinor_order) => {
810                    // Spinor functions. l = two_j.
811                    let so_il = SpinorOrder::increasingm(shl.l, spinor_order.even);
812                    let r2j = r2js[l].clone();
813                    if *spinor_order != so_il {
814                        // `rl` is in increasing-m order by default. See the function `rlmat` for
815                        // the origin of this order.
816                        let perm = spinor_order
817                            .get_perm_of(&so_il)
818                            .expect("Unable to obtain the permutation that maps `spinor_order` to the increasing order.");
819                        if !so_il.even && !symop.is_proper() {
820                            Ok(-r2j.select(Axis(0), &perm.image()).select(Axis(1), &perm.image()))
821                        } else {
822                            Ok(r2j.select(Axis(0), &perm.image()).select(Axis(1), &perm.image()))
823                        }
824                    } else {
825                        if !spinor_order.even && !symop.is_proper() {
826                            Ok(-r2j)
827                        } else {
828                            Ok(r2j)
829                        }
830                    }
831                }
832            }
833        })
834        .collect::<Result<Vec<Array2<Complex<f64>>>, _>>();
835    rmats_res
836}