dos_main Subroutine

  subroutine dos_main(pw90_berry, dis_manifold, pw90_dos, kpoint_dist, kpt_latt, pw90_oper_read, &
                      pw90_band_deriv_degen, pw90_spin, ws_region, w90_system, print_output, &
                      wannier_data, ws_distance, wigner_seitz, HH_R, SS_R, u_matrix, v_matrix, &
                      eigval, real_lattice, scissors_shift, mp_grid, num_bands, num_kpts, &
                      num_wann, effective_model, have_disentangled, spin_decomp, seedname, stdout, &
                      timer, error, comm)
    !================================================!
    !
    !! 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_comms, only: comms_reduce, w90_comm_type, mpirank, mpisize
    use w90_postw90_common, only: pw90common_fourier_R_to_k
    use w90_postw90_types, only: pw90_dos_mod_type, pw90_berry_mod_type, &
                                 pw90_band_deriv_degen_type, pw90_spin_mod_type, pw90_oper_read_type, wigner_seitz_type, &
                                 kpoint_dist_type
    use w90_types, only: print_output_type, wannier_data_type, dis_manifold_type, &
                         ws_region_type, w90_system_type, ws_distance_type, timer_list_type
    use w90_get_oper, only: get_HH_R, get_SS_R
    use w90_io, only: io_date, io_stopwatch_start, io_stopwatch_stop
    use w90_utility, only: utility_diagonalize, utility_recip_lattice_base
    use w90_wan_ham, only: wham_get_eig_deleig

    implicit none

    ! arguments
    type(pw90_berry_mod_type), intent(in)        :: pw90_berry
    type(dis_manifold_type), intent(in)          :: dis_manifold
    type(pw90_dos_mod_type), intent(in)          :: pw90_dos
    type(kpoint_dist_type), intent(in)           :: kpoint_dist
    type(pw90_oper_read_type), intent(in)        :: pw90_oper_read
    type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
    type(pw90_spin_mod_type), intent(in)         :: pw90_spin
    type(ws_region_type), intent(in)             :: ws_region
    type(w90_system_type), intent(in)            :: w90_system
    type(print_output_type), intent(in)          :: print_output
    type(wannier_data_type), intent(in)          :: wannier_data
    type(ws_distance_type), intent(inout)        :: ws_distance
    type(wigner_seitz_type), intent(inout)       :: wigner_seitz
    type(timer_list_type), intent(inout)         :: timer
    type(w90_comm_type), intent(in)               :: comm
    type(w90_error_type), allocatable, intent(out) :: error

    complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :)
    complex(kind=dp), intent(in) :: u_matrix(:, :, :), v_matrix(:, :, :)

    real(kind=dp), intent(in) :: eigval(:, :), real_lattice(3, 3)
    real(kind=dp), intent(in) :: scissors_shift
    real(kind=dp), intent(in) :: kpt_latt(:, :)

    integer, intent(in) :: mp_grid(3)
    integer, intent(in) :: num_bands, num_kpts, num_wann
    integer, intent(in) :: stdout

    character(len=50), intent(in) :: seedname
    logical, intent(in) :: have_disentangled
    logical, intent(in) :: spin_decomp
    logical, intent(in) :: effective_model

    ! local variables
    ! '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) :: recip_lattice(3, 3), volume

    integer :: i, loop_x, loop_y, loop_z, loop_tot, ifreq
    integer :: dos_unit, ndim, ierr
    integer :: my_node_id, num_nodes
    integer :: num_freq !! Number of sampling points

    real(kind=dp), allocatable :: dos_k(:, :)
    real(kind=dp), allocatable :: dos_all(:, :)
    real(kind=dp) :: kweight, kpt(3), omega
    real(kind=dp), allocatable :: dos_energyarray(:)
    real(kind=dp) :: del_eig(num_wann, 3)
    real(kind=dp) :: eig(num_wann), levelspacing_k(num_wann)
    real(kind=dp) :: d_omega !! Step between energies

    complex(kind=dp), allocatable :: HH(:, :)
    complex(kind=dp), allocatable :: delHH(:, :, :)
    complex(kind=dp), allocatable :: UU(:, :)
    logical :: on_root = .false.

    my_node_id = mpirank(comm)
    num_nodes = mpisize(comm)
    if (my_node_id == 0) on_root = .true.

    num_freq = nint((pw90_dos%energy_max - pw90_dos%energy_min)/pw90_dos%energy_step) + 1
    if (num_freq == 1) num_freq = 2
    d_omega = (pw90_dos%energy_max - pw90_dos%energy_min)/(num_freq - 1)

    allocate (dos_energyarray(num_freq), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating dos_energyarray in dos subroutine', comm)
      return
    end if

    do ifreq = 1, num_freq
      dos_energyarray(ifreq) = pw90_dos%energy_min + real(ifreq - 1, dp)*d_omega
    end do

    allocate (HH(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating HH in dos', comm)
      return
    end if
    allocate (delHH(num_wann, num_wann, 3), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating delHH in dos', comm)
      return
    end if
    allocate (UU(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating UU in dos', comm)
      return
    end if

    call get_HH_R(dis_manifold, kpt_latt, print_output, wigner_seitz, HH_R, u_matrix, v_matrix, &
                  eigval, real_lattice, scissors_shift, num_bands, num_kpts, num_wann, &
                  w90_system%num_valence_bands, effective_model, have_disentangled, seedname, ws_distance, ws_region, &
                  stdout, timer, error, comm)

    if (allocated(error)) return

    if (spin_decomp) then
      ndim = 3
      call get_SS_R(dis_manifold, kpt_latt, print_output, pw90_oper_read, SS_R, v_matrix, eigval, &
                    wigner_seitz, ws_distance, ws_region, num_bands, num_kpts, num_wann, &
                    have_disentangled, seedname, stdout, timer, error, comm)

      if (allocated(error)) return

    else
      ndim = 1
    end if

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

    if (print_output%iprint > 0) then

      if (print_output%timing_level > 1) call io_stopwatch_start('dos', timer)

!       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 (pw90_dos%num_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, pw90_dos%num_project
          write (stdout, '(5x,a,2x,i3)') 'n =', pw90_dos%project(i)
        end do
      end if

      write (stdout, '(/,5x,a,f9.4,a,f9.4,a)') &
        'Energy range: [', pw90_dos%energy_min, ',', pw90_dos%energy_max, '] eV'

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

      write (stdout, '(/,/,1x,a20,3(i0,1x))') 'Interpolation grid: ', &
        pw90_dos%kmesh%mesh(1:3)

    end if

    dos_all = 0.0_dp

    call utility_recip_lattice_base(real_lattice, recip_lattice, volume)
    if (pw90_berry%wanint_kpoint_file) then
      !
      ! Unlike for optical properties, this should always work for the DOS
      !
      if (print_output%iprint > 0) 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, kpoint_dist%num_int_kpts_on_node(my_node_id)
        kpt(:) = kpoint_dist%int_kpts(:, loop_tot)
        if (pw90_dos%smearing%use_adaptive) then
          call wham_get_eig_deleig(dis_manifold, kpt_latt, pw90_band_deriv_degen, ws_region, &
                                   print_output, wannier_data, ws_distance, wigner_seitz, delHH, &
                                   HH, HH_R, u_matrix, UU, v_matrix, del_eig, eig, eigval, kpt, &
                                   real_lattice, scissors_shift, mp_grid, num_bands, num_kpts, &
                                   num_wann, w90_system%num_valence_bands, effective_model, &
                                   have_disentangled, seedname, stdout, timer, error, comm)
          if (allocated(error)) return

          call dos_get_levelspacing(del_eig, pw90_dos%kmesh%mesh, levelspacing_k, num_wann, &
                                    recip_lattice)
          call dos_get_k(w90_system%num_elec_per_state, ws_region, kpt, dos_energyarray, eig, &
                         dos_k, num_wann, wannier_data, real_lattice, mp_grid, pw90_dos, &
                         spin_decomp, pw90_spin, ws_distance, wigner_seitz, HH_R, SS_R, &
                         pw90_dos%smearing, error, comm, levelspacing_k=levelspacing_k, UU=UU)
          if (allocated(error)) return

        else
          call pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, HH, &
                                         HH_R, kpt, real_lattice, mp_grid, 0, num_wann, error, comm)
          if (allocated(error)) return

          call utility_diagonalize(HH, num_wann, eig, UU, error, comm)
          if (allocated(error)) return

          call dos_get_k(w90_system%num_elec_per_state, ws_region, kpt, dos_energyarray, eig, &
                         dos_k, num_wann, wannier_data, real_lattice, mp_grid, pw90_dos, &
                         spin_decomp, pw90_spin, ws_distance, wigner_seitz, HH_R, SS_R, &
                         pw90_dos%smearing, error, comm, UU=UU)
          if (allocated(error)) return

        end if
        dos_all = dos_all + dos_k*kpoint_dist%weight(loop_tot)
      end do

    else

      if (print_output%iprint > 0) write (stdout, '(/,1x,a)') 'Sampling the full BZ'

      kweight = 1.0_dp/real(PRODUCT(pw90_dos%kmesh%mesh), kind=dp)
      do loop_tot = my_node_id, PRODUCT(pw90_dos%kmesh%mesh) - 1, num_nodes
        loop_x = loop_tot/(pw90_dos%kmesh%mesh(2)*pw90_dos%kmesh%mesh(3))
        loop_y = (loop_tot - loop_x*(pw90_dos%kmesh%mesh(2) &
                                     *pw90_dos%kmesh%mesh(3)))/pw90_dos%kmesh%mesh(3)
        loop_z = loop_tot - loop_x*(pw90_dos%kmesh%mesh(2)*pw90_dos%kmesh%mesh(3)) &
                 - loop_y*pw90_dos%kmesh%mesh(3)
        kpt(1) = real(loop_x, dp)/real(pw90_dos%kmesh%mesh(1), dp)
        kpt(2) = real(loop_y, dp)/real(pw90_dos%kmesh%mesh(2), dp)
        kpt(3) = real(loop_z, dp)/real(pw90_dos%kmesh%mesh(3), dp)
        if (pw90_dos%smearing%use_adaptive) then
          call wham_get_eig_deleig(dis_manifold, kpt_latt, pw90_band_deriv_degen, ws_region, &
                                   print_output, wannier_data, ws_distance, wigner_seitz, delHH, &
                                   HH, HH_R, u_matrix, UU, v_matrix, del_eig, eig, eigval, kpt, &
                                   real_lattice, scissors_shift, mp_grid, num_bands, num_kpts, &
                                   num_wann, w90_system%num_valence_bands, effective_model, &
                                   have_disentangled, seedname, stdout, timer, error, comm)
          if (allocated(error)) return

          call dos_get_levelspacing(del_eig, pw90_dos%kmesh%mesh, levelspacing_k, num_wann, &
                                    recip_lattice)
          call dos_get_k(w90_system%num_elec_per_state, ws_region, kpt, dos_energyarray, eig, &
                         dos_k, num_wann, wannier_data, real_lattice, mp_grid, pw90_dos, &
                         spin_decomp, pw90_spin, ws_distance, wigner_seitz, HH_R, SS_R, &
                         pw90_dos%smearing, error, comm, levelspacing_k=levelspacing_k, UU=UU)
          if (allocated(error)) return

        else
          call pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, HH, &
                                         HH_R, kpt, real_lattice, mp_grid, 0, num_wann, error, comm)
          if (allocated(error)) return

          call utility_diagonalize(HH, num_wann, eig, UU, error, comm)
          if (allocated(error)) return

          call dos_get_k(w90_system%num_elec_per_state, ws_region, kpt, dos_energyarray, eig, &
                         dos_k, num_wann, wannier_data, real_lattice, mp_grid, pw90_dos, &
                         spin_decomp, pw90_spin, ws_distance, wigner_seitz, HH_R, SS_R, &
                         pw90_dos%smearing, error, comm, UU=UU)
          if (allocated(error)) return

        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', error, comm)
    if (allocated(error)) return

    if (print_output%iprint > 0) then
      write (stdout, '(1x,a)') 'Output data files:'
      write (stdout, '(/,3x,a)') trim(seedname)//'-dos.dat'
      open (newunit=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, :)
      end do
      close (dos_unit)
      if (print_output%timing_level > 1) call io_stopwatch_stop('dos', timer)
    end if

    deallocate (HH, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating HH in dos_main', comm)
      return
    end if
    deallocate (delHH, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating delHH in dos_main', comm)
      return
    end if
    deallocate (UU, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating UU in dos_main', comm)
      return
    end if

  end subroutine dos_main