!
  !----------------------------------------------------------------
  subroutine coulomb ( xq )
  !----------------------------------------------------------------
  ! 
  ! the screened Coulomb interaction W(r,r',w,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
  !
  integer :: ig, igp, ierr, ibnd, ios, recl, unf_recl
  real(dbl) :: g2kin(ngm), et(nbnd_occ, nks)
  complex(dbl) :: vr(nr), psi(ngm, nbnd_occ), dvbare(nr), dvscf(nr)
  !
  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)
  !

  !
  !  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
  !
! write(6,'(/4x,a)') repeat('-',67)
! write(6,'(4x,"Uniform k-point grid for the screened coulomb interaction"/)') 
! do ik = 1, nks
!    write ( 6, '(4x,"k(",i4,") = (",3f12.7,"), wk =",f12.7)') &
!        ik, (xk (ipol, ik) , ipol = 1, 3) , wk (ik)
! enddo
! write(6,'(4x,a/)') repeat('-',67)
  !
  do ik = 1, nks
    !
    !  find wavefunctions for this kpoint, occupied bands
    !
    call eigenstates ( xk(:, ik), vr, g2kin, psi, et(:,ik) ) 
    !
!   write(6,'(4x,i5,4(2x,f6.3))') ik, et(:,ik)*ryd2ev
!   write(6,'(4x,a/)') repeat('-',67)
    !
    !  direct write to file
    !
    write ( iunwfc, rec = ik, iostat = ios) psi 
    !
  enddo

  !
  !  loop on bare perturbations and fixed q point
  !
!@  do ig = 1, ngm
  do ig = 1, 1
    !
    write(6,'(4x,"Screened Coulomb: q =",3f7.3,"  G'' =",3f7.3)') xq, g(:,ig)
    !
    dvbare = 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 = 4pi/Omega / |q+G|^2 \delta(G,G') --> transform in real-space
      ! should we include the 1/Omega ?? @@
      !
      !  dvbare ( nl (ig) ) = dcmplx ( e2 * fpi / (tpiba2 * qg2) , zero )
      !
      ! every Fourier component acts independently:
      !   W(r,r') =   v(r,r') + chi(r,r'')*  W(r'',r') 
      !   W(r,G') =   v(r,G') + chi(r,r'')*  W(r'',G')
      ! c*W         c*v                    c*W                c=cnst
      !
      dvbare ( nl (ig) ) = cone  
      call cfft3 ( dvbare, nr1, nr2, nr3,  1)
      !
      ! solve self-consistently the linear system for this perturbation 
      !
      call solve_linter ( dvbare, dvscf, xq, et, vr )
      !
      ! DEBUG - transform dvbare and dvscf to G space and calculate eps^-1(G,G',q)
      !
      call cfft3 ( dvbare, nr1, nr2, nr3, -1)
      call cfft3 ( dvscf , nr1, nr2, nr3, -1)
      write(6,'(f9.5,5x,2f9.5,5x,2f9.5)') xq(1), dvbare ( nl(1) ), dvscf ( nl(1) )
      !                                                            ^^^^^^^^^^^^^^^
      !                                                             eps^-1(0,0,q)
      !
    endif
    !
  enddo
  
  deallocate ( xk, wk )
  close ( iunwfc )
  close ( iubar )
  close ( iudwf )
  !
  return
  end subroutine coulomb
  !----------------------------------------------------------------
  !