!
  !-----------------------------------------------------------------------
  subroutine green ( ir1, ir2, g2kin, vr, psi, eval)
  !-----------------------------------------------------------------------
  !
  ! Green's function from Haydock's recursion: G(r,r',k,w)
  ! the k-dependence is inside the kinetic energy g2kin 
  !
  ! for some reasons the dot_product(conjg(),()) does not work
  ! need to check this out
  !
  !-----------------------------------------------------------------------
  !
  use gspace
  use constants
  use parameters
  implicit none
  !
  complex(kind=DP) :: vr (nr)
  real(kind=DP) :: g2kin (ngm)
  integer :: ir1, ir2, i, j, k, ir
  !
  ! notation as in Haydock - review (pag 227)
  ! e.g.: bnpunp = b_{n+1} * u_{n+1}, hun = H * u_n
  !
  complex(DP) :: d1 (ngm), d2 (ngm), u (ngm)
  !
  complex(DP) :: psinc (nr), ZDOTC
  real(DP) :: freq(nw), norm
  integer :: ig, iw, ibnd
  real(DP) :: psi (ngm, nbnd), c2(nbnd), eval(nbnd) ! debug
  complex(DP) :: ctmp1, ctmp2, psic(ngm)
  complex(DP) :: gr(nw), gr_xx(nw), gr_exp(nw)
  !
  ! frequency bins for the Green's function
  !
  do iw = 1, nw
    freq(iw) = wmin + ( wmax - wmin ) / float (nw-1) * float(iw-1)
  enddo
  !

  !
  ! DEFINE PSINC'S
  !
  ! first d1 = psinc(r1) ...
  !
  psinc = czero
  psinc (ir1) = cone
  !
  ! go to G-space
  !
  call cfft3 ( psinc, nr1, nr2, nr3, -1)
  do ig = 1, ngm
    d1 (ig) = psinc( nl(ig) )
  enddo
  !
  ! normalize
  !
  norm = DREAL ( ZDOTC (ngm, d1, 1, d1, 1) )
  norm = sqrt(norm)
  call DSCAL (2 * ngm, one / norm, d1, 1)
  !
  ! ... then d2 = psinc(r2) 
  !
  psinc = czero
  psinc (ir2) = cone
  !
  ! go to G-space
  !
  call cfft3 ( psinc, nr1, nr2, nr3, -1)
  do ig = 1, ngm
    d2 (ig) = psinc( nl(ig) )
  enddo
  !
  ! normalize
  !
  norm = DREAL ( ZDOTC (ngm, d2, 1, d2, 1) )
  norm = sqrt(norm)
  call DSCAL (2 * ngm, one / norm, d2, 1)
  !

  !
  ! DEBUG - calculate the Green's function using the sum over states
  ! <d1| G(r,r',k,w) |d2>
  !
  gr_exp = czero
  do ibnd = 1, nbnd
    psic = dcmplx (psi(:,ibnd), zero)
    ctmp1 = ZDOTC (ngm, psic, 1, d1, 1)
    ctmp2 = ZDOTC (ngm, psic, 1, d2, 1)
    gr_exp = gr_exp + conjg(ctmp1)*ctmp2 / ( freq - eval(ibnd)*ryd2ev + ci * delta)
  enddo
  do iw = 1, nw
    write(200,'(3f15.10)') freq(iw), gr_exp(iw) 
  enddo
  !
 
  !
  ! G_12 = G_uu - G_vv - i * (G_ww - G_zz) 
  !
  ! u = 0.5 * (d1+ d2)
  ! v = 0.5 * (d1- d2)
  ! w = 0.5 * (d1+id2)
  ! z = 0.5 * (d1-id2)
  !
  ! NOTE: when ir1=ir2 the procedure at (d1-d2)/2
  ! finds a nihil eigenvalue and the
  ! Green's function is peaked at 0 - spurious
  !
  ! IF I DON'T NORMALIZE IT DOESN NOT WORK - this is due to the fact that I am
  ! using wrong psinc functions - these are not orthogonal!
  !
  gr = czero
  !
  u = 0.5d0 * ( d1 + d2 ) 
  norm = DREAL ( ZDOTC (ngm, u, 1, u, 1) )
  norm = sqrt(norm)
  call DSCAL (2 * ngm, one / norm, u, 1)
  call haydock ( g2kin, vr, freq, u, gr_xx)
  gr_xx = norm * norm * gr_xx
  gr = gr + gr_xx
  !
  if (ir1.ne.ir2) then
    !
    !  trivial solution - spurios peak at 0 in the GF
    !
    u = 0.5d0 * ( d1 - d2 )                    
    norm = DREAL ( ZDOTC (ngm, u, 1, u, 1) )
    norm = sqrt(norm)
    call DSCAL (2 * ngm, one / norm, u, 1)
    call haydock ( g2kin, vr, freq, u, gr_xx)
    gr_xx = norm * norm * gr_xx
    gr = gr - gr_xx
    !
  endif
  !
  u = 0.5d0 * ( d1 + ci * d2 ) 
  norm = DREAL ( ZDOTC (ngm, u, 1, u, 1) )
  norm = sqrt(norm)
  call DSCAL (2 * ngm, one / norm, u, 1)
  call haydock ( g2kin, vr, freq, u, gr_xx)
  gr_xx = norm * norm * gr_xx
  gr = gr - ci * gr_xx
  !
  u = 0.5d0 * ( d1 - ci * d2 )
  norm = DREAL ( ZDOTC (ngm, u, 1, u, 1) )
  norm = sqrt(norm)
  call DSCAL (2 * ngm, one / norm, u, 1)
  call haydock ( g2kin, vr, freq, u, gr_xx)
  gr_xx = norm * norm * gr_xx
  gr = gr + ci * gr_xx
  !
  do iw = 1, nw
    write(201,'(3f15.10)')  freq(iw), gr(iw) 
  enddo
  !
  return
  end subroutine green
  !
  !-----------------------------------------------------------------------
  subroutine haydock ( g2kin, vr, w, un, gr)
  !-----------------------------------------------------------------------
  !
  use gspace
  use constants
  use parameters
  implicit none
  !
  complex(kind=DP) :: vr (nr)
  real(kind=DP) :: g2kin (ngm)
  integer :: i, j, k, ir
  !
  ! notation as in Haydock - review (pag 227)
  ! e.g.: bnpunp = b_{n+1} * u_{n+1}, hun = H * u_n
  complex(DP) :: unp (ngm), un (ngm), unm (ngm), hun (ngm), bnpunp (ngm)
  !
  real(DP) :: w(nw), an, bn, bnp, a(nstep+2), b(nstep+2), bnp2
  complex(DP) :: gr(nw), psinc (nr), ZDOTC
  integer :: ig, istep, iw, nterm 
  complex(DP) :: Z
  real(DP) :: D(nstep+1), E(nstep), WORK
  integer :: INFO, N
  !
  ! variables for the terminator
  !
  real(DP) :: a_term, b2_term, term(nw)
  !

!@@
!@@
!@@ WARNING: FOR SOME REASON I DON'T UNDERSTAND, IF I COMMENT
!@@ ANY OF THE "DEBUG" SECTIONS, THE EXECUTION IS SCREWED UP !!
!@@
!@@

  !
  ! initialize the recursion 
  !
  a = zero
  b = zero
  bn = zero
  !
  do istep = 0, nstep
    !
    ! basic Haydock scheme [Eqs. (5.2)-(5.13) of SSP]
    !
    call h_psi_c ( un, hun, g2kin, vr) 
    an = DREAL ( ZDOTC (ngm, un, 1, hun, 1) )
    !
    ! DEBUG: print out un in real space
    !
    psinc = (0.d0,0.d0)
    do ig = 1, ngm
       psinc ( nl ( ig ) ) = un (ig)
    enddo
    call cfft3 ( psinc, nr1, nr2, nr3,  1)
    do i = 1, nr1
    do j = 1, nr2
    do k = 1, nr3
      ir = i + (j-1) * nr1 + (k-1) * nr1 * nr2
!     write (100+istep,'(3i5,2f9.5)') i, j, k, psinc(ir)
    enddo
    enddo
    enddo
    !
    !  bnpunp = hun - an * un - bn * unm
    !  unp = bunpunp / norm(bnpunp)
    !
    call ZAXPY (ngm, - an, un , 1, hun, 1)
    call ZAXPY (ngm, - bn, unm, 1, hun, 1)
    bnp = DREAL ( ZDOTC (ngm, hun, 1, hun, 1) )
    bnp = sqrt (bnp)
    call DSCAL (2 * ngm, one / bnp, hun, 1)
    call ZCOPY (ngm, hun, 1, unp, 1)
    !
    ! now store and swap
    !
    a (istep+1) = an
    b (istep+2) = bnp 
    !
    call ZCOPY (ngm, un , 1, unm, 1)  ! unm <- un
    call ZCOPY (ngm, unp, 1, un , 1)  ! un  <- unp
    bn = bnp
    !
    if (istep.gt.0) then
     !
     ! DEBUG: eigenvalues of the Hamiltonian at this iteration
     !
     write(6,'(/4x,a)') repeat('-',67)
     write(6,'(4x,"Eigenvalues (eV) of the tridiagonal matrix at iteration",i5)') &
       istep + 1
     write(6,'(4x,a/)') repeat('-',67)
     !
     N = istep + 1
     D = zero
     E = zero
     D(1:N) = a(1:N)
     E(1:N-1) = b(2:N)
     call DSTEV( 'N', N, D, E, Z, 1, WORK, INFO )
     write(6,'(4x,7f9.3)') D(1:N) * ryd2ev
     !
     ! define the terminator using the extrema of the eigenvalue spectrum
     ! (this should be much more stable than trying to fit the fluctuating
     ! coefficients to constant values) - Giannozzi et al, Appl. Num. Math. 
     ! 4, 273 (1988) - Eq. (7.1)
     !
     a_term = 0.5d0 * ( D(1) + D(N) ) * ryd2ev
     b2_term = ( 0.25d0 * ( D(N) - D(1) ) * ryd2ev ) ** 2.d0
     !
    endif
    !  
    ! DEBUG: green's function at each step
    !
    gr = czero
    !
    do i = istep,0,-1
      gr = w + cmplx(zero, delta) - a(i+1) * ryd2ev - ( ryd2ev * b(i+2) ) ** two * gr
      gr = cone / gr
    enddo
    do iw = 1, nw
!     write(100,'(i5,2f15.10)') istep, w(iw), - one/pi * aimag ( gr(iw) )
    enddo
!   write(100,*)
    !
  enddo
  !
  ! Green's function projected on the initial state
  !
  gr = czero
  !
  ! start with a costant chain defined by a_term and b2_term (already in eV)
  !
!
!  @@ THIS IS BAD ! IF I LOOK AT gr GIVEN ONLY BY THE TERMINATOR IT'S A MESS...
!  @@ WHAT THE HELL !!!
!
!  nterm = 500
!  do i = 1, nterm
!    gr = w + cmplx(zero, delta) - a_term - b2_term * gr
!    gr = cone / gr
!  enddo
  !
  do i = nstep,0,-1
    gr = w + cmplx(zero, delta) - a(i+1) * ryd2ev - ( ryd2ev * b(i+2) ) ** two * gr
    gr = cone / gr
  enddo
  !
  ! write down the recursion coefficients
  !  
! write(6,'(/4x,a)') repeat('-',67)
! write(6,'(4x,"Recursion coefficients")') 
! write(6,'(4x,a/)') repeat('-',67)
! do istep = 0, nstep
!   write(6,'(4x,2f15.5)') a(istep+1), b(istep+1)
! enddo
  !
  return
  end subroutine haydock
  !