integral/
operator.rs

1//! One-electron operator-DSL driver (L3).
2//!
3//! [`Basis::int1e`] evaluates any one-electron [`Operator`] — a polynomial in the
4//! position `r` and momentum `p = −i∇` operators — over the basis, returning the
5//! complex matrix as separate real and imaginary parts ([`OperatorMatrix`]). The
6//! decomposition into base overlap integrals lives in [`integral_core::operator`];
7//! this driver loops shell pairs, contracts primitives, applies the
8//! Cartesian→spherical transform to **each** part, and places the blocks.
9//!
10//! ## Additive layer
11//!
12//! This is purely additive: the bespoke fast paths ([`Basis::overlap`],
13//! [`Basis::kinetic`], [`Basis::dipole`]) stay as-is. The DSL exists for
14//! generality (a new integral type is one [`Operator`] declaration, not an engine
15//! change) and is validated *against* those bespoke paths (`tests/operator_dsl`).
16//!
17//! ## Hermiticity / the imaginary unit
18//!
19//! `r`, `r_i r_j`, and `p·p` (kinetic) are real-symmetric → [`OperatorMatrix::imag`]
20//! is ~0; `p` and `r × p` are imaginary/antisymmetric → [`OperatorMatrix::real`]
21//! is ~0. The `i`/sign convention is documented in the module docs.
22
23use integral_core::operator::operator_into;
24use integral_core::os::MAX_L;
25
26pub use integral_core::operator::{Factor, Operator, Term};
27
28use crate::integrals::{place_block, to_func_1e};
29use crate::shell::Basis;
30use crate::IntegralError;
31
32/// The complex matrix of a one-electron operator, stored as separate real and
33/// imaginary parts. Both are dense row-major `nao × nao` (same ordering/strides
34/// as the value builders).
35///
36/// A Hermitian operator with a real Gaussian representation has
37/// [`OperatorMatrix::imag`] ≈ 0 (e.g. `r`, kinetic); an anti-Hermitian one (e.g.
38/// `p`, `r × p`) has [`OperatorMatrix::real`] ≈ 0 and an antisymmetric imaginary
39/// part.
40#[derive(Debug, Clone, PartialEq)]
41pub struct OperatorMatrix {
42    nao: usize,
43    re: Vec<f64>,
44    im: Vec<f64>,
45}
46
47impl OperatorMatrix {
48    /// Matrix dimension `nao`.
49    #[must_use]
50    pub fn nao(&self) -> usize {
51        self.nao
52    }
53
54    /// The real part `Re⟨μ|O|ν⟩` (row-major `nao × nao`).
55    #[must_use]
56    pub fn real(&self) -> &[f64] {
57        &self.re
58    }
59
60    /// The imaginary part `Im⟨μ|O|ν⟩` (row-major `nao × nao`).
61    #[must_use]
62    pub fn imag(&self) -> &[f64] {
63        &self.im
64    }
65
66    /// Largest `|Im⟨μ|O|ν⟩|` — ~0 for a real (Hermitian, real-symmetric)
67    /// operator. A diagnostic for the operator's Hermiticity character.
68    #[must_use]
69    pub fn max_abs_imag(&self) -> f64 {
70        self.im.iter().fold(0.0_f64, |m, &v| m.max(v.abs()))
71    }
72
73    /// Largest `|Re⟨μ|O|ν⟩|` — ~0 for a purely imaginary (anti-Hermitian)
74    /// operator.
75    #[must_use]
76    pub fn max_abs_real(&self) -> f64 {
77        self.re.iter().fold(0.0_f64, |m, &v| m.max(v.abs()))
78    }
79}
80
81impl Basis {
82    /// Evaluate a one-electron [`Operator`] over the basis, returning its complex
83    /// matrix as real + imaginary parts ([`OperatorMatrix`]).
84    ///
85    /// The operator is decomposed into base overlap integrals
86    /// ([`integral_core::operator`]); spherical shells are transformed to their
87    /// `2l+1` components (the transform is applied to each part). This is the
88    /// generic operator-DSL path; the bespoke [`Basis::overlap`] /
89    /// [`Basis::kinetic`] / [`Basis::dipole`] builders remain the fast paths.
90    ///
91    /// # Errors
92    /// [`IntegralError::OperatorMomentumTooHigh`] if any shell has
93    /// `l + op.degree() > MAX_L` (the DSL raises the ket to `l + degree`).
94    pub fn int1e(&self, op: &Operator) -> Result<OperatorMatrix, IntegralError> {
95        let deg = op.degree();
96        for s in self.shells() {
97            if s.l() + deg > MAX_L {
98                return Err(IntegralError::OperatorMomentumTooHigh {
99                    l: s.l(),
100                    degree: deg,
101                    max: MAX_L,
102                });
103            }
104        }
105
106        let n = self.nao();
107        let offs = self.offsets();
108        let mut re = vec![0.0; n * n];
109        let mut im = vec![0.0; n * n];
110        let shells = self.shells();
111
112        for (si, sa) in shells.iter().enumerate() {
113            for (sj, sb) in shells.iter().enumerate() {
114                let (na, nb) = (sa.n_cart(), sb.n_cart());
115                let mut br = vec![0.0; na * nb];
116                let mut bi = vec![0.0; na * nb];
117                for pi in 0..sa.n_prim() {
118                    for pj in 0..sb.n_prim() {
119                        let scale = sa.primitive_coeff(pi) * sb.primitive_coeff(pj);
120                        operator_into(op, sa.prim(pi), sb.prim(pj), scale, &mut br, &mut bi);
121                    }
122                }
123                // c2s each part independently (the transform is real & linear).
124                let br = to_func_1e(br, sa, sb);
125                let bi = to_func_1e(bi, sa, sb);
126                let nbf = sb.n_func();
127                place_block(&mut re, n, offs[si], offs[sj], &br, nbf);
128                place_block(&mut im, n, offs[si], offs[sj], &bi, nbf);
129            }
130        }
131        Ok(OperatorMatrix { nao: n, re, im })
132    }
133}
integral/operator.rs

integral/
operator.rs
