sci_form/forcefield/
params.rs

1/// UFF atom type parameters.
2/// Source: Rappé, A. K. et al. "UFF, a full periodic table force field for molecular mechanics
3/// and molecular dynamics simulations." J. Am. Chem. Soc. 114 (1992): 10024–10035. Table 1.
4#[derive(Debug, Clone)]
5pub struct UffAtomParams {
6    pub label: &'static str,
7    pub element: u8, // Atomic number
8    pub r1: f64,     // Bond radius (Angstroms) — used in bond length formula
9    pub theta0: f64, // Equilibrium valence angle (degrees, convert to radians in builder)
10    pub x1: f64,     // VDW distance parameter (Angstroms)
11    pub d1: f64,     // VDW well depth (kcal/mol)
12    pub zeta: f64,   // GMP electronegativity scale factor
13    pub z_star: f64, // Effective nuclear charge for bond force constant
14    pub chi: f64,    // GMP electronegativity (for bond length correction)
15    pub v_tors: f64, // Torsion barrier Vi (kcal/mol), used as sqrt(Vi*Vj)/2 (eq. 16)
16}
17
18/// Retrieve UFF parameters by atom type label.
19pub fn get_uff_params(label: &str) -> Option<UffAtomParams> {
20    match label {
21        // Hydrogen
22        "H_" => Some(UffAtomParams {
23            label: "H_",
24            element: 1,
25            r1: 0.354,
26            theta0: 180.0,
27            x1: 2.886,
28            d1: 0.044,
29            zeta: 0.712,
30            z_star: 0.712,
31            chi: 4.528,
32            v_tors: 0.0,
33        }),
34        // Carbon
35        "C_3" => Some(UffAtomParams {
36            label: "C_3",
37            element: 6,
38            r1: 0.757,
39            theta0: 109.47,
40            x1: 3.851,
41            d1: 0.105,
42            zeta: 1.912,
43            z_star: 1.912,
44            chi: 5.343,
45            v_tors: 2.119,
46        }),
47        "C_2" => Some(UffAtomParams {
48            label: "C_2",
49            element: 6,
50            r1: 0.732,
51            theta0: 120.0,
52            x1: 3.851,
53            d1: 0.105,
54            zeta: 1.912,
55            z_star: 1.912,
56            chi: 5.343,
57            v_tors: 2.119,
58        }),
59        "C_R" => Some(UffAtomParams {
60            label: "C_R",
61            element: 6,
62            r1: 0.729,
63            theta0: 120.0,
64            x1: 3.851,
65            d1: 0.105,
66            zeta: 1.912,
67            z_star: 1.912,
68            chi: 5.343,
69            v_tors: 2.119,
70        }),
71        "C_1" => Some(UffAtomParams {
72            label: "C_1",
73            element: 6,
74            r1: 0.706,
75            theta0: 180.0,
76            x1: 3.851,
77            d1: 0.105,
78            zeta: 1.912,
79            z_star: 1.912,
80            chi: 5.343,
81            v_tors: 2.119,
82        }),
83        // Nitrogen
84        "N_3" => Some(UffAtomParams {
85            label: "N_3",
86            element: 7,
87            r1: 0.700,
88            theta0: 106.7,
89            x1: 3.660,
90            d1: 0.069,
91            zeta: 2.544,
92            z_star: 2.544,
93            chi: 6.899,
94            v_tors: 0.450,
95        }),
96        "N_2" => Some(UffAtomParams {
97            label: "N_2",
98            element: 7,
99            r1: 0.685,
100            theta0: 111.2,
101            x1: 3.660,
102            d1: 0.069,
103            zeta: 2.544,
104            z_star: 2.544,
105            chi: 6.899,
106            v_tors: 0.450,
107        }),
108        "N_R" => Some(UffAtomParams {
109            label: "N_R",
110            element: 7,
111            r1: 0.699,
112            theta0: 120.0,
113            x1: 3.660,
114            d1: 0.069,
115            zeta: 2.544,
116            z_star: 2.544,
117            chi: 6.899,
118            v_tors: 0.450,
119        }),
120        "N_1" => Some(UffAtomParams {
121            label: "N_1",
122            element: 7,
123            r1: 0.656,
124            theta0: 180.0,
125            x1: 3.660,
126            d1: 0.069,
127            zeta: 2.544,
128            z_star: 2.544,
129            chi: 6.899,
130            v_tors: 0.450,
131        }),
132        // Oxygen
133        "O_3" => Some(UffAtomParams {
134            label: "O_3",
135            element: 8,
136            r1: 0.658,
137            theta0: 104.51,
138            x1: 3.500,
139            d1: 0.060,
140            zeta: 2.300,
141            z_star: 2.300,
142            chi: 8.741,
143            v_tors: 2.100,
144        }),
145        "O_2" => Some(UffAtomParams {
146            label: "O_2",
147            element: 8,
148            r1: 0.634,
149            theta0: 120.0,
150            x1: 3.500,
151            d1: 0.060,
152            zeta: 2.300,
153            z_star: 2.300,
154            chi: 8.741,
155            v_tors: 2.100,
156        }),
157        "O_R" => Some(UffAtomParams {
158            label: "O_R",
159            element: 8,
160            r1: 0.680,
161            theta0: 110.0,
162            x1: 3.500,
163            d1: 0.060,
164            zeta: 2.300,
165            z_star: 2.300,
166            chi: 8.741,
167            v_tors: 2.100,
168        }),
169        "O_1" => Some(UffAtomParams {
170            label: "O_1",
171            element: 8,
172            r1: 0.639,
173            theta0: 180.0,
174            x1: 3.500,
175            d1: 0.060,
176            zeta: 2.300,
177            z_star: 2.300,
178            chi: 8.741,
179            v_tors: 2.100,
180        }),
181        // Halogens
182        "F_" => Some(UffAtomParams {
183            label: "F_",
184            element: 9,
185            r1: 0.668,
186            theta0: 180.0,
187            x1: 3.364,
188            d1: 0.050,
189            zeta: 1.735,
190            z_star: 1.735,
191            chi: 10.874,
192            v_tors: 1.082,
193        }),
194        "Cl" => Some(UffAtomParams {
195            label: "Cl",
196            element: 17,
197            r1: 1.029,
198            theta0: 109.47,
199            x1: 3.947,
200            d1: 0.227,
201            zeta: 2.315,
202            z_star: 2.315,
203            chi: 8.564,
204            v_tors: 1.673,
205        }),
206        "Br" => Some(UffAtomParams {
207            label: "Br",
208            element: 35,
209            r1: 1.137,
210            theta0: 109.47,
211            x1: 4.189,
212            d1: 0.251,
213            zeta: 2.283,
214            z_star: 2.283,
215            chi: 7.790,
216            v_tors: 1.682,
217        }),
218        "I_" => Some(UffAtomParams {
219            label: "I_",
220            element: 53,
221            r1: 1.353,
222            theta0: 109.47,
223            x1: 4.500,
224            d1: 0.339,
225            zeta: 2.650,
226            z_star: 2.650,
227            chi: 6.822,
228            v_tors: 1.530,
229        }),
230        // Sulfur
231        "S_3" => Some(UffAtomParams {
232            label: "S_3",
233            element: 16,
234            r1: 1.020,
235            theta0: 92.1,
236            x1: 4.035,
237            d1: 0.274,
238            zeta: 2.703,
239            z_star: 2.703,
240            chi: 10.000,
241            v_tors: 1.948,
242        }),
243        "S_2" => Some(UffAtomParams {
244            label: "S_2",
245            element: 16,
246            r1: 0.987,
247            theta0: 120.0,
248            x1: 4.035,
249            d1: 0.274,
250            zeta: 2.703,
251            z_star: 2.703,
252            chi: 10.000,
253            v_tors: 1.948,
254        }),
255        "S_R" => Some(UffAtomParams {
256            label: "S_R",
257            element: 16,
258            r1: 0.990,
259            theta0: 120.0,
260            x1: 4.035,
261            d1: 0.274,
262            zeta: 2.703,
263            z_star: 2.703,
264            chi: 10.000,
265            v_tors: 1.948,
266        }),
267        // Phosphorus
268        "P_3" => Some(UffAtomParams {
269            label: "P_3",
270            element: 15,
271            r1: 0.930,
272            theta0: 93.8,
273            x1: 4.147,
274            d1: 0.305,
275            zeta: 2.862,
276            z_star: 2.862,
277            chi: 8.000,
278            v_tors: 2.421,
279        }),
280        // Boron
281        "B_2" => Some(UffAtomParams {
282            label: "B_2",
283            element: 5,
284            r1: 0.828,
285            theta0: 120.0,
286            x1: 3.637,
287            d1: 0.180,
288            zeta: 1.755,
289            z_star: 1.755,
290            chi: 4.961,
291            v_tors: 0.0,
292        }),
293        // Silicon
294        "Si3" => Some(UffAtomParams {
295            label: "Si3",
296            element: 14,
297            r1: 1.117,
298            theta0: 109.47,
299            x1: 4.295,
300            d1: 0.402,
301            zeta: 2.323,
302            z_star: 2.323,
303            chi: 4.168,
304            v_tors: 1.225,
305        }),
306        _ => None,
307    }
308}
309
310/// Compute the UFF equilibrium bond length between two atom types.
311///
312/// Rappé eq. 2–4:
313/// r_ij = r_i + r_j + r_BO - r_EN
314/// r_BO = -λ · (r_i + r_j) · ln(BO)          [λ = 0.1332]
315/// r_EN = r_i · r_j · (√χ_i - √χ_j)² / (χ_i·r_i + χ_j·r_j)
316pub fn get_uff_bond_length(
317    pi: &UffAtomParams,
318    pj: &UffAtomParams,
319    bond_order: &crate::graph::BondOrder,
320) -> f64 {
321    const LAMBDA: f64 = 0.1332;
322
323    let bo_value: f64 = match bond_order {
324        crate::graph::BondOrder::Single => 1.0,
325        crate::graph::BondOrder::Aromatic => 1.5,
326        crate::graph::BondOrder::Double => 2.0,
327        crate::graph::BondOrder::Triple => 3.0,
328        _ => 1.0,
329    };
330
331    let r_bo = if (bo_value - 1.0).abs() < 1e-10 {
332        0.0
333    } else {
334        -LAMBDA * (pi.r1 + pj.r1) * bo_value.ln()
335    };
336
337    let r_en = if (pi.chi - pj.chi).abs() > 1e-6 {
338        let sqrt_diff = pi.chi.sqrt() - pj.chi.sqrt();
339        pi.r1 * pj.r1 * sqrt_diff * sqrt_diff / (pi.chi * pi.r1 + pj.chi * pj.r1)
340    } else {
341        0.0
342    };
343
344    (pi.r1 + pj.r1 + r_bo - r_en).max(0.5)
345}
346
347/// UFF bond stretching force constant K_ij (Rappé eq. 6).
348/// E_bond = ½ K_ij (r - r₀)²    →   K_ij = 664.12 · Z*_i · Z*_j / r₀³
349pub fn get_uff_bond_force_constant(pi: &UffAtomParams, pj: &UffAtomParams, r0: f64) -> f64 {
350    664.12 * pi.z_star * pj.z_star / r0.powi(3)
351}
352
353/// UFF angle bending force constant K_ijk (Rappé eq. 13, simplified).
354/// K_ijk = 332.06 · Z*_i · Z*_k / (r_ij · r_jk · r_ik²)
355/// where r_ik² = r_ij² + r_jk² - 2·r_ij·r_jk·cos(θ₀)
356pub fn get_uff_angle_force_constant(
357    pi: &UffAtomParams,
358    pk: &UffAtomParams,
359    r_ij: f64,
360    r_jk: f64,
361    theta0_rad: f64,
362) -> f64 {
363    let r_ik_sq = r_ij * r_ij + r_jk * r_jk - 2.0 * r_ij * r_jk * theta0_rad.cos();
364    if r_ik_sq < 1e-10 || r_ij < 1e-10 || r_jk < 1e-10 {
365        return 0.0;
366    }
367    332.06 * pi.z_star * pk.z_star / (r_ij * r_jk * r_ik_sq)
368}
369
370/// UFF torsion barriers for bond j–k (Rappé Section 3D, eq. 16).
371/// Returns (V_barrier, periodicity_n, cos_phi0).
372/// V_barrier = ½ · √(V_i · V_j)  for sp3–sp3 bonds.
373pub fn get_uff_torsion_params(
374    hyb_j: &crate::graph::Hybridization,
375    hyb_k: &crate::graph::Hybridization,
376    vj: f64,
377    vk: f64,
378    bond_order: &crate::graph::BondOrder,
379) -> (f64, f64, f64) {
380    use crate::graph::BondOrder::*;
381    use crate::graph::Hybridization::*;
382
383    match (hyb_j, hyb_k, bond_order) {
384        // sp3–sp3: three-fold (staggered → anti preferred)
385        (SP3, SP3, _) => (0.5 * (vj * vk).sqrt(), 3.0, 1.0),
386        // sp3–sp2 mixed: six-fold (nearly flat)
387        (SP3, SP2, _) | (SP2, SP3, _) => (1.0, 6.0, 1.0),
388        // sp2–sp2 conjugated (aromatic or double bond): strong planarity
389        (SP2, SP2, Aromatic) | (SP2, SP2, Double) => (5.0, 2.0, -1.0),
390        // sp2–sp2 single (biaryl/amide): moderate planarity
391        (SP2, SP2, Single) => (10.0, 2.0, -1.0),
392        // sp (linear): no torsion
393        (SP, _, _) | (_, SP, _) => (0.0, 1.0, 1.0),
394        _ => (1.0, 3.0, 1.0),
395    }
396}
sci_form/forcefield/params.rs

sci_form/forcefield/
params.rs
