[Pw_forum] Electron-phonon coupling

[Paolo Giannozzi]

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

-------------- next part --------------

      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

-------------- next part --------------
A non-text attachment was scrubbed...
Name: e-ph.tex
Type: application/x-tex
Size: 2214 bytes
Desc: 
Url : /pipermail/attachments/20030402/c5390192/attachment.bin