! 
  !--------------------------------------------------------------------------------
  subroutine rotate_eigenm ( iq_first, iq, nqc, isym, nsym, s, invs, irt, &
      rtau, xq, isq, cz1, cz2, timerev ) 
  !--------------------------------------------------------------------------------
  !
  ! Here:
  !
  ! 1) determine the unitary matrix which gives the transformed eigenmodes
  !    from the first q in the star to the present q
  !
  ! 2) rotate eigenmodes from the first q in the star (cz1) to the current q (cz2)
  ! 
  ! In rotate_epmat.f90:
  !
  ! 3) rotate the electron-phonon matrix from the cartesian representation
  !    of the first qpoint of the star to the eigenmode representation (using cz1).
  ! 
  ! 4) rotate the electron-phonon matrix from the eigenmode representation
  !    to the cartesian representation of the qpoint iq in the star (with cz2).
  !
  !    Step 3 requires using the rotated eigenmodes to set the phases 
  !    (the diagonalizers of the rotated dyn mat will not
  !    work because of the gages in deg. subspaces and the phases).
  !
  ! W.r.t. the standard method of q2qstar_ph.f90, here I construct the
  ! unitary matrix Gamma defined in Maradudin & Vosko, RMP 1968.
  ! In this way all rotations are just zgemm operations and we throw away  
  ! all nasty indices.
  !
  ! Feliciano Giustino, UCB Feb 2007
  !
  !--------------------------------------------------------------------------
  !
#include "f_defs.h"
  USE kinds, only : DP
  use io_global, only : stdout
  use el_phon, only : dynq, epmatq
  use phcom, only : nmodes
  use control_flags, only : iverbosity
  use pwcom, only : nbnd, nks, nr1, nr2, nr3, at, bg
  use ions_base, ONLY : tau, nat
  USE ions_base, ONLY : amass, tau, nat, ntyp => nsp, ityp
  implicit none
  !
  integer :: iq_first, nqc, isym, isq(48), nsq, nq, nsym, iq
  !  originating q point of the star
  !  total number of qpoints
  !  fractional translation
  !  the symmetry under consideration
  !  ...
  !  degeneracy of the star of q
  !  number of q in the star
  !  number of global (crystal) symmetries
  !  index of this q in the star
  real(kind=DP) :: xq(3), rtau(3,48,nat)
  !  the rotated q vector
  !  the relative position of the rotated atom to the original one
  integer :: s(3,3,48), irt(48,nat), g0vec(3), invs(48)
  !  the symmetry operation for the eigenmodes
  !  the rotated of each atom
  !  the folding G-vector
  !  the index of the inverse operation
  complex(kind=DP) :: gamma(nmodes, nmodes), cz1( nmodes, nmodes), &
       cz2(nmodes, nmodes)
  !  the Gamma matrix for the symmetry operation on the dyn mat
  !  the eigenvectors for the first q in the star
  !  the rotated eigenvectors, for the current q in the star
  logical :: timerev
  !  true if we are using time reversal

  !
  ! variables for lapack zhpevx
  !
  complex(kind=DP) :: cwork( 2*nmodes ), dynp( nmodes*(nmodes+1)/2 )
  integer :: neig, info, ifail( nmodes ), iwork( 5*nmodes )
  real(kind=DP) :: w1( nmodes ), w2(nmodes), rwork( 7*nmodes )
  !  dynp: complex dynmat packed (upper triangular part for zhpevx)
  !  after hermitian-ization
  !
  complex(kind=DP), parameter :: cone = (1.d0, 0.d0), czero = (0.d0, 0.d0)
  real(kind=DP), parameter :: eps = 1.d-10, zero = 0.d0, &
       twopi = 6.28318530717959, rydcm1 = 13.6058d0 * 8065.5d0
  !
  ! work variables 
  !
  integer :: imode, jmode, nu, mu, ipol, jpol, sna, ism1, na, nb, jsym
  real(kind=DP) :: wtmp( nmodes ), v(3), arg, massfac, scart(3,3)
  complex(kind=DP) :: cfac, dyn1(nmodes, nmodes), dyn2( nmodes, nmodes)

  !
  ! ------------------------------------------------------------------
  ! diagonalize dynq(iq_first) --> w1, cz1
  ! ------------------------------------------------------------------
  !
  !  first divide by the square root of masses 
  !
  do na = 1, nat
   do nb = 1, nat
     massfac = 1.d0 / sqrt ( amass(ityp(na)) * amass(ityp(nb)) )
     dyn1(3*(na-1)+1:3*na, 3*(nb-1)+1:3*nb) = &
     dynq(3*(na-1)+1:3*na, 3*(nb-1)+1:3*nb, iq_first) * massfac
   end do
  end do
  !
  do jmode = 1, nmodes
   do imode = 1, jmode
      dynp (imode + (jmode - 1) * jmode/2 ) = &
      ( dyn1 ( imode, jmode) + conjg ( dyn1 ( jmode, imode) ) ) / 2.d0
   enddo
  enddo
  !
  call zhpevx ('V', 'A', 'U', nmodes, dynp , 0.0, 0.0, &
               0, 0,-1.0, neig, w1, cz1, nmodes, cwork, &
               rwork, iwork, ifail, info)
  !
  if (iverbosity.eq.1) then
    !
    !  check the frequencies
    !
    do nu = 1, nmodes
      if ( w1 (nu) .gt. 0.d0 ) then
         wtmp(nu) =  sqrt(abs( w1 (nu) ))
      else
         wtmp(nu) = -sqrt(abs( w1 (nu) ))
      endif
    enddo
    write ( stdout, '(5x,"Frequencies of the matrix for the first q in the star")') 
    write ( stdout, '(6(2x,f10.5))' ) (wtmp(nu)*rydcm1, nu = 1, nmodes)
    !
  endif
  !
  !  here dyn1 is dynq(:,:,iq_first) after dividing by the masses
  !  (the true dyn matrix D_q)
  !
  !  -----------------------------------------------------------------------
  !  the matrix gamma (Maradudin & Vosko, RMP, eq. 2.37)   
  !  -----------------------------------------------------------------------
  !
  !  I have built the matrix by following step-by-step q2star_ph.f90 and
  !  rotate_and_add_dyn.f90
  !
  ism1 = invs (isym)
  !
  !  the symmetry matrix in cartesian coordinates 
  !  (so that we avoid going back and forth with the dynmat)  
  !  note the presence of both at and bg in the transform!
  !
  scart = float ( s ( :, :, ism1) )
  scart = matmul ( matmul ( bg, scart), transpose (at) )
  ! 
  gamma = czero
  do na = 1, nat
    !
    ! the corresponding of na in {S|v}
    sna = irt (isym, na)
    !
    ! cfac = exp[iSq*(tau_K - {S|v} tau_k)]   (Maradudin&Vosko RMP Eq. 2.33)
    ! [v can be ignored since it cancels out, see endnotes. xq is really Sq]
    ! rtau(:,isym,na) = s(:,:,invs(isym)) * tau(:, na) - tau(:,irt(isym,na))) (cartesian)
    !
    arg = twopi * dot_product (xq, rtau (:, isym, na)) 
    cfac = dcmplx (cos(arg),-sin(arg))
    !
    !  the submatrix (sna,na) contains the rotation scart 
    !
    gamma ( 3*(sna-1)+1:3*sna, 3*(na-1)+1:3*na ) = cfac * scart 
    !
  enddo
  !
  !  possibly run some consistency checks
  !
  if ( iverbosity .eq. 1 ) then
    !
    !  D_{Sq} = gamma * D_q * gamma^\dagger (Maradudin & Vosko, RMP, eq. 3.5) 
    ! 
    call zgemm ('n', 'n', nmodes, nmodes, nmodes , cone  , gamma  , &
           nmodes, dyn1 , nmodes, czero , dyn2, nmodes)
    call zgemm ('n', 'c', nmodes, nmodes, nmodes , cone  , dyn2, &
           nmodes, gamma, nmodes, czero , dyn1   , nmodes)
    !
    do jmode = 1, nmodes
     do imode = 1, jmode
       dynp (imode + (jmode - 1) * jmode/2 ) = &
        ( dyn1 ( imode, jmode) + conjg ( dyn1 ( jmode, imode) ) ) / 2.d0
     enddo
    enddo
    !
    call zhpevx ('V', 'A', 'U', nmodes, dynp , 0.0, 0.0, &
                 0, 0,-1.0, neig, w2, cz2, nmodes, cwork, &
                 rwork, iwork, ifail, info)
    !
    !  check the frequencies
    !
    do nu = 1, nmodes
      if ( w2 (nu) .gt. 0.d0 ) then
         wtmp(nu) =  sqrt(abs( w2 (nu) ))
      else
         wtmp(nu) = -sqrt(abs( w2 (nu) ))
      endif
    enddo
    write ( stdout, '(6(2x,f10.5))' ) (wtmp(nu)*rydcm1, nu = 1, nmodes)
    !
  endif
  !
  !
  !  -----------------------------------------------------------------------
  !  transformation of the eigenvectors: e_{Sq} = gamma * e_q  (MV eq. 2.36)
  !  -----------------------------------------------------------------------
  !
  call zgemm ('n', 'n', nmodes, nmodes, nmodes , cone  , gamma  , &
       nmodes, cz1 , nmodes, czero , cz2, nmodes)

  !
  ! in case time reversal is applied we use Eq. (4.4) Maradudin&Vosko, RMP
  !
  if (timerev) cz2 = conjg ( cz2 )
  call errore('rotate_eigenm','rotate_eigenm: time reversal never tested',0)
  !

  !
  !   now check that cz2 diagonalizes dyn2 (after dividing by the masses)
  !
  do na = 1, nat
   do nb = 1, nat
     massfac = 1.d0 / sqrt ( amass(ityp(na)) * amass(ityp(nb)) )
     dyn2(3*(na-1)+1:3*na, 3*(nb-1)+1:3*nb) = &
       dynq(3*(na-1)+1:3*na, 3*(nb-1)+1:3*nb, nqc) * massfac
   end do
  end do
  !
  ! the rotated matrix and the one read from file
  if (iverbosity.eq.1) write (6,'(2f15.10)') dyn2-dyn1 
  !
  ! here I have checked that the matrix rotated with gamma 
  ! is perfectly equal to the one read from file for this q in the star
  !
  call zgemm ('c', 'n', nmodes, nmodes, nmodes, cone, cz2, &
      nmodes, dyn2, nmodes, czero, dyn1, nmodes)
  call zgemm ('n', 'n', nmodes, nmodes, nmodes, cone, dyn1, &
     nmodes, cz2, nmodes, czero, dyn2, nmodes)
  !
  do nu = 1, nmodes
    w2(nu) = abs(dyn2(nu,nu))
    do mu = 1, nmodes
    if ( mu.ne.nu .and. abs(dyn2(mu,nu)).gt.eps ) call errore &
      ('rotate_eigenm','problem with rotated eigenmodes',0)
    enddo
  enddo
  !
  if (iverbosity.eq.1) then
    !
    !  a simple check on the frequencies
    !
    do nu = 1, nmodes
      if ( w2 (nu) .gt. 0.d0 ) then
         wtmp(nu) =  sqrt(abs( w2 (nu) ))
      else
         wtmp(nu) = -sqrt(abs( w2 (nu) ))
      endif
    enddo
    write ( stdout, '(5x,"Frequencies of the matrix for the current q in the star")') 
    write ( stdout, '(6(2x,f10.5))' ) (wtmp(nu)*rydcm1, nu = 1, nmodes)
    !
  endif
  !
  return
  end subroutine rotate_eigenm