gosh_model/
edip.rs

1// [[file:../models.note::178e12ff][178e12ff]]
2use super::*;
3
4use std::collections::{HashMap, HashSet};
5
6use ::edip::EdipParameters;
7use gchemol::neighbors::Neighborhood;
8use vecfx::*;
9// 178e12ff ends here
10
11// [[file:../models.note::6e669f3b][6e669f3b]]
12#[derive(Clone, Debug, Default)]
13pub struct Edip {
14    virial: f64,
15    // for create neighbors
16    nh: Neighborhood,
17}
18
19impl Edip {
20    fn update_nh(&mut self, mol: &Molecule) {
21        self.nh = Neighborhood::new();
22        // use atom index (0-based) for node index
23        self.nh.update(mol.positions().enumerate());
24        if let Some(lat) = mol.get_lattice() {
25            self.nh.set_lattice(lat.matrix().into());
26        }
27    }
28}
29
30impl ChemicalModel for Edip {
31    fn compute(&mut self, mol: &Molecule) -> Result<Computed> {
32        const search_radius: f64 = 4.0;
33
34        // only works for silicon
35        let not_silicon = mol.symbols().any(|x| x != "Si");
36        if not_silicon {
37            bail!("EDIP potential model only works for Silicon");
38        }
39
40        self.update_nh(mol);
41        let n = mol.natoms();
42        let positions = mol.positions().collect_vec();
43        let mut neighbors = vec![];
44        let mut distances = HashMap::new();
45        // FIXME: rewrite for periodic system
46        let lat = mol.get_lattice();
47        for i in 0..n {
48            let mut connected = HashSet::new();
49            for x in self.nh.neighbors(i, search_radius) {
50                // FIXME: avoid recompute pair distance in edip crate
51                let j = x.node;
52                let pi: Vector3f = positions[i].into();
53                let pj: Vector3f = positions[j].into();
54                let d = if let Some(image) = x.image {
55                    // translation periodic image
56                    let t = lat.unwrap().to_cart(image);
57                    pj + t - pi
58                } else {
59                    pj - pi
60                };
61                distances.insert((i, j), d.into());
62                connected.insert(j);
63            }
64
65            neighbors.push(connected);
66        }
67
68        let params = EdipParameters::silicon();
69        let mut forces = vec![[0.0; 3]; n];
70        let (energy, virial) = ::edip::compute_forces(&mut forces, &neighbors, &distances, &params);
71
72        let mut computed = Computed::default();
73        computed.set_energy(energy);
74        computed.set_forces(forces);
75
76        // FIXME: could be removed
77        self.virial = virial;
78
79        Ok(computed)
80    }
81}
82// 6e669f3b ends here
83
84// [[file:../models.note::28122508][28122508]]
85#[test]
86fn test_edip() -> Result<()> {
87    use gchemol::prelude::*;
88    use gchemol::Molecule;
89    use vecfx::*;
90
91    let f = "./tests/files/si10.xyz";
92    let mol = Molecule::from_file(f)?;
93
94    let mut model = Edip::default();
95    let computed = model.compute(&mol)?;
96    let f = computed.get_forces().unwrap();
97    let f_norm = f.as_flat().as_vector_slice().norm();
98    assert!(dbg!(f_norm) <= 0.005);
99    let energy = computed.get_energy().unwrap();
100    approx::assert_relative_eq!(energy, -39.57630354331939, epsilon = 1e-5);
101    approx::assert_relative_eq!(model.virial, 0.00224989660978845, epsilon = 1e-5);
102
103    let f = "./tests/files/si5.xyz";
104    let mol = Molecule::from_file(f)?;
105    let computed = model.compute(&mol)?;
106    let energy = computed.get_energy().unwrap();
107    approx::assert_relative_eq!(energy, -14.566606, epsilon = 1e-5);
108    // FIXME: refactor
109    let virial = model.virial;
110    approx::assert_relative_eq!(virial, -3.643552, epsilon = 1e-5);
111
112    let f = computed.get_forces().unwrap();
113    #[rustfmt::skip]
114    let f_expected = [ -0.19701000,  -0.62522600,   0.02948000,
115                       -0.23330600,   0.43698600,   0.42511400,
116                        0.43297600,   0.18333500,  -0.44864300,
117                       -1.75669300,   0.50149400,  -1.48879300,
118                        1.75403300,  -0.49658900,   1.48284200];
119    approx::assert_relative_eq!(
120        f.as_flat().as_vector_slice(),
121        f_expected.as_vector_slice(),
122        epsilon = 1E-5,
123    );
124
125    // for silicon crystal
126    let mol = Molecule::from_file("./tests/files/si-3x3x3.cif")?;
127    let computed = model.compute(&mol)?;
128    let energy = computed.get_energy().unwrap();
129    let f = computed.get_forces().unwrap();
130    dbg!(f.as_flat().as_vector_slice().norm());
131
132    Ok(())
133}
134// 28122508 ends here
gosh_model/edip.rs

gosh_model/
edip.rs
