!
!----------------------------------------------
SUBROUTINE grid_shuffle()
!----------------------------------------------
!
!  In PH they have a very simple set of meshes:
!     nksq = nks / 2
!     ALLOCATE(ikks(nksq), ikqs(nksq))
!     DO ik=1,nksq
!        ikks(ik) = 2 * ik - 1
!        ikqs(ik) = 2 * ik
!     ENDDO
!And they just generate the eigenfunctions on those grids with an nscf step for each new qpoint.
!I think the time has come to do the full grid shuffle. One NSCF step to generate all the k's.
!This should be sufficient to do a full GW dispersion calculation. Since all we need are all the ks, k+qs
!and k-qs. The current refold routine RELIES on a uniform and complete q-mesh. The trick is now to reformulate
!the folding so that it applies when we become more selective about which k and q points we want to use and we
!start exploiting symmetry reductions. 
!generate uniform {k} and {k \pm q} grids.
!The {k} grid is taken to coincide with the {q} grid generated
!in gwhs.f90, the {k+q} grid is obtained by folding
!In this way we calculate the occupied states only once in gwhs.f90
 
  USE kinds,           ONLY : DP
  USE units_gw,        ONLY : iuwfc, lrwfc
  USE qpoint,          ONLY : xq, npwq, igkq, nksq, ikks, ikqs
  USE disp,            ONLY : nqs, g0vec, eval_occ, gmap, x_q
  USE klist,           ONLY : xk, wk, nkstot, nks
  USE constants,       ONLY : eps8
  USE wvfct,           ONLY : et, nbnd
  USE gvect,           ONLY : ngm

  IMPLICIT NONE

  INTEGER :: iq, count, i, j, k, ik, ipol, ikk, ikq, ig0, igmg0, nw
  INTEGER :: ig, iw, igp, ierr, ibnd, ios, recl, unf_recl, ikf, fold(nkstot), icounter
  REAL    :: g0(3)
  INTEGER, PARAMETER :: ng0vec = 27
  INTEGER :: shift(nkstot)
  !complex(kind=DP) :: aux (ngm, nbnd_occ) , evq (ngm, nbnd_occ)
  COMPLEX(DP) :: aux (ngm, nbnd) , evq (ngm, nbnd)

! SGW integer, parameter :: nksq = nq, nks = 2 * nksq
! In SGW the xk grid is defined to be the q grid here
! However we already have the k grid from the nscf step
! do ik = 1, nksq
!   xk (:, ik) = xq (:, ik)
!   include spin degeneracy
!   wk = two / float ( nksq )
! enddo
! find folding index and G-vector map
! "nkstot" k- and k+q-points in the IBZ calculated for the phonon sym.
!SGW Pilot  integer, parameter :: nksq = nq, nks = 2 * nksq
! HL gmap appears to be fine here...
!  do ig=1,ngm
!     do ig0=1, ng0vec
!         write(6,*), gmap (ig,ig0)
!     enddo
!  enddo
!  do ik = 1, nksq
!  do ik = 1, nqs

   do ik = 1, nks
    call ktokpmq (xk(:,ik), xq, +1, fold(ik) )
   !
   !  the folding G-vector
   !
    g0 = xk(:, fold(ik)) - ( xk(:,ik) + xq )
   !
   write(6,'("folding index ", 3f7.3)') fold(ik)
   write(6,'("xk original ", 3f7.3)') xk(:, ik)
   write(6,'("xq ", 3f7.3)') xq
   write(6,'("xk folded ", 3f7.3)') xk(:, fold(ik))
   write(6,'("g0 ", 3f7.3)') g0
   !
    shift(ik) = 0
   if(ik.eq.1) then
     write(6,'("g0vec and ig0")')
      do ig0 = 1, ng0vec
      if ( ( abs(g0vec(1,ig0)-g0(1)).lt.eps8 ) .and. &
           ( abs(g0vec(2,ig0)-g0(2)).lt.eps8 ) .and. &
           ( abs(g0vec(3,ig0)-g0(3)).lt.eps8 ) ) then
         shift(ik) = ig0
         write(6,'(3f7.3)') g0vec(:,ig0)
         write(6,'(3f7.3)') g0
      endif
      enddo
    endif
!   if (shift(ik).eq.0) call error ('coulomb','folding vector not found',0)
   enddo
!stop
 
! double the grid and add k+q in the even positions
! PH.x set up the list like this in initialize_gw:
!     DO ik=1,nksq
!        ikks(ik) = 2 * ik - 1
!        ikqs(ik) = 2 * ik
!     ENDDO
 
   do ik = nksq, 1, -1
     ikk = 2 * ik - 1
     ikq = 2 * ik
     xk(:,ikk) = xk(:,ik)
     xk(:,ikq) = xk(:,ik) + xq
     wk(ikk) = wk(ik)
     wk(ikq) = 0.d0
   enddo

  do ik = 1, nksq
     !
     ikk = 2 * ik - 1
     ikq = 2 * ik
     !
     !  the folded k+q
     !
     ikf = fold(ik)
     write(6,*) fold(ik)
     !
     !  read occupied wavefunctions for the folded k+q point
     !  c_{k+q} (G) = c_{k+q+G0} (G-G0)

     !read ( iuwfc, rec = 2 * ikf - 1, iostat = ios) aux
     call davcio(aux, lrwfc, iuwfc, 2 * ikf - 1, -1)


     !  set the phase factor of evq for the folding
     !  WARNING: here we loose some information since
     !  the sphere G-G0 has half the surface outside the
     !  cutoff. It is important to make sure that the cutoff
     !  used is enough to guarantee that the lost half-surface
     !  does not create problems. In the EPW code I solved
     !  the issue by translating the G-sphere in the fft
     !  mapping. It's painful, so I will do it only in extremis.
     !  FG Aug 2008

     write(6,*) shift(ik)
     do ibnd = 1, nbnd
        do ig = 1, ngm
         ! problem with shift(ik)=0 not finding G-vector which maps k+q back on to k. 
         ! the mesh may not be self-contained. 
          igmg0 = gmap (ig, shift(ik))
          if (igmg0.eq.0) then
             evq (ig,ibnd) = (0.0d0, 0.0d0) 
          else
             write(6,*) igmg0
             evq (ig,ibnd) = aux ( igmg0, ibnd)
          endif
        enddo
     enddo
    !stop
    !  and write it again in the right place
    !  HL write ( iuwfc, rec = ikq, iostat = ios) evq

    call davcio (evq, lrwfc, iuwfc, ikq, 1)


    ! DEBUG: here I checked that by running
    ! call eigenstates ( xk(:,ikq), vr, g2kin, evq, eval_occ(:,ikf) )
    ! the evq (wfs at k+q) are the same as those obained above
    ! (within a phase factor and gauge in degenerate cases -
    ! I actually checked sum_ibnd |evq(:,ibnd)|^2 )
    ! the eigenvalues

    ! HL need to switch around these indices. In sgw eval_occ contains the eigenvalues on the initially
    ! generated q grid. In gw.x et contains the NSCF eigenvalues.
    ! SGW:
    ! et(:,ikk) = eval_occ(:,ik)
    ! et(:,ikq) = eval_occ(:,ikf)
    ! Maybe it would be best to define a new temporary array which gets passed to solve 
    ! linter and contains the shifted eigenvalues. since the grid shuffle needs to take place for each q.  
    ! i.e. call gridshuffle(iq) -> etq(:,:) -> call solve_linter (evc,eqv, etq()).

    eval_occ(:, ikk) = et(:,ik)
    eval_occ(:, ikq) = et(:,ikf)
    
  enddo !on k
END SUBROUTINE grid_shuffle