oxiphysics_gpu/kernels/md_force/
lennardjoneskernel_traits.rs1use crate::compute::ComputeKernel;
11
12use super::types::LennardJonesKernel;
13
14impl ComputeKernel for LennardJonesKernel {
15 fn name(&self) -> &str {
16 "LennardJonesKernel"
17 }
18 fn execute(&self, inputs: &[&[f64]], outputs: &mut [Vec<f64>], work_size: usize) {
19 if inputs.len() < 2 || outputs.len() < 2 {
20 return;
21 }
22 let pos = inputs[0];
23 let epsilon = inputs[1][0];
24 let sigma = inputs[1][1];
25 let cutoff = inputs[1][2];
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 for j in (i + 1)..n {
33 let xj = [pos[j * 3], pos[j * 3 + 1], pos[j * 3 + 2]];
34 let dx = xi[0] - xj[0];
35 let dy = xi[1] - xj[1];
36 let dz = xi[2] - xj[2];
37 let r2 = dx * dx + dy * dy + dz * dz;
38 if r2 >= cutoff2 || r2 < 1e-30 {
39 continue;
40 }
41 let r2_inv = 1.0 / r2;
42 let sr2 = sigma * sigma * r2_inv;
43 let sr6 = sr2 * sr2 * sr2;
44 let sr12 = sr6 * sr6;
45 potential += 4.0 * epsilon * (sr12 - sr6);
46 let f_mag = 24.0 * epsilon * (2.0 * sr12 - sr6) * r2_inv;
47 forces[i * 3] += f_mag * dx;
48 forces[i * 3 + 1] += f_mag * dy;
49 forces[i * 3 + 2] += f_mag * dz;
50 forces[j * 3] -= f_mag * dx;
51 forces[j * 3 + 1] -= f_mag * dy;
52 forces[j * 3 + 2] -= f_mag * dz;
53 }
54 }
55 outputs[0] = forces;
56 outputs[1] = vec![potential];
57 }
58}