! !----------------------------------------------------------------------- subroutine selfen_elec !----------------------------------------------------------------------- ! ! 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 ! ! Feliciano Giustino, University of California at Berkeley, 2006 ! !----------------------------------------------------------------------- ! #include "f_defs.h" USE kinds, only : DP USE cell_base, ONLY : at, bg USE io_global, ONLY : stdout use phcom, only : nbndsub, lrepmatf, iunepmatf, lgamma, nmodes, & fsthick, eptemp, ngaussw, degaussw, iuetf, wmin, wmax, nw, nbndskip use pwcom, only : nelec, ef use el_phon #ifdef __PARA use mp, only : mp_barrier use para #endif implicit none ! real(kind=DP), parameter :: eps = 20.d0 / 13.6058d0 / 8065.5d0, & ryd2mev = 13605.8, one = 1.d0, ryd2ev = 13.6058, & two = 2.d0, zero = 0.d0, pi = 3.14159265358979 complex(kind = 8), parameter :: ci = (0.d0, 1.d0), cone = (1.d0, 0.d0) integer :: ik, ikk, ikq, ibnd, jbnd, imode, nrec, iq, fermicount complex(kind=DP) epf (ibndmax-ibndmin+1, ibndmax-ibndmin+1) real(kind=DP) :: g2, ekk, ekq, wq, ef0, wgq, wgkq, lambda, & weight, w0g1, w0g2, wgauss, dosef, dos_ef, sigma(nbndsub, nksf), tmp ! ! variables for collecting data from all pools in parallel case ! integer :: nksqtotf real(kind=DP), allocatable :: xkf_all(:,:) , etf_all(:,:), sigma_all (:,:) ! write(6,'(/5x,a)') repeat('=',67) write(6,'(5x,"Electron (Imaginary) Self-Energy in the Migdal Approximation")') write(6,'(5x,a/)') repeat('=',67) ! if ( fsthick.lt.1.d3 ) & write(stdout, '(/5x,a,f10.6,a)' ) & 'Fermi Surface thickness = ', fsthick, ' Ry' write(stdout, '(/5x,a,f10.6,a)' ) & 'Golden Rule strictly enforced with T = ',eptemp, ' Ry' ! ! Fermi level and corresponding DOS ! ! here we take into account that we may skip bands when we wannierize ! (spin-unpolarized) if (nbndskip.gt.0) then nelec = nelec - two * nbndskip write(stdout, '(/5x,"Skipping the first ",i4," bands:")' ) nbndskip write(stdout, '(/5x,"The Fermi level will be determined with ",f9.5," electrons")' ) nelec endif ! call efermig (etf, nbndsub, nksf, nelec, wkf, degaussw, ngaussw, ef0) dosef = dos_ef (ngaussw, degaussw, ef0, etf, wkf, nksf, nbndsub) ! N(Ef) in the equation for lambda is the DOS per spin dosef = dosef / two ! write (6, 100) degaussw, ngaussw write (6, 101) dosef, ef0 * ryd2ev ! ! make sure q point weights add up to 1 wqf = wqf/sum(wqf) ! ! loop over all q points of the fine mesh ! sigma = zero ! do iq = 1, nxqf ! fermicount = 0 ! do ik = 1, nksqf ! if (lgamma) then ikk = ik ikq = ik else ikk = 2 * ik - 1 ikq = ikk + 1 endif ! ! we read the hamiltonian eigenvalues (those at k+q depend on q!) ! nrec = (iq-1) * nksf + ikk call davcio ( etf (ibndmin:ibndmax, ikk), ibndmax-ibndmin+1, iuetf, nrec, - 1) nrec = (iq-1) * nksf + ikq call davcio ( etf (ibndmin:ibndmax, ikq), ibndmax-ibndmin+1, iuetf, nrec, - 1) ! ! here we must have ef, not ef0, to be consistent with ephwann_shuffle if ( ( minval ( abs(etf (:, ikk) - ef) ) .lt. fsthick ) .and. & ( minval ( abs(etf (:, ikq) - ef) ) .lt. fsthick ) ) then ! fermicount = fermicount + 1 ! do imode = 1, nmodes ! ! the phonon frequency and Bose occupation wq = wf (imode, iq) wgq = wgauss( -wq/eptemp, -99) wgq = wgq / ( one - two * wgq ) ! ! we read the e-p matrix ! nrec = (iq-1) * nmodes * nksqf + (imode-1) * nksqf + ik call dasmio ( epf, ibndmax-ibndmin+1, lrepmatf, iunepmatf, nrec, -1) ! do ibnd = 1, ibndmax-ibndmin+1 ! ! the energy of the electron at k ekk = etf (ibndmin-1+ibnd, ikk) - ef0 ! do jbnd = 1, ibndmax-ibndmin+1 ! ! the fermi occupation for k+q ekq = etf (ibndmin-1+jbnd, ikq) - ef0 wgkq = wgauss( -ekq/eptemp, -99) ! ! here we take into account the zero-point sqrt(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 ! ! @@@@ not sure if it is i,j or j,i g2 = abs(epf (jbnd, ibnd))**two / ( two * wq ) ! ! Note the +ci and -ci below, this comes from Eq. (28) of my Notes - FG ! weight = two * wqf(iq) * aimag ( & ( ( wgkq + wgq ) / ( ekk - ( ekq - wq ) - ci * degaussw ) + & ( one - wgkq + wgq ) / ( ekk - ( ekq + wq ) + ci * degaussw ) ) ) ! sigma(ibndmin-1+ibnd,ikk) = sigma(ibndmin-1+ibnd,ikk) + g2 * weight ! enddo enddo ! enddo ! ! endif tfermi endif ! enddo ! ! no collection from pools here, q points are not distributed ! ! end loop on q enddo ! write( stdout, '(/5x,a,i5,a,i5/)' ) & 'Number of (k,k+q) pairs on the Fermi surface: ',fermicount, ' out of ', nksqf ! ! The k points are distributed amond pools: here we collect them ! nksqtotf = nkstotf/kunit ! even-odd for k,k+q ! allocate ( xkf_all ( 3, nkstotf ), & etf_all ( nbndsub, nkstotf ), & sigma_all ( nbndsub, nkstotf) ) ! #ifdef __PARA ! ! note that poolgather works with the doubled grid (k and k+q) ! therefore we need to use the sigma array with both grids, even ! though one of them is useless. This should be fixed be modifying ! poolgather (it's a waste of memory). ! call poolgather ( 3, nkstotf, nksf, xkf, xkf_all ) call poolgather ( nbndsub, nkstotf, nksf, etf, etf_all ) call poolgather ( nbndsub, nkstotf, nksf, sigma, sigma_all) ! #else ! xkf_all = xkf etf_all = etf sigma_all = sigma ! #endif ! write(6,'(5x,"WARNING: only the eigenstates within the Fermi window are meaningful")') ! do ik = 1, nksqtotf ! if (lgamma) then ikk = ik ikq = ik else ikk = 2 * ik - 1 ikq = ikk + 1 endif ! write(stdout,'(/5x,"ik = ",i5," coord.: ", 3f9.5)') ik, xkf_all (:,ikk) write(stdout,'(5x,a)') repeat('-',67) do ibnd = ibndmin, ibndmax ! ! note that ekk does not depend on q (while etf(:,ikq) changes for every q) ekk = etf_all (ibnd, ikk) - ef0 ! write(stdout, 102) ibnd, ryd2ev * ekk, ryd2mev * sigma_all (ibnd,ikk) ! enddo write(stdout,'(5x,a/)') repeat('-',67) ! #ifdef __PARA if (me.eq.1) & #endif write(901,'(6(2x,f9.5))') (ryd2mev*sigma_all(ibnd,ikk),ibnd=ibndmin,ibndmax) ! enddo ! 100 format(5x,'Gaussian Broadening: ',f7.3,' Ry, ngauss=',i4) 101 format(5x,'DOS =',f10.6,' states/spin/Ry/Unit Cell at Ef=',f10.6,' eV') 102 format(5x,'E( ',i3,' )=',f9.3,' eV Im[Sigma]=',f9.3,' meV') ! return end subroutine selfen_elec