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



  // calculate block and grid sizes
  int threadsPerBlock = 256;
  int blocksPerGrid = (int)(numberOfAtomsInsideBox / threadsPerBlock) + 1;

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


  double lambdaCoef = DeviceGetLambdaCoulomb(moleculeID, kindID, box, gpu_isFraction, gpu_molIndex, gpu_kindIndex, gpu_lambdaCoulomb);
  // loop over images
  for(int vectorIndex = offset_vector_index; vectorIndex < offset_vector_index + IMAGES_PER_BLOCK; vectorIndex ++) {
    double dotx = x * gpu_kx[vectorIndex];
    double doty = y * gpu_ky[vectorIndex];
    double dotz = z * gpu_kz[vectorIndex];
