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}