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}