qsym2/basis/
ao.rs

1//! Atomic-orbital basis functions.
2
3use std::cmp::Ordering;
4use std::collections::{HashMap, HashSet};
5use std::convert::TryInto;
6use std::fmt;
7use std::slice::Iter;
8
9use anyhow::{self, ensure, format_err};
10use counter::Counter;
11use derive_builder::Builder;
12use itertools::{izip, Itertools};
13
14use crate::angmom::ANGMOM_LABELS;
15use crate::auxiliary::atom::Atom;
16use crate::auxiliary::misc::ProductRepeat;
17use crate::permutation::{permute_inplace, PermutableCollection, Permutation};
18
19#[cfg(test)]
20#[path = "ao_tests.rs"]
21mod ao_tests;
22
23// -------------------
24// Shell order structs
25// -------------------
26
27// ~~~~~~~~~
28// PureOrder
29// ~~~~~~~~~
30
31/// Structure to contain information about the ordering of pure Gaussians of a certain rank.
32#[derive(Clone, Builder, PartialEq, Eq, Hash)]
33pub struct PureOrder {
34    /// A sequence of $`m_l`$ values giving the ordering of the pure Gaussians.
35    #[builder(setter(custom))]
36    mls: Vec<i32>,
37
38    /// The rank of the pure Gaussians.
39    pub lpure: u32,
40
41    /// The spatial parity of the pure Gaussians: `true` if even under spatial inversion and
42    /// `false` if odd. By default, this is the same as the parity of `Self::lpure`.
43    #[builder(default = "self.default_even()?")]
44    pub even: bool,
45}
46
47impl PureOrderBuilder {
48    fn mls(&mut self, mls: &[i32]) -> &mut Self {
49        let lpure = self.lpure.expect("`lpure` has not been set.");
50        assert!(
51            mls.iter()
52                .map(|m| m.unsigned_abs())
53                .max()
54                .expect("The maximum |m| value could not be determined.")
55                == lpure
56        );
57        assert_eq!(mls.len(), (2 * lpure + 1) as usize);
58        self.mls = Some(mls.to_vec());
59        self
60    }
61
62    fn default_even(&self) -> Result<bool, String> {
63        let lpure = self
64            .lpure
65            .ok_or_else(|| "`lpure` has not been set.".to_string())?;
66        Ok(lpure % 2 == 0)
67    }
68}
69
70impl PureOrder {
71    /// Returns a builder to construct a new [`PureOrder`] structure.
72    fn builder() -> PureOrderBuilder {
73        PureOrderBuilder::default()
74    }
75
76    /// Constructs a new [`PureOrder`] structure from its constituting $`m_l`$ values. The spatial
77    /// inversion parity is determined by the parity of the deduced $`l`$ value.
78    pub fn new(mls: &[i32]) -> Result<Self, anyhow::Error> {
79        let lpure = mls
80            .iter()
81            .map(|m| m.unsigned_abs())
82            .max()
83            .expect("The maximum |m| value could not be determined.");
84        let pure_order = PureOrder::builder()
85            .lpure(lpure)
86            .mls(mls)
87            .build()
88            .map_err(|err| format_err!(err))?;
89        ensure!(pure_order.verify(), "Invalid `PureOrder`.");
90        Ok(pure_order)
91    }
92
93    /// Constructs a new [`PureOrder`] structure for a specified rank with increasing-$`m`$ order.
94    /// The spatial inversion parity is determined by the parity of the deduced $`l`$ value.
95    ///
96    /// # Arguments
97    ///
98    /// * `lpure` - The required pure Gaussian rank.
99    ///
100    /// # Returns
101    ///
102    /// A [`PureOrder`] struct for a specified rank with increasing-$`m`$ order.
103    #[must_use]
104    pub fn increasingm(lpure: u32) -> Self {
105        let lpure_i32 = i32::try_from(lpure).expect("`lpure` cannot be converted to `i32`.");
106        let mls = (-lpure_i32..=lpure_i32).collect_vec();
107        Self::builder()
108            .lpure(lpure)
109            .mls(&mls)
110            .build()
111            .expect("Unable to construct a `PureOrder` structure with increasing-m order.")
112    }
113
114    /// Constructs a new [`PureOrder`] structure for a specified rank with decreasing-$`m`$ order.
115    /// The spatial inversion parity is determined by the parity of the deduced $`l`$ value.
116    ///
117    /// # Arguments
118    ///
119    /// * `lpure` - The required pure Gaussian rank.
120    ///
121    /// # Returns
122    ///
123    /// A [`PureOrder`] struct for a specified rank with decreasing-$`m`$ order.
124    #[must_use]
125    pub fn decreasingm(lpure: u32) -> Self {
126        let lpure_i32 = i32::try_from(lpure).expect("`lpure` cannot be converted to `i32`.");
127        let mls = (-lpure_i32..=lpure_i32).rev().collect_vec();
128        Self::builder()
129            .lpure(lpure)
130            .mls(&mls)
131            .build()
132            .expect("Unable to construct a `PureOrder` structure with decreasing-m order.")
133    }
134
135    /// Constructs a new [`PureOrder`] structure for a specified rank with Molden order. The spatial
136    /// inversion parity is determined by the parity of the deduced $`l`$ value.
137    ///
138    /// # Arguments
139    ///
140    /// * `lpure` - The required pure Gaussian rank.
141    ///
142    /// # Returns
143    ///
144    /// A [`PureOrder`] struct for a specified rank with Molden order.
145    #[must_use]
146    pub fn molden(lpure: u32) -> Self {
147        let lpure_i32 = i32::try_from(lpure).expect("`lpure` cannot be converted to `i32`.");
148        let mls = (0..=lpure_i32)
149            .flat_map(|absm| {
150                if absm == 0 {
151                    vec![0]
152                } else {
153                    vec![absm, -absm]
154                }
155            })
156            .collect_vec();
157        Self::builder()
158            .lpure(lpure)
159            .mls(&mls)
160            .build()
161            .expect("Unable to construct a `PureOrder` structure with Molden order.")
162    }
163
164    /// Verifies if this [`PureOrder`] struct is valid.
165    ///
166    /// # Returns
167    ///
168    /// A boolean indicating if this [`PureOrder`] struct is valid.
169    #[must_use]
170    pub fn verify(&self) -> bool {
171        let mls_set = self.mls.iter().collect::<HashSet<_>>();
172        let lpure = self.lpure;
173        mls_set.len() == self.ncomps() && mls_set.iter().all(|m| m.unsigned_abs() <= lpure)
174    }
175
176    /// Iterates over the constituent $`m_l`$ values.
177    pub fn iter(&self) -> Iter<i32> {
178        self.mls.iter()
179    }
180
181    /// Returns the number of pure components in the shell.
182    pub fn ncomps(&self) -> usize {
183        let lpure = usize::try_from(self.lpure).unwrap_or_else(|_| {
184            panic!(
185                "Unable to convert the pure degree {} to `usize`.",
186                self.lpure
187            )
188        });
189        2 * lpure + 1
190    }
191
192    /// Returns the $`m`$ value with a specified index in this shell.
193    pub fn get_m_with_index(&self, i: usize) -> Option<i32> {
194        self.mls.get(i).cloned()
195    }
196}
197
198impl fmt::Display for PureOrder {
199    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
200        writeln!(f, "Pure rank: {}", self.lpure)?;
201        writeln!(f, "Order:")?;
202        for m in self.iter() {
203            writeln!(f, "  {m}")?;
204        }
205        Ok(())
206    }
207}
208
209impl fmt::Debug for PureOrder {
210    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
211        writeln!(f, "Pure rank: {}", self.lpure)?;
212        writeln!(f, "Order:")?;
213        for m in self.iter() {
214            writeln!(f, "  {m:?}")?;
215        }
216        Ok(())
217    }
218}
219
220impl PermutableCollection for PureOrder {
221    type Rank = usize;
222
223    fn get_perm_of(&self, other: &Self) -> Option<Permutation<Self::Rank>> {
224        let o_mls: HashMap<&i32, usize> = other
225            .mls
226            .iter()
227            .enumerate()
228            .map(|(i, o_m)| (o_m, i))
229            .collect();
230        let image_opt: Option<Vec<Self::Rank>> =
231            self.mls.iter().map(|s_m| o_mls.get(s_m).copied()).collect();
232        image_opt.and_then(|image| Permutation::from_image(image).ok())
233    }
234
235    fn permute(&self, perm: &Permutation<Self::Rank>) -> Result<Self, anyhow::Error> {
236        let mut p_pureorder = self.clone();
237        p_pureorder.permute_mut(perm)?;
238        Ok(p_pureorder)
239    }
240
241    fn permute_mut(&mut self, perm: &Permutation<Self::Rank>) -> Result<(), anyhow::Error> {
242        permute_inplace(&mut self.mls, perm);
243        Ok(())
244    }
245}
246
247// ~~~~~~~~~
248// CartOrder
249// ~~~~~~~~~
250
251/// Structure to contain information about the ordering of Cartesian Gaussians of a certain rank.
252#[derive(Clone, Builder, PartialEq, Eq, Hash)]
253pub struct CartOrder {
254    /// A sequence of $`(l_x, l_y, l_z)`$ tuples giving the ordering of the Cartesian Gaussians.
255    #[builder(setter(custom))]
256    pub cart_tuples: Vec<(u32, u32, u32)>,
257
258    /// The rank of the Cartesian Gaussians.
259    pub lcart: u32,
260}
261
262impl CartOrderBuilder {
263    fn cart_tuples(&mut self, cart_tuples: &[(u32, u32, u32)]) -> &mut Self {
264        let lcart = self.lcart.expect("`lcart` has not been set.");
265        assert!(
266            cart_tuples.iter().all(|(lx, ly, lz)| lx + ly + lz == lcart),
267            "Inconsistent total Cartesian orders between components."
268        );
269        assert_eq!(
270            cart_tuples.len(),
271            ((lcart + 1) * (lcart + 2)).div_euclid(2) as usize,
272            "Unexpected number of components for `lcart` = {}.",
273            lcart
274        );
275        self.cart_tuples = Some(cart_tuples.to_vec());
276        self
277    }
278}
279
280impl CartOrder {
281    /// Returns a builder to construct a new [`CartOrder`] structure.
282    fn builder() -> CartOrderBuilder {
283        CartOrderBuilder::default()
284    }
285
286    /// Constructs a new [`CartOrder`] structure from its constituting tuples, each of which contains
287    /// the $`x`$, $`y`$, and $`z`$ exponents for one Cartesian term.
288    ///
289    /// # Errors
290    ///
291    /// Errors if the Cartesian tuples are invalid (*e.g.* missing components or containing
292    /// inconsistent components).
293    pub fn new(cart_tuples: &[(u32, u32, u32)]) -> Result<Self, anyhow::Error> {
294        let first_tuple = cart_tuples
295            .get(0)
296            .ok_or(format_err!("No Cartesian tuples found."))?;
297        let lcart = first_tuple.0 + first_tuple.1 + first_tuple.2;
298        let cart_order = CartOrder::builder()
299            .lcart(lcart)
300            .cart_tuples(cart_tuples)
301            .build()
302            .map_err(|err| format_err!(err))?;
303        ensure!(cart_order.verify(), "Invalid `CartOrder`.");
304        Ok(cart_order)
305    }
306
307    /// Constructs a new [`CartOrder`] structure for a specified rank with lexicographic order.
308    ///
309    /// # Arguments
310    ///
311    /// * `lcart` - The required Cartesian Gaussian rank.
312    ///
313    /// # Returns
314    ///
315    /// A [`CartOrder`] struct for a specified rank with lexicographic order.
316    #[must_use]
317    pub fn lex(lcart: u32) -> Self {
318        let mut cart_tuples =
319            Vec::with_capacity(((lcart + 1) * (lcart + 2)).div_euclid(2) as usize);
320        for lx in (0..=lcart).rev() {
321            for ly in (0..=(lcart - lx)).rev() {
322                cart_tuples.push((lx, ly, lcart - lx - ly));
323            }
324        }
325        Self::builder()
326            .lcart(lcart)
327            .cart_tuples(&cart_tuples)
328            .build()
329            .expect("Unable to construct a `CartOrder` structure with lexicographic order.")
330    }
331
332    /// Constructs a new [`CartOrder`] structure for a specified rank with Molden order.
333    ///
334    /// # Arguments
335    ///
336    /// * `lcart` - The required Cartesian Gaussian rank up to 4.
337    ///
338    /// # Returns
339    ///
340    /// A [`CartOrder`] struct for a specified rank with Q-Chem order.
341    ///
342    /// # Panics
343    ///
344    /// Panics if `lcart` is greater than 4.
345    #[must_use]
346    pub fn molden(lcart: u32) -> Self {
347        assert!(lcart <= 4, "`lcart` > 4 is not specified by Molden.");
348        let cart_tuples: Vec<(u32, u32, u32)> = match lcart {
349            0 => vec![(0, 0, 0)],
350            1 => vec![(1, 0, 0), (0, 1, 0), (0, 0, 1)],
351            2 => vec![
352                (2, 0, 0),
353                (0, 2, 0),
354                (0, 0, 2),
355                (1, 1, 0),
356                (1, 0, 1),
357                (0, 1, 1),
358            ],
359            3 => vec![
360                (3, 0, 0),
361                (0, 3, 0),
362                (0, 0, 3),
363                (1, 2, 0),
364                (2, 1, 0),
365                (2, 0, 1),
366                (1, 0, 2),
367                (0, 1, 2),
368                (0, 2, 1),
369                (1, 1, 1),
370            ],
371            4 => vec![
372                (4, 0, 0),
373                (0, 4, 0),
374                (0, 0, 4),
375                (3, 1, 0),
376                (3, 0, 1),
377                (1, 3, 0),
378                (0, 3, 1),
379                (1, 0, 3),
380                (0, 1, 3),
381                (2, 2, 0),
382                (2, 0, 2),
383                (0, 2, 2),
384                (2, 1, 1),
385                (1, 2, 1),
386                (1, 1, 2),
387            ],
388            _ => panic!("`lcart` > 4 is not specified by Molden."),
389        };
390        Self::builder()
391            .lcart(lcart)
392            .cart_tuples(&cart_tuples)
393            .build()
394            .expect("Unable to construct a `CartOrder` structure with Molden order.")
395    }
396
397    /// Constructs a new [`CartOrder`] structure for a specified rank with Q-Chem order.
398    ///
399    /// # Arguments
400    ///
401    /// * `lcart` - The required Cartesian Gaussian rank.
402    ///
403    /// # Returns
404    ///
405    /// A [`CartOrder`] struct for a specified rank with Q-Chem order.
406    #[must_use]
407    pub fn qchem(lcart: u32) -> Self {
408        let cart_tuples: Vec<(u32, u32, u32)> = if lcart > 0 {
409            (0..3)
410                .product_repeat(lcart as usize)
411                .filter_map(|tup| {
412                    let mut tup_sorted = tup.clone();
413                    tup_sorted.sort_unstable();
414                    tup_sorted.reverse();
415                    if tup == tup_sorted {
416                        let lcartqns = tup.iter().collect::<Counter<_>>();
417                        Some((
418                            <usize as TryInto<u32>>::try_into(*(lcartqns.get(&0).unwrap_or(&0)))
419                                .expect("Unable to convert Cartesian x-exponent to `u32`."),
420                            <usize as TryInto<u32>>::try_into(*(lcartqns.get(&1).unwrap_or(&0)))
421                                .expect("Unable to convert Cartesian y-exponent to `u32`."),
422                            <usize as TryInto<u32>>::try_into(*(lcartqns.get(&2).unwrap_or(&0)))
423                                .expect("Unable to convert Cartesian z-exponent to `u32`."),
424                        ))
425                    } else {
426                        None
427                    }
428                })
429                .collect()
430        } else {
431            vec![(0, 0, 0)]
432        };
433        Self::builder()
434            .lcart(lcart)
435            .cart_tuples(&cart_tuples)
436            .build()
437            .expect("Unable to construct a `CartOrder` structure with Q-Chem order.")
438    }
439
440    /// Verifies if this [`CartOrder`] struct is valid.
441    ///
442    /// # Returns
443    ///
444    /// A boolean indicating if this [`CartOrder`] struct is valid.
445    #[must_use]
446    pub fn verify(&self) -> bool {
447        let cart_tuples_set = self.cart_tuples.iter().collect::<HashSet<_>>();
448        let lcart = self.lcart;
449        cart_tuples_set.len() == self.ncomps()
450            && cart_tuples_set
451                .iter()
452                .all(|(lx, ly, lz)| lx + ly + lz == lcart)
453    }
454
455    /// Iterates over the constituent tuples.
456    pub fn iter(&self) -> Iter<(u32, u32, u32)> {
457        self.cart_tuples.iter()
458    }
459
460    /// Returns the number of Cartesian components in the shell.
461    pub fn ncomps(&self) -> usize {
462        let lcart = usize::try_from(self.lcart).unwrap_or_else(|_| {
463            panic!(
464                "Unable to convert the Cartesian degree {} to `usize`.",
465                self.lcart
466            )
467        });
468        ((lcart + 1) * (lcart + 2)).div_euclid(2)
469    }
470
471    /// Returns the Cartesian component with a specified index in this shell.
472    pub fn get_cart_tuple_with_index(&self, i: usize) -> Option<(u32, u32, u32)> {
473        self.cart_tuples.get(i).cloned()
474    }
475}
476
477impl fmt::Display for CartOrder {
478    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
479        writeln!(f, "Cartesian rank: {}", self.lcart)?;
480        writeln!(f, "Order:")?;
481        for cart_tuple in self.iter() {
482            writeln!(f, "  {}", cart_tuple_to_str(cart_tuple, true))?;
483        }
484        Ok(())
485    }
486}
487
488impl fmt::Debug for CartOrder {
489    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
490        writeln!(f, "Cartesian rank: {}", self.lcart)?;
491        writeln!(f, "Order:")?;
492        for cart_tuple in self.iter() {
493            writeln!(f, "  {cart_tuple:?}")?;
494        }
495        Ok(())
496    }
497}
498
499impl PermutableCollection for CartOrder {
500    type Rank = usize;
501
502    fn get_perm_of(&self, other: &Self) -> Option<Permutation<Self::Rank>> {
503        let o_cart_tuples: HashMap<&(u32, u32, u32), usize> = other
504            .cart_tuples
505            .iter()
506            .enumerate()
507            .map(|(i, o_cart_tuple)| (o_cart_tuple, i))
508            .collect();
509        let image_opt: Option<Vec<Self::Rank>> = self
510            .cart_tuples
511            .iter()
512            .map(|s_cart_tuple| o_cart_tuples.get(s_cart_tuple).copied())
513            .collect();
514        image_opt.and_then(|image| Permutation::from_image(image).ok())
515    }
516
517    fn permute(&self, perm: &Permutation<Self::Rank>) -> Result<Self, anyhow::Error> {
518        let mut p_cartorder = self.clone();
519        p_cartorder.permute_mut(perm)?;
520        Ok(p_cartorder)
521    }
522
523    fn permute_mut(&mut self, perm: &Permutation<Self::Rank>) -> Result<(), anyhow::Error> {
524        permute_inplace(&mut self.cart_tuples, perm);
525        Ok(())
526    }
527}
528
529/// Translates a Cartesian exponent tuple to a human-understandable string.
530///
531/// # Arguments
532///
533/// * `cart_tuple` - A tuple of $`(l_x, l_y, l_z)`$ specifying the exponents of the Cartesian
534/// components of the Cartesian Gaussian.
535/// * flat - A flag indicating if the string representation is flat (*e.g.* `xxyz`) or compact
536/// (*e.g.* `x^2yz`).
537///
538/// Returns
539///
540/// The string representation of the Cartesian exponent tuple.
541pub(crate) fn cart_tuple_to_str(cart_tuple: &(u32, u32, u32), flat: bool) -> String {
542    if cart_tuple.0 + cart_tuple.1 + cart_tuple.2 == 0u32 {
543        "1".to_string()
544    } else {
545        let cart_array = [cart_tuple.0, cart_tuple.1, cart_tuple.2];
546        let carts = ["x", "y", "z"];
547        Itertools::intersperse(
548            cart_array.iter().enumerate().map(|(i, &l)| {
549                if flat {
550                    carts[i].repeat(l as usize)
551                } else {
552                    match l.cmp(&1) {
553                        Ordering::Greater => format!("{}^{l}", carts[i]),
554                        Ordering::Equal => carts[i].to_string(),
555                        Ordering::Less => String::new(),
556                    }
557                }
558            }),
559            String::new(),
560        )
561        .collect::<String>()
562    }
563}
564
565// ~~~~~~~~~~~
566// SpinorOrder
567// ~~~~~~~~~~~
568
569/// Structure to contain information about the ordering of spinors of a certain rank.
570#[derive(Clone, Builder, PartialEq, Eq, Hash)]
571pub struct SpinorOrder {
572    /// The angular momentum (times two) of the spinor Gaussians. This must be an odd integer.
573    #[builder(setter(custom))]
574    pub two_j: u32,
575
576    /// A sequence of $`2m_j`$ values giving the ordering of the spinor Gaussians.
577    #[builder(setter(custom))]
578    two_mjs: Vec<i32>,
579
580    /// The spatial inversion parity of the spinor Gaussians: `true` if even under spatial inversion
581    /// and `false` if odd.
582    pub even: bool,
583}
584
585impl SpinorOrderBuilder {
586    fn two_j(&mut self, two_j: u32) -> &mut Self {
587        if two_j.rem_euclid(2) != 1 {
588            panic!("`two_j` must be odd.")
589        }
590        self.two_j = Some(two_j);
591        self
592    }
593
594    fn two_mjs(&mut self, two_mjs: &[i32]) -> &mut Self {
595        let two_j = self.two_j.expect("`two_j` has not been set.");
596        let max_two_m = two_mjs
597            .iter()
598            .map(|two_m| two_m.unsigned_abs())
599            .max()
600            .expect("The maximum |2m| value could not be determined.");
601        assert_eq!(
602            max_two_m, two_j,
603            "The maximum |2m| value does not equal 2j = {two_j}."
604        );
605        assert_eq!(
606            two_mjs.len(),
607            (two_j + 1) as usize,
608            "The number of 2m values specified does not equal 2j + 1 = {}",
609            two_j + 1
610        );
611        self.two_mjs = Some(two_mjs.to_vec());
612        self
613    }
614}
615
616impl SpinorOrder {
617    /// Returns a builder to construct a new [`SpinorOrder`] structure.
618    fn builder() -> SpinorOrderBuilder {
619        SpinorOrderBuilder::default()
620    }
621
622    /// Constructs a new [`SpinorOrder`] structure from its constituting $`2m_j`$ values and a
623    /// specified spatial parity.
624    pub fn new(two_mjs: &[i32], even: bool) -> Result<Self, anyhow::Error> {
625        let two_j = two_mjs
626            .iter()
627            .map(|two_m| two_m.unsigned_abs())
628            .max()
629            .ok_or_else(|| format_err!("The maximum |2m| value could not be determined."))?;
630        let spinor_order = SpinorOrder::builder()
631            .two_j(two_j)
632            .two_mjs(two_mjs)
633            .even(even)
634            .build()
635            .map_err(|err| format_err!(err))?;
636        ensure!(spinor_order.verify(), "Invalid `SpinorOrder`.");
637        Ok(spinor_order)
638    }
639
640    /// Constructs a new [`SpinorOrder`] structure for a specified angular momentum with
641    /// increasing-$`m`$ order.
642    ///
643    /// # Arguments
644    ///
645    /// * `two_j` - The required spinor angular momentum (times two).
646    /// * `even` - Boolean indicating whether the spinors are even with respect to spatial
647    /// inversion.
648    ///
649    /// # Returns
650    ///
651    /// A [`SpinorOrder`] struct for a specified angular momentum with increasing-$`m`$ order.
652    #[must_use]
653    pub fn increasingm(two_j: u32, even: bool) -> Self {
654        let two_j_i32 = i32::try_from(two_j).expect("`two_j` cannot be converted to `i32`.");
655        let two_mjs = (-two_j_i32..=two_j_i32).step_by(2).collect_vec();
656        Self::builder()
657            .two_j(two_j)
658            .two_mjs(&two_mjs)
659            .even(even)
660            .build()
661            .expect("Unable to construct a `SpinorOrder` structure with increasing-m order.")
662    }
663
664    /// Constructs a new [`SpinorOrder`] structure for a specified angular momentum with
665    /// decreasing-$`m`$ order.
666    ///
667    /// # Arguments
668    ///
669    /// * `two_j` - The required spinor angular momentum (times two).
670    /// * `even` - Boolean indicating whether the spinors are even with respect to spatial
671    /// inversion.
672    ///
673    /// # Returns
674    ///
675    /// A [`SpinorOrder`] struct for a specified angular momentum with decreasing-$`m`$ order.
676    #[must_use]
677    pub fn decreasingm(two_j: u32, even: bool) -> Self {
678        let two_j_i32 = i32::try_from(two_j).expect("`two_j` cannot be converted to `i32`.");
679        let two_mjs = (-two_j_i32..=two_j_i32).rev().step_by(2).collect_vec();
680        Self::builder()
681            .two_j(two_j)
682            .two_mjs(&two_mjs)
683            .even(even)
684            .build()
685            .expect("Unable to construct a `SpinorOrder` structure with decreasing-m order.")
686    }
687
688    /// Constructs a new [`SpinorOrder`] structure for a specified angular momentum with Molden order.
689    ///
690    /// # Arguments
691    ///
692    /// * `two_j` - The required spinor angular momentum (times two).
693    /// * `even` - Boolean indicating whether the spinors are even with respect to spatial
694    /// inversion.
695    ///
696    /// # Returns
697    ///
698    /// A [`SpinorOrder`] struct for a specified angular momentum with Molden order.
699    #[must_use]
700    pub fn molden(two_j: u32, even: bool) -> Self {
701        let two_j_i32 = i32::try_from(two_j).expect("`two_j` cannot be converted to `i32`.");
702        let two_mjs = (1..=two_j_i32)
703            .step_by(2)
704            .flat_map(|abs_twom| vec![abs_twom, -abs_twom])
705            .collect_vec();
706        Self::builder()
707            .two_j(two_j)
708            .two_mjs(&two_mjs)
709            .even(even)
710            .build()
711            .expect("Unable to construct a `SpinorOrder` structure with Molden order.")
712    }
713
714    /// Verifies if this [`SpinorOrder`] struct is valid.
715    ///
716    /// # Returns
717    ///
718    /// A boolean indicating if this [`SpinorOrder`] struct is valid.
719    #[must_use]
720    pub fn verify(&self) -> bool {
721        let two_mjs_set = self.two_mjs.iter().collect::<HashSet<_>>();
722        let two_j = self.two_j;
723        two_mjs_set.len() == self.ncomps()
724            && two_mjs_set
725                .iter()
726                .all(|two_m| two_m.unsigned_abs() <= two_j)
727    }
728
729    /// Iterates over the constituent $`2m_j`$ values.
730    pub fn iter(&self) -> Iter<i32> {
731        self.two_mjs.iter()
732    }
733
734    /// Returns the number of pure components in the shell.
735    pub fn ncomps(&self) -> usize {
736        let two_j = usize::try_from(self.two_j).unwrap_or_else(|_| {
737            panic!(
738                "Unable to convert the two-j value {} to `usize`.",
739                self.two_j
740            )
741        });
742        two_j + 1
743    }
744
745    /// Returns the $`2m`$ value with a specified index in this shell.
746    pub fn get_two_m_with_index(&self, i: usize) -> Option<i32> {
747        self.two_mjs.get(i).cloned()
748    }
749}
750
751impl fmt::Display for SpinorOrder {
752    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
753        writeln!(
754            f,
755            "Angular momentum: {}/2 ({})",
756            self.two_j,
757            if self.even { "g" } else { "u" }
758        )?;
759        writeln!(f, "Order:")?;
760        for two_m in self.iter() {
761            writeln!(f, "  {two_m}/2")?;
762        }
763        Ok(())
764    }
765}
766
767impl fmt::Debug for SpinorOrder {
768    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
769        writeln!(
770            f,
771            "Angular momentum: {}/2 ({})",
772            self.two_j,
773            if self.even { "g" } else { "u" }
774        )?;
775        writeln!(f, "Order:")?;
776        for two_m in self.iter() {
777            writeln!(f, "  {two_m:?}/2")?;
778        }
779        Ok(())
780    }
781}
782
783impl PermutableCollection for SpinorOrder {
784    type Rank = usize;
785
786    fn get_perm_of(&self, other: &Self) -> Option<Permutation<Self::Rank>> {
787        let o_two_mjs: HashMap<&i32, usize> = other
788            .two_mjs
789            .iter()
790            .enumerate()
791            .map(|(i, o_two_m)| (o_two_m, i))
792            .collect();
793        let image_opt: Option<Vec<Self::Rank>> = self
794            .two_mjs
795            .iter()
796            .map(|s_two_m| o_two_mjs.get(s_two_m).copied())
797            .collect();
798        image_opt.and_then(|image| Permutation::from_image(image).ok())
799    }
800
801    fn permute(&self, perm: &Permutation<Self::Rank>) -> Result<Self, anyhow::Error> {
802        let mut p_pureorder = self.clone();
803        p_pureorder.permute_mut(perm)?;
804        Ok(p_pureorder)
805    }
806
807    fn permute_mut(&mut self, perm: &Permutation<Self::Rank>) -> Result<(), anyhow::Error> {
808        permute_inplace(&mut self.two_mjs, perm);
809        Ok(())
810    }
811}
812
813// ----------
814// ShellOrder
815// ----------
816
817/// Enumerated type to indicate the type of the angular functions in a shell and how they are
818/// ordered.
819#[derive(Clone, PartialEq, Eq, Hash, Debug)]
820pub enum ShellOrder {
821    /// This variant indicates that the angular functions are real solid harmonics. The associated
822    /// value is a [`PureOrder`] struct containing the order of these functions.
823    Pure(PureOrder),
824
825    /// This variant indicates that the angular functions are Cartesian functions. The associated
826    /// value is a [`CartOrder`] struct containing the order of these functions.
827    Cart(CartOrder),
828
829    /// This variant indicates that the angular functions are spinors. The associated
830    /// value is a [`SpinorOrder`] struct containing the order of these functions.
831    Spinor(SpinorOrder),
832}
833
834impl fmt::Display for ShellOrder {
835    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
836        match self {
837            ShellOrder::Pure(pure_order) => write!(
838                f,
839                "Pure ({}) ({})",
840                if pure_order.even { "g" } else { "u" },
841                pure_order.iter().map(|m| m.to_string()).join(", ")
842            ),
843            ShellOrder::Cart(cart_order) => write!(
844                f,
845                "Cart ({}) ({})",
846                if cart_order.lcart % 2 == 0 { "g" } else { "u" },
847                cart_order
848                    .iter()
849                    .map(|cart_tuple| { cart_tuple_to_str(cart_tuple, true) })
850                    .join(", ")
851            ),
852            ShellOrder::Spinor(spinor_order) => write!(
853                f,
854                "Spinor ({}) ({})",
855                if spinor_order.even { "g" } else { "u" },
856                spinor_order
857                    .iter()
858                    .map(|two_m| format!("{two_m}/2"))
859                    .join(", ")
860            ),
861        }
862    }
863}
864
865// ----------
866// BasisShell
867// ----------
868
869/// Structure representing a shell in an atomic-orbital basis set.
870#[derive(Clone, Builder, PartialEq, Eq, Hash, Debug)]
871pub struct BasisShell {
872    /// A non-negative integer indicating the rank of the shell. If this shell is pure, the rank
873    /// equals its angular momentum. If this shell is Cartesian, the rank is the sum of the
874    /// exponents of the Cartesian coordinates. If this shell is a spinor, the rank is twice its
875    /// angular momentum.
876    #[builder(setter(custom))]
877    pub l: u32,
878
879    /// An enum indicating the type of the angular functions in a shell and how they are ordered.
880    #[builder(setter(custom))]
881    pub shell_order: ShellOrder,
882}
883
884impl BasisShellBuilder {
885    fn l(&mut self, l: u32) -> &mut Self {
886        match self.shell_order.as_ref() {
887            Some(ShellOrder::Pure(pure_order)) => assert_eq!(pure_order.lpure, l),
888            Some(ShellOrder::Cart(cart_order)) => assert_eq!(cart_order.lcart, l),
889            Some(ShellOrder::Spinor(spinor_order)) => assert_eq!(spinor_order.two_j, l),
890            None => {}
891        }
892        self.l = Some(l);
893        self
894    }
895
896    fn shell_order(&mut self, shl_ord: ShellOrder) -> &mut Self {
897        match (&shl_ord, self.l) {
898            (ShellOrder::Pure(pure_order), Some(l)) => assert_eq!(pure_order.lpure, l),
899            (ShellOrder::Cart(cart_order), Some(l)) => assert_eq!(cart_order.lcart, l),
900            (ShellOrder::Spinor(spinor_order), Some(l)) => assert_eq!(spinor_order.two_j, l),
901            _ => {}
902        }
903        self.shell_order = Some(shl_ord);
904        self
905    }
906}
907
908impl BasisShell {
909    /// Returns a builder to construct a new [`BasisShell`].
910    ///
911    /// # Returns
912    ///
913    /// A builder to construct a new [`BasisShell`].
914    fn builder() -> BasisShellBuilder {
915        BasisShellBuilder::default()
916    }
917
918    /// Constructs a new [`BasisShell`].
919    ///
920    /// # Arguments
921    ///
922    /// * `l` - The rank of this shell, which is equal to $`l_{\mathrm{pure}}`$ for a pure shell,
923    /// $`l_{\mathrm{cart}}`$ for a Cartesian shell, or $`2j`$ for a spinor shell.
924    /// * `shl_ord` - A [`ShellOrder`] structure specifying the type and ordering of the basis
925    /// functions in this shell.
926    pub fn new(l: u32, shl_ord: ShellOrder) -> Self {
927        match &shl_ord {
928            ShellOrder::Pure(pureorder) => assert_eq!(pureorder.lpure, l),
929            ShellOrder::Cart(cartorder) => assert_eq!(cartorder.lcart, l),
930            ShellOrder::Spinor(spinororder) => assert_eq!(spinororder.two_j, l),
931        }
932        BasisShell::builder()
933            .l(l)
934            .shell_order(shl_ord)
935            .build()
936            .expect("Unable to construct a `BasisShell`.")
937    }
938
939    /// The number of basis functions in this shell.
940    pub fn n_funcs(&self) -> usize {
941        let lsize = self.l as usize;
942        match self.shell_order {
943            ShellOrder::Pure(_) => 2 * lsize + 1,
944            ShellOrder::Cart(_) => ((lsize + 1) * (lsize + 2)).div_euclid(2),
945            ShellOrder::Spinor(_) => lsize + 1,
946        }
947    }
948}
949
950// ---------
951// BasisAtom
952// ---------
953
954/// Structure containing the ordered sequence of the shells for an atom.
955#[derive(Clone, Builder, PartialEq, Eq, Hash, Debug)]
956pub struct BasisAtom<'a> {
957    /// An atom in the basis set.
958    pub(crate) atom: &'a Atom,
959
960    /// The ordered shells associated with this atom.
961    #[builder(setter(custom))]
962    pub(crate) basis_shells: Vec<BasisShell>,
963}
964
965impl<'a> BasisAtomBuilder<'a> {
966    pub(crate) fn basis_shells(&mut self, bss: &[BasisShell]) -> &mut Self {
967        self.basis_shells = Some(bss.to_vec());
968        self
969    }
970}
971
972impl<'a> BasisAtom<'a> {
973    /// Returns a builder to construct a new [`BasisAtom`].
974    ///
975    /// # Returns
976    ///
977    /// A builder to construct a new [`BasisAtom`].
978    pub(crate) fn builder() -> BasisAtomBuilder<'a> {
979        BasisAtomBuilder::default()
980    }
981
982    /// Constructs a new [`BasisAtom`].
983    ///
984    /// # Arguments
985    ///
986    /// * `atm` - A reference to an atom.
987    /// * `bss` - A sequence of [`BasisShell`]s containing the basis functions localised on this
988    /// atom.
989    pub fn new(atm: &'a Atom, bss: &[BasisShell]) -> Self {
990        BasisAtom::builder()
991            .atom(atm)
992            .basis_shells(bss)
993            .build()
994            .expect("Unable to construct a `BasisAtom`.")
995    }
996
997    /// The number of basis functions localised on this atom.
998    fn n_funcs(&self) -> usize {
999        self.basis_shells.iter().map(BasisShell::n_funcs).sum()
1000    }
1001
1002    /// The ordered tuples of 0-based indices indicating the starting (inclusive) and ending
1003    /// (exclusive) positions of the shells on this atom.
1004    fn shell_boundary_indices(&self) -> Vec<(usize, usize)> {
1005        self.basis_shells
1006            .iter()
1007            .scan(0, |acc, basis_shell| {
1008                let start_index = *acc;
1009                *acc += basis_shell.n_funcs();
1010                Some((start_index, *acc))
1011            })
1012            .collect::<Vec<_>>()
1013    }
1014}
1015
1016// -----------------
1017// BasisAngularOrder
1018// -----------------
1019
1020/// Structure containing the angular momentum information of an atomic-orbital basis set that is
1021/// required for symmetry transformation to be performed.
1022#[derive(Clone, Builder, PartialEq, Eq, Hash, Debug)]
1023pub struct BasisAngularOrder<'a> {
1024    /// An ordered sequence of [`BasisAtom`] in the order the atoms are defined in the molecule.
1025    #[builder(setter(custom))]
1026    pub(crate) basis_atoms: Vec<BasisAtom<'a>>,
1027}
1028
1029impl<'a> BasisAngularOrderBuilder<'a> {
1030    pub(crate) fn basis_atoms(&mut self, batms: &[BasisAtom<'a>]) -> &mut Self {
1031        self.basis_atoms = Some(batms.to_vec());
1032        self
1033    }
1034}
1035
1036impl<'a> BasisAngularOrder<'a> {
1037    /// Returns a builder to construct a new [`BasisAngularOrder`].
1038    ///
1039    /// # Returns
1040    ///
1041    /// A builder to construct a new [`BasisAngularOrder`].
1042    #[must_use]
1043    pub(crate) fn builder() -> BasisAngularOrderBuilder<'a> {
1044        BasisAngularOrderBuilder::default()
1045    }
1046
1047    /// Constructs a new [`BasisAngularOrder`] structure from the constituting [`BasisAtom`]s.
1048    ///
1049    /// # Arguments
1050    ///
1051    /// * `batms` - The constituent [`BasisAtom`]s.
1052    pub fn new(batms: &[BasisAtom<'a>]) -> Self {
1053        BasisAngularOrder::builder()
1054            .basis_atoms(batms)
1055            .build()
1056            .expect("Unable to construct a `BasisAngularOrder`.")
1057    }
1058
1059    /// The number of atoms in the basis.
1060    pub fn n_atoms(&self) -> usize {
1061        self.basis_atoms.len()
1062    }
1063
1064    /// The number of basis functions in this basis.
1065    pub fn n_funcs(&self) -> usize {
1066        self.basis_atoms.iter().map(BasisAtom::n_funcs).sum()
1067    }
1068
1069    /// The ordered tuples of 0-based shell indices indicating the starting (inclusive) and ending
1070    /// (exclusive) shell positions of the atoms in this basis.
1071    pub fn atom_boundary_indices(&self) -> Vec<(usize, usize)> {
1072        self.basis_atoms
1073            .iter()
1074            .scan(0, |acc, basis_atom| {
1075                let start_index = *acc;
1076                *acc += basis_atom.n_funcs();
1077                Some((start_index, *acc))
1078            })
1079            .collect::<Vec<_>>()
1080    }
1081
1082    /// The ordered tuples of 0-based function indices indicating the starting (inclusive) and
1083    /// ending (exclusive) positions of the shells in this basis.
1084    pub fn shell_boundary_indices(&self) -> Vec<(usize, usize)> {
1085        let atom_boundary_indices = self.atom_boundary_indices();
1086        self.basis_atoms
1087            .iter()
1088            .zip(atom_boundary_indices)
1089            .flat_map(|(basis_atom, (atom_start, _))| {
1090                basis_atom
1091                    .shell_boundary_indices()
1092                    .iter()
1093                    .map(|(shell_start, shell_end)| {
1094                        (shell_start + atom_start, shell_end + atom_start)
1095                    })
1096                    .collect::<Vec<_>>()
1097            })
1098            .collect::<Vec<_>>()
1099    }
1100
1101    /// An iterator over the constituent [`BasisShell`]s in this basis.
1102    pub fn basis_shells(&self) -> impl Iterator<Item = &BasisShell> + '_ {
1103        self.basis_atoms
1104            .iter()
1105            .flat_map(|basis_atom| basis_atom.basis_shells.iter())
1106    }
1107
1108    /// Determines the permutation of the functions in this [`BasisAngularOrder`] to map `self` to
1109    /// `other`, given that the shells themselves remain unchanged while only the functions in each
1110    /// shell are permuted.
1111    ///
1112    /// For example, consider `self`:
1113    /// ```text
1114    /// S (1)
1115    /// P (x, y, z)
1116    /// D (xx, xy, xz, yy, yz, zz)
1117    /// ```
1118    ///
1119    /// and `other`:
1120    /// ```text
1121    /// S (1)
1122    /// P (y, z, x)
1123    /// D (xx, xy, yy, xz, yz, zz)
1124    /// ```
1125    ///
1126    /// the mapping permutation is given by `π(0, 3, 1, 2, 4, 5, 7, 6, 8, 9)`.
1127    ///
1128    /// # Arguments
1129    ///
1130    /// * `other` - Another [`BasisAngularOrder`] to be compared against.
1131    ///
1132    /// # Returns
1133    ///
1134    /// The mapping permutation, if any.
1135    pub(crate) fn get_perm_of_functions_fixed_shells(
1136        &self,
1137        other: &Self,
1138    ) -> Result<Permutation<usize>, anyhow::Error> {
1139        if self.n_funcs() == other.n_funcs() && self.n_atoms() == other.n_atoms() {
1140            let s_shell_boundaries = self.shell_boundary_indices();
1141            let o_shell_boundaries = other.shell_boundary_indices();
1142            if s_shell_boundaries.len() == o_shell_boundaries.len() {
1143                let image = izip!(
1144                    self.basis_shells(),
1145                    other.basis_shells(),
1146                    s_shell_boundaries.iter(),
1147                    o_shell_boundaries.iter()
1148                )
1149                .map(|(s_bs, o_bs, (s_start, s_end), (o_start, o_end))| {
1150                    if (s_start, s_end) == (o_start, o_end) {
1151                        let s_shl_ord = &s_bs.shell_order;
1152                        let o_shl_ord = &o_bs.shell_order;
1153                        match (s_shl_ord, o_shl_ord) {
1154                            (ShellOrder::Pure(s_po), ShellOrder::Pure(o_po)) => Ok(
1155                                s_po.get_perm_of(&o_po)
1156                                    .unwrap()
1157                                    .image()
1158                                    .iter()
1159                                    .map(|x| s_start + x)
1160                                    .collect_vec(),
1161                            ),
1162                            (ShellOrder::Cart(s_co), ShellOrder::Cart(o_co)) => Ok(
1163                                s_co.get_perm_of(&o_co)
1164                                    .unwrap()
1165                                    .image()
1166                                    .iter()
1167                                    .map(|x| s_start + x)
1168                                    .collect_vec(),
1169                            ),
1170                            _ => Err(format_err!("At least one pair of corresponding shells have mismatched pure/cart.")),
1171                        }
1172                    } else {
1173                        Err(format_err!("At least one pair of corresponding shells have mismatched boundary indices."))
1174                    }
1175                })
1176                .collect::<Result<Vec<_>, _>>()
1177                .and_then(|image_by_shells| {
1178                    let flattened_image = image_by_shells.into_iter().flatten().collect_vec();
1179                    Permutation::from_image(flattened_image)
1180                });
1181                image
1182            } else {
1183                Err(format_err!("Mismatched numbers of shells."))
1184            }
1185        } else {
1186            Err(format_err!("Mismatched numbers of basis functions."))
1187        }
1188    }
1189}
1190
1191impl<'a> PermutableCollection for BasisAngularOrder<'a> {
1192    type Rank = usize;
1193
1194    /// Determines the permutation of [`BasisAtom`]s to map `self` to `other`.
1195    ///
1196    /// # Arguments
1197    ///
1198    /// * `other` - Another [`BasisAngularOrder`] to be compared with `self`.
1199    ///
1200    /// # Returns
1201    ///
1202    /// Returns a permutation that permutes the *ordinary* atoms of `self` to give `other`, or
1203    /// `None` if no such permutation exists.
1204    fn get_perm_of(&self, other: &Self) -> Option<Permutation<Self::Rank>> {
1205        let o_basis_atoms: HashMap<&BasisAtom, usize> = other
1206            .basis_atoms
1207            .iter()
1208            .enumerate()
1209            .map(|(i, basis_atom)| (basis_atom, i))
1210            .collect();
1211        let image_opt: Option<Vec<Self::Rank>> = self
1212            .basis_atoms
1213            .iter()
1214            .map(|s_basis_atom| o_basis_atoms.get(s_basis_atom).copied())
1215            .collect();
1216        image_opt.and_then(|image| Permutation::from_image(image).ok())
1217    }
1218
1219    fn permute(&self, perm: &Permutation<Self::Rank>) -> Result<Self, anyhow::Error> {
1220        let mut p_bao = self.clone();
1221        p_bao.permute_mut(perm)?;
1222        Ok(p_bao)
1223    }
1224
1225    fn permute_mut(&mut self, perm: &Permutation<Self::Rank>) -> Result<(), anyhow::Error> {
1226        permute_inplace(&mut self.basis_atoms, perm);
1227        Ok(())
1228    }
1229}
1230
1231impl<'a> fmt::Display for BasisAngularOrder<'a> {
1232    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
1233        let order_length = self
1234            .basis_shells()
1235            .map(|v| v.shell_order.to_string().chars().count())
1236            .max()
1237            .unwrap_or(20);
1238        let atom_index_length = self.n_atoms().to_string().chars().count();
1239        writeln!(f, "{}", "┈".repeat(17 + atom_index_length + order_length))?;
1240        writeln!(f, " {:>atom_index_length$}  Atom  Shell  Order", "#")?;
1241        writeln!(f, "{}", "┈".repeat(17 + atom_index_length + order_length))?;
1242        for (atm_i, batm) in self.basis_atoms.iter().enumerate() {
1243            let atm = batm.atom;
1244            for (shl_i, bshl) in batm.basis_shells.iter().enumerate() {
1245                if shl_i == 0 {
1246                    writeln!(
1247                        f,
1248                        " {:>atom_index_length$}  {:<4}  {:<5}  {:<order_length$}",
1249                        atm_i,
1250                        atm.atomic_symbol,
1251                        match &bshl.shell_order {
1252                            ShellOrder::Pure(_) | ShellOrder::Cart(_) => ANGMOM_LABELS
1253                                .get(usize::try_from(bshl.l).unwrap_or_else(|err| panic!("{err}")))
1254                                .copied()
1255                                .unwrap_or(&bshl.l.to_string())
1256                                .to_string(),
1257                            ShellOrder::Spinor(spinororder) => format!("{}/2", spinororder.two_j),
1258                        },
1259                        bshl.shell_order
1260                    )?;
1261                } else {
1262                    writeln!(
1263                        f,
1264                        " {:>atom_index_length$}  {:<4}  {:<5}  {:<order_length$}",
1265                        "",
1266                        "",
1267                        match &bshl.shell_order {
1268                            ShellOrder::Pure(_) | ShellOrder::Cart(_) => ANGMOM_LABELS
1269                                .get(usize::try_from(bshl.l).unwrap_or_else(|err| panic!("{err}")))
1270                                .copied()
1271                                .unwrap_or(&bshl.l.to_string())
1272                                .to_string(),
1273                            ShellOrder::Spinor(spinororder) => format!("{}/2", spinororder.two_j),
1274                        },
1275                        bshl.shell_order
1276                    )?;
1277                }
1278            }
1279        }
1280        writeln!(f, "{}", "┈".repeat(17 + atom_index_length + order_length))?;
1281        Ok(())
1282    }
1283}