SUBROUTINE fft6_r2g_coulomb (f_r, f_g) 

! 6-D fourier transform used for W(r,r') to W(G, G')
! which agrees with convention in FG's paper.

  USE io_files,         ONLY : prefix, iunigk
  USE gwsigma,          ONLY : ngmsig, nrsig, nr1sig, nr2sig, nr3sig, nlsig
  USE kinds,            ONLY : DP 
  USE cell_base,        ONLY : omega, alat
  USE control_gw,       ONLY : lgamma
  USE wvfct,            ONLY : npw,npwx, igk
  USE qpoint,           ONLY : xq, npwq, igkq, nksq, ikks, ikqs
  USE gsmooth,          ONLY : nrxxs, nrx1s, nrx2s, nrx3s, nr1s, nr2s, nr3s, nls
  USE gvect,    ONLY : ngm, nrxx, g, nr1, nr2, nr3, nrx1, nrx2, nrx3, nl

IMPLICIT NONE

  INTEGER                          :: ios
  INTEGER                          :: ig, igp, ir, irp
  COMPLEX(DP)                      :: czero

  COMPLEX(DP)                      :: f_g(ngmsig,ngmsig)
  COMPLEX(DP)                      :: aux(nrsig)
  COMPLEX(DP)                      :: f_r (nrsig,nrsig)


!EXPLICIT TEST OF FFT
  do ir = 1, nrsig
     aux = (0.0d0, 0.0d0)
     do irp = 1, nrsig
         aux(irp) = f_r(ir, irp)
     enddo
 !HL forward backward works on full grid
 !   call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1)

     call cft3 (aux, nr1sig, nr2sig, nr3sig, nr1sig, nr2sig, nr3sig, -1)

     do igp = 1, ngmsig
        !f_r (ir,igp) = aux( nl(igp))
        f_r (ir,igp) = aux( nlsig(igp))
     enddo
  enddo

  do igp = 1, ngmsig
     aux = czero
     do ir = 1, nrsig
        aux(ir) = conjg ( f_r(ir,igp) )
     enddo

     !HL forward backward works on smooth grid.
     !call cft3 (aux, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1)

     call cft3 (aux, nr1sig, nr2sig, nr3sig, nr1sig, nr2sig, nr3sig, -1)

     do ig = 1, ngmsig
        !f_r (ig,igp) = conjg ( aux( nl( ig )) ) * omega
        f_r (ig,igp) = conjg ( aux( nlsig( ig )) ) * omega
     enddo
  enddo

  f_g(1:ngmsig,1:ngmsig) = f_r(1:ngmsig, 1:ngmsig)

END SUBROUTINE fft6_r2g_coulomb