m_truncation.f90

Go to the documentation of this file.
00001 module m_truncation 
00002   implicit none
00003 contains
00004   subroutine truncation(NWF,Na1,Na2,Na3,threshold_e,threshold_r,diff_transfers,a1,a2,a3,wcenter_lat,KS_R)
00005     implicit none 
00006     integer,intent(in)::NWF,Na1,Na2,Na3 
00007     real(8),intent(in)::threshold_e,threshold_r,diff_transfers  
00008     real(8),intent(in)::a1(3),a2(3),a3(3) 
00009     real(8),intent(in)::wcenter_lat(3,NWF) 
00010     complex(8),intent(inout)::KS_R(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3)   
00011     complex(8)::tij
00012     complex(8),allocatable::tr(:)
00013     real(8),allocatable::dist(:)
00014     integer,allocatable::ib_map(:)
00015     integer,allocatable::jb_map(:)
00016     integer,allocatable::ia1_map(:)
00017     integer,allocatable::ia2_map(:)
00018     integer,allocatable::ia3_map(:)
00019     integer::ia1,ia2,ia3,ib,jb
00020     integer::Ntr,itr 
00021     real(8)::vec_rij(3),rij 
00022     !
00023     Ntr=NWF*NWF*(2*Na3+1)*(2*Na2+1)*(2*Na1+1) 
00024     allocate(tr(Ntr)); tr=0.0d0  
00025     allocate(dist(Ntr)); dist=0.0d0  
00026     allocate(ib_map(Ntr)); ib_map=0 
00027     allocate(jb_map(Ntr)); jb_map=0 
00028     allocate(ia1_map(Ntr)); ia1_map=0
00029     allocate(ia2_map(Ntr)); ia2_map=0
00030     allocate(ia3_map(Ntr)); ia3_map=0
00031     !
00032     do ia1=-Na1,Na1
00033      do ia2=-Na2,Na2
00034       do ia3=-Na3,Na3
00035        do ib=1,NWF
00036         do jb=1,NWF 
00037          tij=KS_R(ib,jb,ia1,ia2,ia3)
00038          !
00039          if(threshold_e>abs(tij))then 
00040           tij=0.0d0 
00041          endif 
00042          !
00043          vec_rij(:)=(dble(ia1)+wcenter_lat(1,jb)-wcenter_lat(1,ib))*a1(:)&
00044                    +(dble(ia2)+wcenter_lat(2,jb)-wcenter_lat(2,ib))*a2(:)&
00045                    +(dble(ia3)+wcenter_lat(3,jb)-wcenter_lat(3,ib))*a3(:)
00046          rij=dsqrt(vec_rij(1)**2+vec_rij(2)**2+vec_rij(3)**2) 
00047          if(threshold_r<abs(rij))then 
00048           !write(6,*)'rij=',rij 
00049           tij=0.0d0 
00050          endif 
00051          !
00052          KS_R(ib,jb,ia1,ia2,ia3)=tij 
00053          ! 
00054          itr=jb +(ib-1)*NWF +(ia3+Na3)*NWF*NWF +(ia2+Na2)*NWF*NWF*(2*Na3+1) +(ia1+Na1)*NWF*NWF*(2*Na3+1)*(2*Na2+1)  
00055          ib_map(itr)=ib
00056          jb_map(itr)=jb
00057          ia3_map(itr)=ia3
00058          ia2_map(itr)=ia2
00059          ia1_map(itr)=ia1
00060          tr(itr)=tij
00061          dist(itr)=rij
00062          !
00063         enddo!jb 
00064        enddo!ib 
00065       enddo!ia3 
00066      enddo!ia2 
00067     enddo!ia1 
00068     !
00069     call wrt_transfer_analysis(threshold_e,threshold_r,diff_transfers,Ntr,ib_map(1),jb_map(1),ia1_map(1),ia2_map(1),ia3_map(1),dist(1),tr(1)) 
00070     !
00071     deallocate(tr,dist,ib_map,jb_map,ia1_map,ia2_map,ia3_map)
00072     !
00073     return
00074   end subroutine truncation 
00075 
00076   subroutine wrt_transfer_analysis(threshold_e,threshold_r,diff_transfers,Ntr,ib_map,jb_map,ia1_map,ia2_map,ia3_map,dist,tr) 
00077     implicit none 
00078     integer,intent(in)::Ntr
00079     real(8),intent(in)::threshold_e,threshold_r,diff_transfers   
00080     integer,intent(in)::ib_map(Ntr)
00081     integer,intent(in)::jb_map(Ntr)
00082     integer,intent(in)::ia1_map(Ntr)
00083     integer,intent(in)::ia2_map(Ntr)
00084     integer,intent(in)::ia3_map(Ntr)
00085     real(8),intent(in)::dist(Ntr)
00086     complex(8),intent(in)::tr(Ntr)
00087     real(8)::tri,trj 
00088     integer::itr,jtr
00089     integer,allocatable::skip(:)
00090     real(8),parameter::au=27.21151d0 
00091     real(8),parameter::del_zero=(1.d-5)/au  
00092     !
00093     !OPEN(122,W,FILE='./dat.tr') 
00094     !
00095     OPEN(122,FILE='./dat.tr') 
00096     rewind(122)
00097     write(122,'(a30,x,f10.5)')'#Threshold for transfer (eV):',threshold_e*au
00098     write(122,'(a30,x,f10.5)')'#Threshold for distance (AA):',threshold_r 
00099     write(122,'(a30,x,f10.5)')'#difference in transfers (eV):',diff_transfers*au  
00100     write(122,'(a30)')'#i,j,ia1,ia2,ia3,tr,dist:'
00101     allocate(skip(Ntr)); skip=0
00102     do itr=1,Ntr 
00103      tri=abs(dble(tr(itr)))  
00104      if(tri<del_zero)cycle 
00105      if(skip(itr)==1)cycle 
00106      write(122,*) 
00107      write(122,'(a30)')'irreducible transfer:' 
00108      !write(122,'(5i5,2f15.7)') ib_map(itr),jb_map(itr),ia1_map(itr),ia2_map(itr),ia3_map(itr),dble(tr(itr))*au,dist(itr)
00109      write(122,'(5i3,2f10.5)') ib_map(itr),jb_map(itr),ia1_map(itr),ia2_map(itr),ia3_map(itr),dble(tr(itr))*au,dist(itr)
00110      write(122,*) 
00111      do jtr=1,Ntr 
00112       trj=abs(dble(tr(jtr)))  
00113       if(itr==jtr)cycle 
00114       if(skip(jtr)==1)cycle 
00115       if(trj<del_zero)cycle 
00116       if(abs(tri-trj)<diff_transfers)then 
00117        !write(122,'(5i5,2f15.7)') ib_map(jtr),jb_map(jtr),ia1_map(jtr),ia2_map(jtr),ia3_map(jtr),dble(tr(jtr))*au,dist(jtr)
00118        write(122,'(5i3,2f10.5)') ib_map(jtr),jb_map(jtr),ia1_map(jtr),ia2_map(jtr),ia3_map(jtr),dble(tr(jtr))*au,dist(jtr)
00119        !write(122,'(2i10,2f15.7)') itr,jtr,dble(tr(itr))*au,dble(tr(jtr))*au  
00120        skip(jtr)=1
00121       endif 
00122      enddo 
00123      skip(itr)=1 
00124     enddo 
00125     deallocate(skip)
00126     !
00127     return
00128   end subroutine wrt_transfer_analysis 
00129   !
00130   subroutine set_kgrid(dense,kvec)
00131     IMPLICIT NONE
00132     integer,intent(in)::dense(3)!Calculated k-grid 
00133     REAL(8),intent(out)::kvec(3,PRODUCT(dense(1:3))) 
00134     INTEGER::ik,i1,i2,i3,s1,s2,s3 
00135     !
00136     WRITE(6,'(a50,x,3i10)')'Calculated k-grid:',dense(1:3)
00137     !
00138     s1=dense(1)/2; s2=dense(2)/2; s3=dense(3)/2 
00139     ik = 0
00140     DO i1=-(s1-mod(dense(1)+1,2)),s1  
00141        DO i2=-(s2-mod(dense(2)+1,2)),s2  
00142           DO i3=-(s3-mod(dense(3)+1,2)),s3  
00143              ik=ik + 1
00144              kvec(1:3,ik)=DBLE((/i1, i2, i3/))/DBLE(dense(1:3))
00145           END DO
00146        END DO
00147     END DO
00148     !
00149     WRITE(6,'(a50,x,i10)')'Total number of k-grid:',ik 
00150     return
00151   end subroutine set_kgrid 
00152   !
00153   subroutine set_FR(NWF,nkb1,nkb2,nkb3,kgd1,kgd2,kgd3,Na1,Na2,Na3,La1,La2,La3,HR,FR) 
00154     implicit none 
00155     integer,intent(in)::NWF,nkb1,nkb2,nkb3,kgd1,kgd2,kgd3,Na1,Na2,Na3,La1,La2,La3 
00156     complex(8),intent(inout)::HR(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3) 
00157     complex(8),intent(out)::FR(NWF,NWF,-La1:La1,-La2:La2,-La3:La3) 
00158     real(8),allocatable::WR(:,:,:) 
00159     integer::i,j,k 
00160     !integer::ib,jb  
00161     !complex(8)::tij
00162     real(8),parameter::au=27.21151d0 
00163     !
00164     allocate(WR(-Na1:Na1,-Na2:Na2,-Na3:Na3)); WR=1.0 
00165     call get_weightR(nkb1,nkb2,nkb3,Na1,Na2,Na3,WR(-Na1,-Na2,-Na3))      
00166     !
00167     !rm weight factor
00168     !
00169     do i=-Na1,Na1
00170      do j=-Na2,Na2
00171       do k=-Na3,Na3 
00172        HR(:,:,i,j,k)=HR(:,:,i,j,k)/WR(i,j,k) 
00173       enddo
00174      enddo
00175     enddo
00176     deallocate(WR) 
00177     !
00178     !F(R)<--H(R) 
00179     !
00180     do i=-min(Na1,La1),min(Na1,La1)  
00181      do j=-min(Na2,La2),min(Na2,La2)  
00182       do k=-min(Na3,La3),min(Na3,La3)  
00183        FR(:,:,i,j,k)=HR(:,:,i,j,k) 
00184       enddo
00185      enddo
00186     enddo 
00187     !
00188     allocate(WR(-La1:La1,-La2:La2,-La3:La3)); WR=1.0 
00189     call get_weightR(kgd1,kgd2,kgd3,La1,La2,La3,WR(-La1,-La2,-La3))      
00190     ! 
00191     !add weight factor
00192     !
00193     do i=-La1,La1
00194      do j=-La2,La2
00195       do k=-La3,La3 
00196        FR(:,:,i,j,k)=FR(:,:,i,j,k)*WR(i,j,k) 
00197       enddo
00198      enddo
00199     enddo
00200     deallocate(WR) 
00201     !
00202     !output check 
00203     !
00204     !do i=-La1,La1
00205     ! do j=-La2,La2
00206     !  do k=-La3,La3 
00207     !   do ib=1,NWF
00208     !    do jb=1,NWF 
00209     !     tij=FR(ib,jb,i,j,k)
00210     !     write(6,'(5i5,2f15.7)') ib,jb,i,j,k,dble(tij)*au,abs(tij)*au
00211     !    enddo
00212     !   enddo 
00213     !  enddo
00214     ! enddo
00215     !enddo  
00216     ! 
00217     return
00218   end subroutine set_FR 
00219   !
00220   subroutine get_weightR(kgd1,kgd2,kgd3,Na1,Na2,Na3,WEIGHT_R) 
00221     implicit none 
00222     integer,intent(in)::kgd1,kgd2,kgd3,Na1,Na2,Na3 
00223     real(8),intent(out)::WEIGHT_R(-Na1:Na1,-Na2:Na2,-Na3:Na3) 
00224     integer::ia1,ia2,ia3 
00225     integer::ncalck 
00226     real(8)::SUM_REAL
00227     !
00228     !WEIGHT_R BY Y.Nomura NOMURA 
00229     !
00230     SUM_REAL=0.0d0 
00231     WEIGHT_R=1.0d0  
00232     ncalck=kgd1*kgd2*kgd3 
00233     do ia1=-Na1,Na1
00234      do ia2=-Na2,Na2
00235       do ia3=-Na3,Na3
00236        if(abs(ia1)==Na1.and.mod(kgd1,2)==0.and.Na1/=0) WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0 
00237        if(abs(ia2)==Na2.and.mod(kgd2,2)==0.and.Na2/=0) WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0 
00238        if(abs(ia3)==Na3.and.mod(kgd3,2)==0.and.Na3/=0) WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0 
00239        SUM_REAL=SUM_REAL+WEIGHT_R(ia1,ia2,ia3)
00240       enddo!ia3
00241      enddo!ia2
00242     enddo!ia1
00243     write(6,'(a50,f15.8,i8)')'SUM_WEIGHT, ncalck:',SUM_REAL,ncalck 
00244     if(abs(SUM_REAL-dble(ncalck))>1.0d-6)then 
00245      stop 'SUM_WEIGHT/=ncalck'
00246     endif 
00247     return
00248   end subroutine get_weightR    
00249   !
00250 end module m_truncation