sci_form/scf/
gradients.rs

1//! Analytical and numerical nuclear gradients ∂E/∂R.
2//!
3//! Numerical (central finite differences) gradients:
4//!   ∂E/∂R_i ≈ [E(R_i + δ) - E(R_i - δ)] / (2δ)
5//!
6//! Also provides the energy-weighted density matrix W and
7//! nuclear repulsion gradient needed for analytical gradients.
8
9use nalgebra::DMatrix;
10
11use super::scf_loop::{run_scf, ScfConfig};
12use super::types::MolecularSystem;
13
14/// Gradient result for all atoms.
15#[derive(Debug, Clone)]
16pub struct GradientResult {
17    /// Nuclear gradients dE/dR in Hartree/Bohr, shape [n_atoms][3].
18    pub gradients: Vec<[f64; 3]>,
19    /// RMS gradient magnitude.
20    pub rms_gradient: f64,
21    /// Maximum gradient component.
22    pub max_gradient: f64,
23    /// Total energy at the evaluated geometry.
24    pub energy: f64,
25}
26
27impl GradientResult {
28    pub fn norm(&self) -> f64 {
29        self.rms_gradient
30    }
31}
32
33/// Compute nuclear gradients using central finite differences.
34///
35/// Cost: 6N SCF calculations for N atoms.
36pub fn numerical_gradient(
37    system: &MolecularSystem,
38    config: &ScfConfig,
39    step_size: f64,
40) -> GradientResult {
41    let n_atoms = system.n_atoms();
42    let mut gradients = vec![[0.0; 3]; n_atoms];
43
44    let e_ref = run_scf(system, config).total_energy;
45
46    for atom in 0..n_atoms {
47        for coord in 0..3 {
48            let mut sys_plus = system.clone();
49            sys_plus.positions_bohr[atom][coord] += step_size;
50            let e_plus = run_scf(&sys_plus, config).total_energy;
51
52            let mut sys_minus = system.clone();
53            sys_minus.positions_bohr[atom][coord] -= step_size;
54            let e_minus = run_scf(&sys_minus, config).total_energy;
55
56            gradients[atom][coord] = (e_plus - e_minus) / (2.0 * step_size);
57        }
58    }
59
60    let mut sum_sq = 0.0;
61    let mut max_g = 0.0f64;
62    let n_components = (n_atoms * 3) as f64;
63
64    for g in &gradients {
65        for &gc in g {
66            sum_sq += gc * gc;
67            max_g = max_g.max(gc.abs());
68        }
69    }
70
71    GradientResult {
72        gradients,
73        rms_gradient: (sum_sq / n_components).sqrt(),
74        max_gradient: max_g,
75        energy: e_ref,
76    }
77}
78
79/// Compute the energy-weighted density matrix W.
80///
81///   W_μν = 2 Σ_i^occ ε_i C_μi C_νi
82pub fn energy_weighted_density(
83    coefficients: &DMatrix<f64>,
84    orbital_energies: &[f64],
85    n_occupied: usize,
86) -> DMatrix<f64> {
87    let n_basis = coefficients.nrows();
88    let mut w = DMatrix::zeros(n_basis, n_basis);
89
90    for k in 0..n_occupied {
91        let eps_k = orbital_energies[k];
92        for mu in 0..n_basis {
93            for nu in 0..n_basis {
94                w[(mu, nu)] += 2.0 * eps_k * coefficients[(mu, k)] * coefficients[(nu, k)];
95            }
96        }
97    }
98
99    w
100}
101
102/// Compute nuclear repulsion gradient for atom A.
103///
104///   ∂V_nn/∂R_A = Z_A Σ_{B≠A} Z_B (R_A - R_B) / |R_A - R_B|³
105pub fn nuclear_repulsion_gradient(
106    atomic_numbers: &[u8],
107    positions: &[[f64; 3]],
108    atom_a: usize,
109) -> [f64; 3] {
110    let z_a = atomic_numbers[atom_a] as f64;
111    let r_a = positions[atom_a];
112    let mut grad = [0.0; 3];
113
114    for (b, &z_b_u8) in atomic_numbers.iter().enumerate() {
115        if b == atom_a {
116            continue;
117        }
118        let z_b = z_b_u8 as f64;
119        let r_b = positions[b];
120
121        let dx = r_a[0] - r_b[0];
122        let dy = r_a[1] - r_b[1];
123        let dz = r_a[2] - r_b[2];
124        let r2 = dx * dx + dy * dy + dz * dz;
125        let r = r2.sqrt();
126        let r3 = r * r2;
127
128        let factor = z_a * z_b / r3;
129        grad[0] += factor * dx;
130        grad[1] += factor * dy;
131        grad[2] += factor * dz;
132    }
133
134    grad
135}
136
137#[cfg(test)]
138mod tests {
139    use super::*;
140
141    #[test]
142    fn test_nuclear_repulsion_gradient_h2() {
143        let z = vec![1u8, 1];
144        let d = 1.4;
145        let pos = vec![[0.0, 0.0, 0.0], [d, 0.0, 0.0]];
146
147        let g0 = nuclear_repulsion_gradient(&z, &pos, 0);
148        let g1 = nuclear_repulsion_gradient(&z, &pos, 1);
149
150        // Newton's 3rd law
151        for i in 0..3 {
152            assert!((g0[i] + g1[i]).abs() < 1e-12);
153        }
154
155        assert!(g0[0] < 0.0, "Atom 0 should be pushed in -x direction");
156    }
157
158    #[test]
159    fn test_energy_weighted_density_symmetry() {
160        let n = 3;
161        let mut c = DMatrix::zeros(n, n);
162        c[(0, 0)] = 0.7;
163        c[(1, 0)] = 0.5;
164        c[(2, 0)] = 0.3;
165
166        let energies = vec![-0.5, 0.3, 1.2];
167        let w = energy_weighted_density(&c, &energies, 1);
168
169        for i in 0..n {
170            for j in 0..n {
171                assert!((w[(i, j)] - w[(j, i)]).abs() < 1e-12);
172            }
173        }
174    }
175}
sci_form/scf/gradients.rs

sci_form/scf/
gradients.rs
