calc_int_sub.F

Go to the documentation of this file.
00001       subroutine search_kq(NTK,SK0,q1,q2,q3,ik,ikq,shift_G)
00002       implicit none 
00003       integer::NTK
00004       real(8)::SK0(3,NTK)  
00005       real(8)::q1,q2,q3
00006       integer::ik
00007       integer::jk
00008       real(8)::SKQ(3)
00009       integer::ikq,shift_G(3)
00010       real(8),parameter::dlt_BZ=1.0d-6!20170322  
00011 !--
00012       SKQ(1)=SK0(1,ik)+q1
00013       SKQ(2)=SK0(2,ik)+q2 
00014       SKQ(3)=SK0(3,ik)+q3 
00015 !--
00016 !     if(SKQ(1)<-0.5D0)then 
00017 !      SKQ(1)=SKQ(1)+1.0D0 
00018 !      shift_G(1)=-1
00019 !      endif 
00020 !      if(SKQ(2)<-0.5D0)then 
00021 !      SKQ(2)=SKQ(2)+1.0D0 
00022 !      shift_G(2)=-1
00023 !      endif 
00024 !      if(SKQ(3)<-0.5D0)then 
00025 !      SKQ(3)=SKQ(3)+1.0D0 
00026 !      shift_G(3)=-1
00027 !      endif 
00028 !      if(SKQ(1)>0.5D0)then 
00029 !      SKQ(1)=SKQ(1)-1.0D0 
00030 !      shift_G(1)=+1
00031 !      endif 
00032 !      if(SKQ(2)>0.5D0)then 
00033 !      SKQ(2)=SKQ(2)-1.0D0 
00034 !      shift_G(2)=+1
00035 !      endif 
00036 !      if(SKQ(3)>0.5D0)then 
00037 !      SKQ(3)=SKQ(3)-1.0D0 
00038 !      shift_G(3)=+1
00039 !      endif 
00040 !--
00041       if(SKQ(1)>1.50d0+dlt_BZ)then 
00042        SKQ(1)=SKQ(1)-2.0D0 
00043        shift_G(1)=+2 
00044       endif 
00045       if(SKQ(1)>0.5D0+dlt_BZ)then 
00046        SKQ(1)=SKQ(1)-1.0D0 
00047        shift_G(1)=+1
00048       endif 
00049       if(SKQ(1)<=-1.5D0+dlt_BZ)then 
00050        SKQ(1)=SKQ(1)+2.0D0 
00051        shift_G(1)=-2 
00052       endif 
00053       if(SKQ(1)<=-0.5D0+dlt_BZ)then 
00054        SKQ(1)=SKQ(1)+1.0D0 
00055        shift_G(1)=-1
00056       endif 
00057 !--
00058       if(SKQ(2)>1.5D0+dlt_BZ)then 
00059        SKQ(2)=SKQ(2)-2.0D0 
00060        shift_G(2)=+2 
00061       endif 
00062       if(SKQ(2)>0.5D0+dlt_BZ)then 
00063        SKQ(2)=SKQ(2)-1.0D0 
00064        shift_G(2)=+1
00065       endif 
00066       if(SKQ(2)<=-1.5D0+dlt_BZ)then 
00067        SKQ(2)=SKQ(2)+2.0D0 
00068        shift_G(2)=-2 
00069       endif 
00070       if(SKQ(2)<=-0.5D0+dlt_BZ)then 
00071        SKQ(2)=SKQ(2)+1.0D0 
00072        shift_G(2)=-1
00073       endif 
00074 !--
00075       if(SKQ(3)>1.5D0+dlt_BZ)then 
00076        SKQ(3)=SKQ(3)-2.0D0 
00077        shift_G(3)=+2 
00078       endif 
00079       if(SKQ(3)>0.5D0+dlt_BZ)then 
00080        SKQ(3)=SKQ(3)-1.0D0 
00081        shift_G(3)=+1
00082       endif 
00083       if(SKQ(3)<=-1.5D0+dlt_BZ)then 
00084        SKQ(3)=SKQ(3)+2.0D0 
00085        shift_G(3)=-2 
00086       endif 
00087       if(SKQ(3)<=-0.5D0+dlt_BZ)then 
00088        SKQ(3)=SKQ(3)+1.0D0 
00089        shift_G(3)=-1
00090       endif 
00091 !--
00092       do jk=1,NTK
00093        if(ABS(SK0(1,jk)-SKQ(1))<1.D-6.and.
00094      +    ABS(SK0(2,jk)-SKQ(2))<1.D-6.and. 
00095      +    ABS(SK0(3,jk)-SKQ(3))<1.D-6)then 
00096         ikq=jk 
00097        endif 
00098       enddo 
00099       RETURN  
00100       END 
00101 !--
00102       subroutine calc_InterStateMatrix(NTK,NTG,NG0,KG0,C0_K,C0_KQ,
00103      +                                 ik,ikq,nwx2,nwy2,nwz2,
00104      +                                 nfft1,nfft2,Nl123,
00105      +                                 wfunc,fftwk,fs,LG0,NG_for_eps,
00106      +                                 shift_G,m_tmp)
00107       use fft_3d 
00108       implicit none 
00109       type(fft3_struct)::fs 
00110       integer::NTK,NTG 
00111       integer::NG0(NTK),KG0(3,NTG,NTK) 
00112       complex(8)::C0_K(NTG),C0_KQ(NTG) 
00113       integer::ik,ikq 
00114       integer::nwx2,nwy2,nwz2,nfft1,nfft2,Nl123 
00115       real(8)::wfunc(Nl123*2),fftwk(Nl123*2) 
00116       integer::NG_for_eps 
00117       integer::LG0(3,NTG) 
00118       integer::shift_G(3)
00119       integer::ig,igb1,igb2,igb3,ind,igL,igL1,igL2,igL3  
00120       integer::ir,ir1,ir2,ir3
00121       real(8)::dG1,dG2,dG3,phase 
00122       complex(8),allocatable::cell_periodic_k(:)
00123       complex(8),allocatable::cell_periodic_kq(:)
00124       complex(8)::f,pf
00125       complex(kind=8)::phx(nwx2),phy(nwy2),phz(nwz2)
00126       complex(kind=8)::m_tmp(NG_for_eps) 
00127       real(8),parameter::pi=dacos(-1.0d0)
00128       real(8),parameter::tpi=2.0d0*pi 
00129       complex(8),parameter::ci=(0.0D0,1.0D0) 
00130 !--
00131       allocate(cell_periodic_k(Nl123))
00132       allocate(cell_periodic_kq(Nl123))
00133       wfunc(:)=0.0D0
00134       fftwk(:)=0.0D0
00135       do ig=1,NG0(ik) 
00136       igb1=KG0(1,ig,ik) 
00137       igb2=KG0(2,ig,ik) 
00138       igb3=KG0(3,ig,ik) 
00139       igb1=MOD(nwx2+igb1,nwx2)+1
00140       igb2=MOD(nwy2+igb2,nwy2)+1
00141       igb3=MOD(nwz2+igb3,nwz2)+1
00142       ind=igb1+(igb2-1)*nfft1+(igb3-1)*nfft1*nfft2 
00143       wfunc(ind)=dble(C0_K(ig))
00144       wfunc(ind+Nl123)=dimag(C0_K(ig))
00145       enddo 
00146       call fft3_bw(fs,wfunc,fftwk) 
00147       do ir=1,Nl123 
00148       cell_periodic_k(ir)=cmplx(wfunc(ir),wfunc(ir+Nl123)) 
00149       enddo 
00150 !--
00151       wfunc(:)=0.0D0
00152       fftwk(:)=0.0D0
00153       do ig=1,NG0(ikq)
00154       igb1=KG0(1,ig,ikq)
00155       igb2=KG0(2,ig,ikq) 
00156       igb3=KG0(3,ig,ikq) 
00157       igb1=MOD(nwx2+igb1,nwx2)+1
00158       igb2=MOD(nwy2+igb2,nwy2)+1      
00159       igb3=MOD(nwz2+igb3,nwz2)+1 
00160       ind=igb1+(igb2-1)*nfft1+(igb3-1)*nfft1*nfft2 
00161       wfunc(ind)=dble(C0_KQ(ig))
00162       wfunc(ind+Nl123)=dimag(C0_KQ(ig)) 
00163       enddo 
00164       call fft3_bw(fs,wfunc,fftwk)    
00165       do ir=1,Nl123
00166       cell_periodic_kq(ir)=cmplx(wfunc(ir),wfunc(ir+Nl123)) 
00167       enddo 
00168 !--
00169       dG1=dble(shift_G(1))
00170       dG2=dble(shift_G(2))
00171       dG3=dble(shift_G(3))
00172 !YOSHIHIDE YOSHIMOTO 20080620 
00173       do ir1=1,nwx2
00174         phase = tpi*(dble(ir1-1)*dG1/dble(nwx2))
00175         phx(ir1) = exp(-ci*phase)
00176       end do
00177       do ir2=1,nwy2
00178         phase = tpi*(dble(ir2-1)*dG2/dble(nwy2))
00179         phy(ir2) = exp(-ci*phase)
00180       end do
00181       do ir3=1,nwz2
00182         phase = tpi*(dble(ir3-1)*dG3/dble(nwz2))
00183         phz(ir3) = exp(-ci*phase)
00184       end do
00185 !--
00186       do ir3=1,nwz2
00187       do ir2=1,nwy2
00188       do ir1=1,nwx2
00189        ir=ir1+(ir2-1)*nfft1+(ir3-1)*nfft1*nfft2 
00190 !YOSHIHIDE YOSHIMOTO 20080620 
00191          pf=phx(ir1)*phy(ir2)*phz(ir3)
00192          cell_periodic_kq(ir)
00193      + = cell_periodic_kq(ir) * pf 
00194       enddo 
00195       enddo 
00196       enddo 
00197 !--
00198       wfunc(:)=0.0D0
00199       fftwk(:)=0.0D0
00200       do ir=1,Nl123 
00201       f=CONJG(cell_periodic_kq(ir))*cell_periodic_k(ir)
00202       wfunc(ir)=dble(f)
00203       wfunc(ir+Nl123)=dimag(f) 
00204       enddo 
00205       call fft3_fw(fs,wfunc,fftwk) 
00206       do igL=1,NG_for_eps
00207       igL1=-LG0(1,igL)
00208       igL2=-LG0(2,igL)
00209       igL3=-LG0(3,igL) 
00210       igL1=MOD(nwx2+igL1,nwx2)+1
00211       igL2=MOD(nwy2+igL2,nwy2)+1
00212       igL3=MOD(nwz2+igL3,nwz2)+1
00213       ind=igL1+(igL2-1)*nfft1+(igL3-1)*nfft1*nfft2
00214       m_tmp(igL)=cmplx(wfunc(ind),wfunc(ind+Nl123))
00215       enddo 
00216 !--
00217       deallocate(cell_periodic_k)
00218       deallocate(cell_periodic_kq)
00219       RETURN  
00220       END 
00221 !--
00222       subroutine make_eps(NTG,NTGQ,ne,itrs,NG,LGtmp,RWtmp,rginvtmp,
00223      + pgtmp,nnp,L1,L2,L3,packtmp,epsirr,epsmk) 
00224       implicit none 
00225       integer::NTG,NTGQ,itrs,NG,ne,L1,L2,L3,nnp 
00226       integer::LGtmp(3,NTG) 
00227       integer::RWtmp(3) 
00228       real(8)::rginvtmp(3,3) 
00229       integer::pgtmp(3) 
00230       integer::packtmp(-L1:L1,-L2:L2,-L3:L3) 
00231       complex(4)::epsirr(NTGQ,NTGQ,ne) 
00232       integer::ig,jg,i1,i2,i3,j1,j2,j3,k1,k2,k3,iig,jjg 
00233       real(8)::phase 
00234       complex(8)::pf1,pf2 
00235       complex(4)::epsmk(NTGQ,NTGQ,ne) 
00236       real(8),parameter::pi=dacos(-1.0d0)
00237       real(8),parameter::tpi=2.0d0*pi 
00238       complex(8),parameter::ci=(0.0D0,1.0D0) 
00239 !--
00240       epsmk(:,:,:)=0.0d0 
00241       select case(itrs) 
00242       case(1)!=== not time-reversal ===*      
00243       do ig=1,NG 
00244        i1=LGtmp(1,ig); j1=LGtmp(2,ig); k1=LGtmp(3,ig) 
00245        i2=i1+RWtmp(1); j2=j1+RWtmp(2); k2=k1+RWtmp(3) 
00246        i3=int(rginvtmp(1,1))*i2+int(rginvtmp(1,2))*j2 
00247      +   +int(rginvtmp(1,3))*k2 
00248        j3=int(rginvtmp(2,1))*i2+int(rginvtmp(2,2))*j2 
00249      +   +int(rginvtmp(2,3))*k2 
00250        k3=int(rginvtmp(3,1))*i2+int(rginvtmp(3,2))*j2 
00251      +   +int(rginvtmp(3,3))*k2 
00252        iig=packtmp(i3,j3,k3) 
00253        phase=tpi*(dble(i1)*dble(pgtmp(1))
00254      +           +dble(j1)*dble(pgtmp(2))
00255      +           +dble(k1)*dble(pgtmp(3))) 
00256        pf1=exp(-ci*phase/dble(nnp)) 
00257        do jg=1,NG 
00258         i1=LGtmp(1,jg); j1=LGtmp(2,jg); k1=LGtmp(3,jg) 
00259         i2=i1+RWtmp(1); j2=j1+RWtmp(2); k2=k1+RWtmp(3) 
00260         i3=int(rginvtmp(1,1))*i2+int(rginvtmp(1,2))*j2 
00261      +    +int(rginvtmp(1,3))*k2 
00262         j3=int(rginvtmp(2,1))*i2+int(rginvtmp(2,2))*j2 
00263      +    +int(rginvtmp(2,3))*k2 
00264         k3=int(rginvtmp(3,1))*i2+int(rginvtmp(3,2))*j2 
00265      +    +int(rginvtmp(3,3))*k2 
00266         jjg=packtmp(i3,j3,k3) 
00267         phase=tpi*(dble(i1)*dble(pgtmp(1))
00268      +            +dble(j1)*dble(pgtmp(2))
00269      +            +dble(k1)*dble(pgtmp(3))) 
00270         pf2=exp(ci*phase/dble(nnp)) 
00271 !       write(6,'(a,x,2i5)') 'iig jjg',iig,jjg  
00272         epsmk(ig,jg,:)=epsirr(iig,jjg,:)*pf1*pf2 
00273        enddo!jg 
00274       enddo!ig 
00275       case(-1)!=== time-reversal ===*      
00276       do ig=1,NG 
00277        i1=-LGtmp(1,ig); j1=-LGtmp(2,ig); k1=-LGtmp(3,ig) 
00278        i2=i1+RWtmp(1); j2=j1+RWtmp(2); k2=k1+RWtmp(3) 
00279        i3=int(rginvtmp(1,1))*i2+int(rginvtmp(1,2))*j2 
00280      +   +int(rginvtmp(1,3))*k2 
00281        j3=int(rginvtmp(2,1))*i2+int(rginvtmp(2,2))*j2 
00282      +   +int(rginvtmp(2,3))*k2 
00283        k3=int(rginvtmp(3,1))*i2+int(rginvtmp(3,2))*j2 
00284      +   +int(rginvtmp(3,3))*k2 
00285        iig=packtmp(i3,j3,k3) 
00286        phase=tpi*(dble(i1)*dble(pgtmp(1))
00287      +           +dble(j1)*dble(pgtmp(2))
00288      +           +dble(k1)*dble(pgtmp(3))) 
00289        pf1=exp(-ci*phase/dble(nnp)) 
00290        do jg=1,NG 
00291         i1=-LGtmp(1,jg); j1=-LGtmp(2,jg); k1=-LGtmp(3,jg) 
00292         i2=i1+RWtmp(1); j2=j1+RWtmp(2); k2=k1+RWtmp(3) 
00293         i3=int(rginvtmp(1,1))*i2+int(rginvtmp(1,2))*j2 
00294      +    +int(rginvtmp(1,3))*k2 
00295         j3=int(rginvtmp(2,1))*i2+int(rginvtmp(2,2))*j2 
00296      +    +int(rginvtmp(2,3))*k2 
00297         k3=int(rginvtmp(3,1))*i2+int(rginvtmp(3,2))*j2 
00298      +    +int(rginvtmp(3,3))*k2 
00299         jjg=packtmp(i3,j3,k3) 
00300         phase=tpi*(dble(i1)*dble(pgtmp(1))
00301      +            +dble(j1)*dble(pgtmp(2))
00302      +            +dble(k1)*dble(pgtmp(3))) 
00303         pf2=exp(ci*phase/dble(nnp)) 
00304         epsmk(ig,jg,:)=epsirr(jjg,iig,:)*pf1*pf2 
00305        enddo!jg 
00306       enddo!ig 
00307       end select 
00308 !     write(6,*)'finish make_eps'
00309 !--
00310       return 
00311       end 
00312 !--
00313       SUBROUTINE OUTER_PRODUCT(vec_x,vec_y,vec_z)
00314       implicit none 
00315       real(8)::vec_x(3),vec_y(3),vec_z(3) 
00316       vec_z(1)=vec_x(2)*vec_y(3)-vec_x(3)*vec_y(2)
00317       vec_z(2)=vec_x(3)*vec_y(1)-vec_x(1)*vec_y(3) 
00318       vec_z(3)=vec_x(1)*vec_y(2)-vec_x(2)*vec_y(1)
00319       RETURN
00320       END
00321 !--
00322       subroutine search_q(SQI,SQ)
00323       implicit none 
00324       real(8)::SQI(3),SQ(3)
00325       real(8),parameter::dlt_BZ=1.0d-6!20170322  
00326       SQ=SQI 
00327       if(SQ(1)>  1.50d0+dlt_BZ) SQ(1)=SQ(1)-2.0d0 
00328       if(SQ(1)>  0.50d0+dlt_BZ) SQ(1)=SQ(1)-1.0d0 
00329       if(SQ(1)<=-1.50d0+dlt_BZ) SQ(1)=SQ(1)+2.0d0 
00330       if(SQ(1)<=-0.50d0+dlt_BZ) SQ(1)=SQ(1)+1.0d0 
00331       if(SQ(2)>  1.50d0+dlt_BZ) SQ(2)=SQ(2)-2.0d0 
00332       if(SQ(2)>  0.50d0+dlt_BZ) SQ(2)=SQ(2)-1.0d0 
00333       if(SQ(2)<=-1.50d0+dlt_BZ) SQ(2)=SQ(2)+2.0d0 
00334       if(SQ(2)<=-0.50d0+dlt_BZ) SQ(2)=SQ(2)+1.0d0 
00335       if(SQ(3)>  1.50d0+dlt_BZ) SQ(3)=SQ(3)-2.0d0 
00336       if(SQ(3)>  0.50d0+dlt_BZ) SQ(3)=SQ(3)-1.0d0 
00337       if(SQ(3)<=-1.50d0+dlt_BZ) SQ(3)=SQ(3)+2.0d0 
00338       if(SQ(3)<=-0.50d0+dlt_BZ) SQ(3)=SQ(3)+1.0d0 
00339       RETURN  
00340       END 
00341 !--
00342       subroutine kcheck(ktmp,RWtmp) 
00343       implicit none 
00344       real(8),intent(inout)::ktmp(3)
00345       integer,intent(out)::RWtmp(3) 
00346       real(8),parameter::dlt_BZ=1.0d-6!20170321 
00347 !--
00348 !      write(6,'(a8,3f15.10)')'before',ktmp(1),ktmp(2),ktmp(3)
00349 !--
00350       if(ktmp(1)>1.50d0+dlt_BZ)then 
00351        ktmp(1)=ktmp(1)-2.0d0;RWtmp(1)=-2
00352       endif 
00353       if(ktmp(1)>0.50d0+dlt_BZ)then 
00354        ktmp(1)=ktmp(1)-1.0d0;RWtmp(1)=-1
00355       endif 
00356       if(ktmp(1)<=-1.50d0+dlt_BZ)then 
00357        ktmp(1)=ktmp(1)+2.0d0;RWtmp(1)=2 
00358       endif 
00359       if(ktmp(1)<=-0.50d0+dlt_BZ)then 
00360        ktmp(1)=ktmp(1)+1.0d0;RWtmp(1)=1
00361       endif 
00362 !---
00363       if(ktmp(2)>1.50d0+dlt_BZ)then 
00364        ktmp(2)=ktmp(2)-2.0d0;RWtmp(2)=-2 
00365       endif 
00366       if(ktmp(2)>0.50d0+dlt_BZ)then 
00367        ktmp(2)=ktmp(2)-1.0d0;RWtmp(2)=-1
00368       endif 
00369       if(ktmp(2)<=-1.50d0+dlt_BZ)then 
00370        ktmp(2)=ktmp(2)+2.0d0;RWtmp(2)=2 
00371       endif 
00372       if(ktmp(2)<=-0.50d0+dlt_BZ)then 
00373        ktmp(2)=ktmp(2)+1.0d0;RWtmp(2)=1
00374       endif 
00375 !--
00376       if(ktmp(3)>1.50d0+dlt_BZ)then 
00377        ktmp(3)=ktmp(3)-2.0d0;RWtmp(3)=-2 
00378       endif 
00379       if(ktmp(3)>0.50d0+dlt_BZ)then 
00380        ktmp(3)=ktmp(3)-1.0d0;RWtmp(3)=-1
00381       endif 
00382       if(ktmp(3)<=-1.50d0+dlt_BZ)then 
00383        ktmp(3)=ktmp(3)+2.0d0;RWtmp(3)=2 
00384       endif 
00385       if(ktmp(3)<=-0.50d0+dlt_BZ)then 
00386        ktmp(3)=ktmp(3)+1.0d0;RWtmp(3)=1
00387       endif 
00388 !--
00389 !     write(6,'(a8,3f15.10)')'after',ktmp(1),ktmp(2),ktmp(3)
00390 !--
00391       return
00392       end 
00393 !--
00394 !20170322  
00395       subroutine kcheck_trs(ktmp,RWtmp)!20170322 
00396       implicit none 
00397       real(8),intent(inout)::ktmp(3)
00398       integer,intent(out)::RWtmp(3) 
00399       real(8),parameter::dlt_BZ=-1.0d-6 
00400 !--
00401 !      write(6,'(a8,3f15.10)') 'before',ktmp(1),ktmp(2),ktmp(3)
00402 !--
00403       if(ktmp(1)>=1.50d0+dlt_BZ)then 
00404        ktmp(1)=ktmp(1)-2.0d0;RWtmp(1)=-2
00405       endif 
00406       if(ktmp(1)>=0.50d0+dlt_BZ)then 
00407        ktmp(1)=ktmp(1)-1.0d0;RWtmp(1)=-1
00408       endif 
00409       if(ktmp(1)<-1.50d0+dlt_BZ)then 
00410        ktmp(1)=ktmp(1)+2.0d0;RWtmp(1)=2 
00411       endif 
00412       if(ktmp(1)<-0.50d0+dlt_BZ)then 
00413        ktmp(1)=ktmp(1)+1.0d0;RWtmp(1)=1
00414       endif 
00415 !--
00416       if(ktmp(2)>=1.50d0+dlt_BZ)then 
00417        ktmp(2)=ktmp(2)-2.0d0;RWtmp(2)=-2 
00418       endif 
00419       if(ktmp(2)>=0.50d0+dlt_BZ)then 
00420        ktmp(2)=ktmp(2)-1.0d0;RWtmp(2)=-1
00421       endif 
00422       if(ktmp(2)<-1.50d0+dlt_BZ)then 
00423        ktmp(2)=ktmp(2)+2.0d0;RWtmp(2)=2 
00424       endif 
00425       if(ktmp(2)<-0.50d0+dlt_BZ)then 
00426        ktmp(2)=ktmp(2)+1.0d0;RWtmp(2)=1
00427       endif 
00428 !--
00429       if(ktmp(3)>=1.50d0+dlt_BZ)then 
00430        ktmp(3)=ktmp(3)-2.0d0;RWtmp(3)=-2 
00431       endif 
00432       if(ktmp(3)>=0.50d0+dlt_BZ)then 
00433        ktmp(3)=ktmp(3)-1.0d0;RWtmp(3)=-1
00434       endif 
00435       if(ktmp(3)<-1.50d0+dlt_BZ)then 
00436        ktmp(3)=ktmp(3)+2.0d0;RWtmp(3)=2 
00437       endif 
00438       if(ktmp(3)<-0.50d0+dlt_BZ)then 
00439        ktmp(3)=ktmp(3)+1.0d0;RWtmp(3)=1
00440       endif 
00441 !--
00442 !      write(6,'(a8,3f15.10)') 'after',ktmp(1),ktmp(2),ktmp(3)
00443 !--
00444       return
00445       end 
00446 !--
00447       subroutine make_LG0(NTG,b1,b2,b3,Gcut_for_eps,Gcut_for_psi, 
00448      + q1,q2,q3,LG0,NG_for_eps,NG_for_psi)
00449       implicit none 
00450       integer::NTG,igL,igL1,igL2,igL3
00451       integer::NG_for_eps,NG_for_psi            
00452       integer::LG0(3,NTG)    
00453       real(8)::b1(3),b2(3),b3(3) 
00454       real(8)::Gcut_for_eps,Gcut_for_psi  
00455       real(8)::q1,q2,q3,qgL2,qgL(3)
00456       integer,parameter::NGL1=150!100!300 
00457       integer,parameter::NGL2=150!100!300 
00458       integer,parameter::NGL3=150!100!300 
00459 !--
00460       igL=0
00461       do igL1=-NGL1,NGL1 
00462       do igL2=-NGL2,NGL2 
00463       do igL3=-NGL3,NGL3 
00464       qgL(:)=(q1+dble(igL1))*b1(:)
00465      +      +(q2+dble(igL2))*b2(:)
00466      +      +(q3+dble(igL3))*b3(:)    
00467       qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00468       if(qgL2 <= Gcut_for_eps) then 
00469        igL=igL+1 
00470        LG0(1,igL)=igL1
00471        LG0(2,igL)=igL2
00472        LG0(3,igL)=igL3 
00473 !      write(6,*) igL
00474       endif  
00475       enddo 
00476       enddo 
00477       enddo 
00478       NG_for_eps=igL 
00479 !--
00480       do igL1=-NGL1,NGL1 
00481       do igL2=-NGL2,NGL2 
00482       do igL3=-NGL3,NGL3 
00483       qgL(:)=(q1+dble(igL1))*b1(:)
00484      +      +(q2+dble(igL2))*b2(:)
00485      +      +(q3+dble(igL3))*b3(:)    
00486       qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00487       if(qgL2>Gcut_for_eps.and.qgL2<=Gcut_for_psi) then 
00488        igL=igL+1 
00489        LG0(1,igL)=igL1
00490        LG0(2,igL)=igL2
00491        LG0(3,igL)=igL3 
00492 !      write(6,*) igL
00493       endif  
00494       enddo 
00495       enddo 
00496       enddo 
00497       NG_for_psi=igL 
00498 !--
00499 !Do igL=1,NG_for_eps 
00500 !write(6,*) igL,LG0(:,igL)
00501 !enddo 
00502 !--
00503       RETURN 
00504       END 
00505 !--
00506       subroutine make_KG0(NTG,b1,b2,b3,Gcut,q1,q2,q3,KG0,NG) 
00507       implicit none 
00508       integer,intent(in)::NTG
00509       real(8),intent(in)::b1(3),b2(3),b3(3) 
00510       real(8),intent(in)::Gcut 
00511       real(8),intent(in)::q1,q2,q3        
00512       integer,intent(out)::KG0(3,NTG)    
00513       integer,intent(out)::NG 
00514       integer::igL,igL1,igL2,igL3
00515       real(8)::qgL(3),qgL2  
00516       integer,parameter::NGL1=150!100!300
00517       integer,parameter::NGL2=150!100!300 
00518       integer,parameter::NGL3=150!100!300
00519       igL=0
00520       do igL1=-NGL1,NGL1 
00521       do igL2=-NGL2,NGL2 
00522       do igL3=-NGL3,NGL3 
00523       qgL(:)=(q1+dble(igL1))*b1(:)+(q2+dble(igL2))*b2(:)
00524      +      +(q3+dble(igL3))*b3(:)    
00525       qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00526       if(qgL2<=Gcut) then 
00527        igL=igL+1 
00528        KG0(1,igL)=igL1
00529        KG0(2,igL)=igL2
00530        KG0(3,igL)=igL3 
00531       endif  
00532       enddo 
00533       enddo 
00534       enddo 
00535       NG=igL 
00536 !--
00537 !write(6,*) "maxabs KG0_1=",maxval(abs(KG0(1,:))) 
00538 !write(6,*) "maxabs KG0_2=",maxval(abs(KG0(2,:))) 
00539 !write(6,*) "maxabs KG0_3=",maxval(abs(KG0(3,:))) 
00540 !--
00541 !do igL=1,NG 
00542 ! write(6,*) igL,KG0(:,igL)
00543 !enddo 
00544 !--
00545       RETURN 
00546       END