Skip to content

Commit

Permalink
Moved screening related stuff to a new module; implemented the imagin…
Browse files Browse the repository at this point in the history 
…ary part of the bare polarizability.
  • Loading branch information
@nakib
nakib committed May 7, 2024
1 parent 9941fc4 commit af4043d
Show file tree
Hide file tree
Showing 3 changed files with 347 additions and 333 deletions.
3 changes: 1 addition & 2 deletions app/elphbolt.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,8 +31,7 @@ program elphbolt
use phonon_module, only: phonon
use wannier_module, only: wannier
use bte_module, only: bte
use bz_sums, only: calculate_dos, calculate_qTF, calculate_el_dos_fermi, calculate_el_Ws, &
calculate_RPA_dielectric_3d_G0_qpath
use bz_sums, only: calculate_dos, calculate_qTF, calculate_el_dos_fermi, calculate_el_Ws
use interactions, only: calculate_gReq, calculate_gkRp, calculate_3ph_interaction, &
calculate_eph_interaction_ibzq, calculate_eph_interaction_ibzk, &
calculate_echimp_interaction_ibzk, calculate_coarse_grained_3ph_vertex, &
Expand Down
334 changes: 3 additions & 331 deletions src/bz_sums.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ module bz_sums
use misc, only: exit_with_message, print_message, write2file_rank2_real, &
distribute_points, Bose, Fermi, binsearch, mux_vector, twonorm, linspace, &
write2file_rank1_real
use screening_module, only: calculate_qTF
use phonon_module, only: phonon
use electron_module, only: electron
use crystal_module, only: crystal
Expand All @@ -32,345 +33,16 @@ module bz_sums

implicit none

public calculate_dos, calculate_transport_coeff, calculate_qTF, &
public calculate_dos, calculate_transport_coeff, &
calculate_el_dos_fermi, calculate_el_Ws, calculate_cumulative_transport_coeff, &
calculate_el_dos_fermi_gaussian, calculate_el_Ws_gaussian, &
calculate_RPA_dielectric_3d_G0_qpath
calculate_el_dos_fermi_gaussian, calculate_el_Ws_gaussian
private calculate_el_dos, calculate_ph_dos_iso

interface calculate_dos
module procedure :: calculate_el_dos, calculate_ph_dos_iso
end interface calculate_dos

contains

complex(r64) function RPA_polarizability_3d(Omega, q_indvec, el, pcell_vol, T)
!! Polarizability of the 3d Kohn-Sham system in the
!! random-phase approximation (RPA) using Eq. B1 and B2
!! of Knapen, Kozaczuk and Lin Phys. Rev. D 104, 015031 (2021).
!!
!! Here we calculate the diagonal in G-G' space. Moreover,
!! we use the approximation G.r -> 0.
!!
!! Omega Energy of excitation in the electron gas
!! q_indvec Wave vector of excitation in the electron gas (0-based integer triplet)
!! el Electron data type
!! pcell_vol Primitive unit cell volume
!! T Temperature (K)

real(r64), intent(in) :: Omega, pcell_vol, T
integer(i64), intent(in) :: q_indvec(3)
type(electron), intent(in) :: el

!Locals
integer(i64) :: m, n, ik, ikp, ikp_window, k_indvec(3), kp_indvec(3)
real(r64) :: overlap, ek, ekp, realpart, imagpart
complex(r64) :: aux

aux = 0.0
do m = 1, el%numbands
do ik = 1, el%nwv
ek = el%ens(ik, m)

!Apply energy window to initial electron
if(abs(ek - el%enref) > el%fsthick) cycle

!Calculate index vector of final electron: k' = k + q
k_indvec = nint(el%wavevecs(ik, :)*el%wvmesh)
kp_indvec = modulo(k_indvec + q_indvec, el%wvmesh) !0-based index vector

!Multiplex kp_indvec
ikp = mux_vector(kp_indvec, el%wvmesh, 0_i64)

!Check if final electron wave vector is within FBZ blocks
call binsearch(el%indexlist, ikp, ikp_window)
if(ikp_window < 0) cycle

do n = 1, el%numbands
ekp = el%ens(ikp_window, n)

!Apply energy window to final electron
if(abs(ekp - el%enref) > el%fsthick) cycle

!This is |U(k')U^\dagger(k)|_nm squared
!(Recall that U^\dagger(k) is the diagonalizer of the electronic hamiltonian.)
overlap = (abs(dot_product(el%evecs(ikp_window, n, :), el%evecs(ik, m, :))))**2

aux = aux + &
(Fermi(ekp, el%chempot, T) - &
Fermi(ek, el%chempot, T))* &
resolvent(el, m, ik, ekp - Omega)*overlap
end do
end do
end do
!Recall that the resolvent is already normalized in the full wave vector mesh.
!As such, the 1/product(el%wvmesh) is not needed in the expression below.
aux = aux*el%spindeg/pcell_vol

!RPA_polarizability_3d = aux

!Zero out extremely small numbers
realpart = real(aux)
imagpart = imag(aux)
if(abs(real(aux)) < 1.0e-30) realpart = 0.0_r64
if(abs(imag(aux)) < 1.0e-30) imagpart = 0.0_r64

RPA_polarizability_3d = realpart + imagpart*oneI
end function RPA_polarizability_3d

subroutine calculate_RPA_dielectric_3d(el, crys, num)
!! Dielectric function of the 3d Kohn-Sham system in the
!! random-phase approximation (RPA) using Eq. B1 and B2
!! of Knapen, Kozaczuk and Lin Phys. Rev. D 104, 015031 (2021).
!!
!! Here we calculate the diagonal in G-G' space. Moreover,
!! we use the approximation G.r -> 0.
!!
!! el Electron data type
!! crys Crystal data type
!! num Numerics data type

type(electron), intent(in) :: el
type(crystal), intent(in) :: crys
type(numerics), intent(in) :: num

!Locals
real(r64) :: Omega, k(3), kp(3), q(3), ek, ekp, Gplusq_crys(3)
integer(i64) :: ik, ikp, m, n, &
k_indvec(3), kp_indvec(3), q_indvec(3), count, &
start, end, chunk, num_active_images
integer :: ig1, ig2, ig3, igmax, num_gmax
complex(r64) :: polarizability
complex(r64), allocatable :: diel_ik(:)
character(len = 1024) :: filename

print*, el%numbands

!This sets how large the G vectors can be. Choosing 3
!means including -3G to +3G in the calculations. This
!should be a safe range.
igmax = 3
num_gmax = (2*igmax + 1)**3

!Allocate diel_ik to hold maximum possible Omega x G points
allocate(diel_ik(num_gmax*el%nstates_inwindow*el%numbands))

!Distribute points among images
call distribute_points(el%nwv_irred, chunk, start, end, num_active_images)

if(this_image() == 1) then
write(*, "(A, I10)") " #k-vecs = ", el%nwv_irred
write(*, "(A, I10)") " #k-vecs/image <= ", chunk
end if

do ik = start, end !Over IBZ k points
!Initiate counter for Omega x G points
count = 0

!Initial (IBZ blocks) wave vector (crystal coords.)
k = el%wavevecs_irred(ik, :)

!Convert from crystal to 0-based index vector
k_indvec = nint(k*el%wvmesh)

do m = 1, el%numbands
!IBZ electron energy
ek = el%ens_irred(ik, m)

!Check energy window
if(abs(ek - el%enref) > el%fsthick) cycle

do ikp = 1, el%nwv
!Final wave vector (crystal coords.)
kp = el%wavevecs(ikp, :)

!Convert from crystal to 0-based index vector
kp_indvec = nint(kp*el%wvmesh)

!Calculate q_indvec (folded back to the 1BZ)
q_indvec = modulo(kp_indvec - k_indvec, el%wvmesh) !0-based index vector

!Calculate q_indvec (without folding back to the 1BZ)
!q_indvec = kp_indvec - k_indvec !0-based index vector

do n = 1, el%numbands
ekp = el%ens(ikp, n)

!Check energy window
if(abs(el%ens(ikp, n) - el%enref) > el%fsthick) cycle

!Electron-hole pair/plasmon energy
Omega = ekp - ek

!Calculate RPA polarizability
polarizability = &
RPA_polarizability_3d(Omega, q_indvec, el, crys%volume, crys%T)

do ig1 = -igmax, igmax
do ig2 = -igmax, igmax
do ig3 = -igmax, igmax
!Counter for Omega x G points
count = count + 1

!Calculate G+q
Gplusq_crys = ([ig1, ig2, ig3] + q_indvec)/dble(el%wvmesh)

!Calculate RPA dielectric (diagonal in G-G' space)
diel_ik(count) = 1.0_r64 - &
(1.0_r64/twonorm(matmul(crys%reclattvecs, Gplusq_crys)))**2* &
polarizability/perm0*qe*1.0e9_r64
end do
end do
end do

end do
end do
end do

!Change to data output directory
call chdir(trim(adjustl(num%epsilondir)))

!Write data in binary format
!Note: this will overwrite existing data!
write (filename, '(I9)') ik
filename = 'epsilon_RPA.ik'//trim(adjustl(filename))
open(1, file = trim(filename), status = 'replace', access = 'stream')
write(1) count
write(1) diel_ik(1:count)
close(1)
end do

sync all
end subroutine calculate_RPA_dielectric_3d

!DEBUG/TEST
subroutine calculate_RPA_dielectric_3d_G0_qpath(el, crys, num)
!! Dielectric function of the 3d Kohn-Sham system in the
!! random-phase approximation (RPA) using Eq. B1 and B2
!! of Knapen, Kozaczuk and Lin Phys. Rev. D 104, 015031 (2021).
!!
!! Here we calculate the diagonal in G-G' space. Moreover,
!! we use the approximation G.r -> 0.
!!
!! el Electron data type
!! crys Crystal data type
!! num Numerics data type

type(electron), intent(in) :: el
type(crystal), intent(in) :: crys
type(numerics), intent(in) :: num

!Locals
real(r64), allocatable :: energylist(:), qlist(:, :), qmaglist(:)
real(r64) :: qcrys(3)
integer(i64) :: iq, iOmega, numomega, numq, &
start, end, chunk, num_active_images, qxmesh
complex(r64) :: polarizability
complex(r64), allocatable :: diel(:, :)
character(len = 1024) :: filename
real(r64) :: a0, eps0_q0_prefac, omega_plasma

!!$ a0 = 1.0e-9_r64*bohr2nm !Bohr radius in m
!!$ eps0_q0_prefac = (4.0_r64/9.0_r64/pi)/a0* &
!!$ (3.0_r64*abs(sum(el%conc))*1.0e6_r64)**(-1.0_r64/3.0_r64) !
omega_plasma = 0.5025125628E-01 !eV

!TEST
!Material: Si
!Uniform energy mesh from 0-1.0 eV with uniform q-vecs in Gamma-Gamma along x
numomega = 200
numq = 50
qxmesh = 200
!Create qlist in crystal coordinates
allocate(qlist(numq, 3), qmaglist(numq))
do iq = 1, numq
qlist(iq, :) = [(iq - 1.0_r64)/qxmesh, 0.0_r64, 0.0_r64]
qmaglist(iq) = twonorm(matmul(crys%reclattvecs, qlist(iq, :)))
end do

!Create energy grid
allocate(energylist(numomega))
call linspace(energylist, 0.0_r64, 0.5_r64, numomega)

!Allocate diel_ik to hold maximum possible Omega
allocate(diel(numq, numomega))
diel = 0.0_r64

!Distribute points among images
call distribute_points(numq, chunk, start, end, num_active_images)

if(this_image() == 1) then
write(*, "(A, I10)") " #q-vecs = ", numq
write(*, "(A, I10)") " #q-vecs/image <= ", chunk
end if

do iq = start, end !Over IBZ k points
qcrys = qlist(iq, :) !crystal coordinates

if(all(qcrys == 0.0_r64)) then
diel(iq, :) = 1.0_r64 - (omega_plasma/energylist)**2 + 1.0e-3*oneI
else
do iOmega = 1, numomega
!Calculate RPA polarizability
polarizability = &
RPA_polarizability_3d(energylist(iOmega), &
nint(qcrys*numq)+0_i64, el, crys%volume, crys%T)

!Calculate RPA dielectric (diagonal in G-G' space)
diel(iq, iOmega) = 1.0_r64 - &
(1.0_r64/twonorm(matmul(crys%reclattvecs, qcrys)))**2* &
polarizability/perm0*qe*1.0e9_r64
end do
end if
end do

sync all
call co_sum(diel)
sync all

!Print to file
call write2file_rank2_real("RPA_dielectric_3D_G0_qpath", qlist)
call write2file_rank1_real("RPA_dielectric_3D_G0_qmagpath", qmaglist)
call write2file_rank1_real("RPA_dielectric_3D_G0_Omega", energylist)
call write2file_rank2_real("RPA_dielectric_3D_G0_real", real(diel)*1.0_r64)
call write2file_rank2_real("RPA_dielectric_3D_G0_imag", imag(diel)*1.0_r64)
end subroutine calculate_RPA_dielectric_3d_G0_qpath
!!

subroutine calculate_qTF(crys, el)
!! Calculate Thomas-Fermi screening wave vector from the static
!! limit of the Lindhard function.

type(crystal), intent(inout) :: crys
type(electron), intent(in) :: el

!Local variables
real(r64) :: beta, fFD
integer(i64) :: ib, ik

beta = 1.0_r64/kB/crys%T/qe !1/J
crys%qTF=0.d0

if(crys%epsilon0 /= 0) then
call print_message("Calculating Thomas-Fermi screening wave vector...")

do ib = 1, el%numbands
do ik = 1, el%nwv
fFD = Fermi(el%ens(ik, ib), el%chempot, crys%T)
crys%qTF = crys%qTF + fFD*(1.0_r64 - fFD)
end do
end do

!Free-electron gas Thomas-Fermi model
! qTF**2 = spindeg*e^2*beta/nptq/vol_pcell/perm0/epsilon0*Sum_{BZ}f0_{k}(1-f0_{k})
crys%qTF = sqrt(1.0e9_r64*crys%qTF*el%spindeg*beta*qe**2/product(el%wvmesh)&
/crys%volume/perm0/crys%epsilon0) !nm^-1

if(this_image() == 1) then
write(*, "(A, 1E16.8, A)") ' Thomas-Fermi screening wave vector = ', crys%qTF, ' 1/nm'
end if
end if
end subroutine calculate_qTF

subroutine calculate_el_dos_Fermi_Gaussian(el, reclattvecs)
!! Calculate spin-normalized electron density of states at the Fermi level
Expand Down
Loading

0 comments on commit af4043d

Please sign in to comment.