!-*- mode: F90 -*-!
!------------------------------------------------------------!
! This file is distributed as part of the Wannier90 code and !
! under the terms of the GNU General Public License. See the !
! file `LICENSE' in the root directory of the Wannier90      !
! distribution, or http://www.gnu.org/copyleft/gpl.txt       !
!                                                            !
! The webpage of the Wannier90 code is www.wannier.org       !
!                                                            !
! The Wannier90 code is hosted on GitHub:                    !
!                                                            !
! https://github.com/wannier-developers/wannier90            !
!------------------------------------------------------------!

module w90_dos
  !! Compute Density of States
  use w90_constants, only : dp

  implicit none

  private

  public :: dos_main, dos_get_levelspacing, dos_get_k

  integer       :: num_freq
  !! Number of sampling points
  real(kind=dp) :: d_omega
  !! Step between energies

contains

  !=========================================================!
  !                   PUBLIC PROCEDURES                     ! 
  !=========================================================!

  subroutine dos_main
    !=======================================================!
    !                                                       !
    !! Computes the electronic density of states. Can 
    !! resolve into up-spin and down-spin parts, project 
    !! onto selected Wannier orbitals, and use adaptive 
    !! broadening, as in PRB 75, 195121 (2007) [YWVS07].
    !                                                       !
    !=======================================================!

    use w90_io, only             : io_error,io_file_unit,io_date,io_stopwatch,&
                                   seedname,stdout
    use w90_comms, only          : on_root,num_nodes,my_node_id,comms_reduce
    use w90_postw90_common, only : num_int_kpts_on_node,int_kpts,weight,&
                                   pw90common_fourier_R_to_k
    use w90_parameters, only     : num_wann,dos_energy_min,dos_energy_max,&
                                   dos_energy_step,timing_level,&
                                   wanint_kpoint_file,dos_kmesh,&
                                   dos_smr_index,dos_adpt_smr,dos_adpt_smr_fac,&
                                   dos_adpt_smr_max, spin_decomp,&
                                   dos_smr_fixed_en_width,&
                                   dos_project,num_dos_project
    use w90_get_oper, only       : get_HH_R,get_SS_R,HH_R
    use w90_wan_ham, only        : wham_get_eig_deleig
    use w90_utility, only        : utility_diagonalize

    ! 'dos_k' contains contrib. from one k-point, 
    ! 'dos_all' from all nodes/k-points (first summed on one node and 
    ! then reduced (i.e. summed) over all nodes)
    !
    real(kind=dp), allocatable :: dos_k(:,:)
    real(kind=dp), allocatable :: dos_all(:,:)

    real(kind=dp)    :: kweight,kpt(3),omega
    integer          :: i,loop_x,loop_y,loop_z,loop_tot,ifreq
    integer          :: dos_unit,ndim, ierr
    real(kind=dp), dimension(:), allocatable :: dos_energyarray

    complex(kind=dp), allocatable :: HH(:,:)
    complex(kind=dp), allocatable :: delHH(:,:,:)
    complex(kind=dp), allocatable :: UU(:,:)
    real(kind=dp) :: del_eig(num_wann,3)
    real(kind=dp) :: eig(num_wann), levelspacing_k(num_wann)

    num_freq=nint((dos_energy_max-dos_energy_min)/dos_energy_step)+1
    if(num_freq==1) num_freq=2
    d_omega=(dos_energy_max-dos_energy_min)/(num_freq-1)

    allocate(dos_energyarray(num_freq),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating dos_energyarray in '&
         //'dos subroutine')
    do ifreq=1,num_freq
       dos_energyarray(ifreq) = dos_energy_min + real(ifreq-1,dp)*d_omega
    end do

    allocate(HH(num_wann,num_wann),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating HH in dos')
    allocate(delHH(num_wann,num_wann,3),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating delHH in dos')
    allocate(UU(num_wann,num_wann),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating UU in dos')    

    call get_HH_R
    if(spin_decomp) then
       ndim=3
       call get_SS_R
    else
       ndim=1
    end if

    allocate(dos_k(num_freq,ndim))
    allocate(dos_all(num_freq,ndim))

    if(on_root) then

       if (timing_level>1) call io_stopwatch('dos',1)

!       write(stdout,'(/,1x,a)') '============'
!       write(stdout,'(1x,a)')   'Calculating:'
!       write(stdout,'(1x,a)')   '============'

       write(stdout,'(/,/,1x,a)')&
            'Properties calculated in module  d o s'
       write(stdout,'(1x,a)')&
            '--------------------------------------'

       if(num_dos_project==num_wann) then
          write(stdout,'(/,3x,a)') '* Total density of states (_dos)'
       else
          write(stdout,'(/,3x,a)')&
               '* Density of states projected onto selected WFs (_dos)'
          write(stdout,'(3x,a)') 'Selected WFs |Rn> are:'
          do i=1,num_dos_project
             write(stdout,'(5x,a,2x,i3)') 'n =',dos_project(i)
          enddo
       endif

       write(stdout,'(/,5x,a,f9.4,a,f9.4,a)')&
            'Energy range: [',dos_energy_min,',',dos_energy_max,'] eV'

       write(stdout,'(/,5x,a,(f6.3,1x))')&
            'Adaptive smearing width prefactor: ',&
            dos_adpt_smr_fac

       write(stdout, '(/,/,1x,a20,3(i0,1x))') 'Interpolation grid: ',&
            dos_kmesh(1:3) 

    end if

    dos_all=0.0_dp

    if(wanint_kpoint_file) then
       !
       ! Unlike for optical properties, this should always work for the DOS
       !
       if(on_root)  write(stdout,'(/,1x,a)') 'Sampling the irreducible BZ only'

       ! Loop over k-points on the irreducible wedge of the Brillouin zone,
       ! read from file 'kpoint.dat'
       !
       do loop_tot=1,num_int_kpts_on_node(my_node_id)
          kpt(:)=int_kpts(:,loop_tot)
          if (dos_adpt_smr) then
             call wham_get_eig_deleig(kpt,eig,del_eig,HH,delHH,UU)
             call dos_get_levelspacing(del_eig,dos_kmesh,levelspacing_k)
             call dos_get_k(kpt,dos_energyarray,eig,dos_k,&
                  smr_index=dos_smr_index,&
                  adpt_smr_fac=dos_adpt_smr_fac,&
                  adpt_smr_max=dos_adpt_smr_max,&
                  levelspacing_k=levelspacing_k,&
                  UU=UU)
          else
             call pw90common_fourier_R_to_k(kpt,HH_R,HH,0) 
             call utility_diagonalize(HH,num_wann,eig,UU) 
             call dos_get_k(kpt,dos_energyarray,eig,dos_k,&
                  smr_index=dos_smr_index,&
                  smr_fixed_en_width=dos_smr_fixed_en_width,&
                  UU=UU)
          end if
          dos_all=dos_all+dos_k*weight(loop_tot)
       end do

    else
       
       if (on_root) write(stdout,'(/,1x,a)') 'Sampling the full BZ'
       
       kweight = 1.0_dp / real(PRODUCT(dos_kmesh),kind=dp)
       do loop_tot=my_node_id,PRODUCT(dos_kmesh)-1,num_nodes
          loop_x= loop_tot/(dos_kmesh(2)*dos_kmesh(3))
          loop_y=(loop_tot-loop_x*(dos_kmesh(2)&
               *dos_kmesh(3)))/dos_kmesh(3)
          loop_z=loop_tot-loop_x*(dos_kmesh(2)*dos_kmesh(3))&
                -loop_y*dos_kmesh(3)
          kpt(1)=real(loop_x,dp)/real(dos_kmesh(1),dp)
          kpt(2)=real(loop_y,dp)/real(dos_kmesh(2),dp)
          kpt(3)=real(loop_z,dp)/real(dos_kmesh(3),dp)
          if (dos_adpt_smr) then
             call wham_get_eig_deleig(kpt,eig,del_eig,HH,delHH,UU)
             call dos_get_levelspacing(del_eig,dos_kmesh,levelspacing_k)
             call dos_get_k(kpt,dos_energyarray,eig,dos_k,&
                  smr_index=dos_smr_index,&
                  adpt_smr_fac=dos_adpt_smr_fac,&
                  adpt_smr_max=dos_adpt_smr_max,&
                  levelspacing_k=levelspacing_k,&
                  UU=UU)
          else
             call pw90common_fourier_R_to_k(kpt,HH_R,HH,0) 
             call utility_diagonalize(HH,num_wann,eig,UU) 
             call dos_get_k(kpt,dos_energyarray,eig,dos_k,&
                  smr_index=dos_smr_index,&
                  smr_fixed_en_width=dos_smr_fixed_en_width,&
                  UU=UU)             
          end if
          dos_all=dos_all+dos_k*kweight
       end do

    end if

    ! Collect contributions from all nodes
    !
    call comms_reduce(dos_all(1,1),num_freq*ndim,'SUM')

    if(on_root) then
       write(stdout,'(1x,a)')     'Output data files:'
       write(stdout,'(/,3x,a)') trim(seedname)//'-dos.dat'
       dos_unit=io_file_unit()
       open(dos_unit,FILE=trim(seedname)//'-dos.dat',STATUS='UNKNOWN',&
            FORM='FORMATTED')
       do ifreq=1,num_freq
          omega=dos_energyarray(ifreq)
          write(dos_unit,'(4E16.8)') omega,dos_all(ifreq,:)
       enddo
       close(dos_unit)
       if (timing_level>1) call io_stopwatch('dos',2)
    end if

    deallocate(HH,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating HH in dos_main')
    deallocate(delHH,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating delHH in dos_main')
    deallocate(UU,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating UU in dos_main')

  end subroutine dos_main

  ! =========================================================================

! The next routine is commented. It should be working (apart for a
! missing broadcast at the very end, see comments there).  However,
! it should be debugged, and probably the best thing is to avoid to
! resample the BZ, but rather use the calculated DOS (maybe it can be
! something that is done at the end of the DOS routine?)
!~  subroutine find_fermi_level
!~    !==============================================!
!~    !                                              !
!~    ! Finds the Fermi level by integrating the DOS !
!~    !                                              !
!~    !==============================================!
!~
!~    use w90_io, only            : stdout,io_error
!~    use w90_comms
!~    use w90_postw90_common, only : max_int_kpts_on_node,num_int_kpts_on_node,&
!~         int_kpts,weight
!~    use w90_parameters, only    : nfermi,fermi_energy_list,&
!~         num_valence_bands,&
!~         num_wann,dos_num_points,dos_min_energy,&
!~         dos_max_energy,dos_energy_step,&
!~         wanint_kpoint_file
!~
!~#ifdef MPI 
!~    include 'mpif.h'
!~#endif
!~
!~    real(kind=dp) :: kpt(3),sum_max_node,sum_max_all,&
!~         sum_mid_node,sum_mid_all,emin,emax,emid,&
!~         emin_node(0:num_nodes-1),emax_node(0:num_nodes-1),&
!~         ef
!~    integer       :: loop_x,loop_y,loop_z,loop_kpt,loop_nodes,&
!~         loop_iter,ierr,num_int_kpts,ikp
!~
!~    real(kind=dp), allocatable :: eig_node(:,:)
!~    real(kind=dp), allocatable :: levelspacing_node(:,:)
!~
!~    if(on_root) write(stdout,'(/,a)') 'Finding the value of the Fermi level'
!~
!~    if(.not.wanint_kpoint_file) then
!~       !
!~       ! Already done in pw90common_wanint_get_kpoint_file if 
!~       ! wanint_kpoint_file=.true.
!~       !
!~       allocate(num_int_kpts_on_node(0:num_nodes-1))
!~       num_int_kpts=dos_num_points**3
!~       !
!~       ! Local k-point counter on each node (lazy way of doing it, there is
!~       ! probably a smarter way)
!~       !
!~       ikp=0
!~       do loop_kpt=my_node_id,num_int_kpts-1,num_nodes
!~          ikp=ikp+1
!~       end do
!~       num_int_kpts_on_node(my_node_id)=ikp
!~#ifdef MPI
!~       call MPI_reduce(ikp,max_int_kpts_on_node,1,MPI_integer,&
!~            MPI_MAX,0,MPI_COMM_WORLD,ierr)
!~#else
!~       max_int_kpts_on_node=ikp
!~#endif
!~       call comms_bcast(max_int_kpts_on_node,1)
!~    end if
!~
!~    allocate(eig_node(num_wann,max_int_kpts_on_node),stat=ierr)
!~    if (ierr/=0)&
!~         call io_error('Error in allocating eig_node in find_fermi_level')
!~    eig_node=0.0_dp
!~    allocate(levelspacing_node(num_wann,max_int_kpts_on_node),stat=ierr)
!~    if (ierr/=0)&
!~         call io_error('Error in allocating levelspacing in find_fermi_level')
!~    levelspacing_node=0.0_dp
!~
!~    if(wanint_kpoint_file) then
!~       if(on_root) write(stdout,'(/,1x,a)') 'Sampling the irreducible BZ only'
!~       do loop_kpt=1,num_int_kpts_on_node(my_node_id)
!~          kpt(:)=int_kpts(:,loop_kpt)
!~          call dos_get_eig_levelspacing_k(kpt,eig_node(:,loop_kpt),&
!~               levelspacing_node(:,loop_kpt))
!~       end do
!~    else
!~       if (on_root)&
!~            write(stdout,'(/,1x,a)') 'Sampling the full BZ (not using symmetry)'
!~       allocate(weight(max_int_kpts_on_node),stat=ierr)
!~       if (ierr/=0)&
!~            call io_error('Error in allocating weight in find_fermi_level')
!~       weight=0.0_dp
!~       ikp=0
!~       do loop_kpt=my_node_id,num_int_kpts-1,num_nodes
!~          ikp=ikp+1
!~          loop_x=loop_kpt/dos_num_points**2
!~          loop_y=(loop_kpt-loop_x*dos_num_points**2)/dos_num_points
!~          loop_z=loop_kpt-loop_x*dos_num_points**2-loop_y*dos_num_points
!~          kpt(1)=real(loop_x,dp)/dos_num_points
!~          kpt(2)=real(loop_y,dp)/dos_num_points
!~          kpt(3)=real(loop_z,dp)/dos_num_points
!~          weight(ikp)=1.0_dp/dos_num_points**3
!~          call dos_get_eig_levelspacing_k(kpt,eig_node(:,ikp),&
!~               levelspacing_node(:,ikp))
!~       end do
!~    end if
!~
!~    ! Find minimum and maximum band energies within projected subspace
!~    !
!~    emin_node(my_node_id)=&
!~         minval(eig_node(1,1:num_int_kpts_on_node(my_node_id)))
!~    emax_node(my_node_id)=&
!~         maxval(eig_node(num_wann,1:num_int_kpts_on_node(my_node_id)))
!~    if(.not.on_root) then
!~       call comms_send(emin_node(my_node_id),1,root_id)
!~       call comms_send(emax_node(my_node_id),1,root_id)
!~    else
!~       do loop_nodes=1,num_nodes-1
!~          call comms_recv(emin_node(loop_nodes),1,loop_nodes)
!~          call comms_recv(emax_node(loop_nodes),1,loop_nodes)
!~       end do
!~       emin=minval(emin_node)
!~       emax=maxval(emax_node)
!~    end if
!~    call comms_bcast(emin,1)
!~    call comms_bcast(emax,1)
!~
!~    ! Check that the Fermi level lies within the projected subspace
!~    !
!~    sum_max_node=count_states(emax,eig_node,levelspacing_node,&
!~         num_int_kpts_on_node(my_node_id))
!~#ifdef MPI
!~    call MPI_reduce(sum_max_node,sum_max_all,1,MPI_DOUBLE_PRECISION,&
!~         MPI_SUM,0,MPI_COMM_WORLD,ierr)
!~#else
!~    sum_max_all=sum_max_node
!~#endif
!~    if(on_root) then
!~       if(num_valence_bands>sum_max_all) then
!~          write(stdout,*) 'Something wrong in find_fermi_level:'
!~          write(stdout,*)&
!~               '   Fermi level does not lie within projected subspace'
!~          write(stdout,*) 'num_valence_bands= ',num_valence_bands
!~          write(stdout,*) 'sum_max_all= ',sum_max_all
!~          stop 'Stopped: see output file'
!~       end if
!~    end if
!~
!~    ! Now interval search for the Fermi level
!~    !
!~    do loop_iter=1,1000
!~       emid=(emin+emax)/2.0_dp
!~       sum_mid_node=count_states(emid,eig_node,levelspacing_node,&
!~            num_int_kpts_on_node(my_node_id))
!~#ifdef MPI
!~       call MPI_reduce(sum_mid_node,sum_mid_all,1,MPI_DOUBLE_PRECISION,&
!~            MPI_SUM,0,MPI_COMM_WORLD,ierr)
!~#else
!~       sum_mid_all=sum_mid_node
!~#endif
!~       ! This is needed because MPI_reduce only returns sum_mid_all to the 
!~       ! root (To understand: could we use MPI_Allreduce instead?)
!~       !
!~       call comms_bcast(sum_mid_all,1)
!~       if(abs(sum_mid_all-num_valence_bands) < 1.e-10_dp) then
!~          !
!~          ! NOTE: Here should assign a value to an entry in a fermi-level 
!~          !       vector. Then at the end average over adaptive smearing 
!~          !       widths and broadcast the result
!~          !
!~          ef=emid
!~          exit
!~       elseif((sum_mid_all-num_valence_bands) < -1.e-10_dp) then
!~          emin=emid
!~       else
!~          emax=emid
!~       end if
!~    end do
!~    
!~    nfermi=1
!~    allocate(fermi_energy_list(1))
!~    fermi_energy_list(1)=ef
!~    !~~ PROBABLY HERE YOU MAY WANT TO BROADCAST THE ABOVE TWO VARIABLES~!
!~    if(on_root) then
!~       write(stdout,*) ' '
!~       write(stdout,'(1x,a,f10.6,a)')&
!~            'Fermi energy = ',ef, ' eV'
!~       write(stdout,'(1x,a)')&
!~            '---------------------------------------------------------'
!~    end if
!~
!~  end subroutine find_fermi_level


  !> This subroutine calculates the contribution to the DOS of a single k point
  !> 
  !> \todo still to do: adapt spin_get_nk to read in input the UU rotation matrix
  !> 
  !> \note This routine simply provides the dos contribution of a given
  !>       point. This must be externally summed after proper weighting.
  !>       The weight factor (for a full BZ sampling with N^3 points) is 1/N^3 if we want
  !>       the final DOS to be normalized to the total number of electrons.
  !> \note The only factor that is included INSIDE this routine is the spin degeneracy
  !>       factor (=num_elec_per_state variable)
  !> \note The EnergyArray is assumed to be evenly spaced (and the energy spacing
  !>       is taken from EnergyArray(2)-EnergyArray(1))
  !> \note The routine is assuming that EnergyArray has at least two elements.
  !> \note The dos_k array must have dimensions size(EnergyArray) * ndim, where
  !>       ndim=1 if spin_decomp==false, or ndim=3 if spin_decomp==true. This is not checked. 
  !> \note If smearing/binwidth < min_smearing_binwidth_ratio, 
  !>       no smearing is applied (for that k point)
  !>
  !> \param kpt         the three coordinates of the k point vector whose DOS contribution we
  !>                    want to calculate (in relative coordinates)
  !> \param EnergyArray array with the energy grid on which to calculate the DOS (in eV)
  !>                    It must have at least two elements
  !> \param eig_k       array with the eigenvalues at the given k point (in eV)
  !> \param dos_k       array in which the contribution is stored. Three dimensions:
  !>                    dos_k(energyidx, spinidx), where:
  !>                    - energyidx is the index of the energies, corresponding to the one
  !>                      of the EnergyArray array; 
  !>                    - spinidx=1 contains the total dos; if if spin_decomp==.true., then
  !>                      spinidx=2 and spinidx=3 contain the spin-up and spin-down contributions to the DOS
  !> \param smr_index  index that tells the kind of smearing
  !> \param smr_fixed_en_width optional parameter with the fixed energy for smearing, in eV. Can be provided only if the
  !>                    levelspacing_k parameter is NOT given
  !> \param adpt_smr_fac optional parameter with the factor for the adaptive smearing. Can be provided only if the
  !>                    levelspacing_k parameter IS given
  !> \param levelspacing_k optional array with the level spacings, i.e. how much each level changes
  !>                    near a given point of the interpolation mesh, as given by the
  !>                    dos_get_levelspacing() routine
  !>                    If present: adaptive smearing
  !>                    If not present: fixed-energy-width smearing
  subroutine dos_get_k(kpt,EnergyArray,eig_k,dos_k,smr_index,&
       smr_fixed_en_width,adpt_smr_fac,adpt_smr_max,levelspacing_k,UU)
    use w90_io, only            : io_error
    use w90_constants, only     : dp, smearing_cutoff,min_smearing_binwidth_ratio
    use w90_utility, only       : utility_w0gauss
    use w90_parameters, only    : num_wann,spin_decomp,num_elec_per_state,&
                                  num_dos_project,dos_project
    use w90_spin, only          : spin_get_nk

    ! Arguments
    !
    real(kind=dp), dimension(3), intent(in)               :: kpt
    real(kind=dp), dimension(:), intent(in)               :: EnergyArray
    real(kind=dp), dimension(:), intent(in)               :: eig_k
    real(kind=dp), dimension(:,:), intent(out)            :: dos_k
    integer, intent(in)                                   :: smr_index
    real(kind=dp), intent(in),optional                    :: smr_fixed_en_width
    real(kind=dp), intent(in),optional                    :: adpt_smr_fac
    real(kind=dp), intent(in),optional                    :: adpt_smr_max
    real(kind=dp), dimension(:), intent(in),optional      :: levelspacing_k
    complex(kind=dp), dimension(:,:), intent(in),optional :: UU

    ! Adaptive smearing
    !
    real(kind=dp) :: eta_smr,arg

    ! Misc/Dummy
    !
    integer          :: i,j,loop_f,min_f,max_f, num_s_steps
    real(kind=dp)    :: rdum,spn_nk(num_wann),alpha_sq,beta_sq 
    real(kind=dp)    :: binwidth, r_num_elec_per_state
    logical          :: DoSmearing
   
    if (present(levelspacing_k)) then
       if (present(smr_fixed_en_width)) &
            call io_error('Cannot call doskpt with levelspacing_k and '&
            //'with smr_fixed_en_width parameters together')
       if (.not.(present(adpt_smr_fac))) &
            call io_error('Cannot call doskpt with levelspacing_k and '&
            //'without adpt_smr_fac parameter')
       if (.not.(present(adpt_smr_max))) &
            call io_error('Cannot call doskpt with levelspacing_k and '&
            //'without adpt_smr_max parameter')
    else
       if (present(adpt_smr_fac)) &
            call io_error('Cannot call doskpt without levelspacing_k and '&
            //'with adpt_smr_fac parameter')
       if (present(adpt_smr_max)) &
            call io_error('Cannot call doskpt without levelspacing_k and '&
            //'with adpt_smr_max parameter')
       if (.not.(present(smr_fixed_en_width))) &
            call io_error('Cannot call doskpt without levelspacing_k and '&
            //'without smr_fixed_en_width parameter')
    end if

    r_num_elec_per_state = real(num_elec_per_state,kind=dp)

    ! Get spin projections for every band
    !
    if(spin_decomp) call spin_get_nk(kpt,spn_nk)

    binwidth = EnergyArray(2) - EnergyArray(1)
    
    dos_k=0.0_dp
    do i=1,num_wann
       if(spin_decomp) then
          ! Contribution to spin-up DOS of Bloch spinor with component 
          ! (alpha,beta) with respect to the chosen quantization axis
          alpha_sq=(1.0_dp+spn_nk(i))/2.0_dp ! |alpha|^2
          ! Contribution to spin-down DOS 
          beta_sq=1.0_dp-alpha_sq ! |beta|^2 = 1 - |alpha|^2
       end if

       if (.not.present(levelspacing_k)) then
          eta_smr=smr_fixed_en_width
       else
          ! Eq.(35) YWVS07
          eta_smr=min(levelspacing_k(i)*adpt_smr_fac,adpt_smr_max)
!          eta_smr=max(eta_smr,min_smearing_binwidth_ratio) !! No: it would render the next if always false
       end if

       ! Faster optimization: I precalculate the indices
       if (eta_smr/binwidth < min_smearing_binwidth_ratio) then
          min_f= max(nint((eig_k(i) - EnergyArray(1))/&
               (EnergyArray(size(EnergyArray))-EnergyArray(1)) &
               * real(size(EnergyArray)-1,kind=dp)) + 1, 1)
          max_f= min(nint((eig_k(i) - EnergyArray(1))/&
               (EnergyArray(size(EnergyArray))-EnergyArray(1)) &
               * real(size(EnergyArray)-1,kind=dp)) + 1, size(EnergyArray))
          DoSmearing=.false.
       else      
          min_f= max(nint((eig_k(i) - smearing_cutoff * eta_smr - EnergyArray(1))/&
               (EnergyArray(size(EnergyArray))-EnergyArray(1)) &
               * real(size(EnergyArray)-1,kind=dp)) + 1, 1)
          max_f= min(nint((eig_k(i) + smearing_cutoff * eta_smr - EnergyArray(1))/&
               (EnergyArray(size(EnergyArray))-EnergyArray(1)) &
               * real(size(EnergyArray)-1,kind=dp)) + 1, size(EnergyArray))
          DoSmearing=.true.
       end if


       do loop_f=min_f, max_f
          ! kind of smearing read from input (internal smearing_index variable)
          if (DoSmearing) then
             arg=(EnergyArray(loop_f)-eig_k(i))/eta_smr
             rdum=utility_w0gauss(arg,smr_index)/eta_smr
          else
             rdum=1._dp/(EnergyArray(2)-EnergyArray(1))
          end if
          
          !
          ! Contribution to total DOS
          !
          if(num_dos_project==num_wann) then 
             !
             ! Total DOS (default): do not loop over j, to save time
             !
             dos_k(loop_f,1)=dos_k(loop_f,1)+rdum * r_num_elec_per_state
             ! [GP] I don't put num_elec_per_state here below: if we are 
             ! calculating the spin decomposition, we should be doing a 
             ! calcultation with spin-orbit, and thus num_elec_per_state=1!
             if(spin_decomp) then
                ! Spin-up contribution
                dos_k(loop_f,2)=dos_k(loop_f,2)+rdum*alpha_sq
                ! Spin-down contribution
                dos_k(loop_f,3)=dos_k(loop_f,3)+rdum*beta_sq
             end if
          else ! 0<num_dos_project<num_wann
             !
             ! Partial DOS, projected onto the WFs with indices 
             ! n=dos_project(1:num_dos_project)
             !
             do j=1,num_dos_project 
                dos_k(loop_f,1)=dos_k(loop_f,1)+rdum * r_num_elec_per_state&
                     *abs(UU(dos_project(j),i))**2
                if(spin_decomp) then
                   ! Spin-up contribution
                   dos_k(loop_f,2)=dos_k(loop_f,2)&
                        +rdum*alpha_sq*abs(UU(dos_project(j),i))**2
                   ! Spin-down contribution
                   dos_k(loop_f,3)=dos_k(loop_f,3)&
                        +rdum*beta_sq*abs(UU(dos_project(j),i))**2
                end if
             enddo
          endif
       enddo !loop_f
    end do !loop over bands

  end subroutine dos_get_k

  ! =========================================================================

  subroutine dos_get_levelspacing(del_eig,kmesh,levelspacing)
  !! This subroutine calculates the level spacing, i.e. how much the level changes
  !! near a given point of the interpolation mesh
    use w90_parameters, only: num_wann
    use w90_postw90_common, only : pw90common_kmesh_spacing
    
    real(kind=dp), dimension(num_wann,3), intent(in) :: del_eig
    !! Band velocities, already corrected when degeneracies occur
    integer, dimension(3), intent(in)                :: kmesh
    !! array of three integers, giving the number of k points along
    !! each of the three directions defined by the reciprocal lattice vectors
    real(kind=dp), dimension(num_wann), intent(out)  :: levelspacing
    !! On output, the spacing for each of the bands (in eV)

    real(kind=dp) :: Delta_k
    integer :: band
    

    Delta_k=pw90common_kmesh_spacing(kmesh)
    do band=1,num_wann
       levelspacing(band)=&
            sqrt(dot_product(del_eig(band,:),del_eig(band,:)))*Delta_k
    end do

  end subroutine dos_get_levelspacing

! Next routine is commented; it is the older version
!~  subroutine dos_get_eig_levelspacing_k(kpt,eig,levelspacing)
!~
!~    use w90_constants, only     : dp,cmplx_0,cmplx_i,twopi
!~    use w90_io, only            : io_error
!~    use w90_utility, only   : utility_diagonalize
!~    use w90_parameters, only    : num_wann,dos_num_points
!~    use w90_postw90_common, only : pw90common_fourier_R_to_k,pw90common_kmesh_spacing
!~    use w90_get_oper, only      : HH_R
!~    use w90_wan_ham, only   : wham_get_deleig_a
!~
!~    ! Arguments
!~    !
!~    real(kind=dp), intent(in)  :: kpt(3)
!~    real(kind=dp), intent(out) :: eig(num_wann)
!~    real(kind=dp), intent(out) :: levelspacing(num_wann)
!~
!~    complex(kind=dp), allocatable :: HH(:,:)
!~    complex(kind=dp), allocatable :: delHH(:,:,:)
!~    complex(kind=dp), allocatable :: UU(:,:)
!~
!~    ! Adaptive smearing
!~    !
!~    real(kind=dp) :: del_eig(num_wann,3),Delta_k
!~
!~    integer          :: i
!~
!~    allocate(HH(num_wann,num_wann))
!~    allocate(delHH(num_wann,num_wann,3))
!~    allocate(UU(num_wann,num_wann))
!~
!~    call pw90common_fourier_R_to_k(kpt,HH_R,HH,0) 
!~    call utility_diagonalize(HH,num_wann,eig,UU) 
!~    call pw90common_fourier_R_to_k(kpt,HH_R,delHH(:,:,1),1) 
!~    call pw90common_fourier_R_to_k(kpt,HH_R,delHH(:,:,2),2) 
!~    call pw90common_fourier_R_to_k(kpt,HH_R,delHH(:,:,3),3) 
!~    call wham_get_deleig_a(del_eig(:,1),eig,delHH(:,:,1),UU)
!~    call wham_get_deleig_a(del_eig(:,2),eig,delHH(:,:,2),UU)
!~    call wham_get_deleig_a(del_eig(:,3),eig,delHH(:,:,3),UU)
!~
!~    Delta_k=pw90common_kmesh_spacing(dos_num_points)
!~    do i=1,num_wann
!~       levelspacing(i)=&
!~            sqrt(dot_product(del_eig(i,:),del_eig(i,:)))*Delta_k
!~    end do
!~
!~  end subroutine dos_get_eig_levelspacing_k

!~  !=========================================================!
!~  !                   PRIVATE PROCEDURES                    ! 
!~  !=========================================================!
!~
!~
!~  function count_states(energy,eig,levelspacing,npts)
!~
!~    use w90_constants, only     : dp,cmplx_0,cmplx_i,twopi
!~    use w90_utility, only       : utility_wgauss
!~    use w90_postw90_common, only : weight
!~    use w90_parameters, only    : num_wann,dos_adpt_smr_fac
!~
!~    real(kind=dp) :: count_states
!~
!~    ! Arguments
!~    !
!~    real(kind=dp)                  :: energy
!~    real(kind=dp), dimension (:,:) :: eig
!~    real(kind=dp), dimension (:,:) :: levelspacing
!~    integer                        :: npts
!~
!~    ! Misc/Dummy
!~    !
!~    integer       :: loop_k,i
!~    real(kind=dp) :: sum,eta_smr,arg
!~
!~    count_states=0.0_dp
!~    do loop_k=1,npts
!~       sum=0.0_dp
!~       do i=1,num_wann
!~          eta_smr=levelspacing(i,loop_k)*dos_adpt_smr_fac
!~          arg=(energy-eig(i,loop_k))/eta_smr
!~          !
!~          ! For Fe and a 125x125x125 interpolation mesh, E_f=12.6306 with M-P
!~          ! smearing, and E_f=12.6512 with F-D smearing
!~          !
!~          !          sum=sum+utility_wgauss(arg,-99) ! Fermi-Dirac
!~          sum=sum+utility_wgauss(arg,1)    ! Methfessel-Paxton case
!~       end do
!~       count_states=count_states+weight(loop_k)*sum
!~    end do
!~
!~  end function count_states

end module w90_dos

