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}