!
  ! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino 
  !                                                                           
  ! This file is distributed under the terms of the GNU General Public       
  ! License. See the file `LICENSE' in the root directory of the          
  ! present distribution, or http://www.gnu.org/copyleft.gpl.txt .        
  !
  !-----------------------------------------------------------------------
  SUBROUTINE eliashberg_a2f( lambda, lambda_v )
  !-----------------------------------------------------------------------
  !
  !  Compute the Eliasberg spectral function
  !  in the Migdal approximation. 
  !  
  !  If the q-points are not on a uniform grid (i.e. a line)
  !  the function will not be correct
  !  
  !  02/2009 works in serial on ionode at the moment.  can be parallelized
  !  03/2009 added transport spectral function -- this involves a v_k dot v_kq term 
  !          in the quantities coming from selfen_phon.f90.  Not fully implemented  
  !  10/2009 the code is transitioning to 'on-the-fly' phonon selfenergies
  !          and this routine is not currently functional
  !  
  !-----------------------------------------------------------------------
  !
#include "f_defs.h"
  USE kinds,     ONLY : DP
  USE phcom,     ONLY : nmodes
  USE el_phon,   ONLY : nxqf, wf
#ifdef __PARA
  USE mp_global, ONLY : mpime
#endif
  USE io_global, ONLY : stdout, ionode_id
  USE io_files,  ONLY : find_free_unit, prefix
  implicit none
  !
  real(kind=DP), PARAMETER :: eps = 20.d0 / 13.6058d0 / 8065.5d0, &
     ryd2mev = 13605.8, two = 2.d0, zero = 0.d0,  &
     pi = 3.14159265358979, twopi = 6.28318530717958
  complex(kind = 8), PARAMETER :: ci = (0.d0, 1.d0), cone = (1.d0, 0.d0)
  integer, PARAMETER :: nsmear = 10
  integer :: imode, iq, iw, nw, iua2ffil,iua2ftrfil, ismear
  real(kind=DP) :: lambda_v(nmodes, nxqf)
  real(kind=DP)    :: lambda(nmodes, nxqf), iomega, sigma, a2F_tmp, a2F_tr_tmp, &
       om_max, dw, w0, l, l_tr, eps_om, lambda_tot
  real(kind=DP), allocatable :: a2F(:,:), a2F_tr(:,:), l_a2F(:)
  !
  iua2ffil   = find_free_unit()
  iua2ftrfil = find_free_unit()
  !
  CALL start_clock('a2F')
#ifdef __PARA
  IF (mpime.eq.ionode_id) &
#endif
  !
  OPEN (unit = iua2ffil, file = TRIM(prefix)//".a2f", form = 'formatted')
  OPEN (unit = iua2ftrfil, file = TRIM(prefix)//".a2f_tr", form = 'formatted')
  !
  !
  WRITE(stdout,'(/5x,a)') REPEAT('=',67)
  WRITE(stdout,'(5x,"Eliashberg Spectral Function in the Migdal Approximation")') 
  WRITE(stdout,'(5x,a/)') REPEAT('=',67)
  !
  !
  ! size of the grid - eventually we will specify this 
  nw = 500
  lambda_tot = 0.d0
  IF (nw .lt. 10)  nw = 500
  ! smearing on delta functions
  ALLOCATE( a2F(nw+1, nsmear) )
  ALLOCATE( a2F_tr(nw+1, nsmear) )
  !
  ALLOCATE( l_a2F( nsmear) )
  l_a2F(:) = zero
  !
  ! largest frequency
  !
  om_max = MAXVAL ( wf(:,:) ) + 5.d0/ryd2mev
  dw = om_max/float(nw)
  !
  ! value for which below this we ignore contributions to a2F - negative frequencies give strange answers
  !
  eps_om = 1.d0 / ryd2mev
  !
  !
  ! the smearing parameters will need to be specified in the future
  DO ismear = 1, nsmear
     sigma  = 0.05/ryd2mev + (ismear-1) * 0.15 / ryd2mev
!     sigma  = 0.25/ryd2mev + (ismear-1) * 1.0 / ryd2mev
     !
     a2F(:,ismear) = zero
     a2F_tr(:,ismear) = zero
     !
     !
     DO iw = 1, nw+1  ! loop over points on the a2F(w) graph
        !
        iomega = (float(iw)-1.d0) * dw ! step through the frequncies we wish to plot
        !
        DO iq = 1, nxqf ! loop over q-points 
           DO imode = 1, nmodes ! loop over modes
              !
              w0 = wf(imode,iq)
              l  = lambda(imode,iq)
              l_tr = lambda_v(imode, iq) 
              IF (lambda(imode, iq) .lt. 0.d0)  l = 0.d0 ! sanity check
              IF (lambda_v(imode, iq) .lt. 0.d0)  l_tr = 0.d0 ! sanity check
              IF (wf(imode, iq) .lt. eps_om ) w0 = 0.d0
              IF (lambda_v(imode,iq) .lt. 0.d0 ) lambda_v(imode,iq) = 0.d0
              !
              ! a2F(w) = lambda*omega/(2pi sigma) * exp (omega-w)/2sigma^2
              !
              a2F_tmp    =    w0 * l / (two * sigma * nxqf * twopi**0.5)
              a2F_tr_tmp = w0 * l_tr / (two * sigma * nxqf * twopi**0.5)
              IF (a2F_tr_tmp .lt. 0.d0 ) CALL errore &
                   ('a2F','something wrong with a2F_tr ' ,-1)
              !
              !           write (stdout, *) "a2F(iw) ", a2F_tmp * exp(-0.5*(iomega - w0)**2/(sigma/ryd2mev)**2)
              a2F   (iw,ismear) = a2F   (iw,ismear) + a2F_tmp    * exp(-0.5*(iomega - w0)**2/(sigma)**2)
              a2F_tr(iw,ismear) = a2F_tr(iw,ismear) + a2F_tr_tmp * exp(-0.5*(iomega - w0)**2/(sigma)**2)
              !
           ENDDO !imodes
        ENDDO ! iq
        !
        !  output a2F
        !
        IF (ismear .eq. nsmear) WRITE (iua2ffil,   '(f12.6, 15f12.6)') iomega*ryd2mev, a2F(iw,:)
        IF (ismear .eq. nsmear) WRITE (iua2ftrfil, '(f12.6, 15f12.6)') iomega*ryd2mev, a2F_tr(iw,:)
        !
        ! do the integral 2 int (a2F(w)/w dw)
        !
        IF (iomega .gt. zero) l_a2F(ismear) = l_a2F(ismear) + 2 * a2F(iw,ismear)/iomega * dw 
        !
        ENDDO ! iw
     !
  ENDDO ! nsmear
  !
WRITE(iua2ffil,*) " #   Integrated el-ph coupling"
WRITE(iua2ffil,'("  #    ", 15f12.6)') l_a2F(:) 
!
close(iua2ffil)
close(iua2ftrfil)
  !
DO iq = 1, nxqf ! loop over q-points 
   DO imode = 1, nmodes ! loop over modes
      ! for frequencies > 8 meV
      IF (lambda(imode,iq) .gt. 0.d0 .and. wf(imode,iq)*ryd2mev .gt. 8) lambda_tot = lambda_tot + lambda(imode,iq)
   ENDDO
ENDDO
WRITE (6,'(5x,a,f8.3)') "lambda : ", lambda_tot/float(nxqf)
WRITE(6,*)
!
!
!
#ifdef _PARA
  ENDIF
  CALL mp_barrier
#endif
  !
  CALL stop_clock('a2F')
  CALL print_clock('a2F')
  !
  !-----------------------------------------------------------------------
  END SUBROUTINE eliashberg_a2f
  !-----------------------------------------------------------------------