Skip to content

Commit

Permalink
Implementation of hilbert transform
Browse files Browse the repository at this point in the history
  • Loading branch information
@dpaulzc
dpaulzc committed Oct 24, 2024
1 parent 0d49458 commit c97c627
Showing 1 changed file with 89 additions and 15 deletions.
104 changes: 89 additions & 15 deletions src/screening.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,7 @@ subroutine calculate_Hilbert_weights(w_disc, w_cont, zeroplus, Hilbert_weights)
real(r64), intent(in) :: w_disc(:), w_cont(:), zeroplus
!complex(r64), allocatable, intent(out) :: Hilbert_weights(:, :)
real(r64), allocatable, intent(out) :: Hilbert_weights(:, :)
real(r64), allocatable :: Hilbert_weights_i(:, :)

!Locals
integer :: n_disc, n_cont, i, j
Expand All @@ -96,17 +97,18 @@ subroutine calculate_Hilbert_weights(w_disc, w_cont, zeroplus, Hilbert_weights)
n_disc = size(w_disc)
n_cont = size(w_cont)

allocate(Hilbert_weights(n_disc, n_disc))
allocate(Hilbert_weights(n_disc, n_disc),Hilbert_weights_i(n_disc, n_disc))

!Grid spacing of "continuous" variable
dw = w_cont(2) - w_cont(1)
! # CHANGE #
open(1, file = 'hilbertweightcheck')
open(2, file = 'omega-range')
do i = 1, n_disc
w_l = w_disc(i) - dw
w_r = w_disc(i) + dw
!Phi_i = finite_element(w_cont, w_l, w_disc(i), w_r)
Phi_i = 1
Phi_i = finite_element(w_cont, w_l, w_disc(i), w_r)
!Phi_i = 1
do j = 1, n_disc
integrand = Phi_i*(&
(w_disc(j) - w_cont + oneI*zeroplus)/((w_disc(j) - w_cont)**2 + zeroplus**2) - &
Expand All @@ -123,15 +125,79 @@ subroutine calculate_Hilbert_weights(w_disc, w_cont, zeroplus, Hilbert_weights)
!a subroutine.
!call compsimps(integrand, dw, Hilbert_weights(j, i))
call compsimps(real(integrand), dw, realpart)
call compsimps(imag(integrand), dw, imagpart)

Hilbert_weights(j, i) = realpart
Hilbert_weights_i(j, i) = imagpart
end do
end do
! # CHANGE #
write(1,*) Hilbert_weights(:,1)
write(1,*) Hilbert_weights_i(:,1)
write(2,*) w_disc(:)
close(1)
close(2)
end subroutine calculate_Hilbert_weights

subroutine hilbert_transform(Reeps_T, Omegas, spec_eps_T)
!! EQ (4), R. Balito et. al.
!!
!! Reeps_T Real part of bare polarizability
!! Omega Energy of excitation in the electron gas
!! spec_eps_T Spectral head of bare polarizability
!!

real(r64), allocatable, intent(out) :: Reeps_T(:)
real(r64), intent(in) :: Omegas(:)
real(r64), intent(in) :: spec_eps_T(:)

! Local variables
real(r64) :: in1, in2, fin, bl !, domegas
integer(i64) :: nx, kh, np

nx = size(Omegas)
allocate(Reeps_T(nx))

!domegas = Omegas(2)-Omegas(1)

do kh = 1, nx
in1=0.0_r64
in2=0.0_r64
do np=1,nx-1-kh
bl = log(real(np+1.0_r64)/real(np))
if ( (kh+np).lt.nx) then
in1 = in1 - (1.0_r64-(np+1.0_r64)*bl)*spec_eps_T(kh+np) + &
(1.0_r64-np*bl)*spec_eps_T(kh+np+1)
else
!print *,"Hilbert transform: Boundary breached",kh+np,nx
in1 = in1 - (1.0_r64-(np+1.0_r64)*bl)*spec_eps_T(kh+np) + &
(1.0_r64-np*bl)*1e-12_r64
end if
end do

do np=1,kh-1
bl = log(real(np+1.0_r64)/real(np))
if ((kh-np).gt.1) then
in2 = in2 + (1.0_r64-(np+1.0_r64)*bl)*spec_eps_T(kh-np) - &
(1.0_r64-np*bl)*spec_eps_T(kh-np-1)
else
!print *,"Hilbert transform: Boundary breached",kh+np,nx
in2 = in2 + (1.0_r64-(np+1.0_r64)*bl)*spec_eps_T(kh-np) - &
(1.0_r64-np*bl)*1e-12_r64
end if
end do

if (kh.ge.2 .and. kh.le.(nx-1)) then
fin = spec_eps_T(kh+1)-spec_eps_T(kh-1)
else
fin = 1e-12_r64
end if
Reeps_T(kh) = (-1.0_r64/pi)*(fin+in1+in2)

end do

end subroutine hilbert_transform

subroutine head_polarizability_real_3d_T(Reeps_T, Omegas, spec_eps_T, Hilbert_weights_T)
!! Head of the bare real polarizability of the 3d Kohn-Sham system using
!! Hilbert transform for a given set of temperature-dependent quantities.
Expand Down Expand Up @@ -428,12 +494,12 @@ subroutine calculate_RPA_dielectric_3d_G0_scratch(el, crys, num, wann)
complex(r64), allocatable :: diel(:, :)
character(len = 1024) :: filename
real(r64) :: omega_plasma

integer(i64) :: ij
!Silicon
!omega_plasma = 1.0e-9_r64*hbar*sqrt(el%conc_el/perm0/crys%epsiloninf/(0.267*me)) !eV

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

if(this_image() == 1) then
print*, "plasmon energy = ", omega_plasma
Expand Down Expand Up @@ -511,18 +577,26 @@ subroutine calculate_RPA_dielectric_3d_G0_scratch(el, crys, num, wann)

call spectral_head_polarizability_3d_qpath(&
spec_eps, energylist, qcrys, el, wann, crys, num%tetrahedra)

!Calculate re_eps with Hilbert-Kramers-Kronig transform
call calculate_Hilbert_weights(&
w_disc = energylist, &
w_cont = energylist, &
zeroplus = zeroplus, & !Can this magic "small" number be removed?
Hilbert_weights = Hilbert_weights)

call head_polarizability_real_3d_T(Reeps, energylist, spec_eps, Hilbert_weights)

call head_polarizability_imag_3d_T(Imeps, energylist, energylist, spec_eps)
!Calculate re_eps with Hilbert-Kramers-Kronig transform
! call calculate_Hilbert_weights(&
! w_disc = energylist, &
! w_cont = energylist, &
! zeroplus = zeroplus, & !Can this magic "small" number be removed?
! Hilbert_weights = Hilbert_weights)
open(10,file = "hilbertw")
open(20,file = "spec_eps")
!call head_polarizability_real_3d_T(Reeps, energylist, spec_eps, Hilbert_weights)

call hilbert_transform(Reeps, energylist, spec_eps)

call head_polarizability_imag_3d_T(Imeps, energylist, energylist, spec_eps)
do ij=1,numomega
write(10,*) energylist(ij),Reeps(ij),Imeps(ij)
write(20,*) energylist(ij),spec_eps(ij)
end do
close(10)
close(20)
!Calculate RPA dielectric (diagonal in G-G' space)
!!$ diel(iq, :) = 1.0_r64 - &
!!$ (1.0_r64/twonorm(matmul(crys%reclattvecs, qcrys)))**2* &
Expand Down

0 comments on commit c97c627

Please sign in to comment.