!-----------------------------------------------------------------------
subroutine setphases ( kunit, ik0, nng, unimat )
  !---------------------------------------------------------------------
  !
  ! This subroutine gives the rotations which bring the
  ! eigenstates of the Hamiltonian (evc) into a uniquely defined
  ! gauge (independent of the machine or the numerical libraries)
  !
  ! This is done by diagonalizing a fake perturbation within the
  ! degenerate hamiltonian manifold. The eigenstates are assumed
  ! to be labeled in order of increasing energy (as it is usually 
  ! the case).
  !  
  ! The first step is the diagonalization of the perturbation.
  ! This still leaves a phase arbitrariness, which is the fixed
  ! by requiring some c(G) be real.
  !
  ! It should be :  |psi^\prime_i>  =  \sum_j U_{ji} * |psi_j>
  ! and A^\prime = U^\dagger * A * U = matmul(conjg(transpose(U)),matmul(A,U)) 
  ! with U = unimat
  !
  ! Only for 1 proc/pool (This limitation can be removed)
  !
  ! Feliciano Giustino May 2006
  ! 
#include "f_defs.h"
  !
#ifdef __PARA
  USE mp_global, only : my_pool_id, nproc_pool, me_pool
  use mp, only : mp_barrier
#endif
  use io_global, only : stdout
  USE kinds, only : DP
  USE wvfct, ONLY : et
  USE ions_base, ONLY : nat
  use phcom, ONLY : igkq, evq, lrdrho, iudvscf0, fildvscf0
  use control_flags, ONLY : iverbosity
  use pwcom
  use wavefunctions_module, ONLY: evc
  !
  implicit none
  !
  INTEGER, PARAMETER :: nnglen = 100, ndig = 5
  !   reduced size of evc, evq, this is just for setting individual phases 
  !     this may be taken = 1 in principle, but if C(G1) = 0 we are in trouble, so...
  !   number of digits used in integer comparisons to decide the largest c(G) in the set
  COMPLEX(KIND=DP), POINTER :: evtmp(:,:)     
  !   pointer to wavefunctions in the PW basis (either evc or evq)
  INTEGER, POINTER :: igktmp(:)
  !   correspondence k+G <-> G
  INTEGER :: nng, kunit, ik0, nset, ndeg(nbnd)
  !   size of evc, evq ( this is npw or npwq from elphel2) 
  !   granularity of k point distribution
  !   the current k point
  !   number of degenerate subsets for this k point
  !   degeneracy of each subset
  COMPLEX(kind=DP) :: c (nnglen,nbnd) , unimat (nbnd,nbnd), cnew (nnglen,nbnd)
  !   the chunk of evc used for setting the phases
  !   the global rotation matrix
  !   the rotated chunks
  COMPLEX(kind=DP), ALLOCATABLE ::  u (:,:)
  !   the rotation matrix for the degenarate subset
  REAL(kind=DP) :: deltav (nrxxs)
  !   the fake (local) perturbation in real space, it is real to guarantee Hermiticity
  COMPLEX(kind=DP) :: v1 (nrxxs,nspin), v2 (nrxxs,nspin), v3 (nrxxs,nspin)
  !   tmp matrices used to build deltav
  REAL(kind=DP), PARAMETER :: eps = 1.d-5, epsc = 1.d-8
  !   threshold for deciding whether two states are degenerate
  !   threshold for deciding whether c(G) = 0
  !
  ! variables for lapack ZHPEVX 
  !
  integer                       :: neig, info
  integer, allocatable          :: ifail(:), iwork(:)
  real(kind=DP), allocatable    :: w(:), rwork(:)
  complex(kind=DP), allocatable :: up(:), cwork(:), cz(:,:)
  !
  ! work variables
  !
  complex(kind=DP), parameter :: ci = (0.d0, 1.d0), czero = (0.d0, 0.d0)
  complex(kind=DP) :: unitcheck (nbnd, nbnd), ctmp
  real(kind=DP) :: maxnorm, theta, meantheta, tmp, minphase
  integer :: ir, ig, igmax, ibnd, jbnd, mbnd, pbnd, iset, n0, n1, n2, &
     nsmall, ntheta, ibndg, jbndg, itmp, imaxnorm
  logical :: tdegen, exst
  character(len=10) :: errstr
  complex(kind=DP) :: aux1 (nrxxs), deltavpsi (nng)
  complex(kind=DP) :: ZDOTC

#ifdef __PARA
  !
  ! In the case of multiple procs per pool, the fake perturbation defined 
  ! by deltav (ir) = ir will depend on the division of nr1x*nr2x*nr3x
  ! (This can be removed by defined a serious perturbation below...)
  !
  if (nproc_pool>1) call errore &
       ('setphases', 'only one proc per pool to guarantee the same gauge', 1)
#endif
  !
  ! initialize
  !
  ! if we are working with a k (k+q) point, evtmp points to evc (evq)
  ! 
  if ( (kunit.eq.1) .or. ((kunit.eq.2).and.(ik0-2*(ik0/2).eq.1)) ) then
    evtmp  => evc
    igktmp => igk
  else
    evtmp  => evq
    igktmp => igkq
  endif
  c = evtmp(1:nnglen,:)
  !  
  ! build the fake perturbation
  !
  ! This is crap but it works fine. I think it is important to normalize 
  ! the pert to 1 so as to easily keep control on the eigenvalues.
  ! According to wikipedia the word crap is from latin "crappa".
  ! If you do something better you should take into account the 
  ! parallalization on nrxxs.
  !  
!  do ir = 1, nrxxs
!    deltav(ir) = float(ir)/float(nrxxs)
!  enddo
  !
  ! the above is not ok for my BC53 structure... A better choice:
  ! read dvscf (without the bare term) for a few patterns.
  ! The spin index is irrelevant and is kept only for compatibility with davcio_drho.
  ! To get a unique ordering independent of q, the same deltav must be used!
  ! Therefore in input we supply fildvscf0 (the fildvscf calculated, say, at gamma)
  !
  call davcio_drho ( v1,  lrdrho, iudvscf0,       1, -1 )
  call davcio_drho ( v2,  lrdrho, iudvscf0, 3*nat/2, -1 )
  call davcio_drho ( v3,  lrdrho, iudvscf0, 3*nat  , -1 )
  deltav = real ( v1(:,1) + v2(:,1) + v3(:,1) )
  deltav = deltav ** 3.d0
  !
  nset = 1 
  do ibnd = 1, nbnd
    ndeg (ibnd) = 1
    do jbnd = 1, nbnd
       unimat (ibnd, jbnd) = czero
    enddo
    unimat (ibnd, ibnd) = 1.d0 
  enddo
  !
  ! count subsets and their degeneracy
  !
  do ibnd = 2, nbnd
    tdegen = abs( et (ibnd, ik0) - et (ibnd-1, ik0)  ) .lt. eps 
    if (tdegen) then
      ndeg (nset) = ndeg(nset) + 1
    else
      nset = nset + 1
    endif
  enddo
  if (iverbosity.eq.1) write (stdout, *) & 
     ik0, nset, (ndeg (iset) ,iset=1,nset)
  !
  ! -----------------------------------------------------------
  ! determine rotations within each subset 
  ! and copy into the global matrix
  ! -----------------------------------------------------------
  !
  n0 = 0
  do iset = 1, nset
    !
    !
    !  the size of the small rotation
    nsmall = ndeg (iset)
    !
    ! unimat is initialized to the identity matrix,
    ! so when the subset has dimension 1, just do nothing
    !
    if (nsmall.gt.1) then
      !
      ! form the matrix elements of the "perturbation"
      !
      allocate ( u (nsmall, nsmall) )
      u = czero
      do ibnd =1, nsmall
        !
        ! ibnd and jbnd are the indexes of the bands within the subsets
        ! ibndg and jbndg are the global indexes of the bands
        !
        ibndg = ibnd + n0
        !
        aux1(:) = (0.d0, 0.d0)
        do ig = 1, nng
          aux1 (nls (igktmp (ig) ) ) = evtmp (ig, ibndg)
        enddo
        call cft3s (aux1, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, + 2)
        do ir = 1, nrxxs
          aux1 (ir) = aux1 (ir) * deltav (ir) 
        enddo
        call cft3s (aux1, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, - 2)
        deltavpsi (1:nng) = aux1 (nls (igktmp (1:nng) ) )
        ! 
        do jbnd = 1, nsmall
          !
          jbndg = jbnd + n0
          ! 
          u (ibnd, jbnd) = ZDOTC (nng, deltavpsi, 1, evtmp(:, jbndg), 1)
        enddo
        !
      enddo
      !
      ! ok I veryfied that when deltav(ir)=1, u is the unity matrix (othonormality)
      !
#ifdef __PARA
      call reduce (2 * nbnd * nbnd, u)
#endif
      !
      ! check hermiticity
      ! 
      do ibnd = 1, nsmall
       do jbnd = 1, nsmall
         if ( abs( conjg (u (jbnd, ibnd)) - u (ibnd, jbnd) ) .gt. eps ) then
           do mbnd = 1,nsmall
             write(stdout,'(10f15.10)') (u (pbnd, mbnd), pbnd=1,nsmall)
           enddo
             call errore ('setphases','perturbation matrix non hermitian',1)
         endif
       enddo
      enddo
      !
      ! now diagonalize the "perturbation" matrix within the degenerate subset
      !
      allocate ( ifail( nsmall ), iwork( 5*nsmall ), w( nsmall ), rwork( 7*nsmall ),&
                 up( nsmall*(nsmall+1)/2 ), cwork( 2*nsmall ), cz( nsmall, nsmall) )
      !
      !  packed upper triangular part for zhpevx
      do jbnd = 1, nsmall
         do ibnd = 1, nsmall
            up (ibnd + (jbnd - 1) * jbnd/2 ) = u ( ibnd, jbnd)
         enddo
      enddo
      !
      call zhpevx ('V', 'A', 'U', nsmall, up , 0.0, 0.0, &
                   0, 0,-1.0, neig, w, cz, nsmall, cwork, &
                   rwork, iwork, ifail, info)
      if (iverbosity.eq.1) then
         !
         write(stdout, '(5x, "Eigenvalues of fake perturbation: ", 10f9.5)') &
           (w(ibnd),ibnd=1,nsmall)
         write(stdout, *) 
         write(stdout,*) 'before diagonalizing the perturbation:'
         do ibnd = 1, nsmall
           write(stdout,'(10f9.4)') (u (ibnd, jbnd), jbnd=1, nsmall)
         enddo
         !       
         ! recalculate u via similarity transform
         !
         u = matmul ( conjg(transpose( cz)), matmul (u, cz) ) 
         write(stdout,*) 'after diagonalizing the perturbation: via matmul'
         do ibnd = 1, nsmall
           write(stdout,'(10f9.4)') (u (ibnd, jbnd), jbnd=1, nsmall)
         enddo
         write(stdout, '("-----------------"/)') 
         !
      endif
      !
      ! now make sure that all the eigenvalues are nondegenerate
      !
      tdegen = .false.
      do ibnd = 2, nsmall  
        tdegen = tdegen .or. ( abs( w(ibnd) - w(ibnd-1) ) .lt. eps )  
      enddo
      if (tdegen) call errore ('setphases', &
         'eigenvalues of pert matrix degenerate',1) 
      !
      ! ...that they are nonvanishing... 
      !
      do ibnd = 1, nsmall  
        tdegen = tdegen .or. ( abs( w(ibnd) ) .lt. eps )  
      enddo
      if (tdegen) call errore ('setphases', &
         'eigenvalues of pert matrix too small',1)
      !
      ! ...and that they are not too close to 1 (orthonormality...)
      !
      do ibnd = 1, nsmall  
        tdegen = tdegen .or. ( abs( w(ibnd) - 1.d0 ) .lt. eps )  
      enddo
      if (tdegen) call errore ('setphases', &
         'eigenvalues of pert matrix too close to 1',1) 
      !
      ! copy the rotation for the subset into the global rotation
      !
      n1 = n0 + 1
      n2 = n0 + ndeg(iset)
      unimat ( n1:n2, n1:n2) = cz ( 1:nsmall, 1:nsmall)
      !
      deallocate ( ifail, iwork, w, rwork, up, cwork, cz )
      deallocate ( u )
      !
    endif
    !
    ! update global band index
    !
    n0 = n0 + ndeg(iset)
  enddo
  !
  if (iverbosity.eq.1) then
     write(stdout, *) '--- rotations for unitary gauge ----'
     write(stdout,'(i4)') ik0
     write(stdout,'(8f9.5)') (et(ibnd,ik0),ibnd=1,nbnd)
     do ibnd = 1, nbnd
       write(stdout,'(10f9.4)') (unimat (ibnd, jbnd), jbnd=1, nbnd)
     enddo
  endif
  !
  ! -----------------------------------------------------------
  ! now fix the phases and update rotation matrix
  ! -----------------------------------------------------------
  !
  !  rotate the coefficients with the matrix u 
  ! (to set the phase on the evc's *after* the rotation)
  !
  cnew(:,:) = czero
  do ibnd = 1, nbnd
    do jbnd = 1, nbnd
      cnew(1:nnglen, ibnd) = cnew(1:nnglen, ibnd) & 
        + unimat (jbnd, ibnd) * c (1:nnglen, jbnd)
    enddo
  enddo
  !
  ! for every band, find the largest coefficient and determine
  ! the rotation which makes it real and positive 
  !
  do ibnd = 1, nbnd
    ! 
    !  this is to identify the largest c(G) by using *integer* 
    !  comparisons on ndig digits [see remarks at the bottom about this]
    !
    imaxnorm = 0
    do ig = 1, nnglen
      ctmp = cnew (ig, ibnd)
      tmp = conjg ( ctmp ) * ctmp
      itmp = nint (10.d0**float(ndig) * tmp)
      if (itmp.gt.imaxnorm) then
        imaxnorm = itmp
        igmax = ig
      endif
    enddo
    !
    ctmp = cnew (igmax, ibnd)
    tmp = conjg ( ctmp ) * ctmp
    !
    ! ...and the corresponding phase
    !
    if ( abs(tmp) .gt. epsc ) then
      ! note that if x + i * y = rho * cos(theta) + i * rho * sin(theta),
      ! then theta = atan2 ( y, x) (reversed order of x and y!)
      theta = atan2 ( aimag( ctmp ), real ( ctmp ) )
    else
      call errore ('setphases','cnew = 0 for some bands: increase nnglen',1)
    endif
    !  
    if (iverbosity.eq.1) then
      write(stdout, '(3i4,2x,f15.10,2(3x,2f9.5))') ik0, ibnd, igmax, theta/pi*180.d0, &
         ctmp, ctmp * exp ( -ci * theta), exp ( -ci * theta)
    endif
    !
    !  now cancel this phase in the rotation matrix
    !
    unimat (:, ibnd) = unimat (:, ibnd) * exp ( -ci * theta)
    !
  enddo
  !
  if (iverbosity.eq.1) then
     write(stdout, *) '--- rotations including phases  ----'
     do ibnd = 1, nbnd
       write(stdout,'(10f9.4)') (unimat (ibnd, jbnd), jbnd=1, nbnd)
     enddo
  endif
  !
  !  last check: unitarity
  !  (if this test is passed, even with wrong phases the final
  !  results -nonwannierized- should be ok)
  !
  unitcheck = matmul ( conjg( transpose (unimat)), unimat)
  do ibnd = 1, nbnd
    do jbnd = 1, nbnd
      if ( (ibnd.ne.jbnd) .and. ( abs(unitcheck (ibnd, jbnd)) .gt. eps )                  &
     .or.  (ibnd.eq.jbnd) .and. ( abs ( abs(unitcheck (ibnd, jbnd)) - 1.d0 ) .gt. eps ) )&
        call errore ('setphases','final transform not unitary',1)
    enddo
  enddo
  !
  nullify ( evtmp ) 
  nullify ( igktmp ) 
  !
  return
end subroutine setphases
!---------------------------------------------------------------------
!
! NOTA BENE: I truncate the c(G) to a given number of digits in
! order to guarantee that the same norm-ordering of c(G) holds
! for different q-runs, machines, libraries etc:
! run q  = 0 : ig = 3 |c(igmax)| = 0.34.....263  <-- igmax = 3
! run q  = 0 : ig = 4 |c(igmax)| = 0.34.....261
! run q /= 0 : ig = 3 |c(igmax)| = 0.34.....260
! run q /= 0 : ig = 4 |c(igmax)| = 0.34.....265  <-- igmax = 4
! In the situation above, I will have psi(r) with a phase difference
! corresponding to that of c(3) and c(4)...
! (Note that the ordering of the G-vectors is always the same) 
! Mind : ndig should be smaller than machine precision, otherwise the 
! integere comparison is useless
!
!---------------------------------------------------------------------
!