1use std::fmt;
4use std::iter::Sum;
5
6use anyhow::{self, format_err};
7use approx;
8use derive_builder::Builder;
9use itertools::Itertools;
10use log;
11use ndarray::{s, Array1, Array2, Ix2};
12use ndarray_einsum_beta::*;
13use ndarray_linalg::types::Lapack;
14use num::ToPrimitive;
15use num_complex::{Complex, ComplexFloat};
16use num_traits::float::{Float, FloatConst};
17
18use crate::angmom::spinor_rotation_3d::{SpinConstraint, SpinOrbitCoupled, StructureConstraint};
19use crate::auxiliary::molecule::Molecule;
20use crate::basis::ao::BasisAngularOrder;
21use crate::target::density::{DensitiesOwned, Density};
22use crate::target::orbital::MolecularOrbital;
23
24#[cfg(test)]
25mod determinant_tests;
26
27pub mod determinant_analysis;
28mod determinant_transformation;
29
30#[derive(Builder, Clone)]
36#[builder(build_fn(validate = "Self::validate"))]
37pub struct SlaterDeterminant<'a, T, SC>
38where
39 T: ComplexFloat + Lapack,
40 SC: StructureConstraint + fmt::Display,
41{
42 structure_constraint: SC,
44
45 bao: &'a BasisAngularOrder<'a>,
48
49 complex_symmetric: bool,
52
53 #[builder(default = "false")]
56 complex_conjugated: bool,
57
58 mol: &'a Molecule,
60
61 #[builder(setter(custom))]
63 coefficients: Vec<Array2<T>>,
64
65 #[builder(setter(custom))]
67 occupations: Vec<Array1<<T as ComplexFloat>::Real>>,
68
69 #[builder(default = "None")]
71 mo_energies: Option<Vec<Array1<T>>>,
72
73 #[builder(default = "Err(\"Determinant energy not yet set.\".to_string())")]
75 energy: Result<T, String>,
76
77 threshold: <T as ComplexFloat>::Real,
79}
80
81impl<'a, T, SC> SlaterDeterminantBuilder<'a, T, SC>
82where
83 T: ComplexFloat + Lapack,
84 SC: StructureConstraint + fmt::Display,
85{
86 pub fn coefficients(&mut self, cs: &[Array2<T>]) -> &mut Self {
87 self.coefficients = Some(cs.to_vec());
88 self
89 }
90
91 pub fn occupations(&mut self, occs: &[Array1<<T as ComplexFloat>::Real>]) -> &mut Self {
92 self.occupations = Some(occs.to_vec());
93 self
94 }
95
96 fn validate(&self) -> Result<(), String> {
97 let bao = self
98 .bao
99 .ok_or("No `BasisAngularOrder` found.".to_string())?;
100 let nbas = bao.n_funcs();
101 let coefficients = self
102 .coefficients
103 .as_ref()
104 .ok_or("No coefficients found.".to_string())?;
105 let structcons = self
106 .structure_constraint
107 .as_ref()
108 .ok_or("No structure constraints found.".to_string())?;
109 let coefficients_length_check = if coefficients.len() != structcons.n_coefficient_matrices()
110 {
111 log::error!(
112 "Unexpected number of coefficient matrices: {} found, but {} expected for the structure constraint {}.",
113 coefficients.len(),
114 structcons.n_coefficient_matrices(),
115 structcons
116 );
117 false
118 } else {
119 true
120 };
121 let coefficients_shape_check = {
122 let nrows = nbas * structcons.n_explicit_comps_per_coefficient_matrix();
123 if !coefficients.iter().all(|c| c.shape()[0] == nrows) {
124 log::error!(
125 "Unexpected shapes of coefficient matrices: {} {} expected for all coefficient matrices, but {} found.",
126 nrows,
127 if nrows == 1 {"row"} else {"rows"},
128 coefficients.iter().map(|c| c.shape()[0].to_string()).join(", ")
129 );
130 false
131 } else {
132 true
133 }
134 };
135
136 let occupations = self
137 .occupations
138 .as_ref()
139 .ok_or("No occupations found.".to_string())?;
140 let occupations_length_check = if occupations.len() != structcons.n_coefficient_matrices() {
141 log::error!(
142 "Unexpected number of occupation vectors: {} found, but {} expected for the structure constraint {}.",
143 occupations.len(),
144 structcons.n_coefficient_matrices(),
145 structcons
146 );
147 false
148 } else {
149 true
150 };
151 let occupations_shape_check = if !occupations
152 .iter()
153 .zip(coefficients.iter())
154 .all(|(o, c)| o.len() == c.shape()[1])
155 {
156 log::error!(
157 "Mismatched occupations and numbers of orbitals: {}",
158 occupations
159 .iter()
160 .zip(coefficients.iter())
161 .map(|(o, c)| format!("{} vs. {}", o.len(), c.shape()[1]))
162 .join(", ")
163 );
164 false
165 } else {
166 true
167 };
168
169 let mol = self.mol.ok_or("No molecule found.".to_string())?;
170 let natoms_check = mol.atoms.len() == bao.n_atoms();
171 if !natoms_check {
172 log::error!("The number of atoms in the molecule does not match the number of local sites in the basis.");
173 }
174
175 if coefficients_length_check
176 && coefficients_shape_check
177 && occupations_length_check
178 && occupations_shape_check
179 && natoms_check
180 {
181 Ok(())
182 } else {
183 Err("Slater determinant validation failed.".to_string())
184 }
185 }
186}
187
188impl<'a, T, SC> SlaterDeterminant<'a, T, SC>
189where
190 T: ComplexFloat + Clone + Lapack,
191 SC: StructureConstraint + Clone + fmt::Display,
192{
193 pub fn builder() -> SlaterDeterminantBuilder<'a, T, SC> {
195 SlaterDeterminantBuilder::default()
196 }
197}
198
199impl<'a, T, SC> SlaterDeterminant<'a, T, SC>
200where
201 T: ComplexFloat + Clone + Lapack,
202 SC: StructureConstraint + Clone + fmt::Display,
203{
204 pub fn to_orbitals(&self) -> Vec<Vec<MolecularOrbital<'a, T, SC>>> {
212 self.coefficients
213 .iter()
214 .enumerate()
215 .map(|(spini, cs_spini)| {
216 cs_spini
217 .columns()
218 .into_iter()
219 .enumerate()
220 .map(move |(i, c)| {
221 MolecularOrbital::builder()
222 .coefficients(c.to_owned())
223 .energy(self.mo_energies.as_ref().map(|moes| moes[spini][i]))
224 .bao(self.bao)
225 .mol(self.mol)
226 .structure_constraint(self.structure_constraint.clone())
227 .component_index(spini)
228 .complex_symmetric(self.complex_symmetric)
229 .threshold(self.threshold)
230 .build()
231 .expect("Unable to construct a molecular orbital.")
232 })
233 .collect::<Vec<_>>()
234 })
235 .collect::<Vec<_>>()
236 }
237}
238
239impl<'a, T, SC> SlaterDeterminant<'a, T, SC>
240where
241 T: ComplexFloat + Clone + Lapack,
242 SC: StructureConstraint + fmt::Display,
243{
244 pub fn complex_symmetric(&self) -> bool {
246 self.complex_symmetric
247 }
248
249 pub(crate) fn complex_conjugated(&self) -> bool {
251 self.complex_conjugated
252 }
253
254 pub fn structure_constraint(&self) -> &SC {
256 &self.structure_constraint
257 }
258
259 pub fn bao(&self) -> &BasisAngularOrder {
262 self.bao
263 }
264
265 pub fn mol(&self) -> &Molecule {
267 self.mol
268 }
269
270 pub fn energy(&self) -> Result<&T, &String> {
272 self.energy.as_ref()
273 }
274
275 pub fn mo_energies(&self) -> Option<&Vec<Array1<T>>> {
277 self.mo_energies.as_ref()
278 }
279
280 pub fn occupations(&self) -> &Vec<Array1<<T as ComplexFloat>::Real>> {
282 &self.occupations
283 }
284
285 pub fn coefficients(&self) -> &Vec<Array2<T>> {
287 &self.coefficients
288 }
289
290 pub fn threshold(&self) -> <T as ComplexFloat>::Real {
292 self.threshold
293 }
294
295 pub fn nelectrons(&self) -> <T as ComplexFloat>::Real
297 where
298 <T as ComplexFloat>::Real: Sum + From<u16>,
299 {
300 let implicit_factor = self
301 .structure_constraint
302 .implicit_factor()
303 .expect("Unable to retrieve the implicit factor from the structure constraint.");
304 <T as ComplexFloat>::Real::from(
305 implicit_factor
306 .to_u16()
307 .expect("Unable to convert the implicit factor to `u16`."),
308 ) * self
309 .occupations
310 .iter()
311 .map(|occ| occ.iter().copied().sum())
312 .sum()
313 }
314}
315
316impl<'a, T> SlaterDeterminant<'a, T, SpinConstraint>
317where
318 T: ComplexFloat + Clone + Lapack,
319{
320 pub fn to_generalised(&self) -> Self {
328 match self.structure_constraint {
329 SpinConstraint::Restricted(n) => {
330 log::debug!(
331 "Restricted Slater determinant will be augmented to generalised Slater determinant."
332 );
333 let nbas = self.bao.n_funcs();
334
335 let cr = &self.coefficients[0];
336 let occr = &self.occupations[0];
337 let norb = cr.ncols();
338 let mut cg = Array2::<T>::zeros((nbas * usize::from(n), norb * usize::from(n)));
339 let mut occg = Array1::<<T as ComplexFloat>::Real>::zeros((norb * usize::from(n),));
340 (0..usize::from(n)).for_each(|i| {
341 let row_start = nbas * i;
342 let row_end = nbas * (i + 1);
343 let col_start = norb * i;
344 let col_end = norb * (i + 1);
345 cg.slice_mut(s![row_start..row_end, col_start..col_end])
346 .assign(cr);
347 occg.slice_mut(s![col_start..col_end]).assign(occr);
348 });
349 let moeg_opt = self.mo_energies.as_ref().map(|moer| {
350 let mut moeg = Array1::<T>::zeros((norb * usize::from(n),));
351 (0..usize::from(n)).for_each(|i| {
352 let col_start = norb * i;
353 let col_end = norb * (i + 1);
354 moeg.slice_mut(s![col_start..col_end]).assign(&moer[0]);
355 });
356 vec![moeg]
357 });
358 Self::builder()
359 .coefficients(&[cg])
360 .occupations(&[occg])
361 .mo_energies(moeg_opt)
362 .energy(self.energy.clone())
363 .bao(self.bao)
364 .mol(self.mol)
365 .structure_constraint(SpinConstraint::Generalised(n, false))
366 .complex_symmetric(self.complex_symmetric)
367 .threshold(self.threshold)
368 .build()
369 .expect("Unable to spin-generalise a `SlaterDeterminant`.")
370 }
371 SpinConstraint::Unrestricted(n, increasingm) => {
372 log::debug!(
373 "Unrestricted Slater determinant will be augmented to generalised Slater determinant."
374 );
375 let nbas = self.bao.n_funcs();
376 let norb_tot = self.coefficients.iter().map(|c| c.ncols()).sum();
377 let mut cg = Array2::<T>::zeros((nbas * usize::from(n), norb_tot));
378 let mut occg = Array1::<<T as ComplexFloat>::Real>::zeros((norb_tot,));
379
380 let col_boundary_indices = (0..usize::from(n))
381 .scan(0, |acc, ispin| {
382 let start_index = *acc;
383 *acc += self.coefficients[ispin].shape()[1];
384 Some((start_index, *acc))
385 })
386 .collect::<Vec<_>>();
387 (0..usize::from(n)).for_each(|i| {
388 let row_start = nbas * i;
389 let row_end = nbas * (i + 1);
390 let (col_start, col_end) = col_boundary_indices[i];
391 cg.slice_mut(s![row_start..row_end, col_start..col_end])
392 .assign(&self.coefficients[i]);
393 occg.slice_mut(s![col_start..col_end])
394 .assign(&self.occupations[i]);
395 });
396
397 let moeg_opt = self.mo_energies.as_ref().map(|moer| {
398 let mut moeg = Array1::<T>::zeros((norb_tot,));
399 (0..usize::from(n)).for_each(|i| {
400 let (col_start, col_end) = col_boundary_indices[i];
401 moeg.slice_mut(s![col_start..col_end]).assign(&moer[i]);
402 });
403 vec![moeg]
404 });
405
406 Self::builder()
407 .coefficients(&[cg])
408 .occupations(&[occg])
409 .mo_energies(moeg_opt)
410 .energy(self.energy.clone())
411 .bao(self.bao)
412 .mol(self.mol)
413 .structure_constraint(SpinConstraint::Generalised(n, increasingm))
414 .complex_symmetric(self.complex_symmetric)
415 .threshold(self.threshold)
416 .build()
417 .expect("Unable to spin-generalise a `SlaterDeterminant`.")
418 }
419 SpinConstraint::Generalised(_, _) => self.clone(),
420 }
421 }
422}
423
424impl<'a> SlaterDeterminant<'a, f64, SpinConstraint> {
425 pub fn to_densities(
437 &'a self,
438 ) -> Result<DensitiesOwned<'a, f64, SpinConstraint>, anyhow::Error> {
439 let densities = match self.structure_constraint {
440 SpinConstraint::Restricted(_) | SpinConstraint::Unrestricted(_, _) => {
441 self.coefficients().iter().zip(self.occupations().iter()).map(|(c, o)| {
442 let denmat = einsum(
443 "i,mi,ni->mn",
444 &[&o.view(), &c.view(), &c.view()]
445 )
446 .expect("Unable to construct a density matrix from a determinant coefficient matrix.")
447 .into_dimensionality::<Ix2>()
448 .expect("Unable to convert the resultant density matrix to two dimensions.");
449 Density::<f64>::builder()
450 .density_matrix(denmat)
451 .bao(self.bao())
452 .mol(self.mol())
453 .complex_symmetric(self.complex_symmetric())
454 .threshold(self.threshold())
455 .build()
456 })
457 .collect::<Result<Vec<_>, _>>()?
458 }
459 SpinConstraint::Generalised(nspins, _) => {
460 let denmat = einsum(
461 "i,mi,ni->mn",
462 &[&self.occupations[0].view(), &self.coefficients[0].view(), &self.coefficients[0].view()]
463 )
464 .expect("Unable to construct a density matrix from a determinant coefficient matrix.")
465 .into_dimensionality::<Ix2>()
466 .expect("Unable to convert the resultant density matrix to two dimensions.");
467 let nspatial = self.bao.n_funcs();
468 (0..usize::from(nspins)).map(|ispin| {
469 let ispin_denmat = denmat.slice(
470 s![ispin*nspatial..(ispin + 1)*nspatial, ispin*nspatial..(ispin + 1)*nspatial]
471 ).to_owned();
472 Density::<f64>::builder()
473 .density_matrix(ispin_denmat)
474 .bao(self.bao())
475 .mol(self.mol())
476 .complex_symmetric(self.complex_symmetric())
477 .threshold(self.threshold())
478 .build()
479 }).collect::<Result<Vec<_>, _>>()?
480 }
481 };
482 DensitiesOwned::builder()
483 .structure_constraint(self.structure_constraint.clone())
484 .densities(densities)
485 .build()
486 .map_err(|err| format_err!(err))
487 }
488}
489
490impl<'a, T> SlaterDeterminant<'a, Complex<T>, SpinConstraint>
491where
492 T: Float + FloatConst + Lapack,
493 Complex<T>: Lapack,
494{
495 pub fn to_densities(
512 &'a self,
513 ) -> Result<DensitiesOwned<'a, Complex<T>, SpinConstraint>, anyhow::Error> {
514 let densities = match self.structure_constraint {
515 SpinConstraint::Restricted(_) | SpinConstraint::Unrestricted(_, _) => {
516 self.coefficients().iter().zip(self.occupations().iter()).map(|(c, o)| {
517 let denmat = einsum(
518 "i,mi,ni->mn",
519 &[&o.map(Complex::<T>::from).view(), &c.view(), &c.map(Complex::conj).view()]
520 )
521 .expect("Unable to construct a density matrix from a determinant coefficient matrix.")
522 .into_dimensionality::<Ix2>()
523 .expect("Unable to convert the resultant density matrix to two dimensions.");
524 Density::<Complex<T>>::builder()
525 .density_matrix(denmat)
526 .bao(self.bao())
527 .mol(self.mol())
528 .complex_symmetric(self.complex_symmetric())
529 .threshold(self.threshold())
530 .build()
531 })
532 .collect::<Result<Vec<_>, _>>()?
533 }
534 SpinConstraint::Generalised(nspins, _) => {
535 let denmat = einsum(
536 "i,mi,ni->mn",
537 &[
538 &self.occupations[0].map(Complex::<T>::from).view(),
539 &self.coefficients[0].view(),
540 &self.coefficients[0].map(Complex::conj).view()
541 ]
542 )
543 .expect("Unable to construct a density matrix from a determinant coefficient matrix.")
544 .into_dimensionality::<Ix2>()
545 .expect("Unable to convert the resultant density matrix to two dimensions.");
546 let nspatial = self.bao.n_funcs();
547 (0..usize::from(nspins)).map(|ispin| {
548 let ispin_denmat = denmat.slice(
549 s![ispin*nspatial..(ispin + 1)*nspatial, ispin*nspatial..(ispin + 1)*nspatial]
550 ).to_owned();
551 Density::<Complex<T>>::builder()
552 .density_matrix(ispin_denmat)
553 .bao(self.bao())
554 .mol(self.mol())
555 .complex_symmetric(self.complex_symmetric())
556 .threshold(self.threshold())
557 .build()
558 }).collect::<Result<Vec<_>, _>>()?
559 }
560 };
561 DensitiesOwned::builder()
562 .structure_constraint(self.structure_constraint.clone())
563 .densities(densities)
564 .build()
565 .map_err(|err| format_err!(err))
566 }
567}
568
569impl<'a, T> SlaterDeterminant<'a, Complex<T>, SpinOrbitCoupled>
570where
571 T: Float + FloatConst + Lapack,
572 Complex<T>: Lapack,
573{
574 pub fn to_densities(
587 &'a self,
588 ) -> Result<DensitiesOwned<'a, Complex<T>, SpinOrbitCoupled>, anyhow::Error> {
589 let densities = match self.structure_constraint {
590 SpinOrbitCoupled::JAdapted(ncomps) => {
591 let denmat = einsum(
592 "i,mi,ni->mn",
593 &[
594 &self.occupations[0].map(Complex::<T>::from).view(),
595 &self.coefficients[0].view(),
596 &self.coefficients[0].map(Complex::conj).view()
597 ]
598 )
599 .expect("Unable to construct a density matrix from a determinant coefficient matrix.")
600 .into_dimensionality::<Ix2>()
601 .expect("Unable to convert the resultant density matrix to two dimensions.");
602 let nspatial = self.bao.n_funcs();
603 (0..usize::from(ncomps)).map(|icomp| {
604 let icomp_denmat = denmat.slice(
605 s![icomp*nspatial..(icomp + 1)*nspatial, icomp*nspatial..(icomp + 1)*nspatial]
606 ).to_owned();
607 Density::<Complex<T>>::builder()
608 .density_matrix(icomp_denmat)
609 .bao(self.bao())
610 .mol(self.mol())
611 .complex_symmetric(self.complex_symmetric())
612 .threshold(self.threshold())
613 .build()
614 }).collect::<Result<Vec<_>, _>>()?
615 }
616 };
617 DensitiesOwned::builder()
618 .structure_constraint(self.structure_constraint.clone())
619 .densities(densities)
620 .build()
621 .map_err(|err| format_err!(err))
622 }
623}
624
625impl<'a, T, SC> From<SlaterDeterminant<'a, T, SC>> for SlaterDeterminant<'a, Complex<T>, SC>
633where
634 T: Float + FloatConst + Lapack,
635 Complex<T>: Lapack,
636 SC: StructureConstraint + Clone + fmt::Display,
637{
638 fn from(value: SlaterDeterminant<'a, T, SC>) -> Self {
639 SlaterDeterminant::<'a, Complex<T>, SC>::builder()
640 .coefficients(
641 &value
642 .coefficients
643 .into_iter()
644 .map(|coeffs| coeffs.map(Complex::from))
645 .collect::<Vec<_>>(),
646 )
647 .occupations(&value.occupations)
648 .mo_energies(value.mo_energies.map(|moes| {
649 moes.iter()
650 .map(|moe| moe.map(Complex::from))
651 .collect::<Vec<_>>()
652 }))
653 .bao(value.bao)
654 .mol(value.mol)
655 .structure_constraint(value.structure_constraint)
656 .complex_symmetric(value.complex_symmetric)
657 .threshold(value.threshold)
658 .build()
659 .expect("Unable to complexify a `SlaterDeterminant`.")
660 }
661}
662
663impl<'a, T, SC> PartialEq for SlaterDeterminant<'a, T, SC>
667where
668 T: ComplexFloat<Real = f64> + Lapack,
669 SC: StructureConstraint + PartialEq + fmt::Display,
670{
671 fn eq(&self, other: &Self) -> bool {
672 let thresh = (self.threshold * other.threshold).sqrt();
673 let coefficients_eq =
674 self.coefficients.len() == other.coefficients.len()
675 && self.coefficients.iter().zip(other.coefficients.iter()).all(
676 |(scoeffs, ocoeffs)| {
677 approx::relative_eq!(
678 (scoeffs - ocoeffs)
679 .map(|x| ComplexFloat::abs(*x).powi(2))
680 .sum()
681 .sqrt(),
682 0.0,
683 epsilon = thresh,
684 max_relative = thresh,
685 )
686 },
687 );
688 let occs_eq = self.occupations.len() == other.occupations.len()
689 && self
690 .occupations
691 .iter()
692 .zip(other.occupations.iter())
693 .all(|(soccs, ooccs)| {
694 approx::relative_eq!(
695 (soccs - ooccs).map(|x| x.abs().powi(2)).sum().sqrt(),
696 0.0,
697 epsilon = thresh,
698 max_relative = thresh,
699 )
700 });
701 self.structure_constraint == other.structure_constraint
702 && self.bao == other.bao
703 && self.mol == other.mol
704 && coefficients_eq
705 && occs_eq
706 }
707}
708
709impl<'a, T, SC> Eq for SlaterDeterminant<'a, T, SC>
713where
714 T: ComplexFloat<Real = f64> + Lapack,
715 SC: StructureConstraint + Eq + fmt::Display,
716{
717}
718
719impl<'a, T, SC> fmt::Debug for SlaterDeterminant<'a, T, SC>
723where
724 T: fmt::Debug + ComplexFloat + Lapack,
725 <T as ComplexFloat>::Real: Sum + From<u16> + fmt::Debug,
726 SC: StructureConstraint + fmt::Debug + fmt::Display,
727{
728 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
729 write!(
730 f,
731 "SlaterDeterminant[{:?}: {:?} electrons, {} coefficient {} of dimensions {}]",
732 self.structure_constraint,
733 self.nelectrons(),
734 self.coefficients.len(),
735 if self.coefficients.len() != 1 {
736 "matrices"
737 } else {
738 "matrix"
739 },
740 self.coefficients
741 .iter()
742 .map(|coeff| coeff
743 .shape()
744 .iter()
745 .map(|x| x.to_string())
746 .collect::<Vec<_>>()
747 .join("×"))
748 .collect::<Vec<_>>()
749 .join(", ")
750 )?;
751 Ok(())
752 }
753}
754
755impl<'a, T, SC> fmt::Display for SlaterDeterminant<'a, T, SC>
759where
760 T: fmt::Display + ComplexFloat + Lapack,
761 <T as ComplexFloat>::Real: Sum + From<u16> + fmt::Display,
762 SC: StructureConstraint + fmt::Display,
763{
764 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
765 write!(
766 f,
767 "SlaterDeterminant[{}: {} electrons, {} coefficient {} of dimensions {}]",
768 self.structure_constraint,
769 self.nelectrons(),
770 self.coefficients.len(),
771 if self.coefficients.len() != 1 {
772 "matrices"
773 } else {
774 "matrix"
775 },
776 self.coefficients
777 .iter()
778 .map(|coeff| coeff
779 .shape()
780 .iter()
781 .map(|x| x.to_string())
782 .collect::<Vec<_>>()
783 .join("×"))
784 .collect::<Vec<_>>()
785 .join(", ")
786 )?;
787 Ok(())
788 }
789}