!
  !----------------------------------------------------------------
  program gwhs
  !----------------------------------------------------------------
  ! 
  ! pilot code the GW-HS method
  ! empirical pseudopotential method implementation
  ! for nonpolar tetrahedral semiconductors 
  !
  ! 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).
  !
  !----------------------------------------------------------------
  !
  use gspace
  use parameters
  use constants
  implicit none
  !
  ! variables
  !
  real(dbl) :: gcutm, arg, wp
  real(dbl) :: deltag(3), kmq(3), kplusg(3), xk0(3,nk0)
  integer :: ig, ik, ig1, ig2, ideltag, notfound, ierr, ibnd, iq, nk, ik0
  integer :: iw, iwp, iwmwp, nw_sigma
  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(:), wq(:), xq(:,:), g2kin (:)
  complex(dbl), allocatable :: scrcoul (:,:,:), greenf (:,:,:), sigma(:,:,:)
  !
  ! temporary variables - debug
  integer :: i, j, k
  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)
  !
  !----------------------------------------------------------------
  ! READ THE k-POINTS
  !----------------------------------------------------------------
  !
  ! read from file the k-points for the self-energy \Sigma(k)
  !
  do ik = 1, nk0
    read (5,*) xk0(1,ik), xk0(2,ik), xk0(3,ik) 
  enddo
  !
  !----------------------------------------------------------------
  ! CONSTRUCT THE HAMILTONIAN
  !----------------------------------------------------------------
  !
  allocate ( t(ngm, nk0), 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 ik0 = 1, nk0
    do ig = 1, ngm
      kplusg = xk0(:,ik) + g(:,ig)
      t ( ig, ik0) = tpiba2 * dot_product ( kplusg, kplusg )
    enddo
  enddo
  !
  allocate ( h (ngm, ngm, nk0) )
  do ig = 1, ngm
    do ik0 = 1, nk0
        h ( ig, ig, ik0) = t ( ig, ik0)
    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 ik0 = 1, nk0
       if (ideltag.ne.ngm) &
          h ( ig1, ig2, ik0) = ss (ideltag) * vs (ideltag) 
          h ( ig2, ig1, ik0) = h ( ig1, ig2, ik0) 
     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)
  !
  allocate ( e (nbnd), psi (ngm, nbnd), eval(ngm0) )
  allocate ( hk(ngm0,ngm0), u(ngm0,ngm0), fv1(ngm0), fv2(ngm0) )
  allocate ( g2kin (ngm) )
  allocate ( scrcoul (ngm_sigma, ngm_sigma, fbins) )
  allocate ( greenf (ngm_sigma, ngm_sigma, nw) )
!@
nw_sigma = 10
!@
  allocate ( sigma (ngm_sigma, ngm_sigma, nw_sigma) )
  sigma = czero
!@@
! temporary
!@nq = nk0
allocate ( xq(3,nq), wq(nq) )
wq = one/float(nq)
xq = xk0
!@@
  !
  do ik0 = 1, nk0 
    !
    ! generate the uniform {q} grid for the Coulomb interaction
    ! and the grid {k0-q} for the Green's function
    !   
    ! the grids {k} and {k+q} for the dVscf will be obtained
    ! by shuffling the {q} grid
    !
  !....
  !....
    !
    ! now loop over {k0-q} grid for the Green's function
    ! and on {q} for the screened Coulomb interaction
    !
!@    do iq = 1, nq
   do iq = 1, 10
      !
      !  kmq = k0 - q
      !
      kmq = xk0(:,ik0) - xq(:,iq)
      !
!@      do ig1 = 1, ngm_sigma
!@        do ig2 = 1, ngm_sigma
      do ig1 = 1, 1
        do ig2 = 1, 3
         !
         ! the k-dependent kinetic energy in Ry
         ! [Eq. (14) of Ihm,Zunger,Cohen J Phys C 12, 4409 (1979)]
         !
         do ig = 1, ngm
           kplusg = kmq + g(:,ig)
           g2kin ( ig ) = tpiba2 * dot_product ( kplusg, kplusg )
         enddo
         !
write(6,*) iq,ig1,ig2
         call green ( ig1, ig2, g2kin, vr, greenf (ig1, ig2, :) )
         !
       enddo
      enddo
      !
    !... here we should generate all psi_occ on the uniform grid
    !... of q points and carry inside coulomb()
    !...
      !
      call coulomb ( xq(:,iq), scrcoul )
      !
      ! combine Green's function and screened Coulomb ( sum_q wq = 1 )
      !
      do iw = 1, nw_sigma
        !
        do ig1 = 1, ngm_sigma
          do ig2 = 1, ngm_sigma
            !
            do iwp = 1, nw_sigma
              !
           !  iwmwp = ...
              iwmwp = 4 !@ temporary 
              !
           !  wp = ...
              wp = 2.d0 !@ temporary
              !
!
! not sure about the sign of the deltas inside G and W...
!
              sigma(ig1,ig2,iw) = sigma(ig1,ig2,iw) &
                + wq (iq) * ci / twopi * exp ( -ci * delta * wp ) & 
                * greenf(ig1,ig2,iwmwp) * scrcoul (ig1,ig2,iwp)
              !
            enddo
            !
          enddo
        enddo
write(6,*) 'after sigma'
        !
      enddo
      !
    enddo ! end loop on {k0-q} and {q}
    !
  enddo ! end loop on {k0}
  !


  !
  end program gwhs
  !----------------------------------------------------------------
  !