!
   !-----------------------------------------------------------------
   subroutine wigner_seitz (nk1, nk2, nk3, irvec, nrr, ndegen, wslen)
   !-----------------------------------------------------------------
   !
   ! Calculates a grid of points that fall inside of (and eventually 
   ! on the surface of) the Wigner-Seitz supercell centered on the 
   ! origin of the Bravais lattice with primitive translations 
   ! nk1*a_1+nk2*a_2+nk3*a_3
   !  
   ! adapted to f90 by FG from the original f77 version of Marzari, 
   ! Vanderbilt and Souza (remember to stick again the gnu license stuff)
   !
   ! w.r.t. the original version the wigner-seitz vectors are sorted
   ! by increasing lenght in output. In this way the electron and
   ! phonon wigner-seitz vectors should always be the same (even though
   ! the number of them may differ)
   !
   ! BUG FIX: in the case of the tetragonal cell with c>a (LSCO)
   ! the WS points are correctly determined, but there are a few points
   ! with the wrong degeneracies. To avoid this I search for points
   ! in the -2:2 replicas (5^2 replicas). I had the same problem in 
   ! createkmap for the g0vec shift, and also there I have fixed it
   ! by extending the replicas to -2:2 instead of -1:1. @FG May 07
   !
   !-----------------------------------------------------------------
   !
   USE kinds, only : DP
   USE cell_base, ONLY : at, celldm
   implicit none
   !
   integer :: nk1, nk2, nk3, irvec(3,2*nk1*nk2*nk3), ndegen(2*nk1*nk2*nk3), nrr
   real(kind=DP) :: wslen(2*nk1*nk2*nk3)
   !
   ! size of the uniform k mesh (input)
   ! integer components of the ir-th Wigner-Seitz grid point 
   !      in the basis of the lattice vectors (output)
   ! number of Wigner-Seitz grid points (output) 
   ! real-space length, in units of alat
   !
   ! work variables
   integer :: irvec_ (3,2*nk1*nk2*nk3), ndegen_ (2*nk1*nk2*nk3)
   real(kind=DP), parameter :: eps = 1.d-8
   integer :: n1, n2, n3, i1, i2, i3, i, ipol, jpol, ndiff(3), ind(2*nk1*nk2*nk3)
!  real(kind=DP) :: tot, rvec(3), mindist, adot(3,3), dist(27)
   real(kind=DP) :: tot, rvec(3), mindist, adot(3,3), dist(125)
   real(kind=DP), parameter :: bohr2ang = 0.5291772108
   logical :: found
   !
   !  the metric tensor 
   !
   do ipol = 1, 3
    do jpol = 1, 3
      adot (ipol, jpol) = dot_product ( at(:,ipol), at(:,jpol) )
    enddo
   enddo
   !
   ! Loop over grid points r on a unit cell that is 8 times larger than a 
   ! primitive supercell. In the end nrr contains the total number of grids 
   ! points that have been found in the Wigner-Seitz cell
   !
   nrr = 0
!  do n1 = 0, 2*nk1 
!  do n2 = 0, 2*nk2 
!  do n3 = 0, 2*nk3 
   do n1 = 0, 4*nk1 
   do n2 = 0, 4*nk2 
   do n3 = 0, 4*nk3 
     !
!    ! Loop over the 27 points R. R=0 corresponds to i1=i2=i3=1, or icnt=14
     ! Loop over the 5^3 = 125 points R. R=0 corresponds to i1=i2=i3=2, or icnt=63
     !
     i = 0
!    do i1 = 0, 2
!    do i2 = 0, 2
!    do i3 = 0, 2
     do i1 = 0, 4
     do i2 = 0, 4
     do i3 = 0, 4
       i = i + 1
       !
       ! Calculate distance squared |r-R|^2 
       !
       ndiff(1) = n1-i1*nk1
       ndiff(2) = n2-i2*nk2
       ndiff(3) = n3-i3*nk3
       dist(i) = 0.d0
       do ipol = 1, 3
        do jpol = 1, 3
           dist(i) = dist(i) + float(ndiff(ipol))*adot(ipol,jpol)*float(ndiff(jpol))
        enddo
       enddo
       !  
     enddo
     enddo
     enddo
     !
!    ! Sort the  27 vectors R by increasing value of |r-R|^2
     ! Sort the 125 vectors R by increasing value of |r-R|^2
     !
     ! NOTA BENE: hpsort really sorts the dist vector
     ! while the original subroutine by MVS did not. Therefore,
     ! dist(ind(i)) of the original version is here replacerd by
     ! dist(i), while ind(i) is kept.
     !
     ind(1) = 0 ! required for hpsort_eps (see the subroutine)
!    call hpsort_eps(  27, dist, ind, eps)
     call hpsort_eps( 125, dist, ind, eps)
     !
     ! Find all the vectors R with the (same) smallest |r-R|^2;
     ! if R=0 is one of them, then the current point r belongs to 
     ! Wignez-Seitz cell => set found to true
     !
     found = .false.
     i = 1
     mindist = dist(1) 
     do while ( abs(dist(i)-mindist).le.eps ) 
!      if (ind(i).eq.14) found = .true. 
       if (ind(i).eq.63) found = .true. 
       i = i + 1
     enddo
     if (found) then
       nrr = nrr + 1
       ndegen (nrr) = i - 1
!      irvec (1, nrr) = n1 - nk1
!      irvec (2, nrr) = n2 - nk2
!      irvec (3, nrr) = n3 - nk3
       irvec (1, nrr) = n1 - 2*nk1
       irvec (2, nrr) = n2 - 2*nk2
       irvec (3, nrr) = n3 - 2*nk3
     endif
     !
   enddo 
   enddo 
   enddo 
   !
   ! Check the "sum rule"
   !
   tot = 0.d0
   do i = 1, nrr
    tot = tot + 1.d0 / float (ndegen(i)) 
   enddo
   !
   if(abs(tot-float(nk1*nk2*nk3)).gt.eps) call errore &
        ('wigner_seitz','weights do not add up to nk1*nk2*nk3',1)
   !
   ! Hopefully this will never happen, i.e., I think 2*nk1*nk2*nk3 is
   ! an upper bound to the number of lattice points in (or on
   ! the surface of) the Wigner-Seitz supercell
   !
   if(nrr.gt.2*nk1*nk2*nk3) call errore &
     ('wigner_seitz','too many wigseit points, try to increase the bound 2*nk1*nk2*nk3',1)
   !
   ! Now sort the wigner-seitz vectors by increasing magnitude
   !
   do i = 1, nrr
     wslen(i) = 0.d0
     do ipol = 1, 3
      do jpol = 1, 3
         wslen(i) = wslen(i) + float(irvec(ipol,i))*adot(ipol,jpol)*float(irvec(jpol,i))
      enddo
     enddo
     wslen(i) = sqrt(wslen(i))
   enddo 
   !
   ind(1) = 0 ! required for hpsort_eps (see the subroutine)
   call hpsort_eps( nrr, wslen, ind, eps)
   !
   !  now wslen is already sorted, but we still have to sort
   !  irvec and ndegen
   !
   do i = 1, nrr
     ndegen_ (i)  = ndegen(ind(i))
     irvec_ (:,i) = irvec(:,ind(i))
   enddo
   ndegen = ndegen_
   irvec  = irvec_
   !
   return
   end subroutine wigner_seitz
   !---------------------------------------------------------
   !