     if (p < device_settings.num_mpm_markers) {
         const real3 xi = sorted_pos[p];
         real particle_density = 0;
-        LOOP_TWO_RING_GPU(                                              //
-            real weight = N(xi - current_node_location, inv_bin_edge);  //
-            particle_density += grid_mass[current_node] * weight;       //
+        int cx, cy, cz;
+        const real bin_edge = device_settings.bin_edge;
+        const real inv_bin_edge = device_settings.inv_bin_edge;
+        LOOP_TWO_RING_GPUSP(  //
+            real weight = N((xi.x - current_node_locationx) * inv_bin_edge) *
+                          N((xi.y - current_node_locationy) * inv_bin_edge) *
+                          N((xi.z - current_node_locationz) * inv_bin_edge);
+
+            particle_density += grid_mass[current_node] * weight;  //
             )
         // Inverse density to remove division
         particle_density = (bin_edge * bin_edge * bin_edge) / particle_density;

    const int p = blockIdx.x * blockDim.x + threadIdx.x;
    if (p < device_settings.num_mpm_markers) {
        const real3 xi = sorted_pos[p];
        LOOP_TWO_RING_GPU(                                                                     //
            real weight = N(xi - current_node_location, inv_bin_edge) * device_settings.mass;  //
            ATOMIC_ADD(&grid_mass[current_node], weight);                                      //
            )
    }
}

#define a_max 1e13
#define neg_BB1_fallback 0.11
#define neg_BB2_fallback 0.12

#define LOOP_TWO_RING_GPUSP(X)                                                                                         \
    cx = GridCoord(xi.x, inv_bin_edge, system_bounds.minimum[0]);                                                      \
