!
  !----------------------------------------------------------------
  subroutine coulomb ( xxq, scrcoul )
  !----------------------------------------------------------------
  ! 
  ! the screened Coulomb interaction W_q (G,G',w) for a given q
  ! 
  ! we include the spin degeneracy in the kpoint weights
  !
  !----------------------------------------------------------------
  !
  use parameters
  use constants
  use gspace, only : g, nr, nl, nr1, nr2, nr3, ngm, ngm_sigma
  use kspace
  implicit none
  !
  real(dbl) :: xxq(3), ui, uj, uk, qg2, xktmp(3)
  integer :: iq, count, i, j, k, ik, ipol, ikk, ikq
  !
  complex(dbl) :: scrcoul(ngm_sigma, ngm_sigma, fbins)
  integer :: ig, iw, igp, ierr, ibnd, ios, recl, unf_recl, ikf, fold(nq)
  real(dbl) :: g2kin(ngm), et(nbnd_occ, nks), w
  complex(dbl) :: vr(nr), psi(ngm, nbnd_occ), dvbare(nr), dvscf(nr)
  complex(dbl) :: dvscfp (nr), dvscfm (nr) 
  ! these are for +w and -w
  complex(kind=DP) :: evc (ngm, nbnd_occ), evq (ngm, nbnd_occ)
  !
  if (.not.allocated(xk)) allocate ( xk (3,nks), wk(nks) )
  !
  recl = 2 * nbnd_occ * ngm  ! 2 stands for complex
  unf_recl = DIRECT_IO_FACTOR * recl
  open ( iudwf, file = "./silicon.dwf", iostat = ios, form = 'unformatted', &
       status = 'unknown', access = 'direct', recl = unf_recl)
  open ( iubar, file = "./silicon.bar", iostat = ios, form = 'unformatted', &
       status = 'unknown', access = 'direct', recl = unf_recl)
  open ( iudwfp, file = "./silicon.dwfp", iostat = ios, form = 'unformatted', &
       status = 'unknown', access = 'direct', recl = unf_recl)
  open ( iudwfm, file = "./silicon.dwfm", iostat = ios, form = 'unformatted', &
       status = 'unknown', access = 'direct', recl = unf_recl)
  !
  write(6,'(4x,"Screened Coulomb interaction:")')

  !
  !  generate uniform {k} and {k+q} grids 
  !
  !  The {k} grid is taken to coincide with the {q} grid generated
  !  in gwhs.f90, the {k+q} grid is obtained by folding
  !  In this way we calculate the occupied states only once in gwhs.f90
  !
  do ik = 1, nksq
    xk (:, ik) = xq (:, ik) 
    ! include spin degeneracy
    wk = two / float ( nksq )
  enddo
  !
  !  find folding index and phase 
  !
!@@@
xxq = xk(:,4)
!@@@
  do ik = 1, nksq
    call ktokpmq (xk(:,ik), xxq, +1, fold(ik) )
  enddo
  !
  !  double the grid and add k+q in the even positions
  !
  do ik = nksq, 1, -1
    ikk = 2 * ik - 1
    ikq = 2 * ik
    xk(:,ikk) = xk(:,ik)
    xk(:,ikq) = xk(:,ik)+xxq
    wk(ikk) = wk(ik)
    wk(ikq) = 0.d0
  enddo
  !
  do ik = 1, nksq
    !
    ikk = 2 * ik - 1
    ikq = 2 * ik
    !
    !  the folded k+q
    !
    ikf = fold(ik)
    !
    !  read occupied wavefunctions for the folded k+q point
    !
    read ( iunwfc, rec = 2 * ikf - 1, iostat = ios) evq
    !
    !  set the phase factor of evq for the folding
    !
  !@@@
  !  evq = ... ??
  !@@@
    !
    !  and write it again in the right place
    !
    write ( iunwfc, rec = ikq, iostat = ios) evq 
    !
    ! the eigenvalues
    ! 
    et(:,ikk) = eval(:,ik)
    et(:,ikq) = eval(:,ikf) 
    !
!@@@@
!
! NOTA: questo qui giu' e' ok. Per essere sicuro che
!       il folding funziona bene devo assolutamente
!       fare il confronto delle funzioni d'onda QUI
!       cosi' poi sono tranquillo.
!
call eigenstates ( xk(:, ikq), vr, g2kin, psi, et(:,ikq) ) 
write ( iunwfc, rec = ikq, iostat = ios) psi 
!write(6,*) et(:,ikq)
!write(6,'(f12.7)') abs(psi(:,4))**2.d0-abs(evq(:,4))**2.d0
!write(6,'(f12.7)') abs(psi(:,3))**2.d0+abs(psi(:,4))**2.d0-abs(evq(:,3))**2.d0-abs(evq(:,4))**2.d0
!stop
!@@@@
    !
  enddo
  !
  !  loop on bare perturbations ig and fixed q point
  !
  do ig = 1, 3 !@ ngm_sigma
    !
    write(6,'(4x,"ig = ",i5)') ig
    !
    do iw = 1, 2 !@ fbins
      !
      write(6,'(4x,3x,"iw = ",i5)') iw
      !
      w = ( fmin + ( fmax - fmin ) * float(iw-1) / float(fbins-1) ) / ryd2ev
      !
      dvbare = czero
      dvscf  = czero
      qg2 = (g(1,ig)+xxq(1))**2.d0 + (g(2,ig)+xxq(2))**2.d0 + (g(3,ig)+xxq(3))**2.d0
      if (qg2 > 1.d-8) then
        !
        dvbare ( nl ( ig ) ) = cone 
        call cfft3 ( dvbare, nr1, nr2, nr3,  1)
        !
        ! solve self-consistently the linear system for this perturbation 
        ! _dyn is for the dynamical case (w/=0)
        !
        call solve_linter_dyn ( dvbare, dvscf, xxq, et, vr, w )
        !
        ! transform dvscf to G space 
        !
        call cfft3 ( dvscf , nr1, nr2, nr3, -1)
        ! 
      endif
      !
      do igp = 1, ngm_sigma
!
! not sure wheather it's this or the transpose      
!
! should we include the 1/Omega ?? 
!
        scrcoul (ig,igp,iw) = dvscf ( nl (igp) ) * dcmplx ( e2 * fpi / (tpiba2 * qg2) , zero )
      enddo
      !
    enddo
    !
  enddo
  ! 
  close ( iubar, status = 'delete' ) 
  close ( iudwf, status = 'delete' ) 
  close ( iudwfp, status = 'delete' ) 
  close ( iudwfm, status = 'delete' ) 
  !
  return
  end subroutine coulomb
  !----------------------------------------------------------------
  !