!
! Copyright (C) 2001-2003 PWSCF 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 .
!
!--------------------------------------------------------------------
program dos
  !--------------------------------------------------------------------
  !
  ! Input (namelist &inputpp ... &end):                Default value
  !
  !    prefix        prefix of input file produced by pw.x  'pwscf'
  !                    (wavefunctions are not needed)
  !    outdir        directory containing the input file    ./
  !    ngauss1,      gaussian broadening parameters         0
  !     degauss1       if absent, read from file            0.d0
  !    Emin, Emax    min, max energy (eV) for DOS plot      band extrema
  !    DeltaE        energy grid step (eV)                  none
  !    fildos        output file containing DOS(E)          dos.out
  !
  ! In order to use tetrahedron method, the input file must
  ! contain information on tetrahedra (option ltetra=.true.)
  ! if degauss1, ngauss1 are specified they override what is
  ! specified in the input file wrt summation method
  !
  USE io_global,  ONLY : stdout, ionode_id
  USE io_files,   ONLY : nd_nmbr, prefix, tmp_dir
  USE kinds, ONLY : DP
  USE klist,      ONLY : xk, wk, degauss, ngauss, lgauss, nks, nkstot
  USE ktetra,     ONLY : ntetra, tetra, ltetra
  USE wvfct,      ONLY : nbnd, et
  USE lsda_mod,   ONLY : nspin
#ifdef __PARA
  use para,       ONLY : me, mypool
  use mp,         ONLY : mp_bcast
#endif
  implicit none
  character(len=80) :: fildos
  character(len=256) :: outdir
  real(kind=DP) :: E, DOSofE (2), DOSint, Elw, Eup, DeltaE, Emin, Emax, &
       degauss1
  integer :: nrot, ik, n, ndos, ngauss1, ios
  namelist /inputpp/ outdir, prefix, fildos, degauss1,ngauss1,&
       Emin, Emax, DeltaE
  logical :: minus_q

!darek
  integer :: i,j,k,l
  integer :: natwf, nkpts, nbnds
  integer,allocatable :: wfno(:),wfnoat(:),wfc(:),wfm(:),wfl(:),indeks(:,:,:),indeks0(:,:,:)
  character(3),allocatable :: wfat(:),at(:)
  real,allocatable :: kpts(:,:),ebnds(:,:),&
                      psi(:,:,:),psi0(:,:,:),amp(:),&
                      proj1(:),proj2(:),&
                      proj_s(:),proj_p(:),proj_d(:),proj_f(:)
  real :: inter1,inter2
  logical,allocatable :: maska_s(:),maska_p(:),maska_d(:),maska_f(:)
!darek                                        
  
  !
  call start_postproc (nd_nmbr)
#ifdef __PARA
  if (me == 1 .and. mypool == 1) then
#endif
  !
  !   set default values for variables in namelist
  !
  outdir='./'
  prefix ='pwscf'
  fildos =' '
  Emin   =-1000000.
  Emax   = 1000000.
  DeltaE = 0.01
  ngauss1= 0
  degauss1=0.d0
  !
  read (5, inputpp, err=200, iostat=ios )
200 call errore('dos','reading inputpp namelist',abs(ios))
  !
  tmp_dir = trim(outdir)
#ifdef __PARA
  end if
  !
  ! ... Broadcast variables
  !
  CALL mp_bcast( tmp_dir, ionode_id )
  CALL mp_bcast( prefix, ionode_id )
#endif
  !
  call read_file
  !
#ifdef __PARA
  if (me == 1) then
#endif
  if (nks.ne.nkstot) call errore ('dos', 'not implemented', 1)
  !
  if (degauss1.ne.0.d0) then
     degauss=degauss1
     ngauss =ngauss1
     WRITE( stdout,'(/5x,"Gaussian broadening (read from input): ",&
          &        "ngauss,degauss=",i4,f12.6/)') ngauss,degauss
     ltetra=.false.
     lgauss=.true.
  else if (ltetra) then
     WRITE( stdout,'(/5x,"Tetrahedra used"/)')
  else if (lgauss) then
     WRITE( stdout,'(/5x,"Gaussian broadening (read from file): ",&
          &        "ngauss,degauss=",i4,f12.6/)') ngauss,degauss
  else
     if (degauss1.eq.0.d0) call errore('dos','I need a gaussian broadening!',1)
  end if

!darek
!
!reading information from proj.help file
!
	if (ltetra==.true.) then
	    open(10,FILE='proj.help')	
	    1000 FORMAT(5x,"state #",i3,": atom ",i3," (",a3,"), wfc ",i2," (l=",i1," m=",i2,")")
	    read(10,'(I55555)') natwf
	    allocate(wfno(natwf),wfnoat(natwf),wfc(natwf),wfm(natwf),wfl(natwf),wfat(natwf),at(natwf))
	    do i=1, natwf 
		read(10,1000) wfno(i),wfnoat(i),wfat(i),wfc(i),wfl(i),wfm(i)
	    enddo
	    read(10,'(I5,I5)') nkpts,nbnds
	    allocate(kpts(nkpts,3),ebnds(nkpts,nbnds),indeks(nkpts,nbnds,natwf),psi(nkpts,nbnds,natwf))
	    allocate(indeks0(nkpts,nbnds,natwf),psi0(nkpts,nbnds,natwf),amp(natwf),proj1(natwf),proj2(natwf))
	    allocate(maska_s(natwf),maska_p(natwf),maska_d(natwf),maska_f(natwf))
	    allocate(proj_s(wfnoat(natwf)),proj_p(wfnoat(natwf)),proj_d(wfnoat(natwf)),proj_f(wfnoat(natwf)))
	    do i=1,nkpts 
		read(10,'(" k = ",3f14.10)') (kpts (i, j) , j = 1, 3)
		do j=1,nbnds
		    read(10,'(5x,"e = ",f14.10," eV")') ebnds(i,j)
		    do k=1, natwf
			read(10,'(I3,1X,F7.5)') indeks(i,j,k),psi(i,j,k)
		    enddo    
		enddo
	    enddo
	    300 close(10)
!
!sorting
!

	    do i=1,nkpts
		do j=1,nbnds
		    do k=1,natwf
			do l=1,natwf
			    if (indeks(i,j,l)==k) then
				indeks0(i,j,k)=indeks(i,j,l)
				psi0(i,j,k)=psi(i,j,l)
			    endif
			enddo
		    enddo
		enddo
	    enddo
        
	    psi=psi0
	    indeks=indeks0
!
!opening output files
!
	    k=1; j=1
	    at(1)=wfat(1)
	    do i=1, natwf
		if (k/=wfnoat(i)) then
		    k=wfnoat(i)
		    j=j+1
		    at(j)=wfat(i)
		endif
	    enddo
	    open(10,FILE='proj_dos.up')
	    write(10,'(A8)') '# DOS up'
	    write(10,'(A7,<4*wfnoat(natwf)>A12,A12)') &
		    '# E[eV]',(trim(at(i))//'_s',i=1,j),(trim(at(i))//'_p',i=1,j),&
		    (trim(at(i))//'_d',i=1,j),(trim(at(i))//'_f',i=1,j),'Inter'

	    open(20,FILE='proj_dos.dn')
	    write(20,'(A8)') '# DOS dn'
	    write(20,'(A7,<4*wfnoat(natwf)>A12,A12)') &
		    '# E[eV]',(trim(at(i))//'_s',i=1,j),(trim(at(i))//'_p',i=1,j),&
		    (trim(at(i))//'_d',i=1,j),(trim(at(i))//'_f',i=1,j),'Inter'
	endif
!darek

  !
  ! find band extrema
  !
  Elw = et (1, 1)
  Eup = et (nbnd, 1)
  do ik = 2, nks
     Elw = min (Elw, et (1, ik) )
     Eup = max (Eup, et (nbnd, ik) )
  enddo
  if (degauss.ne.0.d0) then
     Eup = Eup + 3d0 * degauss
     Elw = Elw - 3d0 * degauss
  endif
  !
  Emin=max(Emin/13.6058,Elw)
  Emax=min(Emax/13.6058,Eup)
  DeltaE = DeltaE / 13.6058
  ndos = nint ( (Emax - Emin) / DeltaE+0.500001)
  DOSint = 0.0
  !
  open (unit = 4, file = fildos, status = 'unknown', form = 'formatted')
  if (nspin.eq.1) then
     write(4,'("#  E (eV)   dos(E)     Int dos(E)")')
  else
     write(4,'("#  E (eV)   dosup(E)     dosdw(E)   Int dos(E)")')
  end if

  do n= 1, ndos
     E = Emin + (n - 1) * DeltaE
     if (ltetra) then
        call dos_t(et,nspin,nbnd, nks,ntetra,tetra, E, DOSofE)
     else
        call dos_g(et,nspin,nbnd, nks,wk,degauss,ngauss, E, DOSofE)
     endif
     if (nspin.eq.1) then
        DOSint = DOSint + DOSofE (1) * DeltaE
        write (4, '(f7.3,2e12.4)') E * 13.6058, DOSofE(1)/13.6058, DOSint
     else
        DOSint = DOSint + (DOSofE (1) + DOSofE (2) ) * DeltaE
        write (4, '(f7.3,3e12.4)') E * 13.6058, DOSofE/13.6058, DOSint
     endif
     
!darek
!
!for given energy e projection on atomic wfc's are calculated 
!
	if (ltetra==.true.) then
	    e=e*13.6058
	    proj1=0.0; proj2=0.0
	    inter1=0.0; inter2=0.0
	    l=0
	    do i=1,nkpts
		do j=1,nbnds-1
		    if ((ebnds(i,j)<=e).and.(ebnds(i,j+1)>=e)) then
			amp(:)=psi(i,j,:)+(e-ebnds(i,j))/(ebnds(i,j+1)-ebnds(i,j))*(psi(i,j+1,:)-psi(i,j,:))
			proj1=proj1+DOSofE(1)*amp/13.6058
			proj2=proj2+DOSofE(2)*amp/13.6058
			inter1=inter1+DOSofE(1)*(1-sum(amp))/13.6058
			inter2=inter2+DOSofE(2)*(1-sum(amp))/13.6058
			l=l+1
		    endif
		enddo
	    enddo
	    if (l>0) then
		proj1=proj1/l; inter1=inter1/l
		proj2=proj2/l; inter2=inter2/l
	    else
		proj1=0.0; inter1=0.0
		proj2=0.0; inter2=0.0
	    endif
!
!this creates DOS for all atoms and for s,p,d,f states
!
!DOS up
	    proj_s=0; proj_p=0; proj_d=0; proj_f=0
	    do i=1, wfnoat(natwf)
		do j=1,natwf 
		    maska_s(j)=(wfnoat(j)==i.and.wfl(j)==0)
		    maska_p(j)=(wfnoat(j)==i.and.wfl(j)==1)
		    maska_d(j)=(wfnoat(j)==i.and.wfl(j)==2)
		    maska_f(j)=(wfnoat(j)==i.and.wfl(j)==3)
		enddo
		proj_s(i)=sum(proj1,maska_s)
		proj_p(i)=sum(proj1,maska_p)
		proj_d(i)=sum(proj1,maska_d)
		proj_f(i)=sum(proj1,maska_f)
	    enddo
!
!writing DOS to output file in folowing order: e,atom1_s,...,atomN_s,atom1_p,...,atomN_p,...,inter
!where inter means nonlocalized part of DOS 
!
	    write(10,'(F7.3,<4*wfnoat(natwf)>E12.5,E12.5)') &
		    e,(proj_s(i),i=1,wfnoat(natwf)),(proj_p(i),i=1,wfnoat(natwf)),&
		    (proj_d(i),i=1,wfnoat(natwf)),(proj_f(i),i=1,wfnoat(natwf)),inter1

!DOS dn    
	    proj_s=0; proj_p=0; proj_d=0; proj_f=0
	    do i=1, wfnoat(natwf)
		do j=1,natwf 
		    maska_s(j)=(wfnoat(j)==i.and.wfl(j)==0)
		    maska_p(j)=(wfnoat(j)==i.and.wfl(j)==1)
		    maska_d(j)=(wfnoat(j)==i.and.wfl(j)==2)
		    maska_f(j)=(wfnoat(j)==i.and.wfl(j)==3)
		enddo
		proj_s(i)=sum(proj2,maska_s)
		proj_p(i)=sum(proj2,maska_p)
		proj_d(i)=sum(proj2,maska_d)
		proj_f(i)=sum(proj2,maska_f)
	    enddo
    
	    write(20,'(F7.3,<4*wfnoat(natwf)>E12.5,E12.5)') &
		    e,(proj_s(i),i=1,wfnoat(natwf)),(proj_p(i),i=1,wfnoat(natwf)),&
		    (proj_d(i),i=1,wfnoat(natwf)),(proj_f(i),i=1,wfnoat(natwf)),inter2
	endif
!darek     
          
  enddo

  close (unit = 4)
!darek  
	    close(10)
	    close(20)
	    deallocate(wfno,wfnoat,wfc,wfm,wfl,wfat,at)
	    deallocate(kpts,ebnds,indeks,indeks0,psi,psi0,amp,proj1,proj2,maska_s,maska_p,maska_d,maska_f)
!darek
  
#ifdef __PARA
  end if
#endif
  call stop_pp
end program dos