chiqw_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,jk
00007       real(8)::SKQ(3)
00008       integer::ikq,shift_G(3) 
00009       real(8),parameter::dlt_BZ=1.0d-6!20170316 
00010 !--
00011       SKQ(1)=SK0(1,ik)+q1
00012       SKQ(2)=SK0(2,ik)+q2 
00013       SKQ(3)=SK0(3,ik)+q3 
00014 !--
00015       if(SKQ(1)>1.50d0+dlt_BZ)then 
00016        SKQ(1)=SKQ(1)-2.0D0 
00017        shift_G(1)=+2 
00018       endif 
00019       if(SKQ(1)>0.5D0+dlt_BZ)then 
00020        SKQ(1)=SKQ(1)-1.0D0 
00021        shift_G(1)=+1
00022       endif 
00023       if(SKQ(1)<=-1.5D0+dlt_BZ)then 
00024        SKQ(1)=SKQ(1)+2.0D0 
00025        shift_G(1)=-2 
00026       endif 
00027       if(SKQ(1)<=-0.5D0+dlt_BZ)then 
00028        SKQ(1)=SKQ(1)+1.0D0 
00029        shift_G(1)=-1
00030       endif 
00031 !--
00032       if(SKQ(2)>1.5D0+dlt_BZ)then 
00033        SKQ(2)=SKQ(2)-2.0D0 
00034        shift_G(2)=+2 
00035       endif 
00036       if(SKQ(2)>0.5D0+dlt_BZ)then 
00037        SKQ(2)=SKQ(2)-1.0D0 
00038        shift_G(2)=+1
00039       endif 
00040       if(SKQ(2)<=-1.5D0+dlt_BZ)then 
00041        SKQ(2)=SKQ(2)+2.0D0 
00042        shift_G(2)=-2 
00043       endif 
00044       if(SKQ(2)<=-0.5D0+dlt_BZ)then 
00045        SKQ(2)=SKQ(2)+1.0D0 
00046        shift_G(2)=-1
00047       endif 
00048 !--
00049       if(SKQ(3)>1.5D0+dlt_BZ)then 
00050        SKQ(3)=SKQ(3)-2.0D0 
00051        shift_G(3)=+2 
00052       endif 
00053       if(SKQ(3)>0.5D0+dlt_BZ)then 
00054        SKQ(3)=SKQ(3)-1.0D0 
00055        shift_G(3)=+1
00056       endif 
00057       if(SKQ(3)<=-1.5D0+dlt_BZ)then 
00058        SKQ(3)=SKQ(3)+2.0D0 
00059        shift_G(3)=-2 
00060       endif 
00061       if(SKQ(3)<=-0.5D0+dlt_BZ)then 
00062        SKQ(3)=SKQ(3)+1.0D0 
00063        shift_G(3)=-1
00064       endif 
00065 !--
00066       do jk=1,NTK
00067        if(ABS(SK0(1,jk)-SKQ(1))<1.D-6.and.
00068      +    ABS(SK0(2,jk)-SKQ(2))<1.D-6.and. 
00069      +    ABS(SK0(3,jk)-SKQ(3))<1.D-6)then 
00070         ikq=jk 
00071        endif 
00072 !      write(6,'(a6,3f15.7)')'SKQ',SKQ
00073       enddo 
00074       RETURN  
00075       END 
00076 !--
00077       subroutine search_list(SQI,SQ)
00078       implicit none 
00079       real(8)::SQI(3),SQ(3)
00080       real(8),parameter::dlt_BZ=1.0d-6!20170316 
00081       SQ=SQI 
00082       if(SQ(1)>  1.50d0+dlt_BZ) SQ(1)=SQ(1)-2.0d0 
00083       if(SQ(1)>  0.50d0+dlt_BZ) SQ(1)=SQ(1)-1.0d0 
00084       if(SQ(1)<=-1.50d0+dlt_BZ) SQ(1)=SQ(1)+2.0d0 
00085       if(SQ(1)<=-0.50d0+dlt_BZ) SQ(1)=SQ(1)+1.0d0 
00086       if(SQ(2)>  1.50d0+dlt_BZ) SQ(2)=SQ(2)-2.0d0 
00087       if(SQ(2)>  0.50d0+dlt_BZ) SQ(2)=SQ(2)-1.0d0 
00088       if(SQ(2)<=-1.50d0+dlt_BZ) SQ(2)=SQ(2)+2.0d0 
00089       if(SQ(2)<=-0.50d0+dlt_BZ) SQ(2)=SQ(2)+1.0d0 
00090       if(SQ(3)>  1.50d0+dlt_BZ) SQ(3)=SQ(3)-2.0d0 
00091       if(SQ(3)>  0.50d0+dlt_BZ) SQ(3)=SQ(3)-1.0d0 
00092       if(SQ(3)<=-0.50d0+dlt_BZ) SQ(3)=SQ(3)+1.0d0 
00093       if(SQ(3)<=-1.50d0+dlt_BZ) SQ(3)=SQ(3)+2.0d0 
00094       RETURN  
00095       END 
00096 !--
00097       subroutine calc_InterStateMatrix(NTK,NTG,NG0,KG0,C0_K,C0_KQ,
00098      +                                 ik,ikq,nwx2,nwy2,nwz2,
00099      +                                 nfft1,nfft2,Nl123,
00100      +                                 wfunc,fftwk,fs,LG0,NG_for_eps,
00101      +                                 shift_G,m_tmp)
00102 !--
00103       use fft_3d 
00104       implicit none 
00105       type(fft3_struct)::fs 
00106       integer::NTK,NTG 
00107       integer::NG0(NTK),KG0(3,NTG,NTK) 
00108       complex(8)::C0_K(NTG),C0_KQ(NTG) 
00109       integer::ik,ikq 
00110       integer::nwx2,nwy2,nwz2,nfft1,nfft2,Nl123 
00111       real(8)::wfunc(Nl123*2),fftwk(Nl123*2) 
00112       integer::NG_for_eps 
00113       integer::LG0(3,NTG) 
00114       integer::shift_G(3)
00115       integer::ig,igb1,igb2,igb3,ind,igL,igL1,igL2,igL3  
00116       integer::ir,ir1,ir2,ir3 
00117       complex(8)::cell_periodic_k(Nl123)  
00118       complex(8)::cell_periodic_kq(Nl123)  
00119       complex(8)::f
00120       real(8)::dG1,dG2,dG3
00121       real(8)::phase 
00122       complex(8)::pf 
00123       complex(8)::m_tmp(NG_for_eps) 
00124       real(8),parameter::pi=dacos(-1.0d0)
00125       real(8),parameter::tpi=2.0d0*pi 
00126       complex(8),parameter::ci=(0.0D0,1.0D0) 
00127 !--
00128       wfunc(:)=0.0D0
00129       fftwk(:)=0.0D0
00130       do ig=1,NG0(ik) 
00131        igb1=KG0(1,ig,ik) 
00132        igb2=KG0(2,ig,ik) 
00133        igb3=KG0(3,ig,ik) 
00134        igb1=MOD(nwx2+igb1,nwx2)+1
00135        igb2=MOD(nwy2+igb2,nwy2)+1
00136        igb3=MOD(nwz2+igb3,nwz2)+1
00137        ind=igb1+(igb2-1)*nfft1+(igb3-1)*nfft1*nfft2 
00138        wfunc(ind)=dble(C0_K(ig))
00139        wfunc(ind+Nl123)=dimag(C0_K(ig))
00140       enddo 
00141       call fft3_bw(fs,wfunc(1),fftwk(1)) 
00142       do ir=1,Nl123 
00143        cell_periodic_k(ir)=cmplx(wfunc(ir),wfunc(ir+Nl123)) 
00144       enddo 
00145 !     write(6,*) 'cell_k',cell_periodic_k  
00146 !--
00147       wfunc(:)=0.0D0
00148       fftwk(:)=0.0D0
00149       do ig=1,NG0(ikq)
00150        igb1=KG0(1,ig,ikq)
00151        igb2=KG0(2,ig,ikq) 
00152        igb3=KG0(3,ig,ikq) 
00153        igb1=MOD(nwx2+igb1,nwx2)+1
00154        igb2=MOD(nwy2+igb2,nwy2)+1      
00155        igb3=MOD(nwz2+igb3,nwz2)+1 
00156        ind=igb1+(igb2-1)*nfft1+(igb3-1)*nfft1*nfft2 
00157        wfunc(ind)=dble(C0_KQ(ig))
00158        wfunc(ind+Nl123)=dimag(C0_KQ(ig)) 
00159       enddo 
00160       call fft3_bw(fs,wfunc,fftwk)    
00161       do ir=1,Nl123
00162        cell_periodic_kq(ir)=cmplx(wfunc(ir),wfunc(ir+Nl123)) 
00163       enddo 
00164 !    write(6,*) 'cell_kq',cell_periodic_kq  
00165 !--
00166       dG1=dble(shift_G(1))
00167       dG2=dble(shift_G(2))
00168       dG3=dble(shift_G(3))
00169 !     write(6,*) 'dG=',dG1,dG2,dG3
00170 !--
00171       do ir3=1,nwz2
00172        do ir2=1,nwy2
00173         do ir1=1,nwx2
00174          ir=ir1+(ir2-1)*nfft1+(ir3-1)*nfft1*nfft2 
00175          phase=tpi*(dble(ir1-1)*dG1/dble(nwx2)
00176      +             +dble(ir2-1)*dG2/dble(nwy2) 
00177      +             +dble(ir3-1)*dG3/dble(nwz2)) 
00178          pf=exp(-ci*phase) 
00179          cell_periodic_kq(ir)
00180      +  =cell_periodic_kq(ir)*pf 
00181         enddo 
00182        enddo 
00183       enddo 
00184 !--
00185       wfunc(:)=0.0D0
00186       fftwk(:)=0.0D0
00187       do ir=1,Nl123 
00188        f=CONJG(cell_periodic_kq(ir))*cell_periodic_k(ir)
00189        wfunc(ir)=dble(f)
00190        wfunc(ir+Nl123)=dimag(f) 
00191       enddo 
00192       call fft3_fw(fs,wfunc,fftwk) 
00193       do igL=1,NG_for_eps
00194        igL1=-LG0(1,igL)
00195        igL2=-LG0(2,igL)
00196        igL3=-LG0(3,igL) 
00197        igL1=MOD(nwx2+igL1,nwx2)+1
00198        igL2=MOD(nwy2+igL2,nwy2)+1
00199        igL3=MOD(nwz2+igL3,nwz2)+1
00200        ind=igL1+(igL2-1)*nfft1+(igL3-1)*nfft1*nfft2
00201        m_tmp(igL)=cmplx(wfunc(ind),wfunc(ind+Nl123))
00202       enddo 
00203 !--
00204       RETURN  
00205       END 
00206 !--
00207       subroutine calc_VMab(NTK,NTG,NG0,KG0,C0_K,C0_KQ,ik,b1,b2,b3,vm)  
00208       implicit none 
00209       integer::NTK,NTG,ik 
00210       integer::NG0(NTK),KG0(3,NTG,NTK) 
00211       real(8)::b1(3),b2(3),b3(3) 
00212       complex(8)::C0_K(NTG),C0_KQ(NTG) 
00213       integer::ig,igb1,igb2,igb3 
00214       real(8)::G_VEC(3)
00215       complex(8)::SUM_VEC(3)
00216       complex(8)::vm(3) 
00217       complex(8),parameter::ci=(0.0D0,1.0D0) 
00218 !--
00219       SUM_VEC(:)=0.0D0 
00220       do ig=1,NG0(ik) 
00221        igb1=KG0(1,ig,ik) 
00222        igb2=KG0(2,ig,ik) 
00223        igb3=KG0(3,ig,ik) 
00224        G_VEC(:)=dble(igb1)*b1(:)+dble(igb2)*b2(:)+dble(igb3)*b3(:) 
00225        SUM_VEC(:)=SUM_VEC(:)+G_VEC(:)*CONJG(C0_KQ(ig))*C0_K(ig) 
00226       enddo 
00227 !"vm" depends on the direction of q approaching the zero limit
00228 !We employ q=[001],[010],[100] therefore
00229 !Kazuma Nakamura bugfix 2010 10 12
00230       vm(:)=SUM_VEC(:)*ci!ci is multiplied  
00231       RETURN  
00232       END 
00233 !--
00234       subroutine calc_VMaa(NTK,NTG,NG0,KG0,C0_K,C0_KQ,ik,b1,b2,b3,SKT, 
00235      +                     vm)  
00236       implicit none 
00237       integer::NTK,NTG,ik 
00238       integer::NG0(NTK),KG0(3,NTG,NTK) 
00239       real(8)::b1(3),b2(3),b3(3),SKT(3) 
00240       complex(8)::C0_K(NTG),C0_KQ(NTG) 
00241       integer::ig,igb1,igb2,igb3 
00242       real(8)::G_VEC(3)
00243       complex(8)::SUM_VEC(3)
00244       complex(8) :: vm(3) 
00245       complex(8),parameter::ci=(0.0D0,1.0D0) 
00246 !--
00247       SUM_VEC(:)=0.0D0 
00248       do ig=1,NG0(ik) 
00249        igb1=KG0(1,ig,ik) 
00250        igb2=KG0(2,ig,ik) 
00251        igb3=KG0(3,ig,ik) 
00252        G_VEC(:)=(SKT(1)+dble(igb1))*b1(:)
00253      +         +(SKT(2)+dble(igb2))*b2(:)
00254      +         +(SKT(3)+dble(igb3))*b3(:) 
00255        SUM_VEC(:)=SUM_VEC(:)+G_VEC(:)*CONJG(C0_KQ(ig))*C0_K(ig) 
00256       enddo 
00257 !"vm" depends on the direction of q approaching the zero limit 
00258 !We employ q=[001],[010],[100] therefore
00259 !Kazuma Nakamura bugfix 2010 10 12
00260       vm(:)=SUM_VEC(:)*ci!ci is multiplied  
00261       RETURN  
00262       END 
00263 !--
00264       subroutine make_C0(NTG,itrs,NG,KGtmp,RWtmp,rginvtmp,pgtmp,
00265      +           nnp,L1,L2,L3,packtmp,OCCtmp,C0_K) 
00266       implicit none 
00267       integer::NTG,L1,L2,L3,nnp  
00268       integer::itrs 
00269       integer::NG
00270       integer::KGtmp(3,NTG) 
00271       integer::RWtmp(3) 
00272       real(8)::rginvtmp(3,3) 
00273       integer::pgtmp(3) 
00274       integer::packtmp(-L1:L1,-L2:L2,-L3:L3) 
00275       complex(8)::OCCtmp(NTG) 
00276       integer::ig,jg,i1,i2,i3,j1,j2,j3,k1,k2,k3 
00277       real(8)::phase 
00278       complex(8)::pf 
00279       complex(8)::C0_K(NTG) 
00280       real(8),parameter::pi=dacos(-1.0d0)
00281       real(8),parameter::tpi=2.0d0*pi 
00282       complex(8),parameter::ci=(0.0D0,1.0D0) 
00283 !--
00284       C0_K(:)=0.0d0 
00285       if(itrs==1) then 
00286        do ig=1,NG 
00287         i1=KGtmp(1,ig); j1=KGtmp(2,ig); k1=KGtmp(3,ig) 
00288         i2=i1+RWtmp(1); j2=j1+RWtmp(2); k2=k1+RWtmp(3) 
00289         i3=int(rginvtmp(1,1))*i2+int(rginvtmp(1,2))*j2 
00290      +    +int(rginvtmp(1,3))*k2 
00291         j3=int(rginvtmp(2,1))*i2+int(rginvtmp(2,2))*j2 
00292      +    +int(rginvtmp(2,3))*k2 
00293         k3=int(rginvtmp(3,1))*i2+int(rginvtmp(3,2))*j2 
00294      +    +int(rginvtmp(3,3))*k2 
00295         jg=packtmp(i3,j3,k3) 
00296         phase=tpi*(dble(i1)*dble(pgtmp(1))
00297      +            +dble(j1)*dble(pgtmp(2))
00298      +            +dble(k1)*dble(pgtmp(3))) 
00299         pf=exp(-ci*phase/dble(nnp)) 
00300         C0_K(ig)=OCCtmp(jg)*pf 
00301        enddo!ig 
00302       elseif(itrs==-1) then  
00303        do ig=1,NG 
00304         i1=-KGtmp(1,ig); j1=-KGtmp(2,ig); k1=-KGtmp(3,ig) 
00305         i2=i1+RWtmp(1); j2=j1+RWtmp(2); k2=k1+RWtmp(3) 
00306         i3=int(rginvtmp(1,1))*i2+int(rginvtmp(1,2))*j2 
00307      +    +int(rginvtmp(1,3))*k2 
00308         j3=int(rginvtmp(2,1))*i2+int(rginvtmp(2,2))*j2 
00309      +    +int(rginvtmp(2,3))*k2 
00310         k3=int(rginvtmp(3,1))*i2+int(rginvtmp(3,2))*j2 
00311      +    +int(rginvtmp(3,3))*k2 
00312         jg=packtmp(i3,j3,k3) 
00313         phase=tpi*(dble(i1)*dble(pgtmp(1))
00314      +            +dble(j1)*dble(pgtmp(2))
00315      +            +dble(k1)*dble(pgtmp(3))) 
00316         pf=exp(-ci*phase/dble(nnp)) 
00317         C0_K(ig)=OCCtmp(jg)*pf 
00318        enddo!ig 
00319        C0_K(:)=conjg(C0_K(:))  
00320       endif 
00321       return 
00322       end 
00323 !---
00324       SUBROUTINE make_index_kpt(NTK,nkb1,nkb2,nkb3,SK0,index_kpt)      
00325       implicit none 
00326       integer::NTK,nkb1,nkb2,nkb3
00327       real(8)::SK0(3,NTK)  
00328       integer::ik,ix,iy,iz
00329       real(8)::x,y,z
00330       integer::index_kpt(nkb1,nkb2,nkb3)    
00331       ! 
00332       !if(MOD(NTK,2)/=0) then 
00333       ! write(6,*)'i am in make_index for odd'
00334       ! do ik=1,NTK 
00335       !  x=SK0(1,ik)*dble(nkb1) 
00336       !  y=SK0(2,ik)*dble(nkb2)
00337       !  z=SK0(3,ik)*dble(nkb3)  
00338       !  x=x+(dble(nkb1)-1.0d0)/2.0d0 
00339       !  y=y+(dble(nkb2)-1.0d0)/2.0d0
00340       !  z=z+(dble(nkb3)-1.0d0)/2.0d0 
00341       !  ix=idnint(x)+1
00342       !  iy=idnint(y)+1
00343       !  iz=idnint(z)+1
00344       !  index_kpt(ix,iy,iz)=ik
00345       ! enddo 
00346       !else!20170316 
00347       ! write(6,*)'i am in make_index for even'
00348       ! do ik=1,NTK 
00349       !  x=SK0(1,ik)*dble(nkb1) 
00350       !  y=SK0(2,ik)*dble(nkb2)
00351       !  z=SK0(3,ik)*dble(nkb3)  
00352       !  x=x+dble(nkb1)/2.0d0 
00353       !  y=y+dble(nkb2)/2.0d0
00354       !  z=z+dble(nkb3)/2.0d0 
00355       !  ix=idnint(x)
00356       !  iy=idnint(y)
00357       !  iz=idnint(z)
00358       !  index_kpt(ix,iy,iz)=ik
00359       ! enddo 
00360       !endif 
00361       !--
00362       !
00363       !20190520 Kazuma Nakamura
00364       !
00365       do ik=1,NTK 
00366        x=SK0(1,ik)*dble(nkb1) 
00367        y=SK0(2,ik)*dble(nkb2)
00368        z=SK0(3,ik)*dble(nkb3)  
00369        x=x+(dble(nkb1)-dble(mod(nkb1,2)))/2.0d0 
00370        y=y+(dble(nkb2)-dble(mod(nkb2,2)))/2.0d0 
00371        z=z+(dble(nkb3)-dble(mod(nkb3,2)))/2.0d0 
00372        ix=idnint(x)+mod(nkb1,2)
00373        iy=idnint(y)+mod(nkb2,2)
00374        iz=idnint(z)+mod(nkb3,2)
00375        index_kpt(ix,iy,iz)=ik
00376       enddo 
00377       ! 
00378       !--
00379       !do iz=1,nkb3 
00380       ! do iy=1,nkb2
00381       !  do ix=1,nkb1
00382       !   ik=index_kpt(ix,iy,iz)
00383       !   write(6,'(i8,3f15.10)') ik,SK0(:,ik)
00384       !  enddo
00385       ! enddo
00386       !enddo
00387       !--
00388       !
00389       RETURN  
00390       END 
00391 !---
00392       subroutine calc_eps_rpa(NTG,NG_for_eps,LG0,
00393      +                   q1,q2,q3,b1,b2,b3,
00394      +                   chi0,tpi,eps_rpa, 
00395      +                   file_num_chi,file_num_eps)  
00396       implicit none 
00397       integer::NTG,NG_for_eps  
00398       integer::LG0(3,NTG)    
00399       real(8)::q1,q2,q3,b1(3),b2(3),b3(3),tpi 
00400       integer::file_num_chi,file_num_eps 
00401       complex(8)::chi0(NG_for_eps,NG_for_eps)
00402       integer::ig,jg,igL1,igL2,igL3,jgL1,jgL2,jgL3 
00403       real(8)::qgL(3),qgL1,qgL2,qgLi,qgLj   
00404       complex(8)::f
00405       complex(8)::eps_rpa(NG_for_eps,NG_for_eps) 
00406 !
00407       if(file_num_chi==0) write(6,*)'not write chi'
00408 !Factor of 2 is needed for only LDA-polarization calculation, pointed by Yohiro Nohara, 2009 5 27 ---* 
00409       chi0(:,:)=2.0d0*chi0(:,:) 
00410 !cal_eps_rpa 
00411       write(6,*)'calc eps_rpa'
00412       write(6,*)' ' 
00413       eps_rpa(:,:)=0.0D0 
00414       do ig=1,NG_for_eps 
00415        igL1=LG0(1,ig)
00416        igL2=LG0(2,ig)
00417        igL3=LG0(3,ig)
00418        qgL(:)=(q1+dble(igL1))*b1(:)
00419      +       +(q2+dble(igL2))*b2(:)
00420      +       +(q3+dble(igL3))*b3(:)
00421        qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00422        qgLi=dsqrt(qgL2) 
00423       do jg=1,NG_for_eps   
00424        jgL1=LG0(1,jg)
00425        jgL2=LG0(2,jg)
00426        jgL3=LG0(3,jg)
00427        qgL(:)=(q1+dble(jgL1))*b1(:)
00428      +       +(q2+dble(jgL2))*b2(:)
00429      +       +(q3+dble(jgL3))*b3(:)
00430        qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00431        qgLj=dsqrt(qgL2) 
00432 !diag 
00433        if(ig==jg)then 
00434 !       write(file_num_chi,'(3F20.10)') qgLi,chi0(ig,jg)
00435         eps_rpa(ig,jg)=1.0D0-chi0(ig,jg)*2.0D0*tpi/qgLi/qgLj 
00436        endif 
00437 !off diag
00438        if(ig/=jg)then 
00439         eps_rpa(ig,jg)=-chi0(ig,jg)*2.0D0*tpi/qgLi/qgLj 
00440        endif 
00441       enddo 
00442       enddo 
00443 !---
00444       call invZGE(NG_for_eps,eps_rpa(1,1)) 
00445 !---
00446       write(6,*)'calc eps_rpa_inv'
00447       write(6,*)' ' 
00448       do ig=1,NG_for_eps 
00449        igL1=LG0(1,ig)
00450        igL2=LG0(2,ig)
00451        igL3=LG0(3,ig)
00452        qgL(:)=(q1+dble(igL1))*b1(:)
00453      +       +(q2+dble(igL2))*b2(:)
00454      +       +(q3+dble(igL3))*b3(:)
00455        qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00456        qgL1=dsqrt(qgL2) 
00457        f=eps_rpa(ig,ig)
00458 !      write(file_num_eps,'(3F20.10)') qgL1,dble(1.0D0/f),dimag(1.0d0/f)  
00459        write(file_num_eps,'(3F20.10)') qgL1,dble(f),dimag(f)  
00460       enddo 
00461 !---
00462       return 
00463       end 
00464 !---
00465       subroutine calc_eps_rpa_q_0(NTG,NG_for_eps,LG0,
00466      +                   q1,q2,q3,b1,b2,b3,
00467      +                   chi0,tpi,No_G_0,eps_rpa,
00468      +                   wd,delt,wcmplx,   
00469      +                   file_num_chi,file_num_eps)  
00470       implicit none 
00471       integer::NTG,NG_for_eps,No_G_0,file_num_chi,file_num_eps 
00472       integer::LG0(3,NTG)    
00473       real(8)::q1,q2,q3,b1(3),b2(3),b3(3),tpi,wd,delt,w  
00474       complex(8)::chi0(NG_for_eps,NG_for_eps)
00475       complex(8)::wcmplx,f
00476       integer::ig,jg,igL1,igL2,igL3,jgL1,jgL2,jgL3 
00477       real(8)::qgL(3),qgL1,qgL2,qgLi,qgLj   
00478       complex(8)::eps_rpa(NG_for_eps,NG_for_eps) 
00479 !
00480       if(file_num_chi==0) write(6,*)'not write chi'
00481 !Factor of 2 is needed for only LDA-polarization calculation, pointed by Yohiro Nohara, 2009 5 27 ---* 
00482       chi0(:,:)=2.0d0*chi0(:,:) 
00483       !
00484       !20191210 Kazuma Nakamura 
00485       !
00486       w=abs(dble(wcmplx))  
00487       if(w<1.0d-5)then 
00488        w=1.0d-5 
00489        write(6,*)' ' 
00490        write(6,'(a)')'### Since omega=0, change to omega=1.0d-5 ###' 
00491        write(6,*)' ' 
00492       endif 
00493       !
00494 !     wd=0.0d0 
00495 !calc_eps_rpa_q=0 
00496       write(6,*)'calc eps_rpa_q=0'
00497       write(6,*)' ' 
00498       eps_rpa(:,:)=0.0D0 
00499       do ig=1,NG_for_eps 
00500        igL1=LG0(1,ig)
00501        igL2=LG0(2,ig)
00502        igL3=LG0(3,ig)
00503        qgL(:)=(q1+dble(igL1))*b1(:)
00504      +       +(q2+dble(igL2))*b2(:)
00505      +       +(q3+dble(igL3))*b3(:)
00506        qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00507        qgLi=dsqrt(qgL2) 
00508       do jg=1,NG_for_eps   
00509        jgL1=LG0(1,jg)
00510        jgL2=LG0(2,jg)
00511        jgL3=LG0(3,jg)
00512        qgL(:)=(q1+dble(jgL1))*b1(:)
00513      +       +(q2+dble(jgL2))*b2(:)
00514      +       +(q3+dble(jgL3))*b3(:)
00515        qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00516        qgLj=dsqrt(qgL2) 
00517 !head 
00518        if(ig==jg.and.ig==No_G_0)then 
00519 !       write(file_num_chi,'(3F20.10)') qgLi,chi0(ig,jg)
00520         eps_rpa(ig,jg)
00521      + =1.0D0-chi0(ig,jg)*2.0D0*tpi 
00522      + +cmplx(-wd**2/(w**2+delt**2),delt*wd**2/(w**3+w*delt**2))
00523        endif 
00524 !body(diag)
00525        if(ig==jg.and.ig/=No_G_0)then 
00526 !       write(file_num_chi,'(3F20.10)') qgLi,chi0(ig,jg)
00527         eps_rpa(ig,jg)=1.0D0-chi0(ig,jg)*2.0D0*tpi/qgLi/qgLj  
00528        endif 
00529 !wing(column) 
00530        if(ig/=jg.and.ig==No_G_0)then 
00531         eps_rpa(ig,jg)=-chi0(ig,jg)*2.0D0*tpi/qgLj  
00532        endif 
00533 !wing(raw) 
00534        if(ig/=jg.and.jg==No_G_0)then 
00535         eps_rpa(ig,jg)=-chi0(ig,jg)*2.0D0*tpi/qgLi
00536        endif      
00537 !body(off diag)
00538        if(ig/=jg.and.ig/=No_G_0.and.jg/=No_G_0) then 
00539         eps_rpa(ig,jg)=-chi0(ig,jg)*2.0D0*tpi/qgLi/qgLj 
00540        endif 
00541       enddo 
00542       enddo 
00543 !---
00544 !      write(6,*)' '
00545 !      write(6,*)'eps(i,G0)'
00546 !      do ig=1,NG_for_eps
00547 !       write(6,*) ig,eps_rpa(ig,No_G_0)
00548 !      enddo 
00549 !      write(6,*)' '
00550 !      write(6,*)'eps(G0,j)'
00551 !      do ig=1,NG_for_eps
00552 !       write(6,*) ig,eps_rpa(No_G_0,ig)
00553 !      enddo 
00554 !      write(6,*)' '
00555 !---
00556 !      do ig=1,NG_for_eps
00557 !       do jg=1,NG_for_eps
00558 !        write(6,'(2i8,2f15.8)') ig,jg,eps_rpa(ig,jg)
00559 !       enddo 
00560 !      enddo 
00561 !---
00562       call invZGE(NG_for_eps,eps_rpa(1,1)) 
00563 !---
00564       write(6,*)'calc eps_rpa_inv_q=0'
00565       write(6,*)' ' 
00566       do ig=1,NG_for_eps 
00567        igL1=LG0(1,ig)
00568        igL2=LG0(2,ig)
00569        igL3=LG0(3,ig)
00570        qgL(:)=(q1+dble(igL1))*b1(:)
00571      +       +(q2+dble(igL2))*b2(:)
00572      +       +(q3+dble(igL3))*b3(:)
00573        qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00574        qgL1=dsqrt(qgL2) 
00575        f=eps_rpa(ig,ig)
00576 !     write(file_num_eps,'(3F20.10)') qgL1,dble(1.0D0/f),dimag(1.0d0/f) 
00577       write(file_num_eps,'(3F20.10)') qgL1,dble(f),dimag(f)  
00578       enddo 
00579       return 
00580       end 
00581 !--
00582       subroutine judge_FermiInside(NTK,NTB,E_EIG,FermiEnergy,
00583      + WindowInside,Tindx,E_AVE,band_max,band_min) 
00584       implicit none 
00585       integer::NTK,NTB
00586       real(8)::E_EIG(NTB,NTK)  
00587       real(8)::FermiEnergy
00588       integer::ik,ib,fs 
00589       real(8)::e,top,btm 
00590       integer::WindowInside(NTB)    
00591       integer::Tindx(NTB)    
00592       real(8)::E_AVE(NTB)    
00593       real(8)::band_max(NTB)
00594       real(8)::band_min(NTB)
00595 !--
00596       WindowInside(:)=0
00597       Tindx(:)=0
00598       do ib=1,NTB
00599        do ik=1,NTK 
00600         e=E_EIG(ib,ik) 
00601         if(ik==1)then 
00602          top=e
00603          btm=e
00604         endif 
00605         if(e<btm) btm=e
00606         if(e>top) top=e
00607        enddo!ik 
00608 !uocc
00609        if(btm>FermiEnergy) fs=2
00610 !pocc
00611        if(btm<=FermiEnergy.and.top>=FermiEnergy) then
00612         fs=1
00613         Tindx(ib)=ib 
00614        endif 
00615 !docc
00616        if(top<FermiEnergy) fs=0
00617 !
00618        WindowInside(ib)=fs 
00619        E_AVE(ib)=(top+btm)/2.0d0 
00620        band_max(ib)=top
00621        band_min(ib)=btm
00622       enddo!ib 
00623 !--
00624       RETURN 
00625       END 
00626 !--
00627       subroutine judge_WindowInside(NTK,NTB,E_EIG,EL,EU,WindowInside,
00628      + Tindx,E_AVE,band_max,band_min)
00629       implicit none 
00630       integer::NTK,NTB
00631       real(8)::E_EIG(NTB,NTK)  
00632       real(8)::EL,EU
00633       integer::ik,ib,fs 
00634       real(8)::e,top,btm 
00635       integer::WindowInside(NTB)    
00636       integer::Tindx(NTB)    
00637       real(8)::E_AVE(NTB)    
00638       real(8)::band_max(NTB)
00639       real(8)::band_min(NTB)
00640 !--
00641       WindowInside(:)=0
00642       Tindx(:)=0
00643       do ib=1,NTB
00644        do ik=1,NTK 
00645         e=E_EIG(ib,ik) 
00646         if(ik==1)then 
00647          top=e
00648          btm=e
00649         endif 
00650         if(e<btm) btm=e
00651         if(e>top) top=e
00652        enddo!ik 
00653 !case1:docc
00654        if(top<EL) fs=0
00655 !case2:docc
00656        if(EL<=top.and.top<EU.and.btm<EL) fs=0 
00657 !case3:pocc
00658        if(EL<=top.and.top<EU.and.EL<=btm.and.btm<EU)then 
00659         fs=1
00660         Tindx(ib)=ib
00661        endif   
00662 !case4:pocc
00663        if(btm<EL.and.top>EU)then   
00664         fs=1 
00665         Tindx(ib)=ib
00666        endif   
00667 !case5:uocc
00668        if(EL<btm.and.btm<EU.and.top>EU) fs=2 
00669 !case6:uocc 
00670        if(btm>EU) fs=2 
00671 !--
00672        WindowInside(ib)=fs 
00673        E_AVE(ib)=(top+btm)/2.0d0 
00674        band_max(ib)=top
00675        band_min(ib)=btm
00676       enddo!ib 
00677       RETURN 
00678       END 
00679 !--
00680       subroutine kcheck(ktmp,RWtmp) 
00681       implicit none 
00682       real(8),intent(inout)::ktmp(3)
00683       integer,intent(out)::RWtmp(3) 
00684       real(8),parameter::dlt_BZ=1.0d-6!20170316 
00685 !--
00686 !      write(6,'(a8,3f15.10)')'before',ktmp(1),ktmp(2),ktmp(3)
00687 !--
00688       if(ktmp(1)>1.50d0+dlt_BZ)then 
00689        ktmp(1)=ktmp(1)-2.0d0;RWtmp(1)=-2
00690       endif 
00691       if(ktmp(1)>0.50d0+dlt_BZ)then 
00692        ktmp(1)=ktmp(1)-1.0d0;RWtmp(1)=-1
00693       endif 
00694       if(ktmp(1)<=-1.50d0+dlt_BZ)then 
00695        ktmp(1)=ktmp(1)+2.0d0;RWtmp(1)=2 
00696       endif 
00697       if(ktmp(1)<=-0.50d0+dlt_BZ)then 
00698        ktmp(1)=ktmp(1)+1.0d0;RWtmp(1)=1
00699       endif 
00700 !--
00701       if(ktmp(2)>1.50d0+dlt_BZ)then 
00702        ktmp(2)=ktmp(2)-2.0d0;RWtmp(2)=-2 
00703       endif 
00704       if(ktmp(2)>0.50d0+dlt_BZ)then 
00705        ktmp(2)=ktmp(2)-1.0d0;RWtmp(2)=-1
00706       endif 
00707       if(ktmp(2)<=-1.50d0+dlt_BZ)then 
00708        ktmp(2)=ktmp(2)+2.0d0;RWtmp(2)=2 
00709       endif 
00710       if(ktmp(2)<=-0.50d0+dlt_BZ)then 
00711        ktmp(2)=ktmp(2)+1.0d0;RWtmp(2)=1
00712       endif 
00713 !--
00714       if(ktmp(3)>1.50d0+dlt_BZ)then 
00715        ktmp(3)=ktmp(3)-2.0d0;RWtmp(3)=-2 
00716       endif 
00717       if(ktmp(3)>0.50d0+dlt_BZ)then 
00718        ktmp(3)=ktmp(3)-1.0d0;RWtmp(3)=-1
00719       endif 
00720       if(ktmp(3)<=-1.50d0+dlt_BZ)then 
00721        ktmp(3)=ktmp(3)+2.0d0;RWtmp(3)=2 
00722       endif 
00723       if(ktmp(3)<=-0.50d0+dlt_BZ)then 
00724        ktmp(3)=ktmp(3)+1.0d0;RWtmp(3)=1
00725       endif 
00726 !--
00727 !     write(6,'(a8,3f15.10)')'after',ktmp(1),ktmp(2),ktmp(3)
00728 !--
00729       return
00730       end 
00731 !--
00732 !20170316 
00733       subroutine kcheck_trs(ktmp,RWtmp)!20170316 
00734       implicit none 
00735       real(8),intent(inout)::ktmp(3)
00736       integer,intent(out)::RWtmp(3) 
00737       real(8),parameter::dlt_BZ=-1.0d-6 
00738 !--
00739 !      write(6,'(a8,3f15.10)') 'before',ktmp(1),ktmp(2),ktmp(3)
00740 !--
00741       if(ktmp(1)>=1.50d0+dlt_BZ)then 
00742        ktmp(1)=ktmp(1)-2.0d0;RWtmp(1)=-2
00743       endif 
00744       if(ktmp(1)>=0.50d0+dlt_BZ)then 
00745        ktmp(1)=ktmp(1)-1.0d0;RWtmp(1)=-1
00746       endif 
00747       if(ktmp(1)<-1.50d0+dlt_BZ)then 
00748        ktmp(1)=ktmp(1)+2.0d0;RWtmp(1)=2 
00749       endif 
00750       if(ktmp(1)<-0.50d0+dlt_BZ)then 
00751        ktmp(1)=ktmp(1)+1.0d0;RWtmp(1)=1
00752       endif 
00753 !--
00754       if(ktmp(2)>=1.50d0+dlt_BZ)then 
00755        ktmp(2)=ktmp(2)-2.0d0;RWtmp(2)=-2 
00756       endif 
00757       if(ktmp(2)>=0.50d0+dlt_BZ)then 
00758        ktmp(2)=ktmp(2)-1.0d0;RWtmp(2)=-1
00759       endif 
00760       if(ktmp(2)<-1.50d0+dlt_BZ)then 
00761        ktmp(2)=ktmp(2)+2.0d0;RWtmp(2)=2 
00762       endif 
00763       if(ktmp(2)<-0.50d0+dlt_BZ)then 
00764        ktmp(2)=ktmp(2)+1.0d0;RWtmp(2)=1
00765       endif 
00766 !--
00767       if(ktmp(3)>=1.50d0+dlt_BZ)then 
00768        ktmp(3)=ktmp(3)-2.0d0;RWtmp(3)=-2 
00769       endif 
00770       if(ktmp(3)>=0.50d0+dlt_BZ)then 
00771        ktmp(3)=ktmp(3)-1.0d0;RWtmp(3)=-1
00772       endif 
00773       if(ktmp(3)<-1.50d0+dlt_BZ)then 
00774        ktmp(3)=ktmp(3)+2.0d0;RWtmp(3)=2 
00775       endif 
00776       if(ktmp(3)<-0.50d0+dlt_BZ)then 
00777        ktmp(3)=ktmp(3)+1.0d0;RWtmp(3)=1
00778       endif 
00779 !--
00780 !      write(6,'(a8,3f15.10)') 'after',ktmp(1),ktmp(2),ktmp(3)
00781 !--
00782       return
00783       end 
00784 !--
00785       subroutine make_KG0(NTG,b1,b2,b3,Gcut,q1,q2,q3,KG0,NG) 
00786       implicit none 
00787       integer,intent(in)::NTG
00788       real(8),intent(in)::b1(3),b2(3),b3(3) 
00789       real(8),intent(in)::q1,q2,q3,Gcut 
00790       integer::igL,igL1,igL2,igL3
00791       real(8)::qgL(3),qgL2  
00792       integer,intent(out)::KG0(3,NTG)    
00793       integer,intent(out)::NG 
00794       integer,parameter::NGL1=150!100!300
00795       integer,parameter::NGL2=150!100!300 
00796       integer,parameter::NGL3=150!100!300
00797       igL=0
00798       do igL1=-NGL1,NGL1 
00799        do igL2=-NGL2,NGL2 
00800         do igL3=-NGL3,NGL3 
00801          qgL(:)=(q1+dble(igL1))*b1(:)
00802      +         +(q2+dble(igL2))*b2(:)
00803      +         +(q3+dble(igL3))*b3(:)    
00804          qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00805          if(qgL2<=Gcut)then 
00806           igL=igL+1 
00807           KG0(1,igL)=igL1
00808           KG0(2,igL)=igL2
00809           KG0(3,igL)=igL3 
00810          endif  
00811         enddo 
00812        enddo 
00813       enddo 
00814       NG=igL 
00815 !---
00816 !write(6,*)"maxabs KG0_1=",maxval(abs(KG0(1,:))) 
00817 !write(6,*)"maxabs KG0_2=",maxval(abs(KG0(2,:))) 
00818 !write(6,*)"maxabs KG0_3=",maxval(abs(KG0(3,:))) 
00819 !---
00820 !do igL=1,NG 
00821 ! write(6,*) igL,KG0(:,igL)
00822 !enddo 
00823       RETURN 
00824       END 
00825 !---
00826       SUBROUTINE OUTER_PRODUCT(vec_x,vec_y,vec_z)
00827       implicit none 
00828       real(8)::vec_x(3),vec_y(3),vec_z(3) 
00829       vec_z(1)=vec_x(2)*vec_y(3)-vec_x(3)*vec_y(2)
00830       vec_z(2)=vec_x(3)*vec_y(1)-vec_x(1)*vec_y(3) 
00831       vec_z(3)=vec_x(1)*vec_y(2)-vec_x(2)*vec_y(1)
00832       RETURN
00833       END