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            None,
96            self.n_elec,
97            &hf_config,
98        );
99
100        Ok(ScfOutput {
101            energy: result.energy,
102            orbital_energies: result.orbital_energies,
103            iterations: result.iterations,
104            converged: result.converged,
105            density: Some(result.density),
106            mulliken_charges: None,
107        })
108    }
109
110    fn method_name(&self) -> &str {
111        "HF"
112    }
113    fn n_basis(&self) -> usize {
114        self.h_core.nrows()
115    }
116    fn n_electrons(&self) -> usize {
117        self.n_elec
118    }
119}
120
121/// PM3 SCF solver wrapping the existing NDDO implementation.
122pub struct Pm3ScfSolver {
123    pub elements: Vec<u8>,
124    pub positions: Vec<[f64; 3]>,
125}
126
127impl ScfSolver for Pm3ScfSolver {
128    fn solve(&self, _config: &ScfConvergenceConfig) -> Result<ScfOutput, String> {
129        let result = crate::pm3::solver::solve_pm3(&self.elements, &self.positions)?;
130        Ok(ScfOutput {
131            energy: result.total_energy,
132            orbital_energies: result.orbital_energies,
133            iterations: result.scf_iterations,
134            converged: result.converged,
135            density: None,
136            mulliken_charges: Some(result.mulliken_charges),
137        })
138    }
139
140    fn method_name(&self) -> &str {
141        "PM3"
142    }
143    fn n_basis(&self) -> usize {
144        0
145    } // not directly accessible
146    fn n_electrons(&self) -> usize {
147        0
148    }
149}
150
151/// xTB SCC solver wrapping the existing tight-binding implementation.
152pub struct XtbScfSolver {
153    pub elements: Vec<u8>,
154    pub positions: Vec<[f64; 3]>,
155}
156
157impl ScfSolver for XtbScfSolver {
158    fn solve(&self, _config: &ScfConvergenceConfig) -> Result<ScfOutput, String> {
159        let result = crate::xtb::solver::solve_xtb(&self.elements, &self.positions)?;
160        Ok(ScfOutput {
161            energy: result.total_energy,
162            orbital_energies: result.orbital_energies,
163            iterations: result.scc_iterations,
164            converged: result.converged,
165            density: None,
166            mulliken_charges: Some(result.mulliken_charges),
167        })
168    }
169
170    fn method_name(&self) -> &str {
171        "xTB"
172    }
173    fn n_basis(&self) -> usize {
174        0
175    }
176    fn n_electrons(&self) -> usize {
177        0
178    }
179}