!-----------------------------------------------------------------------
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, intra_pool_comm
  use mp, only : mp_barrier, mp_sum
  use random_numbers, only : set_rndm_seed, rndm
#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
  use control_flags, ONLY : iverbosity
  use pwcom
  use wavefunctions_module, ONLY: evc

!$$
  use cell_base
  ! necessary in generating the fake perturbation
!$$
  !
  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

!$$ parameters used in generating deltav(ir) to break the degeneracy
!$$ We basically generate distributed potential sources
!$$ whose center is random but the magnitude is uniform
!$$ (to make things as simple as possible).
!$$ At least for graphene, this random fake perturbation
!$$ successfully breaks degeneracy with >15% difference among
!$$ the eigenvalues. Moreover, if a given perturbation fails
!$$ to break the degeneracy, it repeats the calculation with a
!$$ new random potential until the degeneracy is finally broken.
!$$ However, until now, no second run was necessary.
  integer :: i,j,k
  ! Dummy variables for do loop.
!@JN  real(kind=DP), EXTERNAL :: rndm
  ! rndm: internal random function offered by pwscf:
  !   rand() returns a real number between 0 and 1.
  integer :: nRandCenter, nRandIdx
  ! nRandCenter: The total number of random potential sources
  ! nRandIdx: dummy variable used in a do loop
  real(kind=DP), allocatable :: daRandCoord(:,:)
  ! daRandCoord(3,nRandCenter): cartesian coordinates of
  !   randomly generated centers of potential.
  !   Values of (1,:), (2,:), and (3,:) components
  !   range from 0~nrx1, 0~nrx2, and 0~nrx3, respectively.
  real(kind=DP) :: dSpatialVar, dDistSq, dStdDev, dNorm1, dNorm2
  ! dSpatialVar: Average distance (in atomic units) between
  !   scattering centers.
  ! dDistSq: For a given pixel in a unitcell, the square of
  !   the distance from a specific potential center
  !   (in units of pixel^2).
  ! dStdDev: The standard deviation in real space of the
  !   exponentially decaying potential in units of pixels.
  !   (this length is set equal to dSpatialVar for obvious reasons)
  ! dNorm1: Overall normalization factor of the random potential.
  ! dNorm2: 1/(2*dStdDev*dStdDev).
!$$


  !$$ Random potential initialization.
  !$$ This step is crucial.
  call set_rndm_seed(1);


  !$$ Here, we set the mean distance between potential centers
  !$$ in a.u. This number should be comparable to a.u. in order
  !$$ to break degeneracies of wavefunctions with spatial variation
  !$$ of the order of a.u.'s.
  dSpatialVar = 0.5

  nRandCenter = int(omega/dSpatialVar**3.0)
  allocate(daRandCoord(3,nRandCenter))

  dStdDev = dSpatialVar / (omega/nrxxs)**(1.0/3.0)

  !$$ Just a normalization factor for a Gaussian.
  dNorm1 = 1.0/(2*3.14*dStdDev*dStdDev)**1.5

  !$$ This is to make the eigenvalues of the order 1.
  !$$ (but not crucial)
  dNorm1 = dNorm1 * nrxxs/(nRandCenter+1)

  dNorm2 = 1.0/(2.0*dStdDev*dStdDev)


  if (iverbosity.eq.1) then
     write(stdout,*) 'nRandCenter, nrxxs = ',nRandCenter,nrxxs
     write(stdout,*) 'dStdDev = ',dStdDev
     write(stdout,*) 'dNorm1, dNorm2 = ',dNorm1,dNorm2
     write(stdout,*)
  endif

#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, *) & 
!$$  if (.true.) write (stdout, *) & 
     ik0, nset, (ndeg (iset) ,iset=1,nset)
  !
  ! -----------------------------------------------------------
  ! determine rotations within each subset 
  ! and copy into the global matrix
  ! -----------------------------------------------------------
  !

  !$$ initialize the perturbation
  deltav = czero

  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) )

      !$$ Allocate matrices for rotation outside the do while loop.
      !$$ This step is necessary to iterate until the degeneracy breaks down.
      allocate ( ifail( nsmall ), iwork( 5*nsmall ), w( nsmall ), rwork( 7*nsmall ),&
                 up( nsmall*(nsmall+1)/2 ), cwork( 2*nsmall ), cz( nsmall, nsmall) )

      tdegen = .true.

      !$$
      !$$ repeat until we find a good u matrix
      !$$
      do while(tdegen)

         do nRandIdx=1,nRandCenter
            daRandCoord(1,nRandIdx) = rndm()*nrx1
            daRandCoord(2,nRandIdx) = rndm()*nrx2
            daRandCoord(3,nRandIdx) = rndm()*nrx3
         enddo

         !$$ set up deltav(ir)
         do i=1,nrx1
            do j=1,nrx2
               do k=1,nrx3
                  ir = i + (j - 1) * nrx1 + (k - 1) * nrx1 * nrx2

                  deltav(ir) = 0.0

                  do nRandIdx=1,nRandCenter
                     dDistSq = (i-daRandCoord(1,nRandIdx))*(i-daRandCoord(1,nRandIdx))+&
                              (j-daRandCoord(2,nRandIdx))*(j-daRandCoord(2,nRandIdx))+&
                              (k-daRandCoord(3,nRandIdx))*(k-daRandCoord(3,nRandIdx))
                     ! Add up contributions from all the randomly
                     ! distributed potential centers.
                     deltav(ir) = deltav(ir) + dNorm1*exp(-dDistSq*dNorm2)
                  enddo
               enddo
            enddo
         enddo
         !$$

         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 mp_sum (u, intra_pool_comm)
!@         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
         !
!$$ allocation is done outside the loop
!$$      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
            if (.true.) 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
            endif
            !       
            ! 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) write(stdout,*) 'eigenvalues of pert matrix degenerate'

!$$ in the following, instead of killing the program when the perturbation
!$$ gives degenerate eigenvalues, it changes the perturbation and repeat
!$$ until the degeneracy is broken. The perturbation should not be the same
!$$ for wavefunctions with different wavevector or band indices.
!$$
!$$      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) 

      enddo !$$ repeat until finding a perturbation that gives non-degenerate eigenvalues
      !
      ! 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
  !$$if (.true.) 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
!$$    if (.true.) 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
!$$  if (.true.) 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
!
!---------------------------------------------------------------------
!