!
  !----------------------------------------------------------------
  subroutine coulomb ( xq, 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
  use kspace
  implicit none
  !
  real(dbl) :: xq(3), ui, uj, uk, qg2
  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
  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
  !
  allocate ( xk (3,nks), wk(nks) )
  !

  recl = 2 * nbnd_occ * ngm  ! 2 stands for complex
  unf_recl = DIRECT_IO_FACTOR * recl
  open ( iunwfc, file = "./silicon.wfc", iostat = ios, form = 'unformatted', &
       status = 'unknown', access = 'direct', recl = unf_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)
  !

  !
  !  generate uniform {k} and {k+q} grids - no symmetry-reduction for now
  !  (MP method)
  !
  count = 0
  do i = 1, nk1
    ui = (k1 + 2.d0 * i - nk1 - 1.d0) / (2.d0 * nk1)
    do j = 1, nk2
      uj = (k2 + 2.d0 * j - nk2 - 1.d0) / (2.d0 * nk2)
      do k = 1, nk3
        uk = (k3 + 2.d0 * k - nk3 - 1.d0) / (2.d0 * nk3)
        count = count + 1
        xk (:, count) = ui * bg(:,1) + uj * bg(:,2) + uk * bg(:,3)
      enddo
    enddo
  enddo
  !
  ! include spin degeneracy
  wk = 2.d0 / float ( count )
  !
  !  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)+xq
    wk(ikk) = wk(ik)
    wk(ikq) = 0.d0
  enddo
  !
  do ik = 1, nks
    !
    !  find wavefunctions for this kpoint, occupied bands
    !
    call eigenstates ( xk(:, ik), vr, g2kin, psi, et(:,ik) ) 
    !
    !  direct write to file
    !
    write ( iunwfc, rec = ik, iostat = ios) psi 
    !
  enddo
  !
  !  loop on bare perturbations ig and fixed q point
  !
!@  do ig = 1, ngm_sigma
  do ig = 1, 3
    !
!@    do iw = 1, fbins
    do iw = 2, 2
      !
      w = ( fmin + ( fmax - fmin ) * float(iw-1) / float(fbins-1) ) / ryd2ev
      !
      dvbare = czero
      dvscf  = czero
      qg2 = (g(1,ig)+xq(1))**2.d0 + (g(2,ig)+xq(2))**2.d0 + (g(3,ig)+xq(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, xq, 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
  ! 
  deallocate ( xk, wk )
  close ( iunwfc )
  close ( iubar )
  close ( iudwf )
  close ( iudwfp )
  close ( iudwfm )
  !
  return
  end subroutine coulomb
  !----------------------------------------------------------------
  !