ak_core/calculator/
pair_potential.rs

1use std::iter::zip;
2
3use crate::calculator::{Calculator, CalculatorError};
4use crate::geometry::naive_neighbor_list as build_neighborlist;
5
6pub trait PairPotential {
7    fn cutoff(&self) -> f64;
8    fn pair_energy(&self, r: f64) -> f64;
9    fn pair_force_magnitude(&self, r: f64) -> f64;
10    fn name(&self) -> &str;
11}
12
13impl<T: PairPotential> Calculator for T {
14    fn calculate_energy(&self, view: &crate::StructureView) -> Result<f64, CalculatorError> {
15        let nl = build_neighborlist(view, self.cutoff());
16        let distances = nl.distance.ok_or(CalculatorError::CalculationFailed)?;
17        let energy = distances.iter().map(|r| self.pair_energy(*r)).sum();
18        Ok(energy)
19    }
20
21    fn calculate_forces(
22        &self,
23        view: &crate::StructureView,
24    ) -> Result<Vec<[f64; 3]>, CalculatorError> {
25        let nl = build_neighborlist(view, self.cutoff());
26        let distances = nl.distance.ok_or(CalculatorError::CalculationFailed)?;
27
28        let mut forces: Vec<[f64; 3]> = vec![[0.0_f64; 3]; view.len()];
29
30        for ((i, j), r) in zip(zip(nl.i, nl.j), distances) {
31            let pos_i = view.positions[i];
32            let pos_j = view.positions[j];
33            let r_vec = [
34                pos_j[0] - pos_i[0],
35                pos_j[1] - pos_i[1],
36                pos_j[2] - pos_i[2],
37            ];
38            let f_mag = self.pair_force_magnitude(r);
39            let r_hat = [r_vec[0] / r, r_vec[1] / r, r_vec[2] / r]; // unit vector
40
41            forces[i] = [
42                forces[i][0] + f_mag * r_hat[0],
43                forces[i][1] + f_mag * r_hat[1],
44                forces[i][2] + f_mag * r_hat[2],
45            ];
46            forces[j] = [
47                forces[j][0] - f_mag * r_hat[0],
48                forces[j][1] - f_mag * r_hat[1],
49                forces[j][2] - f_mag * r_hat[2],
50            ];
51        }
52
53        Ok(forces)
54    }
55}
ak_core/calculator/pair_potential.rs

ak_core/calculator/
pair_potential.rs
