[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