sci_form/hf/
basis.rs

1//! Gaussian basis set definitions for HF-3c.
2//!
3//! Minimal basis sets for Hartree-Fock: STO-3G contractions for s and p shells.
4//! Each primitive Gaussian has the form:
5//! $$g(\vec{r}) = N x^l y^m z^n \exp(-\alpha |\vec{r} - \vec{R}|^2)$$
6
7/// Angular momentum type.
8#[derive(Debug, Clone, Copy, PartialEq, Eq)]
9pub enum ShellType {
10    /// s-type (l=0): 1 function.
11    S,
12    /// p-type (l=1): 3 functions (px, py, pz).
13    P,
14}
15
16/// A contracted Gaussian shell centered on an atom.
17#[derive(Debug, Clone)]
18pub struct Shell {
19    /// Center index (atom index).
20    pub center_idx: usize,
21    /// Center position [x, y, z] in Å (converted to Bohr internally).
22    pub center: [f64; 3],
23    /// Shell type (S or P).
24    pub shell_type: ShellType,
25    /// Primitive exponents (α values in Bohr^-2).
26    pub exponents: Vec<f64>,
27    /// Contraction coefficients (normalized).
28    pub coefficients: Vec<f64>,
29}
30
31impl Shell {
32    /// Number of basis functions in this shell.
33    pub fn n_functions(&self) -> usize {
34        match self.shell_type {
35            ShellType::S => 1,
36            ShellType::P => 3,
37        }
38    }
39
40    /// Normalize contraction coefficients by multiplying by primitive norms
41    /// and then normalizing the entire contracted shell.
42    pub fn normalize(&mut self) {
43        // 1. Multiply by primitive norms
44        for (i, &alpha) in self.exponents.iter().enumerate() {
45            let n_prim = match self.shell_type {
46                ShellType::S => (2.0 * alpha / std::f64::consts::PI).powf(0.75),
47                ShellType::P => (128.0 * alpha.powi(5) / std::f64::consts::PI.powi(3)).powf(0.25),
48            };
49            self.coefficients[i] *= n_prim;
50        }
51
52        // 2. Compute self overlap of the unnormalized contracted shell
53        let mut sum_s = 0.0;
54        for (i, &a) in self.exponents.iter().enumerate() {
55            for (j, &b) in self.exponents.iter().enumerate() {
56                let gamma = a + b;
57                let overlap = match self.shell_type {
58                    ShellType::S => (std::f64::consts::PI / gamma).powf(1.5),
59                    ShellType::P => {
60                        let pre = (std::f64::consts::PI / gamma).powf(1.5);
61                        pre * (1.0 / (2.0 * gamma)) // <px | px>
62                    }
63                };
64                sum_s += self.coefficients[i] * self.coefficients[j] * overlap;
65            }
66        }
67
68        // 3. Normalize the final contracted shell
69        let norm_factor = 1.0 / sum_s.sqrt();
70        for c in &mut self.coefficients {
71            *c *= norm_factor;
72        }
73    }
74}
75
76/// Complete basis set for a molecular system.
77#[derive(Debug, Clone)]
78pub struct BasisSet {
79    /// All shells in the basis.
80    pub shells: Vec<Shell>,
81}
82
83impl BasisSet {
84    /// Total number of basis functions.
85    pub fn n_basis(&self) -> usize {
86        self.shells.iter().map(|s| s.n_functions()).sum()
87    }
88}
89
90/// Ångström to Bohr conversion factor.
91pub const ANG_TO_BOHR: f64 = 1.8897259886;
92
93/// Build an STO-3G basis set for the given atoms.
94pub fn build_sto3g_basis(elements: &[u8], positions: &[[f64; 3]]) -> BasisSet {
95    let mut shells = Vec::new();
96    for (idx, (&z, pos)) in elements.iter().zip(positions.iter()).enumerate() {
97        let center = [
98            pos[0] * ANG_TO_BOHR,
99            pos[1] * ANG_TO_BOHR,
100            pos[2] * ANG_TO_BOHR,
101        ];
102        let atom_shells = sto3g_shells(z, idx, center);
103        shells.extend(atom_shells);
104    }
105    BasisSet { shells }
106}
107
108/// STO-3G shell definitions per element.
109///
110/// Exponents and coefficients from Hehre, Stewart, Pople (1969).
111fn sto3g_shells(z: u8, center_idx: usize, center: [f64; 3]) -> Vec<Shell> {
112    let mut shells = match z {
113        1 => vec![Shell {
114            center_idx,
115            center,
116            shell_type: ShellType::S,
117            exponents: vec![3.42525091, 0.62391373, 0.16885540],
118            coefficients: vec![0.15432897, 0.53532814, 0.44463454],
119        }],
120        6 => vec![
121            // 1s core
122            Shell {
123                center_idx,
124                center,
125                shell_type: ShellType::S,
126                exponents: vec![71.6168370, 13.0450960, 3.5305122],
127                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
128            },
129            // 2s valence
130            Shell {
131                center_idx,
132                center,
133                shell_type: ShellType::S,
134                exponents: vec![2.9412494, 0.6834831, 0.2222899],
135                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
136            },
137            // 2p valence
138            Shell {
139                center_idx,
140                center,
141                shell_type: ShellType::P,
142                exponents: vec![2.9412494, 0.6834831, 0.2222899],
143                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
144            },
145        ],
146        7 => vec![
147            Shell {
148                center_idx,
149                center,
150                shell_type: ShellType::S,
151                exponents: vec![99.1061690, 18.0523120, 4.8856602],
152                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
153            },
154            Shell {
155                center_idx,
156                center,
157                shell_type: ShellType::S,
158                exponents: vec![3.7804559, 0.8784966, 0.2857144],
159                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
160            },
161            Shell {
162                center_idx,
163                center,
164                shell_type: ShellType::P,
165                exponents: vec![3.7804559, 0.8784966, 0.2857144],
166                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
167            },
168        ],
169        8 => vec![
170            Shell {
171                center_idx,
172                center,
173                shell_type: ShellType::S,
174                exponents: vec![130.7093200, 23.8088610, 6.4436083],
175                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
176            },
177            Shell {
178                center_idx,
179                center,
180                shell_type: ShellType::S,
181                exponents: vec![5.0331513, 1.1695961, 0.3803890],
182                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
183            },
184            Shell {
185                center_idx,
186                center,
187                shell_type: ShellType::P,
188                exponents: vec![5.0331513, 1.1695961, 0.3803890],
189                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
190            },
191        ],
192        9 => vec![
193            Shell {
194                center_idx,
195                center,
196                shell_type: ShellType::S,
197                exponents: vec![166.6791300, 30.3608120, 8.2168207],
198                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
199            },
200            Shell {
201                center_idx,
202                center,
203                shell_type: ShellType::S,
204                exponents: vec![6.4648032, 1.5022812, 0.4885885],
205                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
206            },
207            Shell {
208                center_idx,
209                center,
210                shell_type: ShellType::P,
211                exponents: vec![6.4648032, 1.5022812, 0.4885885],
212                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
213            },
214        ],
215        _ => {
216            // Fallback: hydrogen-like 1s with scaled exponent
217            vec![Shell {
218                center_idx,
219                center,
220                shell_type: ShellType::S,
221                exponents: vec![
222                    3.42525091 * (z as f64 / 1.0).powi(2),
223                    0.62391373 * (z as f64 / 1.0).powi(2),
224                    0.16885540 * (z as f64 / 1.0).powi(2),
225                ],
226                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
227            }]
228        }
229    };
230
231    for shell in &mut shells {
232        shell.normalize();
233    }
234
235    shells
236}
237
238#[cfg(test)]
239mod tests {
240    use super::*;
241
242    #[test]
243    fn test_h2_basis() {
244        let basis = build_sto3g_basis(&[1, 1], &[[0.0, 0.0, 0.0], [0.0, 0.0, 0.74]]);
245        assert_eq!(basis.n_basis(), 2); // Two 1s orbitals
246        assert_eq!(basis.shells.len(), 2);
247    }
248
249    #[test]
250    fn test_water_basis() {
251        let elements = [8u8, 1, 1];
252        let positions = [
253            [0.0, 0.0, 0.117],
254            [0.0, 0.757, -0.469],
255            [0.0, -0.757, -0.469],
256        ];
257        let basis = build_sto3g_basis(&elements, &positions);
258        // O: 1s + 2s + 2p(×3) = 5, 2×H: 1s each = 2 → total = 7
259        assert_eq!(basis.n_basis(), 7);
260    }
261}
sci_form/hf/basis.rs

sci_form/hf/
basis.rs
