SUBROUTINE gw_product(ik, iq) 

! This routine forms the GW product. 
! SOME GENERAL NOTES: There is some fairly brutal RYTOEV conversion going on here. Should try and approach
! this completely consistently. i.e. convert all frequency related things to rydberg at the start of the 
! calculation and then only convert to eV when printing out to screen to be read by user. 

  USE kinds, ONLY : DP
  USE gvect,    ONLY : ngm, nrxx, g, nr1, nr2, nr3, nrx1, nrx2, nrx3, nl
  USE gsmooth, ONLY : nrxxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, nls, ngms
  USE kinds, ONLY : DP
  USE ions_base, ONLY : nat
  USE lsda_mod, ONLY : nspin
  USE constants, ONLY : e2, fpi, RYTOEV, tpi
  USE disp,            ONLY : nqs
  USE klist,           ONLY: wk

!GW specific variables. Should have a user defined cutoff for the selfenergy, 
  USE freq_gw,    ONLY : fpol, fiu, nfs, nfsmax,&
                         nwcoul, nwgreen, nwalloc, nwsigma, wtmp, wcoul,&
                         wgreen, wsigma, w_ryd, wsigmamin, wsigmamax,&
                         deltaw, wcoulmax, ind_w0mw, ind_w0pw

  USE units_gw,    ONLY : iuncoul, iungreen, iunsigma
  USE gwsigma, ONLY : ngmsig, sigma

  USE qpoint,      ONLY : xq

!Pade arrays
  complex(DP), allocatable :: z(:), u(:), a(:)

!W arrays 
  complex(DP), allocatable :: scrcoul_g (:,:,:)
  complex(DP), allocatable :: scrcoul_pade_g (:,:)
  REAL(DP) :: qg2
  COMPLEX(DP) :: cprefac 
  !various frequency counters 
 
  INTEGER :: iwim, iw 
  INTEGER :: ig, igp, npe, irr, icounter
  INTEGER :: igstart, igstop, igpert
  INTEGER :: iq

! HL Need to think about how all of this is going to be parallelized. 
! #ifdef __PARA
!       scrcoul_g = czero
!       if (me.eq.1.and.mypool.eq.1) then
! #endif
! #ifdef __PARA
!       endif
!       use poolreduce to broadcast the results to every pool
!       call poolreduce ( 2 * ngms * ngms * nwim, scrcoul_g)
! #endif
! Might want to allocate scrcoul_g in gw_allocate

     allocate ( scrcoul_g (ngmsig, ngmsig, nfs) )
     allocate ( scrcoul_pade_g (ngms, ngms) )

     WRITE(6, '("Summing q = ", 3f7.3 )') xq(:)
     scrcoul_g = (0.0d0, 0.0d0)

     read ( iuncoul, rec = iq, iostat = ios) scrcoul_g

!CALCULATE SIGMA^{C}(r,r',w) = \int G(r,r',w + w') W (r, r', w') dw'
DO iq = 1, nqs 
     DO iw = 1, nwcoul
        do ig = 1, ngmsig
          do igp = 1, ngmsig
            do iwim = 1, nfs
               z(iwim) = dcmplx( 0.d0, fiu(iwim))
               a(iwim) = scrcoul_g (ig,igp,iwim)
            enddo
            
            call pade_eval ( nwim, z, a, dcmplx( fiu(iw), eta), scrcoul_pade_g (ig,igp))
            
          enddo
        enddo

        call fft6_g2r ( scrcoul_pade_g, scrcoul)

        !cprefac = deltaw/RYTOEV * wq (iq) * (0.0d0, 1.0d0)/ tpi
        !HL need to think about the weight of each of these q-points
        !and also how i'm meant to write the kpoints to file and how to refold
        !everything into the Brillouin Zone. 


        cprefac = deltaw/RYTOEV * wk (iq) * (0.0d0, 1.0d0)/ tpi

        do iw0 = 1, nwsigma
            !
            iw0mw = ind_w0mw (iw0,iw)
            iw0pw = ind_w0pw (iw0,iw)


!     generate green's function in real space for iw0mw (greenfm)
!
!     HL more parallel related shenanigans 
!       #ifdef __PARA
!             greenf_g = czero
!             if (me.eq.1.and.mypool.eq.1) then
!       #endif

            rec0 = (iw0pw-1) * nk0 * nq + (ik0-1) * nq + (iq-1) + 1
            read ( iungreen, rec = rec0, iostat = ios) greenf_g

!#ifdef __PARA
!          endif
!          use poolreduce to broadcast the results to every pool
!          call poolreduce ( 2 * ngms * ngms, greenf_g )
!#endif

          ! greenf_g is ngms*ngms, greenf is nrs*nrs

        call fft6_g2r ( greenf_g, greenfm )

! generate green's function in real space for iw0pw (greenfp)


!HL more parallel
!         #ifdef __PARA
!                 greenf_g = czero
!                 if (me.eq.1.and.mypool.eq.1) then
!         #endif

          rec0 = (iw0pw-1) * nk0 * nq + (ik0-1) * nq + (iq-1) + 1
          read ( iungreen, rec = rec0, iostat = ios) greenf_g

!#ifdef __PARA
!          endif
!          use poolreduce to broadcast the results to every pool
!          call poolreduce ( 2 * ngms * ngms, greenf_g )
!#endif
!          greenf_g is ngms*ngms, greenf is nrs*nrs
!
          call fft6_g2r ( greenf_g, greenfp )
          sigma (:,:,iw0,1) = sigma (:,:,iw0,1) + cprefac * ( greenfp + greenfm ) * scrcoul

       enddo !on iw0  
     ENDDO ! on frequency convolution over w'
ENDDO ! end looop on {k0-q} and {q} 

!EXCHANGE PART OF THE SELF-ENERGY 
!DO iq = 1, nqs
!ENDDO

END SUBROUTINE gw_product