!
  !----------------------------------------------------------------
  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, 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
  !
! 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_sigma
    !
    write(6,'(4x,"Screened Coulomb: q =",3f7.3,"  G'' =",3f7.3)') xq, g(:,ig)
    !
    do iw = 1, fbins
      !
      w = ( fmin + ( fmax - fmin ) * float(iw-1) / float(fbins-1) ) / ryd2ev
      !
      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 ( 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 ( ig ) = cone 
        call cfft3 ( dvbare, nr1, nr2, nr3,  1)
        !
!        ! DEBUG
!        !
!        ! here I only perform 1 iteration. In output we have
!        ! eps instead of eps^-1 (eps = delta - v_c * chi_0 )
!        ! The value of eps for q->0 agrees with dielec_mat.f90 within 0.05%
!        !
!        call solve_linter_notSCF ( dvbare, dvscfp, xq, et, vr, +w )
!        call solve_linter_notSCF ( dvbare, dvscfm, xq, et, vr, -w )
!        dvscf = 0.5d0 * ( dvscfp + dvscfm )
!        !
!        ! go to G-space, this is  -v_c * chi_0
!        call cfft3 ( dvscf, nr1, nr2, nr3, -1)
!        ! now add the delta function
!        dvscf ( nl (ig) ) = dvscf ( nl (ig) ) + cone
!        !
!        ! END DEBUG
        !
        ! 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 )
!       call solve_linter ( dvbare, dvscf, xq, et, vr, w )
        !
        !
        ! 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,'(50(2x,f9.5))') (dreal(dvscf ( nl(igp) )), igp = 1, 50)
        do igp = 1, 20
         write(70,'(f12.7,3x,f12.7))')  dreal(dvscf ( nl(igp) )), aimag(dvscf(nl(igp)))   
        enddo
!       write(6,'(f9.5,5x,2f9.5,5x,2f9.5)') xq(1), dvbare ( nl(1) ), dvscf ( nl(1) )
        !                                                            ^^^^^^^^^^^^^^^
        !                                                             eps^-1(0,0,q)
      endif
      !
      write(6, *)  w*ryd2ev, dvscf ( nl(1) )
      write(16,*)  w*ryd2ev, real ( dvscf ( nl(1) ) ), aimag ( dvscf ( nl(1) )  )
      !
    enddo
    !
  enddo
  ! 
  deallocate ( xk, wk )
  close ( iunwfc )
  close ( iubar )
  close ( iudwf )
  close ( iudwfp )
  close ( iudwfm )
  !
  return
  end subroutine coulomb
  !----------------------------------------------------------------
  !