sci_form/hf/
d3.rs

1//! Grimme DFT-D3 dispersion correction with Becke-Johnson damping.
2//!
3//! Adds long-range van der Waals interactions missing from HF:
4//! $$E_{disp} = -\sum_{n=6,8} s_n \sum_{A>B} \frac{C_n^{AB}}{R_{AB}^n + f_{damp}^n}$$
5//!
6//! Reference: Grimme, S. et al. J. Chem. Phys. 132 (2010): 154104.
7
8use serde::{Deserialize, Serialize};
9
10/// D3 correction result.
11#[derive(Debug, Clone, Serialize, Deserialize)]
12pub struct D3Result {
13    /// Total D3 dispersion energy (Hartree).
14    pub energy: f64,
15    /// C6 contribution.
16    pub e6: f64,
17    /// C8 contribution.
18    pub e8: f64,
19}
20
21/// D3 parameters for HF-3c (BJ damping).
22const S6: f64 = 1.0;
23const S8: f64 = 0.8777;
24const A1: f64 = 0.4171;
25const A2: f64 = 2.9149;
26
27/// Compute D3-BJ dispersion correction.
28pub fn compute_d3_energy(elements: &[u8], positions_bohr: &[[f64; 3]]) -> D3Result {
29    let n = elements.len();
30    let mut e6 = 0.0;
31    let mut e8 = 0.0;
32
33    for a in 0..n {
34        for b in (a + 1)..n {
35            let r = distance(positions_bohr, a, b);
36            if r < 1e-10 {
37                continue;
38            }
39
40            let c6 = get_c6(elements[a], elements[b]);
41            let c8 = c6_to_c8(c6, elements[a], elements[b]);
42
43            let r0 = (c8 / c6).sqrt();
44            let f6 = bj_damping(r, r0, 6);
45            let f8 = bj_damping(r, r0, 8);
46
47            e6 -= S6 * c6 / (r.powi(6) + f6);
48            e8 -= S8 * c8 / (r.powi(8) + f8);
49        }
50    }
51
52    D3Result {
53        energy: e6 + e8,
54        e6,
55        e8,
56    }
57}
58
59/// Becke-Johnson damping function.
60fn bj_damping(_r: f64, r0: f64, n: i32) -> f64 {
61    let r_cut = A1 * r0 + A2;
62    r_cut.powi(n)
63}
64
65/// Empirical C6 coefficients (Hartree·Bohr⁶).
66fn get_c6(z1: u8, z2: u8) -> f64 {
67    let c6_1 = atomic_c6(z1);
68    let c6_2 = atomic_c6(z2);
69    // Geometric combination rule
70    (c6_1 * c6_2).sqrt()
71}
72
73/// Approximate atomic C6 values.
74fn atomic_c6(z: u8) -> f64 {
75    match z {
76        1 => 6.50,
77        6 => 46.60,
78        7 => 24.20,
79        8 => 15.60,
80        9 => 9.52,
81        15 => 185.0,
82        16 => 134.0,
83        17 => 94.60,
84        _ => 30.0, // rough fallback
85    }
86}
87
88/// Estimate C8 from C6 using Casimir-Polder relation.
89fn c6_to_c8(c6: f64, z1: u8, z2: u8) -> f64 {
90    let q1 = atomic_q(z1);
91    let q2 = atomic_q(z2);
92    3.0 * c6 * (q1 * q2).sqrt()
93}
94
95/// Atomic quadrupole-like parameter for C6→C8.
96fn atomic_q(z: u8) -> f64 {
97    match z {
98        1 => 2.23,
99        6 => 7.40,
100        7 => 6.07,
101        8 => 5.14,
102        9 => 4.09,
103        _ => 5.0,
104    }
105}
106
107fn distance(positions: &[[f64; 3]], a: usize, b: usize) -> f64 {
108    let dx = positions[a][0] - positions[b][0];
109    let dy = positions[a][1] - positions[b][1];
110    let dz = positions[a][2] - positions[b][2];
111    (dx * dx + dy * dy + dz * dz).sqrt()
112}
113
114#[cfg(test)]
115mod tests {
116    use super::super::basis::ANG_TO_BOHR;
117    use super::*;
118
119    #[test]
120    fn test_d3_h2() {
121        let positions = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.74 * ANG_TO_BOHR]];
122        let result = compute_d3_energy(&[1, 1], &positions);
123        // D3 energy should be small and negative
124        assert!(result.energy < 0.0, "D3 should be attractive");
125        assert!(result.energy.abs() < 0.01, "D3 for H2 should be tiny");
126    }
127
128    #[test]
129    fn test_d3_ethane_larger() {
130        // Two carbons + H atoms → more D3
131        let r = 1.54 * ANG_TO_BOHR;
132        let positions = [[0.0, 0.0, 0.0], [0.0, 0.0, r]];
133        let result = compute_d3_energy(&[6, 6], &positions);
134        assert!(result.energy < 0.0);
135    }
136}
sci_form/hf/d3.rs

sci_form/hf/
d3.rs
