      subroutine gvlist(alat,plat,tau,n1,n2,n3,gmax,job,ngmx,ng,kv,gv,
     .  igv,igv2)
C- Create a list of G vectors for a specified uniform mesh and cutoff
C ----------------------------------------------------------------------
Ci Inputs
Ci   alat     Lattice constant
Ci   plat     Real-space primitive lattice vectors
Ci   tau      offset subtracted from gv to measure length
Ci            Used if 100s digit job is 1 or 2
Ci   n1..n3   no. divisions along the three lattice vectors
Ci   job      1s digit
Ci              0 return ng only
Ci              1 return kv and igv
Ci              2 return kv and igv2
Ci              4 return kv and gv
Ci              8 return kv and gv, and sort list
Ci                any combination of 1,2,4,8 is allowed
Ci            100s digit
Ci              0 specify cutoff as |G|<Gmax
Ci              1 specify cutoff as |G+tau|<Gmax
Ci              2 specify cutoff as |G-tau|<Gmax
Ci              Add 4 if to return gv = G+/-tau
Ci            1000s digit
Ci              not used
Ci   gmax     Energy cutoff
Ci   ngmx     Leading dimension of gv,kv,igv; must be at least ng
Ci            if any of these are to be generated.
Co Outputs (which quantities are ouptut depends on job)
Co   ng       Number of lattice vectors
Co   gv       list of reciprocal lattice vectors G
Co   kv       indices for gather/scatter operations.
Co            kv(ig,i=1,2,3) for vector ig=1..ng points to which entry
Co            on the mesh (i1,i2,i3) vector ig belongs:
Co              1<=i1<=n1; 1<=i2<=n2;  1<=i3<=n3
Co   igv      list of reciprocal lattice vectors G, represented as
Co            three integers (multiples of qlat).
Co            gv and igv are related by:
Co              gv(1:3,1:ng) = 2*pi/alat * (qlat * igv(1:ng))
Co   igv2     same as igv except first and second columns are permuted
Cr Remarks
Cr   Collects a list of q + reciprocal lattice vectors (G+q) out of
Cr   a uniform 3D mesh lie within a cutoff gmax.  Sorted by length.
Cr   Criterion for cutoff is |(gv-tau)|<gmax
Cr   Use with the following:
Cr     call gvgetf(ng,1,kv,k1,k2,k3,c,c0)
Cr       to collect elements into list c0 from 3D array c
Cr     call gvputf(ng,1,kv,k1,k2,k3,c0,c)
Cr       to poke elements from list c0 into 3D array c
Cr     call gvaddf(ng,kv,k1,k2,k3,c0,c)
Cr       to add elements from list c0 into 3D array c
Cr
Cr   The following routines are designed to work either with a
Cr   specified uniform mesh of G vectors, specified by primitive
Cr   lattice vectors and numbers of divisions n1,n2,n3 (gvcutof,
Cr   gvlist) or without it (gvlst2).
Cr
Cr     gvctof takes as input a uniform mesh of points n1..n3. gvctof
Cr            generates an appropriate energy cutoff gmax, and counts
Cr            the number of G vectors within the cutoff.  The list of
Cr            vectors is generated by looping over all points
Cr            0..n1,0..n2,0..n3, shortens each vector, and retaining
Cr            those for which G<Gmax.  Because only points on the
Cr            specified mesh are considered, there is no guarantee
Cr            that all vectors G<Gmax will be included in the list.
Cr
Cr     gvlist takes as input uniform mesh of G vectors specified by
Cr            n1..n3 and an energy cutoff Gmax, and creates a list of
Cr            G-vectors that lie within Gmax.  gvlist operates in the
Cr            same way as gvctof, and generates the same list of
Cr            vectors.  The list of vectors is sorted by increasing
Cr            length.
Cr
Cr     gvlst2 is designed to work without any constraint that the list
Cr            of vectors map onto a specified mesh n1..n3.  It takes
Cr            as input an energy cutoff Gmax, and returns a list of
Cr            all G vectors whose length G is G<Gmax.  There is no
Cr            shortening of vectors; the number of divisions n1..n3
Cr            needed to encompass all vectors cannot be specified, but
Cr            is output by gvlst2.  Thus, this routine is only
Cr            suitable in cases where there is no need for the vectors
Cr            to map onto a specified mesh.  Offset q in gvlst2
Cr            corresponds to -tau in gvlist.
Cr
Cu Updates
Cu   20 Jul 02 Bug fix for nonzero tau
Cu   31 May 01 gvlist was significantly altered.  See ogvlst for
Cu             original version.
Cu   07 Feb 01 changed gmax tolerance to be consistent with gvlst2
Cu   10 Jun 00 added extra argument to gvgetf and gvputf
Cu    2 Sep 98 Adapted from nfp gvlist.f
C ----------------------------------------------------------------------
C     implicit none
C ... Passed parameters
      integer job,n1,n2,n3,ng,ngmx,kv(ngmx,3),igv(ngmx,3),igv2(3,*)
      double precision alat,gmax,gv(ngmx,3),plat(3,3),tau(3)
C ... Local parameters
      logical ligv,ligv2,lgv,lsort,lgpq
      integer m,ig,j1,j2,j3,iprint,ogg,okk,oiwk,PRTG,stdo,lgunit,job0
      double precision qlat(3,3),plat1(3,3),qlat1(3,3),qpg(3),gs(3),pi,
     .  tpiba,tol,gg,vol,gmax2,q(3),s1,s2,s3,ddot
      character*256 outs
      parameter (PRTG=30)
C ... Heap
      integer w(1)
      common /w/ w
C ... External calls
      external awrit0,awrit4,awrit6,defrr,dinv33,dpcopy,dpzero,isanrg,
     .         pvglst,rlse,rx,shorbz
cccccccccccccc
      integer,save:: i,iyy
      integer,allocatable:: ixx(:,:)

      integer,parameter :: m_gslist=10
      integer:: n_gslist,i_gslist,if1,if2,if3
      double precision:: gslist(3,m_gslist)
      integer:: gsflags(3)
      logical:: fullmesh
cccccccccccccc

C ... Setup
      stdo = lgunit(1)
      job0  = mod(job,100)
      ligv  = mod(job0,2)   .ne. 0
      ligv2 = mod(job0/2,2) .ne. 0
      lgv   = mod(job0/4,4) .ne. 0
      lsort = mod(job0/8,2) .ne. 0
      pi = 4d0*datan(1d0)
      call dinv33(plat,1,qlat,vol)
      tpiba = 2*pi/alat
      tol = 1d-8
      m = mod(mod(job/100,10),4)
      lgpq = mod(job/100,10) .gt. 4

      call isanrg(job0,0,11,'gvlist','1+10s digit job',.true.)
      call isanrg(m,0,2,'gvlist','100s digit job',.true.)
      call dpzero(q,3)
      if (m .eq. 1) then
        call dpcopy(tau,q,1,3,1d0)
      elseif (m .eq. 2) then
        call dpcopy(tau,q,1,3,-1d0)
      endif

C ... Basis vectors for real-space mesh and recip-space supercell
      do  10  m = 1, 3
        plat1(m,1) = plat(m,1)/n1
        plat1(m,2) = plat(m,2)/n2
        plat1(m,3) = plat(m,3)/n3
        qlat1(m,1) = qlat(m,1)*n1
        qlat1(m,2) = qlat(m,2)*n2
        qlat1(m,3) = qlat(m,3)*n3
   10 continue

C --- Loop through g vectors, shorten, count and keep if within gmax ---
      ig = 0
      gmax2 = (gmax-tol)**2/tpiba**2

      if(fullmesh()) then !kino's path. Something wrong
                          ! ---> Kino think this is OK even for fullmesh()=F. 
                          ! But this path did not go through make check.

       do  21  j1 = 0, n1-1
       do  21  j2 = 0, n2-1
       do  21  j3 = 0, n3-1
       n_gslist=0
       do 21 if1=0,1
       do 21 if2=0,1
       do 21 if3=0,1
       gsflags=(/if1,if2,if3/)
       qpg(1) = j1*qlat(1,1) + j2*qlat(1,2) + j3*qlat(1,3) + q(1)
       qpg(2) = j1*qlat(2,1) + j2*qlat(2,2) + j3*qlat(2,3) + q(2)
       qpg(3) = j1*qlat(3,1) + j2*qlat(3,2) + j3*qlat(3,3) + q(3)
        call shorbzm(gsflags,qpg,gs,qlat1,plat1)
        do i_gslist=1,n_gslist
            if ( gs(1).eq.gslist(1,i_gslist) .and.
     .            gs(2).eq.gslist(2,i_gslist) .and.
     .            gs(3).eq.gslist(3,i_gslist) ) then
              goto 21 
            endif
        enddo
        n_gslist=n_gslist+1
        gslist(1:3,n_gslist)=gs(1:3)
c        if (n_gslist>1)then
c          write(stdo,'(a,3f15.5)')'new gs=',gs
c        endif

C        print 333, j1,j2,j3,gs(1)+tau(1),gs(2)+tau(2),gs(3)+tau(3)
C  333   format(3i4,3f12.6)
        gg = (gs(1)**2+gs(2)**2+gs(3)**2)
        if (gg .le. gmax2) then
          ig = ig+1

C         if (abs(gs(1)+5) .lt.1d-4 .and. gs(3).lt.-10.8) then
C           print *, j1,j2,j3
C           print *, gs
C         endif

        if (job0 .ne. 0) then

          if (ig .gt. ngmx) call rx('gvlist: ng exceeds ngmx')

          kv(ig,1) = j1+1
          kv(ig,2) = j2+1
          kv(ig,3) = j3+1

C     ... Copy lattice vectors to gv,igv,igv2, depending on job
          if (ligv .or. ligv2 .or. .true.) then
            s1 = (gs(1)-q(1))*plat(1,1)
     .         + (gs(2)-q(2))*plat(2,1)
     .         + (gs(3)-q(3))*plat(3,1)
            s2 = (gs(1)-q(1))*plat(1,2)
     .         + (gs(2)-q(2))*plat(2,2)
     .         + (gs(3)-q(3))*plat(3,2)
            s3 = (gs(1)-q(1))*plat(1,3)
     .         + (gs(2)-q(2))*plat(2,3)
     .         + (gs(3)-q(3))*plat(3,3)
c#if DEBUG
c            if ((nint(s1)-s1)**2+(nint(s2)-s2)**2+(nint(s3)-s3)**2
c     .        .gt. 1d-12) call rx('bug in gvlist')
c#endif
            if (ligv) then
              igv(ig,1) = nint(s1)
              igv(ig,2) = nint(s2)
              igv(ig,3) = nint(s3)
            endif
            if (ligv2) then
              igv2(1,ig) = nint(s1)
              igv2(2,ig) = nint(s2)
              igv2(3,ig) = nint(s3)
            endif
          endif

          if (lgv) then
            if (lgpq) then
              gv(ig,1) = gs(1)
              gv(ig,2) = gs(2)
              gv(ig,3) = gs(3)
            else
              gv(ig,1) = gs(1) - q(1)
              gv(ig,2) = gs(2) - q(2)
              gv(ig,3) = gs(3) - q(3)
            endif
          endif

        endif
        endif
   21  continue

      else

       do  20  j1 = 0, n1-1
       do  20  j2 = 0, n2-1
       do  20  j3 = 0, n3-1
        qpg(1) = j1*qlat(1,1) + j2*qlat(1,2) + j3*qlat(1,3) + q(1)
        qpg(2) = j1*qlat(2,1) + j2*qlat(2,2) + j3*qlat(2,3) + q(2)
        qpg(3) = j1*qlat(3,1) + j2*qlat(3,2) + j3*qlat(3,3) + q(3)
        call shorbz(qpg,gs,qlat1,plat1)
C        print 333, j1,j2,j3,gs(1)+tau(1),gs(2)+tau(2),gs(3)+tau(3)
C  333   format(3i4,3f12.6)
        gg = (gs(1)**2+gs(2)**2+gs(3)**2)
        if (gg .le. gmax2) then
          ig = ig+1
C         if (abs(gs(1)+5) .lt.1d-4 .and. gs(3).lt.-10.8) then
C           print *, j1,j2,j3
C           print *, gs
C         endif

        if (job0 .ne. 0) then

          if (ig .gt. ngmx) call rx('gvlist: ng exceeds ngmx')

          kv(ig,1) = j1+1
          kv(ig,2) = j2+1
          kv(ig,3) = j3+1

C     ... Copy lattice vectors to gv,igv,igv2, depending on job
          if (ligv .or. ligv2 .or. .true.) then
            s1 = (gs(1)-q(1))*plat(1,1)
     .         + (gs(2)-q(2))*plat(2,1)
     .         + (gs(3)-q(3))*plat(3,1)
            s2 = (gs(1)-q(1))*plat(1,2)
     .         + (gs(2)-q(2))*plat(2,2)
     .         + (gs(3)-q(3))*plat(3,2)
            s3 = (gs(1)-q(1))*plat(1,3)
     .         + (gs(2)-q(2))*plat(2,3)
     .         + (gs(3)-q(3))*plat(3,3)
c#if DEBUG
c            if ((nint(s1)-s1)**2+(nint(s2)-s2)**2+(nint(s3)-s3)**2
c     .        .gt. 1d-12) call rx('bug in gvlist')
c#endif
            if (ligv) then
              igv(ig,1) = nint(s1)
              igv(ig,2) = nint(s2)
              igv(ig,3) = nint(s3)
            endif
            if (ligv2) then
              igv2(1,ig) = nint(s1)
              igv2(2,ig) = nint(s2)
              igv2(3,ig) = nint(s3)
            endif
          endif

          if (lgv) then
            if (lgpq) then
              gv(ig,1) = gs(1)
              gv(ig,2) = gs(2)
              gv(ig,3) = gs(3)
            else
              gv(ig,1) = gs(1) - q(1)
              gv(ig,2) = gs(2) - q(2)
              gv(ig,3) = gs(3) - q(3)
            endif
          endif

        endif
        endif
   20  continue

      endif
      ng = ig

C ... Printout
      if (iprint() .ge. PRTG) then
#if AWRITE
        s1 = alat*sqrt(ddot(3,plat(1,1),1,plat(1,1),1))/n1
        s2 = alat*sqrt(ddot(3,plat(1,2),1,plat(1,2),1))/n2
        s3 = alat*sqrt(ddot(3,plat(1,3),1,plat(1,3),1))/n3
        call awrit4('%x%N GVLIST: gmax = %,1;4g'//
     .    ' a.u. created %i vectors of %i (%i%%)',
     .    outs,len(outs),0,gmax,ng,n1*n2*n3,(ng*100)/(n1*n2*n3))
        call awrit6('%a%N%9fmesh has %i x %i x %i'//
     .    ' divisions; length %,3;3d, %,3;3d, %,3;3d',outs,
     .    len(outs),0,n1,n2,n3,s1,s2,s3)
        call awrit0('%a',outs,len(outs),-stdo)
#else
      if (iprint() .ge. PRTG) print 861, gmax,ng,n1*n2*n3
  861 format(/' gvlist: cutoff radius ',f7.3,' gives',i7,
     .  '   recips of max',i7)
#endif
      endif

C ... Sort the list of vectors
      if (.not. lsort) return
      call defrr(ogg,   ng)
      call defrr(okk,   ng)
      call defrr(oiwk,  ng)
      if (lgpq) call dpzero(q,3)
      call pvglst(w(ogg),w(okk),w(oiwk),q,ngmx,ng,ligv,ligv2,gv,kv,igv,
     .  igv2)
      call rlse(ogg)

c$$$ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c$$$c      if(job==509) print *,'509q=',q
c$$$c      print *,'ggg: 1 job=',job
c$$$      allocate(ixx(ng,3))
c$$$      do i=1,ng
c$$$         ixx(i,1)=i*1
c$$$         ixx(i,2)=i*2
c$$$         ixx(i,3)=i*3
c$$$      enddo
c$$$c      print *,'ggg: x job=',job
c$$$      if(job==509) then
c$$$        iyy=iyy+1
c$$$        write(6,"('ggg509:',5i12,d13.5)") ngmx,ng,sum(ixx*kv(1:ng,:)),
c$$$     &  sum(ixx*igv(1:ng,1:3)),
c$$$c     &  sum(transpose(ixx)*igv2(1:3,1:ng)),
c$$$     &  sum(ixx*gv(1:ng,:))
c$$$         if(iyy<5) then
c$$$          do i=1,ng
c$$$c         write(6,"('ggg509:',i5,4d13.5)")i,gv(i,:),sqrt(sum(gv(i,:)**2))
c$$$          write(6,"('ggg509:',i5,14d13.5)")i,gv(i,:)-tau,
c$$$     &   sqrt( sum((gv(i,:)-tau)**2))
c$$$          enddo
c$$$         endif
c$$$      endif
c$$$      if(job==8) then
c$$$      write(6,"('ggg008:',3i12,d13.5)") ngmx,ng,sum(ixx*kv(1:ng,:)),
c$$$     &  sum(ixx*gv(1:ng,:))
c$$$      endif
c$$$      deallocate(ixx)
c$$$ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      end

      subroutine pvglst(gg,kk,iwk,p,ngmx,ng,ligv,ligv2,gv,kv,igv,igv2)
C- Kernel called by gvlist to sort gv and kv
C ----------------------------------------------------------------------
Ci Inputs
Ci   gg    :double precision work array of length ng
Ci   kk    :integer work array of length ng
Ci   iwk   :integer work array of length ng : returns permutation of gg
Ci   p     :offset subtracted from g to measure length
Ci   ngmx  :dimensions gv,kv,igv
Ci   ng    :number of G-vectors
Ci   ng    :number of group operations
Ci   ligv  :sort igv also
Ci   ligv2 :sort igv2 also
Cio Inputs/Outputs
Ci         :the following arrays are sorted:
Ci   gv    :list of reciprocal lattice vectors G (gvlist.f)
Ci   kv    :indices for gather/scatter operations (gvlist.f)
Ci   igv   :G vectors in integer format
Ci   igv2  :like igv, but array dimensions inverted
Cr Remarks
Cr
Cu Updates
Cu   20 Jul 02 Bug fix for nonzero tau
Cu   01 Jun 01  Added igv,igv2
C ----------------------------------------------------------------------
      implicit none
      logical ligv,ligv2
      integer ngmx,ng,kk(ng),iwk(ng),kv(ngmx,3),igv(ngmx,3),igv2(3,*)
      double precision gv(ngmx,3),gg(ng),p(3)
C Local variables
      integer ig,m,jg,PRTG,iprint
      parameter (PRTG=80)
c      parameter (PRTG=0)
C ... External calls
      external dvshel

      do  ig = 1, ng
       gg(ig)= (gv(ig,1)+p(1))**2+(gv(ig,2)+p(2))**2+(gv(ig,3)+p(3))**2
ctakao Apr2009. Ordering is clealy fixed---without this, ordering can be ambigous when gg is degenerated.
       gg(ig)= gg(ig) *(1d0 + 1d-15*ig)
      enddo

C     dvshel seems to work a little faster for vectors this size
      call dvshel(1,ng,gg,iwk,1)
C     call dvheap(1,ng,gg,iwk,0d0,11)

C ... Rearrange gv,kv,igv,igv2
      do  20  m = 1, 3
        do  22  ig = 1, ng
          jg = iwk(ig)+1
          gg(ig) = gv(jg,m)
          kk(ig) = kv(jg,m)
   22   continue
        do  24  ig = 1, ng
          gv(ig,m) = gg(ig)
          kv(ig,m) = kk(ig)
   24   continue
        if (ligv) then
          do  25  ig = 1, ng
            jg = iwk(ig)+1
            kk(ig) = igv(jg,m)
   25     continue
          do  26  ig = 1, ng
            gv(ig,m) = gg(ig)
            igv(ig,m) = kk(ig)
   26     continue
        endif
        if (ligv2) then
          do  27  ig = 1, ng
            jg = iwk(ig)+1
            kk(ig) = igv2(m,jg)
   27     continue
          do  28  ig = 1, ng
            gv(ig,m) = gg(ig)
            igv2(m,ig) = kk(ig)
   28     continue
        endif
   20 continue

      if (iprint() .ge. PRTG) then
      print 333, p
  333 format(' gvlist: g vectors after sorting: p=',3f12.6)
c      do  30  ig = 1, min(ng,50)
      do  30  ig = 1, ng
        gg(ig)=(gv(ig,1)+p(1))**2+(gv(ig,2)+p(2))**2+(gv(ig,3)+p(3))**2
C        print 550, ig,gv(ig,1)+p(1),gv(ig,2)+p(2),gv(ig,3)+p(3),
        print 550, ig,gv(ig,1),gv(ig,2),gv(ig,3),
     .     kv(ig,1),kv(ig,2),kv(ig,3),sqrt(gg(ig))
  550   format(i5,3f11.5,3i6,f11.5)
   30 continue
C     pause
      endif

      end

      subroutine gvgetf(ng,n,kv,k1,k2,k3,c,c0)
C- Gathers Fourier coefficients from 3D array c into list c0.
C     implicit none
      integer ng,n,k1,k2,k3,kv(ng,3)
      double complex c0(ng,n),c(k1,k2,k3,n)
      integer ig,i,j1,j2,j3

      do  10  i = 1, n
      do  10  ig = 1, ng
        j1 = kv(ig,1)
        j2 = kv(ig,2)
        j3 = kv(ig,3)
        c0(ig,i) = c(j1,j2,j3,i)
   10 continue
      end
      subroutine gvputf(ng,n,kv,k1,k2,k3,c0,c)
C- Pokes Fourier coefficients from gathered list c0 into 3D array c.
C     implicit none
      integer ng,n,k1,k2,k3,kv(ng,3)
      double complex c0(ng,n),c(k1,k2,k3,n)
      integer ig,i,j1,j2,j3

      call dpzero(c,2*k1*k2*k3*n)
      do  10  i = 1, n
      do  10  ig = 1, ng
        j1 = kv(ig,1)
        j2 = kv(ig,2)
        j3 = kv(ig,3)
        c(j1,j2,j3,i) = c0(ig,i)
   10 continue
      end
      subroutine gvaddf(ng,kv,k1,k2,k3,c0,c)
c  Adds Fourier coefficients from gathered list c0 into 3D array c.
C     implicit none
      integer ng,k1,k2,k3,kv(ng,3)
      double complex c0(ng),c(k1,k2,k3)
      integer ig,j1,j2,j3
      do  10  ig = 1, ng
        j1 = kv(ig,1)
        j2 = kv(ig,2)
        j3 = kv(ig,3)
        c(j1,j2,j3) = c(j1,j2,j3) + c0(ig)
   10 continue
      end

C      subroutine gvlstp(p,ng,gv,sgv)
CC- Adds offset to lattice vectors, and resorts by increasing length
CC ---------------------------------------------------------------
CCi Inputs:
CCi  gv,ng:  lattice vectors, and number
CCi  p:      vector to subtract from gv
CCo  sgv:    lattice vectors, sorted by increasing (p-gv)**2
CCr Remarks
CCr  Shell sorting is apparently faster than heap sorting,
CCr  because the starting vectors are approximately sorted already.
CC ---------------------------------------------------------------
C      implicit none
C      integer ng
C      double precision p(3),gv(3,ng),sgv(3,ng)
C      integer ir,oiwk
C      integer w(1)
C      common /w/ w
C
C      do  10  ir = 1, ng
C   10 sgv(ir,1) =
C     .  (p(1)-gv(1,ir))**2 + (p(2)-gv(2,ir))**2 + (p(3)-gv(3,ir))**2
C      call defi(oiwk, ng*3)
C      call dvshel(1,ng,sgv,w(oiwk),1)
C      call dvperm(3,ng,gv,sgv,w(oiwk),.false.)
C      call rlse(oiwk)
C
C      end
C      subroutine rx(string)
C      character *(*) string
C      print *, string
C      stop
C      end

