oxiphysics_gpu/kernels/md_force/
coulombkernel_traits.rs1#[allow(unused_imports)]
12use super::functions::*;
13use crate::compute::ComputeKernel;
14
15use super::types::CoulombKernel;
16
17#[allow(clippy::needless_range_loop)]
18impl ComputeKernel for CoulombKernel {
19 fn name(&self) -> &str {
20 "CoulombKernel"
21 }
22 fn execute(&self, inputs: &[&[f64]], outputs: &mut [Vec<f64>], work_size: usize) {
23 if inputs.len() < 3 || outputs.len() < 2 {
24 return;
25 }
26 let pos = inputs[0];
27 let charges = inputs[1];
28 let k_e = inputs[2][0];
29 let cutoff = inputs[2][1];
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 let qi = charges[i];
37 for j in (i + 1)..n {
38 let xj = [pos[j * 3], pos[j * 3 + 1], pos[j * 3 + 2]];
39 let dx = xi[0] - xj[0];
40 let dy = xi[1] - xj[1];
41 let dz = xi[2] - xj[2];
42 let r2 = dx * dx + dy * dy + dz * dz;
43 if r2 >= cutoff2 || r2 < 1e-30 {
44 continue;
45 }
46 let r = r2.sqrt();
47 let qj = charges[j];
48 potential += k_e * qi * qj / r;
49 let f_mag = k_e * qi * qj / (r2 * r);
50 forces[i * 3] += f_mag * dx;
51 forces[i * 3 + 1] += f_mag * dy;
52 forces[i * 3 + 2] += f_mag * dz;
53 forces[j * 3] -= f_mag * dx;
54 forces[j * 3 + 1] -= f_mag * dy;
55 forces[j * 3 + 2] -= f_mag * dz;
56 }
57 }
58 outputs[0] = forces;
59 outputs[1] = vec![potential];
60 }
61}