SUBROUTINE sigma_c(ik0) 
! G TIMES W PRODUCT
  USE io_global,     ONLY : stdout, ionode_id, ionode
  USE kinds,         ONLY : DP
  USE lsda_mod,      ONLY : nspin
  USE constants,     ONLY : e2, fpi, RYTOEV, tpi, eps8, pi
  USE disp,          ONLY : nqs, nq1, nq2, nq3, wq, x_q, xk_kpoints
  USE control_gw,    ONLY : lgamma, eta, godbyneeds, padecont, cohsex, modielec, &
                            do_diag_g, do_diag_w, real_freq_w
  USE klist,         ONLY : wk, xk
  USE io_files,      ONLY : iunigk, prefix, tmp_dir
  USE wvfct,         ONLY : nbnd, npw, npwx, igk, g2kin, et
  USE eqv,           ONLY : evq, eprec
  USE freq_gw,       ONLY : fpol, fiu, nfs, nfsmax, &
                            nwcoul, nwgreen, nwalloc, nwsigma, wtmp, wcoul, &
                            wgreen, wsigma, wsigmamin, wsigmamax, &
                            deltaw, wcoulmax, ind_w0mw, ind_w0pw
  USE units_gw,      ONLY : iuncoul, iungreen, iunsigma, lrsigma, lrcoul, lrgrn, iuwfc, lrwfc
  USE qpoint,        ONLY : xq, npwq, igkq, nksq, ikks, ikqs
  USE gwsigma,       ONLY : ngmsco, sigma, sigma_g, nrsco, nlsco, fft6_g2r, ecutsco, ngmsig,&
                            nr1sco, nr2sco, nr3sco, ngmgrn, ngmpol
  USE gvect,         ONLY : g, ngm, ecutwfc, nl
  USE cell_base,     ONLY : tpiba2, tpiba, omega, alat, at
  USE symm_base,     ONLY : nsym, s, time_reversal, t_rev, ftau, invs, nrot
  USE mp_global,     ONLY : inter_pool_comm, intra_pool_comm, mp_global_end, mpime, npool, &
                            nproc_pool, me_pool, my_pool_id, nproc
  USE mp,            ONLY : mp_barrier, mp_bcast, mp_sum
  USE fft_scalar,    ONLY : cfft3d
  USE modes,         ONLY : nsymq, invsymq, gi, gimq, irgq, irotmq, minus_q
  USE wavefunctions_module, ONLY : evc
  USE control_flags,        ONLY : noinv
  IMPLICIT NONE
  COMPLEX(DP)         :: ci, czero
  COMPLEX(DP)         :: phase
  COMPLEX(DP)         :: aux (nrsco)
!For running PWSCF need some variables 
  LOGICAL             :: pade_catch
  LOGICAL             :: found_q
  LOGICAL             :: limit, limq, inv_q, found
!Pade arrays
  COMPLEX(DP), ALLOCATABLE :: z(:), u(:), a(:)
!W arrays 
  COMPLEX(DP), ALLOCATABLE :: scrcoul_g (:,:,:)
  COMPLEX(DP), ALLOCATABLE :: scrcoul_g_R (:,:,:)
  COMPLEX(DP), ALLOCATABLE :: scrcoul_pade_g (:,:)
  COMPLEX(DP), ALLOCATABLE :: scrcoul(:,:)
!v array
!COMPLEX(DP), ALLOCATABLE ::  barcoul(:,:), barcoulr(:,:), barcoul_R(:,:)
  REAL(DP) :: qg2, qg, qxy, qz
!G arrays:
  COMPLEX(DP), ALLOCATABLE :: greenf_g(:,:), greenfr(:,:)
!Integration Variable 
  COMPLEX(DP) :: cprefac
!FREQUENCY GRIDS/COUNTERS
  INTEGER  :: iwim, iw, ikq
  INTEGER  :: iw0, iw0mw, iw0pw
  REAL(DP) :: w_ryd(nwcoul)
!COUNTERS
  INTEGER :: ig, igp, irr, icounter, ir, irp
  INTEGER :: iqstart, iqstop, iqs, nkr
  INTEGER :: iq, ipol, iqrec
  INTEGER :: ikmq, ik0, ik, nkpool
  INTEGER :: rec0, ios
  INTEGER :: counter, ierr
  INTEGER :: inversym, screening
!SYMMETRY
  REAL(DP)              :: xq_ibk(3), xq_ibz(3)
  INTEGER               :: isym, jsym
  INTEGER, ALLOCATABLE  :: gmapsym(:,:)
  COMPLEX(DP), ALLOCATABLE  ::  eigv(:,:)
!For G^NA
  INTEGER     :: igkq_ig(npwx) 
  INTEGER     :: igkq_tmp(npwx) 
  INTEGER     :: ss(3,3)
!q-vector of coulomb potential xq_coul := k_{0} - xk(ik)
  REAL(DP) :: xq_coul(3)
  REAL(DP) :: rcut, spal, zcut
  INTEGER  :: ibnd
!CHECK FOR NAN's
  REAL(DP)     :: ar, ai
!For dirac delta fxn.
  REAL(DP)     :: dirac, x, support
!File related:
  character(len=256) :: tempfile, filename
!Complete file name
  integer*8 :: unf_recl

#define DIRECT_IO_FACTOR 8 
! #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
! iG(W-v)
   ALLOCATE ( scrcoul_g      (ngmpol, ngmpol, nfs)     )
   ALLOCATE ( scrcoul_g_R    (ngmpol, ngmpol, nfs)     )
   ALLOCATE ( scrcoul_pade_g (ngmpol, ngmpol)          )
   ALLOCATE ( greenf_g       (ngmgrn, ngmgrn)          )
!These go on the big grid...
   ALLOCATE ( scrcoul        (nrsco, nrsco)            )
   ALLOCATE ( greenfr        (nrsco, nrsco)            )
!barcoul
!  ALLOCATE ( barcoul     (ngmpol, ngmpol)  )
!  ALLOCATE ( barcoulr    (nrsco,  nrsco)   )
!  ALLOCATE ( barcoul_R   (ngmpol,  ngmpol) )
!Technically only need gmapsym up to ngmpol or ngmgrn...
   ALLOCATE ( gmapsym  (ngm, nrot)   )
   ALLOCATE ( eigv     (ngm, nrot)   )
!This is a memory hog...
   ALLOCATE ( sigma          (nrsco, nrsco, nwsigma)   )
   ALLOCATE  (z(nfs), a(nfs), u(nfs))
!do we need two instances of the frequency?
   w_ryd(:) = wcoul(:)/RYTOEV
   WRITE(6," ")
   WRITE(6,'(4x,"Direct product GW for k0(",i3," ) = (",3f12.7," )")') ik0, (xk(ipol, ik0), ipol=1,3)
   WRITE(6," ")
   WRITE(6,'(4x, "ngmsco, ", i4, " nwsigma, ", i4)') ngmsco, nwsigma
   WRITE(6,'(4x, "nrsco, ", i4, " nfs, ", i4)') nrsco, nfs
   ci = (0.0d0, 1.d0)
   czero = (0.0d0, 0.0d0)
   sigma(:,:,:) = (0.0d0, 0.0d0)
   limit=.false.

   CALL start_clock('sigmac')
   CALL gmap_sym(nrot, s, ftau, gmapsym, eigv, invs)

   if(allocated(sigma)) then
     WRITE(6,'(4x,"Sigma allocated")')
   else
     WRITE(6,'(4x,"Sigma too large!")')
     CALL mp_global_end()
     STOP
   endif

   WRITE(6,'("nsym, nsymq, nsymbrav ", 3i4)'), nsym, nsymq, nrot 

!Set appropriate weights for points in the Brillouin zone.
!Weights of all the k-points are in odd positions in list.
!nksq is number of k points not including k+q.
   wq(:) = 0.0d0
   do iq = 1, nksq
      if(lgamma) then
         write(6, '(" lgamma ")')
         wq(iq) = 0.5d0*wk(iq) 
      else
         wq(iq) = 0.5d0*wk(2*iq-1) 
      endif
   enddo
!Every processor needs access to the files: _gw0si.coul1 and _gw0si.green1
call mp_barrier(inter_pool_comm)
#ifdef __PARA
if(.not.ionode) then
!OPEN coulomb file (only written to by head node).
   filename = trim(prefix)//"."//"coul1"
   tempfile = trim(tmp_dir) // trim(filename)
   unf_recl = DIRECT_IO_FACTOR * int(lrcoul, kind=kind(unf_recl))
   open(iuncoul, file = trim(adjustl(tempfile)), iostat = ios, &
   form = 'unformatted', status = 'OLD', access = 'direct', recl = unf_recl)
!OPEN green file (only written to by head node as well).
   filename = trim(prefix)//"."//"green1"
   tempfile = trim(tmp_dir) // trim(filename)
   unf_recl = DIRECT_IO_FACTOR * int(lrgrn, kind=kind(unf_recl))
   open(iungreen, file = trim(adjustl(tempfile)), iostat = ios, &
   form = 'unformatted', status = 'OLD', access = 'direct', recl = unf_recl)
endif
#endif

#ifdef __PARA
      npool = nproc / nproc_pool
      write(stdout,'("npool", i4, i5)') npool, nksq
      if (npool.gt.1) then
      ! number of q-vec per pool and remainder
        nkpool = nksq / npool
        nkr = nksq - nkpool * npool
      ! the remainder goes to the first nkr pools
        if ( my_pool_id < nkr ) nkpool = nkpool + 1
        iqs = nkpool * my_pool_id + 1
        if ( my_pool_id >= nkr ) iqs = iqs + nkr
      ! the index of the first and the last g vec in this pool
        iqstart = iqs
        iqstop  = iqs - 1 + nkpool
        write (stdout,'(/4x,"Max n. of Kpoint in Sigma_C per pool = ",i5)') iqstop-iqstart+1
      else
#endif
       iqstart = 1
       iqstop = nksq
#ifdef __PARA
      endif
#endif
!ONLY PROCESSORS WITH K points to process: 
IF(iqstop-iqstart+1.ne.0) THEN
   WRITE(1000+mpime, '("mpime ", i4, "  iqstart, iqstop: ", 2i5)')mpime, iqstart, iqstop
   if (nksq.gt.1) rewind (unit = iunigk)
DO iq = iqstart, iqstop
      if (lgamma) then
          ikq = iq
      else
!k+q is in even positions of list (k,k+q)
          ikq = 2*iq
      endif
!q point for convolution \sum_{q \in IBZ_{k}} G_{k+q} W_{-q}
!  q = (k0 + q) - k0
!   xq_ibk(:) = xk(:,ikq) - xk_kpoints(:, ik0)
!  by actually searching for -q and calculating w of -q  
!  we ensure there is no confusion in which symmetry routine is required.
!  SYMMFIX
!  -q =  k0 - (k0 + q)
!    xq_ibk(:) = xk_kpoints(:, ik0) - xk(:, ikq)
     xq_ibk(:) = xk(:, ikq)

   if(ionode)  print*, xk_kpoints(:,ik0)
   if(ionode)  print*, xq_ibk(:)

!Find which symmetry operation rotates xq_ibk back to
!The irreducible brillouin zone and which q \in IBZ it corresponds to.
!q is stored in the list x_q as positive q but all the calculations have
!been done at -q therefore we are just going to calculate \Sum G_{-k+q}W_{-q}
   inv_q=.false.
   call find_q_ibz(xq_ibk, s, iqrec, isym, found_q, inv_q)
!PRL
   iqrec = iq

   if(lgamma) npwq=npw 
   write(6, *)  
   write(6, '("xq_IBK point")')
   write(6, '(3f11.7)') xq_ibk
   write(6, '("equivalent xq_IBZ point, symop, iqrec")')
   write(6, '(3f11.7, 2i4)') x_q(:,iqrec), isym, iqrec
   write(6,*)
!
   write(1000+mpime, *)  
   write(1000+mpime, '("xq_IBK point")')
   write(1000+mpime, '(3f11.7)') xq_ibk
   write(1000+mpime, '("equivalent xq_IBZ point, symop, iqrec")')
   write(1000+mpime, '(3f11.7, 2i4)') x_q(:, iqrec), isym, iqrec
!Need a loop to find all plane waves below ecutsco when 
!igkq takes us outside of this sphere.
!igkq_tmp is gamma centered index up to ngmsco,
!igkq_ig  is the linear index for looping up to npwq.
!Dielectric Function should be written to file at this point
!So we read that in, rotate it, and then apply the coulomb operator.
     scrcoul_g(:,:,:)   = dcmplx(0.0d0, 0.0d0)
     scrcoul_g_R(:,:,:) = dcmplx(0.0d0, 0.0d0)
     if(.not.modielec) CALL davcio(scrcoul_g, lrcoul, iuncoul, iqrec, -1)
!Rotate G_vectors for FFT.
!In EPW FG checked that gmapsym(gmapsym(ig,isym),invs(isym)) = ig
!I have checked that here as well and it works.
!write(6,'(11i4)')gmapsym(gmapsym(:ngmsco,isym),invs(isym))
!Take a look at rotated G-vect map:
!for systems with ftau != 0 should also multiply each matrix element
!by phase, e^{iG'\tau}.
!Another one of these nested loops. 
!Two strategies or alternatives:
!1) Could pad scrcoul_g_R so that all the vectors still fall in 
!   side of it up to ngmsco then trim them off later.
!2) Modify gmapsym so that it only keeps vectors up to ngmsco.
       do igp = 1, ngmpol
          do ig = 1, ngmpol
             if((gmapsym(ig,isym).le.ngmpol).and.(gmapsym(igp,isym).le.ngmpol) &
                  .and.(gmapsym(ig,isym).gt.0).and.(gmapsym(igp,isym).gt.0)) then
                 do iwim = 1, nfs
!Rotating dielectric matrix with phase factor:
                   phase = eigv(ig,isym)*conjg(eigv(igp,isym))
                   scrcoul_g_R(ig, igp, iwim) = scrcoul_g(gmapsym(ig,isym), gmapsym(igp,isym),iwim)*phase
                 enddo
             endif
          enddo
       enddo


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Generate bare coulomb:                      !!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!HL using  sax cutoff
!this should be L_{z}/2
!rcut = 0.50d0*minval(sqrt(sum(at**2,1)))*alat*tpi
!rcut = rcut-rcut/50.0d0
rcut = (float(3)/float(4)/pi*omega*float(nq1*nq2*nq3))**(float(1)/float(3))
!Sax for silicon
!    rcut = 21.329d0
if(.not.modielec) then
DO iw = 1, nfs
   DO ig = 1, ngmpol
!SPHERICAL SCREENING
       qg2 = (g(1,ig) + xq_ibk(1))**2 + (g(2,ig) + xq_ibk(2))**2 + (g(3,ig)+xq_ibk(3))**2
       limq = (qg2.lt.eps8) 
       IF(.not.limq) then
           do igp = 1, ngmpol
              scrcoul_g_R(ig, igp, iw) = scrcoul_g_R(ig,igp,iw)*dcmplx(e2*fpi/(tpiba2*qg2), 0.0d0)
           enddo
       ENDIF
       qg = sqrt(qg2)
       spal = 1.0d0 - cos(rcut*sqrt(tpiba2)*qg)
!@ spherical truncation and enhanced correlation!.
!       spal = 1.0d0
!@ spherical truncation and enhanced correlation!.
!       spal = 4.0d0
!Normal case using truncated coulomb potential.
       if(.not.limq) then
          do igp = 1, ngmpol
              scrcoul_g_R(ig, igp, iw) = scrcoul_g_R(ig,igp,iw)*dcmplx(spal, 0.0d0)
          enddo
       else
!should only occur case iq->0, ig = 0 use vcut (q(0) = (4pi*e2*Rcut^{2})/2
!             write(6,*) (rcut)
!             write(6,*) (fpi*e2*(rcut**2))/2.0d0
!             write(6,*) ig, iw
!             write(6,*) g(:, ig)
!             write(6,*) xq_ibk(:)
!             write(6,*) qg2
!for omega=0,q-->0, G=0 the real part of the head of the dielectric matrix should be real
!we enforce that here:
             if(iw.eq.1) then
                scrcoul_g_R(ig, igp, iw) = real(scrcoul_g_R(ig,igp,iw))
             endif
             do igp = 1, ngmpol
                scrcoul_g_R(ig, igp, iw) = scrcoul_g_R(ig,igp,iw)*dcmplx((fpi*e2*(rcut**2))/2.0d0, 0.0d0)
             enddo
       endif
   ENDDO!ig
ENDDO!nfs
endif
!zeroing wings of W again!
    if(iq.eq.1) then
      do iw = 1, nfs
       do ig = 2, ngmpol
          scrcoul_g_R(ig,1,iw)  = dcmplx(0.0d0, 0.d0)
       enddo
       do igp = 2, ngmpol
          scrcoul_g_R(1,igp,iw) = dcmplx(0.0d0, 0.d0)
       enddo
      enddo
    endif
    if(.not.modielec) then 
        if(godbyneeds) then
        do ig = 1, ngmpol
          do igp = 1, ngmpol 
!         for godby-needs plasmon pole the algebra is done assuming real frequency*i...
!         that is: the calculation is done at i*wp but we pass a real number as the freq. 
               do iw = 1, nfs
                 z(iw) = dcmplx(aimag(fiu(iw)), 0.0d0)
                 u(iw) = scrcoul_g_R(ig, igp, iw)
               enddo
               CALL godby_needs_coeffs(nfs, z, u, a)
               do iw = 1, nfs 
!just overwrite scrcoul_g_R with godby-needs coefficients.
                  scrcoul_g_R (ig, igp, iw) = a(iw)
               enddo
          enddo
        enddo
        else if (padecont) then
          do igp = 1, ngmpol
           do ig = 1, ngmpol
!     Pade input points on the imaginary axis
              do iw = 1, nfs
                 z(iw) = fiu(iw)
                 u(iw) = scrcoul_g_R (ig, igp, iw)
              enddo
! Condition to catch the zero'd qings of the dielectric matrix at q-> 0.
              if(.not.do_diag_w) then
                if(iqrec.ne.1) then
                       call pade_coeff (nfs, z, u, a)
                else if (iqrec.eq.1) then
                     if ((ig.eq.1).xor.(igp.eq.1)) then
                         a=(0.0d0, 0.0d0)
                     else
                         call pade_coeff (nfs, z, u, a)
                     endif
                endif
              else if(do_diag_w) then
                   if(ig.ne.igp) then
                       a=(0.0d0, 0.0d0)
                   else
                       call pade_coeff (nfs, z, u, a)
                   endif
              endif
!     Overwrite scrcoul with Pade coefficients to be passed to pade_eval.
             do iw = 1, nfs 
                scrcoul_g_R (ig, igp, iw) = a(iw)
             enddo
           enddo !enddo on ig
         enddo  !enddo on igp
        else if (real_freq_w) then
             do iw = 1, nfs
                z(iw) = fiu(iw)
                u(iw) = scrcoul_g_R (ig, igp, iw)
             enddo
        else if(.not.padecont.and..not.godbyneeds) then
                 WRITE(6,'("No screening model chosen!")')
        endif
    endif
!Start integration over iw +/- wcoul. 
    WRITE(6,'("Starting Frequency Integration")')
    DO iw = 1, nwcoul
       scrcoul_pade_g(:,:) = (0.0d0, 0.0d0)
    if(.not.modielec) then
        do ig = 1, ngmpol
           do igp = 1, ngmpol
              do iwim = 1, nfs
                  z(iwim) = fiu(iwim)
                  a(iwim) = scrcoul_g_R (ig,igp,iwim)
              enddo
              pade_catch=.false.
              if(padecont) then
                if(.not.do_diag_w) then
                  if(iqrec.ne.1) then
                     call pade_eval ( nfs, z, a, dcmplx(w_ryd(iw), eta), scrcoul_pade_g (ig,igp))
                  else if (iqrec.eq.1) then
                     if ((ig.eq.1).xor.(igp.eq.1)) then
                          scrcoul_pade_g (ig,igp) = (0.0d0, 0.0d0)
                     else
                          call pade_eval ( nfs, z, a, dcmplx(w_ryd(iw), eta), scrcoul_pade_g (ig,igp))
                     endif
                  endif
                else if(do_diag_w) then
                   if(ig.ne.igp) then
                       scrcoul_pade_g (ig,igp) = (0.0d0, 0.0d0)
                   else
                       call pade_eval ( nfs, z, a, dcmplx(w_ryd(iw), eta), scrcoul_pade_g (ig,igp))
                   endif
                endif
              else if(godbyneeds) then
                 scrcoul_pade_g(ig,igp) = a(2)/(dcmplx(w_ryd(iw)**2,0.0d0)-(a(1)-(0.0d0,1.0d0)*eta)**2)
              else if(real_freq_w) then
!all real frequencies are written to file.
                 scrcoul_pade_g(ig,igp) = u(iw)
              else 
                 WRITE(6,'("No screening model chosen!")')
                 STOP
                 call mp_global_end()
              endif
           enddo
        enddo
    else if(modielec)  then
          do ig = 1, ngmpol
             call mod_diel(ig, xq_ibk, w_ryd(iw), scrcoul_pade_g(ig,ig), 1)
            !rewrite as fxns:
            !scrcoul_pade_g(ig,ig) = mod_dielec_(xq,w_ryd,ig)*v_(q+G,truncation)
             qg2 = (g(1,ig) + xq_ibk(1))**2 + (g(2,ig) + xq_ibk(2))**2 + (g(3,ig)+xq_ibk(3))**2
             limq = (qg2.lt.eps8) 
             IF(.not.limq) then
                scrcoul_pade_g(ig, ig) = scrcoul_pade_g(ig,ig)*dcmplx(e2*fpi/(tpiba2*qg2), 0.0d0)
             ENDIF
             qg = sqrt(qg2)
!spal = 1.0d0 - cos(rcut*sqrt(tpiba2)*qg)
!turning off spherical truncation
             spal = 1.0d0
!Normal case using truncated coulomb potential.
             if(.not.limq) then
                   scrcoul_pade_g(ig, ig) = scrcoul_pade_g(ig,ig)*dcmplx(spal, 0.0d0)
             else
                   scrcoul_pade_g(ig, ig) = scrcoul_pade_g(ig,ig)*dcmplx((fpi*e2*(rcut**2))/2.0d0, 0.0d0)
             endif
          enddo
    endif
!W(G,G';w)
!check
        czero = (0.0d0, 0.0d0)
        scrcoul(:,:) = czero
        do ig = 1, ngmpol
           aux(:) = czero
           do igp = 1, ngmpol
                    aux(nlsco(igp)) = scrcoul_pade_g(ig,igp)
           enddo
           call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, +1)
           do irp = 1, nrsco
                    scrcoul(ig, irp) = aux(irp) / omega 
           enddo
        enddo

        do irp = 1, nrsco
           aux = czero
           do ig = 1, ngmpol
              aux(nlsco(ig)) = conjg(scrcoul(ig,irp))
           enddo
           call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, +1)
            scrcoul(1:nrsco,irp) = conjg ( aux )
        enddo
!Now have W(r,r';omega')
       ! simpson quadrature: int_w1^wN f(w)dw = deltaw * 
       ! [ 1/3 ( f1 + fN ) + 4/3 sum_even f_even + 2/3 sum_odd f_odd ]
        cprefac = (deltaw/RYTOEV) * wq(iq) * (0.0d0, 1.0d0)/ tpi
         if ( iw/2*2.eq.iw ) then
            cprefac = cprefac * 4.d0/3.d0
         else
            cprefac = cprefac * 2.d0/3.d0
         endif

! Presumably the fxn should be zero at the cutoff but might as well be consistent...
         if (iw.eq.1) then
             cprefac = (1.0d0/3.0d0)*(deltaw/RYTOEV) * wq(iq) * (0.0d0, 1.0d0)/ tpi
         else if (iw.eq.nwcoul) then
             cprefac = (1.0d0/3.0d0)*(deltaw/RYTOEV) * wq(iq) * (0.0d0, 1.0d0)/ tpi
         endif

         do iw0 = 1, nwsigma
            iw0mw = ind_w0mw (iw0,iw)
            iw0pw = ind_w0pw (iw0,iw)
!rec0 = (iw0mw-1) * 1 * nqs + (ik0-1) * nqs + (iq-1) + 1
            rec0 = (iw0mw-1) * 1 * nksq + (iq-1) + 1

        if (.not.do_diag_g) then
            greenf_g(:,:) = czero
            CALL davcio( greenf_g, lrgrn, iungreen, rec0, -1 )
        else if(do_diag_g) then
            greenf_g(:,:) = czero
            CALL davcio( greenf_g, lrgrn, iungreen, rec0, -1 )
            greenfr(:,:) = czero
            do ig = 1, ngmpol
                greenfr(ig,1) = greenf_g(ig,ig)
            enddo
            greenf_g(:,:) = czero
            do ig = 1, ngmpol 
                greenf_g(ig,ig) = greenfr(ig,1)
            enddo
        endif

            greenfr(:,:) = czero
            do ig = 1, ngmgrn
               aux(:) = czero
               do igp = 1, ngmgrn
                  aux(nlsco(igp)) = greenf_g(ig,igp)
               enddo
               call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, +1)
               do irp = 1, nrsco
                 greenfr(ig, irp) = aux(irp)/omega
               enddo
            enddo
!additionally this means what I've actually calculated in green linsys is G = \psi_k^*(G)\psi_k(G')
!instead of \psi_{k}(G)\psi_{k}^*(G'), in Si inversion symmetry means:  
!G = \psi_{-k}(-G)\psi^*_-k(-G') 
!WHAT IF WE DON'T HAVE INVERSION SYMMETRY?
            do irp = 1, nrsco
               aux = czero
                    do ig = 1, ngmgrn
                          aux(nlsco(ig)) = conjg(greenfr(ig,irp))
                    enddo
                    call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, +1)
                    greenfr(1:nrsco,irp) = conjg ( aux )
            enddo
!!!!!!!!FFT
            sigma (:,:,iw0) = sigma (:,:,iw0) + cprefac * greenfr(:,:)*scrcoul(:,:)
!Now have G(r,r';omega-omega')
!rec0 = (iw0pw-1) * 1 * nqs + (ik0-1) * nqs + (iq-1) + 1
            rec0 = (iw0pw-1) * 1 * nksq + (iq-1) + 1
            if (.not.do_diag_g) then
                greenf_g(:,:) = czero
                CALL davcio( greenf_g, lrgrn, iungreen, rec0, -1 )
            else if(do_diag_g) then 
                greenf_g(:,:) = czero
                CALL davcio( greenf_g, lrgrn, iungreen, rec0, -1 )
                greenfr(:,:) = czero
                do ig = 1, ngmpol 
                   greenfr(ig,1) = greenf_g(ig,ig)
                enddo
                greenf_g(:,:) = czero
                do ig = 1, ngmpol 
                   greenf_g(ig,ig) = greenfr(ig,1)
                enddo
            endif
!Rotate green...
!Inlining FFT:
            greenfr(:,:) = czero
            do ig = 1, ngmgrn
               aux(:) = czero
               do igp = 1, ngmgrn
                  aux(nlsco(igp)) = greenf_g(ig,igp)
               enddo
               call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, +1)
               do irp = 1, nrsco
                  greenfr(ig, irp) = aux(irp) / omega
               enddo
            enddo
            do irp = 1, nrsco
               aux = czero
               do ig = 1, ngmgrn
                       aux(nlsco(ig)) = conjg(greenfr(ig,irp))
               enddo
               call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, +1)
               greenfr(1:nrsco,irp) = conjg(aux) 
            enddo
!Should be Sigma^{*}=G^{*}W^{*}
            sigma (:,:,iw0) = sigma (:,:,iw0) + cprefac * greenfr(:,:)*scrcoul(:,:)
        ENDDO !on iw0  
      ENDDO ! on frequency convolution over w'
    ENDDO ! end loop iqstart, iqstop 
ENDIF

    DEALLOCATE ( gmapsym          )
    DEALLOCATE ( greenfr          )
    DEALLOCATE ( greenf_g         )
    DEALLOCATE ( scrcoul          )
    DEALLOCATE ( scrcoul_pade_g   )
    DEALLOCATE ( scrcoul_g, scrcoul_g_R )
    DEALLOCATE ( z,a,u )

#ifdef __PARA
    CALL mp_barrier(inter_pool_comm)
    CALL mp_sum(sigma, inter_pool_comm)
    CALL mp_barrier(inter_pool_comm)
#endif __PARA


IF (ionode) then
  ALLOCATE ( sigma_g (ngmsco, ngmsco, nwsigma))
  if(allocated(sigma_g)) then
     WRITE(6,'(4x,"Sigma_g allocated")')
  else
     WRITE(6,'(4x,"Sigma_g too large!")')
     CALL mp_global_end()
     STOP
  endif

    WRITE(6,'(4x,"Sigma in G-Space")')
    sigma_g = (0.0d0,0.0d0)
    do iw = 1, nwsigma
      do ir = 1, nrsco
        aux = (0.0d0, 0.0d0)
        do irp = 1, nrsco
           aux(irp) = sigma(ir,irp,iw)
        enddo
        call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, -1)
        do igp = 1, ngmsco
           sigma (ir, igp, iw) = aux(nlsco(igp))
        enddo
      enddo
      do igp = 1, ngmsco
        aux = czero
        do ir = 1, nrsco
           aux(ir) = conjg ( sigma(ir,igp,iw) )
        enddo
        call cfft3d (aux, nr1sco, nr2sco, nr3sco, nr1sco, nr2sco, nr3sco, -1)
        do ig = 1, ngmsco
           sigma (ig,igp,iw) = conjg ( aux( nlsco( ig )) ) * omega
        enddo
      enddo
    enddo
    do iw = 1, nwsigma
      do igp = ngmsco + 1, nrsco
         do ig = ngmsco + 1, nrsco
               sigma (ig, igp, iw) = (0.0d0, 0.0d0)
         enddo 
      enddo
    enddo
!sigma_g = sigma(1:ngmsco,1:ngmsco,:)
    do iw = 1, nwsigma
     do igp = 1, ngmsco
      do ig = 1, ngmsco
           sigma_g(ig,igp,iw)  = sigma(ig,igp,iw)
      enddo
     enddo
    enddo 
!Now write Sigma in G space to file. 
!HL Original:
!Just storing in first record no matter what k-point.
!HLFULLSIG
    CALL davcio (sigma_g, lrsigma, iunsigma, ik0, 1)
!or could store sigma same way green's fxn is stored...
    WRITE(6,'(4x,"Sigma Written to File")')
    CALL stop_clock('sigmac') 
    DEALLOCATE ( sigma_g  )
ENDIF !ionode
    call mp_barrier(inter_pool_comm)
    DEALLOCATE ( sigma    )
    RETURN
END SUBROUTINE sigma_c