sci_form/hf/
scf.rs

1//! Self-Consistent Field (SCF) solver with DIIS acceleration.
2//!
3//! Solves the Roothaan-Hall equations iteratively:
4//! $$\mathbf{FC} = \mathbf{SC}\boldsymbol{\varepsilon}$$
5//!
6//! Uses Löwdin orthogonalization ($S^{-1/2}$) to transform the generalized
7//! eigenvalue problem into a standard one, then applies DIIS for convergence.
8
9use super::fock::{build_fock, electronic_energy};
10use nalgebra::{DMatrix, DVector, SymmetricEigen};
11
12/// SCF convergence result.
13#[derive(Debug, Clone)]
14pub struct ScfResult {
15    /// Converged electronic energy (Hartree).
16    pub energy: f64,
17    /// Orbital energies (eigenvalues).
18    pub orbital_energies: Vec<f64>,
19    /// MO coefficient matrix C (columns are MOs).
20    pub coefficients: DMatrix<f64>,
21    /// Converged density matrix P.
22    pub density: DMatrix<f64>,
23    /// Number of SCF iterations.
24    pub iterations: usize,
25    /// Whether SCF converged.
26    pub converged: bool,
27}
28
29/// SCF solver configuration.
30pub struct ScfConfig {
31    pub max_iter: usize,
32    pub energy_threshold: f64,
33    pub density_threshold: f64,
34    pub diis_size: usize,
35}
36
37impl Default for ScfConfig {
38    fn default() -> Self {
39        ScfConfig {
40            max_iter: 100,
41            energy_threshold: 1e-8,
42            density_threshold: 1e-6,
43            diis_size: 6,
44        }
45    }
46}
47
48/// Run the SCF procedure.
49pub fn solve_scf(
50    h_core: &DMatrix<f64>,
51    s_mat: &DMatrix<f64>,
52    eris: &[f64],
53    n_electrons: usize,
54    config: &ScfConfig,
55) -> ScfResult {
56    let n = h_core.nrows();
57    let n_occ = n_electrons / 2;
58
59    // Löwdin orthogonalization: S^{-1/2}
60    let s_half_inv = lowdin_orthogonalization(s_mat);
61
62    // Initial guess: diagonalize H_core
63    let (mut energies, mut coeffs) = diagonalize_fock(h_core, &s_half_inv);
64    let mut density = build_density(&coeffs, n_occ);
65
66    let mut prev_energy = 0.0;
67    let mut converged = false;
68    let mut iterations = 0;
69
70    // DIIS storage
71    let mut diis_focks: Vec<DMatrix<f64>> = Vec::new();
72    let mut diis_errors: Vec<DMatrix<f64>> = Vec::new();
73
74    for iter in 0..config.max_iter {
75        iterations = iter + 1;
76
77        let fock = build_fock(h_core, &density, eris, n);
78        let energy = electronic_energy(&density, h_core, &fock);
79
80        // DIIS error: FPS - SPF
81        let error = &fock * &density * s_mat - s_mat * &density * &fock;
82        let error_norm = error.iter().map(|v| v * v).sum::<f64>().sqrt();
83
84        // DIIS extrapolation
85        diis_focks.push(fock.clone());
86        diis_errors.push(error);
87        if diis_focks.len() > config.diis_size {
88            diis_focks.remove(0);
89            diis_errors.remove(0);
90        }
91
92        let fock_diis = if diis_focks.len() >= 2 {
93            diis_extrapolate(&diis_focks, &diis_errors)
94        } else {
95            fock
96        };
97
98        let (new_energies, new_coeffs) = diagonalize_fock(&fock_diis, &s_half_inv);
99        let new_density = build_density(&new_coeffs, n_occ);
100
101        let de = (energy - prev_energy).abs();
102
103        if de < config.energy_threshold && error_norm < config.density_threshold {
104            converged = true;
105            energies = new_energies;
106            coeffs = new_coeffs;
107            density = new_density;
108            break;
109        }
110
111        prev_energy = energy;
112        energies = new_energies;
113        coeffs = new_coeffs;
114        density = new_density;
115    }
116
117    let final_energy = electronic_energy(&density, h_core, &build_fock(h_core, &density, eris, n));
118
119    ScfResult {
120        energy: final_energy,
121        orbital_energies: energies.as_slice().to_vec(),
122        coefficients: coeffs,
123        density,
124        iterations,
125        converged,
126    }
127}
128
129fn lowdin_orthogonalization(s: &DMatrix<f64>) -> DMatrix<f64> {
130    let eigen = SymmetricEigen::new(s.clone());
131    let n = s.nrows();
132    let mut s_inv_half = DMatrix::zeros(n, n);
133
134    for i in 0..n {
135        let val = eigen.eigenvalues[i];
136        if val > 1e-10 {
137            let factor = 1.0 / val.sqrt();
138            let col = eigen.eigenvectors.column(i);
139            s_inv_half += factor * col * col.transpose();
140        }
141    }
142    s_inv_half
143}
144
145fn diagonalize_fock(
146    fock: &DMatrix<f64>,
147    s_half_inv: &DMatrix<f64>,
148) -> (DVector<f64>, DMatrix<f64>) {
149    let f_prime = s_half_inv.transpose() * fock * s_half_inv;
150    let eigen = SymmetricEigen::new(f_prime);
151
152    // Sort eigenvalues/vectors by energy
153    let n = eigen.eigenvalues.len();
154    let mut indices: Vec<usize> = (0..n).collect();
155    indices.sort_by(|&a, &b| {
156        eigen.eigenvalues[a]
157            .partial_cmp(&eigen.eigenvalues[b])
158            .unwrap()
159    });
160
161    let sorted_energies = DVector::from_fn(n, |i, _| eigen.eigenvalues[indices[i]]);
162    let sorted_vecs = DMatrix::from_fn(n, n, |r, c| eigen.eigenvectors[(r, indices[c])]);
163
164    // Back-transform: C = S^{-1/2} C'
165    let coeffs = s_half_inv * sorted_vecs;
166    (sorted_energies, coeffs)
167}
168
169fn build_density(coeffs: &DMatrix<f64>, n_occ: usize) -> DMatrix<f64> {
170    let n = coeffs.nrows();
171    let mut density = DMatrix::zeros(n, n);
172    for i in 0..n_occ {
173        let col = coeffs.column(i);
174        density += 2.0 * &col * col.transpose();
175    }
176    density
177}
178
179fn diis_extrapolate(focks: &[DMatrix<f64>], errors: &[DMatrix<f64>]) -> DMatrix<f64> {
180    let m = errors.len();
181    let mut b = DMatrix::zeros(m + 1, m + 1);
182
183    for i in 0..m {
184        for j in 0..=i {
185            let bij: f64 = errors[i]
186                .iter()
187                .zip(errors[j].iter())
188                .map(|(a, b)| a * b)
189                .sum();
190            b[(i, j)] = bij;
191            b[(j, i)] = bij;
192        }
193    }
194    for i in 0..m {
195        b[(m, i)] = -1.0;
196        b[(i, m)] = -1.0;
197    }
198
199    let mut rhs = DVector::zeros(m + 1);
200    rhs[m] = -1.0;
201
202    // Solve B·c = rhs using pseudo-inverse
203    let svd = b.svd(true, true);
204    let c = match svd.solve(&rhs, 1e-10) {
205        Ok(c) => c,
206        Err(_) => {
207            // Fallback: use latest Fock
208            return focks.last().unwrap().clone();
209        }
210    };
211
212    let mut f_diis = DMatrix::zeros(focks[0].nrows(), focks[0].ncols());
213    for i in 0..m {
214        f_diis += c[i] * &focks[i];
215    }
216    f_diis
217}
218
219#[cfg(test)]
220mod tests {
221    use super::*;
222
223    #[test]
224    fn test_lowdin_identity() {
225        let s = DMatrix::identity(3, 3);
226        let s_inv = lowdin_orthogonalization(&s);
227        for i in 0..3 {
228            for j in 0..3 {
229                let expected = if i == j { 1.0 } else { 0.0 };
230                assert!(
231                    (s_inv[(i, j)] - expected).abs() < 1e-10,
232                    "S^{{-1/2}} of identity should be identity"
233                );
234            }
235        }
236    }
237}
sci_form/hf/scf.rs

sci_form/hf/
scf.rs
