qsym2/target/noci/backend/matelem/
hamiltonian.rs

1use std::fmt::{Display, LowerExp};
2use std::marker::PhantomData;
3
4use anyhow::{self, ensure, format_err};
5use derive_builder::Builder;
6use itertools::Itertools;
7use ndarray::{Array2, ArrayView2, ArrayView4, Axis, ScalarOperand};
8use ndarray_linalg::types::Lapack;
9use num::FromPrimitive;
10use num_complex::ComplexFloat;
11
12use crate::angmom::spinor_rotation_3d::StructureConstraint;
13use crate::symmetry::symmetry_transformation::SymmetryTransformable;
14use crate::target::determinant::SlaterDeterminant;
15use crate::target::noci::backend::nonortho::{
16    calc_lowdin_pairing, calc_o0_matrix_element, calc_o1_matrix_element, calc_o2_matrix_element,
17};
18
19use super::OrbitMatrix;
20
21#[cfg(test)]
22#[path = "hamiltonian_tests.rs"]
23mod hamiltonian_tests;
24
25/// Structure for managing the electronic Hamiltonian integrals in an atomic-orbital basis.
26#[derive(Builder)]
27pub struct HamiltonianAO<'a, T, SC>
28where
29    T: ComplexFloat + Lapack,
30    SC: StructureConstraint + Clone,
31{
32    /// The nuclear repulsion energy.
33    enuc: T,
34
35    /// The one-electron integrals in an atomic-orbital basis (with respect to the specified
36    /// structure constraint).
37    onee: ArrayView2<'a, T>,
38
39    /// The two-electron integrals in an atomic-orbital basis (with respect to the specified
40    /// structure constraint).
41    twoe: ArrayView4<'a, T>,
42
43    /// The structure constraint for the wavefunctions described by this Hamiltonian.
44    #[builder(setter(skip), default = "PhantomData")]
45    structure_constraint: PhantomData<SC>,
46}
47
48impl<'a, T, SC> HamiltonianAO<'a, T, SC>
49where
50    T: ComplexFloat + Lapack,
51    SC: StructureConstraint + Clone,
52{
53    /// Returns a builder for [`HamiltonianAO`].
54    pub fn builder() -> HamiltonianAOBuilder<'a, T, SC> {
55        HamiltonianAOBuilder::<T, SC>::default()
56    }
57
58    /// Returns the nuclear repulsion energy.
59    pub fn enuc(&self) -> T {
60        self.enuc
61    }
62
63    /// Returns the one-electron integrals in an atomic-orbital basis (with respect to the specified
64    /// structure constraint).
65    pub fn onee(&'a self) -> &'a ArrayView2<'a, T> {
66        &self.onee
67    }
68
69    /// Returns the two-electron integrals in an atomic-orbital basis (with respect to the specified
70    /// structure constraint).
71    pub fn twoe(&'a self) -> &'a ArrayView4<'a, T> {
72        &self.twoe
73    }
74}
75
76impl<'a, T, SC> HamiltonianAO<'a, T, SC>
77where
78    T: ComplexFloat + Lapack + ScalarOperand + FromPrimitive,
79    <T as ComplexFloat>::Real: LowerExp,
80    SC: StructureConstraint + Display + PartialEq + Clone,
81{
82    /// Calculates the zero-, one-, and two-electron contributions to the matrix element of the
83    /// electronic Hamiltonian between two possibly non-orthogonal Slater determinants.
84    ///
85    /// The matrix element is given by
86    /// ```math
87    ///     \braket{\hat{\iota} ^{w}\Psi | \hat{\mathscr{H}} | ^{x}\Psi}
88    /// ```
89    /// where $`\hat{\iota}`$ is an involutory operator that is either the identity or the
90    /// complex-conjugation operator, which depends on whether the specified determinants have been
91    /// defined with the [`SlaterDeterminant::complex_symmetric`] boolean set to `false` or `true`,
92    /// respectively.
93    ///
94    /// # Arguments
95    ///
96    /// * `det_w` - The determinant $`^{w}\Psi`$.
97    /// * `det_x` - The determinant $`^{x}\Psi`$.
98    /// * `sao` - The atomic-orbital overlap matrix.
99    /// * `thresh_offdiag` - Threshold for determining non-zero off-diagonal elements in the
100    /// orbital overlap matrix between $`^{w}\Psi`$ and $`^{x}\Psi`$ during Löwdin pairing.
101    /// * `thresh_zeroov` - threshold for identifying zero Löwdin overlaps.
102    ///
103    /// # Returns
104    ///
105    /// A tuple containing the zero-, one-, and two-electron contributions to the matrix element.
106    pub fn calc_hamiltonian_matrix_element_contributions(
107        &self,
108        det_w: &SlaterDeterminant<T, SC>,
109        det_x: &SlaterDeterminant<T, SC>,
110        sao: &ArrayView2<T>,
111        thresh_offdiag: <T as ComplexFloat>::Real,
112        thresh_zeroov: <T as ComplexFloat>::Real,
113    ) -> Result<(T, T, T), anyhow::Error> {
114        ensure!(
115            det_w.structure_constraint() == det_x.structure_constraint(),
116            "Inconsistent spin constraints: {} != {}.",
117            det_w.structure_constraint(),
118            det_x.structure_constraint(),
119        );
120        let sc = det_w.structure_constraint();
121
122        if det_w.complex_symmetric() != det_x.complex_symmetric() {
123            return Err(format_err!(
124                "The `complex_symmetric` booleans of the specified determinants do not match: `det_w` (`{}`) != `det_x` (`{}`).",
125                det_w.complex_symmetric(),
126                det_x.complex_symmetric(),
127            ));
128        }
129        let complex_symmetric = det_w.complex_symmetric();
130
131        let lowdin_paired_coefficientss = det_w
132            .coefficients()
133            .iter()
134            .zip(det_w.occupations().iter())
135            .zip(det_x.coefficients().iter().zip(det_x.occupations().iter()))
136            .map(|((cw, occw), (cx, occx))| {
137                let occw_indices = occw
138                    .iter()
139                    .enumerate()
140                    .filter_map(|(i, occ_i)| {
141                        if occ_i.abs() >= det_w.threshold() {
142                            Some(i)
143                        } else {
144                            None
145                        }
146                    })
147                    .collect::<Vec<_>>();
148                let cw_occ = cw.select(Axis(1), &occw_indices);
149                let occx_indices = occx
150                    .iter()
151                    .enumerate()
152                    .filter_map(|(i, occ_i)| {
153                        if occ_i.abs() >= det_x.threshold() {
154                            Some(i)
155                        } else {
156                            None
157                        }
158                    })
159                    .collect::<Vec<_>>();
160                let cx_occ = cx.select(Axis(1), &occx_indices);
161                calc_lowdin_pairing(
162                    &cw_occ.view(),
163                    &cx_occ.view(),
164                    sao,
165                    complex_symmetric,
166                    thresh_offdiag,
167                    thresh_zeroov,
168                )
169            })
170            .collect::<Result<Vec<_>, _>>()?;
171
172        let zeroe_h_wx = calc_o0_matrix_element(&lowdin_paired_coefficientss, self.enuc, sc)?;
173        let onee_h_wx = calc_o1_matrix_element(&lowdin_paired_coefficientss, &self.onee, sc)?;
174        let twoe_h_wx = calc_o2_matrix_element(&lowdin_paired_coefficientss, &self.twoe, sc)?;
175        Ok((zeroe_h_wx, onee_h_wx, twoe_h_wx))
176    }
177
178    /// Calculates the Hamiltonian matrix in the basis of non-orthogonal Slater determinants.
179    ///
180    /// Each matrix element is given by
181    /// ```math
182    ///     \braket{\hat{\iota} ^{w}\Psi | \hat{\mathscr{H}} | ^{x}\Psi}
183    /// ```
184    /// where $`\hat{\iota}`$ is an involutory operator that is either the identity or the
185    /// complex-conjugation operator, which depends on whether the specified determinants have been
186    /// defined with the [`SlaterDeterminant::complex_symmetric`] boolean set to `false` or `true`,
187    /// respectively.
188    ///
189    /// # Arguments
190    ///
191    /// * `dets` - A sequence of Slater determinants to be used as the basis for the Hamiltonian
192    /// matrix.
193    /// * `sao` - The atomic-orbital overlap matrix.
194    /// * `thresh_offdiag` - Threshold for determining non-zero off-diagonal elements in the
195    /// orbital overlap matrix between $`^{w}\Psi`$ and $`^{x}\Psi`$ during Löwdin pairing.
196    /// * `thresh_zeroov` - threshold for identifying zero Löwdin overlaps.
197    ///
198    /// # Returns
199    ///
200    /// The Hamiltonian matrix.
201    pub fn calc_hamiltonian_matrix(
202        &self,
203        dets: &[&SlaterDeterminant<T, SC>],
204        sao: &ArrayView2<T>,
205        thresh_offdiag: <T as ComplexFloat>::Real,
206        thresh_zeroov: <T as ComplexFloat>::Real,
207    ) -> Result<Array2<T>, anyhow::Error> {
208        let dim = dets.len();
209        let mut hmat = Array2::<T>::zeros((dim, dim));
210        for pair in dets.iter().enumerate().combinations_with_replacement(2) {
211            let (w, det_w) = &pair[0];
212            let (x, det_x) = &pair[1];
213            let (zeroe_wx, onee_wx, twoe_wx) = self.calc_hamiltonian_matrix_element_contributions(
214                det_w,
215                det_x,
216                sao,
217                thresh_offdiag,
218                thresh_zeroov,
219            )?;
220            hmat[(*w, *x)] = zeroe_wx + onee_wx + twoe_wx;
221            if *w != *x {
222                let (zeroe_xw, onee_xw, twoe_xw) = self
223                    .calc_hamiltonian_matrix_element_contributions(
224                        det_x,
225                        det_w,
226                        sao,
227                        thresh_offdiag,
228                        thresh_zeroov,
229                    )?;
230                hmat[(*x, *w)] = zeroe_xw + onee_xw + twoe_xw;
231            }
232        }
233        Ok(hmat)
234    }
235}
236
237impl<'a, T, SC> OrbitMatrix<'a, T, SC> for &HamiltonianAO<'a, T, SC>
238where
239    T: ComplexFloat + Lapack + ScalarOperand + FromPrimitive,
240    <T as ComplexFloat>::Real: LowerExp,
241    SC: StructureConstraint + Clone + Display + PartialEq,
242    SlaterDeterminant<'a, T, SC>: SymmetryTransformable,
243{
244    fn calc_matrix_element(
245        &self,
246        det_w: &SlaterDeterminant<T, SC>,
247        det_x: &SlaterDeterminant<T, SC>,
248        sao: &ArrayView2<T>,
249        thresh_offdiag: <T as ComplexFloat>::Real,
250        thresh_zeroov: <T as ComplexFloat>::Real,
251    ) -> Result<T, anyhow::Error> {
252        let (zeroe, onee, twoe) = self.calc_hamiltonian_matrix_element_contributions(
253            det_w,
254            det_x,
255            sao,
256            thresh_offdiag,
257            thresh_zeroov,
258        )?;
259        Ok(zeroe + onee + twoe)
260    }
261}
262
263impl<'a, T, SC> OrbitMatrix<'a, T, SC> for HamiltonianAO<'a, T, SC>
264where
265    T: ComplexFloat + Lapack + ScalarOperand + FromPrimitive,
266    <T as ComplexFloat>::Real: LowerExp,
267    SC: StructureConstraint + Clone + Display + PartialEq,
268    SlaterDeterminant<'a, T, SC>: SymmetryTransformable,
269{
270    fn calc_matrix_element(
271        &self,
272        det_w: &SlaterDeterminant<T, SC>,
273        det_x: &SlaterDeterminant<T, SC>,
274        sao: &ArrayView2<T>,
275        thresh_offdiag: <T as ComplexFloat>::Real,
276        thresh_zeroov: <T as ComplexFloat>::Real,
277    ) -> Result<T, anyhow::Error> {
278        (&self).calc_matrix_element(det_w, det_x, sao, thresh_offdiag, thresh_zeroov)
279    }
280}
qsym2/target/noci/backend/matelem/hamiltonian.rs

qsym2/target/noci/backend/matelem/
hamiltonian.rs
