Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixed a bug in Matrix element module and added r^2 integral #153

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Fixed a bug in Matrix element module and added r^2 integral #153

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
61a120b
r3 integral corrected r2 integral added
Jul 1, 2024
6e01b20
Merge branch 'Matrix_Element' of https://github.com/MuyamiYatara/wann…
Jul 1, 2024
4ad12d8
unfolding.f90 bug fixed
Jul 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 50 additions & 9 deletions src/module.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1460,7 +1460,6 @@ end subroutine factorial

end module element_table


!> module used for matrix element calculation for simulation of ARPES experiment
!-----------------------------------------------------------------------------------------------------!
! implemented functions
Expand Down Expand Up @@ -1848,7 +1847,7 @@ end function Ylm


!> Generate the integrand function
real(dp) function integrand(r, k_f, name, n, l)
real(dp) function integrand_r2(r, k_f, name, n, l)


implicit none
Expand All @@ -1865,8 +1864,31 @@ real(dp) function integrand(r, k_f, name, n, l)

call Z_eff_cal(name, n, l, Z_eff)
res = (r**2)*slater_Rn(r, n, Z_eff)*spherical_bessel(k_f*r, l)
integrand = res
end function integrand
integrand_r2 = res
end function integrand_r2


!> Generate the integrand function
real(dp) function integrand_r3(r, k_f, name, n, l)


implicit none
character(len=10) ,intent(in) :: name
real(dp), intent(in) :: r
real(dp), intent(in) :: k_f
integer, intent(in) :: n
integer, intent(in) :: l


real(dp) :: Z_eff ! effective nuclear charge number under slater rules
real(dp) :: res


call Z_eff_cal(name, n, l, Z_eff)
res = (r**3)*slater_Rn(r, n, Z_eff)*spherical_bessel(k_f*r, l)
integrand_r3 = res
end function integrand_r3



!> Simpson integral
Expand Down Expand Up @@ -1993,15 +2015,16 @@ end subroutine cartesian_to_spherical


!> Main subroutine to calculate the Matrix Element
subroutine get_matrix_element(ele_name, orb_name, k_cart, matrix_element_res)
subroutine get_matrix_element(ele_name, orb_name, k_cart, atom_position_cart, matrix_element_res)


implicit none
character(len=10), intent(in) :: ele_name
character(len=10), intent(in) :: orb_name
real(dp) :: k_cart(3) ! k in cartesian coordinate
real(dp) :: atom_position_cart(3) ! position of calculated atom
real(dp) :: gaunt_value_tmp
complex(dp) :: Ylm_gaunt_vec(3) ! values of gaunt coefficient plus Ylm corresponding to polarization vector
complex(dp) :: Ylm_gaunt_vec(3) ! values of gaunt coefficient multiplied by Ylm corresponding to polarization vector
complex(dp) :: polar_epsilon_vec(3) ! polarization vector


Expand All @@ -2012,7 +2035,7 @@ subroutine get_matrix_element(ele_name, orb_name, k_cart, matrix_element_res)
integer :: dl, m_3 ! dl values from {-1, +1}

real(dp) :: Z_eff
real(dp) :: inte_res, inte_a = 0d0
real(dp) :: inte_res_r3, inte_res_r2, inte_a = 0d0
real(dp) :: k_theta, k_phi, k_f ! k under spherical coordinate, k_f in the unit of angstrom^-1, k_theta and k_phi are in radian
!real(dp) :: matrix_element_res
complex(dp) :: matrix_element_res
Expand Down Expand Up @@ -2042,10 +2065,13 @@ subroutine get_matrix_element(ele_name, orb_name, k_cart, matrix_element_res)
call Z_eff_cal(ele_name, n, l, Z_eff)

matrix_element_res = 0d0


do j = 0, 1
! generate the dl
dl = j*2-1
call integral_simpson(inte_a, inte_b, integrand, k_f, ele_name, n, l+dl, inte_res)
inte_res_r3 = 0d0
call integral_simpson(inte_a, inte_b, integrand_r3, k_f, ele_name, n, l+dl, inte_res_r3)



Expand All @@ -2060,14 +2086,29 @@ subroutine get_matrix_element(ele_name, orb_name, k_cart, matrix_element_res)

Ylm_gaunt_vec(3) = gaunt_value_tmp

!> calculate the term comes from the atom position








!> calculate the matrix element

matrix_element_res = matrix_element_res + (zi)**(l+dl)*4d0*Pi*inte_res*dot_product(polar_epsilon_vec, Ylm_gaunt_vec)
matrix_element_res = matrix_element_res + (zi)**(l+dl)*4d0*Pi*inte_res_r3*dot_product(polar_epsilon_vec, Ylm_gaunt_vec)



end do ! j loop

inte_res_r2 = 0d0
call integral_simpson(inte_a, inte_b, integrand_r2, k_f, ele_name, n, l, inte_res_r2)

matrix_element_res = matrix_element_res + (zi)**l*4d0*Pi*inte_res_r2*dot_product(polar_epsilon_vec, atom_position_cart)&
*Ylm(k_theta, k_phi, l, m)

return
end subroutine get_matrix_element

Expand Down
5 changes: 3 additions & 2 deletions src/unfolding.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -763,6 +763,7 @@ subroutine get_projection_weight_bulk_unfold(ndim, k_SBZ_direct, k_PBZ_direct, p
k_t=k_PBZ_direct-k_SBZ_direct_in_PBZ

allocate(me_values(Folded_cell%NumberofSpinOrbitals))

if ( Matrix_Element_calc == .True. ) then
!@ k_abs is the k_f(3) considering the photon energy
if ( (k_cart_abs**2 - k_cart(1)**2 - k_cart(2)**2 ) .le. 0 ) then
Expand All @@ -783,9 +784,9 @@ subroutine get_projection_weight_bulk_unfold(ndim, k_SBZ_direct, k_PBZ_direct, p

projector_name_PC= adjustl(trim(Folded_cell%proj_name(Folded_cell%spinorbital_to_projector_index(io_PC), Folded_cell%spinorbital_to_atom_index(io_PC))))
me = 0d0 !! initialize the variable
call get_matrix_element(atom_name_PC, projector_name_PC, k_cart, me)
call get_matrix_element(atom_name_PC, projector_name_PC, k_cart, Folded_cell%wannier_centers_direct(:, io_PC), me)
me_values(io_PC) = me

enddo ! io_PC


Expand Down
Loading