!-------------------------------------------------------------------------
      subroutine efield(rhog)
!     subroutine efield(rhog,tfor,tau)
!-----------------------------------------------------------------------
!     rhog input : electronic charge in g space till ng
!
!     note that the average electric field is set equal to zero
!     (G=0 component). The true average is recovered by hand 
!     on the basis of the self-consistent electric field 
!     given as input in the cpe run 


      use ions
      use gvec
      use gvecs
      use parm
      use parms
      use elct
      use cnst
      use ener
      use control
      use work1

      implicit none

      logical                        tfor
      real(kind=8)                   rhor(nnr,nspin),vr(nnr),sigma,sigmaion
      real(kind=8)                   tau(3,nax,nsp),rhogions(ng)
      complex(kind=8)                rhog(ng,nspin)
      integer                        ia,iss,isup,isdw,ix,ig,ir,i,j,k,is,l,m
      complex(kind=8)                fp,fm,ci
      complex(kind=8), pointer::     v(:)
      complex(kind=8), allocatable:: rhotemp(:),vtemp(:)

      ci=(0.,1.)

      v => wrk1
      allocate(rhotemp(ng))
      allocate(vtemp(ng))

      sigma=1.88976 ! 1 angstrom

      if(nspin.eq.1)then
         iss=1
         do i=1,ng
            rhotemp(i)=rhog(i,iss)
         end do
      else
         isup=1
         isdw=2
         do i=1,ng
            rhotemp(i)=rhog(i,isup)+rhog(i,isdw)
         end do
      endif

!      if (tfor) then
!
!        include ions as gaussian distributions
!
!        I AM NOT SURE THAT THE CUTOFF CAN AFFORD GAUSSIANS OF 0.2 ...
!
!        sigmaion=0.2
!        do ig=1,ng    ! to be checked !!!!
!           rhogions(ig)=1./omega*exp(-g(ig)*tpiba2*sigmaion**2./2.)
!        end do
!
!        do is=1,nsp
!           do ia=1,na(is)
!             do ig=1,ng
!               rhotemp(ig)=rhotemp(ig)+rhogions(ig)*&
!                  &zv(is)*exp(-ci*tpiba*dot_product(gx(ig,:),tau(:,ia,is)))
!             end do
!           end do
!        end do
!
!        check wheter the integrated charge is as expected
!
!        if (ng0.eq.2) write(6,*) omega*rhotemp(1)
!
!      end if


      if (ng0.eq.2) vtemp(1)=(0.,0.)
      do ig=ng0,ng
         vtemp(ig)=rhotemp(ig)*fpi/(tpiba2*g(ig))
!        field along x
         vtemp(ig)=vtemp(ig)*gx(ig,1)*ci*tpiba
!        field along y
!         vtemp(ig)=vtemp(ig)*gx(ig,2)*ci*tpiba
!        field along z
!         vtemp(ig)=vtemp(ig)*gx(ig,3)*ci*tpiba
!        gaussian filter
         vtemp(ig)=vtemp(ig)*exp(-(tpiba2*g(ig))*sigma**2./2)

      end do

!     fourier transform of hartree field to r-space (dense grid)

      call zero(2*nnr,v)
      do ig=1,ng
         v(nm(ig))=conjg(vtemp(ig))
         v(np(ig))=vtemp(ig)
      end do

!     v(g) --> v(r)

      call invfft(v,nr1,nr2,nr3,nr1x,nr2x,nr3x)

      do ir=1,nnr
         vr(ir)=real(v(ir))
      end do

      call totrho(vr)

      deallocate(vtemp)
      deallocate(rhotemp)

      return
      end