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}