sub_ksdos.f90

Go to the documentation of this file.
00001 subroutine calc_dos_KS(NTB,NTK,nkb1,nkb2,nkb3,ndosgrd,dosgrd,E_EIG,SK0,delt,dmnr,dmna,b1,b2,b3,dos) 
00002   use m_tetrahedron
00003   implicit none
00004   integer::NTB,NTK,nkb1,nkb2,nkb3,ndosgrd
00005   real(8)::dosgrd(ndosgrd) 
00006   real(8)::E_EIG(NTB,NTK)           
00007   real(8)::SK0(3,NTK)           
00008   real(8)::delt,dmnr,dmna 
00009   real(8)::b1(3),b2(3),b3(3)
00010   real(8)::dos(ndosgrd) 
00011   integer::ie,jb,ik,ikb1,ikb2,ikb3
00012   integer::iomp,omp_get_thread_num  
00013   real(8)::SUM_REAL 
00014   !
00015   !ttrhdrn
00016   !
00017   integer,allocatable::imt1(:)!imt1(4*nkb1*nkb2*nkb3*6)  
00018   integer,allocatable::index_kpt(:,:,:)!index_kpt(nkb1,nkb2,nkb3)    
00019   complex(8),allocatable::fk_1D(:)!fk_1D(NTK)
00020   complex(8),allocatable::gk_1D(:)!gk_1D(NTK)
00021   complex(8),allocatable::fk_3D(:,:,:)!fk_3D(nkb1,nkb2,nkb3) 
00022   complex(8),allocatable::gk_3D(:,:,:)!gk_3D(nkb1,nkb2,nkb3)
00023   complex(8),allocatable::xo(:,:,:)!xo(nkb1,nkb2,nkb3) 
00024   complex(8),allocatable::xow(:,:,:,:)!xow(nsgm,nkb1,nkb2,nkb3)
00025   !
00026   real(8),parameter::pi=dacos(-1.0d0)
00027   real(8),parameter::au=27.21151d0 
00028   !
00029   allocate(imt1(4*nkb1*nkb2*nkb3*6)); imt1=0   
00030   allocate(index_kpt(nkb1,nkb2,nkb3)); index_kpt=0  
00031   call make_index_kpt(NTK,nkb1,nkb2,nkb3,SK0(1,1),index_kpt(1,1,1)) 
00032   call ttrhdrn_mkidx(nkb1,nkb2,nkb3,imt1(1),b1(1),b2(1),b3(1))
00033   !  
00034   allocate(xow(ndosgrd,nkb1,nkb2,nkb3)); xow=0.0d0 
00035 !$OMP PARALLEL PRIVATE(ie,jb,fk_1D,gk_1D,ik,fk_3D,gk_3D,ikb3,ikb2,ikb1,xo,iomp) 
00036   allocate(fk_1D(NTK)); fk_1D=0.0d0 
00037   allocate(gk_1D(NTK)); gk_1D=0.0d0 
00038   allocate(fk_3D(nkb1,nkb2,nkb3)); fk_3D=0.0d0 
00039   allocate(gk_3D(nkb1,nkb2,nkb3)); gk_3D=0.0d0 
00040   allocate(xo(nkb1,nkb2,nkb3)); xo=0.0d0  
00041 !$OMP DO 
00042   do ie=1,ndosgrd 
00043    do jb=1,NTB
00044     fk_1D=0.0d0
00045     gk_1D=0.0d0
00046     do ik=1,NTK 
00047      fk_1D(ik)=1.0d0 
00048      gk_1D(ik)=cmplx(dosgrd(ie)-E_EIG(jb,ik),-delt) 
00049     enddo 
00050     fk_3D=0.0d0 
00051     gk_3D=0.0d0 
00052     do ikb3=1,nkb3
00053      do ikb2=1,nkb2
00054       do ikb1=1,nkb1 
00055        ik=index_kpt(ikb1,ikb2,ikb3) 
00056        fk_3D(ikb1,ikb2,ikb3)=fk_1D(ik)
00057        gk_3D(ikb1,ikb2,ikb3)=gk_1D(ik)
00058       enddo 
00059      enddo 
00060     enddo 
00061     xo=0.0d0 
00062     call ttrhdrn_simple(dmna,dmnr,nkb1,nkb2,nkb3,imt1(1),fk_3D(1,1,1),gk_3D(1,1,1),xo(1,1,1))
00063     xow(ie,:,:,:)=xow(ie,:,:,:)+xo(:,:,:) 
00064    enddo!jb 
00065    iomp=omp_get_thread_num() 
00066    if(iomp.eq.0) then 
00067     write(6,*)'#',ie  
00068    endif 
00069   enddo!ie
00070 !$OMP END DO 
00071   deallocate(fk_1D,gk_1D,fk_3D,gk_3D,xo) 
00072 !$OMP END PARALLEL 
00073   !
00074   dos=0.0d0 
00075   do ie=1,ndosgrd 
00076    SUM_REAL=0.0d0 
00077    do ikb3=1,nkb3
00078     do ikb2=1,nkb2
00079      do ikb1=1,nkb1 
00080       SUM_REAL=SUM_REAL+dabs(dimag(xow(ie,ikb1,ikb2,ikb3)))/pi  
00081      enddo 
00082     enddo 
00083    enddo 
00084    dos(ie)=2.0d0*SUM_REAL/dble(NTK)!2 is spin 
00085   enddo 
00086   !
00087   deallocate(imt1,index_kpt,xow) 
00088 return
00089 end