This subroutine calculates the contribution to the DOS of a single k point
\todo still to do: adapt spin_get_nk to read in input the UU rotation matrix
\note 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 if we want the final DOS to be normalized to the total number of electrons. \note The only factor that is included INSIDE this routine is the spin degeneracy factor (=num_elec_per_state variable) \note The EnergyArray is assumed to be evenly spaced (and the energy spacing is taken from EnergyArray(2)-EnergyArray(1)) \note The routine is assuming that EnergyArray has at least two elements. \note The dos_k array must have dimensions size(EnergyArray) * ndim, where ndim=1 if spin_decomp==false, or ndim=3 if spin_decomp==true. This is not checked. \note If smearing/binwidth < min_smearing_binwidth_ratio, no smearing is applied (for that k point)
\param kpt the three coordinates of the k point vector whose DOS contribution we want to calculate (in relative coordinates) \param EnergyArray array with the energy grid on which to calculate the DOS (in eV) It must have at least two elements \param eig_k array with the eigenvalues at the given k point (in eV) \param dos_k array in which the contribution is stored. Three dimensions: dos_k(energyidx, spinidx), where: - 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 \param smr_index index that tells the kind of smearing \param smr_fixed_en_width optional parameter with the fixed energy for smearing, in eV. Can be provided only if the levelspacing_k parameter is NOT given \param adpt_smr_fac optional parameter with the factor for the adaptive smearing. Can be provided only if the levelspacing_k parameter IS given \param levelspacing_k optional array with the level spacings, i.e. how much each level changes near a given point of the interpolation mesh, as given by the dos_get_levelspacing() routine If present: adaptive smearing If not present: fixed-energy-width smearing
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=dp), | intent(in), | dimension(3) | :: | kpt | ||
| real(kind=dp), | intent(in), | dimension(:) | :: | EnergyArray | ||
| real(kind=dp), | intent(in), | dimension(:) | :: | eig_k | ||
| real(kind=dp), | intent(out), | dimension(:, :) | :: | dos_k | ||
| integer, | intent(in) | :: | smr_index | |||
| real(kind=dp), | intent(in), | optional | :: | smr_fixed_en_width | ||
| real(kind=dp), | intent(in), | optional | :: | adpt_smr_fac | ||
| real(kind=dp), | intent(in), | optional | :: | adpt_smr_max | ||
| real(kind=dp), | intent(in), | optional | dimension(:) | :: | levelspacing_k | |
| complex(kind=dp), | intent(in), | optional | dimension(:, :) | :: | UU |
subroutine dos_get_k(kpt, EnergyArray, eig_k, dos_k, smr_index, &
smr_fixed_en_width, adpt_smr_fac, adpt_smr_max, levelspacing_k, UU)
use w90_io, only: io_error
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, &
num_dos_project, dos_project
use w90_spin, only: spin_get_nk
! Arguments
!
real(kind=dp), dimension(3), intent(in) :: kpt
real(kind=dp), dimension(:), intent(in) :: EnergyArray
real(kind=dp), dimension(:), intent(in) :: eig_k
real(kind=dp), dimension(:, :), intent(out) :: dos_k
integer, intent(in) :: smr_index
real(kind=dp), intent(in), optional :: smr_fixed_en_width
real(kind=dp), intent(in), optional :: adpt_smr_fac
real(kind=dp), intent(in), optional :: adpt_smr_max
real(kind=dp), dimension(:), intent(in), optional :: levelspacing_k
complex(kind=dp), dimension(:, :), intent(in), optional :: UU
! Adaptive smearing
!
real(kind=dp) :: eta_smr, arg
! Misc/Dummy
!
integer :: i, j, loop_f, min_f, max_f, num_s_steps
real(kind=dp) :: rdum, spn_nk(num_wann), alpha_sq, beta_sq
real(kind=dp) :: binwidth, r_num_elec_per_state
logical :: DoSmearing
if (present(levelspacing_k)) then
if (present(smr_fixed_en_width)) &
call io_error('Cannot call doskpt with levelspacing_k and ' &
//'with smr_fixed_en_width parameters together')
if (.not. (present(adpt_smr_fac))) &
call io_error('Cannot call doskpt with levelspacing_k and ' &
//'without adpt_smr_fac parameter')
if (.not. (present(adpt_smr_max))) &
call io_error('Cannot call doskpt with levelspacing_k and ' &
//'without adpt_smr_max parameter')
else
if (present(adpt_smr_fac)) &
call io_error('Cannot call doskpt without levelspacing_k and ' &
//'with adpt_smr_fac parameter')
if (present(adpt_smr_max)) &
call io_error('Cannot call doskpt without levelspacing_k and ' &
//'with adpt_smr_max parameter')
if (.not. (present(smr_fixed_en_width))) &
call io_error('Cannot call doskpt without levelspacing_k and ' &
//'without smr_fixed_en_width parameter')
end if
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)
dos_k = 0.0_dp
do i = 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(i))/2.0_dp ! |alpha|^2
! Contribution to spin-down DOS
beta_sq = 1.0_dp - alpha_sq ! |beta|^2 = 1 - |alpha|^2
end if
if (.not. present(levelspacing_k)) then
eta_smr = smr_fixed_en_width
else
! Eq.(35) YWVS07
eta_smr = min(levelspacing_k(i)*adpt_smr_fac, adpt_smr_max)
! eta_smr=max(eta_smr,min_smearing_binwidth_ratio) !! No: it would render the next if always false
end if
! Faster optimization: I precalculate the indices
if (eta_smr/binwidth < min_smearing_binwidth_ratio) then
min_f = max(nint((eig_k(i) - EnergyArray(1))/ &
(EnergyArray(size(EnergyArray)) - EnergyArray(1)) &
*real(size(EnergyArray) - 1, kind=dp)) + 1, 1)
max_f = min(nint((eig_k(i) - 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(i) - smearing_cutoff*eta_smr - EnergyArray(1))/ &
(EnergyArray(size(EnergyArray)) - EnergyArray(1)) &
*real(size(EnergyArray) - 1, kind=dp)) + 1, 1)
max_f = min(nint((eig_k(i) + smearing_cutoff*eta_smr - 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
! kind of smearing read from input (internal smearing_index variable)
if (DoSmearing) then
arg = (EnergyArray(loop_f) - eig_k(i))/eta_smr
rdum = utility_w0gauss(arg, smr_index)/eta_smr
else
rdum = 1._dp/(EnergyArray(2) - EnergyArray(1))
end if
!
! Contribution to total DOS
!
if (num_dos_project == num_wann) then
!
! Total DOS (default): do not loop over j, to save time
!
dos_k(loop_f, 1) = dos_k(loop_f, 1) + rdum*r_num_elec_per_state
! [GP] 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
dos_k(loop_f, 2) = dos_k(loop_f, 2) + rdum*alpha_sq
! Spin-down contribution
dos_k(loop_f, 3) = dos_k(loop_f, 3) + rdum*beta_sq
end if
else ! 0<num_dos_project<num_wann
!
! Partial DOS, projected onto the WFs with indices
! n=dos_project(1:num_dos_project)
!
do j = 1, num_dos_project
dos_k(loop_f, 1) = dos_k(loop_f, 1) + rdum*r_num_elec_per_state &
*abs(UU(dos_project(j), i))**2
if (spin_decomp) then
! Spin-up contribution
dos_k(loop_f, 2) = dos_k(loop_f, 2) &
+ rdum*alpha_sq*abs(UU(dos_project(j), i))**2
! Spin-down contribution
dos_k(loop_f, 3) = dos_k(loop_f, 3) &
+ rdum*beta_sq*abs(UU(dos_project(j), i))**2
end if
enddo
endif
enddo !loop_f
end do !loop over bands
end subroutine dos_get_k