This subroutine calculates the contribution to the TDF of a single k point
This routine does not use the adaptive smearing; in fact, for non-zero temperatures one often doesn't even need to smear. It simply uses a standard smearing as defined by the variables boltz_TDF_smr_fixed_en_width and boltz_TDF_smr_index
still to do: adapt spin_get_nk to read in input the UU rotation matrix
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/cell_volume if we want to calculate 1/(2*pi)^3 * \int_{BZ} d^3 k The only factor that is included INSIDE this routine is the spin degeneracy factor (=2 if spinors is .false., =1 if spinors is .true.) The EnergyArray is assumed to be evenly spaced (and the energy spacing is taken from EnergyArray(2)-EnergyArray(1)) The routine is assuming that EnergyArray has at least two elements. The meaning of the three indices of the TDF_k array is different with respect to those of the dos_k array returned by the dos_get_k routine The TDF_k array must have dimensions 6 * size(EnergyArray) * ndim, where ndim=1 if spin_decomp==false, or ndim=3 if spin_decomp==true. This is not checked.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=dp), | intent(in), | dimension(3) | :: | kpt | the three coordinates of the k point vector whose DOS contribution we want to calculate (in relative coordinates) |
|
| real(kind=dp), | intent(in), | dimension(:) | :: | EnergyArray | array with the energy grid on which to calculate the DOS (in eV) It must have at least two elements |
|
| real(kind=dp), | intent(in), | dimension(:) | :: | eig_k | array with the eigenvalues at the given k point (in eV) |
|
| real(kind=dp), | intent(in), | dimension(:, :) | :: | deleig_k | array with the band derivatives at the given k point (in eV * angstrom / (2pi) as internally given by the code) already corrected in case of degeneracies, as returned by the wham_get_deleig_a routine |
|
| real(kind=dp), | intent(out), | dimension(:, :, :) | :: | TDF_k | TDF_k array in which the contribution is stored. Three dimensions: TDF_k(ij, energyidx, spinidx), where: - ij indexes the components of the TDF (symmetric) tensor (1=XX, 2=XY, ...); see the global constants defined in the module - 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 |
subroutine TDF_kpt(kpt, EnergyArray, eig_k, deleig_k, TDF_k)
!! This subroutine calculates the contribution to the TDF of a single k point
!!
!! This routine does not use the adaptive smearing; in fact, for non-zero temperatures
!! one often doesn't even need to smear. It simply uses a standard smearing as defined by
!! the variables boltz_TDF_smr_fixed_en_width and boltz_TDF_smr_index
!!
!! still to do: adapt spin_get_nk to read in input the UU rotation matrix
!!
!! 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/cell_volume
!! if we want to calculate 1/(2*pi)^3 * \int_{BZ} d^3 k
!! The only factor that is included INSIDE this routine is the spin degeneracy
!! factor (=2 if spinors is .false., =1 if spinors is .true.)
!! The EnergyArray is assumed to be evenly spaced (and the energy spacing
!! is taken from EnergyArray(2)-EnergyArray(1))
!! The routine is assuming that EnergyArray has at least two elements.
!! The meaning of the three indices of the TDF_k array is different with respect to
!! those of the dos_k array returned by the dos_get_k routine
!! The TDF_k array must have dimensions 6 * size(EnergyArray) * ndim, where
!! ndim=1 if spin_decomp==false, or ndim=3 if spin_decomp==true. This is not checked.
!!
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, &
boltz_TDF_smr_fixed_en_width, boltz_TDF_smr_index, boltz_relax_time
use w90_spin, only: spin_get_nk
! Arguments
!
real(kind=dp), dimension(3), intent(in) :: kpt
!! the three coordinates of the k point vector whose DOS contribution we
!! want to calculate (in relative coordinates)
real(kind=dp), dimension(:), intent(in) :: EnergyArray
!! array with the energy grid on which to calculate the DOS (in eV)
!! It must have at least two elements
real(kind=dp), dimension(:), intent(in) :: eig_k
!! array with the eigenvalues at the given k point (in eV)
real(kind=dp), dimension(:, :), intent(in) :: deleig_k
!! array with the band derivatives at the given k point
!! (in eV * angstrom / (2pi) as internally given by the code)
!! already corrected in case of degeneracies, as returned by the
!! wham_get_deleig_a routine
real(kind=dp), dimension(:, :, :), intent(out) :: TDF_k
!! TDF_k array in which the contribution is stored. Three dimensions:
!! TDF_k(ij, energyidx, spinidx), where:
!! - ij indexes the components of the TDF (symmetric) tensor (1=XX, 2=XY, ...);
!! see the global constants defined in the module
!! - 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
! Adaptive smearing
!
real(kind=dp) :: smear, arg
! Misc/Dummy
!
integer :: BandIdx, loop_f, min_f, max_f
real(kind=dp) :: rdum, spn_nk(num_wann), alpha_sq, beta_sq
real(kind=dp) :: binwidth, r_num_elec_per_state
logical :: DoSmearing
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)
TDF_k = 0.0_dp
do BandIdx = 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(BandIdx))/2.0_dp ! |alpha|^2
! Contribution to spin-down DOS
beta_sq = 1.0_dp - alpha_sq ! |beta|^2 = 1 - |alpha|^2
end if
! Do not use an adaptive smearing here, it would require the knowledge of second derivatives
! And typically, when working at not too small temperatures, smearing is not needed
! Faster optimization: I precalculate the indices
! Value of the smearing in eV; default = 0 eV, i.e. no smearing
smear = boltz_TDF_smr_fixed_en_width
if (smear/binwidth < min_smearing_binwidth_ratio) then
min_f = max(nint((eig_k(BandIdx) - EnergyArray(1))/ &
(EnergyArray(size(EnergyArray)) - EnergyArray(1)) &
*real(size(EnergyArray) - 1, kind=dp)) + 1, 1)
max_f = min(nint((eig_k(BandIdx) - 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(BandIdx) - smearing_cutoff*smear - EnergyArray(1))/ &
(EnergyArray(size(EnergyArray)) - EnergyArray(1)) &
*real(size(EnergyArray) - 1, kind=dp)) + 1, 1)
max_f = min(nint((eig_k(BandIdx) + smearing_cutoff*smear - 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
if (DoSmearing) then
arg = (EnergyArray(loop_f) - eig_k(BandIdx))/smear
rdum = utility_w0gauss(arg, boltz_TDF_smr_index)/smear
else
rdum = 1._dp/(EnergyArray(2) - EnergyArray(1))
end if
!
! Contribution to total DOS
!
TDF_k(XX, loop_f, 1) = TDF_k(XX, loop_f, 1) + rdum* &
r_num_elec_per_state*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 1)
TDF_k(XY, loop_f, 1) = TDF_k(XY, loop_f, 1) + rdum* &
r_num_elec_per_state*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 2)
TDF_k(YY, loop_f, 1) = TDF_k(YY, loop_f, 1) + rdum* &
r_num_elec_per_state*deleig_k(BandIdx, 2)*deleig_k(BandIdx, 2)
TDF_k(XZ, loop_f, 1) = TDF_k(XZ, loop_f, 1) + rdum* &
r_num_elec_per_state*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 3)
TDF_k(YZ, loop_f, 1) = TDF_k(YZ, loop_f, 1) + rdum* &
r_num_elec_per_state*deleig_k(BandIdx, 2)*deleig_k(BandIdx, 3)
TDF_k(ZZ, loop_f, 1) = TDF_k(ZZ, loop_f, 1) + rdum* &
r_num_elec_per_state*deleig_k(BandIdx, 3)*deleig_k(BandIdx, 3)
! 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
TDF_k(XX, loop_f, 2) = TDF_k(XX, loop_f, 2) + rdum* &
alpha_sq*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 1)
TDF_k(XY, loop_f, 2) = TDF_k(XY, loop_f, 2) + rdum* &
alpha_sq*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 2)
TDF_k(YY, loop_f, 2) = TDF_k(YY, loop_f, 2) + rdum* &
alpha_sq*deleig_k(BandIdx, 2)*deleig_k(BandIdx, 2)
TDF_k(XZ, loop_f, 2) = TDF_k(XZ, loop_f, 2) + rdum* &
alpha_sq*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 3)
TDF_k(YZ, loop_f, 2) = TDF_k(YZ, loop_f, 2) + rdum* &
alpha_sq*deleig_k(BandIdx, 2)*deleig_k(BandIdx, 3)
TDF_k(ZZ, loop_f, 2) = TDF_k(ZZ, loop_f, 2) + rdum* &
alpha_sq*deleig_k(BandIdx, 3)*deleig_k(BandIdx, 3)
! Spin-down contribution
TDF_k(XX, loop_f, 3) = TDF_k(XX, loop_f, 3) + rdum* &
beta_sq*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 1)
TDF_k(XY, loop_f, 3) = TDF_k(XY, loop_f, 3) + rdum* &
beta_sq*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 2)
TDF_k(YY, loop_f, 3) = TDF_k(YY, loop_f, 3) + rdum* &
beta_sq*deleig_k(BandIdx, 2)*deleig_k(BandIdx, 2)
TDF_k(XZ, loop_f, 3) = TDF_k(XZ, loop_f, 3) + rdum* &
beta_sq*deleig_k(BandIdx, 1)*deleig_k(BandIdx, 3)
TDF_k(YZ, loop_f, 3) = TDF_k(YZ, loop_f, 3) + rdum* &
beta_sq*deleig_k(BandIdx, 2)*deleig_k(BandIdx, 3)
TDF_k(ZZ, loop_f, 3) = TDF_k(ZZ, loop_f, 3) + rdum* &
beta_sq*deleig_k(BandIdx, 3)*deleig_k(BandIdx, 3)
end if
end do
end do !loop over bands
! I multiply it here, since I am assuming a constant relaxation time, independent of the band index
! (actually, it is also independent of k)
TDF_k = TDF_k*boltz_relax_time
end subroutine TDF_kpt