!
!
! Copyright (C) 2001-2008 Quantum_ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!-----------------------------------------------------------------------

SUBROUTINE coulomb

!-----------------------------------------------------------------------
! This subroutine is the main driver of the COULOMB self consistent cycle
! which gives as output the change of the wavefunctions and the
! change of the self-consistent potential due to the introduction of 
! a charge dvbare(nl(ig)) = 1.00 + i*0.00 at a single fourier component (G). The screened coulomb interaction 
! is given by W_{q}(G,G',iw) = (\delta_{GG'} + drhoscfs_{G,G',iw}) * (e2*fpi)/(q2g*tpiba2)
! W = eps^{-1} v 

  USE kinds, ONLY : DP
  USE ions_base, ONLY : nat
  USE gvect,    ONLY : ngm, nrxx, g, nr1, nr2, nr3, nrx1, nrx2, nrx3, nl
  USE gsmooth, ONLY : nrxxs, nr1s, nr2s, nr3s, nrx1s, nrx2s, nrx3s, nls
  USE constants, ONLY : e2, fpi
  USE cell_base, ONLY : alat, tpiba2
  USE lsda_mod, ONLY : nspin
  USE io_global,  ONLY : stdout, ionode
  USE uspp,  ONLY: okvan
  USE control_gw, ONLY : zue, convt, rec_code
  USE partial,    ONLY : done_irr, comp_irr
  USE modes,      ONLY : nirr, npert, npertx
  USE gwus,       ONLY : int3, int3_nc, int3_paw
  USE uspp_param, ONLY : nhm
  USE eqv,        ONLY : drhoscfs, dvbare
  USE paw_variables, ONLY : okpaw
  USE noncollin_module, ONLY : noncolin, nspin_mag
  USE recover_mod, ONLY : write_rec
  USE mp_global,  ONLY : inter_pool_comm, intra_pool_comm
  USE mp,         ONLY : mp_sum
  USE gwsigma,     ONLY : scrcoul
  USE qpoint,      ONLY : xq
  USE freq_gw,     ONLY : fpol, fiu, nfs, nfsmax

  IMPLICIT NONE
  REAL(DP) :: tcpu, get_clock
! timing variables

  REAL(DP) :: qg2
  INTEGER :: ig, igp, iw, nwim, npe, irr, icounter
  INTEGER :: igstart, igstop, igpert
 
 complex(DP), allocatable :: drhoaux (:,:) 

!
!counters on G and G' imaginary frequencies. 
!should define nwim globally.
!complex(DP) :: dvbare(nrxx).
!array for Bare Perturbation.
!
  LOGICAL :: exst
! used to test the recover file

  EXTERNAL get_clock
  CALL start_clock ('gwqscf')

!DUMMY VARIABLES
  nwim = 1
  npe = 1

!Change in charge density
ALLOCATE (drhoscfs(nrxx , nspin_mag))    
!ALLOCATE (drhoaux(nrxx , nspin_mag))
!drhoscfs(:,1) = (0.000, 0.0000)

!LOOP OVER G

WRITE(stdout,'(4x,"Number of frequencies", i4)') nfs
DO ig = 1, ngm
     igpert = ig
     do iw = 1, nfs

!HL- DEBUG
! Calculating non scf system. One iteration of the linear system solver should print out
! the dielectric constant.
!        CALL solve_linter_nonscf (igpert, drhoscfs)
!        WRITE( stdout, '(/,5x,"End of self-consistent calculation")')
!        IF (convt) THEN
!          WRITE( stdout, '(/,5x,"Convergence has been achieved ")')
!        ELSE
!          WRITE( stdout, '(/,5x,"No convergence has been achieved ")')
!        ENDIF
!HL- DEBUG
!        CALL solve_linter (igpert, drhoscfs)

        CALL solve_linter (igpert, iw, drhoscfs)

        IF (convt) WRITE( stdout, '(/,5x,"Convergence has been achieved ")')

         WRITE(6, "")
         WRITE(stdout,'(4x,"Frequency= ",3f7.3)') fiu(iw) 
         WRITE(stdout,'(4x,"Screened Coulomb: q =",3f7.3,"  G =",3f7.3)') xq(:), g(:,ig)
         call cft3 (drhoscfs, nr1, nr2, nr3, nrx1, nrx2, nrx3, -1)
         WRITE(stdout, '(4x,4x,"inveps_{GG}(q,w) = ", 2f9.5)'), drhoscfs(nl(ig),1) + 1.d0
         WRITE(6,"")


!        HL- Here generate an entire row of the screened coulomb matrix. 
!        qg2 = (g(1,ig)+xq(1))**2 + (g(2,ig)+xq(2))**2 + (g(3,ig)+xq(3))**2
!        do igp = 1, ngm
!        if (qg2 > 1.d-8) then
!        scrcoul(ig, igp, iw)  =  fpi * e2 * drhoscfs( nl(igp), 1) / (tpiba2 * qg2)
!        endif
!        enddo

         rec_code=20
         CALL write_rec('done_drhod',irr,0.0_DP,-1000,.false.,npe,&
                        drhoscfs)
<<<<<<< .mine
         enddo !iw
         !if (ig > 1) stop 
=======
        enddo !iw
        stop
         if (ig.gt.25) stop 

>>>>>>> .r273
         ! HL-   From here on Analytical Continuation to the real axis.
ENDDO !on G
        IF (okvan) THEN
           DEALLOCATE (int3)
           IF (okpaw) DEALLOCATE (int3_paw)
           IF (noncolin) DEALLOCATE(int3_nc)
        ENDIF
        tcpu = get_clock ('GW')
  DEALLOCATE (drhoscfs)
  CALL stop_clock ('gwqscf')
RETURN
END SUBROUTINE coulomb