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/// Map each AO index to its parent atom index.
91pub fn ao_to_atom_map(basis: &BasisSet) -> Vec<usize> {
92    let mut map = Vec::with_capacity(basis.n_basis());
93    for shell in &basis.shells {
94        for _ in 0..shell.n_functions() {
95            map.push(shell.center_idx);
96        }
97    }
98    map
99}
100
101/// Ångström to Bohr conversion factor.
102pub const ANG_TO_BOHR: f64 = 1.8897259886;
103
104/// Build an STO-3G basis set for the given atoms.
105pub fn build_sto3g_basis(elements: &[u8], positions: &[[f64; 3]]) -> BasisSet {
106    let mut shells = Vec::new();
107    for (idx, (&z, pos)) in elements.iter().zip(positions.iter()).enumerate() {
108        let center = [
109            pos[0] * ANG_TO_BOHR,
110            pos[1] * ANG_TO_BOHR,
111            pos[2] * ANG_TO_BOHR,
112        ];
113        let atom_shells = sto3g_shells(z, idx, center);
114        shells.extend(atom_shells);
115    }
116    BasisSet { shells }
117}
118
119/// STO-3G shell definitions per element.
120///
121/// Exponents and coefficients from Hehre, Stewart, Pople (1969).
122fn sto3g_shells(z: u8, center_idx: usize, center: [f64; 3]) -> Vec<Shell> {
123    let mut shells = match z {
124        1 => vec![Shell {
125            center_idx,
126            center,
127            shell_type: ShellType::S,
128            exponents: vec![3.42525091, 0.62391373, 0.16885540],
129            coefficients: vec![0.15432897, 0.53532814, 0.44463454],
130        }],
131        6 => vec![
132            // 1s core
133            Shell {
134                center_idx,
135                center,
136                shell_type: ShellType::S,
137                exponents: vec![71.6168370, 13.0450960, 3.5305122],
138                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
139            },
140            // 2s valence
141            Shell {
142                center_idx,
143                center,
144                shell_type: ShellType::S,
145                exponents: vec![2.9412494, 0.6834831, 0.2222899],
146                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
147            },
148            // 2p valence
149            Shell {
150                center_idx,
151                center,
152                shell_type: ShellType::P,
153                exponents: vec![2.9412494, 0.6834831, 0.2222899],
154                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
155            },
156        ],
157        7 => vec![
158            Shell {
159                center_idx,
160                center,
161                shell_type: ShellType::S,
162                exponents: vec![99.1061690, 18.0523120, 4.8856602],
163                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
164            },
165            Shell {
166                center_idx,
167                center,
168                shell_type: ShellType::S,
169                exponents: vec![3.7804559, 0.8784966, 0.2857144],
170                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
171            },
172            Shell {
173                center_idx,
174                center,
175                shell_type: ShellType::P,
176                exponents: vec![3.7804559, 0.8784966, 0.2857144],
177                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
178            },
179        ],
180        8 => vec![
181            Shell {
182                center_idx,
183                center,
184                shell_type: ShellType::S,
185                exponents: vec![130.7093200, 23.8088610, 6.4436083],
186                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
187            },
188            Shell {
189                center_idx,
190                center,
191                shell_type: ShellType::S,
192                exponents: vec![5.0331513, 1.1695961, 0.3803890],
193                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
194            },
195            Shell {
196                center_idx,
197                center,
198                shell_type: ShellType::P,
199                exponents: vec![5.0331513, 1.1695961, 0.3803890],
200                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
201            },
202        ],
203        9 => vec![
204            Shell {
205                center_idx,
206                center,
207                shell_type: ShellType::S,
208                exponents: vec![166.6791300, 30.3608120, 8.2168207],
209                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
210            },
211            Shell {
212                center_idx,
213                center,
214                shell_type: ShellType::S,
215                exponents: vec![6.4648032, 1.5022812, 0.4885885],
216                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
217            },
218            Shell {
219                center_idx,
220                center,
221                shell_type: ShellType::P,
222                exponents: vec![6.4648032, 1.5022812, 0.4885885],
223                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
224            },
225        ],
226        15 => vec![
227            Shell {
228                center_idx,
229                center,
230                shell_type: ShellType::S,
231                exponents: vec![508.291310, 92.5891370, 25.0571730],
232                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
233            },
234            Shell {
235                center_idx,
236                center,
237                shell_type: ShellType::S,
238                exponents: vec![18.5172490, 4.30422160, 1.39999670],
239                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
240            },
241            Shell {
242                center_idx,
243                center,
244                shell_type: ShellType::P,
245                exponents: vec![18.5172490, 4.30422160, 1.39999670],
246                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
247            },
248            Shell {
249                center_idx,
250                center,
251                shell_type: ShellType::S,
252                exponents: vec![1.56367880, 0.36368650, 0.11828520],
253                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
254            },
255            Shell {
256                center_idx,
257                center,
258                shell_type: ShellType::P,
259                exponents: vec![1.56367880, 0.36368650, 0.11828520],
260                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
261            },
262        ],
263        16 => vec![
264            Shell {
265                center_idx,
266                center,
267                shell_type: ShellType::S,
268                exponents: vec![598.642680, 109.046680, 29.5121090],
269                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
270            },
271            Shell {
272                center_idx,
273                center,
274                shell_type: ShellType::S,
275                exponents: vec![22.1564680, 5.14855900, 1.67441430],
276                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
277            },
278            Shell {
279                center_idx,
280                center,
281                shell_type: ShellType::P,
282                exponents: vec![22.1564680, 5.14855900, 1.67441430],
283                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
284            },
285            Shell {
286                center_idx,
287                center,
288                shell_type: ShellType::S,
289                exponents: vec![1.80579080, 0.41988400, 0.13655330],
290                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
291            },
292            Shell {
293                center_idx,
294                center,
295                shell_type: ShellType::P,
296                exponents: vec![1.80579080, 0.41988400, 0.13655330],
297                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
298            },
299        ],
300        17 => vec![
301            Shell {
302                center_idx,
303                center,
304                shell_type: ShellType::S,
305                exponents: vec![696.408520, 126.888800, 34.3399080],
306                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
307            },
308            Shell {
309                center_idx,
310                center,
311                shell_type: ShellType::S,
312                exponents: vec![25.9670530, 6.03406090, 1.96235810],
313                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
314            },
315            Shell {
316                center_idx,
317                center,
318                shell_type: ShellType::P,
319                exponents: vec![25.9670530, 6.03406090, 1.96235810],
320                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
321            },
322            Shell {
323                center_idx,
324                center,
325                shell_type: ShellType::S,
326                exponents: vec![2.14407210, 0.49841410, 0.16208590],
327                coefficients: vec![-0.09996723, 0.39951283, 0.70011547],
328            },
329            Shell {
330                center_idx,
331                center,
332                shell_type: ShellType::P,
333                exponents: vec![2.14407210, 0.49841410, 0.16208590],
334                coefficients: vec![0.15591627, 0.60768372, 0.39195739],
335            },
336        ],
337        _ => {
338            // Fallback: hydrogen-like 1s with scaled exponent
339            vec![Shell {
340                center_idx,
341                center,
342                shell_type: ShellType::S,
343                exponents: vec![
344                    3.42525091 * (z as f64 / 1.0).powi(2),
345                    0.62391373 * (z as f64 / 1.0).powi(2),
346                    0.16885540 * (z as f64 / 1.0).powi(2),
347                ],
348                coefficients: vec![0.15432897, 0.53532814, 0.44463454],
349            }]
350        }
351    };
352
353    for shell in &mut shells {
354        shell.normalize();
355    }
356
357    shells
358}
359
360#[cfg(test)]
361mod tests {
362    use super::*;
363
364    #[test]
365    fn test_h2_basis() {
366        let basis = build_sto3g_basis(&[1, 1], &[[0.0, 0.0, 0.0], [0.0, 0.0, 0.74]]);
367        assert_eq!(basis.n_basis(), 2); // Two 1s orbitals
368        assert_eq!(basis.shells.len(), 2);
369    }
370
371    #[test]
372    fn test_water_basis() {
373        let elements = [8u8, 1, 1];
374        let positions = [
375            [0.0, 0.0, 0.117],
376            [0.0, 0.757, -0.469],
377            [0.0, -0.757, -0.469],
378        ];
379        let basis = build_sto3g_basis(&elements, &positions);
380        // O: 1s + 2s + 2p(×3) = 5, 2×H: 1s each = 2 → total = 7
381        assert_eq!(basis.n_basis(), 7);
382    }
383}
sci_form/hf/basis.rs

sci_form/hf/
basis.rs
