oxiphysics_gpu/kernels/md_force/
bondforcekernel_traits.rs1#[allow(unused_imports)]
12use super::functions::*;
13use crate::compute::ComputeKernel;
14
15use super::types::BondForceKernel;
16
17impl ComputeKernel for BondForceKernel {
18 fn name(&self) -> &str {
19 "BondForceKernel"
20 }
21 fn execute(&self, inputs: &[&[f64]], outputs: &mut [Vec<f64>], work_size: usize) {
22 if inputs.len() < 2 || outputs.len() < 2 {
23 return;
24 }
25 let pos = inputs[0];
26 let bond_data = inputs[1];
27 let n = work_size;
28 let num_bonds = bond_data.len() / 4;
29 let mut forces = vec![0.0f64; n * 3];
30 let mut energies = vec![0.0f64; num_bonds];
31 for b in 0..num_bonds {
32 let i = bond_data[b * 4] as usize;
33 let j = bond_data[b * 4 + 1] as usize;
34 let k_spring = bond_data[b * 4 + 2];
35 let r0 = bond_data[b * 4 + 3];
36 if i >= n || j >= n {
37 continue;
38 }
39 let dx = pos[j * 3] - pos[i * 3];
40 let dy = pos[j * 3 + 1] - pos[i * 3 + 1];
41 let dz = pos[j * 3 + 2] - pos[i * 3 + 2];
42 let r = (dx * dx + dy * dy + dz * dz).sqrt();
43 if r < 1e-30 {
44 continue;
45 }
46 let delta = r - r0;
47 energies[b] = 0.5 * k_spring * delta * delta;
48 let mag = k_spring * delta / r;
49 forces[i * 3] += mag * dx;
50 forces[i * 3 + 1] += mag * dy;
51 forces[i * 3 + 2] += mag * dz;
52 forces[j * 3] -= mag * dx;
53 forces[j * 3 + 1] -= mag * dy;
54 forces[j * 3 + 2] -= mag * dz;
55 }
56 outputs[0] = forces;
57 outputs[1] = energies;
58 }
59}