oxiphysics_gpu/kernels/md_force/
pairenergyaccumulatekernel_traits.rs1use crate::compute::ComputeKernel;
11
12use super::types::PairEnergyAccumulateKernel;
13
14impl ComputeKernel for PairEnergyAccumulateKernel {
15 fn name(&self) -> &str {
16 "PairEnergyAccumulateKernel"
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 pairs = inputs[1];
24 let epsilon = inputs[2][0];
25 let sigma = inputs[2][1];
26 let cutoff = inputs[2][2];
27 let n = work_size;
28 let cutoff2 = cutoff * cutoff;
29 let mut per_particle_energy = vec![0.0f64; n];
30 let num_pairs = pairs.len() / 2;
31 for p in 0..num_pairs {
32 let i = pairs[p * 2] as usize;
33 let j = pairs[p * 2 + 1] as usize;
34 if i >= n || j >= n {
35 continue;
36 }
37 let dx = pos[i * 3] - pos[j * 3];
38 let dy = pos[i * 3 + 1] - pos[j * 3 + 1];
39 let dz = pos[i * 3 + 2] - pos[j * 3 + 2];
40 let r2 = dx * dx + dy * dy + dz * dz;
41 if r2 >= cutoff2 || r2 < 1e-30 {
42 continue;
43 }
44 let r2_inv = 1.0 / r2;
45 let sr2 = sigma * sigma * r2_inv;
46 let sr6 = sr2 * sr2 * sr2;
47 let sr12 = sr6 * sr6;
48 let energy = 4.0 * epsilon * (sr12 - sr6);
49 per_particle_energy[i] += energy * 0.5;
50 per_particle_energy[j] += energy * 0.5;
51 }
52 let total: f64 = per_particle_energy.iter().sum();
53 outputs[0] = per_particle_energy;
54 outputs[1] = vec![total];
55 }
56}