dreid_forge/model/topology.rs
1//! Force field topology and potential function definitions.
2//!
3//! This module defines the data structures for DREIDING force field
4//! parameters and potential energy functions. After parameterization,
5//! a [`ForgedSystem`] contains all information needed to run molecular
6//! dynamics or energy minimization simulations.
7//!
8//! # Potential Functions
9//!
10//! The DREIDING force field supports several potential function types:
11//!
12//! - **Bonds**: [`BondPotential`] — Harmonic or Morse stretching
13//! - **Angles**: [`AnglePotential`] — Cosine-harmonic or theta-harmonic bending
14//! - **Dihedrals**: [`DihedralPotential`] — Periodic torsion potentials
15//! - **Impropers**: [`ImproperPotential`] — Out-of-plane (planar/umbrella)
16//! - **Van der Waals**: [`VdwPairPotential`] — LJ 12-6 or Exp-6
17//! - **Hydrogen bonds**: [`HBondPotential`] — Directional H-bond terms
18//!
19//! # Output Structure
20//!
21//! The [`ForgedSystem`] struct combines the original molecular system
22//! with computed atom types, partial charges, and all bonded/non-bonded
23//! potential parameters needed for simulation.
24
25use super::system::System;
26
27/// Per-atom force field parameters.
28///
29/// Contains the computed properties for a single atom after
30/// DREIDING parameterization.
31///
32/// # Fields
33///
34/// * `charge` — Partial atomic charge in elementary charge units
35/// * `mass` — Atomic mass in atomic mass units (amu)
36/// * `type_idx` — Index into the atom type table
37#[derive(Debug, Clone, PartialEq)]
38pub struct AtomParam {
39 /// Partial atomic charge (e).
40 pub charge: f64,
41 /// Atomic mass (amu).
42 pub mass: f64,
43 /// Index into the atom type name table.
44 pub type_idx: usize,
45}
46
47/// Bond stretching potential functions.
48#[derive(Debug, Clone, PartialEq)]
49pub enum BondPotential {
50 /// Harmonic bond stretching potential.
51 Harmonic {
52 /// First atom index.
53 i: usize,
54 /// Second atom index.
55 j: usize,
56 /// Force constant (kcal/mol/Ų).
57 k_force: f64,
58 /// Equilibrium bond length (Å).
59 r0: f64,
60 },
61 /// Morse anharmonic bond potential.
62 Morse {
63 /// First atom index.
64 i: usize,
65 /// Second atom index.
66 j: usize,
67 /// Equilibrium bond length (Å).
68 r0: f64,
69 /// Bond dissociation energy (kcal/mol).
70 d0: f64,
71 /// Morse alpha parameter (Å⁻¹).
72 alpha: f64,
73 },
74}
75
76/// Angle bending potential functions.
77#[derive(Debug, Clone, PartialEq)]
78pub enum AnglePotential {
79 /// Cosine-harmonic angle potential (DREIDING default).
80 CosineHarmonic {
81 /// First atom index.
82 i: usize,
83 /// Central atom index.
84 j: usize,
85 /// Third atom index.
86 k: usize,
87 /// Force constant (kcal/mol/rad²).
88 k_force: f64,
89 /// Equilibrium angle (degrees).
90 theta0: f64,
91 },
92 /// Simple theta-harmonic angle potential.
93 ThetaHarmonic {
94 /// First atom index.
95 i: usize,
96 /// Central atom index.
97 j: usize,
98 /// Third atom index.
99 k: usize,
100 /// Force constant (kcal/mol/rad²).
101 k_force: f64,
102 /// Equilibrium angle (degrees).
103 theta0: f64,
104 },
105}
106
107/// Proper dihedral (torsion) potential.
108#[derive(Debug, Clone, PartialEq)]
109pub struct DihedralPotential {
110 /// First atom index.
111 pub i: usize,
112 /// Second atom index (bond axis).
113 pub j: usize,
114 /// Third atom index (bond axis).
115 pub k: usize,
116 /// Fourth atom index.
117 pub l: usize,
118 /// Barrier height (kcal/mol).
119 pub v_barrier: f64,
120 /// Periodicity (number of minima in 360°).
121 pub periodicity: i32,
122 /// Phase offset (degrees).
123 pub phase_offset: f64,
124}
125
126/// Improper dihedral (out-of-plane) potential functions.
127#[derive(Debug, Clone, PartialEq)]
128pub enum ImproperPotential {
129 /// Planar improper for sp² centers.
130 Planar {
131 /// First peripheral atom.
132 i: usize,
133 /// Central atom (sp² center).
134 j: usize,
135 /// Second peripheral atom.
136 k: usize,
137 /// Third peripheral atom.
138 l: usize,
139 /// Force constant (kcal/mol/rad²).
140 k_force: f64,
141 /// Equilibrium out-of-plane angle (degrees).
142 chi0: f64,
143 },
144 /// Umbrella improper for pyramidal centers.
145 Umbrella {
146 /// Central atom (pyramidal center).
147 center: usize,
148 /// First peripheral atom.
149 p1: usize,
150 /// Second peripheral atom.
151 p2: usize,
152 /// Third peripheral atom.
153 p3: usize,
154 /// Force constant (kcal/mol).
155 k_force: f64,
156 /// Equilibrium umbrella angle (degrees).
157 psi0: f64,
158 },
159}
160
161/// Van der Waals non-bonded pair potential functions.
162#[derive(Debug, Clone, PartialEq)]
163pub enum VdwPairPotential {
164 /// Lennard-Jones 12-6 potential.
165 LennardJones {
166 /// First atom type index.
167 type1_idx: usize,
168 /// Second atom type index.
169 type2_idx: usize,
170 /// LJ sigma parameter (Å).
171 sigma: f64,
172 /// LJ epsilon parameter (kcal/mol).
173 epsilon: f64,
174 },
175 /// Exponential-6 potential (Buckingham-like).
176 Exponential6 {
177 /// First atom type index.
178 type1_idx: usize,
179 /// Second atom type index.
180 type2_idx: usize,
181 /// Repulsive prefactor A.
182 a: f64,
183 /// Exponential decay parameter B (Å⁻¹).
184 b: f64,
185 /// Attractive coefficient C (kcal·Å⁶/mol).
186 c: f64,
187 },
188}
189
190/// Hydrogen bond directional potential.
191#[derive(Debug, Clone, PartialEq)]
192pub struct HBondPotential {
193 /// Donor atom type index (D in D-H···A).
194 pub donor_type_idx: usize,
195 /// Hydrogen atom type index (H in D-H···A).
196 pub hydrogen_type_idx: usize,
197 /// Acceptor atom type index (A in D-H···A).
198 pub acceptor_type_idx: usize,
199 /// H-bond equilibrium energy (kcal/mol).
200 pub d0: f64,
201 /// Equilibrium H···A distance (Å).
202 pub r0: f64,
203}
204
205/// Collection of all potential energy functions for a system.
206///
207/// Groups all bonded and non-bonded interaction parameters
208/// computed during DREIDING parameterization.
209#[derive(Debug, Clone, Default)]
210pub struct Potentials {
211 /// Bond stretching potentials.
212 pub bonds: Vec<BondPotential>,
213 /// Angle bending potentials.
214 pub angles: Vec<AnglePotential>,
215 /// Proper dihedral (torsion) potentials.
216 pub dihedrals: Vec<DihedralPotential>,
217 /// Improper dihedral (out-of-plane) potentials.
218 pub impropers: Vec<ImproperPotential>,
219 /// Van der Waals pair potentials between atom types.
220 pub vdw_pairs: Vec<VdwPairPotential>,
221 /// Hydrogen bond potentials.
222 pub h_bonds: Vec<HBondPotential>,
223}
224
225/// A fully parameterized molecular system.
226///
227/// Contains the original [`System`] along with computed DREIDING
228/// force field parameters including atom types, partial charges,
229/// and all potential energy function parameters.
230///
231/// This is the primary output of the [`forge`](crate::forge) function
232/// and contains everything needed to write simulation input files
233/// for molecular dynamics packages.
234///
235/// # Fields
236///
237/// * `system` — Original molecular structure
238/// * `atom_types` — DREIDING atom type names (e.g., "C_3", "O_2")
239/// * `atom_properties` — Per-atom charges, masses, and type indices
240/// * `potentials` — All bonded and non-bonded potential parameters
241#[derive(Debug, Clone)]
242pub struct ForgedSystem {
243 /// The original molecular system with atoms and bonds.
244 pub system: System,
245 /// DREIDING atom type names, indexed by type_idx.
246 pub atom_types: Vec<String>,
247 /// Per-atom force field parameters.
248 pub atom_properties: Vec<AtomParam>,
249 /// All potential energy function parameters.
250 pub potentials: Potentials,
251}
252
253#[cfg(test)]
254mod tests {
255 use super::*;
256 use crate::model::system::System;
257
258 #[test]
259 fn atom_param_fields_and_clone() {
260 let p = AtomParam {
261 charge: -0.34,
262 mass: 12.011,
263 type_idx: 2,
264 };
265 assert!(p.charge < 0.0);
266 assert_eq!(p.mass, 12.011);
267 assert_eq!(p.type_idx, 2);
268 let q = p.clone();
269 assert_eq!(p, q);
270 }
271
272 #[test]
273 fn bond_potential_variants_and_debug() {
274 let h = BondPotential::Harmonic {
275 i: 0,
276 j: 1,
277 k_force: 300.0,
278 r0: 1.23,
279 };
280 let m = BondPotential::Morse {
281 i: 1,
282 j: 2,
283 r0: 1.5,
284 d0: 4.0,
285 alpha: 2.0,
286 };
287 match h {
288 BondPotential::Harmonic { i, j, k_force, r0 } => {
289 assert_eq!(i, 0);
290 assert_eq!(j, 1);
291 assert!(k_force > 0.0);
292 assert!(r0 > 0.0);
293 }
294 _ => panic!("expected Harmonic variant"),
295 }
296 match m {
297 BondPotential::Morse {
298 i,
299 j,
300 r0,
301 d0,
302 alpha,
303 } => {
304 assert_eq!(i, 1);
305 assert_eq!(j, 2);
306 assert!(d0 > 0.0);
307 assert!(alpha > 0.0);
308 assert!(r0 > 0.0);
309 }
310 _ => panic!("expected Morse variant"),
311 }
312 let s = format!("{:?} {:?}", h, m);
313 assert!(s.contains("Harmonic"));
314 assert!(s.contains("Morse"));
315 }
316
317 #[test]
318 fn angle_potential_variants() {
319 let a1 = AnglePotential::CosineHarmonic {
320 i: 0,
321 j: 1,
322 k: 2,
323 k_force: 50.0,
324 theta0: 109.5,
325 };
326 let a2 = AnglePotential::ThetaHarmonic {
327 i: 2,
328 j: 1,
329 k: 0,
330 k_force: 40.0,
331 theta0: 120.0,
332 };
333 match a1 {
334 AnglePotential::CosineHarmonic {
335 k_force, theta0, ..
336 } => {
337 assert_eq!(k_force, 50.0);
338 assert_eq!(theta0, 109.5);
339 }
340 _ => panic!("expected CosineHarmonic"),
341 }
342 match a2 {
343 AnglePotential::ThetaHarmonic {
344 k_force, theta0, ..
345 } => {
346 assert_eq!(k_force, 40.0);
347 assert_eq!(theta0, 120.0);
348 }
349 _ => panic!("expected ThetaHarmonic"),
350 }
351 }
352
353 #[test]
354 fn dihedral_and_improper_variants() {
355 let d = DihedralPotential {
356 i: 0,
357 j: 1,
358 k: 2,
359 l: 3,
360 v_barrier: 2.5,
361 periodicity: 3,
362 phase_offset: 180.0,
363 };
364 assert_eq!(d.periodicity, 3);
365 assert!(d.v_barrier > 0.0);
366
367 let imp1 = ImproperPotential::Planar {
368 i: 0,
369 j: 1,
370 k: 2,
371 l: 3,
372 k_force: 10.0,
373 chi0: 0.0,
374 };
375 let imp2 = ImproperPotential::Umbrella {
376 center: 1,
377 p1: 2,
378 p2: 3,
379 p3: 4,
380 k_force: 5.0,
381 psi0: 180.0,
382 };
383 match imp1 {
384 ImproperPotential::Planar { k_force, chi0, .. } => {
385 assert_eq!(k_force, 10.0);
386 assert_eq!(chi0, 0.0);
387 }
388 _ => panic!("expected Planar"),
389 }
390 match imp2 {
391 ImproperPotential::Umbrella {
392 center, p1, p2, p3, ..
393 } => {
394 assert_eq!(center, 1);
395 assert_eq!(p1, 2);
396 assert_eq!(p2, 3);
397 assert_eq!(p3, 4);
398 }
399 _ => panic!("expected Umbrella"),
400 }
401 }
402
403 #[test]
404 fn vdw_and_hbond_variants_and_potentials_container() {
405 let lj = VdwPairPotential::LennardJones {
406 type1_idx: 0,
407 type2_idx: 1,
408 sigma: 3.5,
409 epsilon: 0.2,
410 };
411 let ex6 = VdwPairPotential::Exponential6 {
412 type1_idx: 1,
413 type2_idx: 2,
414 a: 1000.0,
415 b: 50.0,
416 c: 2.0,
417 };
418 match lj {
419 VdwPairPotential::LennardJones { sigma, epsilon, .. } => {
420 assert_eq!(sigma, 3.5);
421 assert_eq!(epsilon, 0.2);
422 }
423 _ => panic!("expected LennardJones"),
424 }
425 match ex6 {
426 VdwPairPotential::Exponential6 { a, b, c, .. } => {
427 assert!(a > 0.0);
428 assert!(b > 0.0);
429 assert!(c > 0.0);
430 }
431 _ => panic!("expected Exponential6"),
432 }
433
434 let hb = HBondPotential {
435 donor_type_idx: 0,
436 hydrogen_type_idx: 1,
437 acceptor_type_idx: 2,
438 d0: 9.5,
439 r0: 2.75,
440 };
441 assert_eq!(hb.donor_type_idx, 0);
442 assert_eq!(hb.hydrogen_type_idx, 1);
443 assert_eq!(hb.acceptor_type_idx, 2);
444 assert_eq!(hb.d0, 9.5);
445 assert_eq!(hb.r0, 2.75);
446
447 let mut pots = Potentials::default();
448 pots.vdw_pairs.push(lj.clone());
449 pots.vdw_pairs.push(ex6.clone());
450 pots.h_bonds.push(hb.clone());
451 assert_eq!(pots.vdw_pairs.len(), 2);
452 assert_eq!(pots.h_bonds.len(), 1);
453 }
454
455 #[test]
456 fn forged_system_basic_fields_and_debug() {
457 let sys = System::new();
458 let atom_types = vec!["C_3".to_string(), "O_2".to_string()];
459 let atom_props = vec![
460 AtomParam {
461 charge: -0.1,
462 mass: 12.0,
463 type_idx: 0,
464 },
465 AtomParam {
466 charge: -0.2,
467 mass: 16.0,
468 type_idx: 1,
469 },
470 ];
471 let pots = Potentials::default();
472 let fs = ForgedSystem {
473 system: sys.clone(),
474 atom_types: atom_types.clone(),
475 atom_properties: atom_props.clone(),
476 potentials: pots,
477 };
478 assert_eq!(fs.atom_types.len(), 2);
479 assert_eq!(fs.atom_properties.len(), 2);
480 let s = format!("{:?}", fs);
481 assert!(s.contains("ForgedSystem"));
482 assert!(s.contains("atom_types"));
483 assert!(s.contains("atom_properties"));
484 }
485}