!
! 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 nesting_fn
!-----------------------------------------------------------------------
!
! 10/2009 Compute the nesting function. Depends only on electronic
! eigenvalues. Copied from selfen_phon
!
!-----------------------------------------------------------------------
#include "f_defs.h"
USE kinds, only : DP
USE cell_base, ONLY : at, bg
USE io_global, ONLY : stdout
use phcom, only : nmodes
use epwcom, only : nbndsub, lrepmatf, iunepmatf, fsthick, &
eptemp, ngaussw, degaussw, iuetf, &
wmin, wmax, nw, nbndskip, a2f, epf_mem, etf_mem, &
nsmear, delta_smear
use pwcom, only : nelec, ef, isk, nbnd
use el_phon, only : epf17, ibndmax, ibndmin, etf, &
etfq, wkf, xqf, wqf, nksf, nxqf, &
nksqf, wf, nkstotf, xkf, xqf
#ifdef __PARA
use mp, only : mp_barrier,mp_sum
use mp_global, only : me_pool,inter_pool_comm,my_pool_id
#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), czero=(0.d0,0.d0)
integer :: ik, ikk, ikq, ibnd, jbnd, nrec, iq, fermicount, ismear
real(kind=DP) :: ekk, ekq, ef0, wgkk, wgkq, &
weight, w0g1, w0g2, w0gauss, wgauss, dosef, dos_ef, gamma, &
degaussw0, eptemp0
!
real(kind=DP), external :: efermig
logical :: already_skipped
!
!
!
WRITE(6,'(/5x,a)') repeat('=',67)
WRITE(6,'(5x,"Nesting function using double delta approach")')
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(1), ' Ry'
already_skipped = .false.
! here we loop on smearing values
DO ismear = 1, nsmear
!
degaussw0 = (ismear-1)*delta_smear+degaussw
eptemp0 = (ismear-1)*delta_smear+eptemp(1)
!
!
! 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
! this is necessary for a loop on smearings, as each time through
! we do not want to keep subtracting electrons.
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
!
! Note that the weights of k+q points must be set to zero here
! no spin-polarized calculation here
ef0 = efermig(etf,nbndsub,nksf,nelec,wkf,degaussw0,ngaussw,0,isk)
DOsef = dos_ef (ngaussw, degaussw0, ef0, etf, wkf, nksf, nbndsub)
! N(Ef) in the equation for lambda is the DOS per spin
DOsef = dosef / two
!
WRITE (6, 100) degaussw0, ngaussw
WRITE (6, 101) dosef, ef0 * ryd2ev
!
! loop over all q points of the fine mesh (this is in the k-para case
! it should always be k-para for selfen_phon)
!
DO iq = 1, nxqf
!
CALL start_clock('nesting')
!
fermicount = 0
gamma = zero
!
DO ik = 1, nksqf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
! 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
!
fermicount = fermicount + 1
!
!
DO ibnd = 1, ibndmax-ibndmin+1
!
! the fermi occupation for k
ekk = etf (ibndmin-1+ibnd, ikk) - ef0
wgkk = wgauss( -ekk/eptemp0, -99)
!
DO jbnd = 1, ibndmax-ibndmin+1
!
! the fermi occupation for k+q
ekq = etf (ibndmin-1+jbnd, ikq) - ef0
wgkq = wgauss( -ekq/eptemp0, -99)
!
!
! = k-point weight * [f(E_k) - f(E_k+q)]/ [E_k+q - E_k -w_q +id]
! This is the imaginary part of the phonon self-energy, sans the matrix elements
!
! weight = wkf (ikk) * (wgkk - wgkq) * &
! aimag ( cone / ( ekq - ekk - ci * degaussw0 ) )
!
! the below expression is positive-definite, but also an approximation
! which neglects some fine features
!
w0g1 = w0gauss ( ekk / degaussw0, 0) / degaussw0
w0g2 = w0gauss ( ekq / degaussw0, 0) / degaussw0
weight = pi * wkf (ikk) * w0g1 * w0g2
!
gamma = gamma + weight
!
ENDDO ! jbnd
ENDDO ! ibnd
!
!
ENDIF ! endif fsthick
!
CALL stop_clock('nesting')
!
ENDDO ! loop on q
#ifdef __PARA
!
! collect contributions from all pools (sum over k-points)
! this finishes the integral over the BZ (k)
!
CALL mp_sum(gamma,inter_pool_comm)
CALL mp_sum(fermicount, inter_pool_comm)
!
#endif
!
WRITE(6,'(/5x,"iq = ",i5," coord.: ", 3f9.5, " wt: ", f9.5)') iq, xqf(:,iq) , wqf(iq)
WRITE(6,'(5x,a)') repeat('-',67)
!
WRITE(6, 102) gamma
WRITE(6,'(5x,a/)') repeat('-',67)
!
! test only
#ifdef __PARA
! if (me.eq.1) &
IF (me_pool == 0) &
#endif
!
!
WRITE( stdout, '(/5x,a,i8,a,i8/)' ) &
'Number of (k,k+q) pairs on the Fermi surface: ',fermicount, ' out of ', nkstotf/2
ENDDO
!
ENDDO !smears
! ! generate the Eliashberg spectral function
!
! if (a2f) call eliashberg_a2f( gamma_all, lambda_all, gamma_all_v)
!
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,' Nesting function (q)=',f9.3,' units??')
!
end subroutine nesting_fn
!-----------------------------------------------------------------------
subroutine nesting_fn_fly (iq )
!-----------------------------------------------------------------------
!
! compute the imaginary part of the phonon self energy due to electron-
! phonon interaction in the Migdal approximation. This corresponds to
! the phonon 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 routine is similar to the one above
! but it is only called from within ephwann_shuffle and calculates
! the selfenergy for one phonon at a time. Much smaller footprint on
! the disk
!
!-----------------------------------------------------------------------
#include "f_defs.h"
USE kinds, only : DP
USE cell_base, ONLY : at, bg
USE io_global, ONLY : stdout
use phcom, only : nmodes
use epwcom, only : nbndsub, lrepmatf, iunepmatf, fsthick, &
eptemp, ngaussw, degaussw, iuetf, &
wmin, wmax, nw, nbndskip, a2f, epf_mem, etf_mem, &
nsmear, delta_smear
use pwcom, only : nelec, ef, isk, nbnd
use el_phon, only : epf17, ibndmax, ibndmin, etf, &
etfq, wkf, xqf, wqf, nksf, nxqf, &
nksqf, wf, nkstotf, xkf, xqf
#ifdef __PARA
use mp, only : mp_barrier,mp_sum
use mp_global, only : me_pool,inter_pool_comm,my_pool_id
#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), czero=(0.d0,0.d0)
integer :: ik, ikk, ikq, ibnd, jbnd, nrec, iq, fermicount, ismear
real(kind=DP) :: ekk, ekq, ef0, wgkk, wgkq, &
weight, w0g1, w0g2, w0gauss, wgauss, dosef, dos_ef, gamma, &
degaussw0, eptemp0
!
real(kind=DP), external :: efermig
logical :: already_skipped
!
!
!
IF (iq.eq.1) then
WRITE(6,'(/5x,a)') repeat('=',67)
WRITE(6,'(5x,"Nesting Function in the double delta approx")')
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(1), ' Ry'
IF (nbndskip.gt.0) then
nelec = nelec - two * nbndskip
! this is necessary for a loop on smearings, as each time through
! we do not want to keep subtracting electrons.
ENDIF
ENDIF
! here we loop on smearing values
DO ismear = 1, nsmear
!
degaussw0 = (ismear-1)*delta_smear+degaussw
eptemp0 = (ismear-1)*delta_smear+eptemp(1)
!
!
! Fermi level and corresponding DOS
!
! here we take into account that we may skip bands when we wannierize
! (spin-unpolarized)
!
!
! Note that the weights of k+q points must be set to zero here
! no spin-polarized calculation here
ef0 = efermig(etf,nbndsub,nksf,nelec,wkf,degaussw0,ngaussw,0,isk)
DOsef = dos_ef (ngaussw, degaussw0, ef0, etf, wkf, nksf, nbndsub)
! N(Ef) in the equation for lambda is the DOS per spin
DOsef = dosef / two
!
IF (iq.eq.1) then
WRITE (6, 100) degaussw0, ngaussw
WRITE (6, 101) dosef, ef0 * ryd2ev
ENDIF
!
!
CALL start_clock('nesting')
!
fermicount = 0
!
DO ik = 1, nksqf
!
ikk = 2 * ik - 1
ikq = ikk + 1
!
! we read the hamiltonian eigenvalues (those at k+q depend on q!)
!
! when we see references to iq for file readins, it is always = 1 for on the fly calculations
IF (etf_mem) then
etf (ibndmin:ibndmax, ikk) = etfq (ibndmin:ibndmax, ikk, 1)
etf (ibndmin:ibndmax, ikq) = etfq (ibndmin:ibndmax, ikq, 1)
ELSE
nrec = (iq-1) * nksf + ikk
nrec = ikk
CALL davcio ( etf (ibndmin:ibndmax, ikk), ibndmax-ibndmin+1, iuetf, nrec, - 1)
nrec = (iq-1) * nksf + ikq
nrec = 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
!
fermicount = fermicount + 1
!
DO ibnd = 1, ibndmax-ibndmin+1
!
! the fermi occupation for k
ekk = etf (ibndmin-1+ibnd, ikk) - ef0
wgkk = wgauss( -ekk/eptemp0, -99)
!
DO jbnd = 1, ibndmax-ibndmin+1
!
! the fermi occupation for k+q
ekq = etf (ibndmin-1+jbnd, ikq) - ef0
wgkq = wgauss( -ekq/eptemp0, -99)
!
! = k-point weight * [f(E_k) - f(E_k+q)]/ [E_k+q - E_k -w_q +id]
! This is the imaginary part of the phonon self-energy, sans the matrix elements
!
! weight = wkf (ikk) * (wgkk - wgkq) * &
! aimag ( cone / ( ekq - ekk - ci * degaussw ) )
!
! the below expression is positive-definite, but also an approximation
! which neglects some fine features
!
w0g1 = w0gauss ( ekk / degaussw0, 0) / degaussw0
w0g2 = w0gauss ( ekq / degaussw0, 0) / degaussw0
weight = pi * wkf (ikk) * w0g1 * w0g2
!
gamma = gamma + weight
!
ENDDO ! jbnd
ENDDO ! ibnd
!
!
!
ENDIF ! endif fsthick
!
CALL stop_clock('nesting')
!
ENDDO ! loop on k
#ifdef __PARA
!
! collect contributions from all pools (sum over k-points)
! this finishes the integral over the BZ (k)
!
CALL mp_sum(gamma,inter_pool_comm)
CALL mp_sum(fermicount, inter_pool_comm)
!
#endif
!
WRITE(6,'(/5x,"iq = ",i5," coord.: ", 3f9.5, " wt: ", f9.5)') iq, xqf(:,iq) , wqf(iq)
WRITE(6,'(5x,a)') repeat('-',67)
!
WRITE(6, 102) gamma
WRITE(6,'(5x,a/)') repeat('-',67)
CALL flush(6)
!
WRITE( stdout, '(/5x,a,i8,a,i8/)' ) &
'Number of (k,k+q) pairs on the Fermi surface: ',fermicount, ' out of ', nkstotf/2
!
!
ENDDO !smears
!
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,' Nesting function (q)=',f9.3,' units??')
end subroutine nesting_fn_fly
!