This routine calculates the Transport Distribution Function (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.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| 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 , 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 |
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