m_eigenstate.f90

Go to the documentation of this file.
00001 module m_eigenstate  
00002   implicit none
00003 contains
00004   !
00005   subroutine calculate_eigenstate(NWF,Na1,Na2,Na3,HKS,NTK,nkb1,nkb2,nkb3,a1,a2,a3,flg_weight,Ncalck,SK0,EKS,VKS) 
00006     implicit none 
00007     integer,intent(in)::NWF,Na1,Na2,Na3
00008     complex(8),intent(in)::HKS(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3)   
00009     integer,intent(in)::NTK,nkb1,nkb2,nkb3 
00010     real(8),intent(in)::a1(3),a2(3),a3(3)
00011     integer,intent(in)::flg_weight 
00012     integer,intent(in)::Ncalck
00013     real(8),intent(in)::SK0(3,Ncalck) 
00014     real(8),intent(out)::EKS(NWF,Ncalck)           
00015     complex(8),intent(out)::VKS(NWF,NWF,Ncalck)   
00016     real(8),allocatable::WEIGHT_R(:,:,:)!WEIGHT_R(-Na1:Na1,-Na2:Na2,-Na3:Na3)
00017     complex(8),allocatable::pf(:,:,:,:)!pf(-Na1:Na1,-Na2:Na2,-Na3:Na3,NTK)   
00018     integer::ia1min,ia2min,ia3min 
00019     integer::ia1,ia2,ia3,ik 
00020     real(8)::PHASE 
00021     real(8)::SUM_REAL
00022     real(8),parameter::au=27.21151d0
00023     real(8),parameter::tpi=2.0d0*dacos(-1.0d0)
00024     complex(8),parameter::ci=(0.0D0,1.0D0) 
00025     !
00026     !WEIGHT_R BY Y.Nomura NOMURA 
00027     !
00028     allocate(WEIGHT_R(-Na1:Na1,-Na2:Na2,-Na3:Na3)); WEIGHT_R=1.0d0
00029     if(flg_weight.eq.1)then
00030      write(6,'(a40)')'WEIGHT CALCULATED' 
00031      call make_WEIGHT_R(nkb1,nkb2,nkb3,Na1,Na2,Na3,NTK,WEIGHT_R(-Na1,-Na2,-Na3)) 
00032     else
00033      write(6,'(a40)')'WEIGHT NOT CALCULATED' 
00034     endif 
00035     !
00036     !phase factor
00037     !
00038     allocate(pf(-Na1:Na1,-Na2:Na2,-Na3:Na3,Ncalck)); pf=0.0d0 
00039     do ik=1,Ncalck 
00040      do ia3=-Na3,Na3 
00041       do ia2=-Na2,Na2 
00042        do ia1=-Na1,Na1 
00043         !
00044         !NEAREST R SEARCH BY Y.Yoshimoto 
00045         !
00046         call search_Rmin(ia1,ia2,ia3,nkb1,nkb2,nkb3,a1(1),a2(1),a3(1),ia1min,ia2min,ia3min)
00047         PHASE=tpi*(SK0(1,ik)*DBLE(ia1min)+SK0(2,ik)*DBLE(ia2min)+SK0(3,ik)*DBLE(ia3min))           
00048         pf(ia1,ia2,ia3,ik)=EXP(ci*PHASE)*WEIGHT_R(ia1,ia2,ia3)          
00049        enddo!ia1 
00050       enddo!ia2 
00051      enddo!ia3 
00052     enddo!ik  
00053     !
00054     !HKS IN WANNIER BASIS. AND DIAGONALIZE
00055     !
00056     EKS=0.0d0 
00057     VKS=0.0d0 
00058     call eigenvalue(Ncalck,NWF,Na1,Na2,Na3,HKS(1,1,-Na1,-Na2,-Na3),pf(-Na1,-Na2,-Na3,1),EKS(1,1),VKS(1,1,1))  
00059     !
00060     deallocate(WEIGHT_R,pf) 
00061     !
00062     return 
00063   end subroutine calculate_eigenstate 
00064 
00065   subroutine search_Rmin(i,j,k,nkb1,nkb2,nkb3,a1,a2,a3,imin,jmin,kmin)
00066     implicit none
00067     integer::i,j,k,nkb1,nkb2,nkb3
00068     integer::imin,jmin,kmin
00069     real(8)::a1(3),a2(3),a3(3) 
00070     integer::nmin,mmin,lmin
00071     integer::n,m,l 
00072     real(8)::R_pos(3),R_abs,R_min,R_bfr
00073     R_pos(:)=dble(i)*a1(:)+dble(j)*a2(:)+dble(k)*a3(:)
00074     nmin=0;mmin=0;lmin=0
00075     R_bfr=dsqrt(R_pos(1)**2+R_pos(2)**2+R_pos(3)**2) 
00076     R_min=R_bfr 
00077     do n=-3,3
00078      do m=-3,3
00079       do l=-3,3
00080        R_pos(:)=dble(i+n*nkb1)*a1(:)+dble(j+m*nkb2)*a2(:)+dble(k+l*nkb3)*a3(:)
00081        R_abs=dsqrt(R_pos(1)**2+R_pos(2)**2+R_pos(3)**2) 
00082        if(R_min>R_abs)then 
00083         R_min=R_abs
00084         nmin=n
00085         mmin=m
00086         lmin=l
00087        endif 
00088       enddo
00089      enddo
00090     enddo
00091     imin=i+nmin*nkb1
00092     jmin=j+mmin*nkb2
00093     kmin=k+lmin*nkb3
00094     return
00095   end subroutine 
00096 
00097   subroutine eigenvalue(NTK,NWF,Na1,Na2,Na3,HmatR,pf,EKS,VKS)  
00098     implicit none 
00099     integer,intent(in)::NTK,NWF,Na1,Na2,Na3 
00100     complex(8),intent(in)::pf(-Na1:Na1,-Na2:Na2,-Na3:Na3,NTK)   
00101     complex(8),intent(in)::HmatR(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3)   
00102     real(8),intent(inout)::EKS(NWF,NTK)           
00103     complex(8),intent(inout)::VKS(NWF,NWF,NTK)  
00104     complex(8),allocatable::Hin(:,:)!Hin(NWF,NWF)          
00105     real(8),allocatable::E_TMP_R(:)!E_TMP_R(NWF)           
00106     integer::ik,ib,jb,ia1,ia2,ia3
00107     complex(8)::SUM_CMPX 
00108     allocate(Hin(NWF,NWF)); Hin=0.0d0           
00109     allocate(E_TMP_R(NWF)); E_TMP_R=0.0d0            
00110     !
00111     EKS=0.0d0 
00112     VKS=0.0d0 
00113     do ik=1,NTK 
00114      Hin(:,:)=0.0D0               
00115      do ib=1,NWF      
00116       do jb=1,NWF      
00117        SUM_CMPX=0.0D0                    
00118        do ia1=-Na1,Na1 
00119         do ia2=-Na2,Na2 
00120          do ia3=-Na3,Na3 
00121           SUM_CMPX=SUM_CMPX+HmatR(ib,jb,ia1,ia2,ia3)*pf(ia1,ia2,ia3,ik)  
00122          enddo!ia3            
00123         enddo!ia2                       
00124        enddo!ia1
00125        Hin(ib,jb)=SUM_CMPX           
00126       enddo!jb 
00127      enddo!ib 
00128      E_TMP_R(:)=0.0D0 
00129      call diagV(NWF,Hin(1,1),E_TMP_R(1)) 
00130      EKS(:,ik)=E_TMP_R(:) 
00131      !
00132      !H(ib,jb):: ib: basis, jb: eigenvector
00133      !
00134      ![NOTE] trnspose 
00135      !
00136      !VKS(jb,ib):: jb: eigenvector, ib: basis 
00137      !
00138      do ib=1,NWF
00139       do jb=1,NWF
00140        VKS(jb,ib,ik)=Hin(ib,jb) 
00141       enddo
00142      enddo
00143      !
00144     enddo!ik 
00145     deallocate(Hin,E_TMP_R) 
00146     !
00147     return
00148   end subroutine 
00149 
00150   subroutine diagV(nm,mat,eig)
00151     implicit none 
00152     integer,intent(in)::nm
00153     complex(8),intent(inout)::mat(nm,nm)
00154     real(8),intent(out)::eig(nm)
00155     integer::LWORK,LRWORK,LIWORK  
00156     integer,allocatable::iwork_zheevd(:)
00157     real(8),allocatable::rwork_zheevd(:)
00158     complex(8),allocatable::work_zheevd(:)
00159     integer::ind
00160     real(8)::eps 
00161     !
00162     LWORK= 2*nm+nm**2
00163     LRWORK=1+12*nm+3*nm**2
00164     LIWORK=3+10*nm 
00165     allocate(work_zheevd(LWORK));work_zheevd(:)=0.0d0
00166     allocate(rwork_zheevd(LRWORK));rwork_zheevd(:)=0.0d0
00167     allocate(iwork_zheevd(LIWORK));iwork_zheevd(:)=0
00168     eps=1.0d-18
00169     ind=0                 
00170     !
00171     call zheevd("V","U",nm,mat,nm,eig,work_zheevd,LWORK,rwork_zheevd,LRWORK,iwork_zheevd,LIWORK,ind)
00172     !
00173     if(ind/=0)then 
00174      write(6,'(a40,i15)')'ind=',ind 
00175      stop
00176     endif 
00177     !
00178     deallocate(work_zheevd,rwork_zheevd,iwork_zheevd) 
00179     return 
00180   end subroutine
00181   
00182   subroutine diagN(nm,mat,eig)
00183     implicit none 
00184     integer,intent(in)::nm
00185     complex(8),intent(inout)::mat(nm,nm)
00186     real(8),intent(out)::eig(nm)
00187     integer::LWORK,LRWORK,LIWORK  
00188     integer,allocatable::iwork_zheevd(:)
00189     real(8),allocatable::rwork_zheevd(:)
00190     complex(8),allocatable::work_zheevd(:)
00191     integer::ind
00192     real(8)::eps 
00193     !
00194     LWORK= 2*nm+nm**2
00195     LRWORK=1+12*nm+3*nm**2
00196     LIWORK=3+10*nm 
00197     allocate(work_zheevd(LWORK));work_zheevd(:)=0.0d0
00198     allocate(rwork_zheevd(LRWORK));rwork_zheevd(:)=0.0d0
00199     allocate(iwork_zheevd(LIWORK));iwork_zheevd(:)=0
00200     eps=1.0d-18
00201     ind=0                 
00202     !
00203     call zheevd("N","U",nm,mat,nm,eig,work_zheevd,LWORK,rwork_zheevd,LRWORK,iwork_zheevd,LIWORK,ind)
00204     !
00205     if(ind/=0)then 
00206      write(6,'(a40,i15)')'ind=',ind 
00207      stop
00208     endif 
00209     !
00210     deallocate(work_zheevd,rwork_zheevd,iwork_zheevd) 
00211     return 
00212   end subroutine
00213   !
00214   subroutine make_WEIGHT_R(nkb1,nkb2,nkb3,Na1,Na2,Na3,NTK,WEIGHT_R) 
00215     implicit none 
00216     integer,intent(in)::nkb1,nkb2,nkb3,Na1,Na2,Na3,NTK 
00217     real(8),intent(out)::WEIGHT_R(-Na1:Na1,-Na2:Na2,-Na3:Na3)
00218     integer::ia1,ia2,ia3
00219     real(8)::SUM_REAL 
00220     real(8),parameter::dlt_eps=1.0d-6 
00221     !
00222     SUM_REAL=0.0d0 
00223     do ia1=-Na1,Na1
00224      do ia2=-Na2,Na2
00225       do ia3=-Na3,Na3
00226        if(abs(ia1)==Na1.and.mod(nkb1,2)==0.and.Na1/=0) WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0 
00227        if(abs(ia2)==Na2.and.mod(nkb2,2)==0.and.Na2/=0) WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0 
00228        if(abs(ia3)==Na3.and.mod(nkb3,2)==0.and.Na3/=0) WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0 
00229        SUM_REAL=SUM_REAL+WEIGHT_R(ia1,ia2,ia3)
00230       enddo!ia3
00231      enddo!ia2
00232     enddo!ia1
00233     write(6,'(a40,f15.8,i8)')'SUM_WEIGHT,NTK',SUM_REAL,NTK  
00234     if(abs(SUM_REAL-dble(NTK))>dlt_eps)then 
00235      stop 'SUM_WEIGHT/=NTK'
00236     endif 
00237     !
00238     return 
00239   end subroutine make_WEIGHT_R 
00240   !
00241 end module m_eigenstate