Skip to content

Commit

Permalink
First pass through OTF Xee.
Browse files Browse the repository at this point in the history
  • Loading branch information
@nakib
nakib committed Aug 30, 2024
1 parent db12875 commit 6de799d
Show file tree
Hide file tree
Showing 3 changed files with 247 additions and 40 deletions.
174 changes: 170 additions & 4 deletions src/interactions.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ module interactions
#endif

use precision, only: i64, r64
use params, only: pi, twopi, amu, qe, hbar_eVps, perm0, oneI
use params, only: pi, twopi, amu, qe, hbar_eVps, perm0, oneI, kB
use misc, only: exit_with_message, print_message, distribute_points, &
demux_state, mux_vector, mux_state, expi, Bose, binsearch, Fermi, &
twonorm, write2file_rank2_real, demux_vector, interpolate, expm1, &
Expand Down Expand Up @@ -114,6 +114,31 @@ pure real(r64) function gchimp2(el, crys, qcrys, evec_k, evec_kp)

gchimp2 = prefac*overlap*Gsum !ev^2
end function gchimp2

pure real(r64) function gCoul2(el, crys, qcrys, evec_k, evec_kp)
!! Function to calculate the Thomas-Fermi screened
!! squared electron-electron vertex.

type(crystal), intent(in) :: crys
type(electron), intent(in) :: el
real(r64), intent(in) :: qcrys(3)
complex(r64),intent(in) :: evec_k(:), evec_kp(:)

real(r64) :: qcart(3), prefac, overlap, qmag

prefac = 1.0e18_r64*crys%volume**3

!Magnitude of q-vector
qmag = twonorm(matmul(crys%reclattvecs, qcrys))

!This is [U(k')U^\dagger(k)]_nm squared
!(Recall that the electron eigenvectors came out daggered from el_wann_epw.)
overlap = (abs(dot_product(evec_kp, evec_k)))**2

!For G, G' = 0
gCoul2 = prefac*overlap* &
(qe**2/perm0/crys%epsilon0/(qmag**2 + crys%qTF**2))**2 !eV^2/nm^3
end function gCoul2

pure real(r64) function Vm2_3ph(ev1_s1, ev2_s2, ev3_s3, &
Index_i, Index_j, Index_k, ifc3, phases_q2q3, ntrip, nb)
Expand Down Expand Up @@ -1285,9 +1310,9 @@ end function speedup_3ph_interactions

subroutine calculate_W3ph_OTF(ph, num, istate1, T, &
Wm, Wp, istate2_plus, istate3_plus, istate2_minus, istate3_minus)
!! On-the-fly (OTF), serial calcualator of the 3-ph transition probability
!! for all IBZ phonon wave vectors. This subroutine calculates
!! W-(s2q2,s3q3) and W+(s2q2,s3q3) for a given irreducible phonon state s1q1.
!! On-the-fly (OTF), serial calcualator of the 3-ph transition probability.
!! This subroutine calculates W-(s2q2,s3q3) and W+(s2q2,s3q3)
!! for a given irreducible phonon state s1q1.
!! The interaction vertex for the given state is read from disk.

type(phonon), intent(in) :: ph
Expand Down Expand Up @@ -2288,6 +2313,147 @@ subroutine calculate_eph_interaction_ibzk(wann, crys, el, ph, num, key)
sync all
end subroutine calculate_eph_interaction_ibzk

subroutine calculate_Xee_OTF(el, num, istate1, T, crys, X)
!! On-the-fly serial calculator of the e-e transition probability.
!! for a given IBZ electron states within the transport window.

type(electron), intent(in) :: el
type(numerics), intent(in) :: num
type(crystal), intent(in) :: crys
real(r64), intent(in) :: T
integer(i64), intent(in) :: istate1
real(r64), intent(out), allocatable :: X(:)

!Local variables
integer(i64) :: nstates_irred, istate, &
n1, ik1, n2, ik2, n3, ik3, n4, ik4, &
count, nprocs
real(r64) :: const, beta, fermi2, fermi3, fermi4, &
delta_val, occup_fac, en1, en2, en3, en4, g2
!integer(i64), allocatable :: istate_el(:), istate_ph(:)
character(len = 1024) :: filename
procedure(delta_fn), pointer :: delta_fn_ptr => null()
type(vec) :: k1_vec, k2_vec, k3_vec, k4_vec, q_vec

!Inverse temperature energy
beta = 1.0_r64/T/kB

!Associate delta function procedure pointer
delta_fn_ptr => get_delta_fn_pointer(num%tetrahedra)

!Constant factor in transition probability expression
const = twopi/hbar_eVps

!Total number of IBZ blocks states
nstates_irred = el%nwv_irred*el%numbands

!Maxium possible length of transition rates
!Nk**3*Nbands**2
nprocs = el%nstates_inwindow**2*el%numbands

!Allocate transition rates
allocate(X(nprocs))

!Initialize X and and the process tallies
X(:) = 0.0_r64

!Demux state index into band (m) and wave vector (ik) indices
call demux_state(istate1, el%numbands, n1, ik1)

!Electron 1 energy
en1 = el%ens_irred(ik1, n1)

!Initialize process counter for the state (k1, n1)
count = 0

!Apply energy window to electron 1
if(abs(en1 - el%enref) <= el%fsthick) then

!Create initial electron wave vector
k1_vec = vec(el%indexlist_irred(ik1), el%wvmesh, crys%reclattvecs)

!Run over electrons states 3 to 4, eliminating the k4 sum with the
!delta(k1 - k3 + k2 - k4)
do ik3 = 1, el%nwv
!Create 3rd electron wave vector
k3_vec = vec(el%indexlist(ik3), el%wvmesh, crys%reclattvecs)

!q \equiv k1 - k3
q_vec = vec_add(k1_vec, k3_vec, el%wvmesh, crys%reclattvecs)

do n3 = 1, el%numbands
!Electron 3 energy
en3 = el%ens(ik3, n3)

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

!Fermi function of electron 3
fermi3 = Fermi(en3, el%chempot, T)

!Squared matrix element
g2 = gCoul2(el, crys, q_vec%frac, &
el%evecs_irred(ik1, n1, :), el%evecs(ik2, n2, :))

do ik2 = 1, el%nwv
!Create 2nd electron wave vector
k2_vec = vec(el%indexlist(ik2), el%wvmesh, crys%reclattvecs)

!Create final electron wave vector
!delta(q + k2 - k4)
k4_vec = vec_add(q_vec, k2_vec, el%wvmesh, crys%reclattvecs)

do n2 = 1, el%numbands
!Electron 2 energy
en2 = el%ens(ik2, n2)

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

!Fermi function of electron 2
fermi2 = Fermi(en2, el%chempot, T)

do n4 = 1, el%numbands
!Electron 4 energy
en4 = el%ens(ik4, n4)

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

!Increase viable process counter
count = count + 1

!Fermi function of electron 4
fermi4 = Fermi(en4, el%chempot, T)

!Evaulate delta function
delta_val = delta_fn_ptr(en4 + en3 - en1, &
ik2, n2, el%wvmesh, el%simplex_map, &
el%simplex_count, el%simplex_evals)

!Temperature dependent occupation factor
!f1.f2.(1 - f3)(1 - f4)/[f1(1 - f1)]
occup_fac = fermi2*(1.0_r64 - fermi3*exp(beta*(en3 - en1)))* &
(1.0_r64 - fermi4)

!Save transition rate
X(count) = g2*occup_fac*delta_val
end do
end do
end do
end do
end do
end if

if(count > 0) then
!Shrink X
call shrink(X, count)

!Multiply constant/units factor, etc.
X = const*X
end if
end subroutine calculate_Xee_OTF

subroutine calculate_echimp_interaction_ibzk(crys, el, num)
!! Parallel driver of |g_e-chimp(k,k')|^2 over IBZ electron states.
!!
Expand Down
88 changes: 62 additions & 26 deletions src/screening.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,9 +117,14 @@ subroutine calculate_Hilbert_weights(w_disc, w_cont, zeroplus, Hilbert_weights)

!!$ integrand = Phi_i/(w_disc(j)**2 - w_cont**2)

!!$ integrand = Phi_i*(&
!!$ 1.0_r64/(w_disc(j) - w_cont - oneI*zeroplus) - &
!!$ 1.0_r64/(w_disc(j) + w_cont + oneI*zeroplus))

!DBG
integrand = Phi_i*(&
1.0_r64/(w_disc(j) - w_cont - oneI*zeroplus) - &
1.0_r64/(w_disc(j) + w_cont + oneI*zeroplus))
1.0_r64/(-w_disc(j) + w_cont - oneI*zeroplus) - &
1.0_r64/(-w_disc(j) - w_cont + oneI*zeroplus))

!TODO Would be nice to have a compsimps function instead of
!a subroutine.
Expand Down Expand Up @@ -176,18 +181,22 @@ subroutine spectral_head_polarizability_3d(spec_eps, Omegas, q_indvec, el, pcell

!Locals
integer(i64) :: m, n, ik, ikp, ikp_window, iOmega, nOmegas, k_indvec(3), kp_indvec(3)
real(r64) :: overlap, ek, ekp
real(r64) :: overlap, ek, ekp, dOmega, delta, Omega_l, Omega_r
procedure(delta_fn), pointer :: delta_fn_ptr => null()

nOmegas = size(Omegas)

dOmega = Omegas(2) - Omegas(1)

allocate(spec_eps(nOmegas))

!TODO don't have to use tetrahedron since I only need the imag part.
!May use triangles also => easily extensible to the 2D case
!
!Associate delta function procedure pointer
delta_fn_ptr => get_delta_fn_pointer(tetrahedra = .true.)
!delta_fn_ptr => get_delta_fn_pointer(tetrahedra = .true.)
!delta_fn_ptr => get_delta_fn_pointer(tetrahedra = .false.)
!Comment: The delta-fn methods above are giving very different numbers.

spec_eps = 0.0
do iOmega = 1, nOmegas
Expand Down Expand Up @@ -222,13 +231,38 @@ subroutine spectral_head_polarizability_3d(spec_eps, Omegas, q_indvec, el, pcell
!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

!overlap = 1.0

!!$ spec_eps(iOmega) = spec_eps(iOmega) + &
!!$ (Fermi(ek, el%chempot, T) - &
!!$ Fermi(ekp, el%chempot, T))*overlap* &
!!$ delta_fn_ptr(ekp - Omegas(iOmega), ik, m, &
!!$ el%wvmesh, el%simplex_map, &
!!$ el%simplex_count, el%simplex_evals)


!DBG
Omega_l = Omegas(iOmega) - dOmega
Omega_r = Omegas(iOmega) + dOmega

if(Omega_l < ekp - ek .and. ekp - ek < Omega_r) continue

delta = finite_element(ekp - ek, &
Omega_l, Omegas(iOmega), Omega_r)
!!$
!!$ if(iOmega == 1) then
!!$ delta = finite_element(ekp - ek, &
!!$ Omegas_l, Omegas(iOmega), Omegas(iOmega + 1))
!!$ else if(iOmega == nOmegas) then
!!$ delta = finite_element(ekp - ek, &
!!$ Omegas(iOmega - 1), Omegas(iOmega), Omegas(iOmega) + dOmega)
!!$ else
!!$ delta = finite_element(ekp - ek, &
!!$ Omegas(iOmega - 1), Omegas(iOmega), Omegas(iOmega + 1))
!!$ end if
spec_eps(iOmega) = spec_eps(iOmega) + &
(Fermi(ek, el%chempot, T) - &
Fermi(ekp, el%chempot, T))*overlap* &
delta_fn_ptr(ekp - Omegas(iOmega), ik, m, &
el%wvmesh, el%simplex_map, &
el%simplex_count, el%simplex_evals)
Fermi(ekp, el%chempot, T))*overlap*delta/product(el%wvmesh)
end do
end do
end do
Expand All @@ -241,11 +275,7 @@ subroutine spectral_head_polarizability_3d(spec_eps, Omegas, q_indvec, el, pcell
!!$ -spec_eps(iOmega)*pi*sign(1.0_r64, Omegas(iOmega))*el%spindeg/pcell_vol

spec_eps(iOmega) = &
spec_eps(iOmega)*sign(1.0_r64, Omegas(iOmega))*el%spindeg/pcell_vol

!!$ !DBG
!!$ spec_eps(iOmega) = &
!!$ spec_eps(iOmega)*sign(1.0_r64, Omegas(iOmega))*el%spindeg/pcell_vol
-spec_eps(iOmega)*sign(1.0_r64, Omegas(iOmega))*el%spindeg/pcell_vol

!Zero out extremely small numbers !TODO What is small? In what units?
!if(abs(spec_eps(iOmega)) < 1.0e-30_r64) spec_eps(iOmega) = 0.0_r64
Expand All @@ -269,7 +299,7 @@ subroutine calculate_RPA_dielectric_3d_G0_qpath(el, crys, num)

!Locals
real(r64), allocatable :: energylist(:), qlist(:, :), qmaglist(:)
real(r64) :: qcrys(3)
real(r64) :: qcrys(3), zeroplus
integer(i64) :: iq, iOmega, numomega, numq, &
start, end, chunk, num_active_images, qxmesh
real(r64), allocatable :: spec_eps(:)
Expand All @@ -279,9 +309,9 @@ subroutine calculate_RPA_dielectric_3d_G0_qpath(el, crys, num)

!Silicon
!omega_plasma = 1.0e-9*hbar*sqrt(el%conc_el/perm0/crys%epsilon0/(0.267*me)) !eV

!wGaN
omega_plasma = 1.0e-9*hbar*sqrt(el%conc_el/perm0/crys%epsilon0/(0.2*me)) !eV
!omega_plasma = 1.0e-9*hbar*sqrt(el%conc_el/perm0/crys%epsiloninf/(0.2*me)) !eV
omega_plasma = 1.0e-9*hbar*sqrt(el%conc_el/perm0/crys%epsiloninf/(0.22*me)) !eV

if(this_image() == 1) then
print*, "plasmon energy = ", omega_plasma
Expand Down Expand Up @@ -316,15 +346,18 @@ subroutine calculate_RPA_dielectric_3d_G0_qpath(el, crys, num)
allocate(qlist(numq, 3), qmaglist(numq))
do iq = 1, numq
qlist(iq, :) = el%wavevecs_irred(iq, :)
!qmaglist(iq) = twonorm(matmul(crys%reclattvecs, qlist(iq, :)))
qmaglist(iq) = qdist(qlist(iq, :), crys%reclattvecs)
qmaglist(iq) = twonorm(matmul(crys%reclattvecs, qlist(iq, :)))
!qmaglist(iq) = qdist(qlist(iq, :), crys%reclattvecs)
end do

!Create energy grid
numomega = 500
numomega = 200
allocate(energylist(numomega))
call linspace(energylist, 0.0_r64, 0.4_r64, numomega)
call linspace(energylist, 0.0_r64, 0.2_r64, numomega)

!Small number
zeroplus = (energylist(2) - energylist(1))*0.1

!Allocate diel_ik to hold maximum possible Omega
allocate(diel(numq, numomega))
diel = 0.0_r64
Expand Down Expand Up @@ -356,7 +389,7 @@ subroutine calculate_RPA_dielectric_3d_G0_qpath(el, crys, num)
call calculate_Hilbert_weights(&
w_disc = energylist, &
w_cont = energylist, &
zeroplus = 1.0e-6_r64, & !Can this magic "small" number be removed?
zeroplus = zeroplus, & !Can this magic "small" number be removed?
Hilbert_weights = Hilbert_weights)

!call Re_head_polarizability_3d_T(Reeps, energylist, Imeps, Hilbert_weights)
Expand All @@ -380,19 +413,22 @@ subroutine calculate_RPA_dielectric_3d_G0_qpath(el, crys, num)
!energy expression requires a single effective mass.
!Just leaving it here for now to get the plasmon peak in the loss function.
diel(iq, 1:numomega) = 1.0_r64 - &
(omega_plasma/energylist(1:numomega))**2 + oneI*1.0e-9
(omega_plasma/energylist(1:numomega))**2 !+ oneI*1.0e-9
else
!DBG
diel(iq, :) = 1.0_r64 - &
diel(iq, :) = crys%epsiloninf - &
1.0_r64/qmaglist(iq)**2* &
(real(eps) + oneI*pi*(-spec_eps))/perm0*qe*1.0e9_r64
(real(eps) + oneI*pi*(spec_eps))/perm0*qe*1.0e9_r64
!!$
!!$ diel(iq, :) = crys%epsiloninf - &
!!$ 1.0_r64/qmaglist(iq)**2*eps/perm0*qe*1.0e9_r64
end if
end do

call co_sum(diel)

!Handle Omega = 0 case
diel(:, 1) = 1.0_r64 + (crys%qTF/qmaglist(:))**2
!diel(:, 1) = 1.0_r64 + (crys%qTF/qmaglist(:))**2

!Print to file
call write2file_rank2_real("RPA_dielectric_3D_G0_qpath", qlist)
Expand Down
Loading

0 comments on commit 6de799d

Please sign in to comment.