!-----------------------------------------------------------------------
  ! Copyright (C) 2010-2015 Henry Lambert, Feliciano Giustino
  ! 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 do_stern()
  USE io_global,  ONLY : stdout, ionode_id, meta_ionode
  USE kinds,      ONLY : DP
  USE disp,       ONLY : nqs, num_k_pts, xk_kpoints, w_of_q_start
  USE gwsigma,    ONLY : sigma_c_st
  USE gwsymm,     ONLY : ngmunique, ig_unique, use_symm, sym_friend, sym_ig
  USE control_gw, ONLY : done_bands, reduce_io, recover, tmp_dir_gw,&
                          ext_restart, bands_computed, bands_computed, nbnd_occ, lgamma,&
                          do_q0_only, solve_direct, tinvert, lrpa, do_epsil
  USE freq_gw,    ONLY : nfs
  USE units_gw,   ONLY : lrcoul, iuncoul
  USE klist,      ONLY : lgauss
  USE mp_global,  ONLY : inter_image_comm, intra_image_comm, &
                         my_image_id, nimage, root_image
  USE mp,              ONLY : mp_sum, mp_barrier
  USE mp_world,             ONLY : mpime
  USE noncollin_module, ONLY : noncolin, nspin_mag

IMPLICIT NONE

  INTEGER :: iq, ik, ig, igstart, igstop, ios, iq1, iq2
  COMPLEX(DP), ALLOCATABLE :: scrcoul_g(:,:,:,:)
  LOGICAL :: do_band, do_iq, setup_pw, exst, do_matel
  COMPLEX(DP), ALLOCATABLE :: eps_m(:)

  ALLOCATE ( scrcoul_g( sigma_c_st%ngmt, sigma_c_st%ngmt, nfs, nspin_mag))
  ALLOCATE ( ig_unique( sigma_c_st%ngmt) )
  ALLOCATE ( sym_ig(sigma_c_st%ngmt))
  ALLOCATE ( sym_friend(sigma_c_st%ngmt))

  do_iq=.TRUE.
  setup_pw = .TRUE.
  do_band  = .TRUE.
  do_matel = .TRUE.

  if(lgauss) write(stdout, '(//5x,"SYSTEM IS METALLIC")')
  if(.not.do_epsil) then
      iq1 = w_of_q_start
      iq2 = nqs
  else
  ! In case we want to trace a line through the brillouin zone
  ! or get the screening for a particular grid q points (i.e. coulomb matel).
      iq1 = w_of_q_start
      iq2 = num_k_pts
  endif
    
  do iq = iq1, iq2
!Perform head of dielectric matrix calculation.
     CALL start_clock ('epsilq')
     if (iq.eq.1) THEN
        ALLOCATE(eps_m(nfs))
        eps_m(:) = dcmplx(0.0d0,0.0d0)
    !    if(my_image_id.eq.0) THEN
    !      scrcoul_g(:,:,:,:) = (0.0d0, 0.0d0)
    !      CALL prepare_q0(do_band, do_iq, setup_pw, iq)
    !      do_matel = .FALSE.
    !      CALL run_nscf(do_band, do_matel, iq)
    !      CALL initialize_gw()
    !      CALL coulomb_q0G0(iq, eps_m)
    !      CALL clean_pw_gw(iq)
    !    endif
    !    CALL mp_barrier(inter_image_comm)
    endif
    scrcoul_g(:,:,:,:) = (0.0d0, 0.0d0)
    CALL prepare_q(do_band, do_iq, setup_pw, iq)
    do_matel = .FALSE.
    CALL run_nscf(do_band, do_matel, iq)
    CALL initialize_gw()
    if(use_symm) THEN
      write(6,'("")')
      write(6,'(5x, "SYMMETRIZING COULOMB Perturbations")')
      write(6,'("")')
      CALL stern_symm()
    ELSE
      ngmunique = sigma_c_st%ngmt
      do ig = 1, sigma_c_st%ngmt
         ig_unique(ig) = ig
      enddo
    endif
       if(nimage.gt.1) then
          CALL para_img(ngmunique, igstart, igstop)
       else
          igstart = 1
          igstop = ngmunique
       endif
       write(6, '(5x, "iq ",i4, " igstart ", i4, " igstop ", i4)')iq, igstart, igstop
       CALL coulomb(iq, igstart, igstop, scrcoul_g)
       if(nimage.gt.1) THEN
          CALL mp_sum(scrcoul_g, inter_image_comm)
       endif
!Only the meta_image should write to file
       if (meta_ionode) THEN
         CALL unfold_w(scrcoul_g,iq)
         if(solve_direct.and.tinvert) write(1000+mpime, '("UNFOLDING, INVERTING, WRITING W")')
         if(solve_direct.and.tinvert) CALL invert_epsilon(scrcoul_g, iq, eps_m)
         CALL davcio(scrcoul_g, lrcoul, iuncoul, iq, +1, ios)
       endif
       if(allocated(eps_m)) deallocate(eps_m)
       call mp_barrier(inter_image_comm)
       CALL clean_pw_gw(iq)
       if(do_q0_only) GOTO 126
       CALL print_clock ('epsilq')
       CALL stop_clock ('epsilq')
  enddo

126 continue
   write(stdout, '("Finished Calculating Screened Coulomb")')
   deallocate( scrcoul_g )
   deallocate( ig_unique )
   deallocate( sym_ig )
   deallocate( sym_friend )
end SUBROUTINE do_stern