!-----------------------------------------------------------------------------
      subroutine eigfun (nfi,iprint,tbuff,trhor,trhow,c,irb,&
                        &eigrb,bec,rhovan,rhor,rhog,rhos,   &
                        &scale,instate,finstate)
!-----------------------------------------------------------------------------
!     eigenfunction evaluation
!     handle eigenfunction densities as rhor 


      use ions
      use gvec
      use gvecs
      use gvecb, only: ngb
      use cvan
      use parm
      use parms
      use elct
      use cnst
      use pseu
      use ener
      use control
      use work1
#ifdef PARA
      use para_mod
      include 'mpif.h'
#endif

      !implicit none

      real(kind=8)		bec(nhsa,n),rhovan(nat,nhx,nhx,nspin)
      real(kind=8)		rhor(nnr,nspin),rhos(nnrs,nspin),ftrue(nx)
      complex(kind=8)		eigrb(ngb,nas,nsp),c(ngw,nx),rhog(ng,nspin)
      integer			irb(3,nax,nsx),nfi,iprint,ig,iss,ir,i,is
      logical			tbuff,trhor,trhow
      integer           	istart,istop,i1,i2,i3,scale,nar1x,nar2x,nar3x
      integer			istate,instate,finstate
      real(kind=8)      	arg,eigx,rvec(3),rho(2),rhoca(nr1x,nr2x,nr3x,nspin),rhotot
      real(kind=8),allocatable::rhoave(:,:,:,:)

#ifdef PARA
      integer root, proc, ierr, displs(nproc), recvcount(nproc)
      real*8, allocatable:: rhodist(:,:)
      root=0
      allocate(rhodist(nr1x*nr2x*nr3x,nspin)) 
      do proc=1,nproc
        recvcount(proc) =  ncplane*npp(proc)
        if (proc.eq.1) then
           displs(proc)=0
        else
           displs(proc)=displs(proc-1) + recvcount(proc-1)
        end if
      end do
#endif

      nar1x=nr1x/scale
      nar2x=nr2x/scale
      nar3x=nr3x/scale
      allocate(rhoave(nar1x,nar2x,nar3x,nspin))


      write(6,*)
      write(6,*) '---> entering EIGFUN'
      write(6,*)
      if (finstate.lt.instate) call error('eigfun','(finstate.lt.instate)',0)


      ftrue=f

      do istate=instate,finstate

         write(6,*)
         write(6,*) '     > eigenfunction ',istate

         f=(/(0.,i=1,nx)/)
         f(istate)=1.

         call rhoofr (nfi,iprint,tbuff,trhor,trhow,c,irb,&
                     &eigrb,bec,rhovan,rhor,rhog,rhos)

#ifdef PARA
         do is=1,nspin

!          gather the charge density on the first node

           call mpi_barrier ( MPI_COMM_WORLD, ierr)
           call mpi_gatherv (rhor(1,is), recvcount(me), MPI_REAL8,        &
       &                     rhodist(:,is),recvcount, displs, MPI_REAL8,        &
       &                     root, MPI_COMM_WORLD, ierr)
           if (ierr.ne.0) call error('mpi_gatherv','ierr<>0',ierr)
         end do
#endif

!------- RESHAPE RHOR INTO A CARTESIAN MATRIX RHOCA -----------------------------

#ifdef PARA
         if (me.eq.1) then
#endif
           do iss=1,nspin
             do i2=1,nr2x
             do i3=1,nr3x
               istart=1+(i2-1)*nr1x+(i3-1)*nr1x*nr2x
               istop =      i2*nr1x+(i3-1)*nr1x*nr2x
#ifdef PARA
               rhoca(:,i2,i3,iss)=rhodist(istart:istop,iss)
#else
               rhoca(:,i2,i3,iss)=rhor(istart:istop,iss)
#endif
             end do
             end do 
           end do
#ifdef PARA
         endif 
#endif

!------- AVERAGE CARTESIAN DENSITY ----------------------------------------------

#ifdef PARA
      if (me.eq.1) then
#endif
         do iss=1,nspin
           do i1=1,nar1x
           do i2=1,nar2x
           do i3=1,nar3x
              rhoave(i1,i2,i3,iss)=float(scale)**(-float(3))*&
                  &sum(rhoca(scale*(i1-1)+1:scale*i1,&
                            &scale*(i2-1)+1:scale*i2,&
                            &scale*(i3-1)+1:scale*i3,iss))
           end do
           end do 
           end do
         end do
        
         do iss=1,nspin
           do i3=1,nar3x
           do i2=1,nar2x
           do i1=1,nar1x
             write(103,'(e20.10)') rhoave(i1,i2,i3,iss)
           end do
           end do
           end do
         end do
#ifdef PARA
      end if
#endif

      end do

!       reset true values 
      f=ftrue

      write(6,*)
      write(6,*) '<--- exiting EIGFUN'
      write(6,*)
      return
      end