Skip to main content

oxiphysics_gpu/kernels/md_force/
coulombkernel_traits.rs

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