1use 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#[derive(Clone, Builder, PartialEq, Eq, Hash)]
33pub struct PureOrder {
34 #[builder(setter(custom))]
36 mls: Vec<i32>,
37
38 pub lpure: u32,
40
41 #[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 fn builder() -> PureOrderBuilder {
73 PureOrderBuilder::default()
74 }
75
76 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 #[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 #[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 #[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 #[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 pub fn iter(&self) -> Iter<i32> {
178 self.mls.iter()
179 }
180
181 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 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#[derive(Clone, Builder, PartialEq, Eq, Hash)]
253pub struct CartOrder {
254 #[builder(setter(custom))]
256 pub cart_tuples: Vec<(u32, u32, u32)>,
257
258 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 fn builder() -> CartOrderBuilder {
283 CartOrderBuilder::default()
284 }
285
286 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 #[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 #[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 #[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 #[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 pub fn iter(&self) -> Iter<(u32, u32, u32)> {
457 self.cart_tuples.iter()
458 }
459
460 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 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
529pub(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#[derive(Clone, Builder, PartialEq, Eq, Hash)]
571pub struct SpinorOrder {
572 #[builder(setter(custom))]
574 pub two_j: u32,
575
576 #[builder(setter(custom))]
578 two_mjs: Vec<i32>,
579
580 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 fn builder() -> SpinorOrderBuilder {
619 SpinorOrderBuilder::default()
620 }
621
622 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 #[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 #[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 #[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 #[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 pub fn iter(&self) -> Iter<i32> {
731 self.two_mjs.iter()
732 }
733
734 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 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#[derive(Clone, PartialEq, Eq, Hash, Debug)]
820pub enum ShellOrder {
821 Pure(PureOrder),
824
825 Cart(CartOrder),
828
829 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#[derive(Clone, Builder, PartialEq, Eq, Hash, Debug)]
871pub struct BasisShell {
872 #[builder(setter(custom))]
877 pub l: u32,
878
879 #[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 fn builder() -> BasisShellBuilder {
915 BasisShellBuilder::default()
916 }
917
918 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 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#[derive(Clone, Builder, PartialEq, Eq, Hash, Debug)]
956pub struct BasisAtom<'a> {
957 pub(crate) atom: &'a Atom,
959
960 #[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 pub(crate) fn builder() -> BasisAtomBuilder<'a> {
979 BasisAtomBuilder::default()
980 }
981
982 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 fn n_funcs(&self) -> usize {
999 self.basis_shells.iter().map(BasisShell::n_funcs).sum()
1000 }
1001
1002 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#[derive(Clone, Builder, PartialEq, Eq, Hash, Debug)]
1023pub struct BasisAngularOrder<'a> {
1024 #[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 #[must_use]
1043 pub(crate) fn builder() -> BasisAngularOrderBuilder<'a> {
1044 BasisAngularOrderBuilder::default()
1045 }
1046
1047 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 pub fn n_atoms(&self) -> usize {
1061 self.basis_atoms.len()
1062 }
1063
1064 pub fn n_funcs(&self) -> usize {
1066 self.basis_atoms.iter().map(BasisAtom::n_funcs).sum()
1067 }
1068
1069 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 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 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 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 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}