!----------------------------------------------------------------------- subroutine elphsum3 !----------------------------------------------------------------------- ! ! Sum over BZ of the electron-phonon matrix elements el_ph_mat ! and symmetrize over the star of q ! ! Written by Feliciano Giustino for Wannier interpolation. ! Based on the original routine elphsum.f90. See endnote for ! rationale of coefficients and conversion factors. ! ! The summation on k is performed only for those states with energy ! less than fsthick (Fermi Surface thickness) from the FS. ! #include "f_defs.h" USE ions_base, ONLY : nat, tau, ityp USE io_global, ONLY : stdout use pwcom USE kinds, only : DP use phcom use el_phon USE control_flags, ONLY : iswitch, modenum, noinv #ifdef __PARA use para #endif implicit none ! real(kind=DP), parameter :: eps = 20.d0 / 13.6058d0 / 8065.5d0, & pibytwo = 3.1415926 / 2.d0, ryd2ev = 13.6058, & ryd2ghz = 3.289828d6, two = 2.d0, zero = 0.d0 integer :: ik, ikk, ikq, ibnd, jbnd, ipert, jpert, nu, mu, & vu, pu, iuelph, ios, isig, i real(kind=DP) :: weight, w0g1, w0g2, w0gauss, wgauss,dosef, & dos_ef, ef1, phase_space, lambda, gamma, degaussw_ , tmp, & ekk, ekq complex(kind=DP) :: el_ph_sum (3*nat,3*nat), gamma_mat(3*nat,3*nat), & el_ph_mat_f ( nbndsub, nbndsub, 3 * nat) complex(kind = 8), parameter :: czero = (0.d0, 0.d0) ! ! local variables ! integer :: nq, isq (48), imq, na, nb, nt, nta, ntb, imode0, jmode0, irr, jrr ! degeneracy of the star of q ! index of q in the star of a given sym.op. ! index of -q in the star of q (0 if not present) ! counter on atoms ! counter on atomic type ! counter on modes ! counter on representation real(kind=DP) :: sxq (3, 48) ! list of vectors in the star of q ! ! the following are to avoid the mess by star_q.f90 ! integer :: nsym_ , s_ (3, 3, 48), invs_ (48), irt_ (48, nat) real(kind=DP) :: rtau_ (3, 48, nat) ! integer :: nrec, fermicount logical :: tfermikk, tfermikq external dos_ef ! write(stdout, '(/5x,a)' ) 'Entering elphsum3: e-p matrix for the star of q' ! if ( elinterp ) then write(stdout, '(/5x,a)' ) 'Electron interpolation is ON' else write(stdout, '(/5x,a)' ) 'Electron interpolation is OFF' ! ! use the non-interpolated values ! nksf = nks nksqf = nksq ! allocate ( etf ( nbndsub, nksf ), wkf (nksf) ) ! etf = et wkf = wk ! endif ! if ( fsthick.lt.1.d3 ) & write(stdout, '(/5x,a,f10.6,a)' ) & 'Fermi Surface thickness = ', fsthick, ' Ry' write(stdout, '(/5x,a)' ) 'Using Allen approximation' ! if (filelph.ne.' ') then #ifdef __PARA ! parallel case: only first node writes if (me.ne.1) then iuelph = 0 else #endif iuelph = 4 open (unit = iuelph, file = filelph, status = 'unknown', err = & 100, iostat = ios) 100 call errore ('elphon', 'opening file'//filelph, abs (ios) ) rewind (iuelph) write (iuelph, '(3f15.8,f8.3,2i8)') xq, degaussw, ngaussw, 3 * nat write (iuelph, '(6e14.6)') (w2 (nu) , nu = 1, nmodes) #ifdef __PARA end if #endif else iuelph = 0 endif ! ! preliminary check on the nesting function against ! degauss/10 and degauss*10 @ FG ! do isig = 1,3 ! ! just the powers of 10... degaussw_ = degaussw * 10.d0**float(isig-2) ! call efermig (etf, nbndsub, nksf, nelec, wkf, degaussw_ , ngaussw, ef1) dosef = dos_ef (ngaussw, degaussw_ , ef1, etf, wkf, nksf, nbndsub) phase_space = zero ! do ik = 1, nksqf ! if (lgamma) then ikk = ik ikq = ik else ikk = 2 * ik - 1 ikq = ikk + 1 endif ! do ibnd = 1, nbndsub ! ekk = etf (ibnd, ikk) - ef1 if ( abs(ekk) .lt. fsthick ) then ! do jbnd = 1, nbndsub ! ekq = etf (jbnd, ikq) - ef1 if ( abs(ekq) .lt. fsthick ) then ! w0g1 = w0gauss( -ekk/degaussw_ , ngaussw)/degaussw_ w0g2 = w0gauss( -ekq/degaussw_ , ngaussw)/degaussw_ ! phase_space = phase_space + wkf (ikk) * w0g1 * w0g2 endif enddo ! endif enddo ! enddo ! #ifdef __PARA call poolreduce (1, phase_space) #endif ! write (stdout, 9007) degaussw_ , phase_space, phase_space / dosef / dosef ! enddo ! ! Fermi level and corresponding DOS ! ! Note that the weights of k+q points must be set to zero here call efermig (etf, nbndsub, nksf, nelec, wkf, degaussw, ngaussw, ef1) dosef = dos_ef (ngaussw, degaussw, ef1, etf, wkf, nksf, nbndsub) ! N(Ef) in the equation for lambda is the DOS per spin dosef = dosef / two ! ! Sum over bands with gaussian weights ! el_ph_sum = (zero, zero) phase_space = zero ! fermicount = 0 ! do ik = 1, nksqf ! if (lgamma) then ikk = ik ikq = ik else ikk = 2 * ik - 1 ikq = ikk + 1 endif ! tfermikk = .false. tfermikq = .false. ! do ibnd = 1, nbndsub if ( abs(etf (ibnd, ikk) - ef) .lt. fsthick ) tfermikk = .true. if ( abs(etf (ibnd, ikq) - ef) .lt. fsthick ) tfermikq = .true. enddo ! if (tfermikk.and.tfermikq) then ! fermicount = fermicount + 1 ! if ( elinterp ) then ! ! read from file el_ph_mat_f(ibndmin:ibndmax,ibndmin:ibndmax,:) for this ik ! the unread elements are padded with zeros ! write(stdout,'(5x,a,i6)',advance='no') 'Reading el_ph_mat_f from file for ik =',ik ! el_ph_mat_f = czero do ipert = 1, 3 * nat nrec = (ipert - 1) * nksqf + ik call dasmio ( el_ph_mat_f(ibndmin:ibndmax,ibndmin:ibndmax,ipert), & ibndmax-ibndmin+1, lrepmatf, iunepmatf, nrec, -1) enddo ! write(stdout,*) ' ... done' ! else ! ! use noninterpolated matrix ! do ipert = 1, 3 * nat do ibnd = 1, nbnd do jbnd = 1, nbnd el_ph_mat_f(ibnd, jbnd, ipert) = el_ph_mat(ibnd, jbnd, ik, ipert) enddo enddo enddo ! endif ! do ibnd = 1, nbndsub ! ekk = etf (ibnd, ikk) - ef1 ! do jbnd = 1, nbndsub ! ekq = etf (jbnd, ikq) - ef1 ! w0g1 = w0gauss( -ekk/degaussw, ngaussw)/degaussw w0g2 = w0gauss( -ekq/degaussw, ngaussw)/degaussw ! weight = wkf (ikk) * w0g1 * w0g2 ! do jpert = 1, 3 * nat do ipert = 1, 3 * nat el_ph_sum (ipert, jpert) = el_ph_sum (ipert, jpert) & + weight * conjg (el_ph_mat_f (jbnd, ibnd, ipert)) & * el_ph_mat_f (jbnd, ibnd, jpert) enddo enddo ! phase_space = phase_space+weight ! enddo ! enddo ! ! endif tfermi endif ! enddo ! write( stdout, '(/14x,a,i5,a,i5)' ) & 'Number of (k,k+q) pairs on the Fermi surface: ',fermicount, ' out of ', nksqf ! #ifdef __PARA ! ! collect contributions from all pools (sum over k-points) ! call poolreduce (2 * 3 * nat * 3 * nat, el_ph_sum) call poolreduce (1, phase_space) ! #endif ! ! Here el_ph_sum is the matrix in the 'pattern' representation ! I need the symmetrized matrix in cartesian basis for all q ! in the star ! ! Symmetrizes the gamma matrix w.r.t. the small group of q. ! It is a second rank tensor obtained from the tensor product ! of two first rank tensors therefore it transforms as the ! dynamical matrix ! ! The following call also brings ep_ph_sum on the cartesian basis ! call symdyn_munu (el_ph_sum, u, xq, s, invs, rtau, irt, irgq, at, & bg, nsymq, nat, irotmq, minus_q) ! ! Generates the star of q ! call star_q (xq, at, bg, ibrav, symm_type, nat, tau, ityp, nr1, & nr2, nr3, nsym_ , s_ , invs_ , irt_ , rtau_ , nq, sxq, isq, imq, noinv, & modenum) ! ! Rotates and writes the gamma matrices of the star of q ! These are the true matrices that will be used in the interpolation ! call q2qstar_ph (el_ph_sum, at, bg, nat, nsym_ , s_ , invs_ , irt_ , rtau_ , & nq, sxq, isq, imq, iuelph) ! ! !------------------------------------------------------------------------ ! ! ... What follows is only for printing out information needed ! ... when noninterpolated calculations are done (summary and check) ! ! ... go from cartesian to the vibrational eigenmodes ! ... dyn(mu,nu) contains the eigenmodes divided by the masses ! ... i.e. the true eigendisplacements ! gamma_mat = zero ! do nu = 1, nmodes ! do pu = 1, 3 * nat ! do mu = 1, 3 * nat do vu = 1, 3 * nat ! gamma_mat (nu, pu) = gamma_mat (nu, pu) & + conjg (dyn (mu, nu) ) * el_ph_sum (mu, vu) * dyn (vu, pu) ! enddo enddo ! enddo ! enddo ! write (6, 9000) degaussw, ngaussw write (6, 9005) dosef, ef1 * ryd2ev write (6, 9006) phase_space if (iuelph.ne.0) then write (iuelph, *) Write (iuelph, 9000) degaussw, ngaussw write (iuelph, 9005) dosef, ef1 * ryd2ev endif ! gamma_mat = pibytwo * gamma_mat ! do nu = 1, nmodes gamma = real( gamma_mat (nu, nu) ) if (sqrt (abs (w2 (nu) ) ) .gt.eps) then lambda = gamma / pi / w2 (nu) / dosef else lambda = zero endif gamma = gamma * ryd2ev * 1.d3 ! gamma = gamma * ryd2ghz write (6, 9011) nu, lambda, gamma if (iuelph.ne.0) write (iuelph, 9011) nu, lambda, gamma ! enddo ! 9000 format(5x,'Gaussian Broadening: ',f7.3,' Ry, ngauss=',i4) 9005 format(5x,'DOS =',f10.6,' states/spin/Ry/Unit Cell at Ef=',f10.6,' eV') 9006 format(5x,'double delta at Ef =',e10.3) 9007 format(5x,'Gaussian Broadening: ',f10.6, & ' --> double delta at Ef =',e10.3, & ' , double delta / Nf^2 = ', e12.5) 9010 format(5x,'lambda(',i2,')=',f10.6,' gamma=',f10.6,' GHz') 9011 format(5x,'lambda( ',i3,' )=',e20.10,' gamma=',e20.10,' meV') ! if (iuelph.ne.0) close (unit = iuelph) return ! ! --------------------------------------------------------------------- ! eps = 20 cm^-1, in Ry ! ! Note for efermig call: ! Recalculate the Fermi energy Ef=ef1 and the DOS at Ef, dosef = N(Ef) ! for this gaussian broadening ! Note that the weights of k+q points must be set to zero for the ! following call to yield correct results ! ! ! el_ph_sum(mu,nu)=\sum_k\sum_{i,j}[ ! x ! x \delta(e_{k,i}-Ef) \delta(e_{k+q,j} ! ! gamma = \pi \sum_k\sum_{i,j} \delta(e_{k,i}-Ef) \delta(e_{k+q,j}-Ef) ! | \sum_mu z(mu,nu) |^2 ! where z(mu,nu) is the mu component of normal mode nu (z = dyn) ! gamma(nu) is the phonon linewidth of mode nu ! ! the factor 2 comes from the factor sqrt(hbar/2/M/omega) that appears ! in the definition of the electron-phonon matrix element g ! The sqrt(1/M) factor is actually hidden into the normal modes ! ! The factor N(Ef)^2 that appears in most formulations of el-ph interact ! is absent because we sum, not average, over the Fermi surface. ! The factor 2 is provided by the sum over spins (i.e. sum(w_k) = 2 @FG) ! ! lambda is the adimensional el-ph coupling for mode nu: ! lambda(nu)= gamma(nu)/(pi N(Ef) \omega_{q,nu}^2) ! 3.289828x10^6 is the conversion factor from Ry to GHz ! --------------------------------------------------------------------- ! end subroutine elphsum3