vsepr_rs/
lib.rs

1//! vsepr-rs: A lightweight, high-performance molecular geometry optimizer.
2//!
3//! This crate provides a generic engine to refine 3D molecular coordinates 
4//! based on VSEPR (Valence Shell Electron Pair Repulsion) theory.
5//! 
6//! It is designed as a **scaffolder** or **pre-optimizer**: it quickly transforms
7//! raw or overlapping coordinates into a chemically sensible 3D structure that can then
8//! be passed to more rigorous force fields like UFF (Universal Force Field).
9//!
10//! # Example
11//! ```
12//! use vsepr_rs::{VseprOptimizer, AtomTrait, BondTrait};
13//!
14//! struct MyAtom { pos: [f64; 3], element: usize }
15//! impl AtomTrait for MyAtom {
16//!     fn get_position(&self) -> [f64; 3] { self.pos }
17//!     fn set_position(&mut self, pos: [f64; 3]) { self.pos = pos; }
18//!     fn atomic_number(&self) -> usize { self.element }
19//! }
20//!
21//! struct MyBond { pair: (usize, usize) }
22//! impl BondTrait for MyBond {
23//!     fn get_atom_indices(&self) -> (usize, usize) { self.pair }
24//!     fn get_bond_order(&self) -> f32 { 1.0 }
25//! }
26//!
27//! let mut atoms = vec![
28//!     MyAtom { pos: [0.0, 0.0, 0.0], element: 6 }, // C
29//!     MyAtom { pos: [0.1, 0.0, 0.0], element: 1 }, // H
30//!     MyAtom { pos: [0.0, 0.1, 0.0], element: 1 }, // H
31//!     MyAtom { pos: [0.0, 0.0, 0.1], element: 1 }, // H
32//!     MyAtom { pos: [-0.1, 0.0, 0.0], element: 1 }, // H
33//! ];
34//! let bonds = vec![
35//!     MyBond { pair: (0, 1) }, MyBond { pair: (0, 2) },
36//!     MyBond { pair: (0, 3) }, MyBond { pair: (0, 4) },
37//! ];
38//!
39//! let optimizer = VseprOptimizer::default();
40//! optimizer.optimize(&mut atoms, &bonds);
41//! ```
42
43pub mod forcefield;
44pub mod math;
45pub mod optimizer;
46pub mod traits;
47pub mod vsepr;
48
49pub use math::Vec3;
50pub use optimizer::VseprOptimizer;
51pub use traits::{get_covalent_radius, AtomTrait, BondTrait};
52pub use vsepr::{calculate_steric_number, Geometry};
53
54#[cfg(test)]
55mod tests {
56    use super::*;
57
58    struct MyAtom {
59        pos: [f64; 3],
60        element: usize,
61    }
62
63    impl AtomTrait for MyAtom {
64        fn get_position(&self) -> [f64; 3] { self.pos }
65        fn set_position(&mut self, pos: [f64; 3]) { self.pos = pos; }
66        fn atomic_number(&self) -> usize { self.element }
67    }
68
69    struct MyBond {
70        pair: (usize, usize),
71        order: f32,
72    }
73
74    impl BondTrait for MyBond {
75        fn get_atom_indices(&self) -> (usize, usize) { self.pair }
76        fn get_bond_order(&self) -> f32 { self.order }
77    }
78
79    #[test]
80    fn test_methane_tetrahedral() {
81        let mut atoms = vec![
82            MyAtom { pos: [0.0, 0.0, 0.0], element: 6 }, // C
83            MyAtom { pos: [0.0, 0.0, 0.0], element: 1 }, // H
84            MyAtom { pos: [0.0, 0.0, 0.0], element: 1 }, // H
85            MyAtom { pos: [0.0, 0.0, 0.0], element: 1 }, // H
86            MyAtom { pos: [0.0, 0.0, 0.0], element: 1 }, // H
87        ];
88        let bonds = vec![
89            MyBond { pair: (0, 1), order: 1.0 },
90            MyBond { pair: (0, 2), order: 1.0 },
91            MyBond { pair: (0, 3), order: 1.0 },
92            MyBond { pair: (0, 4), order: 1.0 },
93        ];
94
95        let optimizer = VseprOptimizer::default();
96        optimizer.optimize(&mut atoms, &bonds);
97
98        // Check C-H distance (~1.09 A)
99        for i in 1..5 {
100            let dist = Vec3::from(atoms[0].pos).dist(Vec3::from(atoms[i].pos));
101            assert!((dist - 1.09).abs() < 0.1, "C-H bond length too far from ideal: {}", dist);
102        }
103
104        // Check H-C-H angle (ideal ~109.5 deg)
105        let v1 = Vec3::from(atoms[1].pos).sub(Vec3::from(atoms[0].pos));
106        let v2 = Vec3::from(atoms[2].pos).sub(Vec3::from(atoms[0].pos));
107        let angle = v1.angle(v2).to_degrees();
108        assert!((angle - 109.5).abs() < 5.0, "H-C-H angle too far from ideal: {}", angle);
109    }
110
111    #[test]
112    fn test_water_bent() {
113        let mut atoms = vec![
114            MyAtom { pos: [0.0, 0.0, 0.0], element: 8 }, // O
115            MyAtom { pos: [0.0, 0.0, 0.0], element: 1 }, // H
116            MyAtom { pos: [0.0, 0.0, 0.0], element: 1 }, // H
117        ];
118        let bonds = vec![
119            MyBond { pair: (0, 1), order: 1.0 },
120            MyBond { pair: (0, 2), order: 1.0 },
121        ];
122
123        let optimizer = VseprOptimizer::default();
124        optimizer.optimize(&mut atoms, &bonds);
125
126        // O-H distance (~0.96 A)
127        let dist = Vec3::from(atoms[0].pos).dist(Vec3::from(atoms[1].pos));
128        assert!((dist - 0.96).abs() < 0.1);
129
130        // H-O-H angle (VSEPR SN=4, ideal ~109.5, actual ~104.5)
131        // Our simplified model targets the geometry's ideal angle (109.5 for SN=4).
132        let v1 = Vec3::from(atoms[1].pos).sub(Vec3::from(atoms[0].pos));
133        let v2 = Vec3::from(atoms[2].pos).sub(Vec3::from(atoms[0].pos));
134        let angle = v1.angle(v2).to_degrees();
135        assert!((angle - 109.5).abs() < 5.0, "H-O-H angle check failed: {}", angle);
136    }
137
138    #[test]
139    fn test_benzene_from_origin() {
140        let mut atoms = Vec::new();
141        for _ in 0..6 {
142            atoms.push(MyAtom { pos: [0.0, 0.0, 0.0], element: 6 });
143        }
144        for _ in 0..6 {
145             atoms.push(MyAtom { pos: [0.0, 0.0, 0.0], element: 1 });
146        }
147        run_benzene_check(atoms, "Origin Initialization");
148    }
149
150    fn run_benzene_check(mut atoms: Vec<MyAtom>, test_name: &str) {
151        let mut bonds = Vec::new();
152        for i in 0..6 {
153            bonds.push(MyBond { pair: (i, (i + 1) % 6), order: 1.5 });
154        }
155        for i in 0..6 {
156            bonds.push(MyBond { pair: (i, 6 + i), order: 1.0 });
157        }
158
159        let optimizer = VseprOptimizer { iterations: 3000, force_constant: 0.1 };
160        optimizer.optimize(&mut atoms, &bonds);
161
162        println!("--- {} ---", test_name);
163        let p1 = Vec3::from(atoms[0].pos);
164        let p2 = Vec3::from(atoms[1].pos);
165        let dist = p1.dist(p2);
166        println!("C-C dist: {:.3}", dist);
167        assert!((dist - 1.40).abs() < 0.15, "Bond length check failed");
168    }
169}
vsepr_rs/lib.rs

vsepr_rs/
lib.rs
