calc_intJ.F

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