qsym2/target/noci/backend/denmat/mod.rs
1use anyhow::{self, format_err};
2use ndarray::{Array2, Ix2, ScalarOperand};
3use ndarray_einsum::einsum;
4use num_complex::ComplexFloat;
5
6use super::nonortho::LowdinPairedCoefficients;
7
8/// Calculates the unweighted codensity matrix between two Löwdin-paired spin-orbitals $`i`$ in
9/// determinants $`^{w}\Psi`$ and $`^{x}\Psi`$.
10///
11/// The unweighted codensity matrix between the above two spin-orbitals is given by
12/// ```math
13/// ^{wx}\mathbf{P}_i
14/// = \ ^{w}\mathbf{C}_i \otimes
15/// (^{x}\mathbf{G}_i^{\star})^{\lozenge}
16/// = \ ^{w}\mathbf{C}_i
17/// \ ^{x}\mathbf{C}_i^{\dagger\lozenge},
18/// ```
19/// where $`\mathbf{C}_i`$ is the coefficient column vector for spin-orbital $`\chi_i`$ such that
20///
21/// ```math
22/// \chi_i(\mathbf{x}) = \sum_{\mu} \psi_{\mu}(\mathbf{x}) C_{\mu i}.
23/// ```
24///
25/// # Arguments
26///
27/// * `lowdin_paired_coefficients` - Structure containing the Löwdin-paired coefficient matrices.
28/// * `i` - Index of the corresponding Löwdin-paired spin-orbitals for which the unweighted
29/// codensity matrix is to be computed.
30///
31/// # Returns
32///
33/// The unweighted codensity matrix $`^{wx}\mathbf{P}_i`$.
34pub fn calc_unweighted_codensity_matrix<T>(
35 lowdin_paired_coefficients: &LowdinPairedCoefficients<T>,
36 i: usize,
37) -> Result<Array2<T>, anyhow::Error>
38where
39 T: ComplexFloat + 'static,
40{
41 let (cwt, cxt) = lowdin_paired_coefficients.paired_coefficients();
42 let cwi = &cwt.column(i);
43 let cxi = &cxt.column(i);
44 if cwi.shape() != cxi.shape() {
45 Err(format_err!(
46 "Coefficient dimensions mismatched: cwi ({:?}) !~ cxi ({:?}).",
47 cwi.shape(),
48 cxi.shape()
49 ))
50 } else if lowdin_paired_coefficients.complex_symmetric() {
51 einsum("m,n->mn", &[cwi, cxi])
52 .map_err(|err| format_err!(err))
53 .and_then(|p| {
54 p.into_dimensionality::<Ix2>()
55 .map_err(|err| format_err!(err))
56 })
57 } else {
58 einsum("m,n->mn", &[cwi, &cxi.map(|x| x.conj())])
59 .map_err(|err| format_err!(err))
60 .and_then(|p| {
61 p.into_dimensionality::<Ix2>()
62 .map_err(|err| format_err!(err))
63 })
64 }
65}
66
67/// Calculates the weighted codensity matrix between a set of Löwdin-paired spin-orbitals in
68/// determinants $`^{w}\Psi`$ and $`^{x}\Psi`$.
69///
70/// The weighted codensity matrix described above is given by
71/// ```math
72/// ^{wx}\mathbf{W} =
73/// \sum_{\substack{i = 1\\ ^{wx}\lambda_i \neq 0}}^{N_{\mathrm{e}}}
74/// \frac{^{wx}\mathbf{P}_i}{^{wx}\lambda_i}
75/// ```
76/// where $`^{wx}\mathbf{P}_i`$ is the unweighted codensity matrix between Löwdin-paired spin-orbitals
77/// $`i`$ in determinants $`^{w}\Psi`$ and $`^{x}\Psi`$, and the sum runs over up to
78/// :math:`N_{\mathrm{e}}` pairs of Löwdin-paired spin-orbitals of consideration, excluding those
79/// that have zero Löwdin overlaps.
80///
81/// # Arguments
82///
83/// * `lowdin_paired_coefficients` - Structure containing the Löwdin-paired coefficient matrices.
84///
85/// # Returns
86///
87/// The weighted codensity matrix $`^{wx}\mathbf{W}`$.
88pub fn calc_weighted_codensity_matrix<T>(
89 lowdin_paired_coefficients: &LowdinPairedCoefficients<T>,
90) -> Result<Array2<T>, anyhow::Error>
91where
92 T: ComplexFloat + ScalarOperand,
93{
94 let nbasis = lowdin_paired_coefficients.nbasis();
95 lowdin_paired_coefficients.nonzero_indices().iter().fold(
96 Ok::<_, anyhow::Error>(Array2::<T>::zeros((nbasis, nbasis))),
97 |acc_res, &i| {
98 acc_res.and_then(|acc| {
99 Ok(acc
100 + calc_unweighted_codensity_matrix(lowdin_paired_coefficients, i)?
101 / *lowdin_paired_coefficients
102 .lowdin_overlaps()
103 .get(i)
104 .ok_or_else(|| {
105 format_err!("Unable to retrieve the Löwdin overlap with index {i}.")
106 })?)
107 })
108 },
109 )
110}