 #include <vector>
 
 using namespace cub;
+
+#define IMAGES_PER_BLOCK 64
 
 #define FULL_MASK 0xffffffff
 

  int numberOfAtomsInsideBox
)
{

  int particleID =  blockDim.x * blockIdx.x + threadIdx.x;
  if (particleID >= numberOfAtomsInsideBox) return;
  
  double forceX = 0.0, forceY = 0.0, forceZ = 0.0;



  // calculate block and grid sizes
  dim3 threadsPerBlock(256, 1, 1);
  int blocksPerGridX = (int)(numberOfAtomsInsideBox / threadsPerBlock) + 1;
  int blocksPerGridY = (int)(imageSize / IMAGES_PER_BLOCK) + 1;
  dim3 blocksPerGrid(blocksPerGridX, blocksPerGridY, 1);

  CUMALLOC((void **) &gpu_particleCharge, particleCharge.size() * sizeof(double));
  CUMALLOC((void **) &gpu_particleHasNoCharge, particleHasNoCharge.size() * sizeof(bool));
