!
  !----------------------------------------------------------------
  program epm
  !----------------------------------------------------------------
  ! 
  ! test-case for the GW-DFPT proposal 
  ! empirical pseudopotential method implementation
  ! for nonpolar tetrahedral semiconductors 
  !
  ! Fri May 30, 4.26 PM: 
  !  great! it gives exactly the same Si bandstructure as in
  !  Cohen&Bergstresser!
  !
  ! NOTE: I use only one G-grid. The reason is that, while V(G-G')
  ! has in principle components 2G_max, in the EPS scheme the  
  ! largest value happens for G2=11, therefore we can simply use
  ! the smooth grid...
  ! To perform a true scaling test we need to include the double
  ! grid (and count the operations on the smooth and on the dense
  ! grids).
  !
  ! NOTE2: the k-dependent Hamiltonian is real and symmetric 
  !  in reciprocal space
  !  (but not in real space because of the momentum component)
  !  therefore we can find eigenstates with c(G) all real
  !  (the eigenvectors of a rs matrix form a orthogonal matrix)
  !  of course psi(r) does not need to be real
  !
  ! Note that the Hamiltonian in G-space is real for the
  ! diamond structure
  !
  !----------------------------------------------------------------
  !
  use gspace
  use parameters
  use constants
  implicit none
  !
  ! variables
  !
  real(dbl) :: gcutm, arg
  real(dbl) :: xk (3,nk), deltag(3), kplusg(3), xq(3,nq)
  integer :: ig, ik, ig1, ig2, ideltag, notfound, ierr, ibnd, iq
  logical :: allowed, equiv
  real(dbl), allocatable :: t(:,:), ss(:), vs(:), h(:,:,:)
  complex(dbl), allocatable :: vr(:)
  real(kind=DP), parameter :: eps8 = 1.0D-8
  real(kind=DP), allocatable :: eval(:), u(:,:), fv1(:), fv2(:), hk(:,:)
  real(kind=DP), allocatable :: psi(:,:), e(:)
  !
  ! temporary variables - debug
  integer :: i, j, k, ir, ir1, ir2
  real(kind=DP), allocatable :: psincr(:), hpsincr(:)
  complex(kind=DP), allocatable :: psincc(:)
  !

  !
  !----------------------------------------------------------------
  ! DEFINE THE CRYSTAL LATTICE
  !----------------------------------------------------------------
  !
  !  direct lattice vectors (cart. coord. in units of a_0)
  !  [Yu and Cardona, pag 23]
  !
  at( :, 1) = (/ 0.0,  0.5,  0.5 /)  ! a1
  at( :, 2) = (/ 0.5,  0.0,  0.5 /)  ! a2
  at( :, 3) = (/ 0.5,  0.5,  0.0 /)  ! a3
  !     
  !  reciprocal lattice vectors (cart. coord. in units 2 pi/a_0)
  !
  bg( :, 1) = (/ -1.0,  1.0,  1.0 /)  ! b1
  bg( :, 2) = (/  1.0, -1.0,  1.0 /)  ! b2
  bg( :, 3) = (/  1.0,  1.0, -1.0 /)  ! b3
  !
  !  atomic coordinates (cart. coord. in units of a_0)
  !
  tau( :, 1) = (/  0.125,  0.125,  0.125 /)  
  tau( :, 2) = (/ -0.125, -0.125, -0.125 /)  
  !
  !----------------------------------------------------------------
  ! GENERATE THE G-VECTORS
  !----------------------------------------------------------------
  !
  ! the G^2 cutoff in units of 2pi/a_0
  ! Note that in Ry units the kinetic energy is G^2, not G^2/2
  ! (note for the Hamiltonian we need to double the size, 2Gmax, hence the factor 4)
  !
  gcutm = four * ecutwfc / tpiba2      
! gcutm = ecutwfc / tpiba2      
  !
  ! set the fft grid
  !
  ! estimate nr1 and check if it is an allowed value for FFT
  !
  nr1 = 1 + int (2 * sqrt (gcutm) * sqrt( at(1,1)**2 + at(2,1)**2 + at(3,1)**2 ) ) 
  nr2 = 1 + int (2 * sqrt (gcutm) * sqrt( at(1,2)**2 + at(2,2)**2 + at(3,2)**2 ) ) 
  nr3 = 1 + int (2 * sqrt (gcutm) * sqrt( at(1,3)**2 + at(2,3)**2 + at(3,3)**2 ) ) 
  !
  do while (.not.allowed(nr1)) 
    nr1 = nr1 + 1
  enddo
  do while (.not.allowed(nr2)) 
    nr2 = nr2 + 1
  enddo
  do while (.not.allowed(nr3)) 
    nr3 = nr3 + 1
  enddo
  !
  call ggen( gcutm)
!@@
goto 10
!@@
  !
  !----------------------------------------------------------------
  ! GENERATE THE k-POINTS
  !----------------------------------------------------------------
  !
  ! brutally defined by hand - here I should read from file
  !
  do ik = 1, nk
    read (5,*) xk(1,ik), xk(2,ik), xk(3,ik) 
  enddo
  !

  !----------------------------------------------------------------
  ! CONSTRUCT THE HAMILTONIAN
  !----------------------------------------------------------------
  !
  allocate ( t(ngm, nk), ss(ngm), vs(ngm) )
  !
  ! the structure factor
  !
  do ig = 1, ngm
    arg = twopi * ( g(1,ig) * tau( 1, 1) + g(2,ig) * tau( 2, 1) + g(3,ig) * tau( 3, 1) )
    ss (ig) = cos ( arg )
  enddo
  !
  ! the empirical pseudopotential
  !
  vs = zero
  ! integer comparison - careful with other structures
  do ig = 1, ngm
    if     ( int ( gl(igtongl(ig)) ) .eq.  3 ) then
      vs (ig) =  v3
    elseif ( int ( gl(igtongl(ig)) ) .eq.  8 ) then
      vs (ig) =  v8
    elseif ( int ( gl(igtongl(ig)) ) .eq. 11 ) then
      vs (ig) =  v11
    endif
  enddo
  !
  ! the k-dependent kinetic energy in Ry
  ! [Eq. (14) of Ihm,Zunger,Cohen J Phys C 12, 4409 (1979)]
  !
  do ik = 1, nk
    do ig = 1, ngm
      kplusg = xk(:,ik) + g(:,ig)
      t ( ig, ik) = tpiba2 * dot_product ( kplusg, kplusg )
    enddo
  enddo
  !
  allocate ( h (ngm, ngm, nk) )
  do ig = 1, ngm
    do ik = 1, nk
        h ( ig, ig, ik) = t ( ig, ik)
    enddo
  enddo
  !
  do ig1 = 1, ngm
   do ig2 = ig1+1, ngm
     !
     ! define ideltag
     ! 
     deltag = g(:,ig1) - g(:,ig2)
     ideltag = 1 
     do while ( .not. equiv ( deltag, g(:,ideltag) ) .and. ideltag.lt.ngm )
       ideltag = ideltag + 1
     enddo   
     !
     ! the free-electron dispersions look ok (vs=0)
     ! when compared to Yu/Cardona Fig. 2.8
     ! 
     do ik = 1, nk
       if (ideltag.ne.ngm) &
          h ( ig1, ig2, ik) = ss (ideltag) * vs (ideltag) 
          h ( ig2, ig1, ik) = h ( ig1, ig2, ik) 
     enddo
     !
   enddo
  enddo
  !
  ! the empirical pseudopotential in real space 
  ! for further use in h_psi
  !
  allocate ( vr(nr) )
  vr = czero
  do ig = 1, ngm
    vr ( nl ( ig ) ) = dcmplx ( ss (ig) * vs (ig), zero )
  enddo
  call cfft3 ( vr, nr1, nr2, nr3,  1)
  !
  deallocate ( ss, vs)
  !
! !----------------------------------------------------------------
! ! BRUTE-FORCE DIAGONALIZATION
! !----------------------------------------------------------------
! !
! allocate ( eval(ngm), hk(ngm,ngm), u(ngm,ngm), fv1(ngm), fv2(ngm) )
! !
! do ik = 1, nk
!   !
!   hk = h(:,:,ik)
!   !
!   ! 0 = no eigenvectors
!   call rs ( ngm, ngm, hk, eval, 0, u, fv1, fv2, ierr)
!   !
!   write(10,'(i3,10(3x,f13.5))') ik, eval(1:nbnd)*ryd2ev
!   !
! enddo
! !
! deallocate ( eval, hk, u, fv1, fv2 )
  !
  !----------------------------------------------------------------
  ! CONJUGATE GRADIENTS AND GREEN'S FUNCTION FROM RECURSION METHOD
  !----------------------------------------------------------------
  !
  ! I checked that the band-structure from CG is identical to the
  ! one obtained by a full diagonalization -> FFT and H|psi> are ok.
  !
  allocate ( e (nbnd), psi (ngm, nbnd) )
  !
  write(6,'(/4x,a)') repeat('-',67)
  write(6,'(4x,"Eigenvalues (eV) from preconditioned conjugate gradient")') 
  write(6,'(4x,a/)') repeat('-',67)
  !
  allocate ( eval(ngm0), hk(ngm0,ngm0), u(ngm0,ngm0), fv1(ngm0), fv2(ngm0) )
  !
  do ik = 1, nk  
    !
    !----------------------------------------------------------------
    ! CONJUGATE GRADIENTS DIAGONALIZATION
    !----------------------------------------------------------------
    !
    ! trial wavefunctions for CG algorithm, cutoff ecut0 (modules.f90)
    !
    hk = h(1:ngm0,1:ngm0,ik)
    call rs ( ngm0, ngm0, hk, eval, 1, u, fv1, fv2, ierr)
    psi = zero
    do ibnd = 1, nbnd
     do ig = 1, ngm0
        psi(ig,ibnd) = u(ig,ibnd)
     enddo
    enddo
    !
    call cgdiag (nbnd, psi, e, t(1,ik), vr )
    !
    write(6,'(4x,i3,8(2x,f6.3))') ik, e(1:8)*ryd2ev
    write(6,'(4x,a/)') repeat('-',67)
    !
    !----------------------------------------------------------------
    ! GREEN'S FUNCTION G(r,r',k,w) IN THE PSINC BASIS 
    !----------------------------------------------------------------
    !
    ! I put this in the CG loop since I want to compare the GF
    ! calculated by expanding on empty states with the one obtained
    ! by the recursion method
    !
!@    do ir1 = 1, nr
!@      do ir2 = 1, nr
    do ir1 = 134, 134 
      do ir2 = 1000, 1000
       !
!!     call green ( ir1, ir2, t(1,ik), vr )
       call green ( ir1, ir2, t(1,ik), vr, psi, e ) ! @@@ debug
       !
     enddo
    enddo
    !
  enddo
!@@
10 continue
!@@
  !
  !----------------------------------------------------------------
  ! SCREENED COULOMB INTERACTION
  !----------------------------------------------------------------
  !
  !  q point for test purposes - really need a loop here
  !
  xq(:, 1) = (/  0.001, 0.0, 0.0/)
  xq(:, 2) = (/  0.010, 0.0, 0.0/)
  xq(:, 3) = (/  0.050, 0.0, 0.0/)
  xq(:, 4) = (/  0.100, 0.0, 0.0/)
  xq(:, 5) = (/  0.150, 0.0, 0.0/)
  xq(:, 6) = (/  0.200, 0.0, 0.0/) 
  xq(:, 7) = (/  0.300, 0.0, 0.0/) 
  xq(:, 8) = (/  0.500, 0.0, 0.0/) 
  xq(:, 9) = (/  0.750, 0.0, 0.0/) 
  xq(:,10) = (/  1.000, 0.0, 0.0/) 
  xq(:,11) = (/  1.500, 0.0, 0.0/) 
  xq(:,12) = (/  2.000, 0.0, 0.0/) 
  xq(:,13) = (/  6.000, 0.0, 0.0/) 
  xq(:,14) = (/ 10.000, 0.0, 0.0/) 
  xq(:,15) = (/ 2.0, 0.0, 0.0/) 
  xq(:,16) = (/ 1.0, 1.0, 1.0/) 
  !
!@  do iq = 1, nq
  do iq = 1, 14
!@@   call coulomb ( xq(:,iq) )
    call dielec_mat ( xq(:,iq) )
  enddo
  !
  !----------------------------------------------------------------
  ! SELF-ENERGY
  !----------------------------------------------------------------
  !
  end
  !