diff --git a/src/module.f90 b/src/module.f90 index 8369200..132e0a5 100644 --- a/src/module.f90 +++ b/src/module.f90 @@ -1474,7 +1474,6 @@ end subroutine factorial end module element_table - !> module used for matrix element calculation for simulation of ARPES experiment !-----------------------------------------------------------------------------------------------------! ! implemented functions @@ -1862,7 +1861,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 @@ -1879,8 +1878,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 @@ -2007,15 +2029,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 @@ -2026,7 +2049,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 @@ -2056,10 +2079,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) @@ -2074,14 +2100,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 diff --git a/src/sigma_OHE.f90 b/src/sigma_OHE.f90 index f6588f8..0c91b5b 100644 --- a/src/sigma_OHE.f90 +++ b/src/sigma_OHE.f90 @@ -271,7 +271,8 @@ subroutine sigma_ohe_calc_symm(mu_array, KBT_array, BTau_array, Nband_Fermi_Leve !> we get the kpath by Btau_final=-exponent_max*BTauMax, but we only use half of them !> means that we can reach the accuracy as to exp(-exponent_max) - exponent_max= 30 + ! exponent_max= 30 + exponent_max= 15 !> exclude all kpoints with zero velocity x B and large energy away from Fermi level do iband= 1, Nband_Fermi_Level @@ -964,13 +965,17 @@ subroutine cal_sigma_iband_k BTau= BTau_array(ibtau) if (NBTau==1)then - NSlice_Btau_local= 2 + ! NSlice_Btau_local= 2 + NSlice_Btau_local= 1 else - NSlice_Btau_local= (ibtau-1)*NSlice_Btau_inuse/(NBTau-1)/2 + ! NSlice_Btau_local= (ibtau-1)*NSlice_Btau_inuse/(NBTau-1)/2 + NSlice_Btau_local= (ibtau-1)*NSlice_Btau_inuse/(NBTau-1) if (NSlice_Btau_local==0)then - NSlice_Btau_local= 2 + ! NSlice_Btau_local= 2 + NSlice_Btau_local= 1 else - DeltaBtau= exponent_max/2d0/NSlice_Btau_local + ! DeltaBtau= exponent_max/2d0/NSlice_Btau_local + DeltaBtau= exponent_max/NSlice_Btau_local endif endif @@ -980,21 +985,41 @@ subroutine cal_sigma_iband_k !> The core of Chamber formular is to get the average of velocity !> in the relaxation time approximation v_k= klist_iband(iband)%velocity_k(:, 1+ishift) - if (BTau>eps3.and.NSlice_Btau_local>2) then - ! set weight - do it=1, NSlice_Btau_local - weights(it) = exp(-DeltaBtau*(it-1d0)) - enddo - !> trapezoidal integral + ! if (BTau>eps3.and.NSlice_Btau_local>2) then + if (BTau>eps3.and.NSlice_Btau_local>1) then velocity_bar_k= 0d0 - do it=1, NSlice_Btau_local + !> five point integration(n=4) + do it=1, NSlice_Btau_local, 4 !(hj, j=it-1, j mod 4==0, Aj=28/45) velocity_bar_k= velocity_bar_k+ & - DeltaBtau*exp(-(it-1d0)*DeltaBtau)*klist_iband(iband)%velocity_k(:, it+ishift) + 28.0d0/45.0d0*DeltaBtau*exp(-(it-1d0)*DeltaBtau)*klist_iband(iband)%velocity_k(:, it+ishift) + enddo + do it=2, NSlice_Btau_local, 4 !(hj, j=it-1, j mod 4==1, Aj=64/45) + velocity_bar_k= velocity_bar_k+ & + 64.0d0/45.0d0*DeltaBtau*exp(-(it-1d0)*DeltaBtau)*klist_iband(iband)%velocity_k(:, it+ishift) enddo - velocity_bar_k= velocity_bar_k & - - 0.5d0*DeltaBtau*(exp(-(NSlice_Btau_local-1d0)*DeltaBtau)& - * klist_iband(iband)%velocity_k(:, NSlice_Btau_local+ishift) & - + klist_iband(iband)%velocity_k(:, 1+ishift)) + do it=3, NSlice_Btau_local, 4 !(hj, j=it-1, j mod 4==2, Aj=24/45) + velocity_bar_k= velocity_bar_k+ & + 24.0d0/45.0d0*DeltaBtau*exp(-(it-1d0)*DeltaBtau)*klist_iband(iband)%velocity_k(:, it+ishift) + enddo + do it=4, NSlice_Btau_local, 4 !(hj, j=it-1, j mod 4==3, Aj=64/45) + velocity_bar_k= velocity_bar_k+ & + 64.0d0/45.0d0*DeltaBtau*exp(-(it-1d0)*DeltaBtau)*klist_iband(iband)%velocity_k(:, it+ishift) + enddo + velocity_bar_k= velocity_bar_k-14.0d0/45.0d0*DeltaBtau*klist_iband(iband)%velocity_k(:, 1+ishift) + ! set weight + ! do it=1, NSlice_Btau_local + ! weights(it) = exp(-DeltaBtau*(it-1d0)) + ! enddo + ! !> trapezoidal integral + ! velocity_bar_k= 0d0 + ! do it=1, NSlice_Btau_local + ! velocity_bar_k= velocity_bar_k+ & + ! DeltaBtau*exp(-(it-1d0)*DeltaBtau)*klist_iband(iband)%velocity_k(:, it+ishift) + ! enddo + ! velocity_bar_k= velocity_bar_k & + ! - 0.5d0*DeltaBtau*(exp(-(NSlice_Btau_local-1d0)*DeltaBtau)& + ! * klist_iband(iband)%velocity_k(:, NSlice_Btau_local+ishift) & + ! + klist_iband(iband)%velocity_k(:, 1+ishift)) !> Simpson's integral !> add the first and the last term diff --git a/src/unfolding.f90 b/src/unfolding.f90 index c41293d..2533ef0 100644 --- a/src/unfolding.f90 +++ b/src/unfolding.f90 @@ -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 @@ -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