calcTDFandDOS Subroutine

  subroutine calcTDFandDOS(pw90_boltzwann, dis_manifold, pw90_dos, kpt_latt, postw90_oper, &
                           pw90_band_deriv_degen, pw90_spin, ws_region, print_output, &
                           wannier_data, ws_distance, wigner_seitz, HH_R, SS_R, u_matrix, &
                           v_matrix, eigval, real_lattice, TDF, TDFEnergyArray, cell_volume, &
                           scissors_shift, mp_grid, num_bands, num_kpts, num_wann, &
                           num_valence_bands, num_elec_per_state, effective_model, &
                           have_disentangled, spin_decomp, seedname, stdout, timer, error, comm)
    !================================================!
    !! This routine calculates the Transport Distribution Function $$\sigma_{ij}(\epsilon)$$ (TDF)
    !! in units of 1/hbar^2 * eV*fs/angstrom, and possibly the DOS.
    !!
    !! The TDFEnergyArray must be already allocated and initialized with the
    !! energies in eV before calling this routine.
    !!
    !!  The TDF array must be already allocated with dimensions 6 * size(TDFEnergyArray) * ndim, before calling
    !! this routine, where ndim=1 if spin_decomp==false, or ndim=3 if spin_decomp==true. This is not checked.
    !!
    !!  If run in parallel, at the end each processor will have a copy of the full TDF array
    !!
    !!  We assume that the TDFEnergyArray is uniformely spaced and that it has at least
    !!       two elements (no checks are performed on this).
    !!
    !! Note that the order of indices of TDF is different w.r.t. the DOS array (the energy is not the first but
    !!       the second index)
    !!
    !!  If the input flag pw90_boltzwann_bandshift is set to .true., the code will also shift the
    !!       conduction bands by a given amount, as defined by the pw90_boltzwann_bandshift_energyshift
    !!       and pw90_boltzwann_bandshift_firstband input flags.
    !!
    !================================================!

    use w90_constants, only: dp
    use w90_comms, only: comms_bcast, w90_comm_type, mpirank
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_utility, only: utility_recip_lattice_base
    use w90_get_oper, only: get_HH_R, get_SS_R
    use w90_types, only: print_output_type, wannier_data_type, dis_manifold_type, &
                         ws_region_type, ws_distance_type, timer_list_type
    use w90_postw90_types, only: pw90_boltzwann_type, pw90_spin_mod_type, &
                                 pw90_band_deriv_degen_type, pw90_dos_mod_type, pw90_oper_read_type, wigner_seitz_type
    use w90_readwrite, only: w90_readwrite_get_smearing_type
    use w90_wan_ham, only: wham_get_eig_deleig

    implicit none

    ! arguments
    type(pw90_boltzwann_type), intent(in) :: pw90_boltzwann
    type(dis_manifold_type), intent(in) :: dis_manifold
    type(pw90_dos_mod_type), intent(in) :: pw90_dos
    real(kind=dp), intent(in) :: kpt_latt(:, :)
    type(pw90_band_deriv_degen_type), intent(in) :: pw90_band_deriv_degen
    type(pw90_oper_read_type), intent(in) :: postw90_oper
    type(pw90_spin_mod_type), intent(in) :: pw90_spin
    type(print_output_type), intent(in) :: print_output
    type(ws_region_type), intent(in) :: ws_region
    type(w90_comm_type), intent(in) :: comm
    type(wannier_data_type), intent(in) :: wannier_data
    type(wigner_seitz_type), intent(inout) :: wigner_seitz
    type(ws_distance_type), intent(inout) :: ws_distance
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error

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

    real(kind=dp), intent(out) :: TDF(:, :, :) ! (coordinate,Energy,spin)
    !! The TDF(i,EnIdx,spin) output array, where:
    !!        - i is an index from 1 to 6 giving the component of the symmetric tensor
    !!          $$ Sigma_{ij}(\eps) $$,
    !!          where 1=xx, 2=xy, 3=yy, 4=xz, 5=yz, 6=zz
    !!          as defined by the module constants XX, XY, ...
    !!          (in this way the mapping (i,j) -> i+((j-1)*j)/2 is satisfied for the packed storage of the
    !!          upper triangle [i<=j]).
    !!        - EnIdx is the index of the energies; the corresponding energy is given by
    !!          TDFEnergyArray(EndIdx) array (in eV).
    !!        - Spin may be only 1 if spin_decomp=.false. If it is instead true, 1 contains the total TDF,
    !!          2 the spin-up component and 3 the spin-up component
    real(kind=dp), intent(in) :: TDFEnergyArray(:)
    !! TDFEnergyArray The array with the energies for which the TDF is calculated, in eV
    ! Comments:
    ! issue warnings if going outside of the energy window
    ! check that we actually get hbar*velocity in eV*angstrom
    real(kind=dp), intent(in) :: eigval(:, :)
    real(kind=dp), intent(in) :: real_lattice(3, 3)
    real(kind=dp), intent(in) :: cell_volume
    real(kind=dp), intent(in) :: scissors_shift

    complex(kind=dp), intent(in) :: v_matrix(:, :, :), u_matrix(:, :, :)
    complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :) !  <0n|r|Rm>
    complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :) ! <0n|sigma_x,y,z|Rm>

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

    ! local variables
    real(kind=dp) :: recip_lattice(3, 3), volume
    real(kind=dp) :: kpt(3), orig_kpt(3)
    integer :: loop_tot, loop_x, loop_y, loop_z, ierr

    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)

    real(kind=dp), allocatable :: DOS_EnergyArray(:)
    real(kind=dp), allocatable :: DOS_k(:, :), TDF_k(:, :, :)
    real(kind=dp), allocatable :: DOS_all(:, :)
    real(kind=dp) :: kweight
    integer :: ndim, DOS_NumPoints, i, j, k, EnIdx

    character(len=20) :: numfieldsstr
    integer :: boltzdos_unit, NumPtsRefined

    real(kind=dp), parameter :: SPACING_THRESHOLD = 1.e-3
    real(kind=dp) :: min_spacing, max_spacing

    integer :: my_node_id, num_nodes
    logical :: on_root = .false.

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

    if (print_output%iprint > 0 .and. (print_output%timing_level > 0)) call io_stopwatch_start('calcTDF', timer)
    if (print_output%iprint > 0) then
      if (pw90_boltzwann%calc_also_dos) then
        write (stdout, '(3X,A)') "Calculating Transport Distribution function (TDF) and DOS..."
      else
        write (stdout, '(3X,A)') "Calculating Transport Distribution function (TDF)..."
      end if
    end if

    ! I call once the routine to calculate the Hamiltonian in real-space <0n|H|Rm>

    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, &
                  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, postw90_oper, 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

    ! Some initial checks
    if (size(TDF, 1) /= 6 .or. size(TDF, 2) /= size(TDFEnergyArray) .or. size(TDF, 3) /= ndim) then
      call set_error_input(error, 'Wrong size for the TDF array in calcTDF', comm)
      return
    end if

    ! I zero the TDF array before starting
    TDF = 0._dp

    allocate (TDF_k(6, size(TDFEnergyArray), ndim), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating TDF_k in calcTDF', comm)
      return
    end if

    allocate (HH(num_wann, num_wann), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating HH in calcTDF', 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 calcTDF', 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 calcTDF', comm)
      return
    end if

    DOS_NumPoints = int(floor((pw90_boltzwann%dos_energy_max - pw90_boltzwann%dos_energy_min)/pw90_boltzwann%dos_energy_step)) + 1
    if (DOS_NumPoints .eq. 1) DOS_NumPoints = 2
    allocate (DOS_EnergyArray(DOS_NumPoints), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating DOS_EnergyArray in calcTDF', comm)
      return
    end if
    do i = 1, DOS_NumPoints
      DOS_EnergyArray(i) = pw90_boltzwann%dos_energy_min + real(i - 1, dp)*pw90_boltzwann%dos_energy_step
    end do

    allocate (DOS_k(size(DOS_EnergyArray), ndim), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating DOS_k in calcTDF', comm)
      return
    end if
    allocate (DOS_all(size(DOS_EnergyArray), ndim), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating DOS_all in calcTDF', comm)
      return
    end if
    dos_all = 0.0_dp

    ! I open the output files
    if (pw90_boltzwann%calc_also_dos .and. on_root) then
      open (newunit=boltzdos_unit, file=trim(seedname)//'_boltzdos.dat')
    end if

    if (pw90_boltzwann%calc_also_dos .and. on_root .and. (print_output%iprint > 1)) then
      write (stdout, '(5X,A)') "Smearing for DOS: "
      if (pw90_boltzwann%dos_smearing%use_adaptive) then
        write (stdout, '(7X,A)') trim(w90_readwrite_get_smearing_type(pw90_boltzwann%dos_smearing%type_index))//", adaptive"
      else
        if (pw90_boltzwann%dos_smearing%fixed_width/(DOS_EnergyArray(2) - DOS_EnergyArray(1)) < &
            min_smearing_binwidth_ratio) then
          write (stdout, '(7X,A)') "Unsmeared (use smearing width larger than bin width to smear)"
        else
          write (stdout, '(7X,A,G18.10)') trim(w90_readwrite_get_smearing_type(pw90_boltzwann%dos_smearing%type_index))// &
            ", non-adaptive, width (eV) =", pw90_boltzwann%dos_smearing%fixed_width
        end if
      end if
    end if
    if (pw90_boltzwann%calc_also_dos .and. pw90_boltzwann%dos_smearing%use_adaptive .and. &
        (pw90_boltzwann%dos_smearing%fixed_width .ne. 0._dp) .and. on_root) then
      write (stdout, '(5X,A)') "*** WARNING! pw90_boltzwann_dos_smr_fixed_en_width ignored since you chose"
      write (stdout, '(5X,A)') "             an adaptive smearing."
    end if

    if (on_root .and. (print_output%iprint > 1)) then
      if (pw90_boltzwann%tdf_smearing%fixed_width/(TDFEnergyArray(2) - TDFEnergyArray(1)) &
          < min_smearing_binwidth_ratio) then
        write (stdout, '(5X,A)') "Smearing for TDF: "
        write (stdout, '(7X,A)') "Unsmeared (use smearing width larger than bin width to smear)"
      else
        write (stdout, '(5X,A)') "Smearing for TDF: "
        write (stdout, '(7X,A,G18.10)') &
          trim(w90_readwrite_get_smearing_type(pw90_boltzwann%tdf_smearing%type_index))//", non-adaptive, width (eV) =", &
          pw90_boltzwann%tdf_smearing%fixed_width
      end if
    end if

    if (on_root) then
      write (stdout, '(5X,A,I0,A,I0,A,I0)') "k-grid used for band interpolation in BoltzWann: ", &
        pw90_boltzwann%kmesh%mesh(1), 'x', pw90_boltzwann%kmesh%mesh(2), 'x', pw90_boltzwann%kmesh%mesh(3)
      write (stdout, '(5X,A,I1)') "Number of electrons per state: ", num_elec_per_state
      write (stdout, '(5X,A,G18.10)') "Relaxation time (fs): ", pw90_boltzwann%relax_time
      if (print_output%iprint > 1) then
        write (stdout, '(5X,A,G18.10)') "Energy step for TDF (eV): ", pw90_boltzwann%tdf_energy_step
      end if
    end if
    kweight = 1.0_dp/real(PRODUCT(pw90_boltzwann%kmesh%mesh), kind=dp)

    if (pw90_boltzwann%bandshift .and. on_root) then
      write (stdout, '(5X,A,I0,A,G18.10,A)') "Shifting energy bands with index >= ", pw90_boltzwann%bandshift_firstband, " by ", &
        pw90_boltzwann%bandshift_energyshift, " eV."
    end if

    call utility_recip_lattice_base(real_lattice, recip_lattice, volume)
    NumPtsRefined = 0
    min_spacing = 1.e10_dp ! very large initial value
    max_spacing = 0.e0_dp
    ! I loop over all kpoints
    do loop_tot = my_node_id, PRODUCT(pw90_boltzwann%kmesh%mesh) - 1, num_nodes

      ! I get the coordinates for the x,y,z components starting from a single loop variable
      ! (which is better for parallelization purposes)
      ! Important! This works only if loop_tot starts from ZERO and ends with
      !            PRODUCT(pw90_boltzwann_kmesh)-1, so be careful when parallelizing
      loop_x = loop_tot/(pw90_boltzwann%kmesh%mesh(2)*pw90_boltzwann%kmesh%mesh(3))
      loop_y = (loop_tot - loop_x*(pw90_boltzwann%kmesh%mesh(2)*pw90_boltzwann%kmesh%mesh(3)))/pw90_boltzwann%kmesh%mesh(3)
      loop_z = loop_tot - loop_x*(pw90_boltzwann%kmesh%mesh(2)*pw90_boltzwann%kmesh%mesh(3)) - loop_y*pw90_boltzwann%kmesh%mesh(3)

      ! kpt(i) is in in the [0,d-1]/d range, with d=pw90_boltzwann_kmesh(i)
      kpt(1) = (real(loop_x, dp)/real(pw90_boltzwann%kmesh%mesh(1), dp))
      kpt(2) = (real(loop_y, dp)/real(pw90_boltzwann%kmesh%mesh(2), dp))
      kpt(3) = (real(loop_z, dp)/real(pw90_boltzwann%kmesh%mesh(3), dp))

      ! Here I get the band energies and the velocities
      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, num_valence_bands, effective_model, &
                               have_disentangled, seedname, stdout, timer, error, comm)
      if (allocated(error)) return

      call dos_get_levelspacing(del_eig, pw90_boltzwann%kmesh%mesh, levelspacing_k, num_wann, recip_lattice)

      ! Here I apply a scissor operator to the conduction bands, if required in the input
      if (pw90_boltzwann%bandshift) then
        eig(pw90_boltzwann%bandshift_firstband:) = eig(pw90_boltzwann%bandshift_firstband:) + pw90_boltzwann%bandshift_energyshift
      end if

      call TDF_kpt(pw90_boltzwann, ws_region, pw90_spin, wannier_data, ws_distance, wigner_seitz, &
                   HH_R, SS_R, del_eig, eig, TDFEnergyArray, kpt, real_lattice, TDF_k, mp_grid, &
                   num_wann, num_elec_per_state, spin_decomp, error, comm)
      if (allocated(error)) return

      ! As above, the sum of TDF_k * kweight amounts to calculate
      ! spin_degeneracy * V_cell/(2*pi)^3 * \int_BZ d^3k
      ! so that we divide by the cell_volume (in Angstrom^3) to have
      ! the correct integral
      TDF = TDF + TDF_k*kweight/cell_volume

      ! DOS part !

      if (pw90_boltzwann%calc_also_dos) then
        if (pw90_boltzwann%dos_smearing%use_adaptive) then

          ! This may happen if at least one band has zero derivative (along all three directions)
          ! Then I substitute this point with its 8 neighbors (+/- 1/4 of the spacing with the next point on the grid
          ! on each of the three directions)
          min_spacing = min(min_spacing, minval(abs(levelspacing_k)))
          max_spacing = max(max_spacing, maxval(abs(levelspacing_k)))
          if (any(abs(levelspacing_k) < SPACING_THRESHOLD)) then
            orig_kpt = kpt
            NumPtsRefined = NumPtsRefined + 1
            do i = -1, 1, 2
              do j = -1, 1, 2
                do k = -1, 1, 2
                  kpt = orig_kpt + &
                        (/real(i, kind=dp)/real(pw90_boltzwann%kmesh%mesh(1), dp)/4._dp, &
                          real(j, kind=dp)/real(pw90_boltzwann%kmesh%mesh(2), dp)/4._dp, &
                          real(k, kind=dp)/real(pw90_boltzwann%kmesh%mesh(3), dp)/4._dp/)
                  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, &
                                           num_valence_bands, effective_model, have_disentangled, &
                                           seedname, stdout, timer, error, comm)
                  if (allocated(error)) return

                  call dos_get_levelspacing(del_eig, pw90_boltzwann%kmesh%mesh, levelspacing_k, &
                                            num_wann, recip_lattice)
                  call dos_get_k(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_boltzwann%dos_smearing, error, comm, &
                                 levelspacing_k=levelspacing_k)
                  if (allocated(error)) return

                  ! I divide by 8 because I'm substituting a point with its 8 neighbors
                  dos_all = dos_all + dos_k*kweight/8.
                end do
              end do
            end do
          else
            call dos_get_k(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_boltzwann%dos_smearing, error, comm, levelspacing_k=levelspacing_k)
            if (allocated(error)) return

            dos_all = dos_all + dos_k*kweight
          end if
        else
          call dos_get_k(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_boltzwann%dos_smearing, error, comm)
          if (allocated(error)) return

          ! This sum multiplied by kweight amounts to calculate
          ! spin_degeneracy * V_cell/(2*pi)^3 * \int_BZ d^3k
          ! So that the DOS will be in units of 1/eV, normalized so that
          ! \int_{-\infty}^{\infty} DOS(E) dE = Num.Electrons
          dos_all = dos_all + dos_k*kweight
        end if
      end if

    end do

    ! I sum the results of the calculation for the DOS on the root node only
    ! (I have to print the results only)
    if (pw90_boltzwann%calc_also_dos) then
      call comms_reduce(DOS_all(1, 1), size(DOS_all), 'SUM', error, comm)
      if (allocated(error)) return
      call comms_reduce(NumPtsRefined, 1, 'SUM', error, comm)
      if (allocated(error)) return
      call comms_reduce(min_spacing, 1, 'MIN', error, comm)
      if (allocated(error)) return
      call comms_reduce(max_spacing, 1, 'MAX', error, comm)
      if (allocated(error)) return
    end if
    ! I sum the results of the calculation on all nodes, and I store them on all
    ! nodes (because for the following, each node will do a different calculation,
    ! each of which will require the whole knowledge of the TDF array)
    call comms_allreduce(TDF(1, 1, 1), size(TDF), 'SUM', error, comm)
    if (allocated(error)) return

    if (pw90_boltzwann%calc_also_dos .and. on_root) then
      write (boltzdos_unit, '(A)') "# Written by the BoltzWann module of the Wannier90 code."
      write (boltzdos_unit, '(A)') "# The first column."
      if (pw90_boltzwann%dos_smearing%use_adaptive) then
        write (boltzdos_unit, '(A)') '# The second column is the adaptively-smeared DOS'
        write (boltzdos_unit, '(A)') '# (see Yates et al., PRB 75, 195121 (2007)'
        if (spin_decomp) then
          write (boltzdos_unit, '(A)') '# The third column is the spin-up projection of the DOS'
          write (boltzdos_unit, '(A)') '# The fourth column is the spin-down projection of the DOS'
        end if
        write (boltzdos_unit, '(A,1X,G14.6)') '# Smearing coefficient: ', pw90_boltzwann%dos_smearing%adaptive_prefactor
        write (boltzdos_unit, '(A,I0,A,I0)') '# Number of points refined: ', NumPtsRefined, &
          ' out of ', product(pw90_boltzwann%kmesh%mesh)
        write (boltzdos_unit, '(A,G18.10,A,G18.10,A)') '# (Min spacing: ', min_spacing, &
          ', max spacing: ', max_spacing, ')'
      else
        if (pw90_boltzwann%dos_smearing%fixed_width/(DOS_EnergyArray(2) - DOS_EnergyArray(1)) < min_smearing_binwidth_ratio) then
          write (boltzdos_unit, '(A)') '# The second column is the unsmeared DOS.'
        else
          write (boltzdos_unit, '(A,G14.6,A)') '# The second column is the DOS for a fixed smearing of ', &
            pw90_boltzwann%dos_smearing%fixed_width, ' eV.'
        end if
      end if
      write (boltzdos_unit, '(A,1X,G14.6)') '# Cell volume (ang^3): ', cell_volume

      write (boltzdos_unit, '(A)') '# Energy(eV) DOS [DOS DOS ...]'

      ! I save a string with the number of fields to print
      write (numfieldsstr, '(I0)') 1 + ndim
      do EnIdx = 1, size(DOS_EnergyArray)
        write (boltzdos_unit, '(1X,'//trim(numfieldsstr)//'G18.10)') &
          DOS_EnergyArray(EnIdx), dos_all(EnIdx, :)
      end do
    end if

    if (on_root .and. (print_output%timing_level > 0)) call io_stopwatch_stop('calcTDF', timer)
    if (on_root) then
      if (pw90_boltzwann%calc_also_dos) then
        write (stdout, '(3X,A)') "TDF and DOS calculated."
      else
        write (stdout, '(3X,A)') "TDF calculated."
      end if
    end if
    if (on_root) write (stdout, *)

    if (on_root .and. pw90_boltzwann%calc_also_dos) then
      close (boltzdos_unit)
      if (print_output%iprint > 1) write (stdout, '(3X,A)') "DOS written on the "//trim(seedname)//"_boltzdos.dat file."
    end if

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

  end subroutine calcTDFandDOS