Computes the electronic density of states. Can resolve into up-spin and down-spin parts, project onto selected Wannier orbitals, and use adaptive broadening, as in PRB 75, 195121 (2007) [YWVS07].
subroutine dos_main
!=======================================================!
! !
!! Computes the electronic density of states. Can
!! resolve into up-spin and down-spin parts, project
!! onto selected Wannier orbitals, and use adaptive
!! broadening, as in PRB 75, 195121 (2007) [YWVS07].
! !
!=======================================================!
use w90_io, only: io_error, io_file_unit, io_date, io_stopwatch, &
seedname, stdout
use w90_comms, only: on_root, num_nodes, my_node_id, comms_reduce
use w90_postw90_common, only: num_int_kpts_on_node, int_kpts, weight, &
pw90common_fourier_R_to_k
use w90_parameters, only: num_wann, dos_energy_min, dos_energy_max, &
dos_energy_step, timing_level, &
wanint_kpoint_file, dos_kmesh, &
dos_smr_index, dos_adpt_smr, dos_adpt_smr_fac, &
dos_adpt_smr_max, spin_decomp, &
dos_smr_fixed_en_width, &
dos_project, num_dos_project
use w90_get_oper, only: get_HH_R, get_SS_R, HH_R
use w90_wan_ham, only: wham_get_eig_deleig
use w90_utility, only: utility_diagonalize
! 'dos_k' contains contrib. from one k-point,
! 'dos_all' from all nodes/k-points (first summed on one node and
! then reduced (i.e. summed) over all nodes)
!
real(kind=dp), allocatable :: dos_k(:, :)
real(kind=dp), allocatable :: dos_all(:, :)
real(kind=dp) :: kweight, kpt(3), omega
integer :: i, loop_x, loop_y, loop_z, loop_tot, ifreq
integer :: dos_unit, ndim, ierr
real(kind=dp), dimension(:), allocatable :: dos_energyarray
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)
num_freq = nint((dos_energy_max - dos_energy_min)/dos_energy_step) + 1
if (num_freq == 1) num_freq = 2
d_omega = (dos_energy_max - dos_energy_min)/(num_freq - 1)
allocate (dos_energyarray(num_freq), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating dos_energyarray in ' &
//'dos subroutine')
do ifreq = 1, num_freq
dos_energyarray(ifreq) = dos_energy_min + real(ifreq - 1, dp)*d_omega
end do
allocate (HH(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating HH in dos')
allocate (delHH(num_wann, num_wann, 3), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating delHH in dos')
allocate (UU(num_wann, num_wann), stat=ierr)
if (ierr /= 0) call io_error('Error in allocating UU in dos')
call get_HH_R
if (spin_decomp) then
ndim = 3
call get_SS_R
else
ndim = 1
end if
allocate (dos_k(num_freq, ndim))
allocate (dos_all(num_freq, ndim))
if (on_root) then
if (timing_level > 1) call io_stopwatch('dos', 1)
! write(stdout,'(/,1x,a)') '============'
! write(stdout,'(1x,a)') 'Calculating:'
! write(stdout,'(1x,a)') '============'
write (stdout, '(/,/,1x,a)') &
'Properties calculated in module d o s'
write (stdout, '(1x,a)') &
'--------------------------------------'
if (num_dos_project == num_wann) then
write (stdout, '(/,3x,a)') '* Total density of states (_dos)'
else
write (stdout, '(/,3x,a)') &
'* Density of states projected onto selected WFs (_dos)'
write (stdout, '(3x,a)') 'Selected WFs |Rn> are:'
do i = 1, num_dos_project
write (stdout, '(5x,a,2x,i3)') 'n =', dos_project(i)
enddo
endif
write (stdout, '(/,5x,a,f9.4,a,f9.4,a)') &
'Energy range: [', dos_energy_min, ',', dos_energy_max, '] eV'
write (stdout, '(/,5x,a,(f6.3,1x))') &
'Adaptive smearing width prefactor: ', &
dos_adpt_smr_fac
write (stdout, '(/,/,1x,a20,3(i0,1x))') 'Interpolation grid: ', &
dos_kmesh(1:3)
end if
dos_all = 0.0_dp
if (wanint_kpoint_file) then
!
! Unlike for optical properties, this should always work for the DOS
!
if (on_root) write (stdout, '(/,1x,a)') 'Sampling the irreducible BZ only'
! Loop over k-points on the irreducible wedge of the Brillouin zone,
! read from file 'kpoint.dat'
!
do loop_tot = 1, num_int_kpts_on_node(my_node_id)
kpt(:) = int_kpts(:, loop_tot)
if (dos_adpt_smr) then
call wham_get_eig_deleig(kpt, eig, del_eig, HH, delHH, UU)
call dos_get_levelspacing(del_eig, dos_kmesh, levelspacing_k)
call dos_get_k(kpt, dos_energyarray, eig, dos_k, &
smr_index=dos_smr_index, &
adpt_smr_fac=dos_adpt_smr_fac, &
adpt_smr_max=dos_adpt_smr_max, &
levelspacing_k=levelspacing_k, &
UU=UU)
else
call pw90common_fourier_R_to_k(kpt, HH_R, HH, 0)
call utility_diagonalize(HH, num_wann, eig, UU)
call dos_get_k(kpt, dos_energyarray, eig, dos_k, &
smr_index=dos_smr_index, &
smr_fixed_en_width=dos_smr_fixed_en_width, &
UU=UU)
end if
dos_all = dos_all + dos_k*weight(loop_tot)
end do
else
if (on_root) write (stdout, '(/,1x,a)') 'Sampling the full BZ'
kweight = 1.0_dp/real(PRODUCT(dos_kmesh), kind=dp)
do loop_tot = my_node_id, PRODUCT(dos_kmesh) - 1, num_nodes
loop_x = loop_tot/(dos_kmesh(2)*dos_kmesh(3))
loop_y = (loop_tot - loop_x*(dos_kmesh(2) &
*dos_kmesh(3)))/dos_kmesh(3)
loop_z = loop_tot - loop_x*(dos_kmesh(2)*dos_kmesh(3)) &
- loop_y*dos_kmesh(3)
kpt(1) = real(loop_x, dp)/real(dos_kmesh(1), dp)
kpt(2) = real(loop_y, dp)/real(dos_kmesh(2), dp)
kpt(3) = real(loop_z, dp)/real(dos_kmesh(3), dp)
if (dos_adpt_smr) then
call wham_get_eig_deleig(kpt, eig, del_eig, HH, delHH, UU)
call dos_get_levelspacing(del_eig, dos_kmesh, levelspacing_k)
call dos_get_k(kpt, dos_energyarray, eig, dos_k, &
smr_index=dos_smr_index, &
adpt_smr_fac=dos_adpt_smr_fac, &
adpt_smr_max=dos_adpt_smr_max, &
levelspacing_k=levelspacing_k, &
UU=UU)
else
call pw90common_fourier_R_to_k(kpt, HH_R, HH, 0)
call utility_diagonalize(HH, num_wann, eig, UU)
call dos_get_k(kpt, dos_energyarray, eig, dos_k, &
smr_index=dos_smr_index, &
smr_fixed_en_width=dos_smr_fixed_en_width, &
UU=UU)
end if
dos_all = dos_all + dos_k*kweight
end do
end if
! Collect contributions from all nodes
!
call comms_reduce(dos_all(1, 1), num_freq*ndim, 'SUM')
if (on_root) then
write (stdout, '(1x,a)') 'Output data files:'
write (stdout, '(/,3x,a)') trim(seedname)//'-dos.dat'
dos_unit = io_file_unit()
open (dos_unit, FILE=trim(seedname)//'-dos.dat', STATUS='UNKNOWN', &
FORM='FORMATTED')
do ifreq = 1, num_freq
omega = dos_energyarray(ifreq)
write (dos_unit, '(4E16.8)') omega, dos_all(ifreq, :)
enddo
close (dos_unit)
if (timing_level > 1) call io_stopwatch('dos', 2)
end if
deallocate (HH, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating HH in dos_main')
deallocate (delHH, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating delHH in dos_main')
deallocate (UU, stat=ierr)
if (ierr /= 0) call io_error('Error in deallocating UU in dos_main')
end subroutine dos_main