1use std::fmt;
4use std::ops::Mul;
5
6use anyhow::{self, ensure, format_err, Context};
7use approx;
8use derive_builder::Builder;
9use itertools::Itertools;
10use log;
11use ndarray::{Array1, Array2, Array4, Axis, Ix4};
12use ndarray_einsum_beta::*;
13use ndarray_linalg::{
14 eig::Eig,
15 eigh::Eigh,
16 norm::Norm,
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::auxiliary::misc::complex_modified_gram_schmidt;
28use crate::chartab::chartab_group::CharacterProperties;
29use crate::chartab::{DecompositionError, SubspaceDecomposable};
30use crate::symmetry::symmetry_element::symmetry_operation::SpecialSymmetryTransformation;
31use crate::symmetry::symmetry_group::SymmetryGroupProperties;
32use crate::symmetry::symmetry_transformation::{SymmetryTransformable, SymmetryTransformationKind};
33use crate::target::density::Density;
34
35impl<'a, T> Overlap<T, Ix4> for Density<'a, T>
40where
41 T: Lapack
42 + ComplexFloat<Real = <T as Scalar>::Real>
43 + fmt::Debug
44 + Mul<<T as ComplexFloat>::Real, Output = T>,
45 <T as ComplexFloat>::Real: fmt::Debug
46 + approx::RelativeEq<<T as ComplexFloat>::Real>
47 + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
48{
49 fn complex_symmetric(&self) -> bool {
50 self.complex_symmetric
51 }
52
53 fn overlap(
63 &self,
64 other: &Self,
65 metric: Option<&Array4<T>>,
66 metric_h: Option<&Array4<T>>,
67 ) -> Result<T, anyhow::Error> {
68 ensure!(
69 self.density_matrix.shape() == other.density_matrix.shape(),
70 "Inconsistent shapes of density matrices between `self` and `other`."
71 );
72 ensure!(
73 self.bao == other.bao,
74 "Inconsistent basis angular order between `self` and `other`."
75 );
76 let dennorm = self.density_matrix().norm_l2();
77 ensure!(
78 approx::abs_diff_ne!(
79 dennorm,
80 <T as ndarray_linalg::Scalar>::Real::zero(),
81 epsilon = self.threshold()
82 ),
83 "Zero density (density matrix has Frobenius norm {dennorm:.3e})"
84 );
85
86 let sao_4c = metric
87 .ok_or_else(|| format_err!("No atomic-orbital four-centre overlap tensor found for density overlap calculation."))?;
88 let sao_4c_h = metric_h.unwrap_or(sao_4c);
89
90 if self.complex_symmetric() {
91 match (self.complex_conjugated, other.complex_conjugated) {
92 (false, false) => einsum(
93 "ijkl,ji,lk->",
94 &[
95 &sao_4c_h.view(),
96 &self.density_matrix().view(),
97 &other.density_matrix().view(),
98 ],
99 )
100 .map_err(|err| format_err!(err))?
101 .into_iter()
102 .next()
103 .ok_or(format_err!("Unable to extract the density overlap scalar.")),
104 (true, false) => einsum(
105 "ijkl,ji,lk->",
106 &[
107 &sao_4c.view(),
108 &self.density_matrix().view(),
109 &other.density_matrix().view(),
110 ],
111 )
112 .map_err(|err| format_err!(err))?
113 .into_iter()
114 .next()
115 .ok_or(format_err!("Unable to extract the density overlap scalar.")),
116 (false, true) => einsum(
117 "ijkl,ji,lk->",
118 &[
119 &sao_4c.view(),
120 &other.density_matrix().view(),
121 &self.density_matrix().view(),
122 ],
123 )
124 .map_err(|err| format_err!(err))?
125 .into_iter()
126 .next()
127 .ok_or(format_err!("Unable to extract the density overlap scalar.")),
128 (true, true) => einsum(
129 "klij,ji,lk->",
130 &[
131 &sao_4c_h.view(),
132 &self.density_matrix().view(),
133 &other.density_matrix().view(),
134 ],
135 )
136 .map_err(|err| format_err!(err))?
137 .into_iter()
138 .next()
139 .ok_or(format_err!("Unable to extract the density overlap scalar.")),
140 }
141 } else {
142 match (self.complex_conjugated, other.complex_conjugated) {
143 (false, false) => einsum(
144 "ijkl,ji,lk->",
145 &[
146 &sao_4c.view(),
147 &self.density_matrix().mapv(|x| x.conj()).view(),
148 &other.density_matrix().view(),
149 ],
150 )
151 .map_err(|err| format_err!(err))?
152 .into_iter()
153 .next()
154 .ok_or(format_err!("Unable to extract the density overlap scalar.")),
155 (true, false) => einsum(
156 "ijkl,ji,lk->",
157 &[
158 &sao_4c_h.view(),
159 &self.density_matrix().mapv(|x| x.conj()).view(),
160 &other.density_matrix().view(),
161 ],
162 )
163 .map_err(|err| format_err!(err))?
164 .into_iter()
165 .next()
166 .ok_or(format_err!("Unable to extract the density overlap scalar.")),
167 (false, true) => einsum(
168 "ijkl,ji,lk->",
169 &[
170 &sao_4c_h.view(),
171 &other.density_matrix().mapv(|x| x.conj()).view(),
172 &self.density_matrix().view(),
173 ],
174 )
175 .map_err(|err| format_err!(err))?
176 .into_iter()
177 .next()
178 .ok_or(format_err!("Unable to extract the density overlap scalar."))
179 .map(|x| x.conj()),
180 (true, true) => einsum(
181 "klij,ji,lk->",
182 &[
183 &sao_4c.view(),
184 &self.density_matrix().mapv(|x| x.conj()).view(),
185 &other.density_matrix().view(),
186 ],
187 )
188 .map_err(|err| format_err!(err))?
189 .into_iter()
190 .next()
191 .ok_or(format_err!("Unable to extract the density overlap scalar.")),
192 }
193 }
194 }
195
196 fn overlap_definition(&self) -> String {
198 let k = if self.complex_symmetric() { "κ " } else { "" };
199 format!("⟨{k}ρ_1|ρ_2⟩ = ∫ [{k}ρ_1(r)]* ρ_2(r) dr")
200 }
201}
202
203#[derive(Builder, Clone)]
214pub struct DensitySymmetryOrbit<'a, G, T>
215where
216 G: SymmetryGroupProperties,
217 T: ComplexFloat + fmt::Debug + Lapack,
218 Density<'a, T>: SymmetryTransformable,
219{
220 group: &'a G,
222
223 origin: &'a Density<'a, T>,
225
226 pub(crate) linear_independence_threshold: <T as ComplexFloat>::Real,
228
229 integrality_threshold: <T as ComplexFloat>::Real,
232
233 symmetry_transformation_kind: SymmetryTransformationKind,
236
237 #[builder(setter(skip), default = "None")]
239 smat: Option<Array2<T>>,
240
241 #[builder(setter(skip), default = "None")]
244 pub(crate) smat_eigvals: Option<Array1<T>>,
245
246 #[builder(setter(skip), default = "None")]
251 xmat: Option<Array2<T>>,
252
253 eigenvalue_comparison_mode: EigenvalueComparisonMode,
256}
257
258impl<'a, G, T> DensitySymmetryOrbit<'a, G, T>
263where
264 G: SymmetryGroupProperties + Clone,
265 T: ComplexFloat + fmt::Debug + Lapack,
266 Density<'a, T>: SymmetryTransformable,
267{
268 pub fn builder() -> DensitySymmetryOrbitBuilder<'a, G, T> {
270 DensitySymmetryOrbitBuilder::default()
271 }
272}
273
274impl<'a, G> DensitySymmetryOrbit<'a, G, f64>
275where
276 G: SymmetryGroupProperties,
277{
278 fn_calc_xmat_real!(
279 pub calc_xmat
291 );
292}
293
294impl<'a, G, T> DensitySymmetryOrbit<'a, G, Complex<T>>
295where
296 G: SymmetryGroupProperties,
297 T: Float + Scalar<Complex = Complex<T>>,
298 Complex<T>: ComplexFloat<Real = T> + Scalar<Real = T, Complex = Complex<T>> + Lapack,
299 Density<'a, Complex<T>>: SymmetryTransformable + Overlap<Complex<T>, Ix4>,
300{
301 fn_calc_xmat_complex!(
302 pub calc_xmat
314 );
315}
316
317impl<'a, G, T> Orbit<G, Density<'a, T>> for DensitySymmetryOrbit<'a, G, T>
326where
327 G: SymmetryGroupProperties,
328 T: ComplexFloat + fmt::Debug + Lapack,
329 Density<'a, T>: SymmetryTransformable,
330{
331 type OrbitIter = OrbitIterator<'a, G, Density<'a, T>>;
332
333 fn group(&self) -> &G {
334 self.group
335 }
336
337 fn origin(&self) -> &Density<'a, T> {
338 self.origin
339 }
340
341 fn iter(&self) -> Self::OrbitIter {
342 OrbitIterator::new(
343 self.group,
344 self.origin,
345 match self.symmetry_transformation_kind {
346 SymmetryTransformationKind::Spatial => |op, det| {
347 det.sym_transform_spatial(op).with_context(|| {
348 format!("Unable to apply `{op}` spatially on the origin density")
349 })
350 },
351 SymmetryTransformationKind::SpatialWithSpinTimeReversal => |op, det| {
352 det.sym_transform_spatial_with_spintimerev(op).with_context(|| {
353 format!("Unable to apply `{op}` spatially (with spin-including time-reversal) on the origin density")
354 })
355 },
356 SymmetryTransformationKind::Spin => |op, det| {
357 det.sym_transform_spin(op).with_context(|| {
358 format!("Unable to apply `{op}` spin-wise on the origin density")
359 })
360 },
361 SymmetryTransformationKind::SpinSpatial => |op, det| {
362 det.sym_transform_spin_spatial(op).with_context(|| {
363 format!("Unable to apply `{op}` spin-spatially on the origin density")
364 })
365 },
366 },
367 )
368 }
369}
370
371impl<'a, G, T> RepAnalysis<G, Density<'a, T>, T, Ix4> for DensitySymmetryOrbit<'a, G, T>
376where
377 G: SymmetryGroupProperties,
378 G::CharTab: SubspaceDecomposable<T>,
379 T: Lapack
380 + ComplexFloat<Real = <T as Scalar>::Real>
381 + fmt::Debug
382 + Mul<<T as ComplexFloat>::Real, Output = T>,
383 <T as ComplexFloat>::Real: fmt::Debug
384 + Zero
385 + From<u16>
386 + ToPrimitive
387 + approx::RelativeEq<<T as ComplexFloat>::Real>
388 + approx::AbsDiffEq<Epsilon = <T as Scalar>::Real>,
389 Density<'a, T>: SymmetryTransformable,
390{
391 fn set_smat(&mut self, smat: Array2<T>) {
392 self.smat = Some(smat)
393 }
394
395 fn smat(&self) -> Option<&Array2<T>> {
396 self.smat.as_ref()
397 }
398
399 fn xmat(&self) -> &Array2<T> {
400 self.xmat
401 .as_ref()
402 .expect("Orbit overlap orthogonalisation matrix not found.")
403 }
404
405 fn norm_preserving_scalar_map(&self, i: usize) -> Result<fn(T) -> T, anyhow::Error> {
406 if self.origin.complex_symmetric {
407 Err(format_err!("`norm_preserving_scalar_map` is currently not implemented for complex symmetric overlaps."))
408 } else {
409 if self
410 .group
411 .get_index(i)
412 .unwrap_or_else(|| panic!("Group operation index `{i}` not found."))
413 .contains_time_reversal()
414 {
415 Ok(ComplexFloat::conj)
416 } else {
417 Ok(|x| x)
418 }
419 }
420 }
421
422 fn integrality_threshold(&self) -> <T as ComplexFloat>::Real {
423 self.integrality_threshold
424 }
425
426 fn eigenvalue_comparison_mode(&self) -> &EigenvalueComparisonMode {
427 &self.eigenvalue_comparison_mode
428 }
429
430 fn analyse_rep(
443 &self,
444 ) -> Result<
445 <<G as CharacterProperties>::CharTab as SubspaceDecomposable<T>>::Decomposition,
446 DecompositionError,
447 > {
448 log::debug!("Analysing representation symmetry for a density...");
449 let chis = self
450 .calc_characters()
451 .map_err(|err| DecompositionError(err.to_string()))?;
452 log::debug!("Characters calculated.");
453 let res = self.group().character_table().reduce_characters(
454 &chis.iter().map(|(cc, chi)| (cc, *chi)).collect::<Vec<_>>(),
455 self.integrality_threshold(),
456 );
457 log::debug!("Characters reduced.");
458 log::debug!("Analysing representation symmetry for a density... Done.");
459 res
460 }
461}