Skip to main content

oxiphysics_gpu/kernels/md_force/
lennardjoneskernel_traits.rs

1//! # LennardJonesKernel - Trait Implementations
2//!
3//! This module contains trait implementations for `LennardJonesKernel`.
4//!
5//! ## Implemented Traits
6//!
7//! - `ComputeKernel`
8//!
9//! 🤖 Generated with [SplitRS](https://github.com/cool-japan/splitrs)
10
11#[allow(unused_imports)]
12use super::functions::*;
13use crate::compute::ComputeKernel;
14
15use super::types::LennardJonesKernel;
16
17#[allow(clippy::needless_range_loop)]
18impl ComputeKernel for LennardJonesKernel {
19    fn name(&self) -> &str {
20        "LennardJonesKernel"
21    }
22    fn execute(&self, inputs: &[&[f64]], outputs: &mut [Vec<f64>], work_size: usize) {
23        if inputs.len() < 2 || outputs.len() < 2 {
24            return;
25        }
26        let pos = inputs[0];
27        let epsilon = inputs[1][0];
28        let sigma = inputs[1][1];
29        let cutoff = inputs[1][2];
30        let n = work_size;
31        let cutoff2 = cutoff * cutoff;
32        let mut forces = vec![0.0; n * 3];
33        let mut potential = 0.0;
34        for i in 0..n {
35            let xi = [pos[i * 3], pos[i * 3 + 1], pos[i * 3 + 2]];
36            for j in (i + 1)..n {
37                let xj = [pos[j * 3], pos[j * 3 + 1], pos[j * 3 + 2]];
38                let dx = xi[0] - xj[0];
39                let dy = xi[1] - xj[1];
40                let dz = xi[2] - xj[2];
41                let r2 = dx * dx + dy * dy + dz * dz;
42                if r2 >= cutoff2 || r2 < 1e-30 {
43                    continue;
44                }
45                let r2_inv = 1.0 / r2;
46                let sr2 = sigma * sigma * r2_inv;
47                let sr6 = sr2 * sr2 * sr2;
48                let sr12 = sr6 * sr6;
49                potential += 4.0 * epsilon * (sr12 - sr6);
50                let f_mag = 24.0 * epsilon * (2.0 * sr12 - sr6) * r2_inv;
51                forces[i * 3] += f_mag * dx;
52                forces[i * 3 + 1] += f_mag * dy;
53                forces[i * 3 + 2] += f_mag * dz;
54                forces[j * 3] -= f_mag * dx;
55                forces[j * 3 + 1] -= f_mag * dy;
56                forces[j * 3 + 2] -= f_mag * dz;
57            }
58        }
59        outputs[0] = forces;
60        outputs[1] = vec![potential];
61    }
62}