calcTDFandDOS Subroutine

private subroutine calcTDFandDOS(TDF, TDFEnergyArray)

proc~~calctdfanddos~~UsesGraph proc~calctdfanddos calcTDFandDOS module~w90_get_oper w90_get_oper module~w90_get_oper->proc~calctdfanddos module~w90_wan_ham w90_wan_ham module~w90_wan_ham->proc~calctdfanddos module~w90_utility w90_utility module~w90_utility->proc~calctdfanddos module~w90_parameters w90_parameters module~w90_parameters->proc~calctdfanddos module~w90_constants w90_constants module~w90_constants->module~w90_get_oper module~w90_constants->module~w90_wan_ham module~w90_constants->module~w90_utility module~w90_constants->module~w90_parameters module~w90_io w90_io module~w90_constants->module~w90_io module~w90_io->module~w90_parameters
Help

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 boltz_bandshift is set to .true., the code will also shift the conduction bands by a given amount, as defined by the boltz_bandshift_energyshift and boltz_bandshift_firstband input flags.

Arguments

Type IntentOptional AttributesName
real(kind=dp), intent(out), dimension(:,:,:):: TDF

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), dimension(:):: TDFEnergyArray

TDFEnergyArray The array with the energies for which the TDF is calculated, in eV

Calls

proc~~calctdfanddos~~CallsGraph proc~calctdfanddos calcTDFandDOS proc~param_get_smearing_type param_get_smearing_type proc~calctdfanddos->proc~param_get_smearing_type proc~tdf_kpt TDF_kpt proc~calctdfanddos->proc~tdf_kpt interface~comms_allreduce comms_allreduce proc~calctdfanddos->interface~comms_allreduce get_ss_r get_ss_r proc~calctdfanddos->get_ss_r proc~io_error io_error proc~calctdfanddos->proc~io_error get_hh_r get_hh_r proc~calctdfanddos->get_hh_r proc~io_stopwatch io_stopwatch proc~calctdfanddos->proc~io_stopwatch proc~io_file_unit io_file_unit proc~calctdfanddos->proc~io_file_unit proc~wham_get_eig_deleig wham_get_eig_deleig proc~calctdfanddos->proc~wham_get_eig_deleig interface~comms_reduce comms_reduce proc~calctdfanddos->interface~comms_reduce proc~dos_get_k dos_get_k proc~calctdfanddos->proc~dos_get_k proc~dos_get_levelspacing dos_get_levelspacing proc~calctdfanddos->proc~dos_get_levelspacing proc~spin_get_nk spin_get_nk proc~tdf_kpt->proc~spin_get_nk proc~utility_w0gauss utility_w0gauss proc~tdf_kpt->proc~utility_w0gauss proc~comms_allreduce_real comms_allreduce_real interface~comms_allreduce->proc~comms_allreduce_real proc~comms_allreduce_cmplx comms_allreduce_cmplx interface~comms_allreduce->proc~comms_allreduce_cmplx proc~io_stopwatch->proc~io_error proc~wham_get_eig_deleig->get_hh_r pw90common_fourier_r_to_k pw90common_fourier_r_to_k proc~wham_get_eig_deleig->pw90common_fourier_r_to_k proc~utility_diagonalize utility_diagonalize proc~wham_get_eig_deleig->proc~utility_diagonalize proc~wham_get_deleig_a wham_get_deleig_a proc~wham_get_eig_deleig->proc~wham_get_deleig_a proc~comms_reduce_cmplx comms_reduce_cmplx interface~comms_reduce->proc~comms_reduce_cmplx proc~comms_reduce_int comms_reduce_int interface~comms_reduce->proc~comms_reduce_int proc~comms_reduce_real comms_reduce_real interface~comms_reduce->proc~comms_reduce_real proc~dos_get_k->proc~io_error proc~dos_get_k->proc~spin_get_nk proc~dos_get_k->proc~utility_w0gauss interface~pw90common_kmesh_spacing pw90common_kmesh_spacing proc~dos_get_levelspacing->interface~pw90common_kmesh_spacing ss_r ss_r proc~spin_get_nk->ss_r proc~spin_get_nk->pw90common_fourier_r_to_k proc~utility_rotate_diag utility_rotate_diag proc~spin_get_nk->proc~utility_rotate_diag proc~spin_get_nk->proc~utility_diagonalize proc~utility_w0gauss->proc~io_error proc~utility_matmul_diag utility_matmul_diag proc~utility_rotate_diag->proc~utility_matmul_diag proc~utility_diagonalize->proc~io_error zhpevx zhpevx proc~utility_diagonalize->zhpevx proc~wham_get_deleig_a->proc~utility_rotate_diag proc~wham_get_deleig_a->proc~utility_diagonalize proc~utility_rotate utility_rotate proc~wham_get_deleig_a->proc~utility_rotate proc~kmesh_spacing_mesh kmesh_spacing_mesh interface~pw90common_kmesh_spacing->proc~kmesh_spacing_mesh proc~kmesh_spacing_singleinteger kmesh_spacing_singleinteger interface~pw90common_kmesh_spacing->proc~kmesh_spacing_singleinteger
Help

Called By

proc~~calctdfanddos~~CalledByGraph proc~calctdfanddos calcTDFandDOS proc~boltzwann_main boltzwann_main proc~boltzwann_main->proc~calctdfanddos
Help

Source Code


Source Code

  subroutine calcTDFandDOS(TDF,TDFEnergyArray)
    !! 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 boltz_bandshift is set to .true., the code will also shift the
    !!       conduction bands by a given amount, as defined by the boltz_bandshift_energyshift
    !!       and boltz_bandshift_firstband input flags.
    !! 
    use w90_get_oper, only      : get_HH_R, get_SS_R, HH_R
    use w90_parameters, only    : num_wann, boltz_calc_also_dos, &
         boltz_dos_energy_step, boltz_dos_energy_min, boltz_dos_energy_max, &
         boltz_dos_adpt_smr, boltz_dos_smr_fixed_en_width, boltz_dos_adpt_smr_fac, &
         boltz_dos_adpt_smr_max, &
         param_get_smearing_type, boltz_dos_smr_index, boltz_tdf_smr_index
    use w90_utility, only       : utility_diagonalize
    use w90_wan_ham, only       : wham_get_eig_deleig

    real(kind=dp), dimension(:,:,:), 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), dimension(:), 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), dimension(3) :: kpt, orig_kpt
    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, num_s_steps, NumPtsRefined

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

    if (on_root .and. (timing_level>0)) call io_stopwatch('calcTDF',1)
    if (on_root) then
       if (boltz_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
    if(spin_decomp) then
       ndim=3
       call get_SS_R
    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 io_error('Wrong size for the TDF array in calcTDF')
    end if

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

    allocate(TDF_k(6,size(TDFEnergyArray),ndim),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating TDF_k in calcTDF')  

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

    DOS_NumPoints = int(floor((boltz_dos_energy_max - boltz_dos_energy_min)/boltz_dos_energy_step))+1
    if (DOS_NumPoints .eq. 1) DOS_NumPoints = 2
    allocate(DOS_EnergyArray(DOS_NumPoints),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating DOS_EnergyArray in calcTDF')
    do i=1,DOS_NumPoints
       DOS_EnergyArray(i) = boltz_dos_energy_min + real(i-1,dp)*boltz_dos_energy_step
    end do

    allocate(DOS_k(size(DOS_EnergyArray),ndim),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating DOS_k in calcTDF')    
    allocate(DOS_all(size(DOS_EnergyArray),ndim),stat=ierr)
    if (ierr/=0) call io_error('Error in allocating DOS_all in calcTDF')    
    dos_all = 0.0_dp



    ! I open the output files  
    if (boltz_calc_also_dos.and.on_root) then
       boltzdos_unit=io_file_unit()
       open(unit=boltzdos_unit,file=trim(seedname)//'_boltzdos.dat')  
    end if

    if (boltz_calc_also_dos.and.on_root.and.(iprint>1)) then
       write(stdout,'(5X,A)') "Smearing for DOS: "
       if (boltz_dos_adpt_smr) then
          write(stdout,'(7X,A)') trim(param_get_smearing_type(boltz_dos_smr_index)) // ", adaptive"
       else
          if (boltz_dos_smr_fixed_en_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(param_get_smearing_type(boltz_dos_smr_index)) // &
                  ", non-adaptive, width (eV) =", boltz_dos_smr_fixed_en_width
          end if
       end if
    end if
    if(boltz_calc_also_dos.and.boltz_dos_adpt_smr.and.(boltz_dos_smr_fixed_en_width.ne.0._dp).and.on_root) then
       write(stdout,'(5X,A)') "*** WARNING! boltz_dos_smr_fixed_en_width ignored since you chose"
       write(stdout,'(5X,A)') "             an adaptive smearing."
    end if

    if (on_root.and.(iprint>1)) then
       if (boltz_TDF_smr_fixed_en_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(param_get_smearing_type(boltz_TDF_smr_index)) // ", non-adaptive, width (eV) =", &
               boltz_TDF_smr_fixed_en_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: ",&
            boltz_kmesh(1),'x',boltz_kmesh(2),'x',boltz_kmesh(3)
       write(stdout,'(5X,A,I1)') "Number of electrons per state: ", num_elec_per_state
       write(stdout,'(5X,A,G18.10)') "Relaxation time (fs): ", boltz_relax_time
       if (iprint>1) then
          write(stdout,'(5X,A,G18.10)') "Energy step for TDF (eV): ", boltz_tdf_energy_step
       end if
    end if
    kweight = 1.0_dp / real(PRODUCT(boltz_kmesh),kind=dp)

    if (boltz_bandshift.and.on_root) then
       write(stdout,'(5X,A,I0,A,G18.10,A)') "Shifting energy bands with index >= ", boltz_bandshift_firstband, " by ", &
            boltz_bandshift_energyshift, " eV."
    end if

    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(boltz_kmesh)-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(boltz_kmesh)-1, so be careful when parallelizing
       loop_x= loop_tot/(boltz_kmesh(2)*boltz_kmesh(3))
       loop_y=(loop_tot-loop_x*(boltz_kmesh(2)*boltz_kmesh(3)))/boltz_kmesh(3)
       loop_z= loop_tot-loop_x*(boltz_kmesh(2)*boltz_kmesh(3)) -loop_y*boltz_kmesh(3)

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

       ! Here I get the band energies and the velocities
       call wham_get_eig_deleig(kpt, eig, del_eig, HH, delHH, UU)
       call dos_get_levelspacing(del_eig,boltz_kmesh,levelspacing_k)

       ! Here I apply a scissor operator to the conduction bands, if required in the input
       if (boltz_bandshift) then
          eig(boltz_bandshift_firstband:) = eig(boltz_bandshift_firstband:) + boltz_bandshift_energyshift
       end if

       call TDF_kpt(kpt,TDFEnergyArray,eig,del_eig,TDF_k)
       ! 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 (boltz_calc_also_dos) then
          if (boltz_dos_adpt_smr) 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(boltz_kmesh(1),dp)/ 4._dp, &
                              real(j,kind=dp)/real(boltz_kmesh(2),dp)/ 4._dp, &
                              real(k,kind=dp)/real(boltz_kmesh(3),dp)/ 4._dp /)
                         call wham_get_eig_deleig(kpt,eig,del_eig,HH,delHH,UU)
                         call dos_get_levelspacing(del_eig,boltz_kmesh,levelspacing_k)
                         call dos_get_k(kpt,DOS_EnergyArray,eig,dos_k,&
                              smr_index=boltz_dos_smr_index,&
                              adpt_smr_fac=boltz_dos_adpt_smr_fac,&
                              adpt_smr_max=boltz_dos_adpt_smr_max,&
                              levelspacing_k=levelspacing_k)
                         ! 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(kpt,DOS_EnergyArray,eig,dos_k,&
                     smr_index=boltz_dos_smr_index,&
                     adpt_smr_fac=boltz_dos_adpt_smr_fac,&
                     adpt_smr_max=boltz_dos_adpt_smr_max,&
                     levelspacing_k=levelspacing_k)
                dos_all = dos_all + dos_k * kweight                
             end if
          else
             call dos_get_k(kpt,DOS_EnergyArray,eig,dos_k,&
                  smr_index=boltz_dos_smr_index,&
                  smr_fixed_en_width=boltz_dos_smr_fixed_en_width)
             ! 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 (boltz_calc_also_dos) then
       call comms_reduce(DOS_all(1,1),size(DOS_all),'SUM')
       call comms_reduce(NumPtsRefined,1,'SUM')
       call comms_reduce(min_spacing,1,'MIN')
       call comms_reduce(max_spacing,1,'MAX')
    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')

    if (boltz_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 (boltz_dos_adpt_smr) 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: ', boltz_dos_adpt_smr_fac
          write(boltzdos_unit, '(A,I0,A,I0)') '# Number of points refined: ', NumPtsRefined, &
               ' out of ', product(boltz_kmesh)
          write(boltzdos_unit, '(A,G18.10,A,G18.10,A)') '# (Min spacing: ', min_spacing, &
               ', max spacing: ', max_spacing, ')'
       else          
          if (boltz_dos_smr_fixed_en_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 ',&
                  boltz_dos_smr_fixed_en_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. (timing_level>0)) call io_stopwatch('calcTDF',2)
    if (on_root) then
       if (boltz_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. boltz_calc_also_dos) then
       close(boltzdos_unit)
       if (iprint > 1) write(stdout,'(3X,A)') "DOS written on the " // trim(seedname)//"_boltzdos.dat file."
    end if


    deallocate(HH,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating HH in calcTDF')
    deallocate(delHH,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating delHH in calcTDF')
    deallocate(UU,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating UU in calcTDF')
    deallocate(DOS_EnergyArray,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating DOS_EnergyArray in calcTDF')
    deallocate(DOS_k,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating DOS_k in calcTDF')
    deallocate(DOS_all,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating DOS_all in calcTDF')
    deallocate(TDF_k,stat=ierr)
    if (ierr/=0) call io_error('Error in deallocating TDF_k in calcTDF')

  end subroutine calcTDFandDOS


berry_get_imf_klist berry_get_imfgh_klist berry_get_kubo_k berry_main boltzwann_main calcTDFandDOS check_and_sort_similar_centres clean_ws_translate comms_allreduce comms_allreduce_cmplx comms_allreduce_real comms_array_split comms_barrier comms_bcast comms_bcast_char comms_bcast_cmplx comms_bcast_int comms_bcast_logical comms_bcast_real comms_end comms_gatherv comms_gatherv_real comms_recv comms_recv_char comms_recv_cmplx comms_recv_int comms_recv_logical comms_recv_real comms_reduce comms_reduce_cmplx comms_reduce_int comms_reduce_real comms_scatterv comms_scatterv_int_1 comms_scatterv_int_2 comms_scatterv_int_3 comms_scatterv_real comms_send comms_send_char comms_send_cmplx comms_send_int comms_send_logical comms_send_real comms_setup conv_get_seedname conv_read_chkpt conv_read_chkpt_fmt conv_write_chkpt conv_write_chkpt_fmt dis_extract dis_extract_gamma dis_main dis_proj_froz dis_project dis_windows dos_get_k dos_get_levelspacing dos_main fourier_q_to_R gauss_freq geninterp_main get_AA_R get_BB_R get_CC_R get_FF_R get_HH_R get_module_kmesh get_smearing_index get_SS_R get_win_min group hamiltonian_dealloc hamiltonian_get_hr hamiltonian_setup hamiltonian_wigner_seitz hamiltonian_write_hr hamiltonian_write_rmn hamiltonian_write_tb internal_maxloc internal_set_kmesh internal_write_header io_date io_error io_file_unit io_get_seedname io_print_timings io_stopwatch io_time k_path k_slice kmesh_dealloc kmesh_get kmesh_get_bvectors kmesh_shell_automatic kmesh_shell_fixed kmesh_shell_from_file kmesh_spacing_mesh kmesh_spacing_singleinteger kmesh_supercell_sort kmesh_write master_sort_and_group MinusFermiDerivative my_ICOPY orthogonalize_u overlap_dealloc overlap_project overlap_project_gamma overlap_read overlap_rotate param_dealloc param_get_atoms param_get_block_length param_get_keyword param_get_keyword_block param_get_keyword_kpath param_get_keyword_vector param_get_projections param_get_range_vector param_get_smearing_type param_get_vector_length param_in_file param_lib_set_atoms param_memory_estimate param_postw90_write param_read param_read_chkpt param_uppercase param_write param_write_chkpt param_write_header plot_fermi_surface plot_interpolate_bands plot_main plot_u_matrices plot_wannier print_usage pw90common_fourier_R_to_k pw90common_fourier_R_to_k_new pw90common_fourier_R_to_k_vec pw90common_get_occ pw90common_kmesh_spacing pw90common_wanint_data_dist pw90common_wanint_get_kpoint_file pw90common_wanint_param_dist pw90common_wanint_setup qe_erf qe_erfc R_wz_sc R_wz_sc_equiv script_common script_fermi_lines sitesym_dealloc sitesym_dis_extract_symmetry sitesym_read sitesym_replace_d_matrix_band sitesym_slim_d_matrix_band sitesym_symmetrize_gradient sitesym_symmetrize_rotation sitesym_symmetrize_u_matrix sitesym_symmetrize_zmatrix sort spin_get_moment spin_get_moment_k spin_get_nk symmetrize_ukirr TDF_kpt tran_bulk tran_cut_hr_one_dim tran_dealloc tran_find_integral_signatures tran_get_ht tran_green tran_lcr tran_lcr_2c2_build_ham tran_lcr_2c2_sort tran_main tran_parity_enforce tran_read_htC tran_read_htX tran_read_htXY tran_reduce_hr tran_transfer tran_write_xyz utility_cart_to_frac utility_commutator_diag utility_compar utility_diagonalize utility_frac_to_cart utility_im_tr utility_inv2 utility_inv3 utility_lowercase utility_matmul_diag utility_metric utility_re_tr utility_recip_lattice utility_rotate utility_rotate_diag utility_string_to_coord utility_strip utility_translate_home utility_w0gauss utility_wgauss utility_zgemm wann_calc_projection wann_check_unitarity wann_domega wann_main wann_main_gamma wann_omega wann_omega_gamma wann_phases wann_spread_copy wann_svd_omega_i wann_write_r2mn wann_write_vdw_data wann_write_xyz wannier_run wannier_setup wham_get_D_h wham_get_D_h_a wham_get_deleig_a wham_get_eig_deleig wham_get_eig_UU_HH_JJlist wham_get_JJm_list wham_get_JJp_list wham_get_occ_mat_list wigner_seitz ws_translate_dist ws_write_vec