   int molCountInBox = moleculeInvolved.size();
   int *gpu_moleculeInvolved;
   // Each block would handle one molecule
-  int threadsPerBlock = 64;
+  int threadsPerBlock = 32;
   int blocksPerGrid = molCountInBox;
 
   CUMALLOC((void **) &gpu_moleculeInvolved, molCountInBox * sizeof(int));

  }

  __syncthreads();
  // use strid of blockDim.x, which is 64
  // each thread handles one atom rotation
  for(atomIdx = startIdx + threadIdx.x; atomIdx < endIdx; atomIdx += blockDim.x) {
    double3 coor = make_double3(gpu_x[atomIdx], gpu_y[atomIdx], gpu_z[atomIdx]);

  double axisy = roty * (1.0 / rotLen);
  double axisz = rotz * (1.0 / rotLen);
  double matrix[3][3], cross[3][3], tensor[3][3];
  double halfAxx = axx * 0.5;
  double halfAxy = axy * 0.5;
  double halfAxz = axz * 0.5;

  // build cross
  cross[0][0] = 0.0;
