1use std::fmt;
4use std::ops::Mul;
5
6use anyhow::{self, ensure, format_err, Context};
7use approx;
8use derive_builder::Builder;
9use itertools::{izip, Itertools};
10use log;
11use ndarray::{s, Array1, Array2, Axis, Ix2};
12use ndarray_linalg::{
13 eig::Eig,
14 eigh::Eigh,
15 solve::Determinant,
16 types::{Lapack, Scalar},
17 UPLO,
18};
19use num_complex::{Complex, ComplexFloat};
20use num_traits::{Float, ToPrimitive, Zero};
21
22use crate::analysis::{
23 fn_calc_xmat_complex, fn_calc_xmat_real, EigenvalueComparisonMode, Orbit, OrbitIterator,
24 Overlap, RepAnalysis,
25};
26use crate::angmom::spinor_rotation_3d::StructureConstraint;
27use crate::auxiliary::misc::{complex_modified_gram_schmidt, ProductRepeat};
28use crate::chartab::chartab_group::CharacterProperties;
29use crate::chartab::{DecompositionError, SubspaceDecomposable};
30use crate::group::GroupType;
31use crate::symmetry::symmetry_element::symmetry_operation::SpecialSymmetryTransformation;
32use crate::symmetry::symmetry_group::SymmetryGroupProperties;
33use crate::symmetry::symmetry_transformation::{SymmetryTransformable, SymmetryTransformationKind};
34use crate::target::determinant::determinant_analysis::SlaterDeterminantSymmetryOrbit;
35use crate::target::determinant::SlaterDeterminant;
36use crate::target::orbital::MolecularOrbital;
37
38impl<'a, T, SC> Overlap<T, Ix2> for MolecularOrbital<'a, T, SC>
43where
44 T: Lapack
45 + ComplexFloat<Real = <T as Scalar>::Real>
46 + fmt::Debug
47 + Mul<<T as ComplexFloat>::Real, Output = T>,
48 <T as ComplexFloat>::Real: fmt::Debug
49 + approx::RelativeEq<<T as ComplexFloat>::Real>
50 + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
51 SC: StructureConstraint + Eq + fmt::Display,
52{
53 fn complex_symmetric(&self) -> bool {
54 self.complex_symmetric
55 }
56
57 fn overlap(
66 &self,
67 other: &Self,
68 metric: Option<&Array2<T>>,
69 metric_h: Option<&Array2<T>>,
70 ) -> Result<T, anyhow::Error> {
71 ensure!(
72 self.structure_constraint == other.structure_constraint,
73 "Inconsistent spin constraints between `self` and `other`."
74 );
75 ensure!(
76 self.coefficients.len() == other.coefficients.len(),
77 "Inconsistent numbers of coefficient matrices between `self` and `other`."
78 );
79
80 let sao = metric.ok_or_else(|| format_err!("No atomic-orbital metric found."))?;
81 let sao_h = metric_h.unwrap_or(sao);
82 let ov = if self.component_index != other.component_index {
83 T::zero()
84 } else if self.complex_symmetric() {
85 match (self.complex_conjugated, other.complex_conjugated) {
86 (false, false) => self.coefficients.t().dot(sao_h).dot(&other.coefficients),
87 (true, false) => self.coefficients.t().dot(sao).dot(&other.coefficients),
88 (false, true) => other.coefficients.t().dot(sao).dot(&self.coefficients),
89 (true, true) => self
90 .coefficients
91 .t()
92 .dot(&sao_h.t())
93 .dot(&other.coefficients),
94 }
95 } else {
96 match (self.complex_conjugated, other.complex_conjugated) {
97 (false, false) => self
98 .coefficients
99 .t()
100 .mapv(|x| x.conj()) .dot(sao)
102 .dot(&other.coefficients),
103 (true, false) => self
104 .coefficients
105 .t()
106 .mapv(|x| x.conj()) .dot(sao_h)
108 .dot(&other.coefficients),
109 (false, true) => other
110 .coefficients
111 .t()
112 .mapv(|x| x.conj()) .dot(sao_h)
114 .dot(&self.coefficients)
115 .conj(),
116 (true, true) => self
117 .coefficients
118 .t()
119 .mapv(|x| x.conj()) .dot(&sao.t())
121 .dot(&other.coefficients),
122 }
123 };
124 Ok(ov)
125 }
126
127 fn overlap_definition(&self) -> String {
129 let k = if self.complex_symmetric() { "κ " } else { "" };
130 format!("⟨{k}ψ_1|ψ_2⟩ = ∫ [{k}ψ_1(x)]* ψ_2(x) dx")
131 }
132}
133
134#[derive(Builder, Clone)]
145pub struct MolecularOrbitalSymmetryOrbit<'a, G, T, SC>
146where
147 G: SymmetryGroupProperties,
148 T: ComplexFloat + fmt::Debug + Lapack,
149 SC: StructureConstraint + fmt::Display,
150 MolecularOrbital<'a, T, SC>: SymmetryTransformable,
151{
152 group: &'a G,
154
155 origin: &'a MolecularOrbital<'a, T, SC>,
157
158 integrality_threshold: <T as ComplexFloat>::Real,
161
162 pub(crate) linear_independence_threshold: <T as ComplexFloat>::Real,
164
165 symmetry_transformation_kind: SymmetryTransformationKind,
168
169 #[builder(setter(skip), default = "None")]
171 smat: Option<Array2<T>>,
172
173 #[builder(setter(skip), default = "None")]
176 pub(crate) smat_eigvals: Option<Array1<T>>,
177
178 #[builder(setter(skip), default = "None")]
183 xmat: Option<Array2<T>>,
184
185 pub(crate) eigenvalue_comparison_mode: EigenvalueComparisonMode,
188}
189
190impl<'a, G, T, SC> MolecularOrbitalSymmetryOrbit<'a, G, T, SC>
195where
196 G: SymmetryGroupProperties + Clone,
197 T: ComplexFloat + fmt::Debug + Lapack,
198 SC: StructureConstraint + Clone + fmt::Display,
199 MolecularOrbital<'a, T, SC>: SymmetryTransformable,
200{
201 pub fn builder() -> MolecularOrbitalSymmetryOrbitBuilder<'a, G, T, SC> {
203 MolecularOrbitalSymmetryOrbitBuilder::default()
204 }
205
206 pub fn from_orbitals(
221 group: &'a G,
222 orbitals: &'a [Vec<MolecularOrbital<'a, T, SC>>],
223 sym_kind: SymmetryTransformationKind,
224 eig_comp_mode: EigenvalueComparisonMode,
225 integrality_thresh: <T as ComplexFloat>::Real,
226 linear_independence_thresh: <T as ComplexFloat>::Real,
227 ) -> Vec<Vec<Self>> {
228 orbitals
229 .iter()
230 .map(|orbs_spin| {
231 orbs_spin
232 .iter()
233 .map(|orb| {
234 MolecularOrbitalSymmetryOrbit::builder()
235 .group(group)
236 .origin(orb)
237 .integrality_threshold(integrality_thresh)
238 .linear_independence_threshold(linear_independence_thresh)
239 .symmetry_transformation_kind(sym_kind.clone())
240 .eigenvalue_comparison_mode(eig_comp_mode.clone())
241 .build()
242 .expect("Unable to construct a molecular orbital symmetry orbit.")
243 })
244 .collect_vec()
245 })
246 .collect_vec()
247 }
248}
249
250impl<'a, G, SC> MolecularOrbitalSymmetryOrbit<'a, G, f64, SC>
251where
252 G: SymmetryGroupProperties,
253 SC: StructureConstraint + fmt::Display,
254 MolecularOrbital<'a, f64, SC>: SymmetryTransformable,
255{
256 fn_calc_xmat_real!(
257 pub calc_xmat
269 );
270}
271
272impl<'a, G, T, SC> MolecularOrbitalSymmetryOrbit<'a, G, Complex<T>, SC>
273where
274 G: SymmetryGroupProperties,
275 T: Float + Scalar<Complex = Complex<T>>,
276 Complex<T>: ComplexFloat<Real = T> + Scalar<Real = T, Complex = Complex<T>> + Lapack,
277 SC: StructureConstraint + fmt::Display,
278 MolecularOrbital<'a, Complex<T>, SC>: SymmetryTransformable + Overlap<Complex<T>, Ix2>,
279{
280 fn_calc_xmat_complex!(
281 pub calc_xmat
293 );
294}
295
296impl<'a, G, T, SC> Orbit<G, MolecularOrbital<'a, T, SC>>
305 for MolecularOrbitalSymmetryOrbit<'a, G, T, SC>
306where
307 G: SymmetryGroupProperties,
308 T: ComplexFloat + fmt::Debug + Lapack,
309 SC: StructureConstraint + fmt::Display,
310 MolecularOrbital<'a, T, SC>: SymmetryTransformable,
311{
312 type OrbitIter = OrbitIterator<'a, G, MolecularOrbital<'a, T, SC>>;
313
314 fn group(&self) -> &G {
315 self.group
316 }
317
318 fn origin(&self) -> &MolecularOrbital<'a, T, SC> {
319 self.origin
320 }
321
322 fn iter(&self) -> Self::OrbitIter {
323 OrbitIterator::new(
324 self.group,
325 self.origin,
326 match self.symmetry_transformation_kind {
327 SymmetryTransformationKind::Spatial => |op, orb| {
328 orb.sym_transform_spatial(op).with_context(|| {
329 format!("Unable to apply `{op}` spatially on the origin orbital")
330 })
331 },
332 SymmetryTransformationKind::SpatialWithSpinTimeReversal => |op, orb| {
333 orb.sym_transform_spatial_with_spintimerev(op).with_context(|| {
334 format!("Unable to apply `{op}` spatially (with spin-including time reversal) on the origin orbital")
335 })
336 },
337 SymmetryTransformationKind::Spin => |op, orb| {
338 orb.sym_transform_spin(op).with_context(|| {
339 format!("Unable to apply `{op}` spin-wise on the origin orbital")
340 })
341 },
342 SymmetryTransformationKind::SpinSpatial => |op, orb| {
343 orb.sym_transform_spin_spatial(op).with_context(|| {
344 format!("Unable to apply `{op}` spin-spatially on the origin orbital",)
345 })
346 },
347 },
348 )
349 }
350}
351
352impl<'a, G, T, SC> RepAnalysis<G, MolecularOrbital<'a, T, SC>, T, Ix2>
357 for MolecularOrbitalSymmetryOrbit<'a, G, T, SC>
358where
359 G: SymmetryGroupProperties,
360 G::CharTab: SubspaceDecomposable<T>,
361 T: Lapack
362 + ComplexFloat<Real = <T as Scalar>::Real>
363 + fmt::Debug
364 + Mul<<T as ComplexFloat>::Real, Output = T>,
365 <T as ComplexFloat>::Real: fmt::Debug
366 + Zero
367 + approx::RelativeEq<<T as ComplexFloat>::Real>
368 + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
369 SC: StructureConstraint + Eq + fmt::Display,
370 MolecularOrbital<'a, T, SC>: SymmetryTransformable,
371{
372 fn set_smat(&mut self, smat: Array2<T>) {
373 self.smat = Some(smat)
374 }
375
376 fn smat(&self) -> Option<&Array2<T>> {
377 self.smat.as_ref()
378 }
379
380 fn xmat(&self) -> &Array2<T> {
381 self.xmat
382 .as_ref()
383 .expect("Orbit overlap orthogonalisation matrix not found.")
384 }
385
386 fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error> {
387 if self.origin.complex_symmetric {
388 Err(format_err!("`norm_preserving_scalar_map` is currently not implemented for complex symmetric overlaps."))
389 } else {
390 if self
391 .group
392 .get_index(i)
393 .unwrap_or_else(|| panic!("Group operation index `{i}` not found."))
394 .contains_time_reversal()
395 {
396 Ok(ComplexFloat::conj)
397 } else {
398 Ok(|x| x)
399 }
400 }
401 }
402
403 fn integrality_threshold(&self) -> <T as ComplexFloat>::Real {
404 self.integrality_threshold
405 }
406
407 fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode {
408 &self.eigenvalue_comparison_mode
409 }
410
411 fn analyse_rep(
427 &self,
428 ) -> Result<
429 <<G as CharacterProperties>::CharTab as SubspaceDecomposable<T>>::Decomposition,
430 DecompositionError,
431 > {
432 let (valid_symmetry, err_str) = match self.symmetry_transformation_kind {
434 SymmetryTransformationKind::Spatial => (true, String::new()),
435 SymmetryTransformationKind::SpatialWithSpinTimeReversal
436 | SymmetryTransformationKind::Spin
437 | SymmetryTransformationKind::SpinSpatial => {
438 match self.group().group_type() {
439 GroupType::Ordinary(_) => (true, String::new()),
440 GroupType::MagneticGrey(_) | GroupType::MagneticBlackWhite(_) => {
441 (!self.group().unitary_represented(),
442 "Unitary-represented magnetic groups cannot be used for symmetry analysis of a one-electron molecular orbital where spin is treated explicitly.".to_string())
443 }
444 }
445 }
446 };
447 if valid_symmetry {
448 log::debug!("Analysing representation symmetry for an MO...");
449 let chis = self
450 .calc_characters()
451 .map_err(|err| DecompositionError(err.to_string()))?;
452 let res = self.group().character_table().reduce_characters(
453 &chis.iter().map(|(cc, chi)| (cc, *chi)).collect::<Vec<_>>(),
454 self.integrality_threshold(),
455 );
456 log::debug!("Analysing representation symmetry for an MO... Done.");
457 res
458 } else {
459 Err(DecompositionError(err_str))
460 }
461 }
462}
463
464pub fn generate_det_mo_orbits<'a, G, T, SC>(
497 det: &'a SlaterDeterminant<'a, T, SC>,
498 mos: &'a [Vec<MolecularOrbital<'a, T, SC>>],
499 group: &'a G,
500 metric: &Array2<T>,
501 metric_h: Option<&Array2<T>>,
502 integrality_threshold: <T as ComplexFloat>::Real,
503 linear_independence_threshold: <T as ComplexFloat>::Real,
504 symmetry_transformation_kind: SymmetryTransformationKind,
505 eigenvalue_comparison_mode: EigenvalueComparisonMode,
506 use_cayley_table: bool,
507) -> Result<
508 (
509 SlaterDeterminantSymmetryOrbit<'a, G, T, SC>,
510 Vec<Vec<MolecularOrbitalSymmetryOrbit<'a, G, T, SC>>>,
511 ),
512 anyhow::Error,
513>
514where
515 G: SymmetryGroupProperties + Clone,
516 G::CharTab: SubspaceDecomposable<T>,
517 T: Lapack
518 + ComplexFloat<Real = <T as Scalar>::Real>
519 + fmt::Debug
520 + Mul<<T as ComplexFloat>::Real, Output = T>,
521 <T as ComplexFloat>::Real: fmt::Debug
522 + Zero
523 + From<u16>
524 + ToPrimitive
525 + approx::RelativeEq<<T as ComplexFloat>::Real>
526 + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
527 SC: StructureConstraint + Clone + Eq + fmt::Display,
528 SlaterDeterminant<'a, T, SC>: SymmetryTransformable,
529 MolecularOrbital<'a, T, SC>: SymmetryTransformable,
530{
531 log::debug!("Constructing determinant and MO orbits in tandem...");
532 let order = group.order();
533 let mut det_orbit = SlaterDeterminantSymmetryOrbit::<G, T, SC>::builder()
534 .group(group)
535 .origin(det)
536 .integrality_threshold(integrality_threshold)
537 .linear_independence_threshold(linear_independence_threshold)
538 .symmetry_transformation_kind(symmetry_transformation_kind.clone())
539 .eigenvalue_comparison_mode(eigenvalue_comparison_mode.clone())
540 .build()
541 .map_err(|err| format_err!(err))?;
542 let mut mo_orbitss = MolecularOrbitalSymmetryOrbit::from_orbitals(
543 group,
544 mos,
545 symmetry_transformation_kind,
546 eigenvalue_comparison_mode,
547 integrality_threshold,
548 linear_independence_threshold,
549 );
550
551 let mut det_smat = Array2::<T>::zeros((order, order));
552 let mut mo_smatss = mo_orbitss
553 .iter()
554 .map(|mo_orbits| {
555 mo_orbits
556 .iter()
557 .map(|_| Array2::<T>::zeros((order, order)))
558 .collect::<Vec<_>>()
559 })
560 .collect::<Vec<_>>();
561
562 let thresh = det.threshold();
563
564 let sao = metric;
565 let sao_h = metric_h.unwrap_or(sao);
566
567 if let (Some(ctb), true) = (group.cayley_table(), use_cayley_table) {
568 log::debug!("Cayley table available. Group closure will be used to speed up overlap matrix computation.");
569 let mut det_smatw0 = Array1::<T>::zeros(order);
570 let mut mo_smatw0ss = mo_orbitss
571 .iter()
572 .map(|mo_orbits| {
573 mo_orbits
574 .iter()
575 .map(|_| Array1::<T>::zeros(order))
576 .collect::<Vec<_>>()
577 })
578 .collect::<Vec<_>>();
579 let det_0 = det_orbit.origin();
580 for (w, det_w_res) in det_orbit.iter().enumerate() {
581 let det_w = det_w_res?;
582 let w0_ov = izip!(
583 det_w.coefficients(),
584 det_w.occupations(),
585 det_0.coefficients(),
586 det_0.occupations(),
587 )
588 .enumerate()
589 .map(|(ispin, (cw, occw, c0, occ0))| {
590 let nonzero_occ_w = occw.iter().positions(|&occ| occ > thresh).collect_vec();
591 let nonzero_occ_0 = occ0.iter().positions(|&occ| occ > thresh).collect_vec();
592
593 let all_mo_w0_ov_mat = if det.complex_symmetric() {
594 match (det_w.complex_conjugated(), det_0.complex_conjugated()) {
595 (false, false) => cw.t().dot(sao_h).dot(c0),
596 (true, false) => cw.t().dot(sao).dot(c0),
597 (false, true) => c0.t().dot(sao).dot(cw),
598 (true, true) => cw.t().dot(&sao_h.t()).dot(c0),
599 }
600 } else {
601 match (det_w.complex_conjugated(), det_0.complex_conjugated()) {
602 (false, false) => cw.t().mapv(|x| x.conj()).dot(sao).dot(c0),
603 (true, false) => cw.t().mapv(|x| x.conj()).dot(sao_h).dot(c0),
604 (false, true) => c0
605 .t()
606 .mapv(|x| x.conj())
607 .dot(sao_h)
608 .dot(cw)
609 .mapv(|x| x.conj()),
610 (true, true) => cw.t().mapv(|x| x.conj()).dot(&sao.t()).dot(c0),
611 }
612 };
613
614 all_mo_w0_ov_mat
615 .diag()
616 .iter()
617 .enumerate()
618 .for_each(|(imo, mo_w0_ov)| {
619 mo_smatw0ss[ispin][imo][w] = *mo_w0_ov;
620 });
621
622 let occ_mo_w0_ov_mat = all_mo_w0_ov_mat
623 .select(Axis(0), &nonzero_occ_w)
624 .select(Axis(1), &nonzero_occ_0);
625
626 occ_mo_w0_ov_mat
627 .det()
628 .expect("The determinant of the MO overlap matrix `w0` could not be found.")
629 })
630 .fold(T::one(), |acc, x| acc * x);
631
632 let implicit_factor = det.structure_constraint().implicit_factor()?;
633 let w0_ov = if implicit_factor > 1 {
634 let p_i32 = i32::try_from(implicit_factor)?;
635 ComplexFloat::powi(w0_ov, p_i32)
636 } else {
637 w0_ov
638 };
639 det_smatw0[w] = w0_ov;
640 }
641
642 for (i, j) in (0..order).cartesian_product(0..order) {
643 let jinv = ctb
644 .slice(s![.., j])
645 .iter()
646 .position(|&x| x == 0)
647 .ok_or(format_err!(
648 "Unable to find the inverse of group element `{j}`."
649 ))?;
650 let jinv_i = ctb[(jinv, i)];
651
652 for (ispin, mo_smatw0s) in mo_smatw0ss.iter().enumerate() {
653 for (imo, mo_smat_w0) in mo_smatw0s.iter().enumerate() {
654 mo_smatss[ispin][imo][(i, j)] =
655 mo_orbitss[ispin][imo].norm_preserving_scalar_map(jinv)?(mo_smat_w0[jinv_i])
656 }
657 }
658
659 det_smat[(i, j)] = det_orbit.norm_preserving_scalar_map(jinv)?(det_smatw0[jinv_i]);
660 }
661 } else {
662 log::debug!("Cayley table not available or the use of Cayley table not requested. Overlap matrix will be constructed without group-closure speed-up.");
663 let indexed_dets = det_orbit
664 .iter()
665 .map(|det_res| det_res.map_err(|err| err.to_string()))
666 .enumerate()
667 .collect::<Vec<_>>();
668 for det_pair in indexed_dets.iter().product_repeat(2) {
669 let (w, det_w_res) = &det_pair[0];
670 let (x, det_x_res) = &det_pair[1];
671 let det_w = det_w_res
672 .as_ref()
673 .map_err(|err| format_err!(err.to_owned()))
674 .with_context(|| "One of the determinants in the orbit is not available")?;
675 let det_x = det_x_res
676 .as_ref()
677 .map_err(|err| format_err!(err.to_owned()))
678 .with_context(|| "One of the determinants in the orbit is not available")?;
679
680 let wx_ov = izip!(
681 det_w.coefficients(),
682 det_w.occupations(),
683 det_x.coefficients(),
684 det_x.occupations(),
685 )
686 .enumerate()
687 .map(|(ispin, (cw, occw, cx, occx))| {
688 let nonzero_occ_w = occw.iter().positions(|&occ| occ > thresh).collect_vec();
689 let nonzero_occ_x = occx.iter().positions(|&occ| occ > thresh).collect_vec();
690
691 let all_mo_wx_ov_mat = if det.complex_symmetric() {
697 match (det_w.complex_conjugated(), det_x.complex_conjugated()) {
698 (false, false) => cw.t().dot(sao_h).dot(cx),
699 (true, false) => cw.t().dot(sao).dot(cx),
700 (false, true) => cx.t().dot(sao).dot(cw),
701 (true, true) => cw.t().dot(&sao_h.t()).dot(cx),
702 }
703 } else {
704 match (det_w.complex_conjugated(), det_x.complex_conjugated()) {
705 (false, false) => cw.t().mapv(|x| x.conj()).dot(sao).dot(cx),
706 (true, false) => cw.t().mapv(|x| x.conj()).dot(sao_h).dot(cx),
707 (false, true) => cx
708 .t()
709 .mapv(|x| x.conj())
710 .dot(sao_h)
711 .dot(cw)
712 .mapv(|x| x.conj()),
713 (true, true) => cw.t().mapv(|x| x.conj()).dot(&sao.t()).dot(cx),
714 }
715 };
716
717 all_mo_wx_ov_mat
718 .diag()
719 .iter()
720 .enumerate()
721 .for_each(|(imo, mo_wx_ov)| {
722 mo_smatss[ispin][imo][(*w, *x)] = *mo_wx_ov;
723 });
724
725 let occ_mo_wx_ov_mat = all_mo_wx_ov_mat
726 .select(Axis(0), &nonzero_occ_w)
727 .select(Axis(1), &nonzero_occ_x);
728
729 occ_mo_wx_ov_mat
730 .det()
731 .expect("The determinant of the MO overlap matrix could not be found.")
732 })
733 .fold(T::one(), |acc, x| acc * x);
734
735 let implicit_factor = det.structure_constraint().implicit_factor()?;
736 let wx_ov = if implicit_factor > 1 {
737 let p_i32 = i32::try_from(implicit_factor)?;
738 ComplexFloat::powi(wx_ov, p_i32)
739 } else {
740 wx_ov
741 };
742 det_smat[(*w, *x)] = wx_ov;
743 }
744 }
745
746 if det_orbit.origin().complex_symmetric() {
747 det_orbit.set_smat(
748 (det_smat.clone() + det_smat.t().to_owned()).mapv(|x| x / (T::one() + T::one())),
749 );
750 } else {
751 det_orbit.set_smat(
752 (det_smat.clone() + det_smat.t().to_owned().mapv(|x| x.conj()))
753 .mapv(|x| x / (T::one() + T::one())),
754 );
755 };
756
757 mo_orbitss
758 .iter_mut()
759 .enumerate()
760 .for_each(|(ispin, mo_orbits)| {
761 mo_orbits
762 .iter_mut()
763 .enumerate()
764 .for_each(|(imo, mo_orbit)| {
765 if mo_orbit.origin().complex_symmetric() {
766 mo_orbit.set_smat(
767 (mo_smatss[ispin][imo].clone() + mo_smatss[ispin][imo].t().to_owned())
768 .mapv(|x| x / (T::one() + T::one())),
769 )
770 } else {
771 mo_orbit.set_smat(
772 (mo_smatss[ispin][imo].clone()
773 + mo_smatss[ispin][imo].t().to_owned().mapv(|x| x.conj()))
774 .mapv(|x| x / (T::one() + T::one())),
775 )
776 }
777 })
778 });
779
780 log::debug!("Constructing determinant and MO orbits in tandem... Done.");
781 Ok((det_orbit, mo_orbitss))
782}