1use std::collections::HashSet;
4use std::fmt;
5use std::hash::Hash;
6use std::ops::Mul;
7
8use anyhow::{self, ensure, format_err, Context};
9use approx;
10use derive_builder::Builder;
11use itertools::Itertools;
12use log;
13use ndarray::{s, Array1, Array2, Array3, Axis, Ix2};
14use ndarray_linalg::{
15 eig::Eig,
16 eigh::Eigh,
17 types::{Lapack, Scalar},
18 UPLO,
19};
20use num_complex::{Complex, ComplexFloat};
21use num_traits::{Float, ToPrimitive, Zero};
22
23use crate::analysis::{
24 fn_calc_xmat_complex, fn_calc_xmat_real, EigenvalueComparisonMode, Orbit, OrbitIterator,
25 Overlap, RepAnalysis,
26};
27use crate::angmom::spinor_rotation_3d::StructureConstraint;
28use crate::auxiliary::misc::complex_modified_gram_schmidt;
29use crate::chartab::chartab_group::CharacterProperties;
30use crate::chartab::{DecompositionError, SubspaceDecomposable};
31use crate::group::GroupType;
32use crate::io::format::{log_subtitle, qsym2_output, QSym2Output};
33use crate::symmetry::symmetry_element::symmetry_operation::SpecialSymmetryTransformation;
34use crate::symmetry::symmetry_group::SymmetryGroupProperties;
35use crate::symmetry::symmetry_transformation::{SymmetryTransformable, SymmetryTransformationKind};
36use crate::target::determinant::SlaterDeterminant;
37use crate::target::noci::basis::{Basis, OrbitBasis};
38use crate::target::noci::multideterminant::MultiDeterminant;
39
40impl<'a, T, B, SC> Overlap<T, Ix2> for MultiDeterminant<'a, T, B, SC>
45where
46 T: Lapack
47 + ComplexFloat<Real = <T as Scalar>::Real>
48 + fmt::Debug
49 + Mul<<T as ComplexFloat>::Real, Output = T>,
50 <T as ComplexFloat>::Real: fmt::Debug
51 + approx::RelativeEq<<T as ComplexFloat>::Real>
52 + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
53 SC: StructureConstraint + Hash + Eq + Clone + fmt::Display,
54 B: Basis<SlaterDeterminant<'a, T, SC>> + Clone,
55{
56 fn complex_symmetric(&self) -> bool {
57 self.complex_symmetric
58 }
59
60 fn overlap(
79 &self,
80 other: &Self,
81 metric: Option<&Array2<T>>,
82 metric_h: Option<&Array2<T>>,
83 ) -> Result<T, anyhow::Error> {
84 ensure!(
85 self.structure_constraint() == other.structure_constraint(),
86 "Inconsistent structure constraints between `self` and `other`."
87 );
88 ensure!(
89 self.coefficients.len() == other.coefficients.len(),
90 "Inconsistent numbers of coefficients between `self` and `other`."
91 );
92
93 let s_dets = self.basis.iter().collect::<Result<Vec<_>, _>>()?;
94 let o_dets = other.basis.iter().collect::<Result<Vec<_>, _>>()?;
95 let swx_vec = s_dets
96 .iter()
97 .cartesian_product(o_dets.iter())
98 .map(|(w, x)| w.overlap(x, metric, metric_h))
99 .collect::<Result<Vec<_>, _>>()?;
100 let d = self.coefficients.len();
101 let swx = Array2::from_shape_vec((d, d), swx_vec)?;
102
103 if self.complex_symmetric {
104 Ok(self.coefficients.t().dot(&swx).dot(&other.coefficients))
105 } else {
106 Ok(self
107 .coefficients
108 .t()
109 .mapv(|x| x.conj())
110 .dot(&swx)
111 .dot(&other.coefficients))
112 }
113 }
114
115 fn overlap_definition(&self) -> String {
118 let k = if self.complex_symmetric() { "κ " } else { "" };
119 format!("⟨{k}Ψ_1|Ψ_2⟩ = ∫ [{k}Ψ_1(x^Ne)]* Ψ_2(x^Ne) dx^Ne")
120 }
121}
122
123#[derive(Builder, Clone)]
134pub struct MultiDeterminantSymmetryOrbit<'a, 'g, G, T, B, SC>
135where
136 G: SymmetryGroupProperties,
137 T: ComplexFloat + fmt::Debug + Lapack,
138 SC: StructureConstraint + Hash + Eq + fmt::Display,
139 B: 'a + Basis<SlaterDeterminant<'a, T, SC>> + Clone,
140 MultiDeterminant<'a, T, B, SC>: SymmetryTransformable,
141{
142 group: &'g G,
144
145 origin: &'a MultiDeterminant<'a, T, B, SC>,
147
148 pub(crate) linear_independence_threshold: <T as ComplexFloat>::Real,
150
151 integrality_threshold: <T as ComplexFloat>::Real,
154
155 symmetry_transformation_kind: SymmetryTransformationKind,
158
159 #[builder(setter(skip), default = "None")]
162 smat: Option<Array2<T>>,
163
164 #[builder(setter(skip), default = "None")]
167 pub(crate) smat_eigvals: Option<Array1<T>>,
168
169 #[builder(setter(skip), default = "None")]
174 xmat: Option<Array2<T>>,
175
176 eigenvalue_comparison_mode: EigenvalueComparisonMode,
179}
180
181impl<'a, 'g, G, T, B, SC> MultiDeterminantSymmetryOrbit<'a, 'g, G, T, B, SC>
186where
187 G: SymmetryGroupProperties + Clone,
188 T: ComplexFloat + fmt::Debug + Lapack,
189 SC: StructureConstraint + Hash + Eq + Clone + fmt::Display,
190 B: 'a + Basis<SlaterDeterminant<'a, T, SC>> + Clone,
191 MultiDeterminant<'a, T, B, SC>: SymmetryTransformable,
192{
193 pub fn builder() -> MultiDeterminantSymmetryOrbitBuilder<'a, 'g, G, T, B, SC> {
195 MultiDeterminantSymmetryOrbitBuilder::default()
196 }
197
198 pub fn origin(&self) -> &MultiDeterminant<'a, T, B, SC> {
200 self.origin
201 }
202}
203
204impl<'a, 'g, G, B, SC> MultiDeterminantSymmetryOrbit<'a, 'g, G, f64, B, SC>
205where
206 G: SymmetryGroupProperties,
207 SC: StructureConstraint + Hash + Eq + fmt::Display,
208 B: 'a + Basis<SlaterDeterminant<'a, f64, SC>> + Clone,
209 MultiDeterminant<'a, f64, B, SC>: SymmetryTransformable,
210{
211 fn_calc_xmat_real!(
212 pub calc_xmat
224 );
225}
226
227impl<'a, 'g, G, T, B, SC> MultiDeterminantSymmetryOrbit<'a, 'g, G, Complex<T>, B, SC>
228where
229 G: SymmetryGroupProperties,
230 T: Float + Scalar<Complex = Complex<T>>,
231 Complex<T>: ComplexFloat<Real = T> + Scalar<Real = T, Complex = Complex<T>> + Lapack,
232 SC: StructureConstraint + Hash + Eq + fmt::Display,
233 B: 'a + Basis<SlaterDeterminant<'a, Complex<T>, SC>> + Clone,
234 MultiDeterminant<'a, Complex<T>, B, SC>: SymmetryTransformable + Overlap<Complex<T>, Ix2>,
235{
236 fn_calc_xmat_complex!(
237 pub calc_xmat
249 );
250}
251
252impl<'a, 'g, G, T, B, SC> Orbit<G, MultiDeterminant<'a, T, B, SC>>
261 for MultiDeterminantSymmetryOrbit<'a, 'g, G, T, B, SC>
262where
263 G: SymmetryGroupProperties,
264 T: ComplexFloat + fmt::Debug + Lapack,
265 SC: StructureConstraint + Hash + Eq + fmt::Display,
266 B: 'a + Basis<SlaterDeterminant<'a, T, SC>> + Clone,
267 MultiDeterminant<'a, T, B, SC>: SymmetryTransformable,
268{
269 type OrbitIter = OrbitIterator<'a, G, MultiDeterminant<'a, T, B, SC>>;
270
271 fn group(&self) -> &G {
272 self.group
273 }
274
275 fn origin(&self) -> &MultiDeterminant<'a, T, B, SC> {
276 self.origin
277 }
278
279 fn iter(&self) -> Self::OrbitIter {
280 OrbitIterator::new(
281 self.group,
282 self.origin,
283 match self.symmetry_transformation_kind {
284 SymmetryTransformationKind::Spatial => |op, multidet| {
285 multidet.sym_transform_spatial(op).with_context(|| {
286 format!("Unable to apply `{op}` spatially on the origin multi-determinantal wavefunction")
287 })
288 },
289 SymmetryTransformationKind::SpatialWithSpinTimeReversal => |op, multidet| {
290 multidet.sym_transform_spatial_with_spintimerev(op).with_context(|| {
291 format!("Unable to apply `{op}` spatially (with spin-including time reversal) on the origin multi-determinantal wavefunction")
292 })
293 },
294 SymmetryTransformationKind::Spin => |op, multidet| {
295 multidet.sym_transform_spin(op).with_context(|| {
296 format!("Unable to apply `{op}` spin-wise on the origin multi-determinantal wavefunction")
297 })
298 },
299 SymmetryTransformationKind::SpinSpatial => |op, multidet| {
300 multidet.sym_transform_spin_spatial(op).with_context(|| {
301 format!("Unable to apply `{op}` spin-spatially on the origin multi-determinantal wavefunction")
302 })
303 },
304 },
305 )
306 }
307}
308
309impl<'a, 'g, G, T, B, SC> RepAnalysis<G, MultiDeterminant<'a, T, B, SC>, T, Ix2>
314 for MultiDeterminantSymmetryOrbit<'a, 'g, G, T, B, SC>
315where
316 G: SymmetryGroupProperties,
317 G::CharTab: SubspaceDecomposable<T>,
318 T: Lapack
319 + ComplexFloat<Real = <T as Scalar>::Real>
320 + fmt::Debug
321 + Mul<<T as ComplexFloat>::Real, Output = T>,
322 <T as ComplexFloat>::Real: fmt::Debug
323 + Zero
324 + From<u16>
325 + ToPrimitive
326 + approx::RelativeEq<<T as ComplexFloat>::Real>
327 + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
328 SC: StructureConstraint + Hash + Eq + Clone + fmt::Display,
329 B: 'a + Basis<SlaterDeterminant<'a, T, SC>> + Clone,
330 MultiDeterminant<'a, T, B, SC>: SymmetryTransformable,
331{
332 fn set_smat(&mut self, smat: Array2<T>) {
333 self.smat = Some(smat)
334 }
335
336 fn smat(&self) -> Option<&Array2<T>> {
337 self.smat.as_ref()
338 }
339
340 fn xmat(&self) -> &Array2<T> {
341 self.xmat
342 .as_ref()
343 .expect("Orbit overlap orthogonalisation matrix not found.")
344 }
345
346 fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error> {
347 if self.origin.complex_symmetric {
348 Err(format_err!("`norm_preserving_scalar_map` is currently not implemented for complex-symmetric overlaps. This thus precludes the use of the Cayley table to speed up the computation of the orbit overlap matrix."))
349 } else {
350 if self
351 .group
352 .get_index(i)
353 .unwrap_or_else(|| panic!("Group operation index `{i}` not found."))
354 .contains_time_reversal()
355 {
356 Ok(ComplexFloat::conj)
357 } else {
358 Ok(|x| x)
359 }
360 }
361 }
362
363 fn integrality_threshold(&self) -> <T as ComplexFloat>::Real {
364 self.integrality_threshold
365 }
366
367 fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode {
368 &self.eigenvalue_comparison_mode
369 }
370
371 fn analyse_rep(
387 &self,
388 ) -> Result<
389 <<G as CharacterProperties>::CharTab as SubspaceDecomposable<T>>::Decomposition,
390 DecompositionError,
391 > {
392 log::debug!("Analysing representation symmetry for a multi-determinantal wavefunction...");
393 let mut nelectrons_set = self
394 .origin()
395 .basis
396 .iter()
397 .map(|det_res| det_res.and_then(|det| {
398 let det_nelectrons_float = det.nelectrons();
399 if approx::relative_eq!(
400 det_nelectrons_float.round(),
401 det_nelectrons_float,
402 epsilon = self.integrality_threshold,
403 max_relative = self.integrality_threshold
404 ) {
405 det_nelectrons_float.round().to_usize().ok_or(format_err!("Unable to convert the number of electrons `{det_nelectrons_float:.7}` to `usize`."))
406 } else {
407 Err(format_err!("Fractional number of electrons encountered: `{det_nelectrons_float:.7}`"))
408 }
409 }))
410 .collect::<Result<HashSet<_>, _>>()
411 .map_err(|err| DecompositionError(err.to_string()))?;
412 if nelectrons_set.len() != 1 {
413 Err(DecompositionError("Symmetry analysis for multi-determinantal wavefunctions with multiple numbers of electrons is not yet supported.".to_string()))
414 } else {
415 let nelectrons = nelectrons_set.drain().next().ok_or(DecompositionError(
416 "Unable to obtain the number of electrons.".to_string(),
417 ))?;
418 let (valid_symmetry, err_str) = if nelectrons.rem_euclid(2) == 0 {
419 (true, String::new())
421 } else {
422 match self.symmetry_transformation_kind {
424 SymmetryTransformationKind::Spatial => (true, String::new()),
425 SymmetryTransformationKind::SpatialWithSpinTimeReversal
426 | SymmetryTransformationKind::Spin
427 | SymmetryTransformationKind::SpinSpatial => {
428 match self.group().group_type() {
429 GroupType::Ordinary(_) => (true, String::new()),
430 GroupType::MagneticGrey(_) | GroupType::MagneticBlackWhite(_) => {
431 (!self.group().unitary_represented(),
432 "Unitary-represented magnetic groups cannot be used for symmetry analysis of odd-electron systems where spin is treated explicitly.".to_string())
433 }
434 }
435 }
436 }
437 };
438
439 if valid_symmetry {
440 let chis = self
441 .calc_characters()
442 .map_err(|err| DecompositionError(err.to_string()))?;
443 log::debug!("Characters calculated.");
444
445 log_subtitle("Multi-determinantal wavefunction orbit characters");
446 qsym2_output!("");
447 self.characters_to_string(&chis, self.integrality_threshold)
448 .log_output_display();
449 qsym2_output!("");
450
451 let res = self.group().character_table().reduce_characters(
452 &chis.iter().map(|(cc, chi)| (cc, *chi)).collect::<Vec<_>>(),
453 self.integrality_threshold(),
454 );
455 log::debug!("Characters reduced.");
456 log::debug!("Analysing representation symmetry for a multi-determinantal wavefunction... Done.");
457 res
458 } else {
459 Err(DecompositionError(err_str))
460 }
461 }
462 }
463}
464
465impl<'a, 'go, 'g, G, T, SC>
469 MultiDeterminantSymmetryOrbit<
470 'a,
471 'go,
472 G,
473 T,
474 OrbitBasis<'g, G, SlaterDeterminant<'a, T, SC>>,
475 SC,
476 >
477where
478 G: SymmetryGroupProperties + Clone,
479 G::CharTab: SubspaceDecomposable<T>,
480 T: Lapack
481 + ComplexFloat<Real = <T as Scalar>::Real>
482 + fmt::Debug
483 + Mul<<T as ComplexFloat>::Real, Output = T>,
484 <T as ComplexFloat>::Real: fmt::Debug
485 + Zero
486 + From<u16>
487 + ToPrimitive
488 + approx::RelativeEq<<T as ComplexFloat>::Real>
489 + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
490 SC: StructureConstraint + Hash + Eq + Clone + fmt::Display,
491 MultiDeterminant<'a, T, OrbitBasis<'g, G, SlaterDeterminant<'a, T, SC>>, SC>:
492 SymmetryTransformable,
493{
494 pub(crate) fn calc_smat_optimised(
511 &mut self,
512 metric: Option<&Array2<T>>,
513 metric_h: Option<&Array2<T>>,
514 use_cayley_table: bool,
515 ) -> Result<&mut Self, anyhow::Error> {
516 ensure!(
517 self.group.name() == self.origin.basis.group().name(),
518 "Multi-determinantal wavefunction orbit-generating group does not match symmetry analysis group."
519 );
520
521 if let (Some(ctb), true) = (self.group().cayley_table(), use_cayley_table) {
522 log::debug!("Cayley table available. Group closure will be used to speed up overlap matrix computation.");
523
524 let order = self.group.order();
525 let mut smat = Array2::<T>::zeros((order, order));
526 let multidet_0 = self.origin();
527 let det_origins = multidet_0.basis.origins();
528 let n_det_origins = det_origins.len();
529
530 let detov_kpwx_vec = multidet_0
531 .basis
532 .iter()
533 .collect::<Result<Vec<_>, _>>()?
534 .into_iter()
535 .cartesian_product(det_origins.iter())
536 .map(|(w, x)| w.overlap(x, metric, metric_h))
537 .collect::<Result<Vec<_>, _>>()?;
538 let detov_kpwx =
539 Array3::from_shape_vec((order, n_det_origins, n_det_origins), detov_kpwx_vec)?;
540 let ovs = (0..order)
541 .map(|k| {
542 [
543 (0..order),
544 (0..order),
545 (0..n_det_origins),
546 (0..n_det_origins),
547 ]
548 .into_iter()
549 .multi_cartesian_product()
550 .try_fold(T::zero(), |acc, v| {
551 let ip = v[0];
552 let jp = v[1];
553 let w = v[2];
554 let x = v[3];
555 let aipw = multidet_0.coefficients[ip * n_det_origins + w];
556 let ajpx = multidet_0.coefficients[jp * n_det_origins + x];
557
558 let jpinv = ctb.slice(s![.., jp]).iter().position(|&x| x == 0).ok_or(
559 format_err!("Unable to find the inverse of group element `{jp}`."),
560 )?;
561 let kp = ctb[(jpinv, ctb[(k, ip)])];
562 let ukipwjpx = self.norm_preserving_scalar_map(jp)?(detov_kpwx[(kp, w, x)]);
563
564 Ok::<_, anyhow::Error>(
565 acc + self.norm_preserving_scalar_map(k)?(aipw.conj())
566 * ukipwjpx
567 * ajpx,
568 )
569 })
570 })
571 .collect::<Result<Vec<_>, _>>()?;
572 for (i, j) in (0..order).cartesian_product(0..order) {
573 let jinv = ctb
574 .slice(s![.., j])
575 .iter()
576 .position(|&x| x == 0)
577 .ok_or(format_err!(
578 "Unable to find the inverse of group element `{j}`."
579 ))?;
580 let jinv_i = ctb[(jinv, i)];
581 smat[(i, j)] = self.norm_preserving_scalar_map(jinv)?(ovs[jinv_i]);
582 }
583 if self.origin().complex_symmetric() {
584 self.set_smat(
585 (smat.clone() + smat.t().to_owned()).mapv(|x| x / (T::one() + T::one())),
586 )
587 } else {
588 self.set_smat(
589 (smat.clone() + smat.t().to_owned().mapv(|x| x.conj()))
590 .mapv(|x| x / (T::one() + T::one())),
591 )
592 }
593 Ok(self)
594 } else {
595 self.calc_smat(metric, metric_h, use_cayley_table)
596 }
597 }
598}