qsym2/angmom/
spinor_rotation_3d.rs

1//! Three-dimensional rotations of spinors.
2
3use std::cmp;
4use std::fmt;
5
6use anyhow::{self, ensure};
7use approx;
8use factorial::Factorial;
9use nalgebra::Vector3;
10use ndarray::{array, Array2, Axis};
11use num::{BigUint, Complex, Zero};
12use num_traits::ToPrimitive;
13use serde::{Deserialize, Serialize};
14
15use crate::basis::ao::BasisAngularOrder;
16use crate::auxiliary::geometry::normalise_rotation_angle;
17
18#[cfg(test)]
19#[path = "spinor_rotation_3d_tests.rs"]
20mod spinor_rotation_3d_tests;
21
22// ================
23// Trait definition
24// ================
25
26/// Trait for defining the spin-spatial structure of the coefficient matrices.
27pub trait StructureConstraint {
28    // ----------------
29    // Required methods
30    // ----------------
31
32    /// The total number of coefficient matrices expected in this constraint.
33    fn n_coefficient_matrices(&self) -> usize;
34
35    /// The number of isostructural components explicitly specified by each coefficient matrix. The
36    /// basis functions described by any accompanying basis angular order information will pertain to any
37    /// one of the isostructural components.
38    fn n_explicit_comps_per_coefficient_matrix(&self) -> usize;
39
40    /// The number of isostructural components implicitly specified by each coefficient matrix.
41    fn n_implicit_comps_per_coefficient_matrix(&self) -> usize;
42
43    /// The total number of isostructural components given by the product of
44    /// [`Self::n_coefficient_matrices`] and [`Self::n_implicit_comps_per_coefficient_matrix`].
45    fn n_comps(&self) -> usize {
46        self.n_coefficient_matrices() * self.n_implicit_comps_per_coefficient_matrix()
47    }
48
49    // ----------------
50    // Provided methods
51    // ----------------
52    /// The implicit factor given by the ratio between
53    /// [`Self::n_implicit_comps_per_coefficient_matrix`] and
54    /// [`Self::n_explicit_comps_per_coefficient_matrix`]. It will be checked that this ratio
55    /// yields an integer.
56    fn implicit_factor(&self) -> Result<usize, anyhow::Error> {
57        ensure!(
58            self.n_implicit_comps_per_coefficient_matrix()
59                .rem_euclid(self.n_explicit_comps_per_coefficient_matrix())
60                == 0
61        );
62        let implicit_factor = self
63            .n_implicit_comps_per_coefficient_matrix()
64            .div_euclid(self.n_explicit_comps_per_coefficient_matrix());
65        Ok(implicit_factor)
66    }
67}
68
69// ================
70// Enum definitions
71// ================
72
73// --------------
74// SpinConstraint
75// --------------
76
77/// Enumerated type to manage spin constraints and spin space information for the conventional
78/// treatment of spin and spatial degrees of freedom in a decoupled manner.
79#[derive(Clone, Debug, Hash, PartialEq, Eq, Serialize, Deserialize)]
80pub enum SpinConstraint {
81    /// Variant for restricted spin constraint: the spatial parts of all spin spaces are identical.
82    /// The associated value is the number of spin spaces.
83    Restricted(u16),
84
85    /// Variant for unrestricted spin constraint: the spatial parts of different spin spaces are
86    /// different, but spin collinearity is maintained. The associated values are the number of spin
87    /// spaces (*i.e.* the number of different spatial parts that are handled separately) and a
88    /// boolean indicating if the spin spaces are arranged in increasing $`m`$ order.
89    Unrestricted(u16, bool),
90
91    /// Variant for generalised spin constraint: the spatial parts of different spin spaces are
92    /// different, and no spin collinearity is imposed. The associated values are the number of spin
93    /// spaces and a boolean indicating if the spin spaces are arranged in increasing $`m`$ order.
94    Generalised(u16, bool),
95}
96
97impl StructureConstraint for SpinConstraint {
98    fn n_coefficient_matrices(&self) -> usize {
99        match self {
100            SpinConstraint::Restricted(_) => 1,
101            SpinConstraint::Unrestricted(nspins, _) => *nspins as usize,
102            SpinConstraint::Generalised(_, _) => 1,
103        }
104    }
105
106    fn n_explicit_comps_per_coefficient_matrix(&self) -> usize {
107        match self {
108            SpinConstraint::Restricted(_) => 1,
109            SpinConstraint::Unrestricted(_, _) => 1,
110            SpinConstraint::Generalised(nspins, _) => *nspins as usize,
111        }
112    }
113
114    fn n_implicit_comps_per_coefficient_matrix(&self) -> usize {
115        match self {
116            SpinConstraint::Restricted(nspins) => *nspins as usize,
117            SpinConstraint::Unrestricted(_, _) => 1,
118            SpinConstraint::Generalised(nspins, _) => *nspins as usize,
119        }
120    }
121}
122
123impl fmt::Display for SpinConstraint {
124    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
125        match self {
126            Self::Restricted(nspins) => write!(
127                f,
128                "Restricted ({} spin {})",
129                nspins,
130                if *nspins == 1 { "space" } else { "spaces" }
131            ),
132            Self::Unrestricted(nspins, increasingm) => write!(
133                f,
134                "Unrestricted ({} spin {}, {} m)",
135                nspins,
136                if *nspins == 1 { "space" } else { "spaces" },
137                if *increasingm {
138                    "increasing"
139                } else {
140                    "decreasing"
141                }
142            ),
143            Self::Generalised(nspins, increasingm) => write!(
144                f,
145                "Generalised ({} spin {}, {} m)",
146                nspins,
147                if *nspins == 1 { "space" } else { "spaces" },
148                if *increasingm {
149                    "increasing"
150                } else {
151                    "decreasing"
152                }
153            ),
154        }
155    }
156}
157
158// ----------------
159// SpinOrbitCoupled
160// ----------------
161
162/// Enumerated type to manage spin--orbit-coupled layout in the coupled treatment of spin and
163/// spatial degrees of freedom.
164#[derive(Clone, Debug, Hash, PartialEq, Eq, Serialize, Deserialize)]
165pub enum SpinOrbitCoupled {
166    /// Variant for $`j`$-adapted basis functions where each shell consists of $`\ket{j, m_j}`$
167    /// functions. The associated value specifies the total number of duplications of the
168    /// $`j`$-adapted basis (*e.g.* `2` in Dirac--Hartree--Fock). The order of $`m_j`$ in each
169    /// shell should be specified in a [`BasisAngularOrder`] structure elsewhere.
170    JAdapted(u16),
171}
172
173impl StructureConstraint for SpinOrbitCoupled {
174    fn n_coefficient_matrices(&self) -> usize {
175        match self {
176            SpinOrbitCoupled::JAdapted(_) => 1,
177        }
178    }
179
180    fn n_explicit_comps_per_coefficient_matrix(&self) -> usize {
181        match self {
182            SpinOrbitCoupled::JAdapted(ncomps) => *ncomps as usize,
183        }
184    }
185
186    fn n_implicit_comps_per_coefficient_matrix(&self) -> usize {
187        match self {
188            SpinOrbitCoupled::JAdapted(ncomps) => *ncomps as usize,
189        }
190    }
191}
192
193impl fmt::Display for SpinOrbitCoupled {
194    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
195        match self {
196            Self::JAdapted(ncomps) => write!(
197                f,
198                "Spin--orbit-coupled j-adapted ({} {})",
199                ncomps,
200                if *ncomps == 1 {
201                    "component"
202                } else {
203                    "components"
204                },
205            ),
206        }
207    }
208}
209
210// =========
211// Functions
212// =========
213
214/// Returns an element in the Wigner rotation matrix for $`j = 1/2`$ defined by
215///
216/// ```math
217///     \hat{R}(\alpha, \beta, \gamma) \ket{\tfrac{1}{2}m}
218///     = \sum_{m'} \ket{\tfrac{1}{2}m'} D^{(1/2)}_{m'm}(\alpha, \beta, \gamma).
219/// ```
220///
221/// # Arguments
222///
223/// * `mdashi` - Index for $`m'`$ given by $`m'+\tfrac{1}{2}`$.
224/// * `mi` - Index for $`m`$ given by $`m+\tfrac{1}{2}`$.
225/// * `euler_angles` - A triplet of Euler angles $`(\alpha, \beta, \gamma)`$ in radians, following
226/// the Whitaker convention, *i.e.* $`z_2-y-z_1`$ (extrinsic rotations).
227///
228/// # Returns
229///
230/// The element $`D^{(1/2)}_{m'm}(\alpha, \beta, \gamma)`$.
231fn dmat_euler_element(mdashi: usize, mi: usize, euler_angles: (f64, f64, f64)) -> Complex<f64> {
232    assert!(mdashi == 0 || mdashi == 1, "mdashi can only be 0 or 1.");
233    assert!(mi == 0 || mi == 1, "mi can only be 0 or 1.");
234    let (alpha, beta, gamma) = euler_angles;
235    let d = if (mi, mdashi) == (1, 1) {
236        // m = 1/2, mdash = 1/2
237        (beta / 2.0).cos()
238    } else if (mi, mdashi) == (1, 0) {
239        // m = 1/2, mdash = -1/2
240        (beta / 2.0).sin()
241    } else if (mi, mdashi) == (0, 1) {
242        // m = -1/2, mdash = 1/2
243        -(beta / 2.0).sin()
244    } else if (mi, mdashi) == (0, 0) {
245        // m = -1/2, mdash = -1/2
246        (beta / 2.0).cos()
247    } else {
248        panic!("Invalid mi and/or mdashi.");
249    };
250
251    let alpha_basic = alpha.rem_euclid(2.0 * std::f64::consts::PI);
252    let gamma_basic = gamma.rem_euclid(2.0 * std::f64::consts::PI);
253    let i = Complex::<f64>::i();
254    let mut prefactor = (-i
255        * (alpha_basic
256            * (mdashi
257                .to_f64()
258                .unwrap_or_else(|| panic!("Unable to convert `{mdashi}` to `f64`."))
259                - 0.5)
260            + gamma_basic
261                * (mi
262                    .to_f64()
263                    .unwrap_or_else(|| panic!("Unable to convert `{mi}` to `f64`."))
264                    - 0.5)))
265        .exp();
266
267    // Half-integer j = 1/2; double-group behaviours possible.
268    let alpha_double = approx::relative_eq!(
269        alpha.div_euclid(2.0 * std::f64::consts::PI).rem_euclid(2.0),
270        1.0,
271        epsilon = 1e-14,
272        max_relative = 1e-14
273    );
274    let gamma_double = approx::relative_eq!(
275        gamma.div_euclid(2.0 * std::f64::consts::PI).rem_euclid(2.0),
276        1.0,
277        epsilon = 1e-14,
278        max_relative = 1e-14
279    );
280    if alpha_double != gamma_double {
281        prefactor *= -1.0;
282    }
283    prefactor * d
284}
285
286/// Returns the Wigner rotation matrix for $`j = 1/2`$ whose elements are defined by
287///
288/// ```math
289/// \hat{R}(\alpha, \beta, \gamma) \ket{\tfrac{1}{2}m}
290/// = \sum_{m'} \ket{\tfrac{1}{2}m'} D^{(1/2)}_{m'm}(\alpha, \beta, \gamma).
291/// ```
292///
293/// # Arguments
294///
295/// * `euler_angles` - A triplet of Euler angles $`(\alpha, \beta, \gamma)`$ in radians, following
296/// the Whitaker convention, *i.e.* $`z_2-y-z_1`$ (extrinsic rotations).
297/// * `increasingm` - If `true`, the rows and columns of $`\mathbf{D}^{(1/2)}`$ are
298/// arranged in increasing order of $`m_l = -l, \ldots, l`$. If `false`, the order is reversed:
299/// $`m_l = l, \ldots, -l`$. The recommended default is `false`, in accordance with convention.
300///
301/// # Returns
302///
303/// The matrix $`\mathbf{D}^{(1/2)}(\alpha, \beta, \gamma)`$.
304#[must_use]
305pub fn dmat_euler(euler_angles: (f64, f64, f64), increasingm: bool) -> Array2<Complex<f64>> {
306    let mut dmat = Array2::<Complex<f64>>::zeros((2, 2));
307    for mdashi in 0..2 {
308        for mi in 0..2 {
309            dmat[(mdashi, mi)] = dmat_euler_element(mdashi, mi, euler_angles);
310        }
311    }
312    if !increasingm {
313        dmat.invert_axis(Axis(0));
314        dmat.invert_axis(Axis(1));
315    }
316    dmat
317}
318
319/// Returns the Wigner rotation matrix for $`j = 1/2`$ whose elements are defined by
320///
321/// ```math
322/// \hat{R}(\phi\hat{\mathbf{n}}) \ket{\tfrac{1}{2}m}
323/// = \sum_{m'} \ket{\tfrac{1}{2}m'} D^{(1/2)}_{m'm}(\phi\hat{\mathbf{n}}).
324/// ```
325///
326/// The parametrisation of $`\mathbf{D}^{(1/2)}`$ by $`\phi`$ and $`\hat{\mathbf{n}}`$ is given
327/// in (**4**-9.12) of Altmann, S. L. Rotations, Quaternions, and Double Groups. (Dover
328/// Publications, Inc., 2005).
329///
330/// # Arguments
331///
332/// * `angle` - The angle $`\phi`$ of the rotation in radians. A positive rotation is an
333/// anticlockwise rotation when looking down `axis`.
334/// * `axis` - A space-fixed vector defining the axis of rotation. The supplied vector will be
335/// normalised.
336/// * `increasingm` - If `true`, the rows and columns of $`\mathbf{D}^{(1/2)}`$ are
337/// arranged in increasing order of $`m_l = -l, \ldots, l`$. If `false`, the order is reversed:
338/// $`m_l = l, \ldots, -l`$. The recommended default is `false`, in accordance with convention.
339///
340/// # Returns
341///
342/// The matrix $`\mathbf{D}^{(1/2)}(\phi\hat{\mathbf{n}})`$.
343#[must_use]
344pub fn dmat_angleaxis(angle: f64, axis: Vector3<f64>, increasingm: bool) -> Array2<Complex<f64>> {
345    let normalised_axis = axis.normalize();
346    let nx = normalised_axis.x;
347    let ny = normalised_axis.y;
348    let nz = normalised_axis.z;
349
350    let i = Complex::<f64>::i();
351    let double_angle = angle.rem_euclid(4.0 * std::f64::consts::PI);
352    let mut dmat = array![
353        [
354            (double_angle / 2.0).cos() + i * nz * (double_angle / 2.0).sin(),
355            (ny - i * nx) * (double_angle / 2.0).sin()
356        ],
357        [
358            -(ny + i * nx) * (double_angle / 2.0).sin(),
359            (double_angle / 2.0).cos() - i * nz * (double_angle / 2.0).sin(),
360        ]
361    ];
362    if !increasingm {
363        dmat.invert_axis(Axis(0));
364        dmat.invert_axis(Axis(1));
365    }
366    dmat
367}
368
369/// Returns an element in the Wigner rotation matrix for an integral or half-integral
370/// $`j`$, defined by
371/// ```math
372/// \hat{R}(\alpha, \beta, \gamma) \ket{jm}
373/// = \sum_{m'} \ket{jm'} D^{(j)}_{m'm}(\alpha, \beta, \gamma)
374/// ```
375/// where $`-\pi \le \alpha \le \pi`$, $`0 \le \beta \le \pi`$, $`-\pi \le \gamma \le \pi`$.
376///
377/// The explicit expression for the elements of $`\mathbf{D}^{(j)}(\alpha, \beta, \gamma)`$
378/// is given in Professor Anthony Stone's graduate lecture notes on Angular Momentum at the
379/// University of Cambridge in 2006.
380///
381/// # Arguments
382///
383/// * `twoj` - Two times the angular momentum $`2j`$. If this is even, $`j`$ is integral; otherwise,
384/// $`j`$ is half-integral.
385/// * `mdashi` - Index for $`m'`$ given by $`m'+\tfrac{1}{2}`$.
386/// * `mi` - Index for $`m`$ given by $`m+\tfrac{1}{2}`$.
387/// * `euler_angles` - A triplet of Euler angles $`(\alpha, \beta, \gamma)`$ in radians, following
388/// the Whitaker convention, *i.e.* $`z_2-y-z_1`$ (extrinsic rotations).
389///
390/// # Returns
391///
392/// The element $`D^{(j)}_{m'm}(\alpha, \beta, \gamma)`$.
393#[allow(clippy::too_many_lines)]
394pub fn dmat_euler_gen_element(
395    twoj: u32,
396    mdashi: usize,
397    mi: usize,
398    euler_angles: (f64, f64, f64),
399) -> Complex<f64> {
400    assert!(
401        mdashi <= twoj as usize,
402        "`mdashi` must be between 0 and {twoj} (inclusive).",
403    );
404    assert!(
405        mi <= twoj as usize,
406        "`mi` must be between 0 and {twoj} (inclusive).",
407    );
408    let (alpha, beta, gamma) = euler_angles;
409    let j = f64::from(twoj) / 2.0;
410    let mdash = mdashi
411        .to_f64()
412        .unwrap_or_else(|| panic!("Unable to convert `{mdashi}` to `f64`."))
413        - j;
414    let m = mi
415        .to_f64()
416        .unwrap_or_else(|| panic!("Unable to convert `{mi}` to `f64`."))
417        - j;
418
419    let i = Complex::<f64>::i();
420    let alpha_basic = alpha.rem_euclid(2.0 * std::f64::consts::PI);
421    let gamma_basic = gamma.rem_euclid(2.0 * std::f64::consts::PI);
422    let mut prefactor = (-i * (alpha_basic * mdash + gamma_basic * m)).exp();
423
424    if twoj % 2 != 0 {
425        // Half-integer j; double-group behaviours possible.
426        let alpha_double = approx::relative_eq!(
427            alpha.div_euclid(2.0 * std::f64::consts::PI).rem_euclid(2.0),
428            1.0,
429            epsilon = 1e-14,
430            max_relative = 1e-14
431        );
432        let gamma_double = approx::relative_eq!(
433            gamma.div_euclid(2.0 * std::f64::consts::PI).rem_euclid(2.0),
434            1.0,
435            epsilon = 1e-14,
436            max_relative = 1e-14
437        );
438        if alpha_double != gamma_double {
439            prefactor *= -1.0;
440        }
441    }
442
443    // tmax = min(int(j + mdash), int(j - m))
444    // j + mdash = mdashi
445    // j - m = twoj - mi
446    let tmax = cmp::min(mdashi, twoj as usize - mi);
447
448    // tmin = max(0, int(mdash - m))
449    // mdash - m = mdashi - mi
450    let tmin = if mdashi > mi { mdashi - mi } else { 0 };
451
452    let d = (tmin..=tmax).fold(Complex::<f64>::zero(), |acc, t| {
453        // j - m = twoj - mi
454        // j - mdash = twoj - mdashi
455        let num = (BigUint::from(mdashi)
456            .checked_factorial()
457            .unwrap_or_else(|| panic!("Unable to compute the factorial of {mdashi}."))
458            * BigUint::from(twoj as usize - mdashi)
459                .checked_factorial()
460                .unwrap_or_else(|| {
461                    panic!(
462                        "Unable to compute the factorial of {}.",
463                        twoj as usize - mdashi
464                    )
465                })
466            * BigUint::from(mi)
467                .checked_factorial()
468                .unwrap_or_else(|| panic!("Unable to compute the factorial of {mi}."))
469            * BigUint::from(twoj as usize - mi)
470                .checked_factorial()
471                .unwrap_or_else(|| {
472                    panic!("Unable to compute the factorial of {}.", twoj as usize - mi)
473                }))
474        .to_f64()
475        .expect("Unable to convert a `BigUint` value to `f64`.")
476        .sqrt();
477
478        // t <= j + mdash ==> j + mdash - t = mdashi - t >= 0
479        // t <= j - m ==> j - m - t = twoj - mi - t >= 0
480        // t >= 0
481        // t >= mdash - m ==> t - (mdash - m) = t + mi - mdashi >= 0
482        let den = (BigUint::from(mdashi - t)
483            .checked_factorial()
484            .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", mdashi - t))
485            * BigUint::from(twoj as usize - mi - t)
486                .checked_factorial()
487                .unwrap_or_else(|| {
488                    panic!(
489                        "Unable to compute the factorial of {}.",
490                        twoj as usize - mi - t
491                    )
492                })
493            * BigUint::from(t)
494                .checked_factorial()
495                .unwrap_or_else(|| panic!("Unable to compute the factorial of {t}."))
496            * BigUint::from(t + mi - mdashi)
497                .checked_factorial()
498                .unwrap_or_else(|| {
499                    panic!("Unable to compute the factorial of {}.", t + mi - mdashi)
500                }))
501        .to_f64()
502        .expect("Unable to convert a `BigUint` value to `f64`.");
503
504        let trigfactor = (beta / 2.0).cos().powi(
505            i32::try_from(twoj as usize + mdashi - mi - 2 * t).unwrap_or_else(|_| {
506                panic!(
507                    "Unable to convert `{}` to `i32`.",
508                    twoj as usize + mdashi - mi - 2 * t
509                )
510            }),
511        ) * (beta / 2.0).sin().powi(
512            i32::try_from(2 * t + mi - mdashi).unwrap_or_else(|_| {
513                panic!("Unable to convert `{}` to `i32`.", 2 * t + mi - mdashi)
514            }),
515        );
516
517        if t % 2 == 0 {
518            acc + (num / den) * trigfactor
519        } else {
520            acc - (num / den) * trigfactor
521        }
522    });
523
524    prefactor * d
525}
526
527/// Returns the Wigner rotation matrix in the Euler-angle parametrisation for any integral or
528/// half-integral $`j`$ whose elements are defined by
529/// ```math
530/// \hat{R}(\alpha, \beta, \gamma) \ket{jm}
531/// = \sum_{m'} \ket{jm'} D^{(j)}_{m'm}(\alpha, \beta, \gamma)
532/// ```
533/// and given in [`dmat_euler_gen_element`], where $`-\pi \le \alpha \le \pi`$,
534/// $`0 \le \beta \le \pi`$, $`-\pi \le \gamma \le \pi`$.
535///
536/// # Arguments
537///
538/// * `twoj` - Two times the angular momentum $`2j`$. If this is even, $`j`$ is integral; otherwise,
539/// $`j`$ is half-integral.
540/// * `euler_angles` - A triplet of Euler angles $`(\alpha, \beta, \gamma)`$ in radians, following
541/// the Whitaker convention, *i.e.* $`z_2-y-z_1`$ (extrinsic rotations).
542/// * `increasingm` - If `true`, the rows and columns of $`\mathbf{D}^{(j)}`$ are
543/// arranged in increasing order of $`m_l = -l, \ldots, l`$. If `false`, the order is reversed:
544/// $`m_l = l, \ldots, -l`$. The recommended default is `false`, in accordance with convention.
545///
546/// # Returns
547///
548/// The matrix $`\mathbf{D}^{(j)}(\alpha, \beta, \gamma)`$.
549#[must_use]
550pub fn dmat_euler_gen(
551    twoj: u32,
552    euler_angles: (f64, f64, f64),
553    increasingm: bool,
554) -> Array2<Complex<f64>> {
555    let dim = twoj as usize + 1;
556    let mut dmat = Array2::<Complex<f64>>::zeros((dim, dim));
557    for mdashi in 0..dim {
558        for mi in 0..dim {
559            dmat[(mdashi, mi)] = dmat_euler_gen_element(twoj, mdashi, mi, euler_angles);
560        }
561    }
562    if !increasingm {
563        dmat.invert_axis(Axis(0));
564        dmat.invert_axis(Axis(1));
565    }
566    dmat
567}
568
569/// Returns the Wigner rotation matrix in the angle-axis parametrisation for any integral or
570/// half-integral $`j`$  whose elements are defined by
571///
572/// ```math
573/// \hat{R}(\phi\hat{\mathbf{n}}) \ket{jm}
574/// = \sum_{m'} \ket{jm'} D^{(j)}_{m'm}(\phi\hat{\mathbf{n}}).
575/// ```
576///
577/// # Arguments
578///
579/// * `twoj` - Two times the angular momentum $`2j`$. If this is even, $`j`$ is integral; otherwise,
580/// $`j`$ is half-integral.
581/// * `angle` - The angle $`\phi`$ of the rotation in radians. A positive rotation is an
582/// anticlockwise rotation when looking down `axis`.
583/// * `axis` - A space-fixed vector defining the axis of rotation. The supplied vector will be
584/// normalised.
585/// * `increasingm` - If `true`, the rows and columns of $`\mathbf{D}^{(1/2)}`$ are
586/// arranged in increasing order of $`m_l = -l, \ldots, l`$. If `false`, the order is reversed:
587/// $`m_l = l, \ldots, -l`$. The recommended default is `false`, in accordance with convention.
588///
589/// # Returns
590///
591/// The matrix $`\mathbf{D}^{(j)}(\phi\hat{\mathbf{n}})`$.
592#[must_use]
593pub fn dmat_angleaxis_gen_double(
594    twoj: u32,
595    angle: f64,
596    axis: Vector3<f64>,
597    increasingm: bool,
598) -> Array2<Complex<f64>> {
599    let euler_angles = angleaxis_to_euler_double(angle, axis);
600    dmat_euler_gen(twoj, euler_angles, increasingm)
601}
602
603/// Converts an angle and axis of rotation to Euler angles using the equations in Section
604/// (**3**-5.4) in Altmann, S. L. Rotations, Quaternions, and Double Groups. (Dover
605/// Publications, Inc., 2005), but with an extended range,
606///
607/// ```math
608/// 0 \le \alpha \le 2\pi, \quad
609/// 0 \le \beta \le \pi, \quad
610/// 0 \le \gamma \le 4\pi,
611/// ```
612///
613/// such that all angle-axis parametrisations of $`\phi\hat{\mathbf{n}}`$ for
614/// $`0 \le \phi \le 4 \pi`$ are mapped to unique triplets of $`(\alpha, \beta, \gamma)`$,
615/// as explained in Fan, P.-D., Chen, J.-Q., Mcaven, L. & Butler, P. Unique Euler angles and
616/// self-consistent multiplication tables for double point groups. *International Journal of
617/// Quantum Chemistry* **75**, 1–9 (1999),
618/// [DOI](https://doi.org/10.1002/(SICI)1097-461X(1999)75:1<1::AID-QUA1>3.0.CO;2-V).
619///
620/// When $`\beta = 0`$, only the sum $`\alpha+\gamma`$ is determined. Likewise, when
621/// $`\beta = \pi`$, only the difference $`\alpha-\gamma`$ is determined. We thus set
622/// $`\alpha = 0`$ in these cases and solve for $`\gamma`$ without changing the nature of the
623/// results.
624///
625/// # Arguments
626///
627/// * `angle` - The angle $`\phi`$ of the rotation in radians. A positive rotation is an
628/// anticlockwise rotation when looking down `axis`.
629/// * `axis` - A space-fixed vector defining the axis of rotation $`\hat{\mathbf{n}}`$. The supplied
630/// vector will be normalised.
631///
632/// # Returns
633///
634/// The tuple containing the Euler angles $`(\alpha, \beta, \gamma)`$ in radians, following the
635/// Whitaker convention.
636fn angleaxis_to_euler_double(angle: f64, axis: Vector3<f64>) -> (f64, f64, f64) {
637    let normalised_axis = axis.normalize();
638    let nx = normalised_axis.x;
639    let ny = normalised_axis.y;
640    let nz = normalised_axis.z;
641
642    let double_angle = angle.rem_euclid(4.0 * std::f64::consts::PI);
643    let basic_angle = angle.rem_euclid(2.0 * std::f64::consts::PI);
644    let double = approx::relative_eq!(
645        angle.div_euclid(2.0 * std::f64::consts::PI).rem_euclid(2.0),
646        1.0,
647        epsilon = 1e-14,
648        max_relative = 1e-14
649    );
650
651    let cosbeta = 1.0 - 2.0 * (nx.powi(2) + ny.powi(2)) * (basic_angle / 2.0).sin().powi(2);
652    let cosbeta = if cosbeta.abs() > 1.0 {
653        // Numerical errors can cause cosbeta to be outside [-1, 1].
654        approx::assert_relative_eq!(cosbeta.abs(), 1.0, epsilon = 1e-14, max_relative = 1e-14);
655        cosbeta.round()
656    } else {
657        cosbeta
658    };
659
660    // acos gives 0 <= beta <= pi.
661    let beta = cosbeta.acos();
662
663    let (alpha, gamma) =
664        if approx::relative_ne!(cosbeta.abs(), 1.0, epsilon = 1e-14, max_relative = 1e-14) {
665            // cosbeta != 1 or -1, beta != 0 or pi
666            // alpha and gamma are given by Equations (**3**-5.4) to (**3**-5.10)
667            // in Altmann, S. L. Rotations, Quaternions, and Double Groups. (Dover Publications,
668            // Inc., 2005).
669            // These equations yield the same alpha and gamma for phi and phi+2pi.
670            // We therefore account for double-group behaviours separately.
671            let num_alpha =
672                -nx * basic_angle.sin() + 2.0 * ny * nz * (basic_angle / 2.0).sin().powi(2);
673            let den_alpha =
674                ny * basic_angle.sin() + 2.0 * nx * nz * (basic_angle / 2.0).sin().powi(2);
675            let alpha = num_alpha
676                .atan2(den_alpha)
677                .rem_euclid(2.0 * std::f64::consts::PI);
678
679            let num_gamma =
680                nx * basic_angle.sin() + 2.0 * ny * nz * (basic_angle / 2.0).sin().powi(2);
681            let den_gamma =
682                ny * basic_angle.sin() - 2.0 * nx * nz * (basic_angle / 2.0).sin().powi(2);
683            let gamma_raw = num_gamma.atan2(den_gamma);
684            let gamma = if double {
685                (gamma_raw + 2.0 * std::f64::consts::PI).rem_euclid(4.0 * std::f64::consts::PI)
686            } else {
687                gamma_raw.rem_euclid(4.0 * std::f64::consts::PI)
688            };
689
690            (alpha, gamma)
691        } else if approx::relative_eq!(cosbeta, 1.0, epsilon = 1e-14, max_relative = 1e-14) {
692            // cosbeta == 1, beta == 0
693            // cos(0.5(alpha+gamma)) = cos(0.5phi)
694            // We set alpha == 0 by convention.
695            // We then set gamma = phi mod (4*pi).
696            (0.0, double_angle)
697        } else {
698            // cosbeta == -1, beta == pi
699            // sin(0.5phi) must be non-zero, otherwise cosbeta == 1, a
700            // contradiction.
701            // sin(0.5(alpha-gamma)) = -nx*sin(0.5phi)
702            // cos(0.5(alpha-gamma)) = +ny*sin(0.5phi)
703            // We set alpha == 0 by convention.
704            // gamma then lies in [-2pi, 2pi].
705            // We obtain the same gamma for phi and phi+2pi.
706            // We therefore account for double-group behaviours separately.
707            let gamma_raw = 2.0 * nx.atan2(ny);
708            let gamma = if double {
709                (gamma_raw + 2.0 * std::f64::consts::PI).rem_euclid(4.0 * std::f64::consts::PI)
710            } else {
711                gamma_raw.rem_euclid(4.0 * std::f64::consts::PI)
712            };
713
714            (0.0, gamma)
715        };
716
717    (alpha, beta, gamma)
718}
719
720/// Returns the Wigner rotation matrix in the angle-axis parametrisation for any integral or
721/// half-integral $`j`$  whose elements are defined by
722///
723/// ```math
724/// \hat{R}(\phi\hat{\mathbf{n}}) \ket{jm}
725/// = \sum_{m'} \ket{jm'} D^{(j)}_{m'm}(\phi\hat{\mathbf{n}}),
726/// ```
727///
728/// where the angle of rotation is ensured to be in the range $`[-\pi, \pi]`$. In other words, for
729/// half-odd-integer $`j`$, this function only returns Wigner rotation matrices corresponding to
730/// three-dimensional rotations connected to the identity via a homotopy path of class 0.
731///
732/// # Arguments
733///
734/// * `twoj` - Two times the angular momentum $`2j`$. If this is even, $`j`$ is integral; otherwise,
735/// $`j`$ is half-integral.
736/// * `angle` - The angle $`\phi`$ of the rotation in radians. A positive rotation is an
737/// anticlockwise rotation when looking down `axis`.
738/// * `axis` - A space-fixed vector defining the axis of rotation. The supplied vector will be
739/// normalised.
740/// * `increasingm` - If `true`, the rows and columns of $`\mathbf{D}^{(1/2)}`$ are
741/// arranged in increasing order of $`m_l = -l, \ldots, l`$. If `false`, the order is reversed:
742/// $`m_l = l, \ldots, -l`$. The recommended default is `false`, in accordance with convention.
743///
744/// # Returns
745///
746/// The matrix $`\mathbf{D}^{(j)}(\phi\hat{\mathbf{n}})`$.
747#[must_use]
748pub fn dmat_angleaxis_gen_single(
749    twoj: u32,
750    angle: f64,
751    axis: Vector3<f64>,
752    increasingm: bool,
753) -> Array2<Complex<f64>> {
754    let euler_angles = angleaxis_to_euler_single(angle, axis, 1e-14);
755    let mat = dmat_euler_gen(twoj, euler_angles, increasingm);
756    mat
757}
758
759/// Converts an angle and axis of rotation to Euler angles using the equations in Section
760/// (**3**-5.4) in Altmann, S. L. Rotations, Quaternions, and Double Groups. (Dover
761/// Publications, Inc., 2005) such that
762///
763/// ```math
764/// -\pi \le \alpha \le \pi, \quad
765/// 0 \le \beta \le \pi, \quad
766/// -\pi \le \gamma \le \pi.
767/// ```
768///
769/// When $`\beta = 0`$, only the sum $`\alpha+\gamma`$ is determined. Likewise, when
770/// $`\beta = \pi`$, only the difference $`\alpha-\gamma`$ is determined. We thus set
771/// $`\alpha = 0`$ in these cases and solve for $`\gamma`$ without changing the nature of the
772/// results.
773///
774/// However, it is sometimes necessary to modify the above ranges such that the Euler angles
775/// correspond to an $`\mathsf{SU}(2)`$ rotation of homotopy class 0. This is achieved by ensuring
776/// the scalar part of the corresponding quaternion,
777///
778/// ```math
779/// \lambda = \cos\frac{\beta}{2} \sin\frac{\alpha + \gamma}{2},
780/// ```math
781///
782/// is non-negative by adding or subtracting $`2\pi`$ from $`\alpha`$ as necessary, so that the
783/// modified ranges are
784///
785/// ```math
786/// -2\pi \le \alpha \le 2\pi, \quad
787/// 0 \le \beta \le \pi, \quad
788/// -\pi \le \gamma \le \pi.
789/// ```
790///
791/// # Arguments
792///
793/// * `angle` - The angle $`\phi`$ of the rotation in radians. A positive rotation is an
794/// anticlockwise rotation when looking down `axis`.
795/// * `axis` - A space-fixed vector defining the axis of rotation $`\hat{\mathbf{n}}`$. The supplied
796/// vector will be normalised.
797///
798/// # Returns
799///
800/// The tuple containing the Euler angles $`(\alpha, \beta, \gamma)`$ in radians, following the
801/// Whitaker convention.
802fn angleaxis_to_euler_single(angle: f64, axis: Vector3<f64>, thresh: f64) -> (f64, f64, f64) {
803    let normalised_axis = axis.normalize();
804    let nx = normalised_axis.x;
805    let ny = normalised_axis.y;
806    let nz = normalised_axis.z;
807
808    let (normalised_angle, _) = normalise_rotation_angle(angle, thresh);
809
810    let cosbeta = 1.0 - 2.0 * (nx.powi(2) + ny.powi(2)) * (normalised_angle / 2.0).sin().powi(2);
811    let cosbeta = if cosbeta.abs() > 1.0 {
812        // Numerical errors can cause cosbeta to be outside [-1, 1].
813        approx::assert_relative_eq!(cosbeta.abs(), 1.0, epsilon = thresh, max_relative = thresh);
814        cosbeta.round()
815    } else {
816        cosbeta
817    };
818
819    // 0 <= beta <= pi.
820    let beta = cosbeta.acos();
821
822    let (alpha, gamma) =
823        if approx::relative_ne!(cosbeta.abs(), 1.0, epsilon = thresh, max_relative = thresh) {
824            // cosbeta != 1 or -1, beta != 0 or pi, sin(beta) != 0
825            // alpha and gamma are given by Equations (**3**-5.4) to (**3**-5.10)
826            // in Altmann, S. L. Rotations, Quaternions, and Double Groups. (Dover Publications,
827            // Inc., 2005).
828            let sgn = beta.sin().signum();
829            let num_alpha = -nx * normalised_angle.sin()
830                + 2.0 * ny * nz * (normalised_angle / 2.0).sin().powi(2) * sgn;
831            let den_alpha = ny * normalised_angle.sin()
832                + 2.0 * nx * nz * (normalised_angle / 2.0).sin().powi(2) * sgn;
833            // -pi <= alpha <= pi
834            let alpha = num_alpha.atan2(den_alpha);
835
836            let num_gamma = nx * normalised_angle.sin()
837                + 2.0 * ny * nz * (normalised_angle / 2.0).sin().powi(2) * sgn;
838            let den_gamma = ny * normalised_angle.sin()
839                - 2.0 * nx * nz * (normalised_angle / 2.0).sin().powi(2) * sgn;
840            // -pi <= gamma <= pi
841            let gamma = num_gamma.atan2(den_gamma);
842
843            // We shall next determine the scalar part of the quaternion corresponding to the deduced
844            // Euler angles. We need this to ensure that the Euler angles correspond to rotations of
845            // homotopy class 0.
846            // See Eqn. 12-11.4 in Altmann, S. L. Rotations, Quaternions, and Double Groups. (Dover
847            // Publications, Inc., 2005).
848            let lambda = (beta * 0.5).cos() * ((alpha + gamma) * 0.5).cos();
849            if lambda < 0.0 {
850                if approx::relative_eq!(
851                    gamma.abs(),
852                    std::f64::consts::PI,
853                    epsilon = thresh,
854                    max_relative = thresh
855                ) {
856                    (alpha, gamma - gamma.signum() * 2.0 * std::f64::consts::PI)
857                } else if approx::relative_eq!(
858                    alpha.abs(),
859                    std::f64::consts::PI,
860                    epsilon = thresh,
861                    max_relative = thresh
862                ) {
863                    (alpha - alpha.signum() * 2.0 * std::f64::consts::PI, gamma)
864                } else {
865                    (alpha - alpha.signum() * 2.0 * std::f64::consts::PI, gamma)
866                    // panic!("Unable to adjust α = {alpha:+.7} or γ = {gamma:+.7} within [-π, π] so that the resultant Euler angles correspond to a quaternion with a non-negative scalar part (λ = {lambda:+.7}).")
867                }
868            } else {
869                (alpha, gamma)
870            }
871        } else if approx::relative_eq!(cosbeta, 1.0, epsilon = thresh, max_relative = thresh) {
872            // cosbeta == 1, beta == 0
873            // => cos(0.5(alpha+gamma)) = cos(0.5phi)
874            // => alpha + gamma = ±phi
875            // There are two possibilities:
876            //   1. phi == 0, in which case we simply set both alpha and gamma to zero
877            //   2. nx^2 + ny^2 = 0, which means nx = ny = 0, which means we are at either the
878            //      North pole or the South pole. In either case, we set alpha == 0 by convention.
879            //      We then set gamma = phi at the North pole or -phi at the South pole.
880            // In both cases, lambda = cos(±phi/2), and since -π ≤ phi ≤ π, lambda ≥ 0.
881            if approx::relative_eq!(
882                normalised_angle,
883                0.0,
884                epsilon = thresh,
885                max_relative = thresh
886            ) {
887                (0.0, normalised_angle)
888            } else if nz > 0.0 {
889                (0.0, normalised_angle)
890            } else {
891                (0.0, -normalised_angle)
892            }
893        } else {
894            // cosbeta == -1, beta == pi
895            // sin(0.5phi) must be non-zero, otherwise cosbeta == 1, a contradiction.
896            // sin(0.5(alpha-gamma)) = -nx*sin(0.5phi)
897            // cos(0.5(alpha-gamma)) = +ny*sin(0.5phi)
898            // We set alpha == 0 by convention.
899            // gamma then lies in [-pi, pi].
900            // lambda = 0 since beta = pi.
901            let gamma = 2.0 * nx.atan2(ny);
902
903            (0.0, gamma)
904        };
905
906    (alpha, beta, gamma)
907}