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

  gpu_REn[threadID] = 0.0;
  gpu_LJEn[threadID] = 0.0;
  double cutoff = fmax(gpu_rCut[0], gpu_rCutCoulomb[box]);
  if(InRcutGPU(distSq, virX, virY, virZ, gpu_x[gpu_pair1[threadID]],
        gpu_y[gpu_pair1[threadID]], gpu_z[gpu_pair1[threadID]],
        gpu_x[gpu_pair2[threadID]], gpu_y[gpu_pair2[threadID]],
        gpu_z[gpu_pair2[threadID]], 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[gpu_pair1[threadID]];
    int cB = gpu_particleCharge[gpu_pair2[threadID]];
    int kA = gpu_particleKind[gpu_pair1[threadID]];
    int kB = gpu_particleKind[gpu_pair2[threadID]];
    int mA = gpu_particleMol[gpu_pair1[threadID]];
    int mB = gpu_particleMol[gpu_pair2[threadID]];

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

    if(electrostatic) {
      qi_qj_fact = cA * cB * qqFact;
      lambdaCoulomb = DeviceGetLambdaCoulomb(mA, kA, mB, kB, box,
          gpu_isFraction, gpu_molIndex,
          gpu_kindIndex, gpu_lambdaCoulomb);
      gpu_REn[threadID] = CalcCoulombGPU(distSq,
          kA,
          kB,
          qi_qj_fact, gpu_rCutLow[0],
          gpu_ewald[0], gpu_VDW_Kind[0],
          gpu_alpha[box],
          gpu_rCutCoulomb[box],
          gpu_isMartini[0],
          gpu_diElectric_1[0],
          lambdaCoulomb,
          sc_coul,
          sc_sigma_6,
          sc_alpha,
          sc_power,
          gpu_sigmaSq[threadID],
          gpu_count[0]);
    }
    gpu_LJEn[threadID] = CalcEnGPU(distSq,
        kA,
        kB,
        gpu_sigmaSq, gpu_n, gpu_epsilon_Cn,
        gpu_VDW_Kind[0], gpu_isMartini[0],
        gpu_rCut[0], gpu_rOn[0], gpu_count[0],
        lambdaVDW,
        sc_sigma_6, sc_alpha, sc_power, gpu_rMin,
        gpu_rMaxSq, gpu_expConst);
    if(electrostatic) {
      double coulombVir = CalcCoulombForceGPU(distSq, qi_qj_fact,
          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);
      forceRealx = virX * coulombVir;
      forceRealy = virY * coulombVir;
      forceRealz = virZ * coulombVir;
    }
    double 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);
    forceLJx = virX * pVF;
    forceLJy = virY * pVF;
    forceLJz = virZ * pVF;

    atomicAdd(&gpu_aForcex[gpu_pair1[threadID]], forceRealx + forceLJx);
    atomicAdd(&gpu_aForcey[gpu_pair1[threadID]], forceRealy + forceLJy);
    atomicAdd(&gpu_aForcez[gpu_pair1[threadID]], forceRealz + forceLJz);
    atomicAdd(&gpu_aForcex[gpu_pair2[threadID]], -1.0 * (forceRealx + forceLJx));
    atomicAdd(&gpu_aForcey[gpu_pair2[threadID]], -1.0 * (forceRealy + forceLJy));
    atomicAdd(&gpu_aForcez[gpu_pair2[threadID]], -1.0 * (forceRealz + forceLJz));

    atomicAdd(&gpu_mForcex[mA], forceRealx + forceLJx);
    atomicAdd(&gpu_mForcey[mA], forceRealy + forceLJy);
    atomicAdd(&gpu_mForcez[mA], forceRealz + forceLJz);
    atomicAdd(&gpu_mForcex[mB], -1.0 * (forceRealx + forceLJx));
    atomicAdd(&gpu_mForcey[mB], -1.0 * (forceRealy + forceLJy));
    atomicAdd(&gpu_mForcez[mB], -1.0 * (forceRealz + forceLJz));
  }
}


  bool *gpu_isFraction,
  int box)
{
  double distSq;
  double virX, virY, virZ;
  double pRF = 0.0, qi_qj, pVF = 0.0;
  double lambdaVDW = 0.0, lambdaCoulomb = 0.0;
  int threadID = blockIdx.x * blockDim.x + threadIdx.x;
  //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;
