SUBROUTINE sigma_exch(ik0)
  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 ions_base,     ONLY : nat
  USE lsda_mod,      ONLY : nspin
  USE constants,     ONLY : e2, fpi, RYTOEV, tpi, eps8, pi
  USE disp,          ONLY : nqs, nq1, nq2, nq3, wq, x_q 
  USE control_gw,    ONLY : lgamma, eta
  USE klist,         ONLY : wk, xk
  USE io_files,      ONLY : prefix, iunigk
  USE wvfct,         ONLY : nbnd, npw, npwx, igk, g2kin, et
  USE cell_base,     ONLY : omega, tpiba2
  USE eqv,           ONLY : evq, eprec
  USE units_gw,      ONLY : iuncoul, iungreen, iunsigma, lrsigma, lrcoul, lrgrn, iuwfc, lrwfc
  USE qpoint,        ONLY : xq, npwq, igkq, nksq, ikks, ikqs
  USE gwsigma,       ONLY : ngmsig, nr1sex, nr2sex, nr3sex, nrsex, nlsex, ngmsex, fft6_g2r, &
                            sigma_ex, sigma_g_ex, ecutsex, ngmsco

IMPLICIT NONE

!ARRAYS to describe exchange operator.
  LOGICAL :: do_band, do_iq, setup_pw, exst
  COMPLEX(DP), ALLOCATABLE ::  greenf_na (:,:), greenf_nar(:,:)
  COMPLEX(DP), ALLOCATABLE ::  barcoul(:,:), barcoulr(:,:)
  REAL(DP) :: rcut, spal
  INTEGER :: ikmq, ik0, ik
  INTEGER :: ig, igp, npe, irr, icounter, ir, irp
  INTEGER :: iq, ipol, ibnd, jbnd, counter
  INTEGER :: rec0, ios
  REAL(DP) :: qg2, xq0s(3), qg, xxq(3)
  COMPLEX(DP) :: ZDOTC

!Arrays to handle case where nlsco does not contain all G vectors required for |k+G| < ecut
  INTEGER     :: igkq_ig(npwx) 
  INTEGER     :: igkq_tmp(npwx) 

! Self-Energy grid:
! iGv
  ALLOCATE ( barcoul     (ngmsex, ngmsex) )
  ALLOCATE ( barcoulr    (nrsex,  nrsex)  )

  ALLOCATE ( greenf_na   (ngmsex, ngmsex) )
  ALLOCATE ( greenf_nar  (nrsex,  nrsex)  )

  ALLOCATE ( sigma_g_ex  (ngmsex, ngmsex) ) 
  ALLOCATE ( sigma_ex    (nrsex, nrsex)   ) 

DO iq = 1, nqs
!NON-ANALYTIC PART OF THE GREEN'S FUNCTION
      CALL prepare_kmq(do_band, do_iq, setup_pw, iq, ik0)
      CALL run_pwscf(do_band)
      CALL initialize_gw()

      if(lgamma) then
         write(6, '("lgamma=.true.")')
         CALL davcio (evq, lrwfc, iuwfc, 1, -1)
      else
         CALL davcio (evq, lrwfc, iuwfc, 2, -1)
      endif

      if (nksq.gt.1) rewind (unit = iunigk)

      if (nksq.gt.1) then
           read (iunigk, err = 100, iostat = ios) npw, igk
100        CALL errore ('green_linsys', 'reading igk', abs (ios) )
      endif

      if (.not.lgamma.and.nksq.gt.1) then
           read (iunigk, err = 200, iostat = ios) npwq, igkq
200        CALL errore ('green_linsys', 'reading igkq', abs (ios) )
      endif

      if (lgamma) npwq = npw 

!Need a loop to find all plane waves below ecutsco when igkq takes us outside of this sphere.  
       counter  = 0
       igkq_tmp = 0
       igkq_ig  = 0

       do ig = 1, npwx
         if((igkq(ig).le.ngmsco).and.((igkq(ig)).gt.0)) then
             counter = counter + 1
             igkq_tmp (counter) = igkq(ig)
             igkq_ig  (counter) = ig
         endif
       enddo

       greenf_na = (0.0d0, 0.0d0)

!       do ig = 1, npwq
!        do igp = 1, npwq
       do ig = 1, counter
        do igp = 1, counter
         do ibnd = 1, nbnd
         !Working Exchange with if loops
         ! if(((igkq(igp).lt.ngmsig).and.(igkq(ig)).lt.ngmsig).and.&
         !    ((igkq(igp).gt.0).and.(igkq(ig)).gt.0)) then
         !     greenf_na(igkq(ig),igkq(igp)) = greenf_na(igkq(ig), igkq(igp)) + &
         !                                     tpi * (0.0d0, 1.0d0) * conjg((evq(ig,ibnd)))* &
         !                                     evq(igp,ibnd)
         ! endif

         greenf_na(igkq_tmp(ig),igkq_tmp(igp)) =  greenf_na(igkq_tmp(ig), igkq_tmp(igp)) + &
                                                  tpi * (0.0d0, 1.0d0) * conjg((evq(igkq_ig(ig),ibnd)))* &
                                                  evq(igkq_ig(igp), ibnd)
         enddo
        enddo
       enddo

       call fft6_g2r(ngmsex, nrsex, nlsex, greenf_na, greenf_nar, 2, 2)
      !call fft6_g2r(ngmsex, nrsex, nlsex, greenf_na, greenf_nar, 2, 1)

!Form bare coulomb in G-space and then FFT on to real space using nlsex grid. 

      rcut = (float(3)/float(4)/pi*omega*float(nq1*nq2*nq3))**(float(1)/float(3))
      xq0s = (/ 0.01 , 0.00, 0.00 /) ! this should be set from input
      barcoul(:,:) = (0.0d0,0.0d0)

!v(r,r')(q)
    do ig = 1, ngmsex
       qg = sqrt((g(1,ig)+x_q(1,iq))**2.d0 + (g(2,ig)+x_q(2,iq))**2.d0   + (g(3,ig )+x_q(3,iq))**2.d0)
       qg2 = (g(1,ig) + x_q(1,iq))**2.d0   + (g(2,ig) + x_q(2,iq))**2.d0 + ((g(3, ig))+ x_q(3,iq))**2.d0 

       if (qg < eps8) qg =  sqrt((g(1, ig) + xq0s(1))**2.d0 + (g(2, ig) + xq0s(2))**2.d0 & 
                                                            + (g(3, ig) + xq0s(3))**2.d0)

       if (qg2 < eps8) qg2 =  ((g(1, ig) + xq0s(1))**2.d0   + (g(2, ig) + xq0s(2))**2.d0 & 
                                                            + (g(3, ig) + xq0s(3))**2.d0)
       spal = 1.0d0 - cos (rcut * sqrt(tpiba2) * qg)
       barcoul (ig, ig) = e2 * fpi / (tpiba2*qg2) * dcmplx(spal, 0.0d0)
    enddo

!  #ifdef __PARA
!     use poolreduce to bring together the results from each pool
!     call poolreduce ( 2 * nrs * nrs, barcoul)
!  #endif
!  call fft6_g2r(barcoul, barcoulr)
!  write(6,'(26f7.3)') real(barcoul(:,:))

    barcoulr(:,:) = (0.0d0, 0.0d0)

    CALL fft6_g2r(ngmsex, nrsex, nlsex, barcoul, barcoulr, 2, 2)
!   CALL fft6_g2r(ngmsex, nrsex, nlsex, barcoul, barcoulr, 2, 2)
    sigma_ex = sigma_ex + wq(iq)* (0.0d0,1.0d0) / tpi * barcoulr * greenf_nar
    CALL clean_pw_gw(iq)
ENDDO ! on q

!sigma_g_ex operator in G space: won't write to disc unless necessary. 
   CALL sigma_r2g_ex(sigma_ex, sigma_g_ex) 

   DEALLOCATE (sigma_ex)
   DEALLOCATE (greenf_na, greenf_nar)
   DEALLOCATE (barcoul, barcoulr)

END SUBROUTINE sigma_exch