 #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;

    forceZ += factor * gpu_kz[vectorIndex];
  }

  // loop over other particles within the same molecule
  if(blockIdx.y == 0) {
    double intraForce = 0.0, distSq = 0.0, dist = 0.0;
    double distVectX = 0.0, distVectY = 0.0, distVectZ = 0.0;
    int lastParticleWithinSameMolecule = gpu_startMol[particleID] + gpu_lengthMol[particleID];
    for(int otherParticle = gpu_startMol[particleID];
      otherParticle < lastParticleWithinSameMolecule;
      otherParticle++)
    {
      if(particleID != otherParticle) {
        DeviceInRcut(distSq, distVectX, distVectY, distVectZ, gpu_x, gpu_y, gpu_z, particleID, otherParticle, axx, axy, axz, box);
        dist = sqrt(distSq);

        double expConstValue = exp(-1.0 * alphaSq * distSq);
        double qiqj = charge * gpu_particleCharge[otherParticle] * qqFact;
        intraForce = qiqj * lambdaCoef * lambdaCoef / distSq;
        intraForce *= ((erf(alpha * dist) / dist) - constValue * expConstValue);
        forceX -= intraForce * distVectX;
        forceY -= intraForce * distVectY;
        forceZ -= intraForce * distVectZ;
      }
    }
  }
  
  atomicAdd(&gpu_aForceRecx[particleID], forceX);
  atomicAdd(&gpu_aForceRecy[particleID], forceY);
  atomicAdd(&gpu_aForceRecz[particleID], forceZ);
  atomicAdd(&gpu_mForceRecx[moleculeID], forceX);
  atomicAdd(&gpu_mForceRecy[moleculeID], forceY);
  atomicAdd(&gpu_mForceRecz[moleculeID], forceZ);
  
}

