qsym2/angmom/
sh_conversion.rs

1//! Conversions between spherical/solid harmonics and Cartesian functions.
2
3use std::cmp::Ordering;
4
5use factorial::Factorial;
6use ndarray::{Array2, Axis};
7use num::{BigUint, Complex};
8use num_traits::{cast::ToPrimitive, Zero};
9
10use crate::basis::ao::{CartOrder, PureOrder};
11use crate::permutation::PermutableCollection;
12
13#[cfg(test)]
14#[path = "sh_conversion_tests.rs"]
15mod sh_conversion_tests;
16
17/// Calculates the number of combinations of `n` things taken `r` at a time (signed arguments).
18///
19/// If $`n < 0`$ or $`r < 0`$ or $`r > n`$, `0` is returned.
20///
21/// # Arguments
22///
23/// * `n` - Number of things.
24/// * `r` - Number of elements taken.
25///
26/// # Returns
27///
28/// The number of combinations.
29fn comb(n: i32, r: i32) -> BigUint {
30    if n < 0 || r < 0 || r > n {
31        BigUint::zero()
32    } else {
33        let nu = u32::try_from(n).expect("Unable to convert `n` to `u32`.");
34        let ru = u32::try_from(r).expect("Unable to convert `r` to `u32`.");
35        (nu - ru + 1..=nu).product::<BigUint>()
36            / BigUint::from(ru)
37                .checked_factorial()
38                .unwrap_or_else(|| panic!("Unable to compute the factorial of {ru}."))
39    }
40}
41
42/// Calculates the number of combinations of `n` things taken `r` at a time (unsigned arguments).
43///
44/// If $`r > n`$, `0` is returned.
45///
46/// # Arguments
47///
48/// * `n` - Number of things.
49/// * `r` - Number of elements taken.
50///
51/// # Returns
52///
53/// The number of combinations.
54fn combu(nu: u32, ru: u32) -> BigUint {
55    if ru > nu {
56        BigUint::zero()
57    } else {
58        (nu - ru + 1..=nu).product::<BigUint>()
59            / BigUint::from(ru)
60                .checked_factorial()
61                .unwrap_or_else(|| panic!("Unable to compute the factorial of {ru}."))
62    }
63}
64
65/// Calculates the number of permutations of `n` things taken `r` at a time (signed arguments).
66///
67/// If $`n < 0`$ or $`r < 0`$ or $`r > n`$, `0` is returned.
68///
69/// # Arguments
70///
71/// * `n` - Number of things.
72/// * `r` - Number of elements taken.
73///
74/// # Returns
75///
76/// The number of permutations.
77fn perm(n: i32, r: i32) -> BigUint {
78    if n < 0 || r < 0 || r > n {
79        BigUint::zero()
80    } else {
81        let nu = u32::try_from(n).expect("Unable to convert `n` to `u32`.");
82        let ru = u32::try_from(r).expect("Unable to convert `r` to `u32`.");
83        (nu - ru + 1..=nu).product::<BigUint>()
84    }
85}
86
87/// Calculates the number of permutations of `n` things taken `r` at a time (unsigned arguments).
88///
89/// If $`r > n`$, `0` is returned.
90///
91/// # Arguments
92///
93/// * `n` - Number of things.
94/// * `r` - Number of elements taken.
95///
96/// # Returns
97///
98/// The number of permutations.
99fn permu(nu: u32, ru: u32) -> BigUint {
100    if ru > nu {
101        BigUint::zero()
102    } else {
103        (nu - ru + 1..=nu).product::<BigUint>()
104    }
105}
106
107/// Obtains the normalisation constant for a solid harmonic Gaussian, as given in Equation 8 of
108/// Schlegel, H. B. & Frisch, M. J. Transformation between Cartesian and pure spherical harmonic
109/// Gaussians. *International Journal of Quantum Chemistry* **54**, 83–87 (1995),
110/// [DOI](https://doi.org/10.1002/qua.560540202).
111///
112/// The complex solid harmonic Gaussian is defined in Equation 1 of the above reference as
113///
114/// ```math
115/// \tilde{g}(\alpha, l, m, n, \mathbf{r})
116///     = \tilde{N}(n, \alpha) Y_l^m r^n e^{-\alpha r^2},
117/// ```
118///
119/// where $`Y_l^m`$ is a complex spherical harmonic of degree $`l`$ and order $`m`$.
120///
121/// # Arguments
122///
123/// * `n` - The non-negative exponent of the radial part of the solid harmonic Gaussian.
124/// * `alpha` - The coefficient on the exponent of the Gaussian term.
125///
126/// # Returns
127///
128/// The normalisation constant $`\tilde{N}(n, \alpha)`$.
129fn norm_sph_gaussian(n: u32, alpha: f64) -> f64 {
130    let num = (BigUint::from(2u64).pow(2 * n + 3)
131        * BigUint::from(u64::from(n) + 1)
132            .checked_factorial()
133            .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", u64::from(n) + 1)))
134    .to_f64()
135    .expect("Unable to convert a `BigUint` value to `f64`.")
136        * alpha.powf(f64::from(n) + 1.5);
137    let den = BigUint::from(2 * u64::from(n) + 2)
138        .checked_factorial()
139        .unwrap_or_else(|| {
140            panic!(
141                "Unable to compute the factorial of {}.",
142                2 * u64::from(n) + 2
143            )
144        })
145        .to_f64()
146        .expect("Unable to convert a `BigUint` value to `f64`.")
147        * std::f64::consts::PI.sqrt();
148    (num / den).sqrt()
149}
150
151/// Obtains the normalisation constant for a Cartesian Gaussian, as given in Equation 9 of
152/// Schlegel, H. B. & Frisch, M. J. Transformation between Cartesian and pure spherical harmonic
153/// Gaussians. *International Journal of Quantum Chemistry* **54**, 83–87 (1995),
154/// [DOI](https://doi.org/10.1002/qua.560540202).
155///
156/// The Cartesian Gaussian is defined in Equation 2 of the above reference as
157///
158/// ```math
159/// g(\alpha, l_x, l_y, l_z, \mathbf{r})
160///     = N(l_x, l_y, l_z, \alpha) x^{l_x} y^{l_y} z^{l_z} e^{-\alpha r^2}.
161/// ```
162///
163/// # Arguments
164///
165/// * `lcartqns` - A tuple of $`(l_x, l_y, l_z)`$ specifying the non-negative exponents of
166/// the Cartesian components of the Cartesian Gaussian.
167/// * `alpha` - The coefficient on the exponent of the Gaussian term.
168///
169/// # Returns
170///
171/// The normalisation constant $`N(l_x, l_y, l_z, \alpha)`$.
172fn norm_cart_gaussian(lcartqns: (u32, u32, u32), alpha: f64) -> f64 {
173    let (lx, ly, lz) = lcartqns;
174    let lcart = lx + ly + lz;
175    let num = (BigUint::from(2u32).pow(2 * lcart)
176        * BigUint::from(lx)
177            .checked_factorial()
178            .unwrap_or_else(|| panic!("Unable to compute the factorial of {lx}."))
179        * BigUint::from(ly)
180            .checked_factorial()
181            .unwrap_or_else(|| panic!("Unable to compute the factorial of {ly}."))
182        * BigUint::from(lz)
183            .checked_factorial()
184            .unwrap_or_else(|| panic!("Unable to compute the factorial of {lz}.")))
185    .to_f64()
186    .expect("Unable to convert a `BigUint` value to `f64`.")
187        * alpha.powf(f64::from(lcart) + 1.5);
188    let den = (BigUint::from(2 * lx)
189        .checked_factorial()
190        .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", 2 * lx))
191        * BigUint::from(2 * ly)
192            .checked_factorial()
193            .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", 2 * ly))
194        * BigUint::from(2 * lz)
195            .checked_factorial()
196            .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", 2 * lz)))
197    .to_f64()
198    .expect("Unable to convert a `BigUint` value to `f64`.")
199        * std::f64::consts::PI.powi(3).sqrt();
200    (num / den).sqrt()
201}
202
203/// Obtain the complex coefficients $`c(l, m_l, n, l_x, l_y, l_z)`$ based on Equation 15 of
204/// Schlegel, H. B. & Frisch, M. J. Transformation between Cartesian and pure spherical harmonic
205/// Gaussians. *International Journal of Quantum Chemistry* **54**, 83–87 (1995),
206/// [DOI](https://doi.org/10.1002/qua.560540202), but more generalised
207/// for $`l \leq l_{\mathrm{cart}} = l_x + l_y + l_z`$.
208///
209/// Let $`\tilde{g}(\alpha, l, m_l, l_{\mathrm{cart}}, \mathbf{r})`$ be a complex solid
210/// harmonic Gaussian as defined in Equation 1 of
211/// the above reference with $`n = l_{\mathrm{cart}}`$, and let
212/// $`g(\alpha, l_x, l_y, l_z, \mathbf{r})`$ be a Cartesian Gaussian as defined in
213/// Equation 2 of the above reference.
214/// The complex coefficients $`c(l, m_l, n, l_x, l_y, l_z)`$ effect the transformation
215///
216/// ```math
217/// \tilde{g}(\alpha, l, m_l, l_{\mathrm{cart}}, \mathbf{r})
218/// = \sum_{l_x+l_y+l_z=l_{\mathrm{cart}}}
219///     c(l, m_l, l_{\mathrm{cart}}, l_x, l_y, l_z)
220///     g(\alpha, l_x, l_y, l_z, \mathbf{r})
221/// ```
222///
223/// and are given by
224///
225/// ```math
226/// c(l, m_l, l_{\mathrm{cart}}, l_x, l_y, l_z)
227/// = \frac{\tilde{N}(l_{\mathrm{cart}}, \alpha)}{N(l_x, l_y, l_z, \alpha)}
228///     \tilde{c}(l, m_l, l_{\mathrm{cart}}, l_x, l_y, l_z).
229/// ```
230///
231/// The normalisation constants $`\tilde{N}(l_{\mathrm{cart}}, \alpha)`$
232/// and $`N(l_x, l_y, l_z, \alpha)`$ are given in Equations 8 and 9 of
233/// the above reference, and for $`n = l_{\mathrm{cart}}`$, this ratio turns out to be
234/// independent of $`\alpha`$.
235/// The more general form of $`\tilde{c}(l, m_l, l_{\mathrm{cart}}, l_x, l_y, l_z)`$ has been
236/// derived to be
237///
238/// ```math
239/// \tilde{c}(l, m_l, l_{\mathrm{cart}}, l_x, l_y, l_z)
240/// = \frac{\lambda_{\mathrm{cs}}}{2^l l!}
241///     \sqrt{\frac{(2l+1)(l-\lvert m_l \rvert)!}{4\pi(l+\lvert m_l \rvert)!}}
242///     \sum_{i=0}^{(l-\lvert m_l \rvert)/2}
243///         {l\choose i} \frac{(-1)^i(2l-2i)!}{(l-\lvert m_l \rvert -2i)!}\\
244///     \sum_{p=0}^{\lvert m_l \rvert} {{\lvert m_l \rvert} \choose p}
245///         (\pm \mathbb{i})^{\lvert m_l \rvert-p}
246///     \sum_{q=0}^{\Delta l/2} {{\Delta l/2} \choose q} {i \choose j_q}
247///     \sum_{k=0}^{j_q} {q \choose t_{pk}} {j_q \choose k}
248/// ```
249///
250/// where $`+\mathbb{i}`$ applies for $`m_l > 0`$, $`-\mathbb{i}`$
251/// for $`m_l \le 0`$, $`\lambda_{\mathrm{cs}}`$ is the Condon--Shortley
252/// phase given by
253///
254/// ```math
255/// \lambda_{\mathrm{cs}} =
256///     \begin{cases}
257///         (-1)^{m_l} & m_l > 0 \\
258///         1          & m_l \leq 0
259///     \end{cases}
260/// ```
261///
262/// and
263///
264/// ```math
265/// t_{pk} = \frac{l_x-p-2k}{2} \quad \textrm{and} \quad
266/// j_q = \frac{l_x+l_y-\lvert m_l \rvert-2q}{2}.
267/// ```
268///
269/// If $`\Delta l`$ is odd, $`\tilde{c}(l, m_l, l_{\mathrm{cart}}, l_x, l_y, l_z)`$ must vanish.
270/// When  $`t_{pk}`$ or $`j_q`$ is a half-integer, the inner sum in which it is involved
271/// evaluates to zero.
272///
273/// # Arguments
274///
275/// * `lpureqns` - A tuple of $`(l, m_l)`$ specifying the quantum numbers for the spherical
276/// harmonic component of the solid harmonic Gaussian.
277/// * `lcartqns` - A tuple of $`(l_x, l_y, l_z)`$ specifying the exponents of the Cartesian
278/// components of the Cartesian Gaussian.
279/// * `csphase` - If `true`, the Condon--Shortley phase will be used as defined above.
280/// If `false`, this phase will be set to unity.
281///
282/// # Returns
283///
284/// The complex factor $`c(l, m_l, l_{\mathrm{cart}}, l_x, l_y, l_z)`$.
285///
286/// # Panics
287///
288/// Panics when any required factorials cannot be computed.
289#[allow(clippy::too_many_lines)]
290pub fn complexc(lpureqns: (u32, i32), lcartqns: (u32, u32, u32), csphase: bool) -> Complex<f64> {
291    let (l, m) = lpureqns;
292    let li32 = i32::try_from(l).unwrap_or_else(|_| panic!("Cannot convert `{l}` to `i32`."));
293    assert!(
294        m.unsigned_abs() <= l,
295        "m must be between -l and l (inclusive)."
296    );
297    let (lx, ly, lz) = lcartqns;
298    let lxi32 = i32::try_from(lx).unwrap_or_else(|_| panic!("Cannot convert `{lx}` to `i32`."));
299    let lyi32 = i32::try_from(ly).unwrap_or_else(|_| panic!("Cannot convert `{ly}` to `i32`."));
300    let lzi32 = i32::try_from(lz).unwrap_or_else(|_| panic!("Cannot convert `{lz}` to `i32`."));
301    let lcart = lx + ly + lz;
302    let lcarti32 = lxi32 + lyi32 + lzi32;
303    let dl = lcarti32 - li32;
304    if dl % 2 != 0 {
305        return Complex::<f64>::zero();
306    }
307
308    let num = f64::from(
309        (2 * l + 1)
310            * (l - m.unsigned_abs())
311                .checked_factorial()
312                .unwrap_or_else(|| {
313                    panic!(
314                        "Unable to compute the factorial of {}.",
315                        l - m.unsigned_abs()
316                    )
317                }),
318    );
319    let den = 4.0
320        * std::f64::consts::PI
321        * f64::from(
322            (l + m.unsigned_abs())
323                .checked_factorial()
324                .unwrap_or_else(|| {
325                    panic!(
326                        "Unable to compute the factorial of {}.",
327                        l + m.unsigned_abs()
328                    )
329                }),
330        );
331    let mut prefactor =
332        1.0 / f64::from(
333            2u32.pow(l)
334                * l.checked_factorial()
335                    .unwrap_or_else(|| panic!("Unable to compute the factorial of {l}.")),
336        ) * (num / den).sqrt();
337    if csphase && m > 0 {
338        prefactor *=
339            f64::from((-1i32).pow(u32::try_from(m).expect("Unable to convert `m` to `u32`.")));
340    }
341    let ntilde = norm_sph_gaussian(lcart, 1.0);
342    let n = norm_cart_gaussian(lcartqns, 1.0);
343
344    let si =
345        (0..=((l - m.unsigned_abs()).div_euclid(2))).fold(Complex::<f64>::zero(), |acc_si, i| {
346            // i <= (l - |m|) / 2
347            let ii32 =
348                i32::try_from(i).unwrap_or_else(|_| panic!("Cannot convert `{i}` to `i32`."));
349            let mut ifactor = combu(l, i)
350                .to_f64()
351                .expect("Unable to convert a `BigUint` value to `f64`.")
352                * BigUint::from(2 * l - 2 * i)
353                    .checked_factorial()
354                    .unwrap_or_else(|| {
355                        panic!("Unable to compute the factorial of {}.", 2 * l - 2 * i)
356                    })
357                    .to_f64()
358                    .unwrap_or_else(|| {
359                        panic!(
360                            "Unable to convert the factorial of {} to `f64`.",
361                            2 * l - 2 * i
362                        )
363                    })
364                / BigUint::from(l - m.unsigned_abs() - 2 * i)
365                    .checked_factorial()
366                    .unwrap_or_else(|| {
367                        panic!(
368                            "Unable to compute the factorial of {}.",
369                            l - m.unsigned_abs() - 2 * i
370                        )
371                    })
372                    .to_f64()
373                    .unwrap_or_else(|| {
374                        panic!(
375                            "Unable to convert the factorial of {} to `f64`.",
376                            l - m.unsigned_abs() - 2 * i
377                        )
378                    });
379            if i % 2 == 1 {
380                ifactor *= -1.0;
381            };
382            let sp = (0..=(m.unsigned_abs())).fold(Complex::<f64>::zero(), |acc_sp, p| {
383                let pi32 =
384                    i32::try_from(p).unwrap_or_else(|_| panic!("Cannot convert `{p}` to `i32`."));
385                let pfactor = if m > 0 {
386                    combu(m.unsigned_abs(), p)
387                        .to_f64()
388                        .expect("Unable to convert a `BigUint` value to `f64`.")
389                        * Complex::<f64>::i().powu(m.unsigned_abs() - p)
390                } else {
391                    combu(m.unsigned_abs(), p)
392                        .to_f64()
393                        .expect("Unable to convert a `BigUint` value to `f64`.")
394                        * (-1.0 * Complex::<f64>::i()).powu(m.unsigned_abs() - p)
395                };
396                let sq = (0..=(dl.div_euclid(2))).fold(Complex::<f64>::zero(), |acc_sq, q| {
397                    let jq_num = lxi32 + lyi32 - 2 * q - m.abs();
398                    if jq_num.rem_euclid(2) == 0 {
399                        let jq = jq_num.div_euclid(2);
400                        let qfactor = (comb(dl.div_euclid(2), q) * comb(ii32, jq))
401                            .to_f64()
402                            .expect("Unable to convert a `BigUint` value to `f64`.");
403                        let sk = (0..=jq).fold(Complex::<f64>::zero(), |acc_sk, k| {
404                            let tpk_num = lxi32 - pi32 - 2 * k;
405                            if tpk_num.rem_euclid(2) == 0 {
406                                let tpk = tpk_num.div_euclid(2);
407                                let kfactor = (comb(q, tpk) * comb(jq, k))
408                                    .to_f64()
409                                    .expect("Unable to convert a `BigUint` value to `f64`.");
410                                acc_sk + kfactor
411                            } else {
412                                acc_sk
413                            }
414                        });
415                        acc_sq + qfactor * sk
416                    } else {
417                        acc_sq
418                    }
419                });
420                acc_sp + pfactor * sq
421            });
422            acc_si + ifactor * sp
423        });
424    (ntilde / n) * prefactor * si
425}
426
427/// Calculates the overlap between two normalised Cartesian Gaussians of the same order and radial
428/// width, as given in Equation 19 of Schlegel, H. B. & Frisch, M. J. Transformation between
429/// Cartesian and pure spherical harmonic Gaussians. *International Journal of Quantum Chemistry*
430/// **54**, 83–87 (1995), [DOI](https://doi.org/10.1002/qua.560540202).
431///
432/// # Arguments
433///
434/// * `lcartqns1` - A tuple of $`(l_x, l_y, l_z)`$ specifying the exponents of the Cartesian
435/// components of the first Cartesian Gaussian.
436/// * `lcartqns2` - A tuple of $`(l_x, l_y, l_z)`$ specifying the exponents of the Cartesian
437/// components of the first Cartesian Gaussian.
438///
439/// # Returns
440///
441/// The overlap between the two specified normalised Cartesian Gaussians.
442fn cartov(lcartqns1: (u32, u32, u32), lcartqns2: (u32, u32, u32)) -> f64 {
443    let (lx1, ly1, lz1) = lcartqns1;
444    let (lx2, ly2, lz2) = lcartqns2;
445    let lcart1 = lx1 + ly1 + lz1;
446    let lcart2 = lx2 + ly2 + lz2;
447    assert_eq!(
448        lcart1, lcart2,
449        "Only Cartesian Gaussians of the same order are supported."
450    );
451
452    if (lx1 + lx2).rem_euclid(2) == 0
453        && (ly1 + ly2).rem_euclid(2) == 0
454        && (lz1 + lz2).rem_euclid(2) == 0
455    {
456        let num1 = (BigUint::from(lx1 + lx2)
457            .checked_factorial()
458            .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", lx1 + lx2))
459            * BigUint::from(ly1 + ly2)
460                .checked_factorial()
461                .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", ly1 + ly2))
462            * BigUint::from(lz1 + lz2)
463                .checked_factorial()
464                .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", lz1 + lz2)))
465        .to_f64()
466        .expect("Unable to convert a `BigUint` value to `f64`.");
467
468        let den1 = (BigUint::from((lx1 + lx2).div_euclid(2))
469            .checked_factorial()
470            .unwrap_or_else(|| {
471                panic!(
472                    "Unable to compute the factorial of {}.",
473                    (lx1 + lx2).div_euclid(2)
474                )
475            })
476            * BigUint::from((ly1 + ly2).div_euclid(2))
477                .checked_factorial()
478                .unwrap_or_else(|| {
479                    panic!(
480                        "Unable to compute the factorial of {}.",
481                        (ly1 + ly2).div_euclid(2)
482                    )
483                })
484            * BigUint::from((lz1 + lz2).div_euclid(2))
485                .checked_factorial()
486                .unwrap_or_else(|| {
487                    panic!(
488                        "Unable to compute the factorial of {}.",
489                        (lz1 + lz2).div_euclid(2)
490                    )
491                }))
492        .to_f64()
493        .expect("Unable to convert a `BigUint` value to `f64`.");
494
495        let num2 = (BigUint::from(lx1)
496            .checked_factorial()
497            .unwrap_or_else(|| panic!("Unable to compute the factorial of {lx1}."))
498            * BigUint::from(ly1)
499                .checked_factorial()
500                .unwrap_or_else(|| panic!("Unable to compute the factorial of {ly1}."))
501            * BigUint::from(lz1)
502                .checked_factorial()
503                .unwrap_or_else(|| panic!("Unable to compute the factorial of {lz1}."))
504            * BigUint::from(lx2)
505                .checked_factorial()
506                .unwrap_or_else(|| panic!("Unable to compute the factorial of {lx2}."))
507            * BigUint::from(ly2)
508                .checked_factorial()
509                .unwrap_or_else(|| panic!("Unable to compute the factorial of {ly2}."))
510            * BigUint::from(lz2)
511                .checked_factorial()
512                .unwrap_or_else(|| panic!("Unable to compute the factorial of {lz2}.")))
513        .to_f64()
514        .expect("Unable to convert a `BigUint` value to `f64`.");
515
516        let den2 = (BigUint::from(2 * lx1)
517            .checked_factorial()
518            .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", 2 * lx1))
519            * BigUint::from(2 * ly1)
520                .checked_factorial()
521                .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", 2 * ly1))
522            * BigUint::from(2 * lz1)
523                .checked_factorial()
524                .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", 2 * lz1))
525            * BigUint::from(2 * lx2)
526                .checked_factorial()
527                .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", 2 * lx2))
528            * BigUint::from(2 * ly2)
529                .checked_factorial()
530                .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", 2 * ly2))
531            * BigUint::from(2 * lz2)
532                .checked_factorial()
533                .unwrap_or_else(|| panic!("Unable to compute the factorial of {}.", 2 * lz2)))
534        .to_f64()
535        .expect("Unable to convert a `BigUint` value to `f64`.");
536
537        (num1 / den1) * (num2 / den2).sqrt()
538    } else {
539        0.0
540    }
541}
542
543/// Computes the inverse complex coefficients $`c^{-1}(l_x, l_y, l_z, l, m_l, l_{\mathrm{cart}})`$
544/// based on Equation 18 of Schlegel, H. B. & Frisch, M. J. Transformation between
545/// Cartesian and pure spherical harmonic Gaussians. *International Journal of Quantum Chemistry*
546/// **54**, 83–87 (1995), [DOI](https://doi.org/10.1002/qua.560540202), but more generalised for
547/// $`l \leq l_{\mathrm{cart}} = l_x + l_y + l_z`$.
548///
549/// Let $`\tilde{g}(\alpha, l, m_l, l_{\mathrm{cart}}, \mathbf{r})`$ be a complex solid
550/// harmonic Gaussian as defined in Equation 1 of the above reference with
551/// $`n = l_{\mathrm{cart}}`$, and let $`g(\alpha, l_x, l_y, l_z, \mathbf{r})`$ be a Cartesian
552/// Gaussian as defined in Equation 2 of the above reference. The inverse complex coefficients
553/// $`c^{-1}(l_x, l_y, l_z, l, m_l, l_{\mathrm{cart}})`$ effect the inverse transformation
554///
555/// ```math
556/// g(\alpha, l_x, l_y, l_z, \mathbf{r})
557/// = \sum_{l \le l_{\mathrm{cart}} = l_x+l_y+l_z} \sum_{m_l = -l}^{l}
558///     c^{-1}(l_x, l_y, l_z, l, m_l, l_{\mathrm{cart}})
559///     \tilde{g}(\alpha, l, m_l, l_{\mathrm{cart}}, \mathbf{r}).
560/// ```
561///
562/// # Arguments
563///
564/// * `lcartqns` - A tuple of $`(l_x, l_y, l_z)`$ specifying the exponents of the Cartesian
565/// components of the Cartesian Gaussian.
566/// * `lpureqns` - A tuple of $`(l, m_l)`$ specifying the quantum numbers for the spherical
567/// harmonic component of the solid harmonic Gaussian.
568/// * `csphase` - If `true`, the Condon--Shortley phase will be used as defined in
569/// [`complexc`]. If `false`, this phase will be set to unity.
570///
571/// # Returns
572///
573/// $`c^{-1}(l_x, l_y, l_z, l, m_l, l_{\mathrm{cart}})`$.
574pub fn complexcinv(lcartqns: (u32, u32, u32), lpureqns: (u32, i32), csphase: bool) -> Complex<f64> {
575    let (lx, ly, lz) = lcartqns;
576    let lcart = lx + ly + lz;
577    let mut cinv = Complex::<f64>::zero();
578    for lx2 in 0..=lcart {
579        for ly2 in 0..=(lcart - lx2) {
580            let lz2 = lcart - lx2 - ly2;
581            cinv += cartov(lcartqns, (lx2, ly2, lz2))
582                * complexc(lpureqns, (lx2, ly2, lz2), csphase).conj();
583        }
584    }
585    cinv
586}
587
588/// Obtains the transformation matrix $`\boldsymbol{\Upsilon}^{(l)}`$ allowing complex spherical
589/// harmonics to be expressed as linear combinations of real spherical harmonics.
590///
591/// Let $`Y_{lm}`$ be a real spherical harmonic of degree $`l`$. Then, a complex spherical
592/// harmonic of degree $`l`$ and order $`m`$ is given by
593///
594/// ```math
595/// Y_l^m =
596///     \begin{cases}
597///         \frac{\lambda_{\mathrm{cs}}}{\sqrt{2}}
598///         \left(Y_{l\lvert m \rvert}
599///               - \mathbb{i} Y_{l,-\lvert m \rvert}\right)
600///         & \mathrm{if}\ m < 0 \\
601///         Y_{l0} & \mathrm{if}\ m = 0 \\
602///         \frac{\lambda_{\mathrm{cs}}}{\sqrt{2}}
603///         \left(Y_{l\lvert m \rvert}
604///               + \mathbb{i} Y_{l,-\lvert m \rvert}\right)
605///         & \mathrm{if}\ m > 0 \\
606///     \end{cases}
607/// ```
608///
609/// where $`\lambda_{\mathrm{cs}}`$ is the Condon--Shortley phase as defined in [`complexc`].
610/// The linear combination coefficients can then be gathered into a square matrix
611/// $`\boldsymbol{\Upsilon}^{(l)}`$ of dimensions $`(2l+1)\times(2l+1)`$ such that
612///
613/// ```math
614///     Y_l^m = \sum_{m'} Y_{lm'} \Upsilon^{(l)}_{m'm}.
615/// ```
616///
617/// # Arguments
618///
619/// * `l` - The spherical harmonic degree.
620/// * `csphase` - If `true`, $`\lambda_{\mathrm{cs}}`$ is as defined in [`complexc`]. If `false`,
621/// $`\lambda_{\mathrm{cs}} = 1`$.
622/// * `pureorder` - A [`PureOrder`] struct giving the ordering of the components of the pure
623/// Gaussians.
624///
625/// # Returns
626///
627/// The $`\boldsymbol{\Upsilon}^{(l)}`$ matrix.
628pub fn sh_c2r_mat(l: u32, csphase: bool, pureorder: &PureOrder) -> Array2<Complex<f64>> {
629    assert_eq!(pureorder.lpure, l, "Mismatched pure ranks.");
630    let lusize = l as usize;
631    let mut upmat = Array2::<Complex<f64>>::zeros((2 * lusize + 1, 2 * lusize + 1));
632    let po_il = PureOrder::increasingm(l);
633    for &mcomplex in po_il.iter() {
634        let absmreal = mcomplex.unsigned_abs() as usize;
635        match mcomplex.cmp(&0) {
636            Ordering::Less => {
637                // Python-equivalent:
638                // upmat[-absmreal + l, mcomplex + l] = -1.0j / np.sqrt(2)
639                // upmat[+absmreal + l, mcomplex + l] = 1.0 / np.sqrt(2)
640                // mcomplex = -absmreal
641                upmat[(lusize - absmreal, lusize - absmreal)] =
642                    Complex::<f64>::new(0.0, -1.0 / 2.0f64.sqrt());
643                upmat[(lusize + absmreal, lusize - absmreal)] =
644                    Complex::<f64>::new(1.0 / 2.0f64.sqrt(), 0.0);
645            }
646            Ordering::Equal => {
647                upmat[(lusize, lusize)] = Complex::<f64>::from(1.0);
648            }
649            Ordering::Greater => {
650                let lcs = if csphase {
651                    f64::from((-1i32).pow(
652                        u32::try_from(mcomplex).expect("Unable to convert `mcomplex` to `u32`."),
653                    ))
654                } else {
655                    1.0
656                };
657                // Python-equivalent:
658                // upmat[-absmreal + l, mcomplex + l] = lcs * 1.0j / np.sqrt(2)
659                // upmat[+absmreal + l, mcomplex + l] = lcs * 1.0 / np.sqrt(2)
660                // mcomplex = absmreal
661                upmat[(lusize - absmreal, lusize + absmreal)] =
662                    lcs * Complex::<f64>::new(0.0, 1.0 / 2.0f64.sqrt());
663                upmat[(lusize + absmreal, lusize + absmreal)] =
664                    lcs * Complex::<f64>::new(1.0 / 2.0f64.sqrt(), 0.0);
665            }
666        }
667    }
668
669    // upmat is always in increasing-m order. We now permute, if required.
670    if *pureorder != po_il {
671        let perm = pureorder.get_perm_of(&po_il).expect(
672            "Permutation to obtain `pureorder` from the increasing-m order could not be found.",
673        );
674        let image = perm.image();
675        upmat.select(Axis(0), &image).select(Axis(1), &image)
676    } else {
677        upmat
678    }
679}
680
681/// Obtains the matrix $`\boldsymbol{\Upsilon}^{(l)\dagger}`$ allowing real spherical harmonics
682/// to be expressed as linear combinations of complex spherical harmonics.
683///
684/// Let $`Y_l^m`$ be a complex spherical harmonic of degree $`l`$ and order $`m`$.
685/// Then, a real degree-$`l`$ spherical harmonic $`Y_{lm}`$ can be defined as
686///
687/// ```math
688/// Y_{lm} =
689///     \begin{cases}
690///         \frac{\mathbb{i}}{\sqrt{2}}
691///         \left(Y_l^{-\lvert m \rvert}
692///               - \lambda'_{\mathrm{cs}} Y_l^{\lvert m \rvert}\right)
693///         & \mathrm{if}\ m < 0 \\
694///         Y_l^0 & \mathrm{if}\ m = 0 \\
695///         \frac{1}{\sqrt{2}}
696///         \left(Y_l^{-\lvert m \rvert}
697///               + \lambda'_{\mathrm{cs}} Y_l^{\lvert m \rvert}\right)
698///         & \mathrm{if}\ m > 0 \\
699///     \end{cases}
700/// ```
701///
702/// where $`\lambda'_{\mathrm{cs}} = (-1)^{\lvert m \rvert}`$ if the Condon--Shortley phase as
703/// defined in [`complexc`] is employed for the complex spherical harmonics, and
704/// $`\lambda'_{\mathrm{cs}} = 1`$ otherwise. The linear combination coefficients turn out to be
705/// given by the elements of matrix $`\boldsymbol{\Upsilon}^{(l)\dagger}`$ of dimensions
706/// $`(2l+1)\times(2l+1)`$ such that
707///
708/// ```math
709///     Y_{lm} = \sum_{m'} Y_l^{m'} [\Upsilon^{(l)\dagger}]_{m'm}.
710/// ```
711///
712/// It is obvious from the orthonormality of $`Y_{lm}`$ and $`Y_l^m`$ that
713/// $`\boldsymbol{\Upsilon}^{(l)\dagger} = [\boldsymbol{\Upsilon}^{(l)}]^{-1}`$ where
714/// $`\boldsymbol{\Upsilon}^{(l)}`$ is defined in [`sh_c2r_mat`].
715///
716/// # Arguments
717///
718/// * `l` - The spherical harmonic degree.
719/// * `csphase` - If `true`, $`\lambda_{\mathrm{cs}}`$ is as defined in [`complexc`]. If `false`,
720/// $`\lambda_{\mathrm{cs}} = 1`$.
721/// * `pureorder` - A [`PureOrder`] struct giving the ordering of the components of the pure
722/// Gaussians.
723///
724/// # Returns
725///
726/// The $`\boldsymbol{\Upsilon}^{(l)\dagger}`$ matrix.
727pub fn sh_r2c_mat(l: u32, csphase: bool, pureorder: &PureOrder) -> Array2<Complex<f64>> {
728    let mut mat = sh_c2r_mat(l, csphase, pureorder).t().to_owned();
729    mat.par_mapv_inplace(|x| x.conj());
730    mat
731}
732
733/// Obtains the matrix $`\mathbf{U}^{(l_{\mathrm{cart}}, l)}`$ containing linear combination
734/// coefficients of Cartesian Gaussians in the expansion of a complex solid harmonic Gaussian,
735/// *i.e.*, briefly,
736///
737/// ```math
738/// \tilde{\mathbf{g}}^{\mathsf{T}}(l)
739///     = \mathbf{g}^{\mathsf{T}}(l_{\mathrm{cart}})
740///     \ \mathbf{U}^{(l_{\mathrm{cart}}, l)}.
741/// ```
742///
743/// Let $`\tilde{g}(\alpha, \lambda, l_{\mathrm{cart}}, \mathbf{r})`$ be a complex solid harmonic
744/// Gaussian as defined in Equation 1 of Schlegel, H. B. & Frisch, M. J. Transformation between
745/// Cartesian and pure spherical harmonic Gaussians. *International Journal of Quantum Chemistry*
746/// **54**, 83–87 (1995), [DOI](https://doi.org/10.1002/qua.560540202) with
747/// $`n = l_{\mathrm{cart}}`$, and let $`g(\alpha, \lambda_{\mathrm{cart}}, \mathbf{r})`$ be a
748/// Cartesian Gaussian as defined in Equation 2 of the above reference.
749/// Here, $`\lambda`$ is a single index labelling a complex solid harmonic Gaussian of spherical
750/// harmonic degree $`l`$ and order $`m_l`$, and $`\lambda_{\mathrm{cart}}`$ a single index
751/// labelling a Cartesian Gaussian of degrees $`(l_x, l_y, l_z)`$ such that
752/// $`l_x + l_y + l_z = l_{\mathrm{cart}}`$. We can then write
753///
754/// ```math
755/// \tilde{g}(\alpha, \lambda, l_{\mathrm{cart}}, \mathbf{r})
756/// = \sum_{\lambda_{\mathrm{cart}}}
757///     g(\alpha, \lambda_{\mathrm{cart}}, \mathbf{r})
758///     U^{(l_{\mathrm{cart}}, l)}_{\lambda_{\mathrm{cart}}\lambda}
759/// ```
760///
761/// where $`U^{(l_{\mathrm{cart}}, l)}_{\lambda_{\mathrm{cart}}\lambda}`$
762/// is given by the complex coefficients
763///
764/// ```math
765/// U^{(l_{\mathrm{cart}}, l)}_{\lambda_{\mathrm{cart}}\lambda} =
766///     c(l, m_l, l_{\mathrm{cart}}, l_x, l_y, l_z)
767/// ```
768///
769/// defined in [`complexc`].
770///
771/// $`\mathbf{U}^{(l_{\mathrm{cart}}, l)}`$ has dimensions
772/// $`\frac{1}{2}(l_{\mathrm{cart}}+1)(l_{\mathrm{cart}}+2) \times (2l+1)`$ and contains only
773/// zero elements if $`l`$ and $`l_{\mathrm{cart}}`$ have different parities.
774/// It can be verified that
775/// $`\mathbf{V}^{(l,l_{\mathrm{cart}})}
776/// \ \mathbf{U}^{(l_{\mathrm{cart}}, l)} = \boldsymbol{I}_{2l+1}`$, where
777/// $`\mathbf{V}^{(l,l_{\mathrm{cart}})}`$ is given in [`sh_cart2cl_mat`].
778///
779/// # Arguments
780///
781/// * `lcart` - The total Cartesian degree for the Cartesian Gaussians and
782///  also for the radial part of the solid harmonic Gaussian.
783/// * `l` - The degree of the complex spherical harmonic factor in the solid
784///  harmonic Gaussian.
785/// * `cartorder` - A [`CartOrder`] struct giving the ordering of the components of the Cartesian
786/// Gaussians.
787/// * `csphase` - Set to `true` to use the Condon--Shortley phase in the calculations of the $`c`$
788/// coefficients. See [`complexc`] for more details.
789/// * `pureorder` - A [`PureOrder`] struct giving the ordering of the components of the pure
790/// Gaussians.
791///
792/// # Returns
793///
794/// The $`\mathbf{U}^{(l_{\mathrm{cart}}, l)}`$ matrix.
795pub fn sh_cl2cart_mat(
796    lcart: u32,
797    l: u32,
798    cartorder: &CartOrder,
799    csphase: bool,
800    pureorder: &PureOrder,
801) -> Array2<Complex<f64>> {
802    assert_eq!(cartorder.lcart, lcart, "Mismatched Cartesian ranks.");
803    assert_eq!(pureorder.lpure, l, "Mismatched pure ranks.");
804    let mut umat = Array2::<Complex<f64>>::zeros((
805        ((lcart + 1) * (lcart + 2)).div_euclid(2) as usize,
806        2 * l as usize + 1,
807    ));
808    for (i, &m) in pureorder.iter().enumerate() {
809        for (icart, &lcartqns) in cartorder.iter().enumerate() {
810            umat[(icart, i)] = complexc((l, m), lcartqns, csphase);
811        }
812    }
813    umat
814}
815
816/// Obtains the matrix $`\mathbf{V}^{(l, l_{\mathrm{cart}})}`$ containing linear combination
817/// coefficients of complex solid harmonic Gaussians of a specific degree in the expansion of
818/// Cartesian Gaussians, *i.e.*, briefly,
819///
820/// ```math
821/// \mathbf{g}^{\mathsf{T}}(l_{\mathrm{cart}})
822///     = \tilde{\mathbf{g}}^{\mathsf{T}}(l)
823///     \ \mathbf{V}^{(l, l_{\mathrm{cart}})}.
824/// ```
825///
826/// Let $`\tilde{g}(\alpha, \lambda, l_{\mathrm{cart}}, \mathbf{r})`$ be a complex solid harmonic
827/// Gaussian as defined in Equation 1 of Schlegel, H. B. & Frisch, M. J. Transformation between
828/// Cartesian and pure spherical harmonic Gaussians. *International Journal of Quantum Chemistry*
829/// **54**, 83–87 (1995), [DOI](https://doi.org/10.1002/qua.560540202) with
830/// $`n = l_{\mathrm{cart}}`$, and let $`g(\alpha, \lambda_{\mathrm{cart}}, \mathbf{r})`$ be a
831/// Cartesian Gaussian as defined in Equation 2 of the above reference.  Here, $`\lambda`$ is a
832/// single index labelling a complex solid harmonic Gaussian of spherical harmonic degree $`l`$
833/// and order $`m_l`$, and $`\lambda_{\mathrm{cart}}`$ a single index labelling a Cartesian
834/// Gaussian of degrees $`(l_x, l_y, l_z)`$ such that $`l_x + l_y + l_z = l_{\mathrm{cart}}`$.
835/// We can then write
836///
837/// ```math
838/// g(\alpha, \lambda_{\mathrm{cart}}, \mathbf{r})
839/// = \sum_{\substack{\lambda\\ l \leq l_{\mathrm{cart}}}}
840///     \tilde{g}(\alpha, \lambda, l_{\mathrm{cart}}, \mathbf{r})
841///     V^{(l_{\mathrm{cart}})}_{\lambda\lambda_{\mathrm{cart}}}
842/// ```
843///
844/// where $`V^{(l_{\mathrm{cart}})}_{\lambda\lambda_{\mathrm{cart}}}`$ is given by the inverse
845/// complex coefficients
846///
847/// ```math
848/// V^{(l_{\mathrm{cart}})}_{\lambda\lambda_{\mathrm{cart}}} =
849///     c^{-1}(l_x, l_y, l_z, l, m_l, l_{\mathrm{cart}})
850/// ```
851///
852/// defined in [`complexcinv`].
853///
854/// We can order the rows $`\lambda`$ of $`\mathbf{V}^{(l_{\mathrm{cart}})}`$ that have the same
855/// $`l`$ into rectangular blocks of dimensions
856/// $`(2l+1) \times \frac{1}{2}(l_{\mathrm{cart}}+1)(l_{\mathrm{cart}}+2)`$
857/// which give contributions from complex solid harmonic Gaussians of a particular degree $`l`$.
858/// We denote these blocks $`\mathbf{V}^{(l, l_{\mathrm{cart}})}`$.
859/// They contain only zero elements if $`l`$ and $`l_{\mathrm{cart}}`$ have different parities.
860///
861/// # Arguments
862///
863/// * `l` - The degree of the complex spherical harmonic factor in the solid
864///  harmonic Gaussian.
865/// * `lcart` - The total Cartesian degree for the Cartesian Gaussians and
866///  also for the radial part of the solid harmonic Gaussian.
867/// * `cartorder` - A [`CartOrder`] struct giving the ordering of the components of the Cartesian
868/// Gaussians.
869/// * `csphase` - Set to `true` to use the Condon--Shortley phase in the calculations of the
870/// $`c^{-1}`$ coefficients. See [`complexc`] and [`complexcinv`] for more details.
871/// * `pureorder` - A [`PureOrder`] struct giving the ordering of the components of the pure
872/// Gaussians.
873///
874/// # Returns
875///
876/// The $`\mathbf{V}^{(l, l_{\mathrm{cart}})}`$ block.
877pub fn sh_cart2cl_mat(
878    l: u32,
879    lcart: u32,
880    cartorder: &CartOrder,
881    csphase: bool,
882    pureorder: &PureOrder,
883) -> Array2<Complex<f64>> {
884    assert_eq!(pureorder.lpure, l, "Mismatched pure ranks.");
885    assert_eq!(cartorder.lcart, lcart, "Mismatched Cartesian ranks.");
886    let mut vmat = Array2::<Complex<f64>>::zeros((
887        2 * l as usize + 1,
888        ((lcart + 1) * (lcart + 2)).div_euclid(2) as usize,
889    ));
890    for (icart, &lcartqns) in cartorder.iter().enumerate() {
891        for (i, &m) in pureorder.iter().enumerate() {
892            vmat[(i, icart)] = complexcinv(lcartqns, (l, m), csphase);
893        }
894    }
895    vmat
896}
897
898/// Obtain the matrix $`\mathbf{W}^{(l_{\mathrm{cart}}, l)}`$ containing linear combination
899/// coefficients of Cartesian Gaussians in the expansion of a real solid harmonic Gaussian, *i.e.*,
900/// briefly,
901///
902/// ```math
903/// \bar{\mathbf{g}}^{\mathsf{T}}(l)
904///     = \mathbf{g}^{\mathsf{T}}(l_{\mathrm{cart}})
905///     \ \mathbf{W}^{(l_{\mathrm{cart}}, l)}.
906/// ```
907///
908/// Let $`\bar{g}(\alpha, \lambda, l_{\mathrm{cart}}, \mathbf{r})`$ be
909/// a real solid harmonic Gaussian defined in a similar manner to Equation 1 of Schlegel, H. B.
910/// & Frisch, M. J. Transformation between Cartesian and pure spherical harmonic Gaussians.
911/// *International Journal of Quantum Chemistry* **54**, 83–87 (1995),
912/// [DOI](https://doi.org/10.1002/qua.560540202) with $`n = l_{\mathrm{cart}}`$ but with real
913/// rather than complex spherical harmonic factors, and let
914/// $`g(\alpha, \lambda_{\mathrm{cart}}, \mathbf{r})`$ be a Cartesian Gaussian as defined in
915/// Equation 2 of the above reference. Here, $`\lambda`$ is a single index labelling a complex
916/// solid harmonic Gaussian of spherical harmonic degree $`l`$ and order $`m_l`$, and
917/// $`\lambda_{\mathrm{cart}}`$ a single index labelling a Cartesian Gaussian of degrees
918/// $`(l_x, l_y, l_z)`$ such that $`l_x + l_y + l_z = l_{\mathrm{cart}}`$. We can then write
919///
920/// ```math
921/// \bar{g}(\alpha, \lambda, l_{\mathrm{cart}}, \mathbf{r})
922/// = \sum_{\lambda_{\mathrm{cart}}}
923///     g(\alpha, \lambda_{\mathrm{cart}}, \mathbf{r})
924///     W^{(l_{\mathrm{cart}}, l)}_{\lambda_{\mathrm{cart}}\lambda}.
925/// ```
926///
927/// $`\mathbf{W}^{(l_{\mathrm{cart}}, l)}`$ is given by
928///
929/// ```math
930/// \mathbf{W}^{(l_{\mathrm{cart}}, l)}
931/// = \mathbf{U}^{(l_{\mathrm{cart}}, l)}
932///   \boldsymbol{\Upsilon}^{(l)\dagger},
933/// ```
934///
935/// where $`\boldsymbol{\Upsilon}^{(l)\dagger}`$ is defined in [`sh_r2c_mat`] and
936/// $`\mathbf{U}^{(l_{\mathrm{cart}}, l)}`$ in [`sh_cl2cart_mat`].
937/// $`\mathbf{W}^{(l_{\mathrm{cart}}, l)}`$ must be real.
938/// $`\mathbf{W}^{(l_{\mathrm{cart}}, l)}`$ has dimensions
939/// $`\frac{1}{2}(l_{\mathrm{cart}}+1)(l_{\mathrm{cart}}+2) \times (2l+1)`$ and contains only zero
940/// elements if $`l`$ and $`l_{\mathrm{cart}}`$ have different parities. It can be verified that
941/// $`\mathbf{X}^{(l,l_{\mathrm{cart}})}
942/// \ \mathbf{W}^{(l_{\mathrm{cart}}, l)} = \boldsymbol{I}_{2l+1}`$, where
943/// $`\mathbf{X}^{(l,l_{\mathrm{cart}})}`$ is given in
944/// [`sh_cart2rl_mat`].
945///
946/// # Arguments
947///
948/// * lcart - The total Cartesian degree for the Cartesian Gaussians and
949///  also for the radial part of the solid harmonic Gaussian.
950/// * l - The degree of the complex spherical harmonic factor in the solid
951///  harmonic Gaussian.
952/// * cartorder - A [`CartOrder`] struct giving the ordering of the components of the Cartesian
953/// Gaussians.
954/// * `csphase` - Set to `true` to use the Condon--Shortley phase in the calculations of the $`c`$
955/// coefficients. See [`complexc`] for more details.
956/// * `pureorder` - A [`PureOrder`] struct giving the ordering of the components of the pure
957/// Gaussians.
958///
959/// # Returns
960///
961/// The $`\mathbf{W}^{(l_{\mathrm{cart}}, l)}`$ matrix.
962pub fn sh_rl2cart_mat(
963    lcart: u32,
964    l: u32,
965    cartorder: &CartOrder,
966    csphase: bool,
967    pureorder: &PureOrder,
968) -> Array2<f64> {
969    assert_eq!(cartorder.lcart, lcart, "Mismatched Cartesian ranks.");
970    assert_eq!(pureorder.lpure, l, "Mismatched pure ranks.");
971    let upmatdagger = sh_r2c_mat(l, csphase, pureorder);
972    let umat = sh_cl2cart_mat(lcart, l, cartorder, csphase, pureorder);
973    let wmat = umat.dot(&upmatdagger);
974    assert!(
975        wmat.iter()
976            .all(|x| approx::relative_eq!(x.im, 0.0, max_relative = 1e-7, epsilon = 1e-7)),
977        "wmat is not entirely real."
978    );
979    wmat.map(|x| x.re)
980}
981
982/// Obtains the real matrix $`\mathbf{X}^{(l, l_{\mathrm{cart}})}`$ containing linear combination
983/// coefficients of real solid harmonic Gaussians of a specific degree in the expansion of
984/// Cartesian Gaussians, *i.e.*, briefly,
985///
986/// ```math
987/// \mathbf{g}^{\mathsf{T}}(l_{\mathrm{cart}})
988///     = \bar{\mathbf{g}}^{\mathsf{T}}(l)
989///     \ \mathbf{X}^{(l, l_{\mathrm{cart}})}.
990/// ```
991///
992/// Let $`\bar{g}(\alpha, \lambda, l_{\mathrm{cart}}, \mathbf{r})`$ be a real solid harmonic
993/// Gaussian defined in a similar manner to Equation 1 of Schlegel, H. B. & Frisch, M. J.
994/// Transformation between Cartesian and pure spherical harmonic Gaussians. *International
995/// Journal of Quantum Chemistry* **54**, 83–87 (1995),
996/// [DOI](https://doi.org/10.1002/qua.560540202)
997/// with $`n = l_{\mathrm{cart}}`$, but with real rather than complex spherical harmonic factors,
998/// and let $`g(\alpha, \lambda_{\mathrm{cart}}, \mathbf{r})`$ be a Cartesian Gaussian as defined
999/// in Equation 2 of the above reference.  Here, $`\lambda`$ is a single index labelling a real
1000/// solid harmonic Gaussian of spherical harmonic degree $`l`$ and real order $`m_l`$, and
1001/// $`\lambda_{\mathrm{cart}}`$ a single index labelling a Cartesian Gaussian of degrees
1002/// $`(l_x, l_y, l_z)`$ such that $`l_x + l_y + l_z = l_{\mathrm{cart}}`$.
1003/// We can then write
1004///
1005/// ```math
1006/// g(\alpha, \lambda_{\mathrm{cart}}, \mathbf{r})
1007/// = \sum_{\substack{\lambda\\ l \leq l_{\mathrm{cart}}}}
1008///     \bar{g}(\alpha, \lambda, l_{\mathrm{cart}}, \mathbf{r})
1009///     X^{(l_{\mathrm{cart}})}_{\lambda\lambda_{\mathrm{cart}}}.
1010/// ```
1011///
1012/// We can order the rows $`\lambda`$ of $`\mathbf{X}^{(l_{\mathrm{cart}})}`$ that have the same
1013/// $`l`$ into rectangular blocks of dimensions
1014/// $`(2l+1) \times \frac{1}{2}(l_{\mathrm{cart}}+1)(l_{\mathrm{cart}}+2)`$.
1015/// We denote these blocks $`\mathbf{X}^{(l, l_{\mathrm{cart}})}`$ which are given by
1016///
1017/// ```math
1018/// \mathbf{X}^{(l, l_{\mathrm{cart}})}
1019/// = \boldsymbol{\Upsilon}^{(l)} \mathbf{V}^{(l, l_{\mathrm{cart}})},
1020/// ```
1021///
1022/// where $`\boldsymbol{\Upsilon}^{(l)}`$ is defined in
1023/// [`sh_c2r_mat`] and $`\boldsymbol{V}^{(l, l_{\mathrm{cart}})}`$ in [`sh_cart2cl_mat`].
1024/// $`\mathbf{X}^{(l, l_{\mathrm{cart}})}`$ must be real.
1025///
1026/// # Arguments
1027///
1028/// * `l` - The degree of the complex spherical harmonic factor in the solid
1029///  harmonic Gaussian.
1030/// * `lcart` - The total Cartesian degree for the Cartesian Gaussians and
1031///  also for the radial part of the solid harmonic Gaussian.
1032/// * `cartorder` - A [`CartOrder`] struct giving the ordering of the components of the Cartesian
1033/// Gaussians.
1034/// * `csphase` - Set to `true` to use the Condon--Shortley phase in the calculations of the
1035/// $`c^{-1}`$ coefficients. See [`complexc`] and [`complexcinv`] for more details.
1036/// * `pureorder` - A [`PureOrder`] struct giving the ordering of the components of the pure
1037/// Gaussians.
1038///
1039/// # Returns
1040///
1041/// The $`\mathbf{X}^{(l, l_{\mathrm{cart}})}`$ block.
1042pub fn sh_cart2rl_mat(
1043    l: u32,
1044    lcart: u32,
1045    cartorder: &CartOrder,
1046    csphase: bool,
1047    pureorder: &PureOrder,
1048) -> Array2<f64> {
1049    assert_eq!(cartorder.lcart, lcart, "Mismatched Cartesian ranks.");
1050    assert_eq!(pureorder.lpure, l, "Mismatched pure ranks.");
1051    let upmat = sh_c2r_mat(l, csphase, pureorder);
1052    let vmat = sh_cart2cl_mat(l, lcart, cartorder, csphase, pureorder);
1053    let xmat = upmat.dot(&vmat);
1054    assert!(
1055        xmat.iter()
1056            .all(|x| approx::relative_eq!(x.im, 0.0, max_relative = 1e-7, epsilon = 1e-7)),
1057        "xmat is not entirely real."
1058    );
1059    xmat.map(|x| x.re)
1060}
1061
1062/// Returns a list of $`\mathbf{W}^{(l_{\mathrm{cart}}, l)}`$ for
1063/// $`l_{\mathrm{cart}} \ge l \ge 0`$ and $`l \equiv l_{\mathrm{cart}} \mod 2`$.
1064///
1065/// $`\mathbf{W}^{(l_{\mathrm{cart}}, l)}`$ is defined in [`sh_rl2cart_mat`].
1066///
1067/// # Arguments
1068///
1069/// * `lcart` - The total Cartesian degree for the Cartesian Gaussians and
1070///  also for the radial part of the solid harmonic Gaussian.
1071/// * `cartorder` - A [`CartOrder`] struct giving the ordering of the components of the Cartesian
1072/// Gaussians.
1073/// * `csphase` - Set to `true` to use the Condon--Shortley phase in the calculations of the
1074/// $`c`$ coefficients. See [`complexc`] for more details.
1075/// * `pureorder` - A closure to generate a [`PureOrder`] struct giving the ordering of the
1076/// components of the pure Gaussians for a particular value of `l`.
1077///
1078/// # Returns
1079///
1080/// A vector of $`\mathbf{W}^{(l_{\mathrm{cart}}, l)}`$ matrices with
1081/// $`l_{\mathrm{cart}} \ge l \ge 0`$ and $`l \equiv l_{\mathrm{cart}} \mod 2`$ in decreasing
1082/// $`l`$ order.
1083pub fn sh_r2cart(
1084    lcart: u32,
1085    cartorder: &CartOrder,
1086    csphase: bool,
1087    pureorder: fn(u32) -> PureOrder,
1088) -> Vec<Array2<f64>> {
1089    assert_eq!(cartorder.lcart, lcart, "Mismatched Cartesian ranks.");
1090    let lrange = if lcart.rem_euclid(2) == 0 {
1091        #[allow(clippy::range_plus_one)]
1092        (0..lcart + 1).step_by(2).rev()
1093    } else {
1094        #[allow(clippy::range_plus_one)]
1095        (1..lcart + 1).step_by(2).rev()
1096    };
1097    lrange
1098        .map(|l| sh_rl2cart_mat(lcart, l, cartorder, csphase, &pureorder(l)))
1099        .collect()
1100}
1101
1102/// Returns a list of $`\mathbf{X}^{(l, l_{\mathrm{cart}})}`$ for
1103/// $`l_{\mathrm{cart}} \ge l \ge 0`$ and $`l \equiv l_{\mathrm{cart}} \mod 2`$.
1104///
1105/// $`\mathbf{X}^{(l, l_{\mathrm{cart}})}`$ is defined in [`sh_cart2rl_mat`].
1106///
1107/// # Arguments
1108///
1109/// * `lcart` - The total Cartesian degree for the Cartesian Gaussians and
1110///  also for the radial part of the solid harmonic Gaussian.
1111/// * `cartorder` - A [`CartOrder`] struct giving the ordering of the components of the Cartesian
1112/// Gaussians.
1113/// * `csphase` - Set to `true` to use the Condon--Shortley phase in the calculations of the
1114/// $`c^{-1}`$ coefficients. See [`complexc`] and [`complexcinv`] for more details.
1115/// * `pureorder` - A closure to generate a [`PureOrder`] struct giving the ordering of the
1116/// components of the pure Gaussians for a particular value of `l`.
1117///
1118/// # Returns
1119///
1120/// A vector of $`\mathbf{X}^{(l, l_{\mathrm{cart}})}`$ matrices with
1121/// $`l_{\mathrm{cart}} \ge l \ge 0`$ and $`l \equiv l_{\mathrm{cart}} \mod 2`$ in decreasing
1122/// $`l`$ order.
1123pub fn sh_cart2r(
1124    lcart: u32,
1125    cartorder: &CartOrder,
1126    csphase: bool,
1127    pureorder: fn(u32) -> PureOrder,
1128) -> Vec<Array2<f64>> {
1129    assert_eq!(cartorder.lcart, lcart, "Mismatched Cartesian ranks.");
1130    let lrange = if lcart.rem_euclid(2) == 0 {
1131        #[allow(clippy::range_plus_one)]
1132        (0..lcart + 1).step_by(2).rev()
1133    } else {
1134        #[allow(clippy::range_plus_one)]
1135        (1..lcart + 1).step_by(2).rev()
1136    };
1137    lrange
1138        .map(|l| sh_cart2rl_mat(l, lcart, cartorder, csphase, &pureorder(l)))
1139        .collect()
1140}