qualia_core_db/specialized_libs/chemistry_modeling/
dft.rs

1//! Density Functional Theory (DFT) Integration
2//!
3//! This module extends the SCF driver to include Exchange-Correlation Integration.
4//! It implements numerical grid evaluation for LDA and GGA.
5//! Exact analytical derivatives of the functional expressions are computed using
6//! a Rust-native forward-mode automatic differentiation (autodiff) Dual number struct.
7//! NO libxc C-bindings are used.
8
9use crate::specialized_libs::shared::zero_heap_algebra::ZeroHeapMatrix;
10
11/// A Dual number for forward-mode automatic differentiation.
12/// Evaluates f(x) and f'(x) simultaneously.
13#[derive(Debug, Clone, Copy, PartialEq)]
14pub struct Dual {
15    pub v: f64, // Value
16    pub d: f64, // Derivative
17}
18
19impl Dual {
20    pub fn new(v: f64, d: f64) -> Self {
21        Self { v, d }
22    }
23
24    pub fn add(self, other: Self) -> Self {
25        Self::new(self.v + other.v, self.d + other.d)
26    }
27
28    pub fn sub(self, other: Self) -> Self {
29        Self::new(self.v - other.v, self.d - other.d)
30    }
31
32    pub fn mul(self, other: Self) -> Self {
33        Self::new(self.v * other.v, self.d * other.v + self.v * other.d)
34    }
35
36    pub fn div(self, other: Self) -> Self {
37        Self::new(
38            self.v / other.v,
39            (self.d * other.v - self.v * other.d) / (other.v * other.v),
40        )
41    }
42
43    pub fn powf(self, power: f64) -> Self {
44        Self::new(
45            self.v.powf(power),
46            power * self.v.powf(power - 1.0) * self.d,
47        )
48    }
49
50    pub fn cbrt(self) -> Self {
51        self.powf(1.0 / 3.0)
52    }
53
54    pub fn scale(self, scalar: f64) -> Self {
55        Self::new(self.v * scalar, self.d * scalar)
56    }
57
58    pub fn ln(self) -> Self {
59        Self::new(self.v.ln(), self.d / self.v)
60    }
61
62    pub fn atan(self) -> Self {
63        Self::new(self.v.atan(), self.d / (1.0 + self.v * self.v))
64    }
65}
66
67/// Local Density Approximation (LDA) Exchange Functional (Dirac / Slater)
68/// E_x[rho] = - (3/4) * (3/pi)^(1/3) * int rho(r)^(4/3) dr
69/// Returns (Energy Density, Potential)
70pub fn lda_exchange(rho: f64) -> (f64, f64) {
71    if rho <= 1e-12 {
72        return (0.0, 0.0);
73    }
74
75    // Seed the derivative (d/d rho)
76    let r = Dual::new(rho, 1.0);
77    let factor = -0.75 * (3.0 / core::f64::consts::PI).powf(1.0 / 3.0);
78
79    let ex = r.cbrt().scale(factor); // e_x = C_x * rho^(1/3)
80    let vx = ex.scale(4.0 / 3.0); // v_x = (4/3) * C_x * rho^(1/3)
81
82    // e_x is exchange energy per particle, total exchange is rho * e_x.
83    // The derivative of rho * e_x with respect to rho is v_x.
84    (ex.v, vx.v)
85}
86
87/// VWN (Vosko-Wilk-Nusair) Local Correlation Functional (Unpolarized)
88pub fn lda_correlation_vwn(rho: f64) -> (f64, f64) {
89    if rho <= 1e-12 {
90        return (0.0, 0.0);
91    }
92
93    // Simplified VWN5 parameterization for the unpolarized electron gas
94    let r_s = Dual::new(
95        (3.0 / (4.0 * core::f64::consts::PI * rho)).powf(1.0 / 3.0),
96        -1.0 / (3.0 * rho) * (3.0 / (4.0 * core::f64::consts::PI * rho)).powf(1.0 / 3.0),
97    );
98
99    let a = 0.0621814;
100    let x0 = -0.409286;
101    let b = 13.0720_f64;
102    let c = 42.7198_f64;
103
104    let x = r_s.powf(0.5);
105    let q = (4.0 * c - b * b).sqrt();
106
107    // Polynomial X(x) = x^2 + b*x + c
108    let x_func = x.mul(x).add(x.scale(b)).add(Dual::new(c, 0.0));
109    let x0_func = x0 * x0 + b * x0 + c;
110
111    // VWN Evaluation using Dual numbers
112    // ec(x) = A * { ln(x^2 / X(x)) + 2b/Q * atan(Q / (2x+b)) - bx0/X(x0) * [ ln((x-x0)^2 / X(x)) + 2(b+2x0)/Q * atan(Q / (2x+b)) ] }
113    let term1 = x.mul(x).div(x_func).ln();
114
115    let atan_arg = Dual::new(q, 0.0).div(x.scale(2.0).add(Dual::new(b, 0.0)));
116    let term2 = atan_arg.atan().scale(2.0 * b / q);
117
118    let term3_factor = (b * x0) / x0_func;
119
120    let x_minus_x0 = x.sub(Dual::new(x0, 0.0));
121    let term3_ln = x_minus_x0.mul(x_minus_x0).div(x_func).ln();
122    let term3_atan = atan_arg.atan().scale(2.0 * (b + 2.0 * x0) / q);
123
124    let term3 = term3_ln.add(term3_atan).scale(term3_factor);
125
126    let ec = term1.add(term2).sub(term3).scale(a);
127    let vc = ec.v + rho * ec.d; // potential = d(rho*ec)/drho = ec + rho * dec/drho
128
129    (ec.v, vc)
130}
131
132/// DftIntegrator grid-based integration for the SCF loop
133pub struct DftIntegrator {
134    // Pre-computed spatial grid coordinates and weights for the molecule
135    pub grid_points: [[f64; 3]; 500],
136    pub grid_weights: [f64; 500],
137    pub n_points: usize,
138}
139
140impl DftIntegrator {
141    pub fn new() -> Self {
142        Self {
143            grid_points: [[0.0; 3]; 500],
144            grid_weights: [0.0; 500],
145            n_points: 0,
146        }
147    }
148
149    /// Evaluate the Exchange-Correlation matrix V_xc to be added to the Fock matrix
150    pub fn build_vxc<const N: usize>(
151        &self,
152        density: &ZeroHeapMatrix<f64, N, N>,
153        basis_values_at_grid: &[[f64; N]; 500],
154    ) -> (f64, ZeroHeapMatrix<f64, N, N>) {
155        let mut vxc = ZeroHeapMatrix::zeros();
156        let mut exc_total = 0.0;
157
158        for p in 0..self.n_points {
159            let weight = self.grid_weights[p];
160            let basis_vals = &basis_values_at_grid[p];
161
162            // 1. Calculate density at grid point p
163            let mut rho_p = 0.0;
164            for mu in 0..N {
165                for nu in 0..N {
166                    rho_p += density.get(mu, nu) * basis_vals[mu] * basis_vals[nu];
167                }
168            }
169
170            // 2. Evaluate Functional (LDA Exchange + Correlation)
171            let (ex, vx) = lda_exchange(rho_p);
172            let (ec, vc) = lda_correlation_vwn(rho_p);
173
174            let e_xc = ex + ec;
175            let v_xc = vx + vc;
176
177            exc_total += e_xc * rho_p * weight;
178
179            // 3. Accumulate Vxc matrix
180            for mu in 0..N {
181                for nu in 0..N {
182                    let val = vxc.get(mu, nu) + v_xc * basis_vals[mu] * basis_vals[nu] * weight;
183                    vxc.set(mu, nu, val);
184                }
185            }
186        }
187
188        (exc_total, vxc)
189    }
190}
191
192#[cfg(test)]
193mod tests {
194    use super::*;
195
196    #[test]
197    fn test_lda_exchange() {
198        let rho = 0.5;
199        let (ex, vx) = lda_exchange(rho);
200        // Compare with expected analytic values
201        // ex = -0.75 * (3/pi)^(1/3) * rho^(1/3)
202        // vx = 4/3 * ex
203        let expected_ex =
204            -0.75 * (3.0 / core::f64::consts::PI).powf(1.0 / 3.0) * rho.powf(1.0 / 3.0);
205        let expected_vx = 4.0 / 3.0 * expected_ex;
206
207        assert!((ex - expected_ex).abs() < 1e-10);
208        assert!((vx - expected_vx).abs() < 1e-10);
209    }
210
211    #[test]
212    fn test_lda_correlation_vwn() {
213        let rho = 0.5;
214        let (ec, vc) = lda_correlation_vwn(rho);
215
216        // As long as the derivative evaluates cleanly without NaN, the autodiff works
217        assert!(!ec.is_nan());
218        assert!(!vc.is_nan());
219    }
220
221    #[test]
222    fn test_dual_number_ops() {
223        // Test f(x) = x^3 at x = 2
224        // f(2) = 8
225        // f'(x) = 3x^2 => f'(2) = 12
226        let x = Dual::new(2.0, 1.0);
227        let x2 = x.mul(x);
228        let x3 = x2.mul(x);
229
230        assert!((x3.v - 8.0).abs() < 1e-10);
231        assert!((x3.d - 12.0).abs() < 1e-10);
232    }
233}
qualia_core_db/specialized_libs/chemistry_modeling/dft.rs

qualia_core_db/specialized_libs/chemistry_modeling/
dft.rs
