calc_intW.F

Go to the documentation of this file.
00001       PROGRAM MAIN
00002 !--
00003       use m_rdinput       
00004       use fft_3d 
00005       include "config_w3d.h" 
00006       call read_input 
00007 !--
00008 !OPEN(117,R,FILE='dat.bandcalc') 
00009       OPEN(117,FILE='./dir-wfn/dat.bandcalc') 
00010       rewind(117) 
00011       read(117,*) Ecut_for_psi 
00012       read(117,*) FermiEnergy  
00013       read(117,*) Etot
00014       write(6,*)'Ecut_for_psi=',Ecut_for_psi 
00015       write(6,*)'FermiEnergy=',FermiEnergy  
00016       write(6,*)'Etot=',Etot
00017 !--
00018 !OPEN(105,R,FILE='dat.lattice')
00019       OPEN(105,FILE='./dir-wfn/dat.lattice') 
00020       REWIND(105)
00021       READ(105,*) a1(1),a1(2),a1(3)!a1 vector
00022       READ(105,*) a2(1),a2(2),a2(3)!a2 vector
00023       READ(105,*) a3(1),a3(2),a3(3)!a3 vector
00024       write(6,*)'FINISH REDING a1,a2,a3'
00025       CLOSE(105)
00026       call OUTER_PRODUCT(a2(1),a3(1),b1(1))
00027       VOLUME=a1(1)*b1(1)+a1(2)*b1(2)+a1(3)*b1(3)
00028       b1(:)=b1(:)*tpi/VOLUME 
00029       call OUTER_PRODUCT(a3(1),a1(1),b2(1))
00030       b2(:)=b2(:)*tpi/VOLUME 
00031       call OUTER_PRODUCT(a1(1),a2(1),b3(1))
00032       b3(:)=b3(:)*tpi/VOLUME 
00033       write(6,*) 
00034       write(6,*)'======================'
00035       write(6,*)'CHECK OF (A_VEC*B_VEC)'
00036       write(6,*)'======================'
00037       write(6,*) 
00038       write(6,*)'a1*b1=',b1(1)*a1(1)+b1(2)*a1(2)+b1(3)*a1(3)
00039       write(6,*)'a2*b2=',b2(1)*a2(1)+b2(2)*a2(2)+b2(3)*a2(3)
00040       write(6,*)'a3*b3=',b3(1)*a3(1)+b3(2)*a3(2)+b3(3)*a3(3)
00041       write(6,*)'VOLUME OF UNIT CELL=',VOLUME  
00042 !--
00043 !OPEN(100,R,FILE='dat.symmetry') 
00044       OPEN(100,FILE='./dir-wfn/dat.symmetry') 
00045       rewind(100) 
00046       read(100,*) nsymq 
00047       read(100,*) nnp
00048       allocate(rg(3,3,nsymq));rg=0
00049       allocate(pg(3,nsymq));pg=0
00050       allocate(rginv(3,3,nsymq));rginv(:,:,:)=0.0d0 
00051       do iop=1,nsymq
00052        read(100,*)((rg(i,j,iop),i=1,3),j=1,3) 
00053        read(100,*)(pg(i,iop),i=1,3)   
00054       enddo 
00055       CLOSE(100) 
00056       rginv=rg 
00057       do iop=1,nsymq
00058       call invmat(3,rginv(1,1,iop)) 
00059       enddo 
00060       write(6,*)'finish rg'
00061 !check 
00062       do iop=1,nsymq
00063        write(6,*) iop
00064        do i=1,3
00065         write(6,'(3I5,1x,3F15.10)')
00066      +  (rg(i,j,iop),j=1,3),(rginv(i,j,iop),j=1,3)
00067        enddo 
00068       enddo 
00069 !do iop=1,nsymq
00070 !do i=1,3
00071 !do j=1,3
00072 !s=0.0d0 
00073 !do k=1,3
00074 !s=s+rg(i,k,iop)*rginv(k,j,iop)
00075 !enddo 
00076 !write(6,*) i,j,s
00077 !enddo 
00078 !enddo 
00079 !enddo 
00080 !--
00081 !OPEN(101,R,FILE='dat.sample-k') 
00082       OPEN(101,FILE='./dir-wfn/dat.sample-k') 
00083       rewind(101) 
00084       read(101,*) Nk_irr 
00085       allocate(SKI(3,Nk_irr));SKI(:,:)=0.0D0 
00086       do ik=1,Nk_irr 
00087        read(101,*)(SKI(i,IK),i=1,3) 
00088       enddo 
00089       close(101) 
00090       write(6,*)'finish ski'
00091 !20170327 
00092       !call est_NTK(NK_irr,SKI(1,1),NTK,nkb1,nkb2,nkb3)  
00093       !write(6,'(a24,4i10)')'nkb1,nkb2,nkb3,NTK=',nkb1,nkb2,nkb3,NTK  
00094       !Na1=nkb1/2;Na2=nkb2/2;Na3=nkb3/2
00095       !write(6,'(a24,3i10)')'Na1,Na2,Na3=',Na1,Na2,Na3 
00096 !--
00097 !OPEN(132,R,FILE='dat.nkm')20170419  
00098       OPEN(132,FILE='./dir-wfn/dat.nkm') 
00099       allocate(NGI(Nk_irr));NGI(:)=0
00100 !     allocate(NBI(Nk_irr));NBI(:)=0!20170419
00101       rewind(132)
00102       do ik=1,Nk_irr 
00103        read(132,*) NGI(ik)!,NBI(ik)!20170419
00104       enddo 
00105       close(132) 
00106       NTG=maxval(abs(NGI(:))) 
00107 !     NTB=maxval(abs(NBI(:)))!20170419  
00108       write(6,'(a10,i10)')'NTG=',NTG  
00109 !     write(6,'(a10,i10)')'NTB=',NTB!20170419
00110 !--
00111 !OPEN(104,R,FILE='dat.kg') 
00112       OPEN(104,FILE='./dir-wfn/dat.kg') 
00113       rewind(104) 
00114       allocate(KGI(3,NTG,Nk_irr));KGI(:,:,:)=0 
00115       do ik=1,Nk_irr 
00116        read(104,*) NG_for_psi 
00117        NGI(IK)=NG_for_psi 
00118        do ig=1,NG_for_psi 
00119         read(104,*)(KGI(i,ig,ik),i=1,3) 
00120        enddo 
00121       enddo  
00122       close(104)
00123       write(6,*)'finish kgi'
00124 !--
00125       allocate(LKGI(NTG,Nk_irr));LKGI=0.0d0 
00126       do ik=1,Nk_irr 
00127        do ig=1,NGI(ik) 
00128        ktmp(1)=(SKI(1,ik)+dble(KGI(1,ig,ik)))*b1(1)
00129      +        +(SKI(2,ik)+dble(KGI(2,ig,ik)))*b2(1) 
00130      +        +(SKI(3,ik)+dble(KGI(3,ig,ik)))*b3(1) 
00131 
00132        ktmp(2)=(SKI(1,ik)+dble(KGI(1,ig,ik)))*b1(2)
00133      +        +(SKI(2,ik)+dble(KGI(2,ig,ik)))*b2(2) 
00134      +        +(SKI(3,ik)+dble(KGI(3,ig,ik)))*b3(2) 
00135 
00136        ktmp(3)=(SKI(1,ik)+dble(KGI(1,ig,ik)))*b1(3)
00137      +        +(SKI(2,ik)+dble(KGI(2,ig,ik)))*b2(3) 
00138      +        +(SKI(3,ik)+dble(KGI(3,ig,ik)))*b3(3) 
00139        LKGI(ig,ik)=ktmp(1)**2+ktmp(2)**2+ktmp(3)**2
00140        enddo!ig 
00141        write(6,*) maxval(LKGI(:,ik))
00142       enddo!ik 
00143       write(6,*) 
00144       write(6,*) maxval(LKGI(:,:))
00145       Ecut_for_psi=maxval(LKGI(:,:))+1.0d-8
00146       write(6,*)'Ecut_for_psi=',Ecut_for_psi 
00147 !--
00148       L1=maxval(abs(KGI(1,:,:)))+1;write(6,*)'L1=',L1 
00149       L2=maxval(abs(KGI(2,:,:)))+1;write(6,*)'L2=',L2 
00150       L3=maxval(abs(KGI(3,:,:)))+1;write(6,*)'L3=',L3 
00151       allocate(packing(-L1:L1,-L2:L2,-L3:L3,Nk_irr))   
00152       packing(:,:,:,:)=0 
00153       do ik=1,Nk_irr 
00154        do ig=1,NGI(ik) 
00155         i1=KGI(1,ig,ik);j1=KGI(2,ig,ik);k1=KGI(3,ig,ik) 
00156         packing(i1,j1,k1,ik)=ig 
00157        enddo 
00158       enddo 
00159 !--
00160 !fft grid
00161       !20180309 
00162       !m1=maxval(abs(KGI(1,:,:))) 
00163       !m2=maxval(abs(KGI(2,:,:))) 
00164       !m3=maxval(abs(KGI(3,:,:))) 
00165       !nwx2=algn235(2*m1) 
00166       !nwy2=algn235(2*m2) 
00167       !nwz2=algn235(2*m3) 
00168       !
00169       htmp=dsqrt(dot_product(a1,a1))
00170       h1(:)=a1(:)/htmp 
00171       htmp=dsqrt(dot_product(a2,a2))
00172       h2(:)=a2(:)/htmp 
00173       htmp=dsqrt(dot_product(a3,a3)) 
00174       h3(:)=a3(:)/htmp 
00175       d1=abs(dot_product(b1,h1)) 
00176       d2=abs(dot_product(b2,h2)) 
00177       d3=abs(dot_product(b3,h3)) 
00178       qwf=2.0d0*dsqrt(Ecut_for_psi) 
00179       !
00180       !nwx2=algn235(int(qwf/d1)+1,1) 
00181       !nwy2=algn235(int(qwf/d2)+1,1) 
00182       !nwz2=algn235(int(qwf/d3)+1,1) 
00183       !
00184       nwx2=algn235(int(qwf/d1)+1)!20180412 by Y.Nomura 
00185       nwy2=algn235(int(qwf/d2)+1)!20180412 by Y.Nomura 
00186       nwz2=algn235(int(qwf/d3)+1)!20180412 by Y.Nomura 
00187       !
00188       nfft1=nwx2+1
00189       nfft2=nwy2+1
00190       nfft3=nwz2+1
00191       Nl123=nfft1*nfft2*nfft3 
00192       write(6,*)'nwx2=',nwx2 
00193       write(6,*)'nwy2=',nwy2 
00194       write(6,*)'nwz2=',nwz2 
00195       write(6,*)'nfft1=',nfft1 
00196       write(6,*)'nfft2=',nfft2 
00197       write(6,*)'nfft3=',nfft3 
00198       write(6,*)'NL123=',Nl123  
00199       call fft3_init(nwx2,nwy2,nwz2,nfft1,nfft2,nfft3,fs) 
00200 !--
00201 !20180302  
00202       call est_NTK(Nk_irr,Nsymq,SKI(1,1),rg(1,1,1),NTK)
00203       write(6,*)'Estimated NTK=',NTK 
00204 !--
00205 !gen(SK0,numirr,numrot,trs,RW) 
00206       allocate(SK0(3,NTK));SK0(:,:)=0.0d0
00207       allocate(numirr(NTK));numirr(:)=0
00208       allocate(numrot(NTK));numrot(:)=0
00209       allocate(trs(NTK));trs(:)=0
00210       allocate(RW(3,NTK));RW(:,:)=0
00211 !20161207 
00212 !      do ik=1,Nk_irr 
00213 !       SK0(:,ik)=SKI(:,ik) 
00214 !       numirr(ik)=ik;numrot(ik)=1;trs(ik)=1;RW(1:3,ik)=0
00215 !      enddo 
00216 !      jk=Nk_irr 
00217 !--
00218       jk=0
00219       do ik=1,Nk_irr
00220       do iop=1,Nsymq
00221 !sym
00222        ktmp(:)=0.0d0;RWtmp(:)=0  
00223        ktmp(1)=rg(1,1,iop)*SKI(1,ik)
00224      +        +rg(1,2,iop)*SKI(2,ik)
00225      +        +rg(1,3,iop)*SKI(3,ik)
00226        ktmp(2)=rg(2,1,iop)*SKI(1,ik)
00227      +        +rg(2,2,iop)*SKI(2,ik)
00228      +        +rg(2,3,iop)*SKI(3,ik)
00229        ktmp(3)=rg(3,1,iop)*SKI(1,ik)
00230      +        +rg(3,2,iop)*SKI(2,ik)
00231      +        +rg(3,3,iop)*SKI(3,ik)
00232        call kcheck(ktmp(1),RWtmp(1))!rewind check 
00233        do iik=1,jk
00234         if(abs(SK0(1,iik)-ktmp(1))<1.0d-4.and. 
00235      +     abs(SK0(2,iik)-ktmp(2))<1.0d-4.and. 
00236      +     abs(SK0(3,iik)-ktmp(3))<1.0d-4)goto 1000
00237        enddo!iik
00238        jk=jk+1
00239        SK0(:,jk)=ktmp(:)
00240        numirr(jk)=ik;numrot(jk)=iop;trs(jk)=1;RW(:,jk)=RWtmp(:)
00241 !time-reversal
00242 1000   ktmp(:)=0.0d0;RWtmp(:)=0  
00243        ktmp(1)=rg(1,1,iop)*SKI(1,ik)
00244      +        +rg(1,2,iop)*SKI(2,ik)
00245      +        +rg(1,3,iop)*SKI(3,ik)
00246        ktmp(2)=rg(2,1,iop)*SKI(1,ik)
00247      +        +rg(2,2,iop)*SKI(2,ik)
00248      +        +rg(2,3,iop)*SKI(3,ik) 
00249        ktmp(3)=rg(3,1,iop)*SKI(1,ik)
00250      +        +rg(3,2,iop)*SKI(2,ik)
00251      +        +rg(3,3,iop)*SKI(3,ik) 
00252        call kcheck_trs(ktmp(1),RWtmp(1))!rewind check 20170321 
00253        do iik=1,jk
00254         if(abs(SK0(1,iik)-(-ktmp(1)))<1.0d-4.and. 
00255      +     abs(SK0(2,iik)-(-ktmp(2)))<1.0d-4.and. 
00256      +     abs(SK0(3,iik)-(-ktmp(3)))<1.0d-4)goto 2000
00257        enddo!iik
00258        jk=jk+1
00259        SK0(:,jk)=-ktmp(:)
00260        numirr(jk)=ik;numrot(jk)=iop;trs(jk)=-1;RW(:,jk)=RWtmp(:)
00261 2000  enddo 
00262       enddo 
00263 !--
00264       if(NTK/=jk)then 
00265        write(6,*)'ERROR;STOP;NTK should be jk'   
00266        write(6,*)'NTK=',NTK,'jk=',jk;STOP
00267       endif 
00268 !--   
00269        write(6,*) 
00270        write(6,*)'====='
00271        write(6,*)' SK0 '
00272        write(6,*)'====='
00273        write(6,*) 
00274        do ik=1,NTK
00275         write(6,'(I5,3F12.5)') ik,SK0(:,ik) 
00276        enddo 
00277 !--
00278       call est_nkbi(NTK,SK0(1,1),nkb1,nkb2,nkb3)  
00279       write(6,*) 
00280       write(6,'(a24,4i10)')'nkb1,nkb2,nkb3,NTK=',nkb1,nkb2,nkb3,NTK  
00281       Na1=nkb1/2;Na2=nkb2/2;Na3=nkb3/2
00282       write(6,'(a24,3i10)')'Na1,Na2,Na3=',Na1,Na2,Na3 
00283 !--
00284 !gen(NG0,KG0) 
00285       allocate(KGtmp(3,NTG));KGtmp(:,:)=0 
00286       allocate(KG0(3,NTG,NTK));KG0(:,:,:)=0 
00287       allocate(NG0(NTK));NG0(:)=0
00288 !20161207 
00289 !      do jk=1,Nk_irr 
00290 !       NG0(jk)=NGI(jk);KG0(:,:,jk)=KGI(:,:,jk)
00291 !      enddo 
00292 !      do jk=Nk_irr+1,NTK 
00293 !--
00294       do jk=1,NTK 
00295        if(trs(jk)==1)then 
00296         ik=numirr(jk);iop=numrot(jk) 
00297         ktmp(1)=rg(1,1,iop)*SKI(1,ik)
00298      +         +rg(1,2,iop)*SKI(2,ik)
00299      +         +rg(1,3,iop)*SKI(3,ik)+dble(RW(1,jk)) 
00300         ktmp(2)=rg(2,1,iop)*SKI(1,ik)
00301      +         +rg(2,2,iop)*SKI(2,ik)
00302      +         +rg(2,3,iop)*SKI(3,ik)+dble(RW(2,jk))  
00303         ktmp(3)=rg(3,1,iop)*SKI(1,ik)
00304      +         +rg(3,2,iop)*SKI(2,ik)
00305      +         +rg(3,3,iop)*SKI(3,ik)+dble(RW(3,jk))  
00306         call make_KG0(NTG,b1(1),b2(1),b3(1),Ecut_for_psi,
00307      +       ktmp(1),ktmp(2),ktmp(3),KG0(1,1,jk),NG_for_psi)
00308         if(NG_for_psi/=NGI(ik))then 
00309          write(6,*)'ERROR;STOP;NG_for_psi should be NGI(ik)'   
00310          write(6,*)'NG_for_psi=',NG_for_psi,'NGI(ik)=',NGI(ik);STOP
00311         endif 
00312         NG0(jk)=NG_for_psi  
00313        elseif(trs(jk)==-1)then  
00314         ik=numirr(jk); iop=numrot(jk) 
00315         ktmp(1)=rg(1,1,iop)*SKI(1,ik)
00316      +         +rg(1,2,iop)*SKI(2,ik)
00317      +         +rg(1,3,iop)*SKI(3,ik)+dble(RW(1,jk)) 
00318         ktmp(2)=rg(2,1,iop)*SKI(1,ik)
00319      +         +rg(2,2,iop)*SKI(2,ik)
00320      +         +rg(2,3,iop)*SKI(3,ik)+dble(RW(2,jk))  
00321         ktmp(3)=rg(3,1,iop)*SKI(1,ik)
00322      +         +rg(3,2,iop)*SKI(2,ik)
00323      +         +rg(3,3,iop)*SKI(3,ik)+dble(RW(3,jk))  
00324         KGtmp(:,:)=0 
00325         call make_KG0(NTG,b1(1),b2(1),b3(1),Ecut_for_psi,
00326      +       ktmp(1),ktmp(2),ktmp(3),KGtmp(1,1),NG_for_psi)
00327         if(NG_for_psi/=NGI(ik))then 
00328          write(6,*)'ERROR;STOP;NG_for_psi should be NGI(ik)'   
00329          write(6,*)'NG_for_psi=',NG_for_psi,'NGI(ik)=',NGI(ik);STOP
00330         endif 
00331         NG0(jk)=NG_for_psi  
00332         KG0(:,:,jk)=-KGtmp(:,:) 
00333        endif 
00334       enddo 
00335 !--
00336       write(6,*) 
00337       write(6,*)'====='
00338       write(6,*)' NG0 '
00339       write(6,*)'====='
00340       write(6,*) 
00341       do ik=1,NTK 
00342        write(6,'(1I10)') NG0(ik) 
00343       enddo 
00344 !--
00345 !OPEN(113,R,FILE='dat.wan',FORM='unformatted') 
00346       OPEN(113,FILE='./dir-wan/dat.wan',FORM='unformatted') 
00347       REWIND(113)       
00348       read(113) NWF 
00349       allocate(C0(NTG,NWF,NTK));C0(:,:,:)=0.0D0 
00350       do ik=1,NTK
00351        read(113)((C0(ig,iw,ik),ig=1,NG0(ik)),iw=1,NWF)           
00352        !
00353        !20180519 
00354        !do iw=1,NWF
00355        ! read(113)(C0(ig,iw,ik),ig=1,NG0(ik))
00356        !enddo!iw
00357        !
00358       enddo!ik
00359       write(6,*)'FINISH REDING C0'
00360       CLOSE(113)  
00361 !--
00362 !      do ik=1,Nk_irr 
00363 !       do iw=1,NWF
00364 !        do jw=1,NWF 
00365 !         SUM_CMPX=0.0d0 
00366 !         do ig=1,NG0(ik) 
00367 !          SUM_CMPX=SUM_CMPX+CONJG(C0(ig,iw,ik))*C0(ig,jw,ik) 
00368 !         enddo 
00369 !         write(6,'(3i5,x,2f15.10)') ik,iw,jw,SUM_CMPX 
00370 !        enddo 
00371 !       enddo 
00372 !      enddo 
00373 !--
00374 !OPEN(134,R,FILE='dat.wan-center') 
00375       OPEN(134,FILE='./dir-wan/dat.wan-center') 
00376       allocate(coord(3,NWF));coord(:,:)=0.0d0    
00377       rewind(134)
00378       read(134,'(a)') dum_char
00379       read(134,'(a)') dum_char
00380       do iw=1,NWF
00381        read(134,*)(coord(i,iw),i=1,3) 
00382       enddo 
00383 !--
00384 !OPEN(127,R,FILE='dat.sample-q')
00385 !      OPEN(127,FILE='./dir-wfn/dat.sample-q')
00386 !      Nq_irr=Nk_irr 
00387 !      allocate(SQR(3,Nq_irr));SQR(:,:)=0.0D0 
00388 !      rewind(127) 
00389 !      read(127,*) idum 
00390 !      do iq=1,Nq_irr 
00391 !       read(127,*)(SQR(i,iq),i=1,3) 
00392 !      enddo 
00393 !      write(6,*)'FINISH REDING SQR'
00394 !      CLOSE(127) 
00395 !20161207 
00396 !shift[-1/2:1/2]  
00397 !      allocate(SQI(3,Nq_irr));SQI(:,:)=0.0D0 
00398 !      do iq=1,Nq_irr
00399 !       call search_q(SQR(1,iq),SQI(1,iq)) 
00400 !      enddo  
00401 !      write(6,*) 
00402 !      write(6,*)'====================='
00403 !      write(6,*)'shifted SQ (SQR->SQI)'
00404 !      write(6,*)'====================='
00405 !      write(6,*) 
00406 !      do iq=1,Nq_irr 
00407 !       write(6,'(a10,3f10.5,2x,3f10.5)')'SQR->SQI',SQR(:,iq),SQI(:,iq) 
00408 !      enddo 
00409 !--
00410 !OPEN(300,R,FILE='dat.chi_cutoff')!20170402  
00411        OPEN(300,FILE='./dir-eps/dat.chi_cutoff')!20170402 
00412        read(300,*) Ecut_for_eps!20170402  
00413 !--
00414 !OPEN(135,R,FILE='W-grid') 
00415       OPEN(135,FILE='./dir-eps/dat.wgrid') 
00416       rewind(135) 
00417       read(135,*) Num_freq_grid!20170402
00418       ne=Num_freq_grid!20170402
00419       allocate(em(ne));em(:)=0.0d0 
00420       do ie=1,ne 
00421        !read(135,*) em(ie) 
00422         read(135,'(2f15.10)') em(ie)!20171212
00423       enddo 
00424 !--
00425 !OPEN(400,R,FILE='dat.log.400')!20170420  
00426 !OPEN(301,R,FILE='dat.sq')!20170420   
00427       Nq_irr=Nk_irr 
00428       allocate(SQI(3,Nq_irr));SQI(:,:)=0.0D0 
00429       ierr=CHDIR("./dir-eps") 
00430       do iq=1,Nq_irr 
00431        write(dirname,"('q',i3.3)")iq 
00432        ierr=CHDIR(dirname) 
00433        inquire(file='dat.log.400',exist=file_e) 
00434        if(file_e)then 
00435         write(6,*)'dat.log.400 exists in ',trim(dirname) 
00436        else
00437         write(6,*)'error: no dat.log.400 in ',trim(dirname) 
00438        endif 
00439        OPEN(301,FILE='dat.sq') 
00440        rewind(301)
00441        inquire(file='dat.sq',exist=file_e) 
00442        if(file_e)then 
00443         read(301,*)(SQI(i,iq),i=1,3) 
00444        else
00445         write(6,*)'error: no dat.sq in',trim(dirname) 
00446        endif 
00447        close(301) 
00448        ierr=CHDIR("..") 
00449       enddo!iq 
00450       ierr=CHDIR("..") 
00451       call system('pwd') 
00452 !--
00453 !gen(SQ,numirrq,numrotq,trsq,RWq) 
00454       NTQ=NTK!TOTAL NUMBER OF Q POINTS 
00455       allocate(SQ(3,NTQ));SQ(:,:)=0.0d0 
00456       allocate(numirrq(NTQ));numirrq(:)=0 
00457       allocate(numrotq(NTQ));numrotq(:)=0 
00458       allocate(trsq(NTQ));trsq(:)=0
00459       allocate(RWq(3,NTQ));RWq(:,:)=0      
00460 !--
00461       do iq=1,Nq_irr 
00462        SQ(:,iq)=SQI(:,iq) 
00463        numirrq(iq)=iq;numrotq(iq)=1;trsq(iq)=1;RWq(1:3,iq)=0
00464       enddo 
00465       jk=Nq_irr 
00466       do iq=1,Nq_irr
00467       do iop=1,Nsymq
00468 !sym
00469        ktmp(:)=0.0d0;RWtmp(:)=0  
00470        ktmp(1)=rg(1,1,iop)*SQI(1,iq)
00471      +        +rg(1,2,iop)*SQI(2,iq)
00472      +        +rg(1,3,iop)*SQI(3,iq)
00473        ktmp(2)=rg(2,1,iop)*SQI(1,iq)
00474      +        +rg(2,2,iop)*SQI(2,iq)
00475      +        +rg(2,3,iop)*SQI(3,iq)
00476        ktmp(3)=rg(3,1,iop)*SQI(1,iq)
00477      +        +rg(3,2,iop)*SQI(2,iq)
00478      +        +rg(3,3,iop)*SQI(3,iq)
00479        call kcheck(ktmp(1),RWtmp(1))!rewind check 
00480        do iik=1,jk 
00481         if(abs(SQ(1,iik)-ktmp(1))<1.0d-4.and. 
00482      +     abs(SQ(2,iik)-ktmp(2))<1.0d-4.and. 
00483      +     abs(SQ(3,iik)-ktmp(3))<1.0d-4)goto 1100
00484        enddo!iik
00485        jk=jk+1
00486        SQ(:,jk)=ktmp(:)
00487        numirrq(jk)=iq;numrotq(jk)=iop;trsq(jk)=1;RWq(:,jk)=RWtmp(:)
00488 !time-reversal
00489 1100   ktmp(:)=0.0d0;RWtmp(:)=0  
00490        ktmp(1)=rg(1,1,iop)*SQI(1,iq)
00491      +        +rg(1,2,iop)*SQI(2,iq)
00492      +        +rg(1,3,iop)*SQI(3,iq)
00493        ktmp(2)=rg(2,1,iop)*SQI(1,iq)
00494      +        +rg(2,2,iop)*SQI(2,iq)
00495      +        +rg(2,3,iop)*SQI(3,iq) 
00496        ktmp(3)=rg(3,1,iop)*SQI(1,iq)
00497      +        +rg(3,2,iop)*SQI(2,iq)
00498      +        +rg(3,3,iop)*SQI(3,iq) 
00499        call kcheck_trs(ktmp(1),RWtmp(1))!rewind check 20170321 
00500        do iik=1,jk
00501         if(abs(SQ(1,iik)-(-ktmp(1)))<1.0d-4.and. 
00502      +     abs(SQ(2,iik)-(-ktmp(2)))<1.0d-4.and. 
00503      +     abs(SQ(3,iik)-(-ktmp(3)))<1.0d-4)goto 2100
00504        enddo!iik
00505        jk=jk+1
00506        SQ(:,jk)=-ktmp(:)
00507        numirrq(jk)=iq;numrotq(jk)=iop;trsq(jk)=-1;RWq(:,jk)=RWtmp(:)
00508 2100  enddo 
00509       enddo 
00510 !--
00511       if(NTQ/=jk)then 
00512        write(6,*)'ERROR;STOP;NTQ should be jk'   
00513        write(6,*)'NTQ=',NTQ,'jk=',jk;STOP
00514       endif 
00515 !--
00516       write(6,*) 
00517       write(6,*)'============'
00518       write(6,*)' SQ and SK0 '
00519       write(6,*)'============'
00520       write(6,*) 
00521       do iq=1,NTQ 
00522        ik=iq
00523        write(6,'(I5,3F10.5,2x,3F10.5)') iq,SQ(:,iq),SK0(:,ik) 
00524       enddo 
00525       write(6,*) 
00526 !--
00527 !      do iq=1,NTQ 
00528 !       write(6,'(i7,x,i7,x,i7,x,i7,x,3i5)')
00529 !     + iq,numirrq(iq),numrotq(iq),trsq(iq),RWq(:,iq) 
00530 !      enddo 
00531 !      write(6,*) 
00532 !--
00533 !gen(NGQ_eps,NGQ_psi,LG0)
00534       allocate(LG0(3,NTG,NTQ));LG0(:,:,:)=0
00535       allocate(NGQ_eps(NTQ));NGQ_eps(:)=0
00536       allocate(NGQ_psi(NTQ));NGQ_psi(:)=0
00537 !--
00538       do iq=1,Nq_irr 
00539        q1=SQI(1,iq);q2=SQI(2,iq);q3=SQI(3,iq)
00540        call make_LG0(NTG,b1(1),b2(1),b3(1),Ecut_for_eps,
00541      +      Ecut_for_psi,q1,q2,q3,LG0(1,1,iq),NG_for_eps,NG_for_psi) 
00542        NGQ_eps(iq)=NG_for_eps;NGQ_psi(iq)=NG_for_psi  
00543        write(6,'(i8,3f10.5,a8,i8,a8,i10)') 
00544      + iq,q1,q2,q3,'NGeps',NG_for_eps,'NGpsi',NG_for_psi  
00545       enddo 
00546       do iq=Nq_irr+1,NTQ 
00547        if(trsq(iq)==1)then 
00548         iqir=numirrq(iq);iop=numrotq(iq) 
00549         q1=rg(1,1,iop)*SQI(1,iqir)
00550      +    +rg(1,2,iop)*SQI(2,iqir)
00551      +    +rg(1,3,iop)*SQI(3,iqir)+dble(RWq(1,iq)) 
00552         q2=rg(2,1,iop)*SQI(1,iqir)
00553      +    +rg(2,2,iop)*SQI(2,iqir)
00554      +    +rg(2,3,iop)*SQI(3,iqir)+dble(RWq(2,iq))  
00555         q3=rg(3,1,iop)*SQI(1,iqir)
00556      +    +rg(3,2,iop)*SQI(2,iqir)
00557      +    +rg(3,3,iop)*SQI(3,iqir)+dble(RWq(3,iq))  
00558         call make_LG0(NTG,b1(1),b2(1),b3(1),Ecut_for_eps,
00559      +       Ecut_for_psi,q1,q2,q3,LG0(1,1,iq),NG_for_eps,NG_for_psi) 
00560         NGQ_eps(iq)=NG_for_eps;NGQ_psi(iq)=NG_for_psi  
00561         write(6,'(i8,3f10.5,a8,i8,a8,i10)') 
00562      +  iq,q1,q2,q3,'NGeps',NG_for_eps,'NGpsi',NG_for_psi  
00563        elseif(trsq(iq)==-1)then  
00564         iqir=numirrq(iq);iop=numrotq(iq) 
00565         q1=rg(1,1,iop)*SQI(1,iqir)
00566      +    +rg(1,2,iop)*SQI(2,iqir)
00567      +    +rg(1,3,iop)*SQI(3,iqir)+dble(RWq(1,iq)) 
00568         q2=rg(2,1,iop)*SQI(1,iqir)
00569      +    +rg(2,2,iop)*SQI(2,iqir)
00570      +    +rg(2,3,iop)*SQI(3,iqir)+dble(RWq(2,iq))  
00571         q3=rg(3,1,iop)*SQI(1,iqir)
00572      +    +rg(3,2,iop)*SQI(2,iqir)
00573      +    +rg(3,3,iop)*SQI(3,iqir)+dble(RWq(3,iq))  
00574         KGtmp(:,:)=0 
00575         call make_LG0(NTG,b1(1),b2(1),b3(1),Ecut_for_eps,
00576      +       Ecut_for_psi,q1,q2,q3,KGtmp(1,1),NG_for_eps,NG_for_psi) 
00577         NGQ_eps(iq)=NG_for_eps;NGQ_psi(iq)=NG_for_psi  
00578         write(6,'(i8,3f10.5,a8,i8,a8,i10)') 
00579      +  iq,q1,q2,q3,'NGeps',NG_for_eps,'NGpsi',NG_for_psi  
00580         LG0(:,:,iq)=-KGtmp(:,:) 
00581        endif 
00582       enddo 
00583 !--
00584       Lq1=maxval(abs(LG0(1,:,:)))+1;write(6,*)'Lq1=',Lq1 
00585       Lq2=maxval(abs(LG0(2,:,:)))+1;write(6,*)'Lq2=',Lq2 
00586       Lq3=maxval(abs(LG0(3,:,:)))+1;write(6,*)'Lq3=',Lq3 
00587       allocate(packingq(-Lq1:Lq1,-Lq2:Lq2,-Lq3:Lq3,Nq_irr))   
00588       packingq(:,:,:,:)=0 
00589       do iq=1,Nq_irr 
00590        do ig=1,NGQ_eps(iq) 
00591         i1=LG0(1,ig,iq);j1=LG0(2,ig,iq);k1=LG0(3,ig,iq) 
00592         packingq(i1,j1,k1,iq)=ig 
00593         write(6,'(a,5i10)')'iq,ig,i1,j1,k1',iq,ig,i1,j1,k1
00594        enddo 
00595       enddo 
00596 !
00597       write(6,*) 
00598       do iq=1,NTQ 
00599        write(6,'(i5,3f15.10,6i5)')iq,(SQ(i1,iq),i1=1,3),
00600      + numirrq(iq),numrotq(iq),trsq(iq),(RWq(j1,iq),j1=1,3)
00601       enddo  
00602 !-- 
00603       NTGQ=maxval(NGQ_eps(:))!TOTAL NUMBER OF G VECTORS FOR EPSILON
00604       write(6,*)'NTGQ=',NTGQ  
00605 !--
00606 !OPEN(600-,R,FILE='DATA_EPSILON_fRPA',FORM='unformatted') 
00607       allocate(epstmp(NTGQ,NTGQ,ne));epstmp(:,:,:)=0.0D0!real8
00608       allocate(epstmpgm(NTGQ,NTGQ,ne,3));epstmpgm(:,:,:,:)=0.0D0!real8
00609       allocate(epsirr(NTGQ,NTGQ,ne,Nq_irr));epsirr(:,:,:,:)=0.0D0!real4 
00610 !--
00611       ierr=CHDIR("./dir-eps") 
00612       call system('pwd') 
00613       do iq=1,Nq_irr 
00614        if(abs(SKI(1,iq))<1.0d-5.and.abs(SKI(2,iq))<1.0d-5.and.
00615      +    abs(SKI(3,iq))<1.0d-5)then 
00616         write(dirname,"('q',i3.3)")iq
00617         iqgm=iq 
00618         ierr=CHDIR(dirname) 
00619         do ix=1,3 
00620          file_num=600+(ix-1)  
00621          write(filename,'("dat.epsqw.",i3.3)')file_num 
00622          open(file_num,file=filename,form='unformatted') 
00623          write(6,'(a10,a15,a5,a10)') 
00624      +   'read: ',trim(filename),'in ',trim(dirname)  
00625          rewind(file_num) 
00626          NG_for_eps=NGQ_eps(iq)
00627          read(file_num)(((epstmpgm(ig,jg,ie,ix),ig=1,NG_for_eps)
00628      +                          ,jg=1,NG_for_eps),ie=1,ne)
00629         enddo!ix 
00630         epsirr(:,:,:,iq)!epsirr:real4,epstmp:real8 
00631      + =(epstmpgm(:,:,:,1)+epstmpgm(:,:,:,2)+epstmpgm(:,:,:,3))/3.0d0
00632        else
00633         write(dirname,"('q',i3.3)")iq
00634         ierr=CHDIR(dirname) 
00635         file_num=600 
00636         write(filename,'("dat.epsqw.",i3.3)')file_num 
00637         open(file_num,file=filename,form='unformatted') 
00638         write(6,'(a10,a15,a5,a10)') 
00639      +  'read: ',trim(filename),'in ',trim(dirname)  
00640         rewind(file_num) 
00641         NG_for_eps=NGQ_eps(iq)
00642         read(file_num)(((epstmp(ig,jg,ie),ig=1,NG_for_eps)
00643      +                         ,jg=1,NG_for_eps),ie=1,ne)
00644         epsirr(:,:,:,iq)=epstmp(:,:,:)!epsirr:real4,epstmp:real8 
00645        endif!q=0 or not 
00646        ierr=CHDIR("..") 
00647       enddo!iq   
00648       ierr=CHDIR("..") 
00649       call system('pwd') 
00650 !--
00651 !      do iq=1,Nq_irr 
00652 !       NG_for_eps=NGQ_eps(iq)
00653 !       write(6,*) NG_for_eps 
00654 !       do ig=1,NG_for_eps 
00655 !        write(6,*) epsirr(ig,ig,100,iq) 
00656 !       enddo  
00657 !      enddo  
00658 !--
00659       WRITE(6,*) 
00660       WRITE(6,*)'====================================='
00661       WRITE(6,*)'=== INTERACTION CALCULATION START ==='
00662       WRITE(6,*)'====================================='
00663       WRITE(6,*) 
00664 !--
00665 !make rho(G,iw,q/=0)  
00666       allocate(rho(NTG,NWF,NTQ));rho(:,:,:)=0.0D0 
00667       do iq=1,NTQ 
00668        q1=SQ(1,iq);q2=SQ(2,iq);q3=SQ(3,iq)
00669        if(q1.ne.0.0d0.or.q2.ne.0.0d0.or.q3.ne.0.0d0) then 
00670         NG_for_eps=NGQ_eps(iq)
00671         NG_for_psi=NGQ_psi(iq)
00672         allocate(length_qg(NG_for_psi));length_qg(:)=0.0D0 
00673         do igL=1,NG_for_psi   
00674          igL1=LG0(1,igL,iq)
00675          igL2=LG0(2,igL,iq)
00676          igL3=LG0(3,igL,iq)
00677          qgL(:)=(q1+dble(igL1))*b1(:)
00678      +         +(q2+dble(igL2))*b2(:)
00679      +         +(q3+dble(igL3))*b3(:)
00680          qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00681          qgL1=dsqrt(qgL2) 
00682          length_qg(igL)=qgL1
00683         enddo 
00684 !$OMP PARALLEL PRIVATE(ik,iw,igL,shift_G,ikq,C0_K,C0_KQ,
00685 !$OMP&  fftwk,wfunc,rho_tmp,prho)
00686         allocate(prho(NTG,NWF));prho(:,:)=0.0D0 
00687         allocate(fftwk(Nl123*2));fftwk(:)=0.0d0  
00688         allocate(wfunc(Nl123*2));wfunc(:)=0.0d0 
00689         allocate(rho_tmp(NG_for_psi));rho_tmp(:)=0.0D0 
00690         allocate(C0_K(NTG));C0_K(:)=0.0d0     
00691         allocate(C0_KQ(NTG));C0_KQ(:)=0.0d0  
00692 !$OMP DO 
00693         do ik=1,NTK 
00694          do iw=1,NWF 
00695           shift_G(:)=0
00696           call search_kq(NTK,SK0(1,1),q1,q2,q3,ik,ikq,shift_G(1))
00697           C0_K(:)=C0(:,iw,ik) 
00698           C0_KQ(:)=C0(:,iw,ikq) 
00699           rho_tmp(:)=0.0D0 
00700           call calc_InterStateMatrix(NTK,NTG,NG0(1),KG0(1,1,1),
00701      +     C0_K(1),C0_KQ(1),ik,ikq,nwx2,nwy2,nwz2,
00702      +     nfft1,nfft2,Nl123,wfunc(1),fftwk(1),fs,
00703      +     LG0(1,1,iq),NG_for_psi,shift_G(1),rho_tmp(1))
00704           do igL=1,NG_for_psi  
00705            prho(igL,iw)=prho(igL,iw)+rho_tmp(igL)
00706           enddo!igL
00707          enddo!iw  
00708         enddo!ik 
00709 !$OMP END DO 
00710 !$OMP CRITICAL 
00711         rho(:,:,iq)=rho(:,:,iq)+prho(:,:) 
00712 !$OMP END CRITICAL 
00713         deallocate(prho,fftwk,wfunc,rho_tmp,C0_K,C0_KQ) 
00714 !$OMP END PARALLEL 
00715         do iw=1,NWF
00716          do igL=1,NG_for_psi 
00717           rho(igL,iw,iq)=rho(igL,iw,iq)/length_qg(igL)
00718          enddo!igL
00719         enddo!iw 
00720         deallocate(length_qg)
00721         write(6,'(a10,x,2i5)')'FINISH iq',iq
00722 !--
00723 !gen rho(G,iw,q=0); iq=NTQ 
00724        elseif(q1==0.0d0.and.q2==0.0d0.and.q3==0.0d0)then 
00725         write(6,'(a,x,3f10.5)')'q1,q2,q3',q1,q2,q3 
00726         NG_for_eps=NGQ_eps(iq);NG_for_psi=NGQ_psi(iq)
00727         allocate(length_qg(NG_for_psi));length_qg(:)=0.0d0 
00728         do igL=1,NG_for_psi  
00729          igL1=LG0(1,igL,iq)
00730          igL2=LG0(2,igL,iq)
00731          igL3=LG0(3,igL,iq)
00732          if(igL1==0.and.igL2==0.and.igL3==0)then 
00733           NoG0=igL  
00734          else
00735           qgL(:)=(q1+dble(igL1))*b1(:)
00736      +          +(q2+dble(igL2))*b2(:)
00737      +          +(q3+dble(igL3))*b3(:)
00738           qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00739           qgL1=dsqrt(qgL2) 
00740           length_qg(igL)=qgL1
00741          endif  
00742         enddo 
00743 !$OMP PARALLEL PRIVATE(ik,iw,igL,shift_G,ikq,C0_K,C0_KQ,
00744 !$OMP&  fftwk,wfunc,rho_tmp,prho)
00745         allocate(prho(NTG,NWF));prho(:,:)=0.0D0 
00746         allocate(fftwk(Nl123*2));fftwk(:)=0.0d0  
00747         allocate(wfunc(Nl123*2));wfunc(:)=0.0d0 
00748         allocate(rho_tmp(NG_for_psi));rho_tmp(:)=0.0D0 
00749         allocate(C0_K(NTG));C0_K(:)=0.0d0     
00750         allocate(C0_KQ(NTG));C0_KQ(:)=0.0d0  
00751 !$OMP DO 
00752         do ik=1,NTK 
00753          do iw=1,NWF 
00754           ikq=ik 
00755           shift_G(:)=0 
00756           C0_K(:)=C0(:,iw,ik) 
00757           C0_KQ(:)=C0(:,iw,ikq) 
00758           rho_tmp(:)=0.0D0 
00759           call calc_InterStateMatrix(NTK,NTG,NG0(1),KG0(1,1,1),
00760      +     C0_K(1),C0_KQ(1),ik,ikq,nwx2,nwy2,nwz2,
00761      +     nfft1,nfft2,Nl123,wfunc(1),fftwk(1),fs,
00762      +     LG0(1,1,iq),NG_for_psi,shift_G(1),rho_tmp(1))
00763           rho_tmp(NoG0)=0.0d0!rho=0.0 is correct 
00764           do igL=1,NG_for_psi  
00765            prho(igL,iw)=prho(igL,iw)+rho_tmp(igL)
00766           enddo!igL
00767          enddo!iw 
00768         enddo!ik  
00769 !$OMP END DO 
00770 !$OMP CRITICAL 
00771         rho(:,:,iq)=rho(:,:,iq)+prho(:,:) 
00772 !$OMP END CRITICAL 
00773         deallocate(prho,fftwk,wfunc,rho_tmp,C0_K,C0_KQ) 
00774 !$OMP END PARALLEL 
00775         do iw=1,NWF 
00776          do igL=1,NG_for_psi  
00777           if(igL/=NoG0)then 
00778            rho(igL,iw,iq)=rho(igL,iw,iq)/length_qg(igL)
00779           endif 
00780          enddo!igL 
00781         enddo!iw
00782         deallocate(length_qg)
00783         write(6,'(a10,x,2i5)')'FINISH iq',iq
00784        endif 
00785       enddo!iq 
00786       rho(:,:,:)=rho(:,:,:)/dble(NTK)
00787 !--
00788 !     do iq=1,NTQ
00789 !      do iw=1,NWF 
00790 !       do ig=1,NGQ_psi(iq) 
00791 !        write(6,'(i5,2f15.10)') ig,rho(ig,iw,iq) 
00792 !       enddo 
00793 !      enddo 
00794 !     enddo 
00795 !     stop 
00796 !--
00797 !CORRECTION TERM TO DIVERGENCE 
00798 !      ie=calc_ifreq          
00799 !      write(6,*)'ie=',ie 
00800 !--
00801        write(6,*)'HybertsenLouie'!Louie 
00802        write(6,*)'NoG0=',NoG0
00803        write(6,*)'iqgm=',iqgm 
00804        qsz=(6.0D0*(pi**2)/dble(NTQ)/dble(VOLUME))**(1.0D0/3.0D0)  
00805        write(6,*)'qsz=',qsz
00806        write(6,*)'2*qsz/pi=',(2.0d0/pi)*qsz
00807        chead=(2.0D0/pi)*qsz 
00808        write(6,*)'chead(HL) in eV=',chead*au    
00809        allocate(cheadw(ne));cheadw(:)=0.0d0 
00810        do iw=1,ne 
00811         cheadw(iw)=(2.0D0/pi)*epsirr(NoG0,NoG0,iw,iqgm)*qsz 
00812         write(6,*)'cheadw(HL) in eV=',cheadw(iw)*au   
00813         write(6,*)'1/eps=',epsirr(NoG0,NoG0,iw,Nq_irr)
00814        enddo 
00815 !      endif 
00816 !--
00817 !wing 
00818 !      allocate(cwing(NWF,NWF));cwing(:,:)=0.0D0 
00819 !      allocate(wL(NWF));wL(:)=0.0D0 
00820 !      allocate(wR(NWF));wR(:)=0.0D0 
00821 !      do iw=1,NWF 
00822 !       SUM_CMPX_1=0.0D0 
00823 !       SUM_CMPX_2=0.0D0 
00824 !       do igL=1,NGQ_eps(Nq_irr) 
00825 !        if(igL/=NoG0) then 
00826 !         SUM_CMPX_1
00827 !     +  =SUM_CMPX_1
00828 !     +  +rho(igL,iw,Nq_irr) 
00829 !     +  *epsirr(igL,NoG0,ie,Nq_irr) 
00830 !         SUM_CMPX_2
00831 !     +  =SUM_CMPX_2
00832 !     +  +CONJG(rho(igL,iw,Nq_irr))
00833 !     +  *epsirr(NoG0,igL,ie,Nq_irr) 
00834 !        Endif 
00835 !       enddo 
00836 !       wL(iw)=SUM_CMPX_1
00837 !       wR(iw)=SUM_CMPX_2 
00838 !      enddo 
00839 !      do iw=1,NWF
00840 !       write(6,'(4f15.10)') wL(iw),wR(iw)
00841 !      enddo 
00842 !      do iw=1,NWF
00843 !       do jw=1,NWF 
00844 !        cwing(iw,jw)=(qsz**2)*(wL(iw)+wR(jw))/pi 
00845 !        write(6,'(a20,2f15.10)') 'cwing(eV)=',cwing(iw,jw)*au     
00846 !       enddo 
00847 !      enddo 
00848 !      deallocate(cwing,wL,wR) 
00849 !-- 
00850 !make func(\mu,\nu,q)  
00851 !--
00852        allocate(func(NWF,NWF,NTQ));func(:,:,:)=0.0D0
00853        do iq=1,NTQ
00854         do iw=1,NWF
00855          do jw=1,NWF
00856           SUM_CMPX=0.0D0
00857           do igL=1,NGQ_psi(iq)
00858            SUM_CMPX
00859      +    =SUM_CMPX
00860      +    +rho(igL,iw,iq)
00861      +    *CONJG(rho(igL,jw,iq))
00862           enddo 
00863           func(iw,jw,iq)=SUM_CMPX 
00864          enddo 
00865         enddo 
00866         write(6,'(a,x,i5)')'FINISH iq',iq
00867        enddo!iq 
00868        allocate(funcw(NWF,NWF,ne,NTQ));funcw(:,:,:,:)=0.0D0
00869        do iq=1,NTQ
00870         allocate(epsmk(NTGQ,NTGQ,ne))
00871         allocate(packtmp(-Lq1:Lq1,-Lq2:Lq2,-Lq3:Lq3)) 
00872         epsmk(:,:,:)=0.0d0 
00873         iqir=numirrq(iq)
00874         iop=numrotq(iq)
00875         packtmp(:,:,:)=packingq(:,:,:,iqir) 
00876         call make_eps(NTG,NTGQ,ne,trsq(iq),NGQ_eps(iq),LG0(1,1,iq), 
00877      +  RWq(1,iq),rginv(1,1,iop),pg(1,iop),nnp,Lq1,Lq2,Lq3,
00878      +  packtmp(-Lq1,-Lq2,-Lq3),epsirr(1,1,1,iqir),epsmk(1,1,1)) 
00879 !--
00880 !        if(iq.eq.xxx)then 
00881 !         write(6,'(a5,3f15.10)')'SQ',SQ(1,iq),SQ(2,iq),SQ(3,iq)
00882 !         write(6,'(a5,i10)')'iqir',iqir 
00883 !         write(6,'(a5,i10)')'trsq',trsq(iq) 
00884 !         write(6,'(a5,3i10)')'RWq',RWq(1,iq),RWq(2,iq),RWq(3,iq)
00885 !         write(6,'(a5)')'epsirr' 
00886 !         do ie=1,ne 
00887 !          do igL=1,NGQ_eps(iqir)
00888 !           do jgL=1,NGQ_eps(iqir) 
00889 !            write(6,'(i5,i5,2f15.8)')igL,jgL,epsirr(igL,jgL,ie,iqir)
00890 !           enddo
00891 !          enddo 
00892 !         enddo 
00893 !         write(6,'(a5)')'epsmk' 
00894 !         do ie=1,ne 
00895 !          do igL=1,NGQ_eps(iq)
00896 !           do jgL=1,NGQ_eps(iq) 
00897 !            write(6,'(i5,i5,2f15.8)')igL,jgL,epsmk(igL,jgL,ie)
00898 !           enddo
00899 !          enddo 
00900 !         enddo 
00901 !        endif 
00902 !$OMP PARALLEL PRIVATE(ie,iw,jw,SUM_CMPX,igL,jgL)
00903 !$OMP DO 
00904         do ie=1,ne 
00905          do iw=1,NWF
00906           do jw=1,NWF
00907            SUM_CMPX=0.0D0
00908            do igL=1,NGQ_eps(iq)
00909             do jgL=1,NGQ_eps(iq) 
00910               SUM_CMPX
00911      +       =SUM_CMPX
00912      +       +rho(igL,iw,iq)
00913      +       *epsmk(igL,jgL,ie)
00914      +       *CONJG(rho(jgL,jw,iq))
00915             enddo 
00916            enddo 
00917            do igL=NGQ_eps(iq)+1,NGQ_psi(iq)
00918              SUM_CMPX
00919      +      =SUM_CMPX
00920      +      +rho(igL,iw,iq)
00921      +      *CONJG(rho(igL,jw,iq))
00922            enddo 
00923            funcw(iw,jw,ie,iq)=SUM_CMPX 
00924           enddo!jw 
00925          enddo!iw
00926         enddo!ie 
00927 !$OMP END DO 
00928 !$OMP END PARALLEL 
00929         deallocate(epsmk,packtmp)  
00930         write(6,'(a,x,i5)')'FINISH iq',iq
00931        enddo!iq 
00932 !--
00933 !write func  
00934 !      write(6,*) 
00935 !      write(6,*)'----'
00936 !      write(6,*)'func'
00937 !      write(6,*)'----'
00938 !      write(6,*) 
00939 !      Do iq=1,NTQ
00940 !       write(6,*) SQ(1,iq),SQ(2,iq),SQ(3,iq) 
00941 !       Do iw=1,NWF
00942 !        write(6,'(100F12.5)')(func(iw,jw,iq),jw=1,NWF)
00943 !       ENDDO      
00944 !       write(6,*) 
00945 !      ENDDO      
00946 !write funcw  
00947 !      write(6,*) 
00948 !      write(6,*)'-----'
00949 !      write(6,*)'funcw'
00950 !      write(6,*)'-----'
00951 !      write(6,*) 
00952 !      do iq=1,NTQ
00953 !       write(6,*) SQ(1,iq),SQ(2,iq),SQ(3,iq) 
00954 !       do iw=1,NWF
00955 !        write(6,'(100F12.5)')(funcw(iw,jw,1,iq),jw=1,NWF)
00956 !       enddo       
00957 !       write(6,*) 
00958 !      enddo 
00959 !-- 
00960 !make MATRIX_W  
00961        allocate(V_MAT(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3))  
00962        V_MAT(:,:,:,:,:)=0.0D0 
00963        do ia1=-Na1,Na1
00964         do ia2=-Na2,Na2
00965          do ia3=-Na3,Na3
00966           do iw=1,NWF
00967            do jw=1,NWF 
00968             SUM_CMPX=0.0D0 
00969             do iq=1,NTQ 
00970              q1=SQ(1,iq);q2=SQ(2,iq);q3=SQ(3,iq)
00971              phase=tpi*(q1*dble(ia1)+q2*dble(ia2)+q3*dble(ia3)) 
00972              pf=exp(-ci*phase) 
00973              SUM_CMPX=SUM_CMPX+func(iw,jw,iq)*pf
00974             enddo 
00975             V_MAT(iw,jw,ia1,ia2,ia3)
00976      +     =2.0D0*tpi*SUM_CMPX/VOLUME/DBLE(NTQ)+chead 
00977            enddo 
00978           enddo 
00979          enddo 
00980         enddo 
00981        enddo 
00982        allocate(W_MAT(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3,ne))  
00983        W_MAT(:,:,:,:,:,:)=0.0D0 
00984        do ie=1,ne
00985         do ia1=-Na1,Na1
00986          do ia2=-Na2,Na2
00987           do ia3=-Na3,Na3
00988            do iw=1,NWF
00989             do jw=1,NWF 
00990              SUM_CMPX=0.0D0 
00991              do iq=1,NTQ 
00992               q1=SQ(1,iq);q2=SQ(2,iq);q3=SQ(3,iq)
00993               phase=tpi*(q1*dble(ia1)+q2*dble(ia2)+q3*dble(ia3)) 
00994               pf=exp(-ci*phase) 
00995               SUM_CMPX=SUM_CMPX+funcw(iw,jw,ie,iq)*pf
00996              enddo 
00997              W_MAT(iw,jw,ia1,ia2,ia3,ie)
00998      +      =2.0D0*tpi*SUM_CMPX/VOLUME/DBLE(NTQ)+cheadw(ie) 
00999             enddo 
01000            enddo 
01001           enddo 
01002          enddo 
01003         enddo 
01004        enddo 
01005 !--
01006 !WRITE V_MAT  
01007        write(6,*)'# '
01008        write(6,*)'#-----'
01009        write(6,*)'#V_MAT'
01010        write(6,*)'#-----'
01011        write(6,*)'# '
01012        do ia1=-Na1,Na1 
01013         do ia2=-Na2,Na2 
01014          do ia3=-Na3,Na3 
01015           write(6,*) ia1,ia2,ia3           
01016           do iw=1,NWF
01017            write(6,'(100F12.5)')(DBLE(V_MAT(iw,jw,ia1,ia2,ia3)*au), 
01018      +                                         jw=1,NWF)
01019           enddo       
01020           write(6,*) 
01021          enddo 
01022         enddo 
01023        enddo 
01024 !--
01025 !OPEN(200,W,FILE='dat.Vmat')
01026        call system('rm -rf dir-intW') 
01027        call system('mkdir dir-intW') 
01028        ierr=CHDIR("./dir-intW") 
01029        call system("pwd") 
01030        OPEN(200,FILE='dat.Vmat') 
01031        rewind(200)
01032        write(200,'(a)')'#Bare Coulomb: V' 
01033        write(200,'(a)')'#1:R1, 2:R2, 3:R3 (lattice vector)'
01034        write(200,'(a)')'#1:i, 2:j, 3:Re(V_ij) [eV], 4:Im(V_ij) [eV]' 
01035        do ia1=-Na1,Na1
01036         do ia2=-Na2,Na2
01037          do ia3=-Na3,Na3
01038           write(200,*) ia1,ia2,ia3           
01039           do iw=1,NWF 
01040            do jw=1,NWF 
01041             write(200,'(i5,i5,2f20.10)')
01042      +      iw,jw,V_MAT(iw,jw,ia1,ia2,ia3)*au
01043            enddo!jw 
01044           enddo!iw 
01045           write(200,*)
01046          enddo 
01047         enddo 
01048        enddo 
01049        close(200) 
01050 !--
01051 !WRITE W_MAT  
01052 !--
01053        ie=calc_ifreq          
01054        write(6,*)'# '
01055        write(6,*)'#-------------------------'
01056        write(6,*)'# W_MAT(omega=calc_ifreq) '
01057        write(6,*)'#-------------------------'
01058        write(6,*)'# '
01059        do ia1=-Na1,Na1 
01060         do ia2=-Na2,Na2 
01061          do ia3=-Na3,Na3 
01062           write(6,*) ia1,ia2,ia3           
01063           do iw=1,NWF
01064            write(6,'(100F12.5)')(DBLE(W_MAT(iw,jw,ia1,ia2,ia3,ie)*au), 
01065      +                                         jw=1,NWF)
01066           enddo       
01067           write(6,*) 
01068          enddo 
01069         enddo 
01070        enddo 
01071 !--
01072 !OPEN(201,W,FILE='./dir-intW/dat.Wmat')
01073        ie=calc_ifreq          
01074        OPEN(201,FILE='dat.Wmat') 
01075        rewind(201)
01076        write(201,'(a)')'#Screened Coulomb: W at omega=calc_ifreq' 
01077        write(201,'(a)')'#1:R1, 2:R2, 3:R3 (lattice vector)'
01078        write(201,'(a)')'#1:i, 2:j, 3:Re(W_ij) [eV], 4:Im(W_ij) [eV]' 
01079        do ia1=-Na1,Na1
01080         do ia2=-Na2,Na2
01081          do ia3=-Na3,Na3
01082           write(201,*) ia1,ia2,ia3           
01083           do iw=1,NWF 
01084            do jw=1,NWF
01085             write(201,'(i5,i5,2f20.10)')
01086      +      iw,jw,W_MAT(iw,jw,ia1,ia2,ia3,ie)*au 
01087            enddo!jw 
01088           enddo!iw 
01089           write(201,*)
01090          enddo 
01091         enddo 
01092        enddo 
01093        close(201)  
01094 !--
01095 !OPEN(1000,W,FILE='./dir-intW/dat.VvsR.xxx')
01096        do iw=1,NWF 
01097         write(filename,'("dat.VvsR.",i3.3)') iw 
01098         file_num=1000+(iw-1) 
01099         OPEN(file_num,FILE=filename) 
01100         rewind(file_num) 
01101         write(file_num,'(a)')'#R dependence of bare Coulomb' 
01102         write(file_num,'(a)')
01103      +  '#1:R [AA], 2:Re[V(R)] [eV], 3:Im[V(R)] [eV]'
01104         do ia1=-Na1,Na1
01105          do ia2=-Na2,Na2
01106           do ia3=-Na3,Na3
01107            do jw=1,NWF 
01108             rij(:)=dble(ia1)*a1(:)+dble(ia2)*a2(:)+dble(ia3)*a3(:)
01109      +            +(coord(:,jw)-coord(:,iw)) 
01110             rij2=rij(1)**2+rij(2)**2+rij(3)**2
01111             rij1=dsqrt(rij2)
01112             write(file_num,'(3F12.5)')
01113      +      rij1*bohr,V_MAT(iw,jw,ia1,ia2,ia3)*au  
01114            enddo!jw 
01115           enddo!ia3 
01116          enddo!ia2 
01117         enddo!ia1 
01118         close(file_num) 
01119        enddo!iw 
01120 !--
01121 !OPEN(2000,W,FILE='./dir-intW/dat.WvsR')
01122        ie=calc_ifreq          
01123        do iw=1,NWF 
01124         write(filename,'("dat.WvsR.",i3.3)') iw 
01125         file_num=2000+(iw-1) 
01126         OPEN(file_num,FILE=filename) 
01127         rewind(file_num) 
01128         write(file_num,'(a)')
01129      +  '#R dependence of screened Coulomb at omega=calc_ifreq' 
01130         write(file_num,'(a)') 
01131      +  '#1:R [AA], 2:Re[W(R)] [eV], 3:Im[W(R)] [eV]'  
01132         do ia1=-Na1,Na1
01133          do ia2=-Na2,Na2
01134           do ia3=-Na3,Na3
01135             do jw=1,NWF 
01136              rij(:)=dble(ia1)*a1(:)+dble(ia2)*a2(:)+dble(ia3)*a3(:)
01137      +             +(coord(:,jw)-coord(:,iw)) 
01138              rij2=rij(1)**2+rij(2)**2+rij(3)**2
01139              rij1=dsqrt(rij2)
01140              write(file_num,'(3F12.5)')
01141      +       rij1*bohr,W_MAT(iw,jw,ia1,ia2,ia3,ie)*au    
01142             enddo!jw 
01143           enddo!ia3 
01144          enddo!ia2 
01145         enddo!ia1  
01146         close(file_num) 
01147        enddo!iw  
01148 !--
01149 !OPEN(3000,W,FILE='./dir-intW/dat.UvsE')
01150       do iw=1,NWF
01151        do jw=iw,NWF 
01152         write(filename,'("dat.UvsE.",i3.3,"-",i3.3)') iw,jw  
01153         file_num=3000+(iw-1)*NWF+(jw-1)  
01154         OPEN(file_num,FILE=filename) 
01155         rewind(file_num) 
01156         write(file_num,'(a)')
01157      +  '#Omega dependence of onsite Coulomb' 
01158         write(file_num,'(a)') 
01159      +  '#1:Re(w) [eV], 2:Im(w) [eV], 3:Re[U(w)] [eV], 4:Im[U(w)] [eV]'
01160         do ie=1,ne
01161          write(file_num,'(4F12.5)') em(ie)*au,W_MAT(iw,jw,0,0,0,ie)*au 
01162         enddo 
01163         close(file_num) 
01164        enddo
01165       enddo
01166       ierr=CHDIR("..") 
01167       call system("pwd") 
01168 !--
01169 !write model parameters 
01170       ierr=CHDIR("./dir-model") 
01171       call system("pwd") 
01172       ie=calc_ifreq          
01173       !wcut_model=wcut_model/au 
01174       !
01175       !old: 20180604
01176       !
01177       ! call wrt_model(nkb1,nkb2,nkb3,NTK,Na1,Na2,Na3,NWF,wcut_model,
01178       !+     W_MAT(1,1,-Na1,-Na2,-Na3,ie)) 
01179       !
01180       !new: 20181127 
01181       !
01182       !WEIGHT_R 20170317 by Y.Nomura 
01183       !
01184       allocate(WEIGHT_R(-Na1:Na1,-Na2:Na2,-Na3:Na3)); WEIGHT_R=1.0d0
01185       SUM_REAL=0.0d0 
01186       do ia1=-Na1,Na1
01187        do ia2=-Na2,Na2
01188         do ia3=-Na3,Na3
01189          if(abs(ia1)==Na1.and.mod(nkb1,2)==0.and.Na1/=0) 
01190      +    WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0 
01191          if(abs(ia2)==Na2.and.mod(nkb2,2)==0.and.Na2/=0) 
01192      +    WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0 
01193          if(abs(ia3)==Na3.and.mod(nkb3,2)==0.and.Na3/=0) 
01194      +    WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0 
01195          SUM_REAL=SUM_REAL+WEIGHT_R(ia1,ia2,ia3)
01196         enddo!ia3
01197        enddo!ia2
01198       enddo!ia1 
01199       write(6,'(a20,f15.8,i8)')'SUM_WEIGHT,NTK',SUM_REAL,NTK  
01200       if(abs(SUM_REAL-dble(NTK))>1.0d-6)then 
01201        stop 'SUM_WEIGHT/=NTK'
01202       endif 
01203       !
01204       call wrt_model(Na1,Na2,Na3,NWF,W_MAT(1,1,-Na1,-Na2,-Na3,ie),
01205      +     WEIGHT_R(-Na1,-Na2,-Na3)) 
01206       !
01207       ierr=CHDIR("..") 
01208       call system("pwd") 
01209       !
01210       deallocate(C0,epsirr) 
01211       !
01212       STOP
01213       END