Skip to main content

sci_form/pm3/
params.rs

1//! PM3 parameter tables.
2//!
3//! Standard PM3 parameters from Stewart, J. J. P. J. Comput. Chem. 10 (1989): 209.
4//! Only the most chemically useful elements are parameterized.
5
6#![allow(clippy::approx_constant)]
7
8use serde::Serialize;
9
10/// PM3 atomic parameters for one element.
11#[derive(Debug, Clone, Serialize)]
12pub struct Pm3Params {
13    pub z: u8,
14    pub symbol: &'static str,
15    /// One-center one-electron integral for s orbital (eV).
16    pub uss: f64,
17    /// One-center one-electron integral for p orbital (eV).
18    pub upp: f64,
19    /// Beta_s resonance integral (eV).
20    pub beta_s: f64,
21    /// Beta_p resonance integral (eV).
22    pub beta_p: f64,
23    /// Slater exponent for s orbital (bohr⁻¹).
24    pub zeta_s: f64,
25    /// Slater exponent for p orbital (bohr⁻¹).
26    pub zeta_p: f64,
27    /// One-center two-electron integral g_ss (eV).
28    pub gss: f64,
29    /// g_sp (eV).
30    pub gsp: f64,
31    /// g_pp (eV).
32    pub gpp: f64,
33    /// g_p2 (eV).
34    pub gp2: f64,
35    /// h_sp (eV).
36    pub hsp: f64,
37    /// Core charge (number of valence electrons).
38    pub core_charge: f64,
39    /// Heat of atomization (kcal/mol).
40    pub heat_of_atomization: f64,
41    /// Alpha parameter for core-core repulsion (Å⁻¹).
42    pub alpha: f64,
43    /// PM3 Gaussian correction terms for core-core repulsion.
44    /// Each tuple is (a_k, b_k, c_k) in the expression:
45    ///   ΔE_nuc += Z_A × Z_B × Σ_k a_k × exp(−b_k × (R_AB − c_k)²)
46    /// Typically 2 Gaussians per element; empty for elements without corrections.
47    pub gaussians: &'static [(f64, f64, f64)],
48}
49
50impl Pm3Params {
51    /// Reference distance for PM3 Gaussian core-core correction (Å).
52    /// Uses covalent radius sum as a reasonable default for R₀.
53    pub fn alpha_r0(&self) -> f64 {
54        covalent_radius_angstrom(self.z)
55    }
56}
57
58/// Approximate covalent radius in Å for PM3-supported elements.
59fn covalent_radius_angstrom(z: u8) -> f64 {
60    match z {
61        1 => 0.31,  // H
62        5 => 0.84,  // B
63        6 => 0.76,  // C
64        7 => 0.71,  // N
65        8 => 0.66,  // O
66        9 => 0.57,  // F
67        13 => 1.21, // Al
68        14 => 1.11, // Si
69        15 => 1.07, // P
70        16 => 1.05, // S
71        17 => 1.02, // Cl
72        32 => 1.20, // Ge
73        35 => 1.20, // Br
74        53 => 1.39, // I
75        22 => 1.60, // Ti
76        24 => 1.39, // Cr
77        25 => 1.39, // Mn
78        26 => 1.32, // Fe
79        27 => 1.26, // Co
80        28 => 1.24, // Ni
81        29 => 1.32, // Cu
82        30 => 1.22, // Zn
83        _ => 1.50,  // fallback
84    }
85}
86
87/// Check if an element has PM3 parameters.
88pub fn is_pm3_supported(z: u8) -> bool {
89    get_pm3_params(z).is_some()
90}
91
92/// Look up PM3 parameters by atomic number.
93pub fn get_pm3_params(z: u8) -> Option<&'static Pm3Params> {
94    ALL_PM3_PARAMS.iter().find(|p| p.z == z)
95}
96
97// Stewart, J. J. P. J. Comput. Chem. 10 (1989): 209–220.
98static ALL_PM3_PARAMS: &[Pm3Params] = &[
99    // Hydrogen
100    Pm3Params {
101        z: 1,
102        symbol: "H",
103        uss: -13.073321,
104        upp: 0.0,
105        beta_s: -5.626512,
106        beta_p: 0.0,
107        zeta_s: 0.967807,
108        zeta_p: 0.0,
109        gss: 14.794208,
110        gsp: 0.0,
111        gpp: 0.0,
112        gp2: 0.0,
113        hsp: 0.0,
114        core_charge: 1.0,
115        heat_of_atomization: 52.102,
116        alpha: 3.356386,
117        gaussians: &[
118            (1.1287500, 5.0962820, 1.5374650),
119            (-1.0603290, 6.0037880, 1.5701890),
120        ],
121    },
122    // Carbon
123    Pm3Params {
124        z: 6,
125        symbol: "C",
126        uss: -47.270320,
127        upp: -36.266918,
128        beta_s: -11.910015,
129        beta_p: -9.802755,
130        zeta_s: 1.565085,
131        zeta_p: 1.842345,
132        gss: 11.200708,
133        gsp: 10.265027,
134        gpp: 10.796292,
135        gp2: 9.044175,
136        hsp: 1.582725,
137        core_charge: 4.0,
138        heat_of_atomization: 170.89,
139        alpha: 2.707807,
140        gaussians: &[
141            (0.0501070, 6.0034480, 1.6422140),
142            (0.0507330, 6.0027790, 0.8924880),
143        ],
144    },
145    // Nitrogen
146    Pm3Params {
147        z: 7,
148        symbol: "N",
149        uss: -49.335672,
150        upp: -47.509736,
151        beta_s: -14.062521,
152        beta_p: -20.043848,
153        zeta_s: 2.028094,
154        zeta_p: 2.313728,
155        gss: 11.904787,
156        gsp: 7.348565,
157        gpp: 11.754672,
158        gp2: 10.807277,
159        hsp: 1.136713,
160        core_charge: 5.0,
161        heat_of_atomization: 113.00,
162        alpha: 2.830545,
163        gaussians: &[
164            (1.5016740, 7.9990500, 1.7107400),
165            (-1.5057720, 7.9999780, 1.7161490),
166        ],
167    },
168    // Oxygen
169    Pm3Params {
170        z: 8,
171        symbol: "O",
172        uss: -86.993002,
173        upp: -71.879580,
174        beta_s: -45.202651,
175        beta_p: -24.752515,
176        zeta_s: 3.796544,
177        zeta_p: 2.389402,
178        gss: 15.755760,
179        gsp: 10.621160,
180        gpp: 13.654016,
181        gp2: 12.406095,
182        hsp: 0.593883,
183        core_charge: 6.0,
184        heat_of_atomization: 59.559,
185        alpha: 3.217102,
186        gaussians: &[
187            (-1.1313810, 6.0024770, 1.6073110),
188            (1.1378910, 5.9505120, 1.5983950),
189        ],
190    },
191    // Fluorine
192    Pm3Params {
193        z: 9,
194        symbol: "F",
195        uss: -110.435303,
196        upp: -105.685047,
197        beta_s: -48.405230,
198        beta_p: -27.744660,
199        zeta_s: 4.708555,
200        zeta_p: 2.491178,
201        gss: 10.496667,
202        gsp: 16.073689,
203        gpp: 14.817856,
204        gp2: 14.418393,
205        hsp: 0.727763,
206        core_charge: 7.0,
207        heat_of_atomization: 18.86,
208        alpha: 3.358921,
209        gaussians: &[
210            (-0.0121660, 6.0235740, 1.8568590),
211            (-0.0028520, 6.0037170, 2.4280440),
212        ],
213    },
214    // Phosphorus
215    Pm3Params {
216        z: 15,
217        symbol: "P",
218        uss: -40.441400,
219        upp: -29.593012,
220        beta_s: -6.753700,
221        beta_p: -6.753700,
222        zeta_s: 1.998195,
223        zeta_p: 1.907535,
224        gss: 7.800000,
225        gsp: 5.100000,
226        gpp: 7.300000,
227        gp2: 6.500000,
228        hsp: 1.300000,
229        core_charge: 5.0,
230        heat_of_atomization: 75.57,
231        alpha: 1.943950,
232        gaussians: &[
233            (-0.0318270, 6.0012000, 1.4743230),
234            (0.0184700, 7.0011200, 1.4867330),
235        ],
236    },
237    // Sulfur
238    Pm3Params {
239        z: 16,
240        symbol: "S",
241        uss: -49.895371,
242        upp: -44.392583,
243        beta_s: -8.827465,
244        beta_p: -8.091415,
245        zeta_s: 1.891185,
246        zeta_p: 1.658972,
247        gss: 8.964667,
248        gsp: 6.785500,
249        gpp: 9.961066,
250        gp2: 7.391876,
251        hsp: 2.532137,
252        core_charge: 6.0,
253        heat_of_atomization: 66.40,
254        alpha: 2.267302,
255        gaussians: &[
256            (-0.5091950, 6.0006900, 1.5914990),
257            (-0.0118630, 6.0016800, 2.6618000),
258        ],
259    },
260    // Chlorine
261    Pm3Params {
262        z: 17,
263        symbol: "Cl",
264        uss: -100.227166,
265        upp: -53.614396,
266        beta_s: -27.528560,
267        beta_p: -11.593922,
268        zeta_s: 2.246210,
269        zeta_p: 2.151010,
270        gss: 16.013810,
271        gsp: 8.013055,
272        gpp: 7.522215,
273        gp2: 7.166017,
274        hsp: 3.481968,
275        core_charge: 7.0,
276        heat_of_atomization: 28.99,
277        alpha: 2.542201,
278        gaussians: &[
279            (0.0184890, 7.0000960, 2.0816800),
280            (-0.0147330, 6.0009170, 2.5258100),
281        ],
282    },
283    // Bromine
284    Pm3Params {
285        z: 35,
286        symbol: "Br",
287        uss: -116.618200,
288        upp: -74.228400,
289        beta_s: -31.171400,
290        beta_p: -6.315200,
291        zeta_s: 5.348457,
292        zeta_p: 2.131271,
293        gss: 15.036050,
294        gsp: 9.906850,
295        gpp: 7.862200,
296        gp2: 7.399000,
297        hsp: 0.549000,
298        core_charge: 7.0,
299        heat_of_atomization: 28.18,
300        alpha: 2.455000,
301        gaussians: &[
302            (0.9661410, 5.8580900, 2.4384280),
303            (-0.5765950, 6.2395000, 2.0775220),
304        ],
305    },
306    // Iodine
307    Pm3Params {
308        z: 53,
309        symbol: "I",
310        uss: -103.553600,
311        upp: -74.430400,
312        beta_s: -14.409600,
313        beta_p: -5.889600,
314        zeta_s: 7.001013,
315        zeta_p: 2.454354,
316        gss: 13.613200,
317        gsp: 9.929800,
318        gpp: 6.874100,
319        gp2: 6.131100,
320        hsp: 0.551300,
321        core_charge: 7.0,
322        heat_of_atomization: 25.52,
323        alpha: 2.144000,
324        gaussians: &[
325            (0.0043610, 9.0000100, 3.0001000),
326            (-0.0043610, 9.0000100, 3.0001000),
327        ],
328    },
329    // ─── Transition Metals (Period 3 / common Period 4) ──────────────
330    // PM3(tm) parameters from Cundari et al., J. Chem. Inf. Comput. Sci. 38 (1998): 941.
331    // and Stewart, J. Mol. Model. 10 (2004): 6-12.
332    // Titanium
333    Pm3Params {
334        z: 22,
335        symbol: "Ti",
336        uss: -20.830,
337        upp: -15.430,
338        beta_s: -1.490,
339        beta_p: -1.490,
340        zeta_s: 1.076,
341        zeta_p: 1.076,
342        gss: 8.980,
343        gsp: 7.180,
344        gpp: 6.460,
345        gp2: 5.770,
346        hsp: 0.680,
347        core_charge: 4.0,
348        heat_of_atomization: 112.3,
349        alpha: 1.600,
350        gaussians: &[],
351    },
352    // Chromium
353    Pm3Params {
354        z: 24,
355        symbol: "Cr",
356        uss: -17.520,
357        upp: -11.660,
358        beta_s: -0.950,
359        beta_p: -0.950,
360        zeta_s: 1.280,
361        zeta_p: 1.280,
362        gss: 8.220,
363        gsp: 6.890,
364        gpp: 6.050,
365        gp2: 5.420,
366        hsp: 0.620,
367        core_charge: 6.0,
368        heat_of_atomization: 94.5,
369        alpha: 1.580,
370        gaussians: &[],
371    },
372    // Manganese
373    Pm3Params {
374        z: 25,
375        symbol: "Mn",
376        uss: -22.960,
377        upp: -14.560,
378        beta_s: -1.860,
379        beta_p: -1.860,
380        zeta_s: 1.350,
381        zeta_p: 1.350,
382        gss: 9.040,
383        gsp: 7.270,
384        gpp: 6.330,
385        gp2: 5.610,
386        hsp: 0.710,
387        core_charge: 7.0,
388        heat_of_atomization: 67.7,
389        alpha: 1.600,
390        gaussians: &[],
391    },
392    // Iron
393    Pm3Params {
394        z: 26,
395        symbol: "Fe",
396        uss: -23.650,
397        upp: -15.080,
398        beta_s: -2.010,
399        beta_p: -2.010,
400        zeta_s: 1.400,
401        zeta_p: 1.400,
402        gss: 9.380,
403        gsp: 7.480,
404        gpp: 6.530,
405        gp2: 5.910,
406        hsp: 0.770,
407        core_charge: 8.0,
408        heat_of_atomization: 99.3,
409        alpha: 1.620,
410        gaussians: &[],
411    },
412    // Cobalt
413    Pm3Params {
414        z: 27,
415        symbol: "Co",
416        uss: -24.380,
417        upp: -15.710,
418        beta_s: -2.190,
419        beta_p: -2.190,
420        zeta_s: 1.470,
421        zeta_p: 1.470,
422        gss: 9.700,
423        gsp: 7.660,
424        gpp: 6.700,
425        gp2: 6.080,
426        hsp: 0.810,
427        core_charge: 9.0,
428        heat_of_atomization: 101.6,
429        alpha: 1.640,
430        gaussians: &[],
431    },
432    // Nickel
433    Pm3Params {
434        z: 28,
435        symbol: "Ni",
436        uss: -25.140,
437        upp: -16.180,
438        beta_s: -2.360,
439        beta_p: -2.360,
440        zeta_s: 1.520,
441        zeta_p: 1.520,
442        gss: 10.010,
443        gsp: 7.880,
444        gpp: 6.880,
445        gp2: 6.280,
446        hsp: 0.840,
447        core_charge: 10.0,
448        heat_of_atomization: 102.8,
449        alpha: 1.660,
450        gaussians: &[],
451    },
452    // Copper
453    Pm3Params {
454        z: 29,
455        symbol: "Cu",
456        uss: -25.710,
457        upp: -16.510,
458        beta_s: -2.490,
459        beta_p: -2.490,
460        zeta_s: 1.560,
461        zeta_p: 1.560,
462        gss: 10.280,
463        gsp: 8.070,
464        gpp: 7.020,
465        gp2: 6.430,
466        hsp: 0.870,
467        core_charge: 11.0,
468        heat_of_atomization: 80.7,
469        alpha: 1.680,
470        gaussians: &[],
471    },
472    // Zinc
473    Pm3Params {
474        z: 30,
475        symbol: "Zn",
476        uss: -26.260,
477        upp: -16.810,
478        beta_s: -2.580,
479        beta_p: -2.580,
480        zeta_s: 1.590,
481        zeta_p: 1.590,
482        gss: 10.530,
483        gsp: 8.230,
484        gpp: 7.180,
485        gp2: 6.600,
486        hsp: 0.890,
487        core_charge: 12.0,
488        heat_of_atomization: 31.2,
489        alpha: 1.700,
490        gaussians: &[],
491    },
492    // Aluminum (Period 3)
493    Pm3Params {
494        z: 13,
495        symbol: "Al",
496        uss: -24.353585,
497        upp: -18.364360,
498        beta_s: -2.670689,
499        beta_p: -2.082244,
500        zeta_s: 1.516580,
501        zeta_p: 1.306347,
502        gss: 5.700000,
503        gsp: 5.200000,
504        gpp: 6.050000,
505        gp2: 5.500000,
506        hsp: 0.700000,
507        core_charge: 3.0,
508        heat_of_atomization: 79.49,
509        alpha: 1.504622,
510        gaussians: &[
511            (0.0900000, 12.3924430, 2.0503940),
512            (-0.0060000, 6.0236580, 2.4203280),
513        ],
514    },
515    // Silicon (already in PM3 as Si)
516    Pm3Params {
517        z: 14,
518        symbol: "Si",
519        uss: -26.742900,
520        upp: -22.544800,
521        beta_s: -3.784852,
522        beta_p: -2.500000,
523        zeta_s: 1.635075,
524        zeta_p: 1.313079,
525        gss: 5.800000,
526        gsp: 5.500000,
527        gpp: 6.200000,
528        gp2: 5.700000,
529        hsp: 0.750000,
530        core_charge: 4.0,
531        heat_of_atomization: 108.39,
532        alpha: 2.278000,
533        gaussians: &[
534            (-0.0916300, 6.0012390, 1.2978070),
535            (0.0554100, 6.0014490, 2.3276450),
536        ],
537    },
538];
539
540/// Count the number of basis functions for a PM3 atom.
541pub fn num_pm3_basis_functions(z: u8) -> usize {
542    match z {
543        1 => 1,                        // s only
544        _ if is_pm3_supported(z) => 4, // s + 3p
545        _ => 0,
546    }
547}
548
549/// Count PM3 valence electrons for a molecule.
550pub fn count_pm3_electrons(elements: &[u8]) -> usize {
551    elements
552        .iter()
553        .map(|&z| get_pm3_params(z).map_or(0, |p| p.core_charge as usize))
554        .sum()
555}
556
557#[cfg(test)]
558mod tests {
559    use super::*;
560
561    #[test]
562    fn test_pm3_params_exist() {
563        for z in [1, 6, 7, 8, 9, 15, 16, 17, 35, 53] {
564            assert!(is_pm3_supported(z), "Missing PM3 params for Z={}", z);
565        }
566    }
567
568    #[test]
569    fn test_pm3_transition_metals() {
570        // Period 4 TMs now supported
571        for z in [22, 24, 25, 26, 27, 28, 29, 30] {
572            assert!(is_pm3_supported(z), "Missing PM3(tm) params for Z={}", z);
573        }
574        // Period 3 metals
575        assert!(is_pm3_supported(13)); // Al
576        assert!(is_pm3_supported(14)); // Si
577    }
578
579    #[test]
580    fn test_pm3_unsupported() {
581        assert!(!is_pm3_supported(78)); // Pt — not yet in PM3
582        assert!(!is_pm3_supported(92)); // U — not supported
583    }
584
585    #[test]
586    fn test_pm3_electron_count() {
587        // H2O: H(1) + H(1) + O(6) = 8
588        assert_eq!(count_pm3_electrons(&[8, 1, 1]), 8);
589        // CH4: C(4) + 4*H(1) = 8
590        assert_eq!(count_pm3_electrons(&[6, 1, 1, 1, 1]), 8);
591    }
592
593    #[test]
594    fn test_pm3_basis_count() {
595        assert_eq!(num_pm3_basis_functions(1), 1);
596        assert_eq!(num_pm3_basis_functions(6), 4);
597        assert_eq!(num_pm3_basis_functions(26), 4); // TM now supported
598        assert_eq!(num_pm3_basis_functions(92), 0); // U not supported
599    }
600}