sci_form/scf/
core_matrices.rs

1//! Core one-electron matrices: S, T, V, and H⁰ = T + V.
2//!
3//! Groups all one-electron integrals into a single coherent structure
4//! for use throughout the SCF pipeline.
5
6use nalgebra::DMatrix;
7
8use super::basis::BasisSet;
9use super::kinetic_matrix::build_kinetic_matrix;
10use super::nuclear_matrix::build_nuclear_matrix;
11use super::overlap_matrix::build_overlap_matrix;
12
13/// All one-electron matrices needed for an SCF calculation.
14#[derive(Debug, Clone)]
15pub struct CoreMatrices {
16    /// Overlap matrix S (n_basis × n_basis).
17    pub overlap: DMatrix<f64>,
18    /// Kinetic energy matrix T (n_basis × n_basis).
19    pub kinetic: DMatrix<f64>,
20    /// Nuclear attraction matrix V (n_basis × n_basis).
21    pub nuclear: DMatrix<f64>,
22    /// Core Hamiltonian H⁰ = T + V (n_basis × n_basis).
23    pub core_hamiltonian: DMatrix<f64>,
24    /// Number of basis functions.
25    pub n_basis: usize,
26}
27
28impl CoreMatrices {
29    /// Build all one-electron matrices for the molecular system.
30    ///
31    /// Computes S, T, V, and H⁰ = T + V.
32    pub fn build(basis: &BasisSet, elements: &[u8], positions_bohr: &[[f64; 3]]) -> Self {
33        let s = build_overlap_matrix(basis);
34        let t = build_kinetic_matrix(basis);
35        let v = build_nuclear_matrix(basis, elements, positions_bohr);
36        let h_core = &t + &v;
37
38        CoreMatrices {
39            overlap: s,
40            kinetic: t,
41            nuclear: v,
42            core_hamiltonian: h_core,
43            n_basis: basis.n_basis,
44        }
45    }
46}
47
48/// Compute nuclear repulsion energy.
49///
50/// E_nuc = Σ_{A>B} Z_A · Z_B / |R_A - R_B|
51pub fn nuclear_repulsion_energy(elements: &[u8], positions_bohr: &[[f64; 3]]) -> f64 {
52    let n_atoms = elements.len();
53    let mut e_nuc = 0.0;
54
55    for a in 0..n_atoms {
56        for b in (a + 1)..n_atoms {
57            let za = elements[a] as f64;
58            let zb = elements[b] as f64;
59            let dx = positions_bohr[a][0] - positions_bohr[b][0];
60            let dy = positions_bohr[a][1] - positions_bohr[b][1];
61            let dz = positions_bohr[a][2] - positions_bohr[b][2];
62            let r = (dx * dx + dy * dy + dz * dz).sqrt();
63            e_nuc += za * zb / r;
64        }
65    }
66
67    e_nuc
68}
69
70#[cfg(test)]
71mod tests {
72    use super::*;
73
74    #[test]
75    fn test_core_hamiltonian_symmetric() {
76        let elements = [8u8, 1, 1];
77        let positions = [
78            [0.0, 0.0, 0.2216],
79            [0.0, 1.4313, -0.8864],
80            [0.0, -1.4313, -0.8864],
81        ];
82        let basis = BasisSet::sto3g(&elements, &positions);
83        let matrices = CoreMatrices::build(&basis, &elements, &positions);
84
85        let h = &matrices.core_hamiltonian;
86        for i in 0..matrices.n_basis {
87            for j in 0..matrices.n_basis {
88                assert!(
89                    (h[(i, j)] - h[(j, i)]).abs() < 1e-12,
90                    "H⁰ not symmetric at ({}, {})",
91                    i,
92                    j
93                );
94            }
95        }
96    }
97
98    #[test]
99    fn test_h_core_equals_t_plus_v() {
100        let basis = BasisSet::sto3g(&[1, 1], &[[0.0, 0.0, 0.0], [1.4, 0.0, 0.0]]);
101        let matrices = CoreMatrices::build(&basis, &[1, 1], &[[0.0, 0.0, 0.0], [1.4, 0.0, 0.0]]);
102
103        for i in 0..matrices.n_basis {
104            for j in 0..matrices.n_basis {
105                let expected = matrices.kinetic[(i, j)] + matrices.nuclear[(i, j)];
106                assert!(
107                    (matrices.core_hamiltonian[(i, j)] - expected).abs() < 1e-14,
108                    "H⁰ ≠ T + V at ({}, {})",
109                    i,
110                    j
111                );
112            }
113        }
114    }
115
116    #[test]
117    fn test_nuclear_repulsion_h2() {
118        let e_nuc = nuclear_repulsion_energy(&[1, 1], &[[0.0, 0.0, 0.0], [1.4, 0.0, 0.0]]);
119        assert!((e_nuc - 1.0 / 1.4).abs() < 1e-10);
120    }
121
122    #[test]
123    fn test_overlap_diagonal_positive() {
124        let elements = [8u8, 1, 1];
125        let positions = [
126            [0.0, 0.0, 0.2216],
127            [0.0, 1.4313, -0.8864],
128            [0.0, -1.4313, -0.8864],
129        ];
130        let basis = BasisSet::sto3g(&elements, &positions);
131        let matrices = CoreMatrices::build(&basis, &elements, &positions);
132
133        for i in 0..matrices.n_basis {
134            assert!(matrices.overlap[(i, i)] > 0.0, "S[{i},{i}] should be > 0");
135        }
136    }
137
138    #[test]
139    fn test_core_matrices_water_dimensions() {
140        let elements = [8u8, 1, 1];
141        let positions = [[0.0, 0.0, 0.0], [1.43, 1.11, 0.0], [-1.43, 1.11, 0.0]];
142        let basis = BasisSet::sto3g(&elements, &positions);
143        let matrices = CoreMatrices::build(&basis, &elements, &positions);
144
145        assert_eq!(matrices.n_basis, 7);
146        assert_eq!(matrices.overlap.nrows(), 7);
147        assert_eq!(matrices.kinetic.nrows(), 7);
148        assert_eq!(matrices.nuclear.nrows(), 7);
149        assert_eq!(matrices.core_hamiltonian.nrows(), 7);
150    }
151}
sci_form/scf/core_matrices.rs

sci_form/scf/
core_matrices.rs
