[Pw_forum] Electron-phonon coupling

Eyvaz Isaev eyvaz_isaev at yahoo.com
Thu Apr 3 11:58:08 CEST 2003


Dear Paolo, 

Thanks a lot!!! 

Eyvaz.
--- Paolo Giannozzi <giannozz at nest.sns.it> wrote:
> Hi
> 
> > According to 7-th example,  the code is able to
> > calculate \lambda_{nk} (in the example k-point is
> X).
> > I am calculating it for a number of k-points in
> the
> > 1/48 part of the BZ for FCC lattice. The origin of
> the
> > k-points I am using is slightly shifted (so, it is
> > not k=(000)).  I supposed that later \lambda_{nk}
> > could be  integrated using the tetrahedron method.
> 
> you do not need any tetrahedron method: you just
> average over
> the irreducible BZ. You can use gaussian broadening
> to calculate
> alpha^2F(omega). I am attaching the code I used to
> sum over q
> (no warranty that it works) and some notes I wrote.
> 
> Paolo
> --
> Paolo Giannozzi             e-mail: 
> giannozz at nest.sns.it
> Scuola Normale Superiore    Phone:   +39/050509412
> Piazza dei Cavalieri 7      Fax:     +39/050509417,
> 050563513
> I-56126 Pisa, Italy         Office:  Lab. NEST, Via
> della Faggiola 19
> 
> > 
>       program elph
> 
>       implicit none
>       integer npk, nsigx, nmodex, nex
>       parameter(npk=200, nsigx=50, nmodex=100,
> nex=200)
>       integer nks, ios, iuelph, ngauss, ngauss1,
> ngaussq, nsig, nmodes
>       integer ik, ng, mu, nu, i
>       real*8 q(3,npk), wk(npk), degauss(nsigx),
> w2(nmodex),
>      +     dosef(nsigx), ef(nsigx),
> lambdaq(nmodex,nsigx),
>      +     lambda(nsigx), alpha2F(nex,nsigx), logavg
>       real*8 qread(3), dosef1, ef1, degauss1,
> gammaq, lambda2,
>      +     w0gauss, degaussq, emax, deltae, e,
> omega, sum
>       character*50 filelph
>       external w0gauss
> 
> c emax and degaussq in THz
>       read(5,*) emax, degaussq,ngaussq
>       deltae=emax/(nex-1)
>       read(5,*) nks
>       if (nks.gt.npk) call error('lambda','too many
> q-points',npk)
>       sum=0.d0
>       do ik=1,nks
>          read(5,*) q(1,ik), q(2,ik), q(3,ik), wk(ik)
>          sum = sum + wk(ik)
>       end do
>       do ik=1,nks
>          wk(ik)=wk(ik)/sum
>       end do
> 
>       iuelph=4
>       do ik=1,nks
>          read(5,'(a)') filelph
>         
>
open(unit=iuelph,file=filelph,status='old',iostat=ios)
>          read (iuelph,*) qread(1),qread(2),qread(3),
> nsig, nmodes
>          if ( (qread(1)-q(1,ik))**2 +
> (qread(2)-q(2,ik))**2 
>      +       +(qread(3)-q(3,ik))**2 .gt. 1.d-6) 
>      +      call error('lambda','inconsistent q
> read',ik)
>          if (nsig.le.0.or.nsig.gt.nsigx)
>      +        call error('lambda','wrong/too many
> gauss.broad.',nsigx)
>          if (nmodes.le.0.or.nmodes.gt.nmodex)
>      +        call error('lambda','wrong # or too
> many modes',nmodex)
>          if (ik.eq.1) then
>             do ng=1,nsig
>                lambda(ng)=0.d0
>                do i=1,nex                  
>                   alpha2F(i,ng)=0.d0
>                end do
>             end do
>          end if
> c read omega^2
>          read(iuelph,*) (w2(nu),nu=1,nmodes)
> c read data
>          do ng=1,nsig
>             read (iuelph,9000) degauss1, ngauss1
>             if (ik.eq.1) then
>                degauss(ng)=degauss1
>                ngauss =ngauss1
>             else
>                if
> (degauss(ng).ne.degauss1.or.ngauss.ne.ngauss1)
>      +          call error('lambda','inconsistent
> gauss.broad. read',ik)
>             end if
>             read (iuelph,9005) dosef1, ef1
>             if (ik.eq.1) then
>                dosef(ng)=dosef1
>                ef(ng)=ef1
>             else
>                if
> (dosef(ng).ne.dosef1.or.ef(ng).ne.ef1)
>      +          call error('lambda','inconsistent
> DOS(Ef) read',ik)
>             end if
>             do mu=1,nmodes
>                read (iuelph,9010) nu,
> lambdaq(mu,ng), gammaq
>                if (nu.ne.mu) call
> error('lambda','wrong mode read',mu)
> c sum over k-points
>                lambda(ng) = lambda(ng) +
> wk(ik)*lambdaq(mu,ng)
>                do i=1,nex
>                   e=(i-1)*deltae
> c 1 Ry = 3289.828 THz
>                   omega=sqrt(w2(mu))*3289.828
>                   alpha2F(i,ng) = alpha2F(i,ng)
>      +                 +
> wk(ik)*lambdaq(mu,ng)*omega*0.5d0
>      +                 
> *w0gauss((e-omega)/degaussq,ngaussq)/degaussq
>                end do
> c
>             end do
> c
>          end do
>          close(unit=iuelph)
> 
>       end do
> 
>      
> open(unit=iuelph,file='lambda.dat',status='unknown',
>      +     form='formatted')
>       write(iuelph,9014)
>       do ng=1,nsig
> c lambda2 is used as a check
> c logavg is the logarithmic average of omega used in
> McMillan's formula (?)
>          lambda2=0.d0
>          logavg =0.d0
>          do i=2,nex
>             e=(i-1)*deltae
>             lambda2=lambda2 + alpha2F(i,ng)/e
>             logavg =logavg + alpha2F(i,ng)*log(e)/e
>          end do
>          lambda2=lambda2*2.d0*deltae
>          logavg =logavg*2.d0 *deltae
> c 1 THz = 50 K
>          logavg=exp(logavg/lambda2)*50.d0
>          write(6,9015) lambda(ng), lambda2,
> logavg,dosef(ng),degauss(ng)
>          write(iuelph,9016) 
>      +        degauss(ng), lambda(ng), lambda2,
> logavg,dosef(ng)
>       end do
>       close(unit=iuelph)
> 
>      
>
open(unit=iuelph,file='alpha2F.dat',status='unknown',
>      +     form='formatted')
>       write(iuelph,9020) (degauss(ng),ng=1,nsig)
>       do i=1,nex
>          e=(i-1)*deltae
>          write(iuelph,9025)
> e,(alpha2F(i,ng),ng=1,nsig)
>       end do
>       close(unit=iuelph)
> 
>       stop
>  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')
>  9010 format(5x,'lambda(',i2,')=',f10.6,'  
> gamma=',f10.6,' GHz')
>  9014 format('# degauss   lambda    int alpha2F 
> <log w>     N(Ef)')
>  9015 format(5x,'lambda =',f9.6,' (',f10.6,')  <log
> w>=',f9.3,'K  ',
>      +          'N(Ef)=',f9.6,' at degauss=',f5.3)
>  9016 format(f7.3,2f12.6,f10.3,2f12.6)
>  9020 format('# E(THz)',10(f10.3))
>  9025 format(f8.4,10(f10.5))
> 
>       end
> 
> 

> ATTACHMENT part 3 application/x-tex name=e-ph.tex



__________________________________________________
Do you Yahoo!?
Yahoo! Tax Center - File online, calculators, forms, and more
http://tax.yahoo.com



More information about the Pw_forum mailing list