   dim3 dimBlock(block_size[0], block_size[1]);
   dim3 dimGrid(grid_size[0], grid_size[1]);
   
-
   // real part first
   checkCudaErrors(cudaMemcpy( v_ext_d, v_ext_real, sizeof(float) * nprod, cudaMemcpyHostToDevice));
 

  int j = blockIdx.y * blockDim.y + threadIdx.y; //nvirt
  float en=0.0, fn=0.0, em=0.0, fm=0.0, old_re, old_im;
  double d1p, d1pp, d2p, d2pp, alpha, beta;
  int m, index;

  if (i < nfermi)
  {

  checkCudaErrors(cublasSgemm(handle_cublas, CUBLAS_OP_T, CUBLAS_OP_N, nvirt, nfermi, norbs, &alpha, &X4_d[vstart*norbs], norbs, nb2v_d,
        norbs, &beta, nm2v_im_d, nvirt));


  // Normalization!!
  normalize_energy_gpu<<<dimGrid, dimBlock>>>(ksn2e_d, ksn2f_d, omega_re, omega_im, 
      nm2v_re_d, nm2v_im_d, nfermi, norbs, nvirt, vstart);

  // Going back: real part first
  checkCudaErrors(cublasSgemm(handle_cublas, CUBLAS_OP_N, CUBLAS_OP_N, norbs, nfermi, nvirt, &alpha, &X4_d[vstart*norbs], norbs, nm2v_re_d,
