                              uint box)
 {
   int atomNumber = currentCoords.Count();
-  int blocksPerGrid, threadsPerBlock;
   double *gpu_particleCharge;
   double *gpu_final_value;
 

  cudaMemcpy(vars->gpu_dy, currentCOMDiff.y, atomNumber * sizeof(double), cudaMemcpyHostToDevice);
  cudaMemcpy(vars->gpu_dz, currentCOMDiff.z, atomNumber * sizeof(double), cudaMemcpyHostToDevice);
  CUMALLOC((void**) &vars->gpu_rT11, imageSize * sizeof(double));
  CUMALLOC((void**) &vars->gpu_rT12, imageSize * sizeof(double));
  CUMALLOC((void**) &vars->gpu_rT13, imageSize * sizeof(double));
  CUMALLOC((void**) &vars->gpu_rT22, imageSize * sizeof(double));
  CUMALLOC((void**) &vars->gpu_rT23, imageSize * sizeof(double));
  CUMALLOC((void**) &vars->gpu_rT33, imageSize * sizeof(double));
  CUMALLOC((void**) &vars->gpu_vT11, imageSize * sizeof(double));
  CUMALLOC((void**) &vars->gpu_vT12, imageSize * sizeof(double));
  CUMALLOC((void**) &vars->gpu_vT13, imageSize * sizeof(double));
  CUMALLOC((void**) &vars->gpu_vT22, imageSize * sizeof(double));
  CUMALLOC((void**) &vars->gpu_vT23, imageSize * sizeof(double));
  CUMALLOC((void**) &vars->gpu_vT33, imageSize * sizeof(double));
  cudaMemcpy(gpu_particleCharge, &particleCharge[0], particleCharge.size() * sizeof(double), cudaMemcpyHostToDevice);

  // Run the kernel...
  threadsPerBlock = 256;
  blocksPerGrid = (int)(imageSize / threadsPerBlock) + 1;
  VirialReciprocalGPU <<< blocksPerGrid,
                      threadsPerBlock>>>(vars->gpu_x,
                                         vars->gpu_y,

  
  if(blockIdx.y == 0) {
    double constant_part = constVal + 1.0 / gpu_hsqrRef[imageID];
    gpu_rT11[imageID] += factor * (1.0 - 2.0 * constant_part *
                                  gpu_kxRef[imageID] * gpu_kxRef[imageID]);
    gpu_rT12[imageID] += factor * (-2.0 * constant_part *
                                  gpu_kxRef[imageID] * gpu_kyRef[imageID]);
    gpu_rT13[imageID] += factor * (-2.0 * constant_part *
                                  gpu_kxRef[imageID] * gpu_kzRef[imageID]);
    gpu_rT22[imageID] += factor * (1.0 - 2.0 * constant_part *
                                  gpu_kyRef[imageID] * gpu_kyRef[imageID]);
    gpu_rT23[imageID] += factor * (-2.0 * constant_part *
                                  gpu_kyRef[imageID] * gpu_kzRef[imageID]);
    gpu_rT33[imageID] += factor * (1.0 - 2.0 * constant_part *
                                  gpu_kzRef[imageID] * gpu_kzRef[imageID]);
  }

