Hello.
I want to modify some codes for other calculations.
But I don't know much about fortran.
I'm trying to vmef function for vkk, vkq calculation at electron self energy calculation.
but I was confused about the unit and printing result.
some values are xxxxxx+E0x, but some values are xxxxxxx-3xx without E value.
what does it mean -3xx?
And i want to know how to use dnkf functions in eletron self energy module.
last, i want to know about write each vkq component.
Could you help me check the code?
file means ik, ibnd, E, ikk, vkk_1,vkk_2,vkk_3, vkk^2, Im(sigma)
# Electron linewidth = 2*Im(Sigma) (meV)
# ik ibnd E(ibnd) Im(Sigma)(meV)
1 3 0.20912459750849E+00 1 0.11217464632496E+00 0.19475320066878E+00 0.50654454304769249567E-02 -0.22464337352605E+00 0.88781315185286E+02
1 4 0.20912472459197E+00 1 0.13474681982647E+00 0.76596127338991E-01 0.45484610181721646383E+00 -0.12327685904623E+00 0.88781315185286E+02
1 5 0.80002950283306E+01 1 -0.33223079535764E+00 -0.33031245280955E-01 0.20368723164100269174E+00 -0.25186638397400E+00 0.17744735686939E+01
2 3 0.29332887962885E+00 3 0.34092210789731E-01 0.35089491672347E-01 -0.11107861198297262761E+00 0.42410236285429E+00 0.18414035861166E+03
2 4 0.29332900712137E+00 3 -0.87962420861798E-01 0.64408945581074E-01 0.45252503339206739952E+00 0.17048275913375E+00 0.18414035861166E+03
2 5 0.79546643124405E+01 3 0.31036395000744E-01 0.42528305536141E-01 0.42395705933447969116E+01 0.11642600457680E+00 0.43204361296634E+01
3 3 0.46300776118313E+00 5 0.11488646468214E+00 -0.28507430695821E+00 0.92762945650373934692E-01 0.73324063594341E-01 0.12960549021292E+03
3 4 0.46300789093205E+00 5 0.13089499757876E+00 -0.22464337352605E+00 -0.12327685904623308932E+00 -0.47233087182556E+00 0.12960549021292E+03
3 5 0.78471330866398E+01 5 -0.25186638397400E+00 0.14645585997667E-01 -0.85398764031314522427E-01 -0.32348079240387E+00 0.26887337102817E+01
4 3 0.55918222425530E+00 7 -0.32348079240387E+00 -0.17568563379501E+00 0.00000000000000000000E+00 0.00000000000000E+00 0.12124205219744E+03
4 4 0.55918235556153E+00 7 0.00000000000000E+00 0.00000000000000E+00 0.19762625833649861767-321 0.19762625833650-321 0.12124205219744E+03
4 5 0.77944129649239E+01 7 -0.13295284211842E+01 0.51224726160820-319 0.85374543601367402834-319 -0.13295284211842E+01 0.63916610909211E+00
5 3 0.46300776118335E+00 9 0.93412952133953-316 0.83991159793012-322 0.69533558053346369017-309 0.44759381595362E-90 0.12962122780150E+03
5 4 0.46300789093175E+00 9 0.29643938750475-322 0.93511330485353-316 0.49406564584124654418-323 0.44759381595362E-90 0.12962122780150E+03
5 5 0.78471330866460E+01 9 0.98813129168249-322 0.12789802275914-315 0.15000000000000000000E+02 0.44759381595362E-90 0.26887440431509E+01
6 3 0.29332887962925E+00 11 0.26177339447724-309 0.00000000000000E+00 0.00000000000000000000E+00 0.44759381595362E-90 0.18585601832194E+03
6 4 0.29332900712142E+00 11 0.31620201333840-321 0.93414770295530-316 0.19762625833649861767-321 0.44759381595362E-90 0.18585601832194E+03
6 5 0.79546643124653E+01 11 0.93413050947082-316 0.30760289958566-315 0.39525251667299723534-322 0.93412952133953-316 0.43204442780937E+01
7 3 -0.18133857440905E+01 13 0.49406564584125-323 0.30755151675849-315 0.39525251667299723534-322 0.93511330485353-316 0.83475414532130E+02
7 4 -0.71266064322113E+00 13 0.00000000000000E+00
...
Thank you for your help
Hello. I have a question for vmef, dmef function, and I want to implement in electron self energy calculation.
Moderator: stiwari
-
darjeeling
- Posts: 2
- Joined: Fri Mar 25, 2022 5:18 pm
- Affiliation: Kyung Hee University
Re: Hello. I have a question for vmef, dmef function, and I want to implement in electron self energy calculation.
Dear darjeeling:
It is not clear to me which quantities you want to calculate, but dnkf is not related to the electron self-energy.
It helps if you provide your code and you let me know which quantities you want to calculate.
Sincerely,
H. Lee
I guess that there is some formatting issue and in this case, I think that -3xx means E-3xx, that is, 10^(-3xx).some values are xxxxxx+E0x, but some values are xxxxxxx-3xx without E value.
what does it mean -3xx?
It is not clear to me which quantities you want to calculate, but dnkf is not related to the electron self-energy.
It helps if you provide your code and you let me know which quantities you want to calculate.
Sincerely,
H. Lee
-
darjeeling
- Posts: 2
- Joined: Fri Mar 25, 2022 5:18 pm
- Affiliation: Kyung Hee University
Re: Hello. I have a question for vmef, dmef function, and I want to implement in electron self energy calculation.
Dear hlee.
Thank you for your comment.
At first, I thought -3xx was -E3xx.
What I doubt is that when vkk^2 is a combination of -E0x and -3xx of vkk, there are parts that become -3xx.
And I want to calculate momentum relaxed electron-phonon relaxation time value.
There are equations of (4) and (5) in https://doi.org/10.1038/s41567-021-01341-w
So, I'm trying to modifying epw code for this calculation.
In below, i added vkk, vkq, vkksquare and dfnk. but dfnk doesn't work.
Finally, the value you want to calculate is {integral{(vkk^2*tau)dk}/{integral(vkk^2)dk} as shown in Equation (5).
Could you help me put the formula into the internal code?
!
! Copyright (C) 2016-2019 Samuel Ponce', Roxana Margine, Feliciano Giustino
! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino
! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino
!
! This file is distributed under the terms of the GNU General Public
! License. See the file `LICENSE' in the root directory of the
! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
!
!----------------------------------------------------------------------
MODULE selfen
!----------------------------------------------------------------------
!!
!! This module contains the various self-energy routines
!!
IMPLICIT NONE
!
CONTAINS
!
!-----------------------------------------------------------------------
SUBROUTINE selfen_elec_q(iqq, iq, totq, first_cycle)
!-----------------------------------------------------------------------
!!
!! Compute the imaginary part of the electron self energy due to electron-
!! phonon interaction in the Migdal approximation. This corresponds to
!! the electron linewidth (half width). The phonon frequency is taken into
!! account in the energy selection rule.
!!
!! Use matrix elements, electronic eigenvalues and phonon frequencies
!! from ep-wannier interpolation
!!
!! This routines computes the contribution from phonon iq to all k-points
!! The outer loop in ephwann_shuffle.f90 will loop over all iq points
!! The contribution from each iq is summed at the end of this subroutine
!! for iqq=totq to recover the per-ik electron self energy
!!
!! RM 24/02/2014
!! Redefined the size of sigmar_all, sigmai_all, and zi_all within the fermi windwow
!!
!-----------------------------------------------------------------------
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE io_var, ONLY : linewidth_elself
USE modes, ONLY : nmodes
USE epwcom, ONLY : nstemp, nbndsub, shortrange, fsthick, ngaussw, degaussw, &
eps_acustic, efermi_read, fermi_energy, restart, restart_step
USE pwcom, ONLY : ef
USE elph2, ONLY : etf, ibndmin, ibndmax, nkqf, xqf, eta, nbndfst, &
nkf, epf17, wf, wqf, xkf, nkqtotf, adapt_smearing, &
sigmar_all, sigmai_all, sigmai_mode, zi_all, efnew, &
nktotf, lower_bnd, gtemp, &
vmef, dmef
USE control_flags, ONLY : iverbosity
USE constants_epw, ONLY : kelvin2eV, ryd2mev, ryd2ev, one, two, zero, ci, eps4, eps6, eps8, hbar
USE constants, ONLY : pi
USE mp, ONLY : mp_barrier, mp_sum
USE mp_global, ONLY : inter_pool_comm
USE mp_world, ONLY : mpime
USE io_global, ONLY : ionode_id
USE io_selfen, ONLY : selfen_el_write
USE poolgathering, ONLY : poolgather2
!
IMPLICIT NONE
!
CHARACTER(LEN = 20) :: tp
!! string for temperatures
CHARACTER(LEN = 256) :: fileselfen
!! file name for self energy
LOGICAL, INTENT(inout) :: first_cycle
!! Use to determine weather this is the first cycle after restart
INTEGER, INTENT(in) :: iqq
!! Q-point index from selecq.fmt window
INTEGER, INTENT(in) :: iq
!! Q-point index from full grid
INTEGER, INTENT(in) :: totq
!! Total number of q-points from the selecq.fmt grid.
!
! Local variables
INTEGER :: n
!! Integer for the degenerate average over eigenstates
INTEGER :: ik
!! Counter on the k-point index
INTEGER :: ikk
!! k-point index
INTEGER :: ikq
!! q-point index
INTEGER :: ibnd
!! Counter on bands at k
INTEGER :: jbnd
!! Counter on bands at k+q
INTEGER :: imode
!! Counter on mode
INTEGER :: fermicount
!! Number of states on the Fermi surface
INTEGER :: itemp
!! Counter on temperatures
INTEGER :: ierr
!! Error status
!
REAL(KIND = DP) :: g2
!! Electron-phonon matrix elements squared in Ry^2
REAL(KIND = DP) :: ef0
!! Fermi energy level
REAL(KIND = DP) :: ekk
!! Eigen energy at k on the fine grid relative to the Fermi level
REAL(KIND = DP) :: ekk1
!! Temporary variable to the eigenenergies at k for the degenerate average
REAL(KIND = DP) :: ekq
!! Eigen energy at k+q on the fine grid relative to the Fermi level
REAL(KIND = DP) :: etmp1
!! Temporary variable to store etmp1 = ekk - (ekq - wq)
REAL(KIND = DP) :: etmp2
!! Temporary variable to strore etmp2 = ekk - (ekq + wq)
REAL(KIND = DP) :: sq_etmp1
!! Temporary variable to store etmp1^2
REAL(KIND = DP) :: sq_etmp2
!! Temporary variable to store etmp2^2
REAL(KIND = DP) :: wgkq
!! Fermi-Dirac occupation factor $f_{mk+q}(T)$
REAL(KIND = DP) :: fact1
!! Temporary variable to store $f_{mk+q}(T) + n_{q\nu}(T)$
REAL(KIND = DP) :: fact2
!! Temporary variable to store $1 - f_{mk+q}(T) + n_{q\nu}(T)$
REAL(KIND = DP) :: weight
!! Self-energy factor
!!$$ N_q \Re(\frac{f_{mk+q}(T) + n_{q\nu}(T)}{\varepsilon_{nk} - \varepsilon_{mk+q} + \omega_{q\nu} - i\delta}) $$
!!$$ + N_q \Re(\frac{1 - f_{mk+q}(T) + n_{q\nu}(T)}{\varepsilon_{nk} - \varepsilon_{mk+q} - \omega_{q\nu} - i\delta}) $$
REAL(KIND = DP) :: w0g1
!! Dirac delta at k for the imaginary part of $\Sigma$
REAL(KIND = DP) :: w0g2
!! Dirac delta at k+q for the imaginary part of $\Sigma$
REAL(KIND = DP) :: inv_degaussw
!! Inverse of degaussw define for efficiency reasons
REAL(KIND = DP) :: inv_eptemp
!! Inverse of temperature define for efficiency reasons
REAL(KIND = DP) :: eta_tmp
!! Temporary variable eta2
REAL(KIND = DP) :: sq_eta_tmp
!! Temporary eta2^2
REAL(KIND = DP) :: inv_eta_tmp
!! Temporary varialbe inv_eta
REAL(KIND = DP) :: tmp1
!! Temporary variable to store real part of Sigma for the degenerate average
REAL(KIND = DP) :: tmp2
!! Temporary variable to store imag part of Sigma for the degenerate average
REAL(KIND = DP) :: tmp3
!! Temporary variable to store Z for the degenerate average
REAL(KIND = DP) :: DDOT
!! Dot product function
REAL(KIND = DP) :: sigmar_tmp(nbndfst)
!! Temporary array to store the real-part of Sigma
REAL(KIND = DP) :: sigmai_tmp(nbndfst)
!! Temporary array to store the imag-part of Sigma
REAL(KIND = DP) :: zi_tmp(nbndfst)
!! Temporary array to store the Z
REAL(KIND = DP), EXTERNAL :: wgauss
!! Fermi-Dirac distribution function (when -99)
REAL(KIND = DP), EXTERNAL :: w0gauss
!! This function computes the derivative of the Fermi-Dirac function
!! It is therefore an approximation for a delta function
REAL(KIND = DP) :: wq(nmodes)
!! Phonon frequency on the fine grid
REAL(KIND = DP) :: inv_wq(nmodes)
!! $frac{1}{2\omega_{q\nu}}$ defined for efficiency reasons
REAL(KIND = DP) :: g2_tmp(nmodes)
!! If the phonon frequency is too small discart g
REAL(KIND = DP) :: wgq(nmodes)
!! Bose occupation factor $n_{q\nu}(T)$
REAL(KIND = DP) :: eta2(nbndfst, nmodes, nktotf)
!! Temporary array to store the current smearing eta
REAL(KIND = DP) :: inv_eta(nbndfst, nmodes, nktotf)
!! Temporary array to store the inverse of the eta for speed purposes
REAL(KIND = DP), ALLOCATABLE :: xkf_all(:,
!! Collect k-point coordinate from all pools in parallel case
REAL(KIND = DP), ALLOCATABLE :: etf_all(:,
!! Collect eigenenergies from all pools in parallel case
REAL(KIND = DP) :: vkk(3, nbndfst)
!! Electronic velocity $v_{nk}$
REAL(KIND = DP) :: vkq(3, nbndfst)
!! Electronic velocity $v_{nk+q}$
REAL(KIND = DP) :: coskkq(nbndfst, nbndfst)
!! $$(v_k \cdot v_{k+q}) / |v_k|^2$$
REAL(KIND = DP) :: vkkk(3, nbndfst, nkf)
!! vkk ik
REAL(KIND = DP) :: vkqq(3, nbndfst, nbndfst, nkf)
!! vkk ik
REAL(KIND = DP) :: vkksquare(nbndfst, nkf)
!! $$DDOT(vkk,vkk)
REAL(KIND = DP) :: dfnk
!! Derivative Fermi distribution $$-df_{nk}/dE_{nk}$$
REAL(KIND = DP) :: dfnkk(nbndfst, nkf)
!! Derivative Fermi distribution $$-df_{nk}/dE_{nk}$$
REAL(KIND = DP) :: etemp
!! Temperature in Ry (this includes division by kb)
!
! SP: Define the inverse so that we can efficiently multiply instead of dividing
inv_degaussw = one /degaussw
! To avoid if branching in the loop
inv_eta(:, :, = zero
IF (adapt_smearing) THEN
DO ik = 1, nkf
DO ibnd = 1, nbndfst
DO imode = 1, nmodes
inv_eta(ibnd, imode, ik) = one / (DSQRT(two) * eta(imode, ibnd, ik))
eta2(ibnd, imode, ik) = DSQRT(two) * eta(imode, ibnd, ik)
ENDDO
ENDDO
ENDDO
ELSE
DO ik = 1, nkf
DO ibnd = 1, nbndfst
DO imode = 1, nmodes
inv_eta(ibnd, imode, ik) = inv_degaussw
eta2(ibnd, imode, ik) = degaussw
ENDDO
ENDDO
ENDDO
ENDIF
DO itemp = 1, nstemp ! loop over temperatures
inv_eptemp = one / gtemp(itemp)
etemp = gtemp(itemp)
!
! Now pre-treat phonon modes for efficiency
! Treat phonon frequency and Bose occupation
wq(:) = zero
DO imode = 1, nmodes
wq(imode) = wf(imode, iq)
IF (wq(imode) > eps_acustic) THEN
g2_tmp(imode) = one
wgq(imode) = wgauss(-wq(imode) * inv_eptemp, -99)
wgq(imode) = wgq(imode) / (one - two * wgq(imode))
inv_wq(imode) = one / (two * wq(imode))
ELSE
g2_tmp(imode) = zero
wgq(imode) = zero
inv_wq(imode) = zero
ENDIF
ENDDO
!
IF (iqq == 1) THEN
!
WRITE(stdout, '(/5x, a)') REPEAT('=', 67)
WRITE(stdout, '(5x, "Electron (Imaginary) Self-Energy in the Migdal Approximation")')
WRITE(stdout, '(5x, a/)') REPEAT('=', 67)
!
IF (fsthick < 1.d3) WRITE(stdout, '(/5x, a, f10.6, a)' ) 'Fermi Surface thickness = ', fsthick * ryd2ev, ' eV'
WRITE(stdout, '(/5x, a, f10.6, a)') 'Golden Rule strictly enforced with T = ', gtemp(itemp) * ryd2ev, ' eV'
!
ENDIF
!
! Fermi level
!
IF (efermi_read) THEN
ef0 = fermi_energy
ELSE
ef0 = efnew
ENDIF
!
IF ((iqq == 1) .AND. .NOT. adapt_smearing) THEN
WRITE(stdout, 100) degaussw * ryd2ev, ngaussw
WRITE(stdout, '(a)') ' '
ENDIF
!
IF (restart) THEN
! Make everythin 0 except the range of k-points we are working on
sigmar_all(:, 1:lower_bnd - 1, = zero
sigmar_all(:, lower_bnd + nkf:nktotf, = zero
sigmai_all(:, 1:lower_bnd - 1, = zero
sigmai_all(:, lower_bnd + nkf:nktotf, = zero
zi_all(:, 1:lower_bnd - 1, = zero
zi_all(:, lower_bnd + nkf:nktotf, = zero
!
ENDIF
!
! In the case of a restart do not add the first step
IF (first_cycle .and. itemp == nstemp) THEN
first_cycle = .FALSE.
ELSE
!
! loop over all k points of the fine mesh
!
fermicount = 0
DO ik = 1, nkf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
! here we must have ef, not ef0, to be consistent with ephwann_shuffle
! (but in this case they are the same)
coskkq = zero
! this is coskkq = (vk dot vkq) / |vk||vkq|
! In principle the only coskkq contributing to lambda_tr are both near the
! Fermi surface and the magnitudes will not differ greatly between vk and vkq
! we may implement the approximation to the angle between k and k+q
! vectors also listed in Grimvall
!
IF ((MINVAL(ABS(etf(:, ikk) - ef)) < fsthick) .AND. &
(MINVAL(ABS(etf(:, ikq) - ef)) < fsthick)) THEN
!
fermicount = fermicount + 1
DO imode = 1, nmodes
DO ibnd = 1, nbndfst
!
! the energy of the electron at k (relative to Ef)
ekk = etf(ibndmin - 1 + ibnd, ikk) - ef0
!
eta_tmp = eta2(ibnd, imode, ik)
sq_eta_tmp = eta_tmp**two
inv_eta_tmp = inv_eta(ibnd, imode, ik)
dfnk = w0gauss(ekk / etemp, -99) / etemp
!dfnkk(ibnd, ikk) = w0gauss(ekk / etemp, -99) / etemp
!
DO jbnd = 1, nbndfst
!
!
! in below, vkk, vkq, coskkq are added
! for calculating the momentum relaxed el-ph self energy calculation
! v_(k,i) = 1/m <ki|p|ki> = 2 * dmef (:, i,i,k)
! 1/m = 2 in Rydberg atomic units
!
vkk(:, ibnd) = REAL(vmef(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk))
vkkk(:, ibnd, ikk) = REAL(vmef(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk))
vkq(:, jbnd) = REAL(vmef(:, ibndmin - 1 + jbnd, ibndmin - 1 + jbnd, ikq))
vkqq(:, ibnd, jbnd, ikq) = REAL(vmef(:, ibndmin - 1 + jbnd, ibndmin - 1 + jbnd, ikq))
IF (ABS(vkk(1, ibnd)**two + vkk(2, ibnd)**two + vkk(3, ibnd)**two) > eps4) &
coskkq(ibnd, jbnd) = DDOT(3, vkk(:, ibnd), 1, vkq(:, jbnd), 1) / &
((DSQRT(vkk(1, ibnd)**two + vkk(2, ibnd)**two + vkk(3, ibnd)**two)) * (DSQRT(vkq(1, ibnd)**two + vkq(2, ibnd)**two + vkq(3, ibnd)**two)))
IF (ABS(vkk(1, ibnd)**two + vkk(2, ibnd)**two + vkk(3, ibnd)**two) > eps4) THEN
vkksquare(ibnd, ikk) = DDOT(3, vkk(:, ibnd), 1, vkk(:, ibnd), 1)
ELSE
vkksquare(ibnd, ikk) = 0
END IF
!
! the energy of the electron at k+q (relative to Ef)
ekq = etf(ibndmin - 1 + jbnd, ikq) - ef0
! the Fermi occupation at k+q
wgkq = wgauss(-ekq * inv_eptemp, -99)
!
! here we take into account the zero-point DSQRT(hbar/2M\omega)
! with hbar = 1 and M already contained in the eigenmodes
! g2 is Ry^2, wkf must already account for the spin factor
!
! SP: Shortrange is disabled for efficiency reasons
!IF (shortrange .AND. ( ABS(xqf(1, iq))> eps8 .OR. ABS(xqf(2, iq))> eps8 &
! .OR. ABS(xqf(3, iq))> eps8 )) THEN
! ! SP: The abs has to be removed. Indeed the epf17 can be a pure imaginary
! ! number, in which case its square will be a negative number.
! g2 = REAL((epf17(jbnd, ibnd, imode, ik)**two) * inv_wq(imode) * g2_tmp(imode))
!ELSE
! g2 = (ABS(epf17(jbnd, ibnd, imode, ik))**two) * inv_wq(imode) * g2_tmp(imode)
!ENDIF
!
g2 = (ABS(epf17(jbnd, ibnd, imode, ik))**two) * inv_wq(imode) * g2_tmp(imode)
!
! There is a sign error for wq in Eq. 9 of Comp. Phys. Comm. 181, 2140 (2010). - RM
! The sign was corrected according to Eq. (7.282) page 489 from Mahan's book
! (Many-Particle Physics, 3rd edition)
!
fact1 = wgkq + wgq(imode)
fact2 = one - wgkq + wgq(imode)
etmp1 = ekk - (ekq - wq(imode))
etmp2 = ekk - (ekq + wq(imode))
sq_etmp1 = etmp1 * etmp1
sq_etmp2 = etmp2 * etmp2
!
weight = wqf(iq) * REAL(fact1 / (etmp1 - ci * eta_tmp) + fact2 / (etmp2 - ci * eta_tmp)) * (1.0d0 - coskkq(ibnd, jbnd))
!
! \Re\Sigma [Eq. 3 in Comput. Phys. Commun. 209, 116 (2016)]
sigmar_all(ibnd, ik + lower_bnd - 1, itemp) = sigmar_all(ibnd, ik + lower_bnd - 1, itemp) + g2 * weight
!
! Logical implementation
! weight = wqf(iq) * aimag( &
! ( ( wgkq + wgq ) / ( ekk - ( ekq - wq ) - ci * degaussw ) + &
! ( one - wgkq + wgq ) / ( ekk - ( ekq + wq ) - ci * degaussw ) ) )
!
! Delta implementation
w0g1 = w0gauss(etmp1 * inv_eta_tmp, 0) * inv_eta_tmp
w0g2 = w0gauss(etmp2 * inv_eta_tmp, 0) * inv_eta_tmp
!
weight = pi * wqf(iq) * (fact1 * w0g1 + fact2 * w0g2) * (1.0d0 - coskkq(ibnd, jbnd))
!
! \Im\Sigma using delta approx. [Eq. 8 in Comput. Phys. Commun. 209, 116 (2016)]
sigmai_all(ibnd, ik + lower_bnd - 1, itemp) = sigmai_all(ibnd, ik + lower_bnd - 1, itemp) + g2 * weight
!
! Mode-resolved
IF (iverbosity == 3) THEN
sigmai_mode(ibnd, imode, ik + lower_bnd - 1, itemp) = sigmai_mode(ibnd, imode, ik + lower_bnd - 1, itemp) + &
g2 * weight
ENDIF
!
! Z FACTOR: -\frac{\partial\Re\Sigma}{\partial\omega}
!
weight = wqf(iq) * &
(fact1 * (sq_etmp1 - sq_eta_tmp) / (sq_etmp1 + sq_eta_tmp)**two + &
fact2 * (sq_etmp2 - sq_eta_tmp) / (sq_etmp2 + sq_eta_tmp)**two)
!
zi_all(ibnd, ik + lower_bnd - 1, itemp) = zi_all(ibnd, ik + lower_bnd - 1, itemp) + g2 * weight
!
ENDDO !jbnd
!need to add integral (dfnk * vkk^2 * elselftau) / integral (dfnk * vkk^2)
ENDDO !ibnd
ENDDO !imode
ENDIF ! endif fsthick
ENDDO ! end loop on k
!
! Creation of a restart point
IF (restart) THEN
IF (MOD(iqq, restart_step) == 0 .and. itemp == nstemp) THEN
WRITE(stdout, '(5x, a, i10)' ) 'Creation of a restart point at ', iqq
CALL mp_sum(sigmar_all, inter_pool_comm)
CALL mp_sum(sigmai_all, inter_pool_comm)
CALL mp_sum(zi_all, inter_pool_comm)
CALL mp_sum(fermicount, inter_pool_comm)
CALL mp_barrier(inter_pool_comm)
CALL selfen_el_write(iqq, totq, nktotf, sigmar_all, sigmai_all, zi_all)
ENDIF
ENDIF
ENDIF ! in case of restart, do not do the first one
ENDDO ! itemp
!
! The k points are distributed among pools: here we collect them
!
IF (iqq == totq) THEN
!
ALLOCATE(xkf_all(3, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_elec_q', 'Error allocating xkf_all', 1)
ALLOCATE(etf_all(nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_elec_q', 'Error allocating etf_all', 1)
xkf_all(:, = zero
etf_all(:, = zero
!
#if defined(__MPI)
!
! note that poolgather2 works with the doubled grid (k and k+q)
!
CALL poolgather2(3, nkqtotf, nkqf, xkf, xkf_all)
CALL poolgather2(nbndsub, nkqtotf, nkqf, etf, etf_all)
CALL mp_sum(sigmar_all, inter_pool_comm)
CALL mp_sum(sigmai_all, inter_pool_comm)
IF (iverbosity == 3) CALL mp_sum(sigmai_mode, inter_pool_comm)
CALL mp_sum(zi_all, inter_pool_comm)
CALL mp_sum(fermicount, inter_pool_comm)
CALL mp_barrier(inter_pool_comm)
!
#else
!
xkf_all = xkf
etf_all = etf
!
#endif
!
DO itemp = 1, nstemp
! Average over degenerate eigenstates:
WRITE(stdout, '(5x,"Average over degenerate eigenstates is performed")')
WRITE(stdout, '(5x, a, f8.3, a)') "Temperature: ", gtemp(itemp) * ryd2ev / kelvin2eV, "K"
!
DO ik = 1, nktotf
ikk = 2 * ik - 1
ikq = ikk + 1
!
DO ibnd = 1, nbndfst
ekk = etf_all(ibndmin - 1 + ibnd, ikk)
n = 0
tmp1 = zero
tmp2 = zero
tmp3 = zero
DO jbnd = 1, nbndfst
ekk1 = etf_all(ibndmin - 1 + jbnd, ikk)
IF (ABS(ekk1 - ekk) < eps6) THEN
n = n + 1
tmp1 = tmp1 + sigmar_all(jbnd, ik, itemp)
tmp2 = tmp2 + sigmai_all(jbnd, ik, itemp)
tmp3 = tmp3 + zi_all(jbnd, ik, itemp)
ENDIF
!
ENDDO ! jbnd
sigmar_tmp(ibnd) = tmp1 / FLOAT(n)
sigmai_tmp(ibnd) = tmp2 / FLOAT(n)
zi_tmp(ibnd) = tmp3 / FLOAT(n)
!
ENDDO ! ibnd
sigmar_all(:, ik, itemp) = sigmar_tmp(:)
sigmai_all(:, ik, itemp) = sigmai_tmp(:)
zi_all(:, ik, itemp) = zi_tmp(:)
!
ENDDO ! nktotf
!
! Output electron SE here after looping over all q-points (with their contributions
! summed in sigmar_all, etc.)
!
WRITE(stdout, '(5x,"WARNING: only the eigenstates within the Fermi window are meaningful")')
!
IF (mpime == ionode_id) THEN
! Write to file
WRITE(tp, "(f8.3)") gtemp(itemp) * ryd2ev / kelvin2eV
fileselfen = 'linewidth.elself.' // trim(adjustl(tp)) // 'K'
OPEN(UNIT = linewidth_elself, FILE = fileselfen)
WRITE(linewidth_elself, '(a)') '# Electron linewidth = 2*Im(Sigma) (meV)'
IF (iverbosity == 3) THEN
WRITE(linewidth_elself, '(a)') '# ik ibnd E(ibnd) imode Im(Sigma)(meV)'
ELSE
WRITE(linewidth_elself, '(a)') '# ik ibnd E(ibnd) Im(Sigma)(meV)'
ENDIF
!
DO ik = 1, nktotf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
WRITE(stdout, '(/5x, "ik = ", i7," coord.: ", 3f12.7)') ik, xkf_all(:, ikk)
WRITE(stdout, '(5x, a)') REPEAT('-', 67)
!
DO ibnd = 1, nbndfst
!
! note that ekk does not depend on q
ekk = etf_all(ibndmin - 1 + ibnd, ikk) - ef0
!
! calculate Z = 1 / ( 1 -\frac{\partial\Sigma}{\partial\omega} )
zi_all(ibnd, ik, itemp) = one / (one + zi_all(ibnd, ik, itemp))
!
WRITE(stdout, 102) ibndmin - 1 + ibnd, ryd2ev * ekk, ryd2mev * sigmar_all(ibnd, ik, itemp), &
ryd2mev * sigmai_all(ibnd,ik, itemp), zi_all(ibnd, ik, itemp), &
one / zi_all(ibnd, ik, itemp) - one
IF (iverbosity == 3) THEN
DO imode = 1, nmodes
WRITE(linewidth_elself, '(i9, 2x)', ADVANCE = 'no') ik
WRITE(linewidth_elself, '(i9, 2x)', ADVANCE = 'no') ibndmin - 1 + ibnd
WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') ryd2ev * ekk
WRITE(linewidth_elself, '(i9, 2x)', ADVANCE = 'no') imode
WRITE(linewidth_elself, '(E22.14, 2x)') ryd2mev * sigmai_mode(ibnd, imode, ik, itemp)
ENDDO
ELSE
WRITE(linewidth_elself, '(i9, 2x)', ADVANCE = 'no') ik
WRITE(linewidth_elself, '(i9, 2x)', ADVANCE = 'no') ibndmin - 1 + ibnd
WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') ryd2ev * ekk
!WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') ryd2ev * dfnkk(ibnd, ikk)
WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkkk(1, ibnd, ikk)
WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkkk(2, ibnd, ikk)
WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkkk(3, ibnd, ikk)
WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkksquare(ibnd, ikk)
!DO jbnd = 1, nbndfst
! WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkqq(1, ibnd, jbnd, ikq)
! WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkqq(2, ibnd, jbnd, ikq)
! WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkqq(3, ibnd, jbnd, ikq)
!ENDDO
WRITE(linewidth_elself, '(E22.14, 2x)') ryd2mev * sigmai_all(ibnd, ik, itemp)
ENDIF
!
ENDDO
WRITE(stdout, '(5x, a/)') REPEAT('-', 67)
!
ENDDO
CLOSE(linewidth_elself)
ENDIF
!
DO ibnd = 1, nbndfst
DO ik = 1, nktotf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
! note that ekk does not depend on q
ekk = etf_all(ibndmin - 1 + ibnd, ikk) - ef0
!
! calculate Z = 1 / (1 - \frac{\partial\Sigma}{\partial\omega})
!zi_all(ibnd,ik) = one / (one + zi_all(ibnd,ik))
!
WRITE(stdout, '(2i9, 5f12.4)') ik, ibndmin - 1 + ibnd, ryd2ev * ekk, ryd2mev * sigmar_all(ibnd, ik, itemp), &
ryd2mev * sigmai_all(ibnd, ik, itemp), zi_all(ibnd, ik, itemp), &
one / zi_all(ibnd, ik, itemp) - one
!
ENDDO
!
WRITE(stdout, '(a)') ' '
!
ENDDO
!
ENDDO ! itemp
DEALLOCATE(xkf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_elec_q', 'Error deallocating xkf_all', 1)
DEALLOCATE(etf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_elec_q', 'Error deallocating etf_all', 1)
!
ENDIF
!
100 FORMAT(5x, 'Gaussian Broadening: ', f10.6, ' eV, ngauss=', i4)
102 FORMAT(5x, 'E( ', i3, ' )=', f9.4, ' eV Re[Sigma]=', f15.6, ' meV Im[Sigma]=', &
f15.6, ' meV Z=', f15.6, ' lam=', f15.6)
!
RETURN
!
!-----------------------------------------------------------------------
END SUBROUTINE selfen_elec_q
!-----------------------------------------------------------------------
Thank you again for your kind help.
Sincerely,
Thank you for your comment.
At first, I thought -3xx was -E3xx.
What I doubt is that when vkk^2 is a combination of -E0x and -3xx of vkk, there are parts that become -3xx.
And I want to calculate momentum relaxed electron-phonon relaxation time value.
There are equations of (4) and (5) in https://doi.org/10.1038/s41567-021-01341-w
So, I'm trying to modifying epw code for this calculation.
In below, i added vkk, vkq, vkksquare and dfnk. but dfnk doesn't work.
Finally, the value you want to calculate is {integral{(vkk^2*tau)dk}/{integral(vkk^2)dk} as shown in Equation (5).
Could you help me put the formula into the internal code?
!
! Copyright (C) 2016-2019 Samuel Ponce', Roxana Margine, Feliciano Giustino
! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino
! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino
!
! This file is distributed under the terms of the GNU General Public
! License. See the file `LICENSE' in the root directory of the
! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
!
!----------------------------------------------------------------------
MODULE selfen
!----------------------------------------------------------------------
!!
!! This module contains the various self-energy routines
!!
IMPLICIT NONE
!
CONTAINS
!
!-----------------------------------------------------------------------
SUBROUTINE selfen_elec_q(iqq, iq, totq, first_cycle)
!-----------------------------------------------------------------------
!!
!! Compute the imaginary part of the electron self energy due to electron-
!! phonon interaction in the Migdal approximation. This corresponds to
!! the electron linewidth (half width). The phonon frequency is taken into
!! account in the energy selection rule.
!!
!! Use matrix elements, electronic eigenvalues and phonon frequencies
!! from ep-wannier interpolation
!!
!! This routines computes the contribution from phonon iq to all k-points
!! The outer loop in ephwann_shuffle.f90 will loop over all iq points
!! The contribution from each iq is summed at the end of this subroutine
!! for iqq=totq to recover the per-ik electron self energy
!!
!! RM 24/02/2014
!! Redefined the size of sigmar_all, sigmai_all, and zi_all within the fermi windwow
!!
!-----------------------------------------------------------------------
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE io_var, ONLY : linewidth_elself
USE modes, ONLY : nmodes
USE epwcom, ONLY : nstemp, nbndsub, shortrange, fsthick, ngaussw, degaussw, &
eps_acustic, efermi_read, fermi_energy, restart, restart_step
USE pwcom, ONLY : ef
USE elph2, ONLY : etf, ibndmin, ibndmax, nkqf, xqf, eta, nbndfst, &
nkf, epf17, wf, wqf, xkf, nkqtotf, adapt_smearing, &
sigmar_all, sigmai_all, sigmai_mode, zi_all, efnew, &
nktotf, lower_bnd, gtemp, &
vmef, dmef
USE control_flags, ONLY : iverbosity
USE constants_epw, ONLY : kelvin2eV, ryd2mev, ryd2ev, one, two, zero, ci, eps4, eps6, eps8, hbar
USE constants, ONLY : pi
USE mp, ONLY : mp_barrier, mp_sum
USE mp_global, ONLY : inter_pool_comm
USE mp_world, ONLY : mpime
USE io_global, ONLY : ionode_id
USE io_selfen, ONLY : selfen_el_write
USE poolgathering, ONLY : poolgather2
!
IMPLICIT NONE
!
CHARACTER(LEN = 20) :: tp
!! string for temperatures
CHARACTER(LEN = 256) :: fileselfen
!! file name for self energy
LOGICAL, INTENT(inout) :: first_cycle
!! Use to determine weather this is the first cycle after restart
INTEGER, INTENT(in) :: iqq
!! Q-point index from selecq.fmt window
INTEGER, INTENT(in) :: iq
!! Q-point index from full grid
INTEGER, INTENT(in) :: totq
!! Total number of q-points from the selecq.fmt grid.
!
! Local variables
INTEGER :: n
!! Integer for the degenerate average over eigenstates
INTEGER :: ik
!! Counter on the k-point index
INTEGER :: ikk
!! k-point index
INTEGER :: ikq
!! q-point index
INTEGER :: ibnd
!! Counter on bands at k
INTEGER :: jbnd
!! Counter on bands at k+q
INTEGER :: imode
!! Counter on mode
INTEGER :: fermicount
!! Number of states on the Fermi surface
INTEGER :: itemp
!! Counter on temperatures
INTEGER :: ierr
!! Error status
!
REAL(KIND = DP) :: g2
!! Electron-phonon matrix elements squared in Ry^2
REAL(KIND = DP) :: ef0
!! Fermi energy level
REAL(KIND = DP) :: ekk
!! Eigen energy at k on the fine grid relative to the Fermi level
REAL(KIND = DP) :: ekk1
!! Temporary variable to the eigenenergies at k for the degenerate average
REAL(KIND = DP) :: ekq
!! Eigen energy at k+q on the fine grid relative to the Fermi level
REAL(KIND = DP) :: etmp1
!! Temporary variable to store etmp1 = ekk - (ekq - wq)
REAL(KIND = DP) :: etmp2
!! Temporary variable to strore etmp2 = ekk - (ekq + wq)
REAL(KIND = DP) :: sq_etmp1
!! Temporary variable to store etmp1^2
REAL(KIND = DP) :: sq_etmp2
!! Temporary variable to store etmp2^2
REAL(KIND = DP) :: wgkq
!! Fermi-Dirac occupation factor $f_{mk+q}(T)$
REAL(KIND = DP) :: fact1
!! Temporary variable to store $f_{mk+q}(T) + n_{q\nu}(T)$
REAL(KIND = DP) :: fact2
!! Temporary variable to store $1 - f_{mk+q}(T) + n_{q\nu}(T)$
REAL(KIND = DP) :: weight
!! Self-energy factor
!!$$ N_q \Re(\frac{f_{mk+q}(T) + n_{q\nu}(T)}{\varepsilon_{nk} - \varepsilon_{mk+q} + \omega_{q\nu} - i\delta}) $$
!!$$ + N_q \Re(\frac{1 - f_{mk+q}(T) + n_{q\nu}(T)}{\varepsilon_{nk} - \varepsilon_{mk+q} - \omega_{q\nu} - i\delta}) $$
REAL(KIND = DP) :: w0g1
!! Dirac delta at k for the imaginary part of $\Sigma$
REAL(KIND = DP) :: w0g2
!! Dirac delta at k+q for the imaginary part of $\Sigma$
REAL(KIND = DP) :: inv_degaussw
!! Inverse of degaussw define for efficiency reasons
REAL(KIND = DP) :: inv_eptemp
!! Inverse of temperature define for efficiency reasons
REAL(KIND = DP) :: eta_tmp
!! Temporary variable eta2
REAL(KIND = DP) :: sq_eta_tmp
!! Temporary eta2^2
REAL(KIND = DP) :: inv_eta_tmp
!! Temporary varialbe inv_eta
REAL(KIND = DP) :: tmp1
!! Temporary variable to store real part of Sigma for the degenerate average
REAL(KIND = DP) :: tmp2
!! Temporary variable to store imag part of Sigma for the degenerate average
REAL(KIND = DP) :: tmp3
!! Temporary variable to store Z for the degenerate average
REAL(KIND = DP) :: DDOT
!! Dot product function
REAL(KIND = DP) :: sigmar_tmp(nbndfst)
!! Temporary array to store the real-part of Sigma
REAL(KIND = DP) :: sigmai_tmp(nbndfst)
!! Temporary array to store the imag-part of Sigma
REAL(KIND = DP) :: zi_tmp(nbndfst)
!! Temporary array to store the Z
REAL(KIND = DP), EXTERNAL :: wgauss
!! Fermi-Dirac distribution function (when -99)
REAL(KIND = DP), EXTERNAL :: w0gauss
!! This function computes the derivative of the Fermi-Dirac function
!! It is therefore an approximation for a delta function
REAL(KIND = DP) :: wq(nmodes)
!! Phonon frequency on the fine grid
REAL(KIND = DP) :: inv_wq(nmodes)
!! $frac{1}{2\omega_{q\nu}}$ defined for efficiency reasons
REAL(KIND = DP) :: g2_tmp(nmodes)
!! If the phonon frequency is too small discart g
REAL(KIND = DP) :: wgq(nmodes)
!! Bose occupation factor $n_{q\nu}(T)$
REAL(KIND = DP) :: eta2(nbndfst, nmodes, nktotf)
!! Temporary array to store the current smearing eta
REAL(KIND = DP) :: inv_eta(nbndfst, nmodes, nktotf)
!! Temporary array to store the inverse of the eta for speed purposes
REAL(KIND = DP), ALLOCATABLE :: xkf_all(:,
!! Collect k-point coordinate from all pools in parallel case
REAL(KIND = DP), ALLOCATABLE :: etf_all(:,
!! Collect eigenenergies from all pools in parallel case
REAL(KIND = DP) :: vkk(3, nbndfst)
!! Electronic velocity $v_{nk}$
REAL(KIND = DP) :: vkq(3, nbndfst)
!! Electronic velocity $v_{nk+q}$
REAL(KIND = DP) :: coskkq(nbndfst, nbndfst)
!! $$(v_k \cdot v_{k+q}) / |v_k|^2$$
REAL(KIND = DP) :: vkkk(3, nbndfst, nkf)
!! vkk ik
REAL(KIND = DP) :: vkqq(3, nbndfst, nbndfst, nkf)
!! vkk ik
REAL(KIND = DP) :: vkksquare(nbndfst, nkf)
!! $$DDOT(vkk,vkk)
REAL(KIND = DP) :: dfnk
!! Derivative Fermi distribution $$-df_{nk}/dE_{nk}$$
REAL(KIND = DP) :: dfnkk(nbndfst, nkf)
!! Derivative Fermi distribution $$-df_{nk}/dE_{nk}$$
REAL(KIND = DP) :: etemp
!! Temperature in Ry (this includes division by kb)
!
! SP: Define the inverse so that we can efficiently multiply instead of dividing
inv_degaussw = one /degaussw
! To avoid if branching in the loop
inv_eta(:, :,
= zeroIF (adapt_smearing) THEN
DO ik = 1, nkf
DO ibnd = 1, nbndfst
DO imode = 1, nmodes
inv_eta(ibnd, imode, ik) = one / (DSQRT(two) * eta(imode, ibnd, ik))
eta2(ibnd, imode, ik) = DSQRT(two) * eta(imode, ibnd, ik)
ENDDO
ENDDO
ENDDO
ELSE
DO ik = 1, nkf
DO ibnd = 1, nbndfst
DO imode = 1, nmodes
inv_eta(ibnd, imode, ik) = inv_degaussw
eta2(ibnd, imode, ik) = degaussw
ENDDO
ENDDO
ENDDO
ENDIF
DO itemp = 1, nstemp ! loop over temperatures
inv_eptemp = one / gtemp(itemp)
etemp = gtemp(itemp)
!
! Now pre-treat phonon modes for efficiency
! Treat phonon frequency and Bose occupation
wq(:) = zero
DO imode = 1, nmodes
wq(imode) = wf(imode, iq)
IF (wq(imode) > eps_acustic) THEN
g2_tmp(imode) = one
wgq(imode) = wgauss(-wq(imode) * inv_eptemp, -99)
wgq(imode) = wgq(imode) / (one - two * wgq(imode))
inv_wq(imode) = one / (two * wq(imode))
ELSE
g2_tmp(imode) = zero
wgq(imode) = zero
inv_wq(imode) = zero
ENDIF
ENDDO
!
IF (iqq == 1) THEN
!
WRITE(stdout, '(/5x, a)') REPEAT('=', 67)
WRITE(stdout, '(5x, "Electron (Imaginary) Self-Energy in the Migdal Approximation")')
WRITE(stdout, '(5x, a/)') REPEAT('=', 67)
!
IF (fsthick < 1.d3) WRITE(stdout, '(/5x, a, f10.6, a)' ) 'Fermi Surface thickness = ', fsthick * ryd2ev, ' eV'
WRITE(stdout, '(/5x, a, f10.6, a)') 'Golden Rule strictly enforced with T = ', gtemp(itemp) * ryd2ev, ' eV'
!
ENDIF
!
! Fermi level
!
IF (efermi_read) THEN
ef0 = fermi_energy
ELSE
ef0 = efnew
ENDIF
!
IF ((iqq == 1) .AND. .NOT. adapt_smearing) THEN
WRITE(stdout, 100) degaussw * ryd2ev, ngaussw
WRITE(stdout, '(a)') ' '
ENDIF
!
IF (restart) THEN
! Make everythin 0 except the range of k-points we are working on
sigmar_all(:, 1:lower_bnd - 1,
= zerosigmar_all(:, lower_bnd + nkf:nktotf,
= zerosigmai_all(:, 1:lower_bnd - 1,
= zerosigmai_all(:, lower_bnd + nkf:nktotf,
= zerozi_all(:, 1:lower_bnd - 1,
= zerozi_all(:, lower_bnd + nkf:nktotf,
= zero!
ENDIF
!
! In the case of a restart do not add the first step
IF (first_cycle .and. itemp == nstemp) THEN
first_cycle = .FALSE.
ELSE
!
! loop over all k points of the fine mesh
!
fermicount = 0
DO ik = 1, nkf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
! here we must have ef, not ef0, to be consistent with ephwann_shuffle
! (but in this case they are the same)
coskkq = zero
! this is coskkq = (vk dot vkq) / |vk||vkq|
! In principle the only coskkq contributing to lambda_tr are both near the
! Fermi surface and the magnitudes will not differ greatly between vk and vkq
! we may implement the approximation to the angle between k and k+q
! vectors also listed in Grimvall
!
IF ((MINVAL(ABS(etf(:, ikk) - ef)) < fsthick) .AND. &
(MINVAL(ABS(etf(:, ikq) - ef)) < fsthick)) THEN
!
fermicount = fermicount + 1
DO imode = 1, nmodes
DO ibnd = 1, nbndfst
!
! the energy of the electron at k (relative to Ef)
ekk = etf(ibndmin - 1 + ibnd, ikk) - ef0
!
eta_tmp = eta2(ibnd, imode, ik)
sq_eta_tmp = eta_tmp**two
inv_eta_tmp = inv_eta(ibnd, imode, ik)
dfnk = w0gauss(ekk / etemp, -99) / etemp
!dfnkk(ibnd, ikk) = w0gauss(ekk / etemp, -99) / etemp
!
DO jbnd = 1, nbndfst
!
!
! in below, vkk, vkq, coskkq are added
! for calculating the momentum relaxed el-ph self energy calculation
! v_(k,i) = 1/m <ki|p|ki> = 2 * dmef (:, i,i,k)
! 1/m = 2 in Rydberg atomic units
!
vkk(:, ibnd) = REAL(vmef(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk))
vkkk(:, ibnd, ikk) = REAL(vmef(:, ibndmin - 1 + ibnd, ibndmin - 1 + ibnd, ikk))
vkq(:, jbnd) = REAL(vmef(:, ibndmin - 1 + jbnd, ibndmin - 1 + jbnd, ikq))
vkqq(:, ibnd, jbnd, ikq) = REAL(vmef(:, ibndmin - 1 + jbnd, ibndmin - 1 + jbnd, ikq))
IF (ABS(vkk(1, ibnd)**two + vkk(2, ibnd)**two + vkk(3, ibnd)**two) > eps4) &
coskkq(ibnd, jbnd) = DDOT(3, vkk(:, ibnd), 1, vkq(:, jbnd), 1) / &
((DSQRT(vkk(1, ibnd)**two + vkk(2, ibnd)**two + vkk(3, ibnd)**two)) * (DSQRT(vkq(1, ibnd)**two + vkq(2, ibnd)**two + vkq(3, ibnd)**two)))
IF (ABS(vkk(1, ibnd)**two + vkk(2, ibnd)**two + vkk(3, ibnd)**two) > eps4) THEN
vkksquare(ibnd, ikk) = DDOT(3, vkk(:, ibnd), 1, vkk(:, ibnd), 1)
ELSE
vkksquare(ibnd, ikk) = 0
END IF
!
! the energy of the electron at k+q (relative to Ef)
ekq = etf(ibndmin - 1 + jbnd, ikq) - ef0
! the Fermi occupation at k+q
wgkq = wgauss(-ekq * inv_eptemp, -99)
!
! here we take into account the zero-point DSQRT(hbar/2M\omega)
! with hbar = 1 and M already contained in the eigenmodes
! g2 is Ry^2, wkf must already account for the spin factor
!
! SP: Shortrange is disabled for efficiency reasons
!IF (shortrange .AND. ( ABS(xqf(1, iq))> eps8 .OR. ABS(xqf(2, iq))> eps8 &
! .OR. ABS(xqf(3, iq))> eps8 )) THEN
! ! SP: The abs has to be removed. Indeed the epf17 can be a pure imaginary
! ! number, in which case its square will be a negative number.
! g2 = REAL((epf17(jbnd, ibnd, imode, ik)**two) * inv_wq(imode) * g2_tmp(imode))
!ELSE
! g2 = (ABS(epf17(jbnd, ibnd, imode, ik))**two) * inv_wq(imode) * g2_tmp(imode)
!ENDIF
!
g2 = (ABS(epf17(jbnd, ibnd, imode, ik))**two) * inv_wq(imode) * g2_tmp(imode)
!
! There is a sign error for wq in Eq. 9 of Comp. Phys. Comm. 181, 2140 (2010). - RM
! The sign was corrected according to Eq. (7.282) page 489 from Mahan's book
! (Many-Particle Physics, 3rd edition)
!
fact1 = wgkq + wgq(imode)
fact2 = one - wgkq + wgq(imode)
etmp1 = ekk - (ekq - wq(imode))
etmp2 = ekk - (ekq + wq(imode))
sq_etmp1 = etmp1 * etmp1
sq_etmp2 = etmp2 * etmp2
!
weight = wqf(iq) * REAL(fact1 / (etmp1 - ci * eta_tmp) + fact2 / (etmp2 - ci * eta_tmp)) * (1.0d0 - coskkq(ibnd, jbnd))
!
! \Re\Sigma [Eq. 3 in Comput. Phys. Commun. 209, 116 (2016)]
sigmar_all(ibnd, ik + lower_bnd - 1, itemp) = sigmar_all(ibnd, ik + lower_bnd - 1, itemp) + g2 * weight
!
! Logical implementation
! weight = wqf(iq) * aimag( &
! ( ( wgkq + wgq ) / ( ekk - ( ekq - wq ) - ci * degaussw ) + &
! ( one - wgkq + wgq ) / ( ekk - ( ekq + wq ) - ci * degaussw ) ) )
!
! Delta implementation
w0g1 = w0gauss(etmp1 * inv_eta_tmp, 0) * inv_eta_tmp
w0g2 = w0gauss(etmp2 * inv_eta_tmp, 0) * inv_eta_tmp
!
weight = pi * wqf(iq) * (fact1 * w0g1 + fact2 * w0g2) * (1.0d0 - coskkq(ibnd, jbnd))
!
! \Im\Sigma using delta approx. [Eq. 8 in Comput. Phys. Commun. 209, 116 (2016)]
sigmai_all(ibnd, ik + lower_bnd - 1, itemp) = sigmai_all(ibnd, ik + lower_bnd - 1, itemp) + g2 * weight
!
! Mode-resolved
IF (iverbosity == 3) THEN
sigmai_mode(ibnd, imode, ik + lower_bnd - 1, itemp) = sigmai_mode(ibnd, imode, ik + lower_bnd - 1, itemp) + &
g2 * weight
ENDIF
!
! Z FACTOR: -\frac{\partial\Re\Sigma}{\partial\omega}
!
weight = wqf(iq) * &
(fact1 * (sq_etmp1 - sq_eta_tmp) / (sq_etmp1 + sq_eta_tmp)**two + &
fact2 * (sq_etmp2 - sq_eta_tmp) / (sq_etmp2 + sq_eta_tmp)**two)
!
zi_all(ibnd, ik + lower_bnd - 1, itemp) = zi_all(ibnd, ik + lower_bnd - 1, itemp) + g2 * weight
!
ENDDO !jbnd
!need to add integral (dfnk * vkk^2 * elselftau) / integral (dfnk * vkk^2)
ENDDO !ibnd
ENDDO !imode
ENDIF ! endif fsthick
ENDDO ! end loop on k
!
! Creation of a restart point
IF (restart) THEN
IF (MOD(iqq, restart_step) == 0 .and. itemp == nstemp) THEN
WRITE(stdout, '(5x, a, i10)' ) 'Creation of a restart point at ', iqq
CALL mp_sum(sigmar_all, inter_pool_comm)
CALL mp_sum(sigmai_all, inter_pool_comm)
CALL mp_sum(zi_all, inter_pool_comm)
CALL mp_sum(fermicount, inter_pool_comm)
CALL mp_barrier(inter_pool_comm)
CALL selfen_el_write(iqq, totq, nktotf, sigmar_all, sigmai_all, zi_all)
ENDIF
ENDIF
ENDIF ! in case of restart, do not do the first one
ENDDO ! itemp
!
! The k points are distributed among pools: here we collect them
!
IF (iqq == totq) THEN
!
ALLOCATE(xkf_all(3, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_elec_q', 'Error allocating xkf_all', 1)
ALLOCATE(etf_all(nbndsub, nkqtotf), STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_elec_q', 'Error allocating etf_all', 1)
xkf_all(:,
= zeroetf_all(:,
= zero!
#if defined(__MPI)
!
! note that poolgather2 works with the doubled grid (k and k+q)
!
CALL poolgather2(3, nkqtotf, nkqf, xkf, xkf_all)
CALL poolgather2(nbndsub, nkqtotf, nkqf, etf, etf_all)
CALL mp_sum(sigmar_all, inter_pool_comm)
CALL mp_sum(sigmai_all, inter_pool_comm)
IF (iverbosity == 3) CALL mp_sum(sigmai_mode, inter_pool_comm)
CALL mp_sum(zi_all, inter_pool_comm)
CALL mp_sum(fermicount, inter_pool_comm)
CALL mp_barrier(inter_pool_comm)
!
#else
!
xkf_all = xkf
etf_all = etf
!
#endif
!
DO itemp = 1, nstemp
! Average over degenerate eigenstates:
WRITE(stdout, '(5x,"Average over degenerate eigenstates is performed")')
WRITE(stdout, '(5x, a, f8.3, a)') "Temperature: ", gtemp(itemp) * ryd2ev / kelvin2eV, "K"
!
DO ik = 1, nktotf
ikk = 2 * ik - 1
ikq = ikk + 1
!
DO ibnd = 1, nbndfst
ekk = etf_all(ibndmin - 1 + ibnd, ikk)
n = 0
tmp1 = zero
tmp2 = zero
tmp3 = zero
DO jbnd = 1, nbndfst
ekk1 = etf_all(ibndmin - 1 + jbnd, ikk)
IF (ABS(ekk1 - ekk) < eps6) THEN
n = n + 1
tmp1 = tmp1 + sigmar_all(jbnd, ik, itemp)
tmp2 = tmp2 + sigmai_all(jbnd, ik, itemp)
tmp3 = tmp3 + zi_all(jbnd, ik, itemp)
ENDIF
!
ENDDO ! jbnd
sigmar_tmp(ibnd) = tmp1 / FLOAT(n)
sigmai_tmp(ibnd) = tmp2 / FLOAT(n)
zi_tmp(ibnd) = tmp3 / FLOAT(n)
!
ENDDO ! ibnd
sigmar_all(:, ik, itemp) = sigmar_tmp(:)
sigmai_all(:, ik, itemp) = sigmai_tmp(:)
zi_all(:, ik, itemp) = zi_tmp(:)
!
ENDDO ! nktotf
!
! Output electron SE here after looping over all q-points (with their contributions
! summed in sigmar_all, etc.)
!
WRITE(stdout, '(5x,"WARNING: only the eigenstates within the Fermi window are meaningful")')
!
IF (mpime == ionode_id) THEN
! Write to file
WRITE(tp, "(f8.3)") gtemp(itemp) * ryd2ev / kelvin2eV
fileselfen = 'linewidth.elself.' // trim(adjustl(tp)) // 'K'
OPEN(UNIT = linewidth_elself, FILE = fileselfen)
WRITE(linewidth_elself, '(a)') '# Electron linewidth = 2*Im(Sigma) (meV)'
IF (iverbosity == 3) THEN
WRITE(linewidth_elself, '(a)') '# ik ibnd E(ibnd) imode Im(Sigma)(meV)'
ELSE
WRITE(linewidth_elself, '(a)') '# ik ibnd E(ibnd) Im(Sigma)(meV)'
ENDIF
!
DO ik = 1, nktotf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
WRITE(stdout, '(/5x, "ik = ", i7," coord.: ", 3f12.7)') ik, xkf_all(:, ikk)
WRITE(stdout, '(5x, a)') REPEAT('-', 67)
!
DO ibnd = 1, nbndfst
!
! note that ekk does not depend on q
ekk = etf_all(ibndmin - 1 + ibnd, ikk) - ef0
!
! calculate Z = 1 / ( 1 -\frac{\partial\Sigma}{\partial\omega} )
zi_all(ibnd, ik, itemp) = one / (one + zi_all(ibnd, ik, itemp))
!
WRITE(stdout, 102) ibndmin - 1 + ibnd, ryd2ev * ekk, ryd2mev * sigmar_all(ibnd, ik, itemp), &
ryd2mev * sigmai_all(ibnd,ik, itemp), zi_all(ibnd, ik, itemp), &
one / zi_all(ibnd, ik, itemp) - one
IF (iverbosity == 3) THEN
DO imode = 1, nmodes
WRITE(linewidth_elself, '(i9, 2x)', ADVANCE = 'no') ik
WRITE(linewidth_elself, '(i9, 2x)', ADVANCE = 'no') ibndmin - 1 + ibnd
WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') ryd2ev * ekk
WRITE(linewidth_elself, '(i9, 2x)', ADVANCE = 'no') imode
WRITE(linewidth_elself, '(E22.14, 2x)') ryd2mev * sigmai_mode(ibnd, imode, ik, itemp)
ENDDO
ELSE
WRITE(linewidth_elself, '(i9, 2x)', ADVANCE = 'no') ik
WRITE(linewidth_elself, '(i9, 2x)', ADVANCE = 'no') ibndmin - 1 + ibnd
WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') ryd2ev * ekk
!WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') ryd2ev * dfnkk(ibnd, ikk)
WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkkk(1, ibnd, ikk)
WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkkk(2, ibnd, ikk)
WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkkk(3, ibnd, ikk)
WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkksquare(ibnd, ikk)
!DO jbnd = 1, nbndfst
! WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkqq(1, ibnd, jbnd, ikq)
! WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkqq(2, ibnd, jbnd, ikq)
! WRITE(linewidth_elself, '(E22.14, 2x)', ADVANCE = 'no') vkqq(3, ibnd, jbnd, ikq)
!ENDDO
WRITE(linewidth_elself, '(E22.14, 2x)') ryd2mev * sigmai_all(ibnd, ik, itemp)
ENDIF
!
ENDDO
WRITE(stdout, '(5x, a/)') REPEAT('-', 67)
!
ENDDO
CLOSE(linewidth_elself)
ENDIF
!
DO ibnd = 1, nbndfst
DO ik = 1, nktotf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
! note that ekk does not depend on q
ekk = etf_all(ibndmin - 1 + ibnd, ikk) - ef0
!
! calculate Z = 1 / (1 - \frac{\partial\Sigma}{\partial\omega})
!zi_all(ibnd,ik) = one / (one + zi_all(ibnd,ik))
!
WRITE(stdout, '(2i9, 5f12.4)') ik, ibndmin - 1 + ibnd, ryd2ev * ekk, ryd2mev * sigmar_all(ibnd, ik, itemp), &
ryd2mev * sigmai_all(ibnd, ik, itemp), zi_all(ibnd, ik, itemp), &
one / zi_all(ibnd, ik, itemp) - one
!
ENDDO
!
WRITE(stdout, '(a)') ' '
!
ENDDO
!
ENDDO ! itemp
DEALLOCATE(xkf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_elec_q', 'Error deallocating xkf_all', 1)
DEALLOCATE(etf_all, STAT = ierr)
IF (ierr /= 0) CALL errore('selfen_elec_q', 'Error deallocating etf_all', 1)
!
ENDIF
!
100 FORMAT(5x, 'Gaussian Broadening: ', f10.6, ' eV, ngauss=', i4)
102 FORMAT(5x, 'E( ', i3, ' )=', f9.4, ' eV Re[Sigma]=', f15.6, ' meV Im[Sigma]=', &
f15.6, ' meV Z=', f15.6, ' lam=', f15.6)
!
RETURN
!
!-----------------------------------------------------------------------
END SUBROUTINE selfen_elec_q
!-----------------------------------------------------------------------
Thank you again for your kind help.
Sincerely,
Re: Hello. I have a question for vmef, dmef function, and I want to implement in electron self energy calculation.
Dear darjeeling:
In addition, I think that instead of the subroutine of selfen_elec_q, you had better start from the subroutines in the transport modules, for instance, io_transport.f90, transport.f90, transport_iter.f90, etc.
Other developers would add more details.
Sincerely,
H. Lee
As I said, I think that -3xx means E-3xx, that is, 10^(-3xx).At first, I thought -3xx was -E3xx.
What I doubt is that when vkk^2 is a combination of -E0x and -3xx of vkk, there are parts that become -3xx.
In addition, I think that instead of the subroutine of selfen_elec_q, you had better start from the subroutines in the transport modules, for instance, io_transport.f90, transport.f90, transport_iter.f90, etc.
Other developers would add more details.
Sincerely,
H. Lee