oxiphysics_gpu/kernels/md_force/
ewaldrealspacekernel_traits.rs

1//! # EwaldRealSpaceKernel - Trait Implementations
2//!
3//! This module contains trait implementations for `EwaldRealSpaceKernel`.
4//!
5//! ## Implemented Traits
6//!
7//! - `ComputeKernel`
8//!
9//! 🤖 Generated with [SplitRS](https://github.com/cool-japan/splitrs)
10
11use crate::compute::ComputeKernel;
12
13use super::functions::erfc_approx;
14#[allow(unused_imports)]
15use super::functions::*;
16use super::types::EwaldRealSpaceKernel;
17
18impl ComputeKernel for EwaldRealSpaceKernel {
19    fn name(&self) -> &str {
20        "EwaldRealSpaceKernel"
21    }
22    #[allow(clippy::needless_range_loop)]
23    fn execute(&self, inputs: &[&[f64]], outputs: &mut [Vec<f64>], work_size: usize) {
24        if inputs.len() < 3 || outputs.len() < 2 {
25            return;
26        }
27        let pos = inputs[0];
28        let charges = inputs[1];
29        let alpha = inputs[2][0];
30        let r_cutoff = inputs[2][1];
31        let box_len = inputs[2][2];
32        let n = work_size;
33        let r_cut2 = r_cutoff * r_cutoff;
34        let half_box = box_len * 0.5;
35        let mut forces = vec![0.0f64; n * 3];
36        let mut energy = 0.0f64;
37        for i in 0..n {
38            let xi = [pos[i * 3], pos[i * 3 + 1], pos[i * 3 + 2]];
39            let qi = charges[i];
40            for j in (i + 1)..n {
41                let xj = [pos[j * 3], pos[j * 3 + 1], pos[j * 3 + 2]];
42                let qj = charges[j];
43                let mut dx = xi[0] - xj[0];
44                let mut dy = xi[1] - xj[1];
45                let mut dz = xi[2] - xj[2];
46                if dx > half_box {
47                    dx -= box_len;
48                } else if dx < -half_box {
49                    dx += box_len;
50                }
51                if dy > half_box {
52                    dy -= box_len;
53                } else if dy < -half_box {
54                    dy += box_len;
55                }
56                if dz > half_box {
57                    dz -= box_len;
58                } else if dz < -half_box {
59                    dz += box_len;
60                }
61                let r2 = dx * dx + dy * dy + dz * dz;
62                if r2 >= r_cut2 || r2 < 1e-30 {
63                    continue;
64                }
65                let r = r2.sqrt();
66                let ar = alpha * r;
67                let erfc_ar = erfc_approx(ar);
68                energy += qi * qj * erfc_ar / r;
69                let two_alpha_over_sqrt_pi = 2.0 * alpha / std::f64::consts::PI.sqrt();
70                let deriv = -(erfc_ar / r + two_alpha_over_sqrt_pi * (-ar * ar).exp()) / r;
71                let f_mag = -qi * qj * deriv / r;
72                forces[i * 3] += f_mag * dx;
73                forces[i * 3 + 1] += f_mag * dy;
74                forces[i * 3 + 2] += f_mag * dz;
75                forces[j * 3] -= f_mag * dx;
76                forces[j * 3 + 1] -= f_mag * dy;
77                forces[j * 3 + 2] -= f_mag * dz;
78            }
79        }
80        outputs[0] = forces;
81        outputs[1] = vec![energy];
82    }
83}
oxiphysics_gpu/kernels/md_force/ewaldrealspacekernel_traits.rs

oxiphysics_gpu/kernels/md_force/
ewaldrealspacekernel_traits.rs
