   int box)
 {
   int threadID = blockIdx.x * blockDim.x + threadIdx.x;
-  if(threadID >= pairSize)
-    return;
+
   double distSq;
   double qi_qj_fact;
   double qqFact = 167000.0;

  bool *gpu_isFraction,
  int box)
{
  int threadID = blockIdx.x * blockDim.x + threadIdx.x;
  if(threadID >= pairSize)
    return;

  double distSq;
  double virX, virY, virZ;
  double pRF = 0.0, qi_qj, pVF = 0.0;
  double lambdaVDW = 0.0, lambdaCoulomb = 0.0;
  //tensors for VDW and real part of electrostatic
  gpu_vT11[threadID] = 0.0, gpu_vT22[threadID] = 0.0, gpu_vT33[threadID] = 0.0;
  gpu_rT11[threadID] = 0.0, gpu_rT22[threadID] = 0.0, gpu_rT33[threadID] = 0.0;

  double diff_comx, diff_comy, diff_comz;
  double cutoff = fmax(gpu_rCut[0], gpu_rCutCoulomb[box]);

  int currentCell = blockIdx.x / 27;
  int nCellIndex = blockIdx.x;
  int neighborCell = gpu_neighborList[nCellIndex];

  // calculate number of particles inside neighbor Cell
  int particlesInsideCurrentCell, particlesInsideNeighboringCells;
  int endIndex = neighborCell != numberOfCells - 1 ?
    gpu_cellStartIndex[neighborCell+1] : atomNumber;
  particlesInsideNeighboringCells = endIndex - gpu_cellStartIndex[neighborCell];

  // Calculate number of particles inside current Cell
  endIndex = currentCell != numberOfCells - 1 ?
    gpu_cellStartIndex[currentCell+1] : atomNumber;
  particlesInsideCurrentCell = endIndex - gpu_cellStartIndex[currentCell];

  // total number of pairs
  int numberOfPairs = particlesInsideCurrentCell * particlesInsideNeighboringCells;

  for(int pairIndex = threadIdx.x; pairIndex < numberOfPairs; pairIndex += blockDim.x) {
    int neighborParticleIndex = pairIndex / particlesInsideCurrentCell;
    int currentParticleIndex = pairIndex % particlesInsideCurrentCell;

    int currentParticle = gpu_cellVector[gpu_cellStartIndex[currentCell] + currentParticleIndex];
    int neighborParticle = gpu_cellVector[gpu_cellStartIndex[neighborCell] + neighborParticleIndex];

    if(currentParticle < neighborParticle) {
      if(InRcutGPU(distSq, virX, virY, virZ, gpu_x[currentParticle],
            gpu_y[currentParticle], gpu_z[currentParticle],
            gpu_x[neighborParticle], gpu_y[neighborParticle],
            gpu_z[neighborParticle], xAxes, yAxes, zAxes, xAxes / 2.0,
            yAxes / 2.0, zAxes / 2.0, cutoff, gpu_nonOrth[0],
            gpu_cell_x, gpu_cell_y, gpu_cell_z, gpu_Invcell_x, gpu_Invcell_y,
            gpu_Invcell_z)) {
        int cA = gpu_particleCharge[currentParticle];
        int cB = gpu_particleCharge[neighborParticle];
        int kA = gpu_particleKind[currentParticle];
        int kB = gpu_particleKind[neighborParticle];
        int mA = gpu_particleMol[currentParticle];
        int mB = gpu_particleMol[neighborParticle];

        lambdaVDW = DeviceGetLambdaVDW(mA, kA, mB, kB, box, gpu_isFraction,
            gpu_molIndex, gpu_kindIndex, gpu_lambdaVDW);

        diff_comx = gpu_comx[mA] - gpu_comx[mB];
        diff_comy = gpu_comy[mA] - gpu_comy[mB];
        diff_comz = gpu_comz[mA] - gpu_comz[mB];

        diff_comx = MinImageSignedGPU(diff_comx, xAxes, xAxes / 2.0);
        diff_comy = MinImageSignedGPU(diff_comy, yAxes, yAxes / 2.0);
        diff_comz = MinImageSignedGPU(diff_comz, zAxes, zAxes / 2.0);

        if(electrostatic) {
          qi_qj = cA * cB;
          lambdaCoulomb = DeviceGetLambdaCoulomb(mA, kA, mB, kB, box,
              gpu_isFraction, gpu_molIndex,
              gpu_kindIndex, gpu_lambdaCoulomb);
          pRF = CalcCoulombForceGPU(distSq, qi_qj, gpu_VDW_Kind[0], gpu_ewald[0],
              gpu_isMartini[0], gpu_alpha[box],
              gpu_rCutCoulomb[box], gpu_diElectric_1[0],
              gpu_sigmaSq, sc_coul, sc_sigma_6, sc_alpha,
              sc_power, lambdaCoulomb, gpu_count[0],
              kA, kB);

          gpu_rT11[threadID] = pRF * (virX * diff_comx);
          gpu_rT22[threadID] = pRF * (virY * diff_comy);
          gpu_rT33[threadID] = pRF * (virZ * diff_comz);

          //extra tensor calculations
          gpu_rT12[threadID] = pRF * (0.5 * (virX * diff_comy + virY * diff_comx));
          gpu_rT13[threadID] = pRF * (0.5 * (virX * diff_comz + virZ * diff_comx));
          gpu_rT23[threadID] = pRF * (0.5 * (virY * diff_comz + virZ * diff_comy));
        }

        pVF = CalcEnForceGPU(distSq, kA, kB,
            gpu_sigmaSq, gpu_n, gpu_epsilon_Cn, gpu_rCut[0],
            gpu_rOn[0], gpu_isMartini[0], gpu_VDW_Kind[0],
            gpu_count[0], lambdaVDW, sc_sigma_6, sc_alpha,
            sc_power, gpu_rMin, gpu_rMaxSq, gpu_expConst);

        gpu_vT11[threadID] = pVF * (virX * diff_comx);
        gpu_vT22[threadID] = pVF * (virY * diff_comy);
        gpu_vT33[threadID] = pVF * (virZ * diff_comz);

        //extra tensor calculations
        gpu_vT12[threadID] = pVF * (0.5 * (virX * diff_comy + virY * diff_comx));
        gpu_vT13[threadID] = pVF * (0.5 * (virX * diff_comz + virZ * diff_comx));
        gpu_vT23[threadID] = pVF * (0.5 * (virY * diff_comz + virZ * diff_comy));
      }
    }
  }
}

