Skip to main content

sci_form/hf/
scf_trait.rs

1//! Unified SCF solver trait for HF, PM3, and xTB backends.
2//!
3//! Provides a common interface for self-consistent field procedures
4//! across different quantum chemistry methods.
5
6use nalgebra::DMatrix;
7
8/// Common SCF convergence result.
9#[derive(Debug, Clone)]
10pub struct ScfOutput {
11    /// Converged total energy (method-specific units).
12    pub energy: f64,
13    /// Orbital energies (eigenvalues).
14    pub orbital_energies: Vec<f64>,
15    /// Number of SCF iterations.
16    pub iterations: usize,
17    /// Whether SCF converged.
18    pub converged: bool,
19    /// Converged density matrix.
20    pub density: Option<DMatrix<f64>>,
21    /// Mulliken charges (if available).
22    pub mulliken_charges: Option<Vec<f64>>,
23}
24
25/// Configuration for SCF convergence.
26#[derive(Debug, Clone)]
27pub struct ScfConvergenceConfig {
28    /// Maximum number of SCF iterations.
29    pub max_iter: usize,
30    /// Energy convergence threshold.
31    pub energy_threshold: f64,
32    /// Density convergence threshold.
33    pub density_threshold: f64,
34    /// DIIS history size (0 to disable).
35    pub diis_size: usize,
36    /// Level shift for virtual orbitals (0.0 to disable).
37    pub level_shift: f64,
38    /// Enable ADIIS for initial iterations before switching to DIIS.
39    pub use_adiis: bool,
40    /// Iteration threshold for switching from ADIIS to DIIS.
41    pub adiis_switch_iter: usize,
42}
43
44impl Default for ScfConvergenceConfig {
45    fn default() -> Self {
46        ScfConvergenceConfig {
47            max_iter: 100,
48            energy_threshold: 1e-8,
49            density_threshold: 1e-6,
50            diis_size: 8,
51            level_shift: 0.0,
52            use_adiis: false,
53            adiis_switch_iter: 5,
54        }
55    }
56}
57
58/// Unified trait for SCF solvers (HF, PM3, xTB).
59pub trait ScfSolver {
60    /// Run the SCF procedure to convergence.
61    fn solve(&self, config: &ScfConvergenceConfig) -> Result<ScfOutput, String>;
62
63    /// Get the method name for display.
64    fn method_name(&self) -> &str;
65
66    /// Get the number of basis functions.
67    fn n_basis(&self) -> usize;
68
69    /// Get the number of electrons.
70    fn n_electrons(&self) -> usize;
71}
72
73/// HF SCF solver wrapping the existing Roothaan-Hall implementation.
74pub struct HfScfSolver {
75    pub h_core: DMatrix<f64>,
76    pub s_mat: DMatrix<f64>,
77    pub eris: Vec<f64>,
78    pub n_elec: usize,
79}
80
81impl ScfSolver for HfScfSolver {
82    fn solve(&self, config: &ScfConvergenceConfig) -> Result<ScfOutput, String> {
83        let hf_config = super::scf::ScfConfig {
84            max_iter: config.max_iter,
85            energy_threshold: config.energy_threshold,
86            density_threshold: config.density_threshold,
87            diis_size: config.diis_size,
88            level_shift: config.level_shift,
89        };
90
91        let result = super::scf::solve_scf(
92            &self.h_core,
93            &self.s_mat,
94            &self.eris,
95            self.n_elec,
96            &hf_config,
97        );
98
99        Ok(ScfOutput {
100            energy: result.energy,
101            orbital_energies: result.orbital_energies,
102            iterations: result.iterations,
103            converged: result.converged,
104            density: Some(result.density),
105            mulliken_charges: None,
106        })
107    }
108
109    fn method_name(&self) -> &str {
110        "HF"
111    }
112    fn n_basis(&self) -> usize {
113        self.h_core.nrows()
114    }
115    fn n_electrons(&self) -> usize {
116        self.n_elec
117    }
118}
119
120/// PM3 SCF solver wrapping the existing NDDO implementation.
121pub struct Pm3ScfSolver {
122    pub elements: Vec<u8>,
123    pub positions: Vec<[f64; 3]>,
124}
125
126impl ScfSolver for Pm3ScfSolver {
127    fn solve(&self, _config: &ScfConvergenceConfig) -> Result<ScfOutput, String> {
128        let result = crate::pm3::solver::solve_pm3(&self.elements, &self.positions)?;
129        Ok(ScfOutput {
130            energy: result.total_energy,
131            orbital_energies: result.orbital_energies,
132            iterations: result.scf_iterations,
133            converged: result.converged,
134            density: None,
135            mulliken_charges: Some(result.mulliken_charges),
136        })
137    }
138
139    fn method_name(&self) -> &str {
140        "PM3"
141    }
142    fn n_basis(&self) -> usize {
143        0
144    } // not directly accessible
145    fn n_electrons(&self) -> usize {
146        0
147    }
148}
149
150/// xTB SCC solver wrapping the existing tight-binding implementation.
151pub struct XtbScfSolver {
152    pub elements: Vec<u8>,
153    pub positions: Vec<[f64; 3]>,
154}
155
156impl ScfSolver for XtbScfSolver {
157    fn solve(&self, _config: &ScfConvergenceConfig) -> Result<ScfOutput, String> {
158        let result = crate::xtb::solver::solve_xtb(&self.elements, &self.positions)?;
159        Ok(ScfOutput {
160            energy: result.total_energy,
161            orbital_energies: result.orbital_energies,
162            iterations: result.scc_iterations,
163            converged: result.converged,
164            density: None,
165            mulliken_charges: Some(result.mulliken_charges),
166        })
167    }
168
169    fn method_name(&self) -> &str {
170        "xTB"
171    }
172    fn n_basis(&self) -> usize {
173        0
174    }
175    fn n_electrons(&self) -> usize {
176        0
177    }
178}