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::{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;
28use crate::chartab::chartab_group::CharacterProperties;
29use crate::chartab::{DecompositionError, SubspaceDecomposable};
30use crate::group::GroupType;
31use crate::io::format::{log_subtitle, qsym2_output, QSym2Output};
32use crate::symmetry::symmetry_element::symmetry_operation::SpecialSymmetryTransformation;
33use crate::symmetry::symmetry_group::SymmetryGroupProperties;
34use crate::symmetry::symmetry_transformation::{SymmetryTransformable, SymmetryTransformationKind};
35use crate::target::determinant::SlaterDeterminant;
36
37impl<'a, T, SC> Overlap<T, Ix2> for SlaterDeterminant<'a, T, SC>
42where
43 T: Lapack
44 + ComplexFloat<Real = <T as Scalar>::Real>
45 + fmt::Debug
46 + Mul<<T as ComplexFloat>::Real, Output = T>,
47 <T as ComplexFloat>::Real: fmt::Debug
48 + approx::RelativeEq<<T as ComplexFloat>::Real>
49 + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
50 SC: StructureConstraint + Eq + fmt::Display,
51{
52 fn complex_symmetric(&self) -> bool {
53 self.complex_symmetric
54 }
55
56 fn overlap(
75 &self,
76 other: &Self,
77 metric: Option<&Array2<T>>,
78 metric_h: Option<&Array2<T>>,
79 ) -> Result<T, anyhow::Error> {
80 ensure!(
81 self.structure_constraint == other.structure_constraint,
82 "Inconsistent structure constraints between `self` and `other`."
83 );
84 ensure!(
85 self.coefficients.len() == other.coefficients.len(),
86 "Inconsistent numbers of coefficient matrices between `self` and `other`."
87 );
88 ensure!(
89 self.bao == other.bao,
90 "Inconsistent basis angular order between `self` and `other`."
91 );
92
93 let thresh = Float::sqrt(self.threshold * other.threshold);
94 ensure!(self
95 .occupations
96 .iter()
97 .chain(other.occupations.iter())
98 .all(|occs| occs.iter().all(|&occ| approx::relative_eq!(
99 occ,
100 occ.round(),
101 epsilon = thresh,
102 max_relative = thresh
103 ))),
104 "Overlaps between determinants with fractional occupation numbers are currently not supported."
105 );
106
107 let sao = metric.ok_or_else(|| format_err!("No atomic-orbital metric found."))?;
108 let sao_h = metric_h.unwrap_or(sao);
109
110 let ov = izip!(
111 &self.coefficients,
112 &self.occupations,
113 &other.coefficients,
114 &other.occupations
115 )
116 .map(|(cw, occw, cx, occx)| {
117 let nonzero_occ_w = occw.iter().positions(|&occ| occ > thresh).collect_vec();
118 let cw_o = cw.select(Axis(1), &nonzero_occ_w);
119 let nonzero_occ_x = occx.iter().positions(|&occ| occ > thresh).collect_vec();
120 let cx_o = cx.select(Axis(1), &nonzero_occ_x);
121
122 let mo_ov_mat = if self.complex_symmetric() {
123 match (self.complex_conjugated, other.complex_conjugated) {
124 (false, false) => cw_o.t().dot(sao_h).dot(&cx_o),
125 (true, false) => cw_o.t().dot(sao).dot(&cx_o),
126 (false, true) => cx_o.t().dot(sao).dot(&cw_o),
127 (true, true) => cw_o.t().dot(&sao_h.t()).dot(&cx_o),
128 }
129 } else {
130 match (self.complex_conjugated, other.complex_conjugated) {
131 (false, false) => cw_o.t().mapv(|x| x.conj()).dot(sao).dot(&cx_o),
132 (true, false) => cw_o.t().mapv(|x| x.conj()).dot(sao_h).dot(&cx_o),
133 (false, true) => cx_o
134 .t()
135 .mapv(|x| x.conj())
136 .dot(sao_h)
137 .dot(&cw_o)
138 .mapv(|x| x.conj()),
139 (true, true) => cw_o.t().mapv(|x| x.conj()).dot(&sao.t()).dot(&cx_o),
140 }
141 };
142 mo_ov_mat
143 .det()
144 .expect("The determinant of the MO overlap matrix could not be found.")
145 })
146 .fold(T::one(), |acc, x| acc * x);
147
148 let implicit_factor = self.structure_constraint.implicit_factor()?;
149 if implicit_factor > 1 {
150 let p_i32 = i32::try_from(implicit_factor)?;
151 Ok(ComplexFloat::powi(ov, p_i32))
152 } else {
153 Ok(ov)
154 }
155 }
156
157 fn overlap_definition(&self) -> String {
159 let k = if self.complex_symmetric() { "κ " } else { "" };
160 format!("⟨{k}Ψ_1|Ψ_2⟩ = ∫ [{k}Ψ_1(x^Ne)]* Ψ_2(x^Ne) dx^Ne")
161 }
162}
163
164#[derive(Builder, Clone)]
175pub struct SlaterDeterminantSymmetryOrbit<'a, G, T, SC>
176where
177 G: SymmetryGroupProperties,
178 T: ComplexFloat + fmt::Debug + Lapack,
179 SC: StructureConstraint + fmt::Display,
180 SlaterDeterminant<'a, T, SC>: SymmetryTransformable,
181{
182 group: &'a G,
184
185 origin: &'a SlaterDeterminant<'a, T, SC>,
187
188 linear_independence_threshold: <T as ComplexFloat>::Real,
190
191 integrality_threshold: <T as ComplexFloat>::Real,
194
195 symmetry_transformation_kind: SymmetryTransformationKind,
198
199 #[builder(setter(skip), default = "None")]
201 smat: Option<Array2<T>>,
202
203 #[builder(setter(skip), default = "None")]
206 pub(crate) smat_eigvals: Option<Array1<T>>,
207
208 #[builder(setter(skip), default = "None")]
213 xmat: Option<Array2<T>>,
214
215 eigenvalue_comparison_mode: EigenvalueComparisonMode,
218}
219
220impl<'a, G, T, SC> SlaterDeterminantSymmetryOrbit<'a, G, T, SC>
225where
226 G: SymmetryGroupProperties + Clone,
227 T: ComplexFloat + fmt::Debug + Lapack,
228 SC: StructureConstraint + Clone + fmt::Display,
229 SlaterDeterminant<'a, T, SC>: SymmetryTransformable,
230{
231 pub fn builder() -> SlaterDeterminantSymmetryOrbitBuilder<'a, G, T, SC> {
233 SlaterDeterminantSymmetryOrbitBuilder::default()
234 }
235}
236
237impl<'a, G, SC> SlaterDeterminantSymmetryOrbit<'a, G, f64, SC>
238where
239 G: SymmetryGroupProperties,
240 SC: StructureConstraint + fmt::Display,
241 SlaterDeterminant<'a, f64, SC>: SymmetryTransformable,
242{
243 fn_calc_xmat_real!(
244 pub calc_xmat
256 );
257}
258
259impl<'a, G, T, SC> SlaterDeterminantSymmetryOrbit<'a, G, Complex<T>, SC>
260where
261 G: SymmetryGroupProperties,
262 T: Float + Scalar<Complex = Complex<T>>,
263 Complex<T>: ComplexFloat<Real = T> + Scalar<Real = T, Complex = Complex<T>> + Lapack,
264 SC: StructureConstraint + fmt::Display,
265 SlaterDeterminant<'a, Complex<T>, SC>: SymmetryTransformable + Overlap<Complex<T>, Ix2>,
266{
267 fn_calc_xmat_complex!(
268 pub calc_xmat
280 );
281}
282
283impl<'a, G, T, SC> Orbit<G, SlaterDeterminant<'a, T, SC>>
292 for SlaterDeterminantSymmetryOrbit<'a, G, T, SC>
293where
294 G: SymmetryGroupProperties,
295 T: ComplexFloat + fmt::Debug + Lapack,
296 SC: StructureConstraint + fmt::Display,
297 SlaterDeterminant<'a, T, SC>: SymmetryTransformable,
298{
299 type OrbitIter = OrbitIterator<'a, G, SlaterDeterminant<'a, T, SC>>;
300
301 fn group(&self) -> &G {
302 self.group
303 }
304
305 fn origin(&self) -> &SlaterDeterminant<'a, T, SC> {
306 self.origin
307 }
308
309 fn iter(&self) -> Self::OrbitIter {
310 OrbitIterator::new(
311 self.group,
312 self.origin,
313 match self.symmetry_transformation_kind {
314 SymmetryTransformationKind::Spatial => |op, det| {
315 det.sym_transform_spatial(op).with_context(|| {
316 format!("Unable to apply `{op}` spatially on the origin determinant")
317 })
318 },
319 SymmetryTransformationKind::SpatialWithSpinTimeReversal => |op, det| {
320 det.sym_transform_spatial_with_spintimerev(op).with_context(|| {
321 format!("Unable to apply `{op}` spatially (with spin-including time reversal) on the origin determinant")
322 })
323 },
324 SymmetryTransformationKind::Spin => |op, det| {
325 det.sym_transform_spin(op).with_context(|| {
326 format!("Unable to apply `{op}` spin-wise on the origin determinant")
327 })
328 },
329 SymmetryTransformationKind::SpinSpatial => |op, det| {
330 det.sym_transform_spin_spatial(op).with_context(|| {
331 format!("Unable to apply `{op}` spin-spatially on the origin determinant")
332 })
333 },
334 },
335 )
336 }
337}
338
339impl<'a, G, T, SC> RepAnalysis<G, SlaterDeterminant<'a, T, SC>, T, Ix2>
344 for SlaterDeterminantSymmetryOrbit<'a, G, T, SC>
345where
346 G: SymmetryGroupProperties,
347 G::CharTab: SubspaceDecomposable<T>,
348 T: Lapack
349 + ComplexFloat<Real = <T as Scalar>::Real>
350 + fmt::Debug
351 + Mul<<T as ComplexFloat>::Real, Output = T>,
352 <T as ComplexFloat>::Real: fmt::Debug
353 + Zero
354 + From<u16>
355 + ToPrimitive
356 + approx::RelativeEq<<T as ComplexFloat>::Real>
357 + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
358 SC: StructureConstraint + Eq + fmt::Display,
359 SlaterDeterminant<'a, T, SC>: SymmetryTransformable,
360{
361 fn set_smat(&mut self, smat: Array2<T>) {
362 self.smat = Some(smat)
363 }
364
365 fn smat(&self) -> Option<&Array2<T>> {
366 self.smat.as_ref()
367 }
368
369 fn xmat(&self) -> &Array2<T> {
370 self.xmat
371 .as_ref()
372 .expect("Orbit overlap orthogonalisation matrix not found.")
373 }
374
375 fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error> {
376 if self.origin.complex_symmetric {
377 Err(format_err!("`norm_preserving_scalar_map` is currently not implemented for complex symmetric overlaps."))
378 } else {
379 if self
380 .group
381 .get_index(i)
382 .unwrap_or_else(|| panic!("Group operation index `{i}` not found."))
383 .contains_time_reversal()
384 {
385 Ok(ComplexFloat::conj)
386 } else {
387 Ok(|x| x)
388 }
389 }
390 }
391
392 fn integrality_threshold(&self) -> <T as ComplexFloat>::Real {
393 self.integrality_threshold
394 }
395
396 fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode {
397 &self.eigenvalue_comparison_mode
398 }
399
400 fn analyse_rep(
416 &self,
417 ) -> Result<
418 <<G as CharacterProperties>::CharTab as SubspaceDecomposable<T>>::Decomposition,
419 DecompositionError,
420 > {
421 log::debug!("Analysing representation symmetry for a Slater determinant...");
422 let nelectrons_float = self.origin().nelectrons();
423 if approx::relative_eq!(
424 nelectrons_float.round(),
425 nelectrons_float,
426 epsilon = self.integrality_threshold,
427 max_relative = self.integrality_threshold
428 ) {
429 let nelectrons_usize = nelectrons_float.round().to_usize().unwrap_or_else(|| {
430 panic!(
431 "Unable to convert the number of electrons `{nelectrons_float:.7}` to `usize`."
432 );
433 });
434 let (valid_symmetry, err_str) = if nelectrons_usize.rem_euclid(2) == 0 {
435 (true, String::new())
437 } else {
438 match self.symmetry_transformation_kind {
440 SymmetryTransformationKind::Spatial => (true, String::new()),
441 SymmetryTransformationKind::SpatialWithSpinTimeReversal
442 | SymmetryTransformationKind::Spin
443 | SymmetryTransformationKind::SpinSpatial => {
444 match self.group().group_type() {
445 GroupType::Ordinary(_) => (true, String::new()),
446 GroupType::MagneticGrey(_) | GroupType::MagneticBlackWhite(_) => {
447 (!self.group().unitary_represented(),
448 "Unitary-represented magnetic groups cannot be used for symmetry analysis of odd-electron systems where spin is treated explicitly.".to_string())
449 }
450 }
451 }
452 }
453 };
454
455 if valid_symmetry {
456 let chis = self
457 .calc_characters()
458 .map_err(|err| DecompositionError(err.to_string()))?;
459 log::debug!("Characters calculated.");
460
461 log_subtitle("Determinant orbit characters");
462 qsym2_output!("");
463 self.characters_to_string(&chis, self.integrality_threshold)
464 .log_output_display();
465 qsym2_output!("");
466
467 let res = self.group().character_table().reduce_characters(
468 &chis.iter().map(|(cc, chi)| (cc, *chi)).collect::<Vec<_>>(),
469 self.integrality_threshold(),
470 );
471 log::debug!("Characters reduced.");
472 log::debug!("Analysing representation symmetry for a Slater determinant... Done.");
473 res
474 } else {
475 Err(DecompositionError(err_str))
476 }
477 } else {
478 Err(DecompositionError(format!(
479 "Symmetry analysis for determinant with non-integer number of electrons `{nelectrons_float:.7}` (threshold = {:.3e}) not supported.",
480 self.integrality_threshold
481 )))
482 }
483 }
484}