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}