Skip to main content

sci_form/forcefield/etkdg_3d/
mod.rs

1//! ETKDG 3D force field matching RDKit's construct3DForceField.
2//! Uses flat-bottom distance constraints (zero gradient within tolerance)
3//! instead of harmonic bond stretch + angle bend.
4
5pub mod builder;
6pub mod energy;
7pub mod gradient;
8pub mod optimizer;
9
10pub use builder::*;
11pub use energy::*;
12pub use gradient::*;
13pub use optimizer::*;
14
15const KNOWN_DIST_TOL: f64 = 0.01;
16const KNOWN_DIST_K: f64 = 100.0;
17
18/// A pre-computed flat-bottom distance constraint
19#[derive(Clone)]
20pub struct DistConstraint {
21    pub i: usize,
22    pub j: usize,
23    pub min_len: f64,
24    pub max_len: f64,
25    pub k: f64,
26}
27
28/// Pre-computed M6 torsion contribution (Fourier series)
29#[derive(Clone)]
30pub struct M6TorsionContrib {
31    pub i: usize,
32    pub j: usize,
33    pub k: usize,
34    pub l: usize,
35    pub signs: [f64; 6],
36    pub v: [f64; 6],
37}
38
39/// Pre-computed UFF inversion contribution (out-of-plane)
40/// Matches RDKit's InversionContrib exactly:
41///   E = K * (C0 + C1*sinY + C2*cos2W)
42/// where Y is the Wilson out-of-plane angle.
43/// For SP2 C/N/O: C0=1, C1=-1, C2=0 → E = K*(1-sinY)
44#[derive(Clone)]
45pub struct UFFInversionContrib {
46    /// Neighbor atom 1 (idx1 in RDKit)
47    pub at1: usize,
48    /// Central atom (idx2 in RDKit)
49    pub at2: usize,
50    /// Neighbor atom 2 (idx3 in RDKit)
51    pub at3: usize,
52    /// Neighbor atom 3 (idx4 in RDKit)
53    pub at4: usize,
54    pub force_constant: f64,
55    pub c0: f64,
56    pub c1: f64,
57    pub c2: f64,
58}
59
60/// Flat-bottom angle constraint (degrees) for SP (linear) centers.
61/// E = k * angleTerm^2, where angleTerm = clamp(angle - [min,max], 0)
62#[derive(Clone)]
63pub struct AngleConstraint {
64    pub i: usize, // outer atom 1
65    pub j: usize, // central atom
66    pub k: usize, // outer atom 2
67    pub min_deg: f64,
68    pub max_deg: f64,
69    pub force_k: f64,
70}
71
72/// Pre-computed 3D ETKDG force field terms.
73/// Distance constraints are split into three groups matching RDKit's construct3DForceField
74/// ordering: 1-2 bonds, 1-3 angles, and long-range. This ensures identical
75/// floating-point accumulation order in energy/gradient, which is critical
76/// for BFGS trajectory reproducibility.
77pub struct Etkdg3DFF {
78    pub dist_12: Vec<DistConstraint>,
79    pub dist_13: Vec<DistConstraint>,
80    pub dist_long: Vec<DistConstraint>,
81    pub angle_constraints: Vec<AngleConstraint>,
82    pub torsion_contribs: Vec<M6TorsionContrib>,
83    pub inversion_contribs: Vec<UFFInversionContrib>,
84    pub oop_k: f64,
85    pub torsion_k_omega: f64,
86    pub bounds_force_scaling: f64,
87    pub use_m6_torsions: bool,
88}