m_dos.f90

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