m_dos.f90

Go to the documentation of this file.
00001 module m_dos 
00002   implicit none
00003 contains
00004   !
00005   subroutine calculate_dos(NTB,NTK,nkb1,nkb2,nkb3,electron_number,FermiEnergy,filename,b1,b2,b3,SK0,EIG,VKS)
00006     implicit none 
00007     integer,intent(in)::NTB,NTK,nkb1,nkb2,nkb3
00008     real(8),intent(in)::electron_number
00009     real(8),intent(in)::b1(3),b2(3),b3(3)
00010     real(8),intent(in)::SK0(3,NTK)
00011     real(8),intent(in)::EIG(NTB,NTK) 
00012     real(8),intent(inout)::FermiEnergy
00013     character(99),intent(in)::filename 
00014     !
00015     complex(8),intent(in)::VKS(NTB,NTB,NTK)
00016     real(8),allocatable::pdos(:,:)!pdos(ndosgrd,NTB) 
00017     real(8)::flg_pdos
00018     !
00019     real(8)::emax!=maxval(EIG)
00020     real(8)::emin!=minval(EIG)
00021     integer::ndosgrd!=int(2.0d0*diff/dlt)+1
00022     real(8),allocatable::dosgrd(:)!dosgrd(ndosgrd) 
00023     real(8),allocatable::dos(:)!dos(ndosgrd) 
00024     real(8),allocatable::efline(:)!efline(ndosgrd)   
00025     !
00026     real(8),parameter::au=27.21151d0
00027     !real(8),parameter::delt=0.005d0/au!Greens function delt in au 
00028     real(8),parameter::delt=0.01d0/au!Greens function delt in au 
00029     real(8),parameter::dmna=1.0d-3!Ttrhdrn parameter dmna in au 
00030     real(8),parameter::dmnr=1.0d-3!Ttrhdrn parameter dmnr in au 
00031     real(8),parameter::delw=2.0d0*delt!Grid width in au 
00032     ! 
00033     !dos-grid
00034     !
00035     call est_ndosgrd(NTB,NTK,EIG(1,1),delw,emin,emax,ndosgrd)  
00036     allocate(dosgrd(ndosgrd));dosgrd=0.0d0 
00037     call make_dosgrd(emin,delw,ndosgrd,dosgrd(1))  
00038     !
00039     !calc dos 
00040     !
00041     allocate(dos(ndosgrd));dos=0.0d0 
00042     call calc_dos(NTB,NTK,nkb1,nkb2,nkb3,ndosgrd,dosgrd(1),EIG(1,1),SK0(1,1),delt,dmnr,dmna,b1(1),b2(1),b3(1),dos(1)) 
00043     !
00044     !calc pdos 
00045     !
00046     flg_pdos=maxval(abs(VKS)) 
00047     if(flg_pdos/=0.0d0)then
00048       write(6,'(a40,f15.8,a)')'flg_pdos/=0.0:',flg_pdos,': PDOS calculate:' 
00049       allocate(pdos(ndosgrd,NTB));pdos=0.0d0 
00050       call calc_pdos(NTB,NTK,nkb1,nkb2,nkb3,ndosgrd,dosgrd(1),EIG(1,1),VKS(1,1,1),SK0(1,1),delt,dmnr,dmna,b1(1),b2(1),b3(1),pdos(1,1)) 
00051     endif 
00052     !
00053     !estimate ef 
00054     !
00055     if(electron_number/=0.0d0)then
00056      call est_ef(ndosgrd,delw,electron_number,dosgrd(1),dos(1),FermiEnergy) 
00057     endif 
00058     allocate(efline(ndosgrd));efline=0.0d0 
00059     call make_efline(ndosgrd,FermiEnergy,delw,dosgrd(1),dos(1),efline(1)) 
00060     !
00061     !wrt dat.dos 
00062     !
00063     call wrt_dos(filename,ndosgrd,dosgrd(1),dos(1),efline(1)) 
00064     !
00065     !wrt dat.pdos 
00066     !
00067     if(flg_pdos/=0.0d0)then
00068       call wrt_pdos(filename,NTB,ndosgrd,dosgrd(1),pdos(1,1),efline(1)) 
00069     endif 
00070     !
00071     return
00072   end subroutine calculate_dos 
00073   ! 
00074   subroutine est_ndosgrd(NTB,NTK,E_EIG,deltw,emin,emax,ndosgrd)  
00075     implicit none
00076     integer,intent(in)::NTB,NTK
00077     real(8),intent(in)::E_EIG(NTB,NTK)
00078     real(8),intent(in)::deltw 
00079     real(8),intent(out)::emax,emin
00080     integer,intent(out)::ndosgrd
00081     real(8)::diff 
00082     real(8),parameter::au=27.21151d0
00083     !
00084     emax=maxval(E_EIG)
00085     emin=minval(E_EIG)
00086     diff=emax-emin 
00087     !
00088     !define grid range
00089     !
00090     ndosgrd=int(2.0d0*diff/deltw)+1
00091     emax=emax+0.50d0*diff
00092     emin=emin-0.50d0*diff
00093     !
00094     write(6,*)
00095     write(6,'(a40)')'+++ m_dos: est_ndosgrd +++'
00096     write(6,'(a40)')'GRID DATA FOR DOS'
00097     write(6,'(a40,f15.8)')'emin(eV)',emin*au 
00098     write(6,'(a40,f15.8)')'emax(eV)',emax*au 
00099     write(6,'(a40,f15.8)')'deltw(eV)',deltw*au 
00100     write(6,'(a40,i15)')'ndosgrd',ndosgrd 
00101     write(6,*)
00102     return
00103   end subroutine
00104   !
00105   subroutine make_dosgrd(emin,deltw,ndosgrd,dosgrd)  
00106     implicit none
00107     real(8),intent(in)::emin,deltw 
00108     integer,intent(in)::ndosgrd
00109     real(8),intent(out)::dosgrd(ndosgrd) 
00110     integer::ie 
00111     !
00112     dosgrd=0.0d0 
00113     do ie=1,ndosgrd
00114      dosgrd(ie)=emin+dble(ie)*deltw 
00115     enddo 
00116     return
00117   end subroutine 
00118   !
00119   subroutine est_ef(ndosgrd,deltw,electron_number,dosgrd,dos,FermiEnergy)
00120     implicit none
00121     integer,intent(in)::ndosgrd
00122     real(8),intent(in)::deltw,electron_number 
00123     real(8),intent(in)::dosgrd(ndosgrd) 
00124     real(8),intent(in)::dos(ndosgrd) 
00125     real(8),intent(out)::FermiEnergy 
00126     real(8)::SUM_REAL
00127     integer::ie 
00128     real(8),parameter::au=27.21151d0 
00129     !
00130     SUM_REAL=0.0d0 
00131     do ie=1,ndosgrd-1
00132      !--
00133      !SUM_REAL=SUM_REAL+deltw*dos(ie) 
00134      !--
00135      !Trapezoidal 
00136      SUM_REAL=SUM_REAL+deltw*(dos(ie)+dos(ie+1))/2.0d0  
00137     enddo 
00138     !
00139     write(6,'(a40)')'+++ m_dos: est_ef +++'
00140     write(6,'(a40,f15.8)')'SUM of DOS',SUM_REAL
00141     write(6,'(a40,f15.8)')'FermiEnergy(before)=',FermiEnergy*au  
00142     SUM_REAL=0.0d0 
00143     do ie=1,ndosgrd-1
00144      if(SUM_REAL>=electron_number)goto 3000 
00145        !--
00146        !SUM_REAL=SUM_REAL+deltw*dos(ie) 
00147        !--
00148        !Trapezoidal 
00149        SUM_REAL=SUM_REAL+deltw*(dos(ie)+dos(ie+1))/2.0d0  
00150     enddo 
00151   3000 FermiEnergy=dosgrd(ie) 
00152     write(6,'(a40,f15.8)')'FermiEnergy(after)=',dosgrd(ie)*au  
00153     write(6,'(a40,f15.8)')'electron_number',SUM_REAL
00154     return
00155   end subroutine 
00156   !
00157   subroutine make_efline(ndosgrd,FermiEnergy,deltw,dosgrd,dos,efline) 
00158     implicit none 
00159     integer,intent(in)::ndosgrd
00160     real(8),intent(in)::FermiEnergy,deltw
00161     real(8),intent(in)::dosgrd(ndosgrd) 
00162     real(8),intent(in)::dos(ndosgrd) 
00163     real(8),intent(out)::efline(ndosgrd) 
00164     real(8)::diff,diff_old,dosmax,N0ttr 
00165     integer::ie,ief 
00166     real(8)::SUM_REAL 
00167     real(8),parameter::au=27.21151d0
00168     !
00169     !efline(ndosgrd) 
00170     !
00171     diff_old=1.0d0 
00172     do ie=1,ndosgrd
00173      diff=dabs(dosgrd(ie)-FermiEnergy) 
00174      if(diff_old>diff) then 
00175       ief=ie
00176       diff_old=diff 
00177      endif  
00178     enddo!ie  
00179     dosmax=1.2d0*maxval(abs(dos(:)))/au  
00180     dosmax=dble(nint(dosmax*10.0d0))/10.0d0 
00181     do ie=1,ndosgrd 
00182      if(ie<=ief) then 
00183       efline(ie)=dosmax
00184      else 
00185       efline(ie)=0.0d0 
00186      endif 
00187     enddo!ie 
00188     !
00189     !Integral dos(w) dw
00190     !
00191     SUM_REAL=0.0d0 
00192     do ie=1,ndosgrd-1 
00193      !--
00194      !SUM_REAL=SUM_REAL+dos(ie)*deltw 
00195      !--
00196      !Trapezoidal 
00197      SUM_REAL=SUM_REAL+deltw*(dos(ie)+dos(ie+1))/2.0d0  
00198      !write(6,'(2f15.10)') dosgrd(ie)*au,dos(ie)/au  
00199     enddo 
00200     N0ttr=dos(ief)/2.0d0 !2 is spin
00201     write(6,*)
00202     write(6,'(a40)')'+++ m_dos: est_efline +++'
00203     write(6,'(a40,f15.8)')'Integral dos(w) dw',SUM_REAL 
00204     write(6,'(a40,f15.8)')'N(0) in au per calc cell',N0ttr 
00205     write(6,'(a40,f15.8)')'N(0) in eV per calc cell',N0ttr/au   
00206     write(6,*)
00207     return
00208   end subroutine 
00209   !
00210   subroutine calc_dos(NTB,NTK,nkb1,nkb2,nkb3,ndosgrd,dosgrd,E_EIG,SK0,delt,dmnr,dmna,b1,b2,b3,dos) 
00211     use m_tetrahedron
00212     implicit none
00213     integer,intent(in)::NTB,NTK,nkb1,nkb2,nkb3,ndosgrd
00214     real(8),intent(in)::dosgrd(ndosgrd) 
00215     real(8),intent(in)::E_EIG(NTB,NTK)           
00216     real(8),intent(in)::SK0(3,NTK)           
00217     real(8),intent(in)::delt,dmnr,dmna 
00218     real(8),intent(in)::b1(3),b2(3),b3(3)
00219     real(8),intent(out)::dos(ndosgrd) 
00220     integer::index_kpt(nkb1,nkb2,nkb3) 
00221     integer::imt1(4*nkb1*nkb2*nkb3*6)  
00222     integer::ie,jb,ik,ikb1,ikb2,ikb3
00223     integer::iomp,omp_get_thread_num  
00224     real(8)::SUM_REAL 
00225     complex(8)::fk_1D(NTK)
00226     complex(8)::fk_3D(nkb1,nkb2,nkb3) 
00227     complex(8)::gk_1D(NTK)
00228     complex(8)::gk_3D(nkb1,nkb2,nkb3)
00229     complex(8)::xo(nkb1,nkb2,nkb3) 
00230     complex(8)::xow(ndosgrd,nkb1,nkb2,nkb3)
00231     real(8),parameter::pi=dacos(-1.0d0)
00232     real(8),parameter::au=27.21151d0 
00233     !
00234     call make_index_kpt(NTK,nkb1,nkb2,nkb3,SK0(1,1),index_kpt(1,1,1)) 
00235     call ttrhdrn_mkidx(nkb1,nkb2,nkb3,imt1(1),b1(1),b2(1),b3(1))
00236     xow=0.0d0 
00237 !$OMP PARALLEL PRIVATE(ie,jb,fk_1D,gk_1D,ik,fk_3D,gk_3D,&
00238 !$OMP ikb3,ikb2,ikb1,xo,iomp) 
00239 !$OMP DO 
00240     do ie=1,ndosgrd 
00241      do jb=1,NTB
00242       fk_1D=0.0d0
00243       gk_1D=0.0d0
00244       do ik=1,NTK 
00245        fk_1D(ik)=1.0d0 
00246        gk_1D(ik)=cmplx(dosgrd(ie)-E_EIG(jb,ik),-delt) 
00247       enddo!ik 
00248       fk_3D=0.0d0 
00249       gk_3D=0.0d0 
00250       do ikb3=1,nkb3
00251        do ikb2=1,nkb2
00252         do ikb1=1,nkb1 
00253          ik=index_kpt(ikb1,ikb2,ikb3) 
00254          fk_3D(ikb1,ikb2,ikb3)=fk_1D(ik)
00255          gk_3D(ikb1,ikb2,ikb3)=gk_1D(ik)
00256         enddo!ikb1 
00257        enddo!ikb2 
00258       enddo!ikb3 
00259       xo=0.0d0 
00260       call ttrhdrn_simple(dmna,dmnr,nkb1,nkb2,nkb3,imt1(1),& 
00261        fk_3D(1,1,1),gk_3D(1,1,1),xo(1,1,1))
00262       xow(ie,:,:,:)=xow(ie,:,:,:)+xo(:,:,:) 
00263      enddo!jb 
00264      iomp=omp_get_thread_num() 
00265      !if(iomp.eq.0) write(6,*) '#',ie  
00266     enddo!ie
00267 !$OMP END DO 
00268 !$OMP END PARALLEL 
00269     dos=0.0d0 
00270     do ie=1,ndosgrd 
00271      SUM_REAL=0.0d0 
00272      do ikb3=1,nkb3
00273       do ikb2=1,nkb2
00274        do ikb1=1,nkb1 
00275         SUM_REAL=SUM_REAL+dabs(dimag(xow(ie,ikb1,ikb2,ikb3)))/pi  
00276        enddo 
00277       enddo 
00278      enddo 
00279      dos(ie)=2.0d0*SUM_REAL/dble(NTK)!2 is spin 
00280     enddo 
00281     !
00282     return
00283   end subroutine calc_dos
00284   !
00285   subroutine calc_pdos(NTB,NTK,nkb1,nkb2,nkb3,ndosgrd,dosgrd,E_EIG,VKS,SK0,delt,dmnr,dmna,b1,b2,b3,pdos) 
00286     use m_tetrahedron
00287     implicit none
00288     integer,intent(in)::NTB,NTK,nkb1,nkb2,nkb3,ndosgrd
00289     real(8),intent(in)::dosgrd(ndosgrd) 
00290     real(8),intent(in)::E_EIG(NTB,NTK)           
00291     real(8),intent(in)::SK0(3,NTK)           
00292     real(8),intent(in)::delt,dmnr,dmna 
00293     real(8),intent(in)::b1(3),b2(3),b3(3)
00294     !
00295     complex(8),intent(in)::VKS(NTB,NTB,NTK)
00296     complex(8)::pxow(ndosgrd,NTB,nkb1,nkb2,nkb3)
00297     real(8),intent(out)::pdos(ndosgrd,NTB) 
00298     !
00299     integer::index_kpt(nkb1,nkb2,nkb3) 
00300     integer::imt1(4*nkb1*nkb2*nkb3*6)  
00301     integer::ie,ib,jb,ik,ikb1,ikb2,ikb3
00302     integer::iomp,omp_get_thread_num  
00303     real(8)::SUM_REAL 
00304     complex(8)::fk_1D(NTK)
00305     complex(8)::fk_3D(nkb1,nkb2,nkb3) 
00306     complex(8)::gk_1D(NTK)
00307     complex(8)::gk_3D(nkb1,nkb2,nkb3)
00308     complex(8)::xo(nkb1,nkb2,nkb3) 
00309     real(8),parameter::pi=dacos(-1.0d0)
00310     real(8),parameter::au=27.21151d0 
00311     !
00312     call make_index_kpt(NTK,nkb1,nkb2,nkb3,SK0(1,1),index_kpt(1,1,1)) 
00313     call ttrhdrn_mkidx(nkb1,nkb2,nkb3,imt1(1),b1(1),b2(1),b3(1))
00314     pxow=0.0d0 
00315 !$OMP PARALLEL PRIVATE(ie,jb,fk_1D,gk_1D,ik,fk_3D,gk_3D,&
00316 !$OMP ikb3,ikb2,ikb1,xo,iomp) 
00317 !$OMP DO 
00318     do ie=1,ndosgrd 
00319      do jb=1,NTB
00320       fk_1D=0.0d0
00321       gk_1D=0.0d0
00322       do ik=1,NTK 
00323        fk_1D(ik)=1.0d0 
00324        gk_1D(ik)=cmplx(dosgrd(ie)-E_EIG(jb,ik),-delt) 
00325       enddo!ik 
00326       fk_3D=0.0d0 
00327       gk_3D=0.0d0 
00328       do ikb3=1,nkb3
00329        do ikb2=1,nkb2
00330         do ikb1=1,nkb1 
00331          ik=index_kpt(ikb1,ikb2,ikb3) 
00332          fk_3D(ikb1,ikb2,ikb3)=fk_1D(ik)
00333          gk_3D(ikb1,ikb2,ikb3)=gk_1D(ik)
00334         enddo!ikb1 
00335        enddo!ikb2 
00336       enddo!ikb3 
00337       xo=0.0d0 
00338       call ttrhdrn_simple(dmna,dmnr,nkb1,nkb2,nkb3,imt1(1),& 
00339        fk_3D(1,1,1),gk_3D(1,1,1),xo(1,1,1))
00340       pxow(ie,jb,:,:,:)=xo(:,:,:) 
00341      enddo!jb 
00342      iomp=omp_get_thread_num() 
00343      !if(iomp.eq.0) write(6,*) '#',ie  
00344     enddo!ie
00345 !$OMP END DO 
00346 !$OMP END PARALLEL 
00347     pdos=0.0d0 
00348     do ie=1,ndosgrd 
00349      do ib=1,NTB
00350       SUM_REAL=0.0d0 
00351       do jb=1,NTB 
00352        do ikb3=1,nkb3
00353         do ikb2=1,nkb2
00354          do ikb1=1,nkb1 
00355           ik=index_kpt(ikb1,ikb2,ikb3) 
00356           SUM_REAL=SUM_REAL+(abs(VKS(jb,ib,ik))**2)*dabs(dimag(pxow(ie,jb,ikb1,ikb2,ikb3)))/pi  
00357          enddo!ikb1 
00358         enddo!ikb2 
00359        enddo!ikb3 
00360       enddo!jb  
00361       pdos(ie,ib)=2.0d0*SUM_REAL/dble(NTK)!2 is spin 
00362      enddo!ib 
00363     enddo!ie 
00364     !
00365     return
00366   end subroutine calc_pdos 
00367   !
00368   subroutine make_index_kpt(NTK,nkb1,nkb2,nkb3,SK0,index_kpt)      
00369     implicit none 
00370     integer::NTK,nkb1,nkb2,nkb3
00371     real(8)::SK0(3,NTK)  
00372     integer::ik,ix,iy,iz
00373     real(8)::x,y,z
00374     integer::index_kpt(nkb1,nkb2,nkb3)    
00375     ! 
00376     !if(MOD(NTK,2)/=0)then 
00377     ! !write(6,*)'i am in make_index for odd'
00378     ! do ik=1,NTK 
00379     !  x=SK0(1,ik)*dble(nkb1) 
00380     !  y=SK0(2,ik)*dble(nkb2)
00381     !  z=SK0(3,ik)*dble(nkb3)  
00382     !  x=x+(dble(nkb1)-1.0d0)/2.0d0 
00383     !  y=y+(dble(nkb2)-1.0d0)/2.0d0
00384     !  z=z+(dble(nkb3)-1.0d0)/2.0d0 
00385     !  ix=idnint(x)+1
00386     !  iy=idnint(y)+1
00387     !  iz=idnint(z)+1
00388     !  index_kpt(ix,iy,iz)=ik
00389     ! enddo 
00390     !else!20170316 
00391     ! !write(6,*)'i am in make_index for even'
00392     ! do ik=1,NTK 
00393     !  x=SK0(1,ik)*dble(nkb1) 
00394     !  y=SK0(2,ik)*dble(nkb2)
00395     !  z=SK0(3,ik)*dble(nkb3)  
00396     !  x=x+dble(nkb1)/2.0d0 
00397     !  y=y+dble(nkb2)/2.0d0
00398     !  z=z+dble(nkb3)/2.0d0 
00399     !  ix=idnint(x)
00400     !  iy=idnint(y)
00401     !  iz=idnint(z)
00402     !  index_kpt(ix,iy,iz)=ik
00403     ! enddo 
00404     !endif 
00405     !--
00406     !
00407     !20190520 Kazuma Nakamura
00408     !
00409     do ik=1,NTK 
00410      x=SK0(1,ik)*dble(nkb1) 
00411      y=SK0(2,ik)*dble(nkb2)
00412      z=SK0(3,ik)*dble(nkb3)  
00413      x=x+(dble(nkb1)-dble(mod(nkb1,2)))/2.0d0 
00414      y=y+(dble(nkb2)-dble(mod(nkb2,2)))/2.0d0 
00415      z=z+(dble(nkb3)-dble(mod(nkb3,2)))/2.0d0 
00416      ix=idnint(x)+mod(nkb1,2)
00417      iy=idnint(y)+mod(nkb2,2)
00418      iz=idnint(z)+mod(nkb3,2)
00419      index_kpt(ix,iy,iz)=ik
00420     enddo 
00421     ! 
00422     return 
00423   end subroutine  
00424   !
00425   subroutine wrt_dos(filename,ndosgrd,dosgrd,dos,efline) 
00426     implicit none
00427     integer,intent(in)::ndosgrd
00428     real(8),intent(in)::dosgrd(ndosgrd)
00429     real(8),intent(in)::dos(ndosgrd)
00430     real(8),intent(in)::efline(ndosgrd)
00431     character(99),intent(in)::filename
00432     character(99)::fname,fname_arg 
00433     integer::ie 
00434     real(8),parameter::au=27.21151d0 
00435     !
00436     !OPEN(300,W,file='./dir-wan/dat.dos.xxx-total')
00437     !
00438     fname_arg='-total'
00439     !
00440     write(fname,'(a,a)') trim(filename),trim(fname_arg) 
00441     OPEN(300,file=trim(fname)) 
00442     rewind(300) 
00443     do ie=1,ndosgrd
00444      write(300,'(3f15.10)') dosgrd(ie)*au,dos(ie)/au,efline(ie)  
00445     enddo!ie 
00446     close(300)
00447     return
00448   end subroutine 
00449   !
00450   subroutine wrt_pdos(filename,NTB,ndosgrd,dosgrd,pdos,efline) 
00451     implicit none
00452     integer,intent(in)::NTB
00453     integer,intent(in)::ndosgrd
00454     real(8),intent(in)::dosgrd(ndosgrd)
00455     real(8),intent(in)::pdos(ndosgrd,NTB)
00456     real(8),intent(in)::efline(ndosgrd)
00457     character(99),intent(in)::filename
00458     character(99)::fname,fname_arg 
00459     integer::ie,ib 
00460     real(8),parameter::au=27.21151d0 
00461     !
00462     !OPEN(400,W,file='./dir-wan/dat.dos.wannier-.ib')
00463     !
00464     do ib=1,NTB 
00465      fname_arg='-'
00466      write(fname,'(a,a,i3.3)') trim(filename),trim(fname_arg),ib
00467      OPEN(400,FILE=TRIM(fname)) 
00468      rewind(400)  
00469      do ie=1,ndosgrd
00470       write(400,'(3f15.10)') dosgrd(ie)*au,pdos(ie,ib)/au,efline(ie)  
00471      enddo!ie 
00472      close(400)
00473     enddo!ib 
00474     return
00475   end subroutine wrt_pdos 
00476   !
00477 end module m_dos