! !----------------------------------------------------------------------- subroutine spectral_func !----------------------------------------------------------------------- ! ! compute the electron spectral function including the electron- ! phonon interaction in the Migdal approximation. ! ! We take the trace of the spectral function to simulate the photoemission ! intensity. I do not consider the c-axis average for the time being. ! The main approximation is constant dipole matrix element and diagonal ! selfenergy. The diagonality can be checked numerically. ! ! Use matrix elements, electronic eigenvalues and phonon frequencies ! from ep-wannier interpolation ! ! Feliciano Giustino, University of California at Berkeley, 2007 ! !----------------------------------------------------------------------- ! #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 :: iw, ik, ikk, ikq, ibnd, jbnd, imode, nrec, iq, fermicount complex(kind=DP) epf (ibndmax-ibndmin+1, ibndmax-ibndmin+1), weight real(kind=DP) :: g2, ekk, ekq, wq, ef0, wgq, wgkq, lambda, ww, dw, & w0g1, w0g2, wgauss, dosef, dos_ef, sigmar(nbndsub, nksf, nw), & sigmai(nbndsub, nksf, nw), tmp, a(nw, nksf) ! ! variables for collecting data from all pools in parallel case ! integer :: nksqtotf real(kind=DP), allocatable :: xkf_all(:,:) , etf_all(:,:), sigmar_all (:,:), & sigmai_all (:,:), a_all(:,:) ! write(6,'(/5x,a)') repeat('=',67) write(6,'(5x,"Electron Spectral Function in the Migdal Approximation")') write(6,'(5x,a/)') repeat('=',67) ! ! energy range and spacing for spectral function ! wmin = -fsthick wmax = fsthick dw = ( wmax - wmin ) / float (nw-1) ! 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 ! sigmar = zero sigmai = 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 ! ! the loop over the energy variable ! do iw = 1, nw ! ww = wmin + float (iw-1) * dw ! weight = two * wqf(iq) * & ( ( wgkq + wgq ) / ( ww - ( ekq - wq ) - ci * degaussw ) + & ( one - wgkq + wgq ) / ( ww - ( ekq + wq ) + ci * degaussw ) ) ! sigmar(ibndmin-1+ibnd,ikk,iw) = sigmar(ibndmin-1+ibnd,ikk,iw) + g2 * real ( weight ) sigmai(ibndmin-1+ibnd,ikk,iw) = sigmai(ibndmin-1+ibnd,ikk,iw) + g2 * aimag ( weight ) ! enddo ! enddo enddo ! enddo ! ! endif tfermi endif ! enddo ! ! no collection from pools here, q points are not distributed ! ! end loop on q enddo ! ! construct the trace of the spectral function (assume diagonal selfenrgy and constant ! matrix elements for dipole transitions) ! do ik = 1, nksqf ! if (lgamma) then ikk = ik ikq = ik else ikk = 2 * ik - 1 ikq = ikk + 1 endif ! do iw = 1, nw ! ww = wmin + float (iw-1) * dw a(iw,ikk) = zero ! do ibnd = 1, ibndmax-ibndmin+1 ! ! the energy of the electron at k ekk = etf (ibndmin-1+ibnd, ikk) - ef0 ! a(iw,ikk) = a(iw,ikk) + abs( sigmai(ibndmin-1+ibnd,ikk,iw) ) / pi / & ( ( ww - ekk - sigmar(ibndmin-1+ibnd,ikk,iw) )**two + ( sigmai(ibndmin-1+ibnd,ikk,iw) )**two ) ! enddo ! enddo ! 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 ), & a_all ( nw, 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 ( nw, nkstotf, nksf, a, a_all) ! #else ! xkf_all = xkf etf_all = etf a_all = a ! #endif ! 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 iw = 1, nw ww = wmin + float (iw-1) * dw write(stdout, 103) ik, ryd2ev * ww, a_all(iw,ikk) / ryd2ev #ifdef __PARA if (me.eq.1) & #endif write(903,'(2x,i7,2x,f10.5,2x,e12.5)') ik, ryd2ev * ww, a_all(iw,ikk) / ryd2ev enddo write(903,*) write(stdout,'(5x,a/)') repeat('-',67) ! ! 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') 103 format(5x,'ik = ',i7,' w = ',f10.5,' eV A(k,w) = ',e12.5,' eV^-1') ! return end subroutine spectral_func