! ! 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 . ! !----------------------------------------------------------------------- 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 ! ! !----------------------------------------------------------------------- #include "f_defs.h" USE kinds, only : DP USE cell_base, ONLY : at, bg USE io_global, ONLY : stdout use phcom, only : lgamma, nmodes use epwcom, only : nbndsub, lrepmatf, iunepmatf, & fsthick, eptemp, ngaussw, degaussw, iuetf, & wmin, wmax, nw, nbndskip, ecutse, parallel_k, & parallel_q, epf_mem, etf_mem use pwcom, only : nelec, ef, isk use el_phon, only : etf, ibndmin, ibndmax, nksf, etfq, & epf17, wkf, nksqf, nxqf, wf, wqf, xkf, nkstotf use mp_global, only : me_pool, inter_pool_comm #ifdef __PARA use mp #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, & weight, wgauss, dosef, dos_ef, sigmar2(nbndsub, nksf), & sigmai2(nbndsub, nksf), zi2(nbndsub, nksf) logical :: already_skipped real(kind=DP), external :: efermig ! ! variables for collecting data from all pools in parallel case ! integer :: nksqtotf real(kind=DP), allocatable :: xkf_all(:,:) , etf_all(:,:), sigma_all (:,:), & sigmar2_all (:,:), sigmai2_all (:,:), zi2_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' already_skipped = .false. ! ! 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 IF (.not. already_skipped) then nelec = nelec - two * nbndskip already_skipped = .true. WRITE(stdout, '(/5x,"Skipping the first ",i4," bands:")' ) nbndskip WRITE(stdout, '(/5x,"The Fermi level will be determined with ",f9.5," electrons")' ) nelec ENDIF ENDIF ! ef0 = efermig(etf,nbndsub,nksf,nelec,wkf,degaussw,ngaussw,0,isk) 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 ! ! loop over all k points of the fine mesh ! sigmar2 = zero sigmai2 = zero zi2 = zero ! DO ik = 1, nksqf ! IF (lgamma) then ikk = ik ikq = ik ELSE ikk = 2 * ik - 1 ikq = ikk + 1 ENDIF ! ! ! loop over the q-points ! DO iq = 1, nxqf ! ! ! we read the hamiltonian eigenvalues (those at k+q depend on q!) ! IF (etf_mem) then etf (ibndmin:ibndmax, ikk) = etfq (ibndmin:ibndmax, ikk, iq) etf (ibndmin:ibndmax, ikq) = etfq (ibndmin:ibndmax, ikq, iq) ELSE 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) ENDIF ! ! ! 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 ! 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 ! IF (etf_mem) then epf(:,:) = epf17 ( ik, iq, :, :,imode) ELSE nrec = (iq-1) * nmodes * nksqf + (imode-1) * nksqf + ik CALL dasmio ( epf, ibndmax-ibndmin+1, lrepmatf, iunepmatf, nrec, -1) ENDIF ! ! ! DO ibnd = 1, ibndmax-ibndmin+1 ! ! the energy of the electron at k (relative to Ef) 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 ! 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 ) ) ) !@ if ( abs(ekq-ekk) .gt. ecutse ) weight = 0.d0 ! sigmar2(ibndmin-1+ibnd,ikk) = sigmar2(ibndmin-1+ibnd,ikk) + g2 * weight ! weight = wqf(iq) * aimag ( & ( ( wgkq + wgq ) / ( ekk - ( ekq - wq ) - ci * degaussw ) + & ( one - wgkq + wgq ) / ( ekk - ( ekq + wq ) - ci * degaussw ) ) ) ! weight = wqf(iq) * aimag ( one / ( ekk - ekq - ci * degaussw ) ) ! just a DOS !@ if ( abs(ekq-ekk) .gt. ecutse ) weight = 0.d0 ! sigmai2(ibndmin-1+ibnd,ikk) = sigmai2(ibndmin-1+ibnd,ikk) + g2 * weight ! Z FACTOR: -\frac{\partial\Re\Sigma}{\partial\omega} ! weight = wqf(iq) * & ( & ( wgkq + wgq ) * ( (ekk - ( ekq - wq ))**two - degaussw**two ) / & ( (ekk - ( ekq - wq ))**two + degaussw**two )**two & + ( one - wgkq + wgq ) * ( (ekk - ( ekq + wq ))**two - degaussw**two ) / & ( (ekk - ( ekq + wq ))**two + degaussw**two )**two & ) !@ if ( abs(ekq-ekk) .gt. ecutse ) weight = 0.d0 ! zi2(ibndmin-1+ibnd,ikk) = zi2(ibndmin-1+ibnd,ikk) + g2 * weight ! ENDDO !jbnd ENDDO !ibnd ! ENDDO !imodes ! ! endif tfermi ENDIF ! ENDDO ! iq's ! WRITE (6,*) ik WRITE (6,'(f10.4)') sigmar2*ryd2mev ! ENDDO ! end loop on k ! 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/2 ! even-odd for k,k+q ! allocate ( xkf_all ( 3, nkstotf ), & etf_all ( nbndsub, nkstotf ), & sigmar2_all ( nbndsub, nkstotf), & sigmai2_all ( nbndsub, nkstotf), & sigma_all ( nbndsub, nkstotf), & zi2_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, sigmar2, sigmar2_all) CALL poolgather ( nbndsub, nkstotf, nksf, sigmai2, sigmai2_all) CALL poolgather ( nbndsub, nkstotf, nksf, zi2, zi2_all) ! #else ! xkf_all = xkf etf_all = etf sigmar2_all = sigmar2 sigmai2_all = sigmai2 sigma_all = sigmar2 !???? zi2_all = zi2 ! #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_pool == 0) & #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') ! end subroutine selfen_elec