ak_core/calculator/
lennard_jones.rs

1use crate::calculator::pair_potential::PairPotential;
2
3pub struct LennardJones {
4    epsilon: f64,
5    sigma: f64,
6    cutoff: f64,
7}
8
9impl LennardJones {
10    pub fn new(epsilon: f64, sigma: f64, cutoff: f64) -> Self {
11        Self {
12            epsilon,
13            sigma,
14            cutoff,
15        }
16    }
17}
18
19impl PairPotential for LennardJones {
20    fn pair_energy(&self, r: f64) -> f64 {
21        let x6 = (self.sigma / r).powi(6);
22        let x12 = x6.powi(2);
23        4.0 * self.epsilon * (x12 - x6)
24    }
25
26    fn pair_force_magnitude(&self, r: f64) -> f64 {
27        let x6 = (self.sigma / r).powi(6);
28        let x12 = x6.powi(2);
29        24.0 * self.epsilon * (2.0 * x12 - x6) / r
30    }
31
32    fn cutoff(&self) -> f64 {
33        self.cutoff
34    }
35
36    fn name(&self) -> &str {
37        "Lennard Jones"
38    }
39}
40
41#[cfg(test)]
42mod test {
43    use crate::Structure;
44    use crate::calculator::{Calculator, LennardJones, PairPotential};
45
46    fn test_structure_2atoms(sigma: f64) -> Structure {
47        let rmin: f64 = 2.0_f64.powf(1.0 / 6.0) * sigma;
48
49        let positions = [[0.0, 0.0, 0.0], [rmin, 0.0, 0.0]].to_vec();
50        let numbers = [1, 1].to_vec();
51        let cell = [[10.0, 0.0, 0.0], [0.0, 10.0, 0.0], [0.0, 0.0, 10.0]];
52        let pbc = [false, false, false];
53        Structure::new(positions.clone(), numbers.clone(), cell, pbc)
54    }
55
56    fn test_structure_3atoms(sigma: f64) -> Structure {
57        let rmin: f64 = 2.0_f64.powf(1.0 / 6.0) * sigma;
58
59        let positions = [[0.0, 0.0, 0.0], [rmin, 0.0, 0.0], [2.0 * rmin, 0.0, 0.0]].to_vec();
60        let numbers = [1, 1, 1].to_vec();
61        let cell = [[10.0, 0.0, 0.0], [0.0, 10.0, 0.0], [0.0, 0.0, 10.0]];
62        let pbc = [false, false, false];
63        Structure::new(positions.clone(), numbers.clone(), cell, pbc)
64    }
65
66    #[test]
67    fn pair_energy_test() {
68        let lj = LennardJones::new(1.0, 1.0, 5.0);
69        let rmin: f64 = 2.0_f64.powf(1.0 / 6.0) * lj.sigma;
70
71        let pair_energy = lj.pair_energy(rmin);
72
73        assert_eq!(pair_energy, -1.0);
74    }
75
76    #[test]
77    fn calculate_test_two_atoms() {
78        let lj = LennardJones::new(1.5, 1.0, 5.0);
79        let structure = test_structure_2atoms(lj.sigma);
80        let result = lj.calculate(&structure.view()).unwrap();
81        assert_eq!(result.energy, Some(-lj.epsilon));
82        assert!(
83            result
84                .forces
85                .unwrap()
86                .iter()
87                .flatten()
88                .all(|&f| f.abs() < 1e-10)
89        );
90    }
91
92    #[test]
93    fn calculate_test_three_atoms() {
94        let lj = LennardJones::new(1.25, 1.0, 2.0);
95        let structure = test_structure_3atoms(lj.sigma);
96        let result = lj.calculate(&structure.view()).unwrap();
97        assert_eq!(result.energy, Some(-2.0 * lj.epsilon));
98        assert!(
99            result
100                .forces
101                .unwrap()
102                .iter()
103                .flatten()
104                .all(|&f| f.abs() < 1e-10)
105        );
106    }
107
108    #[test]
109    fn calculate_test_long_cutoff_three_atoms() {
110        let lj = LennardJones::new(1.25, 1.0, 5.0);
111        let structure = test_structure_3atoms(lj.sigma);
112        let result = lj.calculate(&structure.view()).unwrap();
113        assert!(result.energy.unwrap() < -2.0 * lj.epsilon)
114    }
115}
ak_core/calculator/lennard_jones.rs

ak_core/calculator/
lennard_jones.rs
