qsym2/angmom/sh_rotation_3d.rs
1//! Three-dimensional rotations of spherical harmonics.
2
3use std::cmp::Ordering;
4
5use approx;
6use nalgebra::{Rotation3, Unit, Vector3};
7use ndarray::{Array2, Axis, ShapeBuilder};
8use num_traits::ToPrimitive;
9
10#[cfg(test)]
11#[path = "sh_rotation_3d_tests.rs"]
12mod sh_rotation_3d_tests;
13
14/// Returns the generalised Kronecker delta $`\delta_{ij}`$ for any $`i`$ and $`j`$ that have a
15/// partial equivalence relation.
16///
17/// # Returns
18///
19/// `0` if $`i \ne j`$, `1` if $`i = j`$.
20fn kdelta<T: PartialEq>(i: &T, j: &T) -> u8 {
21 u8::from(i == j)
22}
23
24/// Returns the function $`_iP^l_{\mu m'}`$ as defined in Table 2 of Ivanic, J. & Ruedenberg, K.
25/// Rotation Matrices for Real Spherical Harmonics. Direct Determination by Recursion.
26/// *The Journal of Physical Chemistry* **100**, 9099–9100 (1996),
27/// [DOI](https://doi.org/10.1021/jp953350u).
28///
29/// # Arguments
30///
31/// * `i` - The index $`i`$ satisfying $`-1 \le i \le 1`$.
32/// * `l` - The spherical harmonic order $`l \ge 2`$.
33/// * `mu` - The index $`\mu`$ satisfying $`-l+1 \le \mu \le l-1`$.
34/// * `mdash` - The index $`m'`$ satisfying $`-l \le m' \le l`$.
35/// * `rmat` - The representation matrix of the transformation of interest in the basis of
36/// coordinate *functions* $`(y, z, x)`$, which are isosymmetric to the real spherical harmonics
37/// $`(Y_{1, -1}, Y_{1, 0}, Y_{1, 1})`$.
38/// * `rlm1` - The representation matrix of the transformation of interest in the basis of the real
39/// spherical harmonics $`Y_{l-1, m}`$ ordered by increasing $`m`$.
40///
41/// # Returns
42///
43/// The value of $`_iP^l_{\mu m'}`$.
44fn func_p(i: i8, l: u32, mu: i64, mdash: i64, rmat: &Array2<f64>, rlm1: &Array2<f64>) -> f64 {
45 assert!(i.abs() <= 1, "`i` must be between -1 and 1 (inclusive).");
46 assert!(l >= 2, "`l` must be at least 2.");
47 let li64 = i64::from(l);
48 assert!(
49 mu.abs() < li64,
50 "Index `mu` = {mu} lies outside [{}, {}].",
51 -(li64) + 1,
52 l - 1
53 );
54 assert!(
55 mdash.abs() <= li64,
56 "Index `mdash` = {mdash} lies outside [-{l}, {l}].",
57 );
58 assert_eq!(rmat.shape(), &[3, 3], "`rmat` must be a 3 × 3 matrix.");
59 assert_eq!(
60 rlm1.shape(),
61 &[2 * l as usize - 1, 2 * l as usize - 1],
62 "`rlm1` must be a {} × {} matrix.",
63 2 * l as usize - 1,
64 2 * l as usize - 1
65 );
66
67 let ii = usize::try_from(i + 1).expect("Unable to convert `i + 1` to `usize`.");
68 let mui =
69 usize::try_from(mu + (li64 - 1)).expect("Unable to convert `mu + (l - 1)` to `usize`.");
70 let mdashi = usize::try_from(mdash + li64).expect("Unable to convert `mdash + l` to `usize`.");
71 let lusize = usize::try_from(l).expect("Unable to convert `l` to `usize`.");
72 if mdash == li64 {
73 // Easier-to-read expression:
74 // R[i + 1, 1 + 1] * Rlm1[mu + (l - 1), l - 1 + (l - 1)]
75 // - R[i + 1, -1 + 1] * Rlm1[mu + (l - 1), -l + 1 + (l - 1)]
76 rmat[(ii, 2)] * rlm1[(mui, 2 * (lusize - 1))] - rmat[(ii, 0)] * rlm1[(mui, 0)]
77 } else if mdash == -li64 {
78 // Easier-to-read expression:
79 // R[i + 1, 1 + 1] * Rlm1[mu + (l - 1), -l + 1 + (l - 1)]
80 // + R[i + 1, -1 + 1] * Rlm1[mu + (l - 1), l - 1 + (l - 1)]
81 rmat[(ii, 2)] * rlm1[(mui, 0)] + rmat[(ii, 0)] * rlm1[(mui, 2 * (lusize - 1))]
82 } else {
83 // Easier-to-read expression:
84 // R[i + 1, 0 + 1] * Rlm1[mu + (l - 1), mdash + (l - 1)]
85 rmat[(ii, 1)] * rlm1[(mui, mdashi - 1)]
86 }
87}
88
89/// Returns the function $`U^l_{mm'}`$ as defined in Table 2 of Ivanic, J. & Ruedenberg, K.
90/// Rotation Matrices for Real Spherical Harmonics. Direct Determination by Recursion.
91/// *The Journal of Physical Chemistry* **100**, 9099–9100 (1996),
92/// [DOI](https://doi.org/10.1021/jp953350u).
93///
94/// # Arguments
95///
96/// * `l` - The spherical harmonic order $`l \ge 2`$.
97/// * `m` - The index $`m`$ satisfying $`-l+1 \le \mu \le l-1`$.
98/// * `mdash` - The index $`m'`$ satisfying $`-l \le m' \le l`$.
99/// * `rmat` - The representation matrix of the transformation of interest in the basis of
100/// coordinate *functions* $`(y, z, x)`$, which are isosymmetric to the real spherical harmonics
101/// $`(Y_{1, -1}, Y_{1, 0}, Y_{1, 1})`$.
102/// * rlm1 - The representation matrix of the transformation of interest in the basis of the real
103/// spherical harmonics $`Y_{l-1, m}`$ ordered by increasing $`m`$.
104///
105/// # Returns
106///
107/// The value of $`U^l_{mm'}`$.
108fn func_u(l: u32, m: i64, mdash: i64, rmat: &Array2<f64>, rlm1: &Array2<f64>) -> f64 {
109 func_p(0, l, m, mdash, rmat, rlm1)
110}
111
112/// Returns the function $`V^l_{mm'}`$ as defined in Table 2 of Ivanic, J. & Ruedenberg, K.
113/// Rotation Matrices for Real Spherical Harmonics. Direct Determination by Recursion.
114/// *The Journal of Physical Chemistry* **100**, 9099–9100 (1996),
115/// [DOI](https://doi.org/10.1021/jp953350u).
116///
117/// # Arguments
118///
119/// * `l` - The spherical harmonic order $`l \ge 2`$.
120/// * `m` - The index $`m`$ satisfying $`-l \le m \le l`$.
121/// * `mdash` - The index $`m'`$ satisfying $`-l \le m' \le l`$.
122/// * `rmat` - The representation matrix of the transformation of interest in the basis of
123/// coordinate *functions* $`(y, z, x)`$, which are isosymmetric to the real spherical harmonics
124/// $`(Y_{1, -1}, Y_{1, 0}, Y_{1, 1})`$.
125/// * `rlm1` - The representation matrix of the transformation of interest in the basis of the real
126/// spherical harmonics $`Y_{l-1, m}`$ ordered by increasing $`m`$.
127///
128/// # Returns
129///
130/// The value of $`V^l_{mm'}`$.
131fn func_v(l: u32, m: i64, mdash: i64, rmat: &Array2<f64>, rlm1: &Array2<f64>) -> f64 {
132 assert!(l >= 2, "`l` must be at least 2.");
133 let li64 = i64::from(l);
134 assert!(m.abs() <= li64, "Index `m` = {m} lies outside [-{l}, {l}].",);
135 assert!(
136 mdash.abs() <= li64,
137 "Index `mdash` = {mdash} lies outside [-{l}, {l}].",
138 );
139 assert_eq!(rmat.shape(), &[3, 3], "`rmat` must be a 3 × 3 matrix.");
140 assert_eq!(
141 rlm1.shape(),
142 &[2 * l as usize - 1, 2 * l as usize - 1],
143 "`rlm1` must be a {} × {} matrix.",
144 2 * l as usize - 1,
145 2 * l as usize - 1
146 );
147
148 match m.cmp(&0) {
149 Ordering::Greater => {
150 func_p(1, l, m - 1, mdash, rmat, rlm1) * (f64::from(1 + kdelta(&m, &1))).sqrt()
151 - func_p(-1, l, -m + 1, mdash, rmat, rlm1) * (f64::from(1 - kdelta(&m, &1)))
152 }
153 Ordering::Less => {
154 func_p(1, l, m + 1, mdash, rmat, rlm1) * (f64::from(1 - kdelta(&m, &(-1))))
155 + func_p(-1, l, -m - 1, mdash, rmat, rlm1)
156 * (f64::from(1 + kdelta(&m, &(-1)))).sqrt()
157 }
158 Ordering::Equal => {
159 func_p(1, l, 1, mdash, rmat, rlm1) + func_p(-1, l, -1, mdash, rmat, rlm1)
160 }
161 }
162}
163
164/// Returns the function $`W^l_{mm'}`$ as defined in Table 2 of Ivanic, J. & Ruedenberg, K.
165/// Rotation Matrices for Real Spherical Harmonics. Direct Determination by Recursion.
166/// *The Journal of Physical Chemistry* **100**, 9099–9100 (1996),
167/// [DOI](https://doi.org/10.1021/jp953350u).
168///
169/// # Arguments
170///
171/// * `l` - The spherical harmonic order $`l \ge 2`$.
172/// * `m` - The index $`m`$ satisfying $`-l+2 \le m \le l-2`$ and $`m \ne 0`$.
173/// * `mdash` - The index $`m'`$ satisfying $`-l \le m' \le l`$.
174/// * `rmat` - The representation matrix of the transformation of interest in the basis of
175/// coordinate *functions* $`(y, z, x)`$, which are isosymmetric to the real spherical harmonics
176/// $`(Y_{1, -1}, Y_{1, 0}, Y_{1, 1})`$.
177/// * `rlm1` - The representation matrix of the transformation of interest in the basis of the real
178/// spherical harmonics $`Y_{l-1, m}`$ ordered by increasing $`m`$.
179///
180/// # Returns
181///
182/// The value of $`W^l_{mm'}`$.
183fn func_w(l: u32, m: i64, mdash: i64, rmat: &Array2<f64>, rlm1: &Array2<f64>) -> f64 {
184 assert!(l >= 2, "`l` must be at least 2.");
185 let li64 = i64::from(l);
186 assert!(
187 m.abs() <= li64 - 2,
188 "Index `m` = {m} lies outside [{}, {}].",
189 -li64 + 2,
190 li64 - 2
191 );
192 assert_ne!(m, 0, "`m` cannot be zero.");
193 assert!(
194 mdash.abs() <= li64,
195 "Index `mdash` = {mdash} lies outside [-{l}, {l}].",
196 );
197 assert_eq!(rmat.shape(), &[3, 3], "`rmat` must be a 3 × 3 matrix.");
198 assert_eq!(
199 rlm1.shape(),
200 &[2 * l as usize - 1, 2 * l as usize - 1],
201 "`rlm1` must be a {} × {} matrix.",
202 2 * l as usize - 1,
203 2 * l as usize - 1
204 );
205 if m > 0 {
206 func_p(1, l, m + 1, mdash, rmat, rlm1) + func_p(-1, l, -m - 1, mdash, rmat, rlm1)
207 } else {
208 func_p(1, l, m - 1, mdash, rmat, rlm1) - func_p(-1, l, -m + 1, mdash, rmat, rlm1)
209 }
210}
211
212/// Returns the coefficient $`u^l_{mm'}`$ as defined in Table 1 of Ivanic, J. & Ruedenberg, K.
213/// Rotation Matrices for Real Spherical Harmonics. Direct Determination by Recursion.
214/// *The Journal of Physical Chemistry* **100**, 9099–9100 (1996),
215/// [DOI](https://doi.org/10.1021/jp953350u).
216///
217/// # Arguments
218///
219/// * `l` - The spherical harmonic order $`l`$.
220/// * `m` - The index $`m`$ satisfying $`-l \le m \le l`$.
221/// * `mdash` - The index $`m'`$ satisfying $`-l \le m' \le l`$.
222///
223/// # Returns
224///
225/// The value of $`u^l_{mm'}`$.
226fn coeff_u(l: u32, m: i64, mdash: i64) -> f64 {
227 let li64 = i64::from(l);
228 assert!(m.abs() <= li64, "Index `m` = {m} lies outside [-{l}, {l}].",);
229 assert!(
230 mdash.abs() <= li64,
231 "Index `mdash` = {mdash} lies outside [-{l}, {l}].",
232 );
233
234 let num = (li64 + m) * (li64 - m);
235 let den = if mdash.abs() < li64 {
236 (li64 + mdash) * (li64 - mdash)
237 } else {
238 (2 * li64) * (2 * li64 - 1)
239 };
240 (num.to_f64()
241 .unwrap_or_else(|| panic!("Unable to convert `{num}` to `f64`."))
242 / den
243 .to_f64()
244 .unwrap_or_else(|| panic!("Unable to convert `{den}` to `f64`.")))
245 .sqrt()
246}
247
248/// Returns the coefficient $`v^l_{mm'}`$ as defined in Table 1 of Ivanic, J. & Ruedenberg, K.
249/// Rotation Matrices for Real Spherical Harmonics. Direct Determination by Recursion.
250/// *The Journal of Physical Chemistry* **100**, 9099–9100 (1996),
251/// [DOI](https://doi.org/10.1021/jp953350u).
252///
253/// # Arguments
254///
255/// * `l` - The spherical harmonic order $`l`$.
256/// * `m` - The index $`m`$ satisfying $`-l \le m \le l`$.
257/// * `mdash` - The index $`m'`$ satisfying $`-l \le m' \le l`$.
258///
259/// # Returns
260///
261/// The value of $`v^l_{mm'}`$.
262fn coeff_v(l: u32, m: i64, mdash: i64) -> f64 {
263 let li64 = i64::from(l);
264 assert!(m.abs() <= li64, "Index `m` = {m} lies outside [-{l}, {l}].",);
265 assert!(
266 mdash.abs() <= li64,
267 "Index `mdash` = {mdash} lies outside [-{l}, {l}].",
268 );
269
270 let num = (1 + i64::from(kdelta(&m, &0))) * (li64 + m.abs() - 1) * (li64 + m.abs());
271 let den = if mdash.abs() < li64 {
272 (li64 + mdash) * (li64 - mdash)
273 } else {
274 (2 * li64) * (2 * li64 - 1)
275 };
276 0.5 * (num
277 .to_f64()
278 .unwrap_or_else(|| panic!("Unable to convert `{num}` to `f64`."))
279 / den
280 .to_f64()
281 .unwrap_or_else(|| panic!("Unable to convert `{den}` to `f64`.")))
282 .sqrt()
283 * f64::from(1 - 2 * i16::from(kdelta(&m, &0)))
284}
285
286/// Returns the coefficient $`w^l_{mm'}`$ as defined in Table 1 of Ivanic, J. & Ruedenberg, K.
287/// Rotation Matrices for Real Spherical Harmonics. Direct Determination by Recursion.
288/// *The Journal of Physical Chemistry* **100**, 9099–9100 (1996),
289/// [DOI](https://doi.org/10.1021/jp953350u).
290///
291/// # Arguments
292///
293/// * `l` - The spherical harmonic order $`l`$.
294/// * `m` - The index $`m`$ satisfying $`-l \le m \le l`$.
295/// * `mdash` - The index $`m'`$ satisfying $`-l \le m' \le l`$.
296///
297/// # Returns
298///
299/// The value of $`w^l_{mm'}`$.
300fn coeff_w(l: u32, m: i64, mdash: i64) -> f64 {
301 let li64 = i64::from(l);
302 assert!(m.abs() <= li64, "Index `m` = {m} lies outside [-{l}, {l}].",);
303 assert!(
304 mdash.abs() <= li64,
305 "Index `mdash` = {mdash} lies outside [-{l}, {l}].",
306 );
307
308 let num = (li64 - m.abs() - 1) * (li64 - m.abs());
309 let den = if mdash.abs() < li64 {
310 (li64 + mdash) * (li64 - mdash)
311 } else {
312 (2 * li64) * (2 * li64 - 1)
313 };
314 -0.5 * ((num
315 .to_f64()
316 .unwrap_or_else(|| panic!("Unable to convert `{num}` to `f64`.")))
317 / (den
318 .to_f64()
319 .unwrap_or_else(|| panic!("Unable to convert `{den}` to `f64`."))))
320 .sqrt()
321 * f64::from(1 - kdelta(&m, &0))
322}
323
324/// Returns the representation matrix $`\mathbf{R}`$ for a rotation in the basis of the coordinate
325/// *functions* $`(y, z, x)`$.
326///
327/// Let $`\hat{R}(\phi, \hat{\mathbf{n}})`$ be a rotation parametrised by the angle $`\phi`$ and
328/// axis $`\hat{\mathbf{n}}`$. The corresponding representation matrix
329/// $`\mathbf{R}(\phi, \hat{\mathbf{n}})`$ is defined as
330///
331/// ```math
332/// \hat{R}(\phi, \hat{\mathbf{n}})\ (y, z, x)
333/// = (y, z, x) \mathbf{R}(\phi, \hat{\mathbf{n}})
334/// ```
335///
336/// See Section **2**-4 of Altmann, S. L. Rotations, Quaternions, and Double Groups. (Dover
337/// Publications, Inc., 2005) for a detailed discussion on how $`(y, z, x)`$ should be considered
338/// as coordinate *functions*.
339///
340/// # Arguments
341///
342/// * `angle` - The angle $`\phi`$ of the rotation in radians. A positive rotation is an
343/// anticlockwise rotation when looking down `axis`.
344/// * `axis` - A space-fixed vector defining the axis of rotation. The supplied vector will be
345/// normalised.
346
347/// # Returns
348///
349/// The representation matrix $`\mathbf{R}(\phi, \hat{\mathbf{n}})`$.
350///
351/// # Panics
352///
353/// Panics when a three-dimensional rotation matrix cannot be constructed for `angle` and `axis`.
354#[must_use]
355pub fn rmat(angle: f64, axis: Vector3<f64>) -> Array2<f64> {
356 if axis.norm() < 1e-14 {
357 Array2::eye(3)
358 } else {
359 let normalised_axis = Unit::new_normalize(axis);
360 let rot = Rotation3::from_axis_angle(&normalised_axis, angle);
361 // nalgebra matrix iter is column-major.
362 let rot_array = Array2::<f64>::from_shape_vec(
363 (3, 3).f(),
364 rot.into_inner().iter().copied().collect::<Vec<_>>(),
365 )
366 .unwrap_or_else(
367 |_| panic!(
368 "Unable to construct a three-dimensional rotation matrix for angle {angle} and axis {axis}."
369 )
370 );
371 rot_array
372 .select(Axis(0), &[1, 2, 0])
373 .select(Axis(1), &[1, 2, 0])
374 }
375}
376
377/// Computes the representation matrix $`\mathbf{R}^l`$ for a transformation of interest in the
378/// basis of real spherical harmonics $`Y_{lm}`$ ordered by increasing $`m`$, as defined in
379/// Equation 5.8 and given in Equation 8.1 of Ivanic, J. & Ruedenberg, K. Rotation Matrices for
380/// Real Spherical Harmonics. Direct Determination by Recursion. *The Journal of Physical
381/// Chemistry* **100**, 9099–9100 (1996), [DOI](https://doi.org/10.1021/jp953350u).
382///
383/// # Arguments
384///
385/// * `l` - The spherical harmonic order $`l \ge 2`$.
386/// * `rmat` - The representation matrix of the transformation of interest in the basis of coordinate
387/// *functions* $`(y, z, x)`$, which are isosymmetric to the real spherical harmonics
388/// $`(Y_{1, -1}, Y_{1, 0}, Y_{1, 1})`$.
389/// * `rlm1` - The representation matrix of the transformation of interest in the basis of the real
390/// spherical harmonics $`Y_{l-1, m}`$ ordered by increasing $`m`$.
391///
392/// # Returns
393///
394/// The required representation matrix $`\mathbf{R}^l`$.
395///
396/// # Panics
397///
398/// Panics when `l` is less than `2`.
399#[must_use]
400pub fn rlmat(l: u32, rmat: &Array2<f64>, rlm1: &Array2<f64>) -> Array2<f64> {
401 assert!(l >= 2, "`l` must be at least 2.");
402 let li64 = i64::from(l);
403 let mut rl = Array2::<f64>::zeros((2 * l as usize + 1, 2 * l as usize + 1));
404 for (mi, m) in (-li64..=li64).enumerate() {
405 for (mdashi, mdash) in (-li64..=li64).enumerate() {
406 let cu = coeff_u(l, m, mdash);
407 let f_u = if approx::relative_ne!(cu.abs(), 0.0, epsilon = 1e-14, max_relative = 1e-14)
408 {
409 func_u(l, m, mdash, rmat, rlm1)
410 } else {
411 0.0
412 };
413
414 let cv = coeff_v(l, m, mdash);
415 let f_v = if approx::relative_ne!(cv.abs(), 0.0, epsilon = 1e-14, max_relative = 1e-14)
416 {
417 func_v(l, m, mdash, rmat, rlm1)
418 } else {
419 0.0
420 };
421
422 let cw = coeff_w(l, m, mdash);
423 let f_w = if approx::relative_ne!(cw.abs(), 0.0, epsilon = 1e-14, max_relative = 1e-14)
424 {
425 func_w(l, m, mdash, rmat, rlm1)
426 } else {
427 0.0
428 };
429
430 rl[(mi, mdashi)] = cu * f_u + cv * f_v + cw * f_w;
431 }
432 }
433 rl
434}