!
  !-----------------------------------------------------------------------
  subroutine selfen_phon
  !-----------------------------------------------------------------------
  !
  !  compute the imaginary part of the phonon self energy due to electron-
  !  phonon interaction in the Migdal approximation. This corresponds to 
  !  the phonon linewidth (half width). The phonon frequency is taken into
  !  account in the energy selection rule.
  !
  !  Use matrix elements, electronic eigenvalues and phonon frequencies
  !  from ep-wannier interpolation
  !
  !  Feliciano Giustino, University of California at Berkeley, 2006
  !
  !-----------------------------------------------------------------------
  !
#include "f_defs.h"
  USE kinds, only : DP
  USE cell_base, ONLY : at, bg
  USE io_global, ONLY : stdout
  use phcom, only : nbndsub, lrepmatf, iunepmatf, lgamma, nmodes, &
      fsthick, eptemp, ngaussw, degaussw, iuetf, wmin, wmax, nw, nbndskip
  use pwcom, only : nelec, ef
  use el_phon
#ifdef __PARA
  use para
  use mp, only : mp_barrier
#endif
  implicit none
  !
  real(kind=DP), parameter :: eps = 20.d0 / 13.6058d0 / 8065.5d0, &
     ryd2mev = 13605.8, one = 1.d0, ryd2ev = 13.6058,         &
     two = 2.d0, zero = 0.d0, pi = 3.14159265358979
  complex(kind = 8), parameter :: ci = (0.d0, 1.d0), cone = (1.d0, 0.d0)
  integer :: ik, ikk, ikq, ibnd, jbnd, imode, nrec, iq, fermicount
  complex(kind=DP) epf (ibndmax-ibndmin+1, ibndmax-ibndmin+1)
  real(kind=DP) :: g2, ekk, ekq, wq, ef0, wgkk, wgkq, lambda, &
     weight, w0g1, w0g2, wgauss, dosef, dos_ef, gamma(nmodes)
  !

  write(6,'(/5x,a)') repeat('=',67)
  write(6,'(5x,"Phonon (Imaginary) Self-Energy in the Migdal Approximation")') 
  write(6,'(5x,a/)') repeat('=',67)
  !
  if ( fsthick.lt.1.d3 ) &
    write(stdout, '(/5x,a,f10.6,a)' ) &
      'Fermi Surface thickness = ', fsthick, ' Ry'
  write(stdout, '(/5x,a,f10.6,a)' ) &
    'Golden Rule strictly enforced with T = ',eptemp, ' Ry'
  !
  ! Fermi level and corresponding DOS
  !
  ! here we take into account that we may skip bands when we wannierize
  ! (spin-unpolarized)
  if (nbndskip.gt.0) then
     nelec = nelec - two * nbndskip
     write(stdout, '(/5x,"Skipping the first ",i4," bands:")' ) nbndskip
     write(stdout, '(/5x,"The Fermi level will be determined with ",f9.5," electrons")' ) nelec
  endif
  !
  !   Note that the weights of k+q points must be set to zero here
  call efermig (etf, nbndsub, nksf, nelec, wkf, degaussw, ngaussw, ef0)
  dosef = dos_ef (ngaussw, degaussw, ef0, etf, wkf, nksf, nbndsub)
  !   N(Ef) in the equation for lambda is the DOS per spin
  dosef = dosef / two
  !
  write (6, 100) degaussw, ngaussw
  write (6, 101) dosef, ef0 * ryd2ev
  !
  ! loop over all q points of the fine mesh
  ! 
  do iq = 1, nxqf
    !
    fermicount = 0
    gamma = zero
    !
    do ik = 1, nksqf
      !
      if (lgamma) then
        ikk = ik
        ikq = ik
      else
        ikk = 2 * ik - 1
        ikq = ikk + 1
      endif
      !
      ! we read the hamiltonian eigenvalues (those at k+q depend on q!) 
      !
      nrec = (iq-1) * nksf + ikk
      call davcio ( etf (ibndmin:ibndmax, ikk), ibndmax-ibndmin+1, iuetf, nrec, - 1)
      nrec = (iq-1) * nksf + ikq
      call davcio ( etf (ibndmin:ibndmax, ikq), ibndmax-ibndmin+1, iuetf, nrec, - 1)
      !
      ! here we must have ef, not ef0, to be consistent with ephwann_shuffle
      if ( ( minval ( abs(etf (:, ikk) - ef) ) .lt. fsthick ) .and. &
           ( minval ( abs(etf (:, ikq) - ef) ) .lt. fsthick ) ) then
        !
        fermicount = fermicount + 1
        !
        do imode = 1, nmodes
          !
          ! the phonon frequency
          wq = wf (imode, iq)
          !
          !  we read the e-p matrix
          !
          nrec = (iq-1) * nmodes * nksqf + (imode-1) * nksqf + ik
          call dasmio ( epf, ibndmax-ibndmin+1, lrepmatf, iunepmatf, nrec, -1)
          !
          do ibnd = 1, ibndmax-ibndmin+1
            !
            !  the fermi occupation for k
            ekk = etf (ibndmin-1+ibnd, ikk) - ef0
            wgkk = wgauss( -ekk/eptemp, -99)
            !
            do jbnd = 1, ibndmax-ibndmin+1
              !
              !  the fermi occupation for k+q
              ekq = etf (ibndmin-1+jbnd, ikq) - ef0
              wgkq = wgauss( -ekq/eptemp, -99)  
              !
              ! here we take into account the zero-point sqrt(hbar/2M\omega)
              ! with hbar = 1 and M already contained in the eigenmodes
              ! g2 is Ry^2, wkf must already account for the spin factor
              !

! NON FUNZIONA SCAMBIANDO i,j

   !@         g2 = abs(epf (ibnd, jbnd))**two / ( two * wq )
              g2 = abs(epf (jbnd, ibnd))**two / ( two * wq )
              !
              weight = wkf (ikk) * (wgkk - wgkq) * &
                 aimag ( cone / ( ekq - ekk - wq - ci * degaussw ) ) 
              gamma (imode) = gamma (imode) + weight * g2 
              !
            enddo
          enddo
          !
        enddo
        !
        ! endif tfermi 
      endif 
      !
    enddo
#ifdef __PARA
    ! collect contributions from all pools (sum over k-points)
    call poolreduce ( nmodes, gamma)
#endif
    !
    write(6,'(/5x,"iq = ",i5," coord.: ", 3f9.5)') iq, xqf(:,iq) 
    write(6,'(5x,a)') repeat('-',67)
    do imode = 1, nmodes
      ! 
      wq = wf (imode, iq)
      lambda = zero
      if ( sqrt(abs(wq)) .gt. eps ) lambda = gamma(imode) / pi / wq**two / dosef
      !
      write(6, 102) imode, lambda, ryd2mev * gamma(imode), ryd2mev * wq
    enddo
    write(6,'(5x,a/)') repeat('-',67)
    !
    ! test only
#ifdef __PARA
    if (me.eq.1) & 
#endif
    write(900,'(6(2x,f9.5))') (ryd2mev*gamma(imode),imode=1,nmodes)
    !
    ! end loop on q
  enddo
  !
  write( stdout, '(/5x,a,i5,a,i5/)' ) &
    'Number of (k,k+q) pairs on the Fermi surface: ',fermicount, ' out of ', nksqf
  !
100 format(5x,'Gaussian Broadening: ',f7.3,' Ry, ngauss=',i4)
101 format(5x,'DOS =',f10.6,' states/spin/Ry/Unit Cell at Ef=',f10.6,' eV')
102 format(5x,'lambda( ',i3,' )=',f9.3,'   gamma=',f9.3,' meV','   omega=',f9.3,' meV')
  !
  return
  end subroutine selfen_phon