!
  !---------------------------------------------------------------------------------
  subroutine hambloch2wan ( nbnd, nbndsub, nks, nksq, nkbl, lgamma, et, xk, cu, &
     nrr, irvec, wslen, chw )
  !---------------------------------------------------------------------------------
  !
  !  From the Hamiltonian in Bloch representationi (coarse mesh), 
  !  find the corresponding Hamiltonian in Wannier representation 
  !
  !  input  : 
  !
  !  output : 
  !
  !  Feliciano Giustino, UCB
  !
  !--------------------------------------------------------------------------------
  !
#include "f_defs.h"
  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
  !
  integer :: nbnd, nbndsub, nks, nksq, nkbl, nrr, irvec (3, nrr)
  ! number of bands 
  ! number of bands in the optimal subspace 
  ! number of kpoints
  ! number of kpoint blocks, in the pool
  ! number of kpoint blocks, total 
  ! number of WS points and coordinates
  logical :: lgamma
  ! true if we are working with q=0
  real(kind=DP) :: et (nbnd, nks), xk (3, nks), wslen (nrr) 
  ! hamiltonian eigenvalues, coarse mesh
  ! kpoint coordinates (cartesian in units of 2piba)
  ! WS vectors length (alat units)
  complex(kind=DP) :: cu (nbnd, nbndsub, nks)
  ! rotation matrix from wannier code
  !
  !  output variables
  !
  complex(kind=DP) :: chw (nbndsub, nbndsub, nrr)  
  !
  ! work variables 
  !
  complex(kind=DP) :: chs(nbndsub, nbndsub, nks)
  ! Hamiltonian in smooth Bloch basis, coarse mesh 
  real(kind=DP), parameter :: bohr2ang = 0.5291772108, twopi = 6.28318530717959
  complex(kind=DP), parameter :: ci = (0.d0,1.d0), czero = (0.d0, 0.d0)
  integer :: ik, ikk, ikq, ibnd, jbnd, imode0, ipol, nsize, rest, &
             ik0, ind, ix, jx, ir, mbnd, pbnd, i
  real(kind=DP) :: rdotk, tmp
  complex(kind=DP) :: cfac, ctmp
  !
  !----------------------------------------------------------
  !    STEP 1: rotation to optimally smooth Bloch states
  !----------------------------------------------------------
  !
  call start_clock ( 'Ham: step 1' )
  !
  !  H~ (k) = U(k)^\dagger * H(k) * U(k)
  !  H~ (k) is chs( nbndsub, nbndsub, ik )
  !
  do ik = 1, nksq
   !
   if (lgamma) then
      ikk = ik
   else
      ikk = 2 * ik - 1
   endif
   !
   do jbnd = 1, nbndsub
    do ibnd = 1, jbnd
       !
       ctmp = czero
       !
       do mbnd = 1, nbnd
         ctmp = ctmp + conjg(cu (mbnd,ibnd,ikk)) * et (mbnd,ikk) * cu (mbnd,jbnd,ikk)
       enddo
       !
       chs (ibnd , jbnd , ikk) = ctmp 
       chs (jbnd , ibnd , ikk) = conjg(ctmp)
       !
    enddo
   enddo
  enddo
  !
  call stop_clock ( 'Ham: step 1' )
  !
  !----------------------------------------------------------
  !  STEP 2: Fourier transform to go into Wannier basis
  !----------------------------------------------------------
  !
  !  H (R) = (1/nk) sum_k e^{-ikR} H~ (k)
  !  chw (nbndsub, nbndsub, ir) is H (R)
  !
  call start_clock ( 'Ham: step 2' )
  !
  ! bring xk in crystal coordinates
  !
  call cryst_to_cart (nks, xk, at, -1)
  !
  chw ( :, :, :) = czero 
  ! 
  do ir = 1, nrr
    !
    do ik = 1, nksq
       !
       if (lgamma) then
         ikk = ik
       else
         ikk = 2 * ik - 1
       endif
       !
       rdotk = twopi * dot_product( xk ( :, ikk), float(irvec( :, ir) ))
       cfac = exp( -ci*rdotk ) / float(nkbl)
       chw ( :, :, ir ) = chw ( :, :, ir ) + cfac * chs ( :, :, ikk )
       !
    enddo
    !
  enddo
#ifdef __PARA
  call poolreduce (2 * nrr * nbndsub * nbndsub , chw )
#endif
  !
  ! bring xk back into cart coord
  !
  call cryst_to_cart (nks, xk, bg, 1)
  !
!@  if (iverbosity.eq.1) then
    !
    !  check spatial decay of Hamiltonian in Wannier basis
    !  the unit in r-space is angstrom, and I am plotting
    !  the matrix for the first mode only
    !
#ifdef __PARA
    if (mpime.eq.ionode_id) then
#endif
!@      write(300, '(/3x,a/)') 'Spatial decay of Hamiltonian in Wannier basis'
      do ir = 1, nrr
        !
        tmp =  maxval ( abs( chw (:,:,ir)) )
        write(300, *) wslen(ir) * celldm (1) * bohr2ang, tmp
        !
      enddo
#ifdef __PARA
    endif
    call mp_barrier()
#endif
!@  endif
  !
  call stop_clock ( 'Ham: step 2' )
  !
  return
  end subroutine hambloch2wan
  !-----------------------------------------------------