sci_form/forcefield/
mmff94.rs

1use super::traits::ForceFieldContribution;
2
3// ─── MMFF94 Atom Type Assignment ─────────────────────────────────────────────
4
5/// MMFF94 atom type (simplified subset for common organic elements).
6/// Full MMFF94 has 99 types; we cover the most frequent organic types.
7#[derive(Debug, Clone, Copy, PartialEq, Eq)]
8#[repr(u8)]
9pub enum Mmff94AtomType {
10    CR = 1,    // Alkyl carbon, SP3
11    CSp2 = 2,  // Vinyl carbon, SP2
12    CSp = 3,   // Acetylenic carbon, SP
13    CO = 4,    // Carbonyl carbon C=O
14    CR4R = 20, // Carbon in 4-membered ring
15    CR3R = 22, // Carbon in 3-membered ring
16    CB = 37,   // Aromatic carbon
17    NR = 8,    // Amine nitrogen SP3
18    N2 = 9,    // Imine nitrogen SP2
19    NC = 10,   // Isonitrile nitrogen SP
20    NAm = 40,  // Amide nitrogen
21    NR2 = 39,  // Aromatic nitrogen (pyridine)
22    OR = 6,    // Alcohol/ether oxygen
23    O2 = 7,    // Carbonyl oxygen
24    OX = 32,   // Carboxylate oxygen
25    F = 11,    // Fluorine
26    Cl = 12,   // Chlorine
27    Br = 13,   // Bromine
28    I = 14,    // Iodine
29    S = 15,    // Thiol sulfur
30    SO = 17,   // Sulfoxide S=O
31    SO2 = 18,  // Sulfone S(=O)2
32    P = 25,    // Phosphorus SP3
33    HC = 5,    // Hydrogen on carbon
34    HO = 21,   // Hydrogen on oxygen
35    HN = 23,   // Hydrogen on nitrogen
36    HS = 71,   // Hydrogen on sulfur
37    Unknown = 0,
38}
39
40/// Assign MMFF94 atom type from element number and hybridization.
41pub fn assign_mmff94_type(
42    element: u8,
43    hyb: &crate::graph::Hybridization,
44    is_aromatic: bool,
45    is_amide_n: bool,
46) -> Mmff94AtomType {
47    use crate::graph::Hybridization::*;
48    match element {
49        1 => Mmff94AtomType::HC,   // Default; caller should refine based on neighbor
50        5 => Mmff94AtomType::CSp2, // Approximate: boron as sp2 carbon
51        6 => {
52            if is_aromatic {
53                Mmff94AtomType::CB
54            } else {
55                match hyb {
56                    SP => Mmff94AtomType::CSp,
57                    SP2 => Mmff94AtomType::CSp2,
58                    _ => Mmff94AtomType::CR,
59                }
60            }
61        }
62        7 => {
63            if is_aromatic {
64                Mmff94AtomType::NR2
65            } else if is_amide_n {
66                Mmff94AtomType::NAm
67            } else {
68                match hyb {
69                    SP => Mmff94AtomType::NC,
70                    SP2 => Mmff94AtomType::N2,
71                    _ => Mmff94AtomType::NR,
72                }
73            }
74        }
75        8 => match hyb {
76            SP2 => Mmff94AtomType::O2,
77            _ => Mmff94AtomType::OR,
78        },
79        9 => Mmff94AtomType::F,
80        15 => Mmff94AtomType::P,
81        16 => Mmff94AtomType::S,
82        17 => Mmff94AtomType::Cl,
83        35 => Mmff94AtomType::Br,
84        53 => Mmff94AtomType::I,
85        _ => Mmff94AtomType::Unknown,
86    }
87}
88
89// ─── MMFF94 Bond Stretching ──────────────────────────────────────────────────
90
91/// MMFF94 bond stretching (quartic form):
92/// E = 0.5 · k_b · Δr² · (1 + cs · Δr + 7/12 · cs² · Δr²)
93/// where cs = -2.0 Å⁻¹ (cubic stretch constant)
94pub struct Mmff94BondStretch {
95    pub atom_i: usize,
96    pub atom_j: usize,
97    pub k_b: f64, // Force constant (md/Å)
98    pub r0: f64,  // Equilibrium bond length (Å)
99}
100
101const MMFF94_CUBIC_STRETCH: f64 = -2.0;
102
103impl ForceFieldContribution for Mmff94BondStretch {
104    fn evaluate_energy_and_inject_gradient(&self, coords: &[f64], grad: &mut [f64]) -> f64 {
105        let ri = self.atom_i * 3;
106        let rj = self.atom_j * 3;
107        let dx = coords[ri] - coords[rj];
108        let dy = coords[ri + 1] - coords[rj + 1];
109        let dz = coords[ri + 2] - coords[rj + 2];
110        let dist = (dx * dx + dy * dy + dz * dz).sqrt().max(1e-8);
111        let dr = dist - self.r0;
112        let cs = MMFF94_CUBIC_STRETCH;
113        let cs2 = cs * cs;
114
115        // Energy: 143.9325 * 0.5 * kb * dr^2 * (1 + cs*dr + 7/12*cs^2*dr^2)
116        let energy =
117            143.9325 * 0.5 * self.k_b * dr * dr * (1.0 + cs * dr + (7.0 / 12.0) * cs2 * dr * dr);
118
119        // Gradient: dE/dr
120        let de_dr = 143.9325 * self.k_b * dr * (1.0 + 1.5 * cs * dr + (7.0 / 6.0) * cs2 * dr * dr);
121        let scale = de_dr / dist;
122        grad[ri] += scale * dx;
123        grad[ri + 1] += scale * dy;
124        grad[ri + 2] += scale * dz;
125        grad[rj] -= scale * dx;
126        grad[rj + 1] -= scale * dy;
127        grad[rj + 2] -= scale * dz;
128
129        energy
130    }
131}
132
133// ─── MMFF94 Angle Bending ────────────────────────────────────────────────────
134
135/// MMFF94 angle bending:
136/// E = 0.5 · k_a · (Δθ)² · (1 + cb · Δθ)
137/// where cb = -0.014 deg⁻¹ (cubic bend constant)
138pub struct Mmff94AngleBend {
139    pub atom_i: usize,
140    pub atom_j: usize, // central
141    pub atom_k: usize,
142    pub k_a: f64,    // Force constant (md·Å/rad²)
143    pub theta0: f64, // Equilibrium angle (radians)
144}
145
146const MMFF94_CUBIC_BEND: f64 = -0.014;
147
148impl ForceFieldContribution for Mmff94AngleBend {
149    fn evaluate_energy_and_inject_gradient(&self, coords: &[f64], grad: &mut [f64]) -> f64 {
150        let ri = self.atom_i * 3;
151        let rj = self.atom_j * 3;
152        let rk = self.atom_k * 3;
153
154        let rji = [
155            coords[ri] - coords[rj],
156            coords[ri + 1] - coords[rj + 1],
157            coords[ri + 2] - coords[rj + 2],
158        ];
159        let rjk = [
160            coords[rk] - coords[rj],
161            coords[rk + 1] - coords[rj + 1],
162            coords[rk + 2] - coords[rj + 2],
163        ];
164
165        let d_ji = (rji[0] * rji[0] + rji[1] * rji[1] + rji[2] * rji[2])
166            .sqrt()
167            .max(1e-8);
168        let d_jk = (rjk[0] * rjk[0] + rjk[1] * rjk[1] + rjk[2] * rjk[2])
169            .sqrt()
170            .max(1e-8);
171
172        let cos_theta = (rji[0] * rjk[0] + rji[1] * rjk[1] + rji[2] * rjk[2]) / (d_ji * d_jk);
173        let cos_theta_clamped = cos_theta.clamp(-1.0, 1.0);
174        let theta = cos_theta_clamped.acos();
175        let dt = (theta - self.theta0) * 180.0 / std::f64::consts::PI; // In degrees for MMFF94
176
177        let cb = MMFF94_CUBIC_BEND;
178        let energy = 0.043844 * 0.5 * self.k_a * dt * dt * (1.0 + cb * dt);
179
180        // Gradient (simplified: project along angle bisector normal)
181        let de_dtheta =
182            0.043844 * self.k_a * dt * (1.0 + 1.5 * cb * dt) * (180.0 / std::f64::consts::PI);
183        let sin_theta = (1.0 - cos_theta_clamped * cos_theta_clamped)
184            .sqrt()
185            .max(1e-8);
186        let dcos = -1.0 / sin_theta;
187        let pref = de_dtheta * dcos;
188
189        for dim in 0..3 {
190            let term_i = pref * (rjk[dim] / (d_ji * d_jk) - cos_theta * rji[dim] / (d_ji * d_ji))
191                / d_ji
192                * d_ji;
193            let term_k = pref * (rji[dim] / (d_ji * d_jk) - cos_theta * rjk[dim] / (d_jk * d_jk))
194                / d_jk
195                * d_jk;
196            grad[ri + dim] += term_i;
197            grad[rk + dim] += term_k;
198            grad[rj + dim] -= term_i + term_k;
199        }
200
201        energy
202    }
203}
204
205// ─── MMFF94 Torsion ──────────────────────────────────────────────────────────
206
207/// MMFF94 torsion:
208/// E = 0.5 · (V1·(1+cos φ) + V2·(1-cos 2φ) + V3·(1+cos 3φ))
209pub struct Mmff94Torsion {
210    pub atom_i: usize,
211    pub atom_j: usize,
212    pub atom_k: usize,
213    pub atom_l: usize,
214    pub v1: f64,
215    pub v2: f64,
216    pub v3: f64,
217}
218
219impl ForceFieldContribution for Mmff94Torsion {
220    fn evaluate_energy_and_inject_gradient(&self, coords: &[f64], grad: &mut [f64]) -> f64 {
221        let p = |idx: usize| -> [f64; 3] {
222            [coords[idx * 3], coords[idx * 3 + 1], coords[idx * 3 + 2]]
223        };
224        let pi = p(self.atom_i);
225        let pj = p(self.atom_j);
226        let pk = p(self.atom_k);
227        let pl = p(self.atom_l);
228
229        // Compute dihedral using standard atan2 method
230        let b1 = [pj[0] - pi[0], pj[1] - pi[1], pj[2] - pi[2]];
231        let b2 = [pk[0] - pj[0], pk[1] - pj[1], pk[2] - pj[2]];
232        let b3 = [pl[0] - pk[0], pl[1] - pk[1], pl[2] - pk[2]];
233
234        let cross = |a: [f64; 3], b: [f64; 3]| -> [f64; 3] {
235            [
236                a[1] * b[2] - a[2] * b[1],
237                a[2] * b[0] - a[0] * b[2],
238                a[0] * b[1] - a[1] * b[0],
239            ]
240        };
241        let dot = |a: [f64; 3], b: [f64; 3]| -> f64 { a[0] * b[0] + a[1] * b[1] + a[2] * b[2] };
242
243        let n1 = cross(b1, b2);
244        let n2 = cross(b2, b3);
245        let m1 = cross(n1, b2);
246
247        let b2_len = dot(b2, b2).sqrt().max(1e-8);
248        let x = dot(n1, n2);
249        let y = dot(m1, n2) / b2_len;
250        let phi = (-y).atan2(x);
251
252        let energy = 0.5
253            * (self.v1 * (1.0 + phi.cos())
254                + self.v2 * (1.0 - (2.0 * phi).cos())
255                + self.v3 * (1.0 + (3.0 * phi).cos()));
256
257        // Numerical gradient for torsion (analytical is complex; use central differences)
258        let eps = 1e-5;
259        for atom_idx in [self.atom_i, self.atom_j, self.atom_k, self.atom_l] {
260            for dim in 0..3 {
261                let idx = atom_idx * 3 + dim;
262                let orig = coords[idx];
263                let mut c_plus = coords.to_vec();
264                let mut c_minus = coords.to_vec();
265                c_plus[idx] = orig + eps;
266                c_minus[idx] = orig - eps;
267
268                let phi_p = {
269                    let pi = [
270                        c_plus[self.atom_i * 3],
271                        c_plus[self.atom_i * 3 + 1],
272                        c_plus[self.atom_i * 3 + 2],
273                    ];
274                    let pj = [
275                        c_plus[self.atom_j * 3],
276                        c_plus[self.atom_j * 3 + 1],
277                        c_plus[self.atom_j * 3 + 2],
278                    ];
279                    let pk = [
280                        c_plus[self.atom_k * 3],
281                        c_plus[self.atom_k * 3 + 1],
282                        c_plus[self.atom_k * 3 + 2],
283                    ];
284                    let pl = [
285                        c_plus[self.atom_l * 3],
286                        c_plus[self.atom_l * 3 + 1],
287                        c_plus[self.atom_l * 3 + 2],
288                    ];
289                    let b1 = [pj[0] - pi[0], pj[1] - pi[1], pj[2] - pi[2]];
290                    let b2 = [pk[0] - pj[0], pk[1] - pj[1], pk[2] - pj[2]];
291                    let b3 = [pl[0] - pk[0], pl[1] - pk[1], pl[2] - pk[2]];
292                    let nn1 = cross(b1, b2);
293                    let nn2 = cross(b2, b3);
294                    let mm1 = cross(nn1, b2);
295                    let b2l = dot(b2, b2).sqrt().max(1e-8);
296                    (-dot(mm1, nn2) / b2l).atan2(dot(nn1, nn2))
297                };
298                let phi_m = {
299                    let pi = [
300                        c_minus[self.atom_i * 3],
301                        c_minus[self.atom_i * 3 + 1],
302                        c_minus[self.atom_i * 3 + 2],
303                    ];
304                    let pj = [
305                        c_minus[self.atom_j * 3],
306                        c_minus[self.atom_j * 3 + 1],
307                        c_minus[self.atom_j * 3 + 2],
308                    ];
309                    let pk = [
310                        c_minus[self.atom_k * 3],
311                        c_minus[self.atom_k * 3 + 1],
312                        c_minus[self.atom_k * 3 + 2],
313                    ];
314                    let pl = [
315                        c_minus[self.atom_l * 3],
316                        c_minus[self.atom_l * 3 + 1],
317                        c_minus[self.atom_l * 3 + 2],
318                    ];
319                    let b1 = [pj[0] - pi[0], pj[1] - pi[1], pj[2] - pi[2]];
320                    let b2 = [pk[0] - pj[0], pk[1] - pj[1], pk[2] - pj[2]];
321                    let b3 = [pl[0] - pk[0], pl[1] - pk[1], pl[2] - pk[2]];
322                    let nn1 = cross(b1, b2);
323                    let nn2 = cross(b2, b3);
324                    let mm1 = cross(nn1, b2);
325                    let b2l = dot(b2, b2).sqrt().max(1e-8);
326                    (-dot(mm1, nn2) / b2l).atan2(dot(nn1, nn2))
327                };
328
329                let e_p = 0.5
330                    * (self.v1 * (1.0 + phi_p.cos())
331                        + self.v2 * (1.0 - (2.0 * phi_p).cos())
332                        + self.v3 * (1.0 + (3.0 * phi_p).cos()));
333                let e_m = 0.5
334                    * (self.v1 * (1.0 + phi_m.cos())
335                        + self.v2 * (1.0 - (2.0 * phi_m).cos())
336                        + self.v3 * (1.0 + (3.0 * phi_m).cos()));
337                grad[idx] += (e_p - e_m) / (2.0 * eps);
338            }
339        }
340
341        energy
342    }
343}
344
345// ─── MMFF94 Buffered 14-7 Van der Waals ──────────────────────────────────────
346
347/// Dispersión estérica repulsiva/atractiva regida por el Potencial Amortiguado 14-7 (Buffered 14-7) de Halgren.
348pub struct Mmff94BufferedVanDerWaals {
349    pub atom_i_idx: usize,
350    pub atom_j_idx: usize,
351    pub radius_star: f64,   // Parámetro dimensional empírico cruzado R*ij
352    pub epsilon_depth: f64, // Factor de profundidad termodinámica eps_ij
353}
354
355impl ForceFieldContribution for Mmff94BufferedVanDerWaals {
356    fn evaluate_energy_and_inject_gradient(&self, coords: &[f64], grad: &mut [f64]) -> f64 {
357        let root_i = self.atom_i_idx * 3;
358        let root_j = self.atom_j_idx * 3;
359
360        let delta_x = coords[root_i] - coords[root_j];
361        let delta_y = coords[root_i + 1] - coords[root_j + 1];
362        let delta_z = coords[root_i + 2] - coords[root_j + 2];
363
364        let dist_squared = delta_x * delta_x + delta_y * delta_y + delta_z * delta_z;
365        let mut dist_r = dist_squared.sqrt();
366
367        // Tope asintótico absoluto inferior para colisiones
368        if dist_r < 1e-8 {
369            dist_r = 1e-8;
370        }
371
372        // Algebra fraccionaria amortiguadora: E_vdW = eps * (1.07 R* / (R + 0.07 R*))^7 * ((1.12 R*^7 / (R^7 + 0.12 R*^7)) - 2)
373        let r_star_powered_7 = self.radius_star.powi(7);
374        let dist_r_powered_7 = dist_r.powi(7);
375
376        let repulsive_denominator = dist_r + 0.07 * self.radius_star;
377        let repulsive_term = (1.07 * self.radius_star / repulsive_denominator).powi(7);
378
379        let attractive_denominator = dist_r_powered_7 + 0.12 * r_star_powered_7;
380        let attractive_term = (1.12 * r_star_powered_7 / attractive_denominator) - 2.0;
381
382        let vdw_total_energy = self.epsilon_depth * repulsive_term * attractive_term;
383
384        // Derivación espacial analítica
385        let gradient_rep_term = -7.0 * repulsive_term / repulsive_denominator;
386        let gradient_attr_term = -7.0 * dist_r.powi(6) * (1.12 * r_star_powered_7)
387            / (attractive_denominator * attractive_denominator);
388
389        let force_scalar_magnitude = self.epsilon_depth
390            * (gradient_rep_term * attractive_term + repulsive_term * gradient_attr_term);
391
392        // Factorización cartesiana
393        let vector_prefactor = force_scalar_magnitude / dist_r;
394        let grad_x = vector_prefactor * delta_x;
395        let grad_y = vector_prefactor * delta_y;
396        let grad_z = vector_prefactor * delta_z;
397
398        grad[root_i] += grad_x;
399        grad[root_i + 1] += grad_y;
400        grad[root_i + 2] += grad_z;
401
402        grad[root_j] -= grad_x;
403        grad[root_j + 1] -= grad_y;
404        grad[root_j + 2] -= grad_z;
405
406        vdw_total_energy
407    }
408}
409
410// ─── Gradient Validation ─────────────────────────────────────────────────────
411
412/// Validate analytical gradients against numerical (central-difference) gradients.
413/// Returns max absolute error across all coordinates.
414pub fn validate_gradients(term: &dyn ForceFieldContribution, coords: &[f64], eps: f64) -> f64 {
415    let n = coords.len();
416    let mut analytical_grad = vec![0.0; n];
417    term.evaluate_energy_and_inject_gradient(coords, &mut analytical_grad);
418
419    let mut max_err = 0.0f64;
420    for i in 0..n {
421        let mut c_plus = coords.to_vec();
422        let mut c_minus = coords.to_vec();
423        c_plus[i] += eps;
424        c_minus[i] -= eps;
425
426        let mut g_dummy = vec![0.0; n];
427        let e_plus = term.evaluate_energy_and_inject_gradient(&c_plus, &mut g_dummy);
428        g_dummy.fill(0.0);
429        let e_minus = term.evaluate_energy_and_inject_gradient(&c_minus, &mut g_dummy);
430
431        let numerical = (e_plus - e_minus) / (2.0 * eps);
432        let err = (analytical_grad[i] - numerical).abs();
433        max_err = max_err.max(err);
434    }
435    max_err
436}
437
438// ─── MMFF94 Builder (assemble terms for a molecule) ──────────────────────────
439
440/// Simple bond/angle/torsion lookup parameters for building MMFF94 terms.
441/// Uses fallback empirical rules when proper MMFF94 parameter tables are not available.
442pub struct Mmff94Builder;
443
444impl Mmff94Builder {
445    /// Estimate MMFF94 bond stretching parameters from elements and bond order.
446    fn bond_params(elem_i: u8, elem_j: u8, _bond_order: u8) -> (f64, f64) {
447        // r0 (Å) from covalent radii sum, kb from empirical rule
448        let r_cov = |e: u8| -> f64 {
449            match e {
450                1 => 0.31,
451                6 => 0.76,
452                7 => 0.71,
453                8 => 0.66,
454                9 => 0.57,
455                15 => 1.07,
456                16 => 1.05,
457                17 => 1.02,
458                35 => 1.20,
459                53 => 1.39,
460                _ => 1.0,
461            }
462        };
463        let r0 = r_cov(elem_i) + r_cov(elem_j);
464        let kb = 5.0; // Fallback; real MMFF94 uses typed parameters
465        (kb, r0)
466    }
467
468    /// Build all MMFF94 force field terms for a parsed molecule.
469    ///
470    /// `elements`: atomic numbers.
471    /// `bonds`: list of (atom_i, atom_j, bond_order).
472    /// `coords`: flat xyz coordinates.
473    ///
474    /// Returns boxed force field contributions.
475    pub fn build(
476        elements: &[u8],
477        bonds: &[(usize, usize, u8)],
478    ) -> Vec<Box<dyn ForceFieldContribution>> {
479        let n_atoms = elements.len();
480        let mut terms: Vec<Box<dyn ForceFieldContribution>> = Vec::new();
481
482        // Bond stretching terms
483        for &(i, j, order) in bonds {
484            let (kb, r0) = Self::bond_params(elements[i], elements[j], order);
485            terms.push(Box::new(Mmff94BondStretch {
486                atom_i: i,
487                atom_j: j,
488                k_b: kb,
489                r0,
490            }));
491        }
492
493        // Angle bending: find all i-j-k where (i,j) and (j,k) are bonded
494        let mut neighbors: Vec<Vec<usize>> = vec![Vec::new(); n_atoms];
495        for &(i, j, _) in bonds {
496            neighbors[i].push(j);
497            neighbors[j].push(i);
498        }
499        for j in 0..n_atoms {
500            let nbrs = &neighbors[j];
501            for a in 0..nbrs.len() {
502                for b in (a + 1)..nbrs.len() {
503                    let i = nbrs[a];
504                    let k = nbrs[b];
505                    terms.push(Box::new(Mmff94AngleBend {
506                        atom_i: i,
507                        atom_j: j,
508                        atom_k: k,
509                        k_a: 0.5,                       // Fallback force constant
510                        theta0: 109.5_f64.to_radians(), // SP3 default
511                    }));
512                }
513            }
514        }
515
516        // Torsion terms: find all i-j-k-l where (i,j), (j,k), (k,l) are bonded
517        for &(j, k, _) in bonds {
518            for &i in &neighbors[j] {
519                if i == k {
520                    continue;
521                }
522                for &l in &neighbors[k] {
523                    if l == j || l == i {
524                        continue;
525                    }
526                    terms.push(Box::new(Mmff94Torsion {
527                        atom_i: i,
528                        atom_j: j,
529                        atom_k: k,
530                        atom_l: l,
531                        v1: 0.0,
532                        v2: 1.0,
533                        v3: 0.0, // Generic 2-fold barrier
534                    }));
535                }
536            }
537        }
538
539        // 1-4 vdW terms (atoms separated by 3 bonds)
540        for &(j, k, _) in bonds {
541            for &i in &neighbors[j] {
542                if i == k {
543                    continue;
544                }
545                for &l in &neighbors[k] {
546                    if l == j || l == i {
547                        continue;
548                    }
549                    let r_star = 3.5; // Generic
550                    let eps = 0.05;
551                    terms.push(Box::new(Mmff94BufferedVanDerWaals {
552                        atom_i_idx: i,
553                        atom_j_idx: l,
554                        radius_star: r_star,
555                        epsilon_depth: eps,
556                    }));
557                }
558            }
559        }
560
561        terms
562    }
563
564    /// Compute total energy from all MMFF94 terms.
565    pub fn total_energy(
566        terms: &[Box<dyn ForceFieldContribution>],
567        coords: &[f64],
568    ) -> (f64, Vec<f64>) {
569        let n = coords.len();
570        let mut grad = vec![0.0; n];
571        let mut total = 0.0;
572        for term in terms {
573            total += term.evaluate_energy_and_inject_gradient(coords, &mut grad);
574        }
575        (total, grad)
576    }
577}
578
579#[cfg(test)]
580mod tests {
581    use super::*;
582
583    #[test]
584    fn test_mmff94_vdw_energy() {
585        let term = Mmff94BufferedVanDerWaals {
586            atom_i_idx: 0,
587            atom_j_idx: 1,
588            radius_star: 3.6,
589            epsilon_depth: 0.1,
590        };
591        let coords = vec![0.0, 0.0, 0.0, 3.6, 0.0, 0.0];
592        let mut grad = vec![0.0; 6];
593        let e = term.evaluate_energy_and_inject_gradient(&coords, &mut grad);
594        // At equilibrium distance, energy should be near -epsilon
595        assert!(e < 0.0 && e > -0.2, "vdW energy at R*: {e}");
596    }
597
598    #[test]
599    fn test_mmff94_vdw_gradient_validation() {
600        let term = Mmff94BufferedVanDerWaals {
601            atom_i_idx: 0,
602            atom_j_idx: 1,
603            radius_star: 3.6,
604            epsilon_depth: 0.1,
605        };
606        let coords = vec![0.0, 0.0, 0.0, 4.0, 0.0, 0.0];
607        let max_err = validate_gradients(&term, &coords, 1e-5);
608        assert!(max_err < 1e-4, "vdW gradient error: {max_err}");
609    }
610
611    #[test]
612    fn test_mmff94_bond_stretch() {
613        let term = Mmff94BondStretch {
614            atom_i: 0,
615            atom_j: 1,
616            k_b: 5.0,
617            r0: 1.5,
618        };
619        // At equilibrium: energy ≈ 0
620        let coords_eq = vec![0.0, 0.0, 0.0, 1.5, 0.0, 0.0];
621        let mut grad = vec![0.0; 6];
622        let e_eq = term.evaluate_energy_and_inject_gradient(&coords_eq, &mut grad);
623        assert!(e_eq.abs() < 1e-10, "bond stretch at r0: {e_eq}");
624
625        // Stretched: energy > 0
626        let coords_stretch = vec![0.0, 0.0, 0.0, 2.0, 0.0, 0.0];
627        grad.fill(0.0);
628        let e_str = term.evaluate_energy_and_inject_gradient(&coords_stretch, &mut grad);
629        assert!(
630            e_str > 0.0,
631            "bond stretch energy should be positive: {e_str}"
632        );
633    }
634
635    #[test]
636    fn test_mmff94_bond_stretch_gradient_validation() {
637        let term = Mmff94BondStretch {
638            atom_i: 0,
639            atom_j: 1,
640            k_b: 5.0,
641            r0: 1.5,
642        };
643        let coords = vec![0.0, 0.0, 0.0, 2.0, 0.1, 0.0];
644        let max_err = validate_gradients(&term, &coords, 1e-5);
645        assert!(max_err < 1e-3, "bond stretch gradient error: {max_err}");
646    }
647
648    #[test]
649    fn test_mmff94_torsion_energy() {
650        let term = Mmff94Torsion {
651            atom_i: 0,
652            atom_j: 1,
653            atom_k: 2,
654            atom_l: 3,
655            v1: 0.0,
656            v2: 5.0,
657            v3: 0.0,
658        };
659        // Planar trans: phi ≈ 180°
660        let coords = vec![-1.5, 1.0, 0.0, 0.0, 0.0, 0.0, 1.5, 0.0, 0.0, 3.0, 1.0, 0.0];
661        let mut grad = vec![0.0; 12];
662        let e = term.evaluate_energy_and_inject_gradient(&coords, &mut grad);
663        assert!(e.is_finite(), "torsion energy should be finite: {e}");
664    }
665
666    #[test]
667    fn test_mmff94_atom_typing() {
668        use crate::graph::Hybridization;
669        let t = assign_mmff94_type(6, &Hybridization::SP3, false, false);
670        assert_eq!(t, Mmff94AtomType::CR);
671        let t = assign_mmff94_type(6, &Hybridization::SP2, true, false);
672        assert_eq!(t, Mmff94AtomType::CB);
673        let t = assign_mmff94_type(7, &Hybridization::SP3, false, true);
674        assert_eq!(t, Mmff94AtomType::NAm);
675    }
676
677    #[test]
678    fn test_mmff94_builder_ethane() {
679        // Ethane: C-C with 6 hydrogens
680        let elements = vec![6, 6, 1, 1, 1, 1, 1, 1]; // C, C, H×6
681        let bonds = vec![
682            (0, 1, 1), // C-C
683            (0, 2, 1),
684            (0, 3, 1),
685            (0, 4, 1), // C-H
686            (1, 5, 1),
687            (1, 6, 1),
688            (1, 7, 1), // C-H
689        ];
690        // Staggered ethane coordinates (approximate)
691        let coords = vec![
692            0.0, 0.0, 0.0, // C0
693            1.54, 0.0, 0.0, // C1
694            -0.5, 0.9, 0.0, // H
695            -0.5, -0.9, 0.0, // H
696            -0.5, 0.0, 0.9, // H
697            2.04, 0.9, 0.0, // H
698            2.04, -0.9, 0.0, // H
699            2.04, 0.0, 0.9, // H
700        ];
701        let terms = Mmff94Builder::build(&elements, &bonds);
702        assert!(!terms.is_empty(), "should produce force field terms");
703
704        let (energy, grad) = Mmff94Builder::total_energy(&terms, &coords);
705        assert!(
706            energy.is_finite(),
707            "total energy should be finite: {energy}"
708        );
709        assert!(
710            grad.iter().all(|g| g.is_finite()),
711            "all gradients should be finite"
712        );
713    }
714
715    #[test]
716    fn test_mmff94_builder_gradient_consistency() {
717        // Verify total gradient is consistent with numerical for a simple system
718        let elements = vec![6, 6];
719        let bonds = vec![(0, 1, 1)];
720        let coords = vec![0.0, 0.0, 0.0, 1.6, 0.1, 0.0];
721        let terms = Mmff94Builder::build(&elements, &bonds);
722        let (_, analytical_grad) = Mmff94Builder::total_energy(&terms, &coords);
723
724        let eps = 1e-5;
725        for i in 0..coords.len() {
726            let mut cp = coords.clone();
727            let mut cm = coords.clone();
728            cp[i] += eps;
729            cm[i] -= eps;
730            let (ep, _) = Mmff94Builder::total_energy(&terms, &cp);
731            let (em, _) = Mmff94Builder::total_energy(&terms, &cm);
732            let numerical = (ep - em) / (2.0 * eps);
733            let err = (analytical_grad[i] - numerical).abs();
734            assert!(
735                err < 0.1,
736                "gradient mismatch at coord {i}: anal={:.6} num={:.6} err={:.6}",
737                analytical_grad[i],
738                numerical,
739                err
740            );
741        }
742    }
743}
sci_form/forcefield/mmff94.rs

sci_form/forcefield/
mmff94.rs
