calc_gwakw.f90

Go to the documentation of this file.
00001 subroutine calc_gwakw(NWF,NTK,nsgm,Na1,Na2,Na3,nkb1,nkb2,nkb3,NSK_BAND_DISP,idlt,shift_ef,SK_BAND_DISP,a1,a2,a3,&
00002   sgmw,KS_R,XC_R,SX_R,SC_R,AKW,gw_sigma_kw)
00003   !
00004   implicit none 
00005   integer::NWF,NTK,nsgm,Na1,Na2,Na3,nkb1,nkb2,nkb3,NSK_BAND_DISP
00006   real(8)::idlt,shift_ef 
00007   real(8)::SK_BAND_DISP(3,NSK_BAND_DISP)
00008   real(8)::a1(3),a2(3),a3(3)
00009   real(8)::sgmw(nsgm)  
00010   complex(8)::KS_R(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3)   
00011   complex(8)::XC_R(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3)   
00012   complex(8)::SX_R(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3)   
00013   complex(4)::SC_R(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3,nsgm)   
00014   !
00015   real(8),allocatable::WEIGHT_R(:,:,:)!WEIGHT_R(-Na1:Na1,-Na2:Na2,-Na3:Na3)
00016   complex(8),allocatable::pf(:,:,:,:)!pf(-Na1:Na1,-Na2:Na2,-Na3:Na3,NSK_BAND_DISP)   
00017   complex(8),allocatable::EMK(:,:,:)!EMK(NWF,NSK_BAND_DISP,nsgm)           
00018   complex(4),allocatable::GW_R(:,:,:,:,:,:)!GW_R(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3,nsgm)   
00019   !
00020   integer::ik,ib,jb,ia1,ia2,ia3,ie,ia1min,ia2min,ia3min 
00021   real(8)::PHASE           
00022   real(8)::SUM_REAL 
00023   complex(8)::w,d,en 
00024   !
00025   real(8),parameter::au=27.21151d0
00026   real(8),parameter::tpi=2.0d0*dacos(-1.0d0)
00027   complex(8),parameter::ci=(0.0D0,1.0D0) 
00028   !
00029   real(8),intent(out)::AKW(NSK_BAND_DISP,nsgm)           
00030   real(8),intent(out)::gw_sigma_kw(NSK_BAND_DISP,nsgm)           
00031   ! 
00032   !1. GW_R 
00033   !
00034   allocate(GW_R(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3,nsgm));GW_R=0.0d0   
00035   do ia1=-Na1,Na1 
00036    do ia2=-Na2,Na2 
00037     do ia3=-Na3,Na3 
00038      do ib=1,NWF          
00039       do jb=1,NWF          
00040        do ie=1,nsgm
00041         GW_R(ib,jb,ia1,ia2,ia3,ie)&
00042        =KS_R(ib,jb,ia1,ia2,ia3)&
00043        -XC_R(ib,jb,ia1,ia2,ia3)&
00044        -SX_R(ib,jb,ia1,ia2,ia3)&
00045        +SC_R(ib,jb,ia1,ia2,ia3,ie) 
00046        enddo 
00047       enddo
00048      enddo
00049     enddo
00050    enddo
00051   enddo
00052   write(6,*)'# finish make H(R)'
00053   !
00054   !2. WEIGHT_R
00055   !
00056   allocate(WEIGHT_R(-Na1:Na1,-Na2:Na2,-Na3:Na3)); WEIGHT_R=1.0d0
00057   SUM_REAL=0.0d0 
00058   do ia1=-Na1,Na1
00059    do ia2=-Na2,Na2
00060     do ia3=-Na3,Na3
00061      if(abs(ia1)==Na1.and.mod(nkb1,2)==0.and.Na1/=0) WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0 
00062      if(abs(ia2)==Na2.and.mod(nkb2,2)==0.and.Na2/=0) WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0 
00063      if(abs(ia3)==Na3.and.mod(nkb3,2)==0.and.Na3/=0) WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0 
00064      SUM_REAL=SUM_REAL+WEIGHT_R(ia1,ia2,ia3)
00065     enddo
00066    enddo
00067   enddo 
00068   write(6,'(a20,f15.8,i8)')'SUM_WEIGHT,NTK',SUM_REAL,NTK  
00069   if(abs(SUM_REAL-dble(NTK))>1.0d-6)then 
00070    stop 'SUM_WEIGHT/=NTK'
00071   endif 
00072   !
00073   !3. PHASE FACTOR 
00074   !
00075   allocate(pf(-Na1:Na1,-Na2:Na2,-Na3:Na3,NSK_BAND_DISP)); pf=0.0d0 
00076   do ik=1,NSK_BAND_DISP          
00077    do ia3=-Na3,Na3 
00078     do ia2=-Na2,Na2 
00079      do ia1=-Na1,Na1 
00080       !
00081       !NEAREST R SEARCH BY Y.Yoshimoto 
00082       !
00083       call search_Rmin(ia1,ia2,ia3,nkb1,nkb2,nkb3,a1(1),a2(1),a3(1),ia1min,ia2min,ia3min)
00084       PHASE=tpi*(SK_BAND_DISP(1,ik)*DBLE(ia1min)+SK_BAND_DISP(2,ik)*DBLE(ia2min)+SK_BAND_DISP(3,ik)*DBLE(ia3min))           
00085       pf(ia1,ia2,ia3,ik)=EXP(ci*PHASE)*WEIGHT_R(ia1,ia2,ia3)          
00086      enddo 
00087     enddo 
00088    enddo 
00089   enddo 
00090   write(6,*)'# finish make pf'
00091   !
00092   !4. HGW IN WANNIER BASIS. AND DIAGONALIZE
00093   !
00094   allocate(EMK(NWF,NSK_BAND_DISP,nsgm)); EMK=0.0d0 
00095   call make_emk(NSK_BAND_DISP,NWF,nsgm,Na1,Na2,Na3,GW_R(1,1,-Na1,-Na2,-Na3,1),pf(-Na1,-Na2,-Na3,1),EMK(1,1,1))  
00096   !
00097   !5. CALC SPECTRAL FUNCTION 
00098   !
00099   AKW=0.0d0
00100   do ik=1,NSK_BAND_DISP 
00101    do jb=1,NWF 
00102     do ie=1,nsgm 
00103      w=cmplx(sgmw(ie)+shift_ef,-idlt) 
00104      en=EMK(jb,ik,ie)
00105      d=1.0d0/(w-en) 
00106      AKW(ik,ie)=AKW(ik,ie)+abs(imag(d)) 
00107     enddo 
00108    enddo 
00109   enddo 
00110   !
00111   !6. CALC SIGMA_GW SPECTRAL FUNCTION 
00112   !
00113   gw_sigma_kw=0.0d0
00114   do ie=1,nsgm 
00115    do ik=1,NSK_BAND_DISP 
00116     do jb=1,NWF 
00117      gw_sigma_kw(ik,ie)=gw_sigma_kw(ik,ie)+abs(imag(EMK(jb,ik,ie))) 
00118     enddo!jb 
00119    enddo!ik 
00120   enddo!ie  
00121   !
00122   deallocate(WEIGHT_R,pf,EMK,GW_R) 
00123   !
00124 return
00125 end