!
  !---------------------------------------------------------------------------------
  subroutine ephbloch2wanp ( nbnd, nmodes, xk, nq, dynq, irvec, wslen, &
    nrk, nrr, epmatwe, epmatwp )
  !---------------------------------------------------------------------------------
  !
  !  From the EP Matrix in Electron Bloch representation (coarse mesh), 
  !  find the corresponding matrix in Phonon Wannier representation 
  !
  !  input  : 
  !
  !  output : 
  !
  !  Feliciano Giustino, UCB
  !
  !--------------------------------------------------------------------------------
  !
#include "f_defs.h"
  use io_global, only : stdout
  USE kinds, only : DP
  use pwcom, only : at, bg, celldm
  use control_flags, ONLY : iverbosity
#ifdef __PARA
  use para
  USE io_global, ONLY : ionode_id
  USE mp_global, ONLY : mpime
  USE mp, ONLY : mp_barrier
#endif
  implicit none
  !
  !  input variables - note irvec is dimensioned with nrr_k 
  !                    (which is assumed to be larger than nrr_q)
  !
  integer :: nbnd, nrk, nmodes, nq, nrr, irvec (3, nrk)
  ! number of electronic bands
  ! number of electronic WS points
  ! number of branches
  ! number of qpoints
  ! number of WS points and coordinates
  complex(kind=DP) :: dynq (nmodes, nmodes, nq), epmatwe (nbnd, nbnd, nrk, nmodes, nq)
  ! dynamical matrix in bloch representation (Cartesian coordinates)
  ! EP matrix in electron-wannier representation and phonon bloch representation
  !   (Cartesian coordinates)
  real(kind=DP) :: xk (3, nq), wslen (nrr) 
  ! kpoint coordinates (cartesian in units of 2piba)
  ! WS vectors length (alat units)
  !
  !  output variables
  !
  complex(kind=DP) :: epmatwp (nbnd, nbnd, nrk, nmodes, nrr)
  ! EP matrix in electron-wannier representation and phonon-wannier representation
  !
  !
  ! variables for lapack ZHPEVX
  !
  integer :: neig, info, ifail( nmodes ), iwork( 5*nmodes )
  real(kind=DP) :: w2( nmodes ), rwork( 7*nmodes )
  complex(kind=DP) :: champ( nmodes*(nmodes+1)/2 ), &
    cwork( 2*nmodes ), u( nmodes, nmodes)
  !
  ! work variables 
  !
  real(kind=DP), parameter :: bohr2ang = 0.5291772108, twopi = 6.28318530717959
  complex(kind=DP), parameter :: ci = (0.d0,1.d0), czero = (0.d0, 0.d0), cone = (1.d0,0.d0)
  integer :: ik, ikk, ikq, ibnd, jbnd, imode0, ipol, nsize, rest, &
             ik0, ind, ix, jx, ir, mbnd, pbnd, i, nu, ire, irp
  real(kind=DP) :: rdotk, tmp, dist, rvec1(3), rvec2(3), len1, len2
  complex(kind=DP) :: cfac, ctmp
  real(kind=DP), parameter :: rydcm1 = 13.6058d0 * 8065.5d0
  !

  !
  !----------------------------------------------------------
  !  Fourier transform to go into Wannier basis
  !----------------------------------------------------------
  !
  !  D (R) = (1/nk) sum_k e^{-ikR} D (k)
  !
  ! bring xk in crystal coordinates
  !
  call cryst_to_cart (nq, xk, at, -1)
  !
  epmatwp (:,:,:,:,:) = czero
  ! 
  do ir = 1, nrr
    !
    do ik = 1, nq
       !
       rdotk = twopi * dot_product( xk ( :, ik), float(irvec( :, ir) ))
       cfac = exp( -ci*rdotk ) / float(nq)
       epmatwp (:, :, :, :, ir) = epmatwp (:, :, :, :, ir) + cfac * epmatwe(:, :, :, :, ik) 
       !
    enddo
    !
  enddo
  !
  ! bring xk back into cart coord
  !
  call cryst_to_cart (nq, xk, bg, 1)
  !
!@  if (iverbosity.eq.1) then
    !
    !  check spatial decay of dynamical matrix in Wannier basis
    !  the unit in r-space is angstrom
    !
    !  we plot: 
    !
    !  R_e, R_p, max_{m,n,nu} |g(m,n,nu;R_e,R_p)|
    !
#ifdef __PARA
    if (mpime.eq.ionode_id) then
#endif
      open(unit=303,file='303.dat',status='unknown')
!@    write(303, '(/3x,a/)') 'Spatial decay of EP matrix in phonon Wannier basis'
      do ire = 1, nrk
       do irp = 1, nrr
        !
        rvec1 = float(irvec(1,ire))*at(:,1) + &
                float(irvec(2,ire))*at(:,2) + &
                float(irvec(3,ire))*at(:,3) 
        rvec2 = float(irvec(1,irp))*at(:,1) + &
                float(irvec(2,irp))*at(:,2) + &
                float(irvec(3,irp))*at(:,3) 
        len1 = sqrt(rvec1(1)**2.d0+rvec1(2)**2.d0+rvec1(3)**2.d0)
        len2 = sqrt(rvec2(1)**2.d0+rvec2(2)**2.d0+rvec2(3)**2.d0)
        tmp =  maxval ( abs( epmatwp (:, :, ire, :, irp) ) )
        !
        ! rvec1 : electron-electron0 distance
        ! rvec2 : phonon - electron0 distance 
        !
        write(303, '(5f15.10)') len1 * celldm (1) * bohr2ang, &
                                len2 * celldm (1) * bohr2ang, tmp
        !
       enddo
      enddo
      close(303) 
#ifdef __PARA
    endif
    call mp_barrier()
#endif
!@  endif
  !
  return
  end subroutine ephbloch2wanp
  !-----------------------------------------------------