sci_form/scf/
scf_loop.rs

1//! Main SCF iteration driver.
2//!
3//! Orchestrates the complete Roothaan-Hall SCF procedure:
4//!
5//! 1. Build one-electron matrices
6//! 2. Compute orthogonalization matrix
7//! 3. Initial guess (core Hamiltonian)
8//! 4. SCF loop with DIIS acceleration
9//! 5. Return converged result
10
11use nalgebra::DMatrix;
12
13use super::basis::BasisSet;
14use super::constants::{
15    DIIS_SUBSPACE_SIZE, HARTREE_TO_EV, SCF_DENSITY_THRESHOLD, SCF_ENERGY_THRESHOLD, SCF_MAX_ITER,
16};
17use super::core_matrices::{nuclear_repulsion_energy, CoreMatrices};
18use super::density_matrix::{build_density_matrix, density_rms_change};
19use super::diis::DiisAccelerator;
20use super::energy;
21use super::fock_matrix::build_fock_matrix;
22use super::mulliken::mulliken_analysis;
23use super::orthogonalization::{back_transform, lowdin_orthogonalization, transform_to_orthogonal};
24use super::two_electron::TwoElectronIntegrals;
25use super::types::{MolecularSystem, ScfResult};
26
27/// Configuration for the SCF calculation.
28#[derive(Debug, Clone)]
29pub struct ScfConfig {
30    /// Maximum number of SCF iterations.
31    pub max_iterations: usize,
32    /// Energy convergence threshold (Hartree).
33    pub energy_threshold: f64,
34    /// Density matrix RMS convergence threshold.
35    pub density_threshold: f64,
36    /// DIIS subspace size (0 = no DIIS).
37    pub diis_size: usize,
38    /// Level shift for virtual orbitals (Hartree).
39    pub level_shift: f64,
40    /// Damping factor for density mixing (0 = no damping, 0.5 = 50% old).
41    pub damping: f64,
42    /// Use rayon-parallelized two-electron integral computation.
43    pub use_parallel_eri: bool,
44    /// Basis set size threshold for parallel ERI.
45    pub parallel_threshold: usize,
46}
47
48impl Default for ScfConfig {
49    fn default() -> Self {
50        Self {
51            max_iterations: SCF_MAX_ITER,
52            energy_threshold: SCF_ENERGY_THRESHOLD,
53            density_threshold: SCF_DENSITY_THRESHOLD,
54            diis_size: DIIS_SUBSPACE_SIZE,
55            level_shift: 0.0,
56            damping: 0.0,
57            use_parallel_eri: false,
58            parallel_threshold: 20,
59        }
60    }
61}
62
63impl ScfConfig {
64    /// Return a config with rayon-parallel ERI enabled.
65    pub fn parallel() -> Self {
66        Self {
67            use_parallel_eri: true,
68            parallel_threshold: 0,
69            ..Self::default()
70        }
71    }
72}
73
74/// Run a complete Hartree-Fock SCF calculation.
75pub fn run_scf(system: &MolecularSystem, config: &ScfConfig) -> ScfResult {
76    let n_electrons = system.n_electrons();
77    let n_occupied = n_electrons / 2;
78
79    // Build basis set
80    let basis = BasisSet::sto3g(&system.atomic_numbers, &system.positions_bohr);
81    let n_basis = basis.n_basis;
82
83    // Build one-electron matrices
84    let core_matrices = CoreMatrices::build(&basis, &system.atomic_numbers, &system.positions_bohr);
85    let s = &core_matrices.overlap;
86    let h_core = &core_matrices.core_hamiltonian;
87
88    // Two-electron integrals
89    let use_parallel = config.use_parallel_eri && basis.n_basis >= config.parallel_threshold;
90    let eris = if use_parallel {
91        #[cfg(feature = "parallel")]
92        {
93            TwoElectronIntegrals::compute_parallel(&basis)
94        }
95        #[cfg(not(feature = "parallel"))]
96        {
97            TwoElectronIntegrals::compute(&basis)
98        }
99    } else {
100        TwoElectronIntegrals::compute(&basis)
101    };
102
103    // Nuclear repulsion
104    let e_nuc = nuclear_repulsion_energy(&system.atomic_numbers, &system.positions_bohr);
105
106    // Löwdin orthogonalization
107    let (x, _n_independent) = lowdin_orthogonalization(s, 1e-8);
108
109    // Initial guess — diagonalize H_core in orthogonal basis
110    let h_ortho = transform_to_orthogonal(h_core, &x);
111    let eigen = h_ortho.symmetric_eigen();
112
113    let mut indices: Vec<usize> = (0..n_basis).collect();
114    indices.sort_by(|&a, &b| {
115        eigen.eigenvalues[a]
116            .partial_cmp(&eigen.eigenvalues[b])
117            .unwrap()
118    });
119
120    let mut c_ortho = DMatrix::zeros(n_basis, n_basis);
121    let mut orbital_energies = vec![0.0; n_basis];
122    for (new_idx, &old_idx) in indices.iter().enumerate() {
123        orbital_energies[new_idx] = eigen.eigenvalues[old_idx];
124        for i in 0..n_basis {
125            c_ortho[(i, new_idx)] = eigen.eigenvectors[(i, old_idx)];
126        }
127    }
128
129    let mut c = back_transform(&c_ortho, &x);
130    let mut p = build_density_matrix(&c, n_occupied);
131    let mut p_old = p.clone();
132
133    // SCF loop
134    let mut diis = DiisAccelerator::new(config.diis_size);
135    let mut e_total = 0.0;
136    let mut converged = false;
137    let mut scf_iter = 0;
138    let mut fock = h_core.clone();
139
140    for iteration in 0..config.max_iterations {
141        scf_iter = iteration + 1;
142
143        fock = build_fock_matrix(h_core, &p, &eris);
144
145        // Level shift
146        if config.level_shift > 0.0 {
147            let n = fock.nrows();
148            for i in 0..n {
149                fock[(i, i)] += config.level_shift;
150            }
151        }
152
153        // DIIS extrapolation
154        if config.diis_size > 0 {
155            diis.add_iteration(&fock, &p, s);
156            if let Some(f_diis) = diis.extrapolate() {
157                fock = f_diis;
158            }
159        }
160
161        // Diagonalize in orthogonal basis
162        let f_ortho = transform_to_orthogonal(&fock, &x);
163        let eigen = f_ortho.symmetric_eigen();
164
165        let mut indices: Vec<usize> = (0..n_basis).collect();
166        indices.sort_by(|&a, &b| {
167            eigen.eigenvalues[a]
168                .partial_cmp(&eigen.eigenvalues[b])
169                .unwrap()
170        });
171
172        for (new_idx, &old_idx) in indices.iter().enumerate() {
173            orbital_energies[new_idx] = eigen.eigenvalues[old_idx];
174            for i in 0..n_basis {
175                c_ortho[(i, new_idx)] = eigen.eigenvectors[(i, old_idx)];
176            }
177        }
178
179        c = back_transform(&c_ortho, &x);
180        let p_new = build_density_matrix(&c, n_occupied);
181
182        if config.damping > 0.0 {
183            p = &p_new * (1.0 - config.damping) + &p_old * config.damping;
184        } else {
185            p = p_new;
186        }
187
188        let e_new = energy::total_energy(&p, h_core, &fock, e_nuc);
189        let delta_e = (e_new - e_total).abs();
190        let delta_p = density_rms_change(&p, &p_old);
191
192        if delta_e < config.energy_threshold && delta_p < config.density_threshold && iteration > 0
193        {
194            converged = true;
195            e_total = e_new;
196            break;
197        }
198
199        e_total = e_new;
200        p_old = p.clone();
201    }
202
203    let e_elec = energy::electronic_energy(&p, h_core, &fock);
204
205    let homo_energy = if n_occupied > 0 {
206        orbital_energies[n_occupied - 1]
207    } else {
208        0.0
209    };
210    let lumo_energy = if n_occupied < n_basis {
211        Some(orbital_energies[n_occupied])
212    } else {
213        None
214    };
215    let gap_ev = lumo_energy
216        .map(|lumo| (lumo - homo_energy) * HARTREE_TO_EV)
217        .unwrap_or(0.0);
218
219    let mulliken = mulliken_analysis(&p, s, &basis.function_to_atom, &system.atomic_numbers);
220
221    ScfResult {
222        orbital_energies,
223        mo_coefficients: c,
224        density_matrix: p,
225        electronic_energy: e_elec,
226        nuclear_repulsion: e_nuc,
227        total_energy: e_total,
228        homo_energy,
229        lumo_energy,
230        gap_ev,
231        mulliken_charges: mulliken.charges,
232        scf_iterations: scf_iter,
233        converged,
234        n_basis,
235        n_electrons,
236        overlap_matrix: s.clone(),
237        fock_matrix: fock,
238    }
239}
240
241#[cfg(test)]
242mod tests {
243    use super::super::constants::ANGSTROM_TO_BOHR;
244    use super::*;
245
246    #[test]
247    fn test_scf_h2() {
248        let system = MolecularSystem {
249            atomic_numbers: vec![1, 1],
250            positions_bohr: vec![[0.0, 0.0, 0.0], [0.74 * ANGSTROM_TO_BOHR, 0.0, 0.0]],
251            charge: 0,
252            multiplicity: 1,
253        };
254
255        let config = ScfConfig::default();
256        let result = run_scf(&system, &config);
257
258        assert!(result.total_energy < 0.0, "Total energy should be negative");
259        assert_eq!(result.n_basis, 2);
260        assert_eq!(result.n_electrons, 2);
261    }
262
263    #[test]
264    fn test_scf_converges() {
265        let system = MolecularSystem {
266            atomic_numbers: vec![1, 1],
267            positions_bohr: vec![[0.0, 0.0, 0.0], [1.4, 0.0, 0.0]],
268            charge: 0,
269            multiplicity: 1,
270        };
271
272        let result = run_scf(&system, &ScfConfig::default());
273        assert!(result.converged, "SCF should converge for H2");
274        assert!(result.scf_iterations < 50);
275    }
276}
sci_form/scf/scf_loop.rs

sci_form/scf/
scf_loop.rs
