!
!----------------------------------------------
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 there will come a time when we need 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
!FG: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)
!SGW:
!FG: 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