sci_form/scf/
orthogonalization.rs

1//! Löwdin and canonical orthogonalization of the overlap matrix.
2//!
3//! The generalized eigenvalue problem Fc = Scε is converted to a
4//! standard eigenvalue problem by orthogonalizing the basis:
5//!
6//!   F' = X†FX,  F'c' = c'ε,  c = Xc'
7//!
8//! where X = S^{-1/2} is the (symmetric) orthogonalization matrix.
9//!
10//! # Löwdin Orthogonalization
11//!
12//! S = UΛU†  →  S^{-1/2} = UΛ^{-1/2}U†
13//!
14//! This preserves the character of the original basis functions
15//! as much as possible (least-distortion orthogonalization).
16
17use nalgebra::DMatrix;
18
19/// Compute the symmetric orthogonalization matrix X = S^{-1/2}.
20///
21/// Uses eigendecomposition of S:
22///   S = U Λ U^T
23///   X = U Λ^{-1/2} U^T
24///
25/// Eigenvalues below `threshold` are discarded (linear dependency).
26/// Returns `(X, n_independent)`.
27pub fn lowdin_orthogonalization(s: &DMatrix<f64>, threshold: f64) -> (DMatrix<f64>, usize) {
28    let n = s.nrows();
29    let eigen = s.clone().symmetric_eigen();
30
31    // Count valid eigenvalues (above threshold)
32    let n_independent = eigen.eigenvalues.iter().filter(|&&e| e > threshold).count();
33
34    // Build X = U Λ^{-1/2} U^T
35    let mut x = DMatrix::zeros(n, n);
36
37    for k in 0..n {
38        let eigenval = eigen.eigenvalues[k];
39        if eigenval > threshold {
40            let inv_sqrt = 1.0 / eigenval.sqrt();
41            // X += inv_sqrt * u_k * u_k^T
42            for i in 0..n {
43                for j in 0..n {
44                    x[(i, j)] += inv_sqrt * eigen.eigenvectors[(i, k)] * eigen.eigenvectors[(j, k)];
45                }
46            }
47        }
48    }
49
50    (x, n_independent)
51}
52
53/// Transform a matrix to the orthogonal basis: M' = X† M X.
54pub fn transform_to_orthogonal(m: &DMatrix<f64>, x: &DMatrix<f64>) -> DMatrix<f64> {
55    x.transpose() * m * x
56}
57
58/// Back-transform coefficients from orthogonal to original basis: C = X C'.
59pub fn back_transform(c_prime: &DMatrix<f64>, x: &DMatrix<f64>) -> DMatrix<f64> {
60    x * c_prime
61}
62
63#[cfg(test)]
64mod tests {
65    use super::*;
66
67    #[test]
68    fn test_lowdin_identity() {
69        // S = I → X = I
70        let s = DMatrix::identity(3, 3);
71        let (x, n_indep) = lowdin_orthogonalization(&s, 1e-10);
72        assert_eq!(n_indep, 3);
73
74        for i in 0..3 {
75            for j in 0..3 {
76                let expected = if i == j { 1.0 } else { 0.0 };
77                assert!((x[(i, j)] - expected).abs() < 1e-10);
78            }
79        }
80    }
81
82    #[test]
83    fn test_lowdin_produces_orthonormal_basis() {
84        // Create a simple 2×2 overlap matrix
85        let mut s = DMatrix::identity(2, 2);
86        s[(0, 1)] = 0.5;
87        s[(1, 0)] = 0.5;
88
89        let (x, _) = lowdin_orthogonalization(&s, 1e-10);
90
91        // X†SX should equal I
92        let identity_check = x.transpose() * &s * &x;
93        for i in 0..2 {
94            for j in 0..2 {
95                let expected = if i == j { 1.0 } else { 0.0 };
96                assert!(
97                    (identity_check[(i, j)] - expected).abs() < 1e-10,
98                    "X†SX[{},{}] = {}, expected {}",
99                    i,
100                    j,
101                    identity_check[(i, j)],
102                    expected
103                );
104            }
105        }
106    }
107
108    #[test]
109    fn test_lowdin_singular_overlap() {
110        // Create a nearly-singular overlap matrix
111        let mut s = DMatrix::identity(3, 3);
112        s[(0, 1)] = 0.99;
113        s[(1, 0)] = 0.99;
114        s[(0, 2)] = 0.5;
115        s[(2, 0)] = 0.5;
116        s[(1, 2)] = 0.5;
117        s[(2, 1)] = 0.5;
118
119        let (x, n_indep) = lowdin_orthogonalization(&s, 1e-6);
120
121        // Should handle near-singularity gracefully
122        assert!(n_indep <= 3);
123        assert!(n_indep >= 2);
124
125        // X†SX should still be (approximately) identity for retained dimensions
126        let check = x.transpose() * &s * &x;
127        for i in 0..check.nrows() {
128            // Diagonal should be close to 1
129            assert!(
130                (check[(i, i)] - 1.0).abs() < 1e-8 || check[(i, i)].abs() < 1e-8,
131                "Diagonal element [{},{}] = {} unexpected",
132                i,
133                i,
134                check[(i, i)]
135            );
136        }
137    }
138
139    #[test]
140    fn test_transform_and_back_transform() {
141        let m = DMatrix::from_row_slice(2, 2, &[1.0, 0.3, 0.3, 2.0]);
142        let x = DMatrix::from_row_slice(2, 2, &[0.8, 0.2, 0.2, 0.9]);
143
144        let m_prime = transform_to_orthogonal(&m, &x);
145        let c_prime = DMatrix::identity(2, 2);
146        let c = back_transform(&c_prime, &x);
147
148        // c should equal x for identity c'
149        assert!((c - &x).norm() < 1e-14);
150
151        // m' should be X†MX
152        let expected = x.transpose() * &m * &x;
153        assert!((m_prime - expected).norm() < 1e-14);
154    }
155}
sci_form/scf/orthogonalization.rs

sci_form/scf/
orthogonalization.rs
