gw.f90

Go to the documentation of this file.
00001 PROGRAM GW
00002 !
00003 use mpi
00004 use fft_3d 
00005 use m_rd_input       
00006 use m_rd_dat_wfn
00007 use m_rd_dat_wan 
00008 use m_rd_dat_eps 
00009 include "config.h" 
00010 !
00011 !mpi 
00012 !
00013 comm=MPI_COMM_WORLD
00014 call MPI_INIT(ierr)
00015 call MPI_COMM_RANK(comm,myrank,ierr)
00016 call MPI_COMM_SIZE(comm,nproc,ierr)
00017 write(6,'(a10,i10,a10,i10)')'myrank=',myrank,'nproc=',nproc
00018 !file_id=9000+myrank 
00019 call MPI_BARRIER(comm,ierr)
00020 start_time=MPI_Wtime() 
00021 !
00022 !read input
00023 !
00024 if(myrank.eq.0)then 
00025  call read_input 
00026  gw_grid_separation=gw_grid_separation/au
00027 endif 
00028 call MPI_BARRIER(comm,ierr)
00029 !
00030 call MPI_Bcast(gw_grid_separation,1,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00031 call MPI_Bcast(Rc_range_spacing,1,MPI_INTEGER,0,comm,ierr) 
00032 call MPI_Bcast(Ncalc,1,MPI_INTEGER,0,comm,ierr) 
00033 call MPI_Bcast(calc_SC,1,MPI_LOGICAL,0,comm,ierr) 
00034 !
00035 call MPI_Bcast(N_sym_points,1,MPI_INTEGER,0,comm,ierr) 
00036 call MPI_Bcast(Ndiv,1,MPI_INTEGER,0,comm,ierr) 
00037 call MPI_Bcast(reading_sk_format,1,MPI_INTEGER,0,comm,ierr) 
00038 call MPI_BARRIER(comm,ierr)
00039 !
00040 if(myrank/=0) allocate(dense(3))
00041 if(myrank/=0) allocate(SK_sym_pts(3,N_sym_points))  
00042 call MPI_Bcast(dense,3,MPI_INTEGER,0,comm,ierr)
00043 call MPI_Bcast(SK_sym_pts,3*N_sym_points,MPI_DOUBLE_PRECISION,0,comm,ierr)
00044 call MPI_BARRIER(comm,ierr)
00045 !
00046 !read wfn 
00047 !
00048 if(myrank.eq.0)then 
00049  call rd_dat_symmetry 
00050  call rd_dat_bandcalc 
00051  call rd_dat_lattice 
00052  call rd_dat_sample_k 
00053  call rd_dat_nkm 
00054  call rd_dat_kg 
00055  call rd_dat_eigenvalue  
00056  call rd_dat_wavefunction 
00057 endif 
00058 call MPI_BARRIER(comm,ierr)
00059 !
00060 !sym(100)
00061 !
00062 call MPI_Bcast(nsymq,1,MPI_INTEGER,0,comm,ierr) 
00063 call MPI_Bcast(nnp,1,MPI_INTEGER,0,comm,ierr) 
00064 !
00065 !bandcalc(117) 
00066 !
00067 call MPI_Bcast(Ecut_for_psi,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00068 call MPI_Bcast(FermiEnergy,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00069 call MPI_Bcast(Etot,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00070 !
00071 !avec(105)
00072 !
00073 call MPI_Bcast(a1,3,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00074 call MPI_Bcast(a2,3,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00075 call MPI_Bcast(a3,3,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00076 call MPI_Bcast(b1,3,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00077 call MPI_Bcast(b2,3,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00078 call MPI_Bcast(b3,3,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00079 call MPI_Bcast(VOLUME,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00080 call MPI_Bcast(a,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00081 call MPI_Bcast(b,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00082 call MPI_Bcast(c,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00083 call MPI_Bcast(alp,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00084 call MPI_Bcast(bet,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00085 call MPI_Bcast(gmm,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00086 !
00087 !sample-k(101)
00088 !
00089 call MPI_Bcast(Nk_irr,1,MPI_INTEGER,0,comm,ierr) 
00090 call MPI_Bcast(NTK,1,MPI_INTEGER,0,comm,ierr) 
00091 call MPI_Bcast(nkb1,1,MPI_INTEGER,0,comm,ierr) 
00092 call MPI_Bcast(nkb2,1,MPI_INTEGER,0,comm,ierr) 
00093 call MPI_Bcast(nkb3,1,MPI_INTEGER,0,comm,ierr) 
00094 call MPI_Bcast(Na1,1,MPI_INTEGER,0,comm,ierr) 
00095 call MPI_Bcast(Na2,1,MPI_INTEGER,0,comm,ierr) 
00096 call MPI_Bcast(Na3,1,MPI_INTEGER,0,comm,ierr) 
00097 !
00098 !kg(104)
00099 !
00100 call MPI_Bcast(L1,1,MPI_INTEGER,0,comm,ierr) 
00101 call MPI_Bcast(L2,1,MPI_INTEGER,0,comm,ierr) 
00102 call MPI_Bcast(L3,1,MPI_INTEGER,0,comm,ierr) 
00103 !
00104 !fft 
00105 !
00106 call MPI_Bcast(nwx2,1,MPI_INTEGER,0,comm,ierr) 
00107 call MPI_Bcast(nwy2,1,MPI_INTEGER,0,comm,ierr) 
00108 call MPI_Bcast(nwz2,1,MPI_INTEGER,0,comm,ierr) 
00109 !
00110 !nkm(132)  
00111 !
00112 call MPI_Bcast(NTG,1,MPI_INTEGER,0,comm,ierr) 
00113 !
00114 !eigenvalue(111) 
00115 !
00116 call MPI_Bcast(NTB,1,MPI_INTEGER,0,comm,ierr) 
00117 !
00118 !cir(102) 
00119 !
00120 call MPI_Bcast(ncomp,1,MPI_INTEGER,0,comm,ierr) 
00121 call MPI_BARRIER(comm,ierr)
00122 !
00123 if(myrank/=0) allocate(rg(3,3,nsymq))
00124 if(myrank/=0) allocate(pg(3,nsymq))
00125 if(myrank/=0) allocate(rginv(3,3,nsymq))  
00126 if(myrank/=0) allocate(SKI(3,Nk_irr))  
00127 if(myrank/=0) allocate(SK0(3,NTK))  
00128 if(myrank/=0) allocate(numirr(NTK)) 
00129 if(myrank/=0) allocate(numrot(NTK)) 
00130 if(myrank/=0) allocate(trs(NTK)) 
00131 if(myrank/=0) allocate(RW(3,NTK))
00132 if(myrank/=0) allocate(numMK(Nk_irr))
00133 if(myrank/=0) allocate(NGI(Nk_irr)) 
00134 if(myrank/=0) allocate(NG0(NTK))
00135 if(myrank/=0) allocate(KGI(3,NTG,Nk_irr)) 
00136 if(myrank/=0) allocate(KG0(3,NTG,NTK))
00137 if(myrank/=0) allocate(packing(-L1:L1,-L2:L2,-L3:L3,Nk_irr))
00138 if(myrank/=0) allocate(E_EIGI(NTB,Nk_irr))  
00139 if(myrank/=0) allocate(CIR(NTG,NTB,Nk_irr)) 
00140 call MPI_BARRIER(comm,ierr)
00141 !
00142 !20180922
00143 !
00144 !mem_size=dble(NTG*NTB*Nk_irr)*16.0d0/1024.0d0/1024.0d0/1024.0d0 
00145  mem_size=dble(NTG*NTB*Nk_irr)*8.0d0/1024.0d0/1024.0d0/1024.0d0 
00146 !
00147 if(myrank.eq.0)then 
00148  write(6,'(a35,f20.15)')'mem size CIR (GB)',mem_size
00149 endif 
00150 ! 
00151 call MPI_Bcast(rg,3*3*nsymq,MPI_INTEGER,0,comm,ierr) 
00152 call MPI_Bcast(pg,3*nsymq,MPI_INTEGER,0,comm,ierr) 
00153 call MPI_Bcast(rginv,3*3*nsymq,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00154 call MPI_Bcast(SKI,3*Nk_irr,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00155 call MPI_Bcast(SK0,3*NTK,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00156 call MPI_Bcast(numirr,NTK,MPI_INTEGER,0,comm,ierr) 
00157 call MPI_Bcast(numrot,NTK,MPI_INTEGER,0,comm,ierr) 
00158 call MPI_Bcast(trs,NTK,MPI_INTEGER,0,comm,ierr) 
00159 call MPI_Bcast(RW,3*NTK,MPI_INTEGER,0,comm,ierr) 
00160 call MPI_Bcast(numMK,Nk_irr,MPI_INTEGER,0,comm,ierr)
00161 call MPI_Bcast(NGI,Nk_irr,MPI_INTEGER,0,comm,ierr) 
00162 call MPI_Bcast(NG0,NTK,MPI_INTEGER,0,comm,ierr) 
00163 call MPI_Bcast(KGI,3*NTG*Nk_irr,MPI_INTEGER,0,comm,ierr) 
00164 call MPI_Bcast(KG0,3*NTG*NTK,MPI_INTEGER,0,comm,ierr) 
00165 call MPI_Bcast(packing(-L1,-L2,-L3,1),(2*L1+1)*(2*L2+1)*(2*L3+1)*Nk_irr,MPI_INTEGER,0,comm,ierr) 
00166 call MPI_Bcast(E_EIGI,NTB*Nk_irr,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00167 !
00168 !20180922
00169 !
00170 !call MPI_Bcast(CIR,NTG*NTB*Nk_irr,MPI_DOUBLE_COMPLEX,0,comm,ierr) 
00171  call MPI_Bcast(CIR,NTG*NTB*Nk_irr,MPI_COMPLEX,0,comm,ierr) 
00172 !
00173 call MPI_BARRIER(comm,ierr)
00174 !
00175 !read wan 
00176 !
00177 if(myrank.eq.0)then 
00178  call rd_dat_ns_nb 
00179  call rd_dat_umat 
00180  call rd_dat_wan 
00181  call rd_dat_hmatr 
00182 endif 
00183 call MPI_BARRIER(comm,ierr)
00184 !
00185 !ns-nb(149)  
00186 !
00187 call MPI_Bcast(Mb,1,MPI_INTEGER,0,comm,ierr) 
00188 call MPI_Bcast(Mt,1,MPI_INTEGER,0,comm,ierr) 
00189 !
00190 !umat(150)  
00191 !
00192 call MPI_Bcast(NWF,1,MPI_INTEGER,0,comm,ierr) 
00193 call MPI_BARRIER(comm,ierr)
00194 !
00195 if(myrank/=0) allocate(Ns(NTK))
00196 if(myrank/=0) allocate(Nb(NTK))
00197 if(myrank/=0) allocate(Nt(NTK))
00198 if(myrank/=0) allocate(UNT(Mb,NWF,NTK))
00199 if(myrank/=0) allocate(C0_WN(NTG,NWF,NTK))
00200 if(myrank/=0) allocate(H_MAT_R(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3))
00201 call MPI_BARRIER(comm,ierr)
00202 !
00203 call MPI_Bcast(Ns,NTK,MPI_INTEGER,0,comm,ierr) 
00204 call MPI_Bcast(Nb,NTK,MPI_INTEGER,0,comm,ierr) 
00205 call MPI_Bcast(Nt,NTK,MPI_INTEGER,0,comm,ierr) 
00206 call MPI_Bcast(UNT,Mb*NWF*NTK,MPI_DOUBLE_COMPLEX,0,comm,ierr) 
00207 call MPI_Bcast(C0_WN,NTG*NWF*NTK,MPI_DOUBLE_COMPLEX,0,comm,ierr) 
00208 call MPI_Bcast(H_MAT_R(1,1,-Na1,-Na2,-Na3),NWF*NWF*(2*Na1+1)*(2*Na2+1)*(2*Na3+1),MPI_DOUBLE_COMPLEX,0,comm,ierr) 
00209 call MPI_BARRIER(comm,ierr)
00210 !
00211 !read eps 
00212 !
00213 if(myrank.eq.0)then 
00214  call rd_dat_chi_cutoff
00215  call rd_dat_ttrhdrn 
00216  call rd_dat_wgrid 
00217  call rd_dat_sq 
00218  call rd_dat_eps
00219  ! 
00220  !espinv=epsinv-1 
00221  !
00222  do iq=1,Nq_irr 
00223   NG_for_eps=NGQ_eps(iq) 
00224   do ie=1,ne 
00225    do igL=1,NG_for_eps 
00226     epsirr(igL,igL,ie,iq)=epsirr(igL,igL,ie,iq)-1.0d0    
00227    enddo 
00228   enddo 
00229  enddo 
00230  !
00231 endif!myrank.eq.0 
00232 !
00233 !chi_cutoff(300)  
00234 !
00235 call MPI_Bcast(Ecut_for_eps,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00236 !
00237 !ttrhdrn(302) 
00238 !
00239 call MPI_Bcast(idlt,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00240 call MPI_Bcast(dmna,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00241 call MPI_Bcast(dmnr,1,MPI_DOUBLE_PRECISION,0,comm,ierr)
00242 !
00243 !wgrid(135)  
00244 !
00245 call MPI_Bcast(Num_freq_grid,1,MPI_INTEGER,0,comm,ierr) 
00246 call MPI_Bcast(ne,1,MPI_INTEGER,0,comm,ierr) 
00247 !
00248 !sq(301)  
00249 !
00250 call MPI_Bcast(Nq_irr,1,MPI_INTEGER,0,comm,ierr) 
00251 call MPI_Bcast(NTQ,1,MPI_INTEGER,0,comm,ierr) 
00252 call MPI_Bcast(NTGQ,1,MPI_INTEGER,0,comm,ierr) 
00253 call MPI_Bcast(Lq1,1,MPI_INTEGER,0,comm,ierr) 
00254 call MPI_Bcast(Lq2,1,MPI_INTEGER,0,comm,ierr) 
00255 call MPI_Bcast(Lq3,1,MPI_INTEGER,0,comm,ierr) 
00256 call MPI_BARRIER(comm,ierr)
00257 !
00258 if(myrank/=0) allocate(em(ne))
00259 if(myrank/=0) allocate(pole_of_chi(ne))
00260 if(myrank/=0) allocate(mat_b(ne,ne))
00261 if(myrank/=0) allocate(SQI(3,Nq_irr)) 
00262 if(myrank/=0) allocate(SQ(3,NTQ)) 
00263 if(myrank/=0) allocate(numirrq(NTQ)) 
00264 if(myrank/=0) allocate(numrotq(NTQ))  
00265 if(myrank/=0) allocate(trsq(NTQ))  
00266 if(myrank/=0) allocate(RWq(3,NTQ))
00267 if(myrank/=0) allocate(LG0(3,NTG,NTQ))    
00268 if(myrank/=0) allocate(NGQ_eps(NTQ))
00269 if(myrank/=0) allocate(NGQ_psi(NTQ)) 
00270 if(myrank/=0) allocate(packingq(-Lq1:Lq1,-Lq2:Lq2,-Lq3:Lq3,Nq_irr))
00271 if(myrank/=0) allocate(epsirr(NTGQ,NTGQ,ne,Nq_irr))!real4 
00272 call MPI_BARRIER(comm,ierr)
00273 !
00274 mem_size=dble(NTGQ*NTGQ*ne*Nq_irr)*8.0d0/1024.0d0/1024.0d0/1024.0d0 
00275 if(myrank.eq.0)then 
00276  write(6,'(a35,f20.15)')'mem size epsirr (GB)',mem_size
00277 endif 
00278 ! 
00279 call MPI_Bcast(em,ne,MPI_DOUBLE_COMPLEX,0,comm,ierr) 
00280 call MPI_Bcast(pole_of_chi,ne,MPI_DOUBLE_COMPLEX,0,comm,ierr)
00281 call MPI_Bcast(mat_b,ne*ne,MPI_DOUBLE_COMPLEX,0,comm,ierr)
00282 call MPI_Bcast(SQI,3*Nq_irr,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00283 call MPI_Bcast(SQ,3*NTQ,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00284 call MPI_Bcast(numirrq,NTQ,MPI_INTEGER,0,comm,ierr) 
00285 call MPI_Bcast(numrotq,NTQ,MPI_INTEGER,0,comm,ierr) 
00286 call MPI_Bcast(trsq,NTQ,MPI_INTEGER,0,comm,ierr) 
00287 call MPI_Bcast(RWq,3*NTQ,MPI_INTEGER,0,comm,ierr)
00288 call MPI_Bcast(LG0,3*NTG*NTQ,MPI_INTEGER,0,comm,ierr) 
00289 call MPI_Bcast(NGQ_eps,NTQ,MPI_INTEGER,0,comm,ierr) 
00290 call MPI_Bcast(NGQ_psi,NTQ,MPI_INTEGER,0,comm,ierr) 
00291 call MPI_Bcast(packingq(-Lq1,-Lq2,-Lq3,1),(2*Lq1+1)*(2*Lq2+1)*(2*Lq3+1)*Nq_irr,MPI_INTEGER,0,comm,ierr) 
00292 call MPI_Bcast(epsirr,NTGQ*NTGQ*ne*Nq_irr,MPI_COMPLEX,0,comm,ierr)
00293 call MPI_BARRIER(comm,ierr)
00294 !
00295 !if(myrank.eq.0)then 
00296 ! write(6,*) 
00297 ! write(6,*)'====================================='
00298 ! write(6,*)'dir-eps: FREQUENCY GRID IN EPSQW (eV)'
00299 ! write(6,*)'====================================='
00300 ! write(6,*) 
00301 ! do ie=1,ne 
00302 !  write(6,'(2f15.7)') em(ie)*au 
00303 ! enddo 
00304 ! !write(6,*) 
00305 ! !write(6,*)'==========================='
00306 ! !write(6,*)'dir-eps: POLE OF MODEL GRID'
00307 ! !write(6,*)'==========================='
00308 ! !write(6,*) 
00309 ! !do ie=1,ne 
00310 ! ! write(6,'(2f15.10)') pole_of_chi(ie)
00311 ! !enddo 
00312 !endif!myrank.eq.0  
00313 !
00314 !Rc, idlt, Ncalc
00315 !
00316 avec_length(1)=dsqrt(a1(1)**2+a1(2)**2+a1(3)**2)
00317 avec_length(2)=dsqrt(a2(1)**2+a2(2)**2+a2(3)**2)
00318 avec_length(3)=dsqrt(a3(1)**2+a3(2)**2+a3(3)**2)
00319 Rc=Rc_range_spacing*maxval(avec_length) 
00320 !
00321 !Ncalc is set to NTB (defalut)
00322 !
00323 if(Ncalc.eq.0) Ncalc=NTB 
00324 !
00325 if(myrank.eq.0)then 
00326  write(6,*) 
00327  write(6,'(a40,x,i10)')'Ncalc (considered bands in self-energy)=',Ncalc 
00328  write(6,'(a40,x,f15.10)')'idlt (eV)=',idlt*au 
00329  write(6,'(a40,x,f15.10)')'Rc (AA)=',Rc*bohr   
00330 endif 
00331 !
00332 !fbk 
00333 !
00334 allocate(fbk(NTB,NTK));fbk(:,:)=0.0d0 
00335 do ik=1,NTK
00336  ikir=numirr(ik) 
00337  do ib=1,NTB
00338   FACTOR=(FermiEnergy-E_EIGI(ib,ikir))/idlt    
00339   fbk(ib,ik)=atan(FACTOR)/pi+0.5d0 
00340  enddo 
00341 enddo 
00342 SUM_REAL=0.0d0 
00343 do ik=1,NTK
00344  ikir=numirr(ik) 
00345  do ib=1,NTB
00346   SUM_REAL=SUM_REAL+fbk(ib,ik) 
00347   en=(E_EIGI(ib,ikir)-FermiEnergy)*au  
00348  enddo
00349 enddo
00350 if(myrank.eq.0)then
00351  write(6,*) 
00352  write(6,'(a40,x,f15.7)')'Nele=',2.0d0*SUM_REAL/dble(NTK)  
00353 endif 
00354 !
00355 !make sgmw
00356 !
00357 if(gw_grid_separation<idlt) gw_grid_separation=idlt
00358 !
00359 if(myrank.eq.0)then 
00360  bandmin=minval(E_EIGI(1+minval(Ns),:)) 
00361  bandmax=maxval(E_EIGI(minval(Nb)+minval(Ns),:))
00362  diff_band_energy=bandmax-bandmin 
00363  chiqw_grd_size=dble(em(ne)) 
00364  emax=bandmax+diff_band_energy 
00365  emin=bandmin-diff_band_energy 
00366  ecmax=bandmax+diff_band_energy+chiqw_grd_size 
00367  ecmin=bandmin-diff_band_energy-chiqw_grd_size  
00368  !
00369  call estimate_nsgm(ecmin,emin,emax,ecmax,gw_grid_separation,nproc,nsgm)
00370  allocate(sgmw(nsgm));sgmw=0.0d0 
00371  call make_sgmw(ecmin,emin,emax,gw_grid_separation,nsgm,sgmw(1))
00372  !
00373  write(6,*) 
00374  write(6,'(a40,x,i10)')'number of chiqw_grd =',Num_freq_grid 
00375  write(6,'(a40,x,f15.7)')'bandmax (eV)=',bandmax*au
00376  write(6,'(a40,x,f15.7)')'bandmin (eV)=',bandmin*au
00377  write(6,'(a40,x,f15.7)')'diff_band_energy (eV)=',diff_band_energy*au  
00378  write(6,'(a40,x,f15.7)')'chiqw_grd_range (eV)=',chiqw_grd_size*au  
00379  write(6,*) 
00380  write(6,'(a40,x,i10)')'number of gw_grd =',nsgm 
00381  write(6,'(a40,x,f15.7)')'emax (eV)=',emax*au
00382  write(6,'(a40,x,f15.7)')'emin (eV)=',emin*au
00383  write(6,'(a40,x,f15.7)')'ecmax (eV)=',ecmax*au
00384  write(6,'(a40,x,f15.7)')'ecmin (eV)=',ecmin*au
00385  write(6,'(a40,x,f15.7)')'gw_grid_separation (eV)=',gw_grid_separation*au
00386  write(6,*) 
00387  !write(6,*)'============================'
00388  !write(6,*)'FREQUENCY GRID IN EPSQW (eV)'
00389  !write(6,*)'============================'
00390  !write(6,*) 
00391  !do ie=1,ne 
00392  ! write(6,'(2f15.7)') em(ie)*au 
00393  !enddo 
00394  !write(6,*) 
00395  !write(6,*)'==========================='
00396  !write(6,*)'dir-eps: POLE OF MODEL GRID'
00397  !write(6,*)'==========================='
00398  !write(6,*) 
00399  !do ie=1,ne 
00400  ! write(6,'(2f15.10)') pole_of_chi(ie)
00401  !enddo 
00402  !write(6,*) 
00403  !write(6,*)'========================='
00404  !write(6,*)'FREQUENCY GRID IN GW (eV)'
00405  !write(6,*)'========================='
00406  !write(6,*) 
00407  !do ie=1,nsgm  
00408  ! write(6,'(f15.7)') sgmw(ie)*au 
00409  !enddo 
00410  !write(6,*) 
00411 endif!myrank.eq.0
00412 !
00413 call MPI_Bcast(nsgm,1,MPI_INTEGER,0,comm,ierr) 
00414 call MPI_BARRIER(comm,ierr)
00415 if(myrank/=0) allocate(sgmw(nsgm))
00416 call MPI_BARRIER(comm,ierr)
00417 call MPI_Bcast(sgmw,nsgm,MPI_DOUBLE_PRECISION,0,comm,ierr) 
00418 call MPI_BARRIER(comm,ierr)
00419 !
00420 !mkdir dir-gw
00421 !
00422 if(myrank.eq.0)then!master  
00423  call system('mkdir -p dir-gw')!20180509
00424  ierr=CHDIR("./dir-gw") 
00425  call system('pwd') 
00426  ! 
00427  !OPEN(136,R,FILE='./dat.gwgrid') 
00428  ! 
00429  OPEN(136,FILE='./dat.gwgrid') 
00430  rewind(136) 
00431  write(136,*) nsgm 
00432  do ie=1,nsgm 
00433   write(136,'(f20.10)') sgmw(ie) 
00434  enddo 
00435  ierr=CHDIR("..") 
00436  call system('pwd') 
00437 endif!myrankeq0 
00438 flush(6) 
00439 !
00440 !MPI 
00441 !
00442 pnq=NTQ/nproc 
00443 if(mod(NTQ,nproc).ne.0) then
00444  nbufq=pnq+1
00445 else
00446  nbufq=pnq
00447 end if
00448 if(myrank.lt.mod(NTQ,nproc)) then
00449  pnq=pnq+1
00450  bnq=pnq*myrank+1 
00451  enq=bnq+pnq-1
00452 else
00453  bnq=(pnq+1)*mod(NTQ,nproc)+pnq*(myrank-mod(NTQ,nproc))+1 
00454  enq=bnq+pnq-1
00455 end if
00456 !
00457 !######
00458 !  SC 
00459 !######
00460 !
00461 if(calc_sc)then!.true.=default 
00462  !
00463  !write(file_id,'(a)')'## SC calc start ##'
00464  !write(file_id,'(a,3i7)')'bnq,enq,pnq',bnq,enq,pnq 
00465  !
00466  if(myrank.eq.0)then 
00467   write(6,*)'## SC calc start ##' 
00468   write(6,*) 
00469  endif 
00470  !
00471  !fft
00472  !
00473  nfft1=nwx2+1; nfft2=nwy2+1; nfft3=nwz2+1; Nl123=nfft1*nfft2*nfft3 
00474  call fft3_init(nwx2,nwy2,nwz2,nfft1,nfft2,nfft3,fs) 
00475  !
00476  allocate(pSC(nsgm,Mb,Mb,Nk_irr)); pSC(:,:,:,:)=0.0d0 
00477  do iq=1,pnq
00478   !
00479   q1=SQ(1,bnq+iq-1);q2=SQ(2,bnq+iq-1);q3=SQ(3,bnq+iq-1)
00480   !
00481   !q.ne.0 
00482   !
00483   if(q1/=0.0d0.or.q2/=0.0d0.or.q3/=0.0d0)then 
00484    NG_for_eps=NGQ_eps(bnq+iq-1)
00485    allocate(length_qg(NG_for_eps)); length_qg=0.0D0 
00486    allocate(atten_factor(NG_for_eps)); atten_factor=0.0D0 
00487    allocate(epsmk(NTGQ,NTGQ,ne)); epsmk=0.0d0 
00488    iqir=numirrq(bnq+iq-1)
00489    iop=numrotq(bnq+iq-1)
00490    !
00491    call make_eps(NTG,NTGQ,ne,trsq(bnq+iq-1),NGQ_eps(bnq+iq-1),LG0(1,1,bnq+iq-1),RWq(1,bnq+iq-1),&
00492    rginv(1,1,iop),pg(1,iop),nnp,Lq1,Lq2,Lq3,packingq(-Lq1,-Lq2,-Lq3,iqir),epsirr(1,1,1,iqir),&
00493    epsmk(1,1,1)) 
00494    !
00495    length_qg(:)=0.0d0 
00496    atten_factor(:)=0.0d0 
00497    do igL=1,NG_for_eps   
00498     igL1=LG0(1,igL,bnq+iq-1)
00499     igL2=LG0(2,igL,bnq+iq-1)
00500     igL3=LG0(3,igL,bnq+iq-1)
00501     qgL(:)=(q1+dble(igL1))*b1(:)+(q2+dble(igL2))*b2(:)+(q3+dble(igL3))*b3(:)
00502     qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00503     qgL1=dsqrt(qgL2) 
00504     length_qg(igL)=qgL1
00505     atten_factor(igL)=dsqrt(1.0d0-dcos(qgL1*Rc))   
00506    enddo!igL 
00507    !
00508 !$OMP PARALLEL PRIVATE(ib,ik,jb,shift_G,ikq,C0_K,C0_KmQ,rho,wfunc,fftwk,rho_tmp,ikir,iop,kb,ie,SUM_CMPX,&
00509 !$OMP&         ikqir,ikqop,igL,jgL,vecf,je,veca,delta,de,sgn,dnm,pSComp) 
00510    allocate(rho(NG_for_eps,Mb,Nk_irr));rho=0.0d0
00511    allocate(fftwk(Nl123*2));fftwk=0.0d0 
00512    allocate(wfunc(Nl123*2));wfunc=0.0d0 
00513    allocate(rho_tmp(NG_for_eps));rho_tmp=0.0d0
00514    allocate(pSComp(nsgm,Mb,Mb,Nk_irr)); pSComp=0.0d0  
00515    allocate(C0_K(NTG));C0_K=0.0d0
00516    allocate(C0_KmQ(NTG));C0_KmQ=0.0d0
00517    allocate(vecf(ne));vecf=0.0d0
00518    allocate(veca(ne));veca=0.0d0
00519 !$OMP DO  
00520    do ib=1,Ncalc
00521     if(myrank.eq.0) write(6,*)'ib=',ib
00522     !
00523     rho(:,:,:)=0.0D0 
00524     !
00525     do ikir=1,Nk_irr
00526      !
00527      ik=numMK(ikir) 
00528      iop=numrot(ik) 
00529      !
00530      do jb=1,Nb(ik)!Mb  
00531       !
00532       call make_C0_for_given_band(NTG,trs(ik),NG0(ik),KG0(1,1,ik),RW(1,ik),rginv(1,1,iop),pg(1,iop),nnp,&
00533       L1,L2,L3,packing(-L1,-L2,-L3,ikir),CIR(1,jb+Ns(ik),ikir),C0_K(1)) 
00534       ! 
00535       shift_G(:)=0
00536       call search_kq(NTK,SK0(1,1),-q1,-q2,-q3,ik,ikq,shift_G(1))
00537       shift_G(:)=-shift_G(:)
00538       !
00539       ikqir=numirr(ikq)
00540       ikqop=numrot(ikq) 
00541       !
00542       call make_C0_for_given_band(NTG,trs(ikq),NG0(ikq),KG0(1,1,ikq),RW(1,ikq),rginv(1,1,ikqop),pg(1,ikqop),nnp,&
00543       L1,L2,L3,packing(-L1,-L2,-L3,ikqir),CIR(1,ib,ikqir),C0_KmQ(1)) 
00544       !
00545       call calc_InterStateMatrix(NTK,NTG,NG0(1),KG0(1,1,1),C0_KmQ(1),C0_K(1),ikq,ik,nwx2,nwy2,nwz2,&
00546       nfft1,nfft2,Nl123,wfunc(1),fftwk(1),fs,LG0(1,1,bnq+iq-1),NG_for_eps,shift_G(1),rho_tmp(1))
00547       !
00548       rho(:,jb,ikir)=rho_tmp(:)
00549       !
00550      enddo!jb  
00551     enddo!ikir 
00552     !
00553     do ikir=1,Nk_irr 
00554      !
00555      ik=numMK(ikir) 
00556      !
00557      shift_G(:)=0
00558      call search_kq(NTK,SK0(1,1),-q1,-q2,-q3,ik,ikq,shift_G(1))
00559      !
00560      do jb=1,Nb(ik)!Mb 
00561       do kb=1,Nb(ik)!Mb 
00562        vecf=0.0d0
00563        do ie=1,ne 
00564         SUM_CMPX=0.0d0
00565         do igL=1,NG_for_eps  
00566          do jgL=1,NG_for_eps 
00567           !
00568           SUM_CMPX=SUM_CMPX+rho(igL,jb,ikir)/length_qg(igL)*atten_factor(igL)*epsmk(igL,jgL,ie)&
00569           *CONJG(rho(jgL,kb,ikir))/length_qg(jgL)*atten_factor(jgL) 
00570           !
00571          enddo!jgL 
00572         enddo!igL 
00573         vecf(ie)=2.0d0*tpi*SUM_CMPX/dble(NTQ)/VOLUME 
00574        enddo!ie 
00575        !
00576        veca=0.0d0 
00577        do je=1,ne 
00578         SUM_CMPX=0.0d0 
00579         do ie=1,ne 
00580          SUM_CMPX=SUM_CMPX+mat_b(je,ie)*vecf(ie) 
00581         enddo!ie 
00582         veca(je)=SUM_CMPX
00583        enddo!je 
00584        !
00585        do ie=1,nsgm 
00586         delta=sgmw(ie) 
00587         !
00588         ikqir=numirr(ikq) 
00589         de=delta-E_EIGI(ib,ikqir)
00590         !
00591         if(E_EIGI(ib,ikqir)>FermiEnergy)then 
00592          sgn=1.0d0
00593         else
00594          sgn=-1.0d0
00595         endif 
00596         !
00597         do je=1,ne 
00598          !
00599          dnm=de-(pole_of_chi(je)-ci*idlt)*sgn
00600          !
00601          pSComp(ie,jb,kb,ikir)=pSComp(ie,jb,kb,ikir)+veca(je)/dnm 
00602          !
00603         enddo!je 
00604        enddo!ie  
00605       enddo!kb 
00606      enddo!jb 
00607     enddo!ikir  
00608    enddo!ib 
00609 !$OMP END DO
00610 !$OMP CRITICAL
00611    pSC=pSC+pSComp
00612 !$OMP END CRITICAL
00613    deallocate(fftwk,wfunc,rho_tmp,rho,pSComp,C0_K,C0_KmQ,vecf,veca)
00614 !$OMP END PARALLEL
00615    if(myrank.eq.0) write(6,*)'FINISH iq',iq 
00616    deallocate(length_qg,atten_factor,epsmk) 
00617    !
00618   elseif(q1==0.0d0.and.q2==0.0d0.and.q3==0.0d0) then 
00619    !
00620    !q.eq.0 
00621    !
00622    NG_for_eps=NGQ_eps(bnq+iq-1)
00623    allocate(length_qg(NG_for_eps));length_qg=0.0d0 
00624    allocate(atten_factor(NG_for_eps));atten_factor=0.0d0  
00625    allocate(rho(NG_for_eps,Mb,Nk_irr));rho=0.0d0 
00626    allocate(epsmk(NTGQ,NTGQ,ne));epsmk(:,:,:)=0.0d0 
00627    iqir=numirrq(bnq+iq-1)
00628    iop=numrotq(bnq+iq-1)
00629    !
00630    call make_eps(NTG,NTGQ,ne,trsq(bnq+iq-1),NGQ_eps(bnq+iq-1),LG0(1,1,bnq+iq-1),RWq(1,bnq+iq-1),&
00631    rginv(1,1,iop),pg(1,iop),nnp,Lq1,Lq2,Lq3,packingq(-Lq1,-Lq2,-Lq3,iqir),epsirr(1,1,1,iqir),&
00632    epsmk(1,1,1)) 
00633    !
00634    length_qg(:)=0.0D0 
00635    do igL=1,NG_for_eps   
00636     igL1=LG0(1,igL,bnq+iq-1)
00637     igL2=LG0(2,igL,bnq+iq-1)
00638     igL3=LG0(3,igL,bnq+iq-1)
00639     if(igL1==0.and.igL2==0.and.igL3==0)then 
00640      No_G_0=igL  
00641      !20180822
00642      length_qg(igL)=1.0d0 
00643      atten_factor(igL)=0.0d0 
00644     else
00645      qgL(:)=(q1+dble(igL1))*b1(:)+(q2+dble(igL2))*b2(:)+(q3+dble(igL3))*b3(:)
00646      qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00647      qgL1=dsqrt(qgL2) 
00648      length_qg(igL)=qgL1
00649      atten_factor(igL)=dsqrt(1.0d0-dcos(qgL1*Rc))   
00650     endif 
00651    enddo!igL 
00652    !
00653    do ib=1,Ncalc
00654     if(myrank.eq.0) write(6,*)'ib=',ib
00655     rho(:,:,:)=0.0D0 
00656 !$OMP PARALLEL PRIVATE(ik,jb,ikqir,ikqop,shift_G,ikq,C0_K,C0_KmQ,wfunc,fftwk,rho_tmp,ikir,iop) 
00657     allocate(fftwk(Nl123*2));fftwk=0.0d0 
00658     allocate(wfunc(Nl123*2));wfunc=0.0d0  
00659     allocate(rho_tmp(NG_for_eps));rho_tmp=0.0d0 
00660     allocate(C0_K(NTG));C0_K=0.0d0
00661     allocate(C0_KmQ(NTG));C0_KmQ=0.0d0
00662 !$OMP DO
00663     do ikir=1,Nk_irr 
00664      !
00665      ik=numMK(ikir) 
00666      iop=numrot(ik) 
00667      !
00668      do jb=1,Nb(ik)!Mb 
00669       !
00670       call make_C0_for_given_band(NTG,trs(ik),NG0(ik),KG0(1,1,ik),RW(1,ik),rginv(1,1,iop),pg(1,iop),nnp,&
00671       L1,L2,L3,packing(-L1,-L2,-L3,ikir),CIR(1,jb+Ns(ik),ikir),C0_K(1)) 
00672       !
00673       shift_G(:)=0
00674       ikq=ik 
00675       !
00676       ikqir=numirr(ikq)
00677       ikqop=numrot(ikq) 
00678       !
00679       call make_C0_for_given_band(NTG,trs(ikq),NG0(ikq),KG0(1,1,ikq),RW(1,ikq),rginv(1,1,ikqop),pg(1,ikqop),nnp,&
00680       L1,L2,L3,packing(-L1,-L2,-L3,ikqir),CIR(1,ib,ikqir),C0_KmQ(1)) 
00681       !
00682       call calc_InterStateMatrix(NTK,NTG,NG0(1),KG0(1,1,1),C0_KmQ(1),C0_K(1),ikq,ik,nwx2,nwy2,nwz2,&
00683       nfft1,nfft2,Nl123,wfunc(1),fftwk(1),fs,LG0(1,1,bnq+iq-1),NG_for_eps,shift_G(1),rho_tmp(1))
00684       !
00685       rho(:,jb,ikir)=rho_tmp(:)
00686       !
00687      enddo!jb
00688     enddo!ikir 
00689 !$OMP END DO
00690     deallocate(fftwk,wfunc,rho_tmp,C0_K,C0_KmQ)
00691 !$OMP END PARALLEL
00692 !--
00693 !$OMP PARALLEL PRIVATE(ik,jb,kb,ie,SUM_CMPX,igL,jgL,vecf,je,veca,ikq,delta,ikir,ikqir,de,sgn,dnm) 
00694     allocate(vecf(ne));vecf=0.0d0
00695     allocate(veca(ne));veca=0.0d0
00696 !$OMP DO
00697     do ikir=1,Nk_irr 
00698      !
00699      ik=numMK(ikir)
00700      !
00701      do jb=1,Nb(ik)!Mb 
00702       do kb=1,Nb(ik)!Mb 
00703        vecf=0.0d0 
00704        do ie=1,ne 
00705         SUM_CMPX=0.0d0
00706         do igL=1,NG_for_eps 
00707          !if(igL==No_G_0) cycle 
00708          do jgL=1,NG_for_eps 
00709           !if(jgL==No_G_0) cycle 
00710           !
00711           SUM_CMPX=SUM_CMPX+rho(igL,jb,ikir)/length_qg(igL)*atten_factor(igL)*epsmk(igL,jgL,ie)&
00712           *CONJG(rho(jgL,kb,ikir))/length_qg(jgL)*atten_factor(jgL) 
00713           !
00714          enddo!jgL 
00715         enddo!igL  
00716         vecf(ie)=2.0d0*tpi*SUM_CMPX/VOLUME/dble(NTQ)  
00717        enddo!ie 
00718        veca=0.0d0
00719        do je=1,ne 
00720         SUM_CMPX=0.0d0 
00721         do ie=1,ne 
00722          SUM_CMPX=SUM_CMPX+mat_b(je,ie)*vecf(ie) 
00723         enddo!ie  
00724         veca(je)=SUM_CMPX 
00725        enddo!je  
00726        do ie=1,nsgm 
00727         ikq=ik 
00728         delta=sgmw(ie) 
00729         !
00730         ikqir=numirr(ikq) 
00731         de=delta-E_EIGI(ib,ikqir)
00732         !
00733         if(E_EIGI(ib,ikqir)>FermiEnergy)then 
00734          sgn=1.0d0
00735         else
00736          sgn=-1.0d0
00737         endif 
00738         !
00739         do je=1,ne 
00740          dnm=de-(pole_of_chi(je)-ci*idlt)*sgn
00741          !
00742          pSC(ie,jb,kb,ikir)=pSC(ie,jb,kb,ikir)+veca(je)/dnm 
00743          !
00744         enddo!je  
00745        enddo!ie  
00746       enddo!kb 
00747      enddo!jb 
00748     enddo!ikir 
00749 !$OMP END DO
00750     deallocate(vecf,veca)
00751 !$OMP END PARALLEL
00752    enddo!ib 
00753    deallocate(length_qg,atten_factor,rho,epsmk) 
00754    !
00755    !<head contribution> 
00756    !
00757    !(i) Spencer-Alabi
00758    !
00759    chead=(tpi/dble(NTQ)/VOLUME)*Rc*Rc   
00760    !
00761    !(ii) Hybertsen-Louie
00762    !
00763    !qsz=(6.0D0*(pi**2)/dble(NTQ)/dble(VOLUME))**(1.0D0/3.0D0)  
00764    !chead=(2.0D0/pi)*qsz 
00765    !
00766    !write(file_id,*)'bnq+iq-1=',bnq+iq-1 
00767    !write(file_id,*)'correction_head=',chead   
00768    !
00769    !cGW
00770    !
00771    !chead=0.0d0    
00772    !write(file_id,*)'cGW, then correction_head=0'
00773    !
00774    !write(file_id,*)'No_G_0=',No_G_0 
00775    allocate(vecf(ne));vecf=0.0d0
00776    allocate(veca(ne));veca=0.0d0
00777    do ikir=1,Nk_irr 
00778     !
00779     ik=numMK(ikir) 
00780     !
00781     do jb=1,Nb(ik)!Mb 
00782      !
00783      vecf=0.0d0 
00784      do ie=1,ne 
00785       vecf(ie)=epsirr(No_G_0,No_G_0,ie,bnq+iq-1)
00786      enddo!ie   
00787      !
00788      veca=0.0d0 
00789      do je=1,ne 
00790       SUM_CMPX=0.0d0 
00791       do ie=1,ne 
00792        SUM_CMPX=SUM_CMPX+mat_b(je,ie)*vecf(ie) 
00793       enddo!ie 
00794       veca(je)=SUM_CMPX 
00795      enddo!je  
00796      !
00797      do ie=1,nsgm  
00798       delta=sgmw(ie) 
00799       !
00800       de=delta-E_EIGI(jb+Ns(ik),ikir) 
00801       !
00802       if(E_EIGI(jb+Ns(ik),ikir)>FermiEnergy)then 
00803        sgn=1.0d0
00804       else
00805        sgn=-1.0d0
00806       endif 
00807       do je=1,ne 
00808        dnm=de-(pole_of_chi(je)-ci*idlt)*sgn
00809        !
00810        pSC(ie,jb,jb,ikir)=pSC(ie,jb,jb,ikir)+veca(je)/dnm*chead  
00811        !
00812       enddo!je
00813      enddo!ie  
00814     enddo!jb 
00815    enddo!ikir 
00816    deallocate(vecf,veca) 
00817    if(myrank.eq.0) write(6,*)'FINISH iq',iq 
00818   endif 
00819  enddo!iq 
00820  call MPI_BARRIER(comm,ierr)
00821  !write(file_id,*)'I finished pSC calc' 
00822  allocate(SCirr(nsgm,Mb,Mb,Nk_irr)); SCirr(:,:,:,:)=0.0d0 
00823  !
00824  !20180922
00825  !
00826  !mem_size=dble(nsgm*Mb*Mb*Nk_irr)*16.0d0/1024.0d0/1024.0d0/1024.0d0 
00827   mem_size=dble(nsgm*Mb*Mb*Nk_irr)*8.0d0/1024.0d0/1024.0d0/1024.0d0 
00828  if(myrank.eq.0)then 
00829   write(6,'(a35,f20.15)')'mem size pSC or SCirr (GB)',mem_size
00830  endif 
00831  !call MPI_REDUCE(pSC,SCirr,nsgm*Mb*Mb*Nk_irr,MPI_DOUBLE_COMPLEX,MPI_SUM,0,comm,ierr)
00832  call MPI_REDUCE(pSC,SCirr,nsgm*Mb*Mb*Nk_irr,MPI_COMPLEX,MPI_SUM,0,comm,ierr)
00833  call MPI_BARRIER(comm,ierr)
00834  !call MPI_Bcast(SCirr,nsgm*Mb*Mb*Nk_irr,MPI_DOUBLE_COMPLEX,0,comm,ierr) 
00835  call MPI_Bcast(SCirr,nsgm*Mb*Mb*Nk_irr,MPI_COMPLEX,0,comm,ierr) 
00836  call MPI_BARRIER(comm,ierr)
00837  !
00838  !MPI for nsgm 
00839  !
00840  pnw=nsgm/nproc 
00841  if(mod(nsgm,nproc).ne.0) then
00842   nbufw=pnw+1
00843  else
00844   nbufw=pnw 
00845  end if
00846  if(myrank.lt.mod(nsgm,nproc)) then
00847   pnw=pnw+1
00848   bnw=pnw*myrank+1 
00849   enw=bnw+pnw-1
00850  else
00851   bnw=(pnw+1)*mod(nsgm,nproc)+pnw*(myrank-mod(nsgm,nproc))+1 
00852   enw=bnw+pnw-1
00853  end if
00854  !write(file_id,'(a,3i7)')'bnw,enw,pnw',bnw,enw,pnw            
00855  !
00856  allocate(SC(pnw,Mb,Mb,NTK)); SC=0.0d0
00857  ! 
00858  !mem_size=dble(pnw*Mb*Mb*NTK)*16.0d0/1024.0d0/1024.0d0/1024.0d0 
00859  !if(myrank.eq.0)then 
00860  ! write(6,'(a35,f20.15)')'mem size SC (GB)',mem_size
00861  !endif 
00862  !
00863  do ik=1,NTK 
00864   ikir=numirr(ik) 
00865   if(trs(ik)==1)then 
00866    do ie=1,pnw!nsgm
00867     do ib=1,Nb(ik)!Mb
00868      do jb=1,Nb(ik)!Mb 
00869       SC(ie,ib,jb,ik)=SCirr(bnw+ie-1,ib,jb,ikir) 
00870      enddo 
00871     enddo 
00872    enddo 
00873   elseif(trs(ik)==-1)then 
00874    do ie=1,pnw!nsgm  
00875     do ib=1,Nb(ik)!Mb
00876      do jb=1,Nb(ik)!Mb 
00877       SC(ie,ib,jb,ik)=SCirr(bnw+ie-1,jb,ib,ikir) 
00878      enddo 
00879     enddo 
00880    enddo 
00881   endif 
00882  enddo 
00883  deallocate(pSC)
00884  !
00885  allocate(pf(NTK,-Na1:Na1,-Na2:Na2,-Na3:Na3)); pf=0.0d0     
00886  do ik=1,NTK 
00887   do ia3=-Na3,Na3 
00888    do ia2=-Na2,Na2 
00889     do ia1=-Na1,Na1 
00890      phase=tpi*(SK0(1,ik)*dble(ia1)+SK0(2,ik)*dble(ia2)+SK0(3,ik)*dble(ia3)) 
00891      pf(ik,ia1,ia2,ia3)=exp(-ci*phase) 
00892     enddo  
00893    enddo  
00894   enddo  
00895  enddo  
00896  !
00897  !make pSCR
00898  !
00899  if(myrank.eq.0)then
00900   write(6,*) 
00901   write(6,*)'start SCR' 
00902   write(6,*) 
00903  endif 
00904  !
00905  allocate(pSCR(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3,pnw)); pSCR=0.0d0 
00906  !--
00907  !OMP CRITICAL: off-diagonal case (default)
00908  !--
00909  do ie=1,pnw!nsgm  
00910   do ia1=-Na1,Na1
00911    do ia2=-Na2,Na2
00912     do ia3=-Na3,Na3
00913      do iw=1,NWF
00914       do jw=1,NWF 
00915        SUM_CMPX=0.0D0 
00916 !$OMP PARALLEL private(psum,ik,jb,kb)  
00917         psum=0.0d0 
00918 !$OMP DO 
00919         do ik=1,NTK 
00920          do jb=1,Nb(ik) 
00921           do kb=1,Nb(ik) 
00922            psum=psum+CONJG(UNT(jb,iw,ik))*SC(ie,jb,kb,ik)*UNT(kb,jw,ik)*pf(ik,ia1,ia2,ia3) 
00923           enddo!kb 
00924          enddo!jb 
00925         enddo!ik 
00926 !$OMP END DO 
00927 !$OMP CRITICAL 
00928         SUM_CMPX=SUM_CMPX+psum 
00929 !$OMP END CRITICAL 
00930 !$OMP END PARALLEL 
00931        pSCR(iw,jw,ia1,ia2,ia3,ie)=SUM_CMPX/DBLE(NTK)
00932       enddo!jw
00933      enddo!iw
00934     enddo!ia3
00935    enddo!ia2 
00936   enddo!ia1 
00937   !write(file_id,'(a,i7)')'ie=',ie 
00938   !if(myrank.eq.0) write(6,*)'FINISH ie',ie  
00939   write(6,*)'FINISH ie myrank',ie,myrank   
00940  enddo!ie
00941  !--
00942  !OMP CRITICAL: diagonal case (not default)
00943  !--
00944  if(.false.)then 
00945   do ie=1,pnw!nsgm  
00946    do ia1=-Na1,Na1
00947     do ia2=-Na2,Na2
00948      do ia3=-Na3,Na3
00949       do iw=1,NWF
00950        do jw=1,NWF 
00951         SUM_CMPX=0.0D0 
00952 !$OMP PARALLEL private(psum,ik,jb)  
00953          psum=0.0d0 
00954 !$OMP DO 
00955          do ik=1,NTK 
00956           do jb=1,Nb(ik) 
00957            psum=psum+CONJG(UNT(jb,iw,ik))*SC(ie,jb,jb,ik)*UNT(jb,jw,ik)*pf(ik,ia1,ia2,ia3) 
00958           enddo!jb 
00959          enddo!ik 
00960 !$OMP END DO 
00961 !$OMP CRITICAL 
00962         SUM_CMPX=SUM_CMPX+psum 
00963 !$OMP END CRITICAL 
00964 !$OMP END PARALLEL 
00965         pSCR(iw,jw,ia1,ia2,ia3,ie)=SUM_CMPX/DBLE(NTK)
00966        enddo!jw
00967       enddo!iw
00968      enddo!ia3
00969     enddo!ia2 
00970    enddo!ia1 
00971    !write(file_id,'(a,i7)')'ie=',ie 
00972    if(myrank.eq.0) write(6,*)'FINISH ie',ie  
00973   enddo!ie
00974  endif
00975  !
00976  deallocate(SC,pf) 
00977  call MPI_BARRIER(comm,ierr)
00978  !write(file_id,*)'I finished pSCR calc' 
00979  !
00980  !make MAT_SC_R 
00981  !
00982  if(myrank.eq.0)then 
00983   allocate(MAT_SC_R(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3,nsgm)); MAT_SC_R=0.0d0 
00984  endif 
00985  !
00986  ! 
00987  mem_size=dble(pnw*NWF*NWF*(2*Na1+1)*(2*Na2+1)*(2*Na3+1))*8.0d0/1024.0d0/1024.0d0/1024.0d0 
00988  if(myrank.eq.0)then 
00989   write(6,'(a35,f20.15)')'mem size pSCR (GB)',mem_size
00990  endif 
00991  ! 
00992  !
00993  call MPI_Gather(pSCR,pnw*NWF*NWF*(2*Na1+1)*(2*Na2+1)*(2*Na3+1),MPI_COMPLEX,&
00994                  MAT_SC_R,pnw*NWF*NWF*(2*Na1+1)*(2*Na2+1)*(2*Na3+1),MPI_COMPLEX,0,comm,ierr)
00995  !
00996  !
00997  if(myrank.eq.0)then 
00998   !write(5004) MAT_SC_R 
00999   write(6,*)'finish MAT_SC_R' 
01000  endif!myrank.eq.0 
01001  !
01002 else
01003  write(6,*)'skip SC calc.' 
01004  if(myrank.eq.0)then 
01005   allocate(MAT_SC_R(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3,nsgm)); MAT_SC_R=0.0d0 
01006  endif 
01007 endif!calc SC or not 
01008 !
01009 !#############
01010 !  OUTPUT SC 
01011 !#############
01012 !
01013 if(calc_sc)then!.true.=default 
01014  if(myrank.eq.0)then 
01015   ierr=CHDIR("./dir-gw") 
01016   call system('pwd') 
01017   !
01018   !OPEN(156,FILE='./dat.energy_vs_sc') 
01019   !
01020   OPEN(156,FILE='./dat.energy_vs_sc') 
01021   rewind(156) 
01022   do ikir=1,Nk_irr 
01023    ik=numMK(ikir)
01024    do jb=1,Nb(ik)!Mb 
01025     min_diff=1.0d0 
01026     do ie=1,nsgm
01027      en=abs(sgmw(ie)-E_EIGI(jb+Ns(ik),ikir))  
01028      if(en<min_diff)then 
01029       min_diff=en 
01030       min_ie=ie
01031      endif  
01032     enddo!ie 
01033     en=E_EIGI(jb+Ns(ik),ikir)-FermiEnergy 
01034     write(156,'(3F20.10)') en*au,SCirr(min_ie,jb,jb,ikir)*au 
01035    enddo!jb 
01036   enddo!ikir 
01037   close(156) 
01038   !
01039   !OPEN(157,W,FILE='MAT_SC_R') 
01040   !
01041   !rewind(157)
01042   !do ia1=-Na1,Na1
01043   ! do ia2=-Na2,Na2 
01044   !  do ia3=-Na3,Na3 
01045   !   do jw=1,NWF
01046   !    do iw=1,NWF 
01047   !     write(157)(MAT_SC_R(ie,iw,jw,ia1,ia2,ia3),ie=1,nsgm) 
01048   !    enddo!iw  
01049   !   enddo!jw  
01050   !  enddo!ia3 
01051   ! enddo!ia2 
01052   !enddo!ia1 
01053   !close(157) 
01054   !
01055   ierr=CHDIR("..") 
01056   call system('pwd') 
01057  endif!myrank.eq.0
01058 endif!calc SC or not 
01059 !
01060 deallocate(SCirr) 
01061 !
01062 !######
01063 !  SX
01064 !######
01065 !
01066 !write(file_id,'(a)')'## SX calc start ##'
01067 !
01068 if(myrank.eq.0)then 
01069  write(6,*)'## SX calc start ##' 
01070  write(6,*) 
01071 endif 
01072 !
01073 !fft
01074 !
01075 nfft1=nwx2+1; nfft2=nwy2+1; nfft3=nwz2+1; Nl123=nfft1*nfft2*nfft3 
01076 call fft3_init(nwx2,nwy2,nwz2,nfft1,nfft2,nfft3,fs) 
01077 !
01078 allocate(SXirr(Mb,Mb,Nk_irr));SXirr(:,:,:)=0.0d0 
01079 allocate(pSX(Mb,Mb,Nk_irr));pSX(:,:,:)=0.0d0 
01080 do iq=1,pnq 
01081  !
01082  q1=SQ(1,bnq+iq-1); q2=SQ(2,bnq+iq-1); q3=SQ(3,bnq+iq-1)
01083  !
01084  !q.ne.0 
01085  !
01086  if(q1/=0.0d0.or.q2/=0.0d0.or.q3/=0.0d0)then 
01087   NG_for_psi=NGQ_psi(bnq+iq-1)
01088   allocate(atten_factor(NG_for_psi));atten_factor(:)=0.0D0 
01089   allocate(length_qg(NG_for_psi));length_qg(:)=0.0D0 
01090   allocate(rho(NG_for_psi,Mb,Nk_irr))
01091   do igL=1,NG_for_psi   
01092    igL1=LG0(1,igL,bnq+iq-1)
01093    igL2=LG0(2,igL,bnq+iq-1)
01094    igL3=LG0(3,igL,bnq+iq-1)
01095    qgL(:)=(q1+dble(igL1))*b1(:)+(q2+dble(igL2))*b2(:)+(q3+dble(igL3))*b3(:)
01096    qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
01097    qgL1=dsqrt(qgL2) 
01098    length_qg(igL)=qgL1
01099    atten_factor(igL)=1.0d0-dcos(qgL1*Rc)  
01100   enddo 
01101   do ib=1,Ncalc!Mt
01102    rho(:,:,:)=0.0D0 
01103 !$OMP PARALLEL PRIVATE(ik,shift_G,ikq,C0_K,C0_KmQ,jb,wfunc,fftwk,rho_tmp,ikir,iop,ikqir,ikqop) 
01104    allocate(fftwk(Nl123*2));fftwk=0.0d0 
01105    allocate(wfunc(Nl123*2));wfunc=0.0d0 
01106    allocate(rho_tmp(NG_for_psi));rho_tmp=0.0d0 
01107    allocate(C0_K(NTG));C0_K=0.0d0
01108    allocate(C0_KmQ(NTG));C0_KmQ=0.0d0
01109 !$OMP DO
01110    do ikir=1,Nk_irr 
01111     !
01112     ik=numMK(ikir) 
01113     iop=numrot(ik) 
01114     !
01115     do jb=1,Nb(ik)!Mb 
01116      !
01117      call make_C0_for_given_band(NTG,trs(ik),NG0(ik),KG0(1,1,ik),RW(1,ik),rginv(1,1,iop),pg(1,iop),nnp,&
01118      L1,L2,L3,packing(-L1,-L2,-L3,ikir),CIR(1,jb+Ns(ik),ikir),C0_K(1)) 
01119      ! 
01120      shift_G(:)=0
01121      call search_kq(NTK,SK0(1,1),-q1,-q2,-q3,ik,ikq,shift_G(1))
01122      shift_G(:)=-shift_G(:)  
01123      !
01124      ikqir=numirr(ikq)
01125      ikqop=numrot(ikq) 
01126      !
01127      call make_C0_for_given_band(NTG,trs(ikq),NG0(ikq),KG0(1,1,ikq),RW(1,ikq),rginv(1,1,ikqop),pg(1,ikqop),nnp,&
01128      L1,L2,L3,packing(-L1,-L2,-L3,ikqir),CIR(1,ib,ikqir),C0_KmQ(1)) 
01129      !
01130      !
01131      call calc_InterStateMatrix(NTK,NTG,NG0(1),KG0(1,1,1),C0_KmQ(1),C0_K(1),ikq,ik,nwx2,nwy2,nwz2,&
01132      nfft1,nfft2,Nl123,wfunc(1),fftwk(1),fs,LG0(1,1,bnq+iq-1),NG_for_psi,shift_G(1),rho_tmp(1))
01133      !
01134      rho(:,jb,ikir)=rho_tmp(:)
01135      !
01136     enddo!jb  
01137    enddo!ik 
01138 !$OMP END DO 
01139    deallocate(fftwk,wfunc,rho_tmp,C0_K,C0_KmQ)
01140 !$OMP END PARALLEL 
01141 !$OMP PARALLEL PRIVATE(ik,ikir,shift_G,ikq,jb,kb,SUM_CMPX,igL)  
01142 !$OMP DO
01143    do ikir=1,Nk_irr 
01144     !
01145     ik=numMK(ikir) 
01146     ! 
01147     shift_G(:)=0
01148     call search_kq(NTK,SK0(1,1),-q1,-q2,-q3,ik,ikq,shift_G(1))
01149     !
01150     do jb=1,Nb(ik)!Mb 
01151      do kb=1,Nb(ik)!Mb 
01152       SUM_CMPX=0.0d0
01153       do igL=1,NG_for_psi 
01154        !
01155        SUM_CMPX=SUM_CMPX+fbk(ib,ikq)*rho(igL,jb,ikir)&
01156        *CONJG(rho(igL,kb,ikir))/((length_qg(igL))**2)*atten_factor(igL)
01157        !
01158       enddo 
01159       !
01160       pSX(jb,kb,ikir)=pSX(jb,kb,ikir)+SUM_CMPX 
01161       !
01162      enddo!kb 
01163     enddo!jb 
01164    enddo!ik 
01165 !$OMP END DO 
01166 !$OMP END PARALLEL 
01167   enddo!ib 
01168   !write(file_id,*)'FINISH iq',iq
01169   if(myrank.eq.0) write(6,*)'FINISH iq',iq 
01170   deallocate(length_qg,atten_factor,rho) 
01171  elseif(q1==0.0d0.and.q2==0.0d0.and.q3==0.0d0)then 
01172   !
01173   !q.eq.0 
01174   !
01175   !write(file_id,'(a,x,3f10.5)')'q1,q2,q3=',q1,q2,q3  
01176   NG_for_psi=NGQ_psi(bnq+iq-1)
01177   allocate(length_qg(NG_for_psi));length_qg(:)=0.0D0 
01178   allocate(atten_factor(NG_for_psi));atten_factor(:)=0.0D0 
01179   allocate(rho(NG_for_psi,Mb,Nk_irr))
01180   do igL=1,NG_for_psi 
01181    igL1=LG0(1,igL,bnq+iq-1)
01182    igL2=LG0(2,igL,bnq+iq-1)
01183    igL3=LG0(3,igL,bnq+iq-1)
01184    if(igL1==0.and.igL2==0.and.igL3==0)then 
01185     No_G_0=igL  
01186     !20180822
01187     length_qg(igL)=1.0d0 
01188     atten_factor(igL)=0.0d0 
01189    else
01190     qgL(:)=(q1+dble(igL1))*b1(:)+(q2+dble(igL2))*b2(:)+(q3+dble(igL3))*b3(:)
01191     qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
01192     qgL1=dsqrt(qgL2) 
01193     length_qg(igL)=qgL1
01194     atten_factor(igL)=1.0d0-dcos(qgL1*Rc)  
01195    endif 
01196   enddo!igL 
01197   do ib=1,Ncalc!Mt 
01198    rho(:,:,:)=0.0D0 
01199 !$OMP PARALLEL PRIVATE(ik,shift_G,ikq,C0_K,C0_KmQ,jb,wfunc,fftwk,rho_tmp,ikir,iop,ikqir,ikqop) 
01200    allocate(fftwk(Nl123*2));fftwk=0.0d0 
01201    allocate(wfunc(Nl123*2));wfunc=0.0d0 
01202    allocate(rho_tmp(NG_for_psi));rho_tmp=0.0d0
01203    allocate(C0_K(NTG));C0_K=0.0d0
01204    allocate(C0_KmQ(NTG));C0_KmQ=0.0d0
01205 !$OMP DO
01206    do ikir=1,Nk_irr 
01207     !
01208     ik=numMK(ikir) 
01209     iop=numrot(ik) 
01210     !
01211     do jb=1,Nb(ik)!Mb 
01212      !
01213      call make_C0_for_given_band(NTG,trs(ik),NG0(ik),KG0(1,1,ik),RW(1,ik),rginv(1,1,iop),pg(1,iop),nnp,&
01214      L1,L2,L3,packing(-L1,-L2,-L3,ikir),CIR(1,jb+Ns(ik),ikir),C0_K(1)) 
01215      !
01216      shift_G(:)=0
01217      ikq=ik 
01218      !
01219      ikqir=numirr(ikq)
01220      ikqop=numrot(ikq) 
01221      !
01222      call make_C0_for_given_band(NTG,trs(ikq),NG0(ikq),KG0(1,1,ikq),RW(1,ikq),rginv(1,1,ikqop),pg(1,ikqop),nnp,&
01223      L1,L2,L3,packing(-L1,-L2,-L3,ikqir),CIR(1,ib,ikqir),C0_KmQ(1)) 
01224      !
01225      call calc_InterStateMatrix(NTK,NTG,NG0(1),KG0(1,1,1),C0_KmQ(1),C0_K(1),ikq,ik,nwx2,nwy2,nwz2,&
01226      nfft1,nfft2,Nl123,wfunc(1),fftwk(1),fs,LG0(1,1,bnq+iq-1),NG_for_psi,shift_G(1),rho_tmp(1))
01227      !
01228      rho(:,jb,ikir)=rho_tmp(:)
01229      !
01230     enddo!jb 
01231    enddo!ik  
01232 !$OMP END DO
01233    deallocate(fftwk,wfunc,rho_tmp,C0_K,C0_KmQ)
01234 !$OMP END PARALLEL 
01235 !$OMP PARALLEL PRIVATE(ikir,ik,ikq,jb,kb,SUM_CMPX,igL) 
01236 !$OMP DO
01237    do ikir=1,Nk_irr 
01238     !
01239     ik=numMK(ikir)
01240     !
01241     ikq=ik 
01242     !
01243     do jb=1,Nb(ik)!Mb 
01244      do kb=1,Nb(ik)!Mb 
01245       SUM_CMPX=0.0d0
01246       do igL=1,NG_for_psi 
01247        !if(igL.ne.No_G_0)then 
01248        !
01249        SUM_CMPX=SUM_CMPX+fbk(ib,ikq)*rho(igL,jb,ikir)*CONJG(rho(igL,kb,ikir))/((length_qg(igL))**2)*atten_factor(igL)
01250        !
01251        !endif 
01252       enddo!igL 
01253       !
01254       pSX(jb,kb,ikir)=pSX(jb,kb,ikir)+SUM_CMPX 
01255       !
01256      enddo!kb 
01257     enddo!jb 
01258    enddo!ik 
01259 !$OMP END DO
01260 !$OMP END PARALLEL 
01261   enddo ! ib 
01262   !write(file_id,*)'FINISH iq',iq
01263   !if(myrank.eq.0) write(6,*)'FINISH iq',iq 
01264   write(6,*)'FINISH iq myrank',iq,myrank  
01265   deallocate(length_qg,atten_factor,rho) 
01266  endif 
01267 enddo!iq 
01268 !
01269 call MPI_BARRIER(comm,ierr)
01270 !write(file_id,*)'I finished pSX calc' 
01271 call MPI_REDUCE(pSX,SXirr,Mb*Mb*Nk_irr,MPI_DOUBLE_COMPLEX,MPI_SUM,0,comm,ierr)
01272 !
01273 if(myrank.eq.0)then 
01274  SXirr(:,:,:)=2.0d0*tpi*SXirr(:,:,:)/dble(NTQ)/VOLUME 
01275  !
01276  !<head contribution> iq=NTQ
01277  !
01278  !(i) Spencer-Alabi
01279  !
01280  chead=(tpi/dble(NTQ)/VOLUME)*Rc*Rc   
01281  !
01282  !(ii) Hybertsen-Louie
01283  !
01284  !qsz=(6.0D0*(pi**2)/dble(NTQ)/dble(VOLUME))**(1.0D0/3.0D0)  
01285  !chead=(2.0D0/pi)*qsz 
01286  !
01287  write(6,*)'correction_head=',chead   
01288  !
01289  allocate(SX_DIV(Mb,Mb,Nk_irr));SX_DIV(:,:,:)=0.0d0 
01290  do ikir=1,Nk_irr 
01291   ik=numMK(ikir)
01292   do jb=1,Nb(ik)!Mb 
01293    SX_DIV(jb,jb,ikir)=fbk(jb+Ns(ik),ik)*chead 
01294   enddo!jb 
01295  enddo!ik 
01296  SXirr=SXirr+SX_DIV!fGW
01297  !--
01298  !SXirr=SXirr!cGW
01299  !--
01300  allocate(SX(Mb,Mb,NTK));SX(:,:,:)=0.0d0 
01301  do ik=1,NTK 
01302   ikir=numirr(ik) 
01303   if(trs(ik)==1) then 
01304    SX(:,:,ik)=SXirr(:,:,ikir) 
01305   elseif(trs(ik)==-1) then 
01306    do ib=1,Nb(ik)!Mb
01307     do jb=1,Nb(ik)!Mb
01308      SX(ib,jb,ik)=SXirr(jb,ib,ikir) 
01309     enddo!jb
01310    enddo!ib 
01311   endif 
01312  enddo!ik 
01313  !
01314  allocate(pf(-Na1:Na1,-Na2:Na2,-Na3:Na3,NTK));pf=0.0d0     
01315  do ik=1,NTK 
01316   do ia3=-Na3,Na3 
01317    do ia2=-Na2,Na2 
01318     do ia1=-Na1,Na1 
01319      phase=tpi*(SK0(1,ik)*dble(ia1)+SK0(2,ik)*dble(ia2)+SK0(3,ik)*dble(ia3)) 
01320      pf(ia1,ia2,ia3,ik)=exp(-ci*phase) 
01321     enddo  
01322    enddo  
01323   enddo  
01324  enddo  
01325  write(6,*)'#finish make pf'
01326  !
01327  allocate(MAT_SX_R(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3));MAT_SX_R(:,:,:,:,:)=0.0d0 
01328  do ia3=-Na3,Na3
01329   do ia2=-Na2,Na2
01330    do ia1=-Na1,Na1
01331     do iw=1,NWF
01332      do jw=1,NWF 
01333       SUM_CMPX=0.0D0 
01334       if(.true.)then!.true.=off-diag, .false.=diag
01335        do ik=1,NTK 
01336         do ib=1,Nb(ik) 
01337          do jb=1,Nb(ik) 
01338           SUM_CMPX=SUM_CMPX+CONJG(UNT(ib,iw,ik))*SX(ib,jb,ik)*UNT(jb,jw,ik)*pf(ia1,ia2,ia3,ik) 
01339          enddo 
01340         enddo 
01341        enddo 
01342       else 
01343        do ik=1,NTK 
01344         do ib=1,Nb(ik) 
01345          SUM_CMPX=SUM_CMPX+CONJG(UNT(ib,iw,ik))*SX(ib,ib,ik)*UNT(ib,jw,ik)*pf(ia1,ia2,ia3,ik) 
01346         enddo 
01347        enddo 
01348       endif 
01349       MAT_SX_R(iw,jw,ia1,ia2,ia3)=SUM_CMPX/DBLE(NTK)
01350      enddo 
01351     enddo 
01352    enddo 
01353   enddo 
01354  enddo 
01355  deallocate(SX,pf) 
01356  write(6,*)'finish MAT_SX_R' 
01357  !
01358  !write(5003) MAT_SX_R 
01359  !
01360  write(6,*) 
01361  write(6,*)'================'
01362  write(6,*)' MAT_SX_R in eV '
01363  write(6,*)'================'
01364  write(6,*) 
01365  do ia1=-Na1,Na1 
01366   do ia2=-Na2,Na2 
01367    do ia3=-Na3,Na3 
01368     write(6,*) ia1,ia2,ia3           
01369     do iw=1,NWF
01370      write(6,'(300F15.8)')(dble(MAT_SX_R(iw,jw,ia1,ia2,ia3)*au),jw=1,NWF)
01371     enddo!iw  
01372     write(6,*)  
01373    enddo!ia3       
01374   enddo!ia2 
01375  enddo!ia1 
01376  !
01377 endif!myrank.eq.0
01378 !
01379 !#############
01380 !  OUTPUT SX 
01381 !#############
01382 !
01383 if(myrank.eq.0)then 
01384  ierr=CHDIR("./dir-gw") 
01385  call system('pwd') 
01386  !
01387  !OPEN(154,FILE='./dat.energy_vs_sx') 
01388  !
01389  OPEN(154,FILE='./dat.energy_vs_sx') 
01390  rewind(154) 
01391  do ikir=1,Nk_irr 
01392   ik=numMK(ikir)
01393   do jb=1,Nb(ik)!Mb 
01394    en=(E_EIGI(jb+Ns(ik),ikir)-FermiEnergy)*au 
01395    write(154,'(3F20.10)') en,SXirr(jb,jb,ikir)*au 
01396   enddo!jb 
01397  enddo!ikir 
01398  close(154)
01399  !
01400  deallocate(SXirr) 
01401  !
01402  !OPEN(155,FILE='./dat.sx_mat_r') 
01403  !
01404  OPEN(155,FILE='./dat.sx_mat_r') 
01405  rewind(155)
01406  write(155,'(a)')'#sx_mat_r'
01407  write(155,'(a)')'#1:R1, 2:R2, 3:R3 (lattice vector)'
01408  write(155,'(a)')'#1:i, 2:j, 3:Re(Sx_ij) [eV], 4:Im(Sx_ij) [eV]' 
01409  do ia1=-Na1,Na1 
01410   do ia2=-Na2,Na2 
01411    do ia3=-Na3,Na3
01412     write(155,*) ia1,ia2,ia3           
01413     do iw=1,NWF
01414      do jw=1,NWF 
01415       write(155,'(i5,i5,2f20.10)') iw,jw,MAT_SX_R(iw,jw,ia1,ia2,ia3)*au 
01416      enddo!jw       
01417     enddo!iw       
01418     write(155,*) 
01419    enddo!ia3
01420   enddo!ia2
01421  enddo!ia1
01422  close(155) 
01423  ierr=CHDIR("..") 
01424  call system('pwd') 
01425 endif!myrank.eq.0
01426 !
01427 !#######
01428 !  VXC
01429 !#######
01430 !
01431 if(myrank.eq.0)then
01432  !
01433  !write(file_id,'(a)')'## Vxc calc start ##'
01434  !
01435  write(6,*)'## Vxc calc start ##' 
01436  write(6,*) 
01437  !
01438  !
01439  !OPEN(151,FILE='./dat.vxc',FORM='unformatted') 
01440  !
01441  OPEN(151,FILE='./dat.vxc',FORM='unformatted') 
01442  REWIND(151)       
01443  read(151) header 
01444  read(151) aa,nsnd 
01445  if(nsnd>1)then
01446   write(6,*)'nspin*ndens>1; the program not supported; then stop'
01447   stop
01448  endif 
01449  read(151) nrx2,nry2,nrz2; write(6,*)'nrx2,nry2,nrz2=',nrx2,nry2,nrz2 
01450  allocate(vxcr(nrx2,nry2,nrz2));vxcr(:,:,:)=0.0d0 
01451  read(151) vxcr(:,:,:) 
01452  close(151)  
01453  write(6,*)'FINISH REDING VXCR'
01454  !
01455  !fft
01456  !
01457  nfft1=nrx2+1; nfft2=nry2+1; nfft3=nrz2+1; Nl123=nfft1*nfft2*nfft3 
01458  call fft3_init(nrx2,nry2,nrz2,nfft1,nfft2,nfft3,fs) 
01459  allocate(fftwk(Nl123*2),stat=err) 
01460  allocate(wfunc(Nl123*2),stat=err) 
01461  !
01462  allocate(VXCirr(Mb,Mb,Nk_irr));VXCirr(:,:,:)=0.0d0 
01463  allocate(MAT_VXC(Mb,Mb));MAT_VXC=0.0d0             
01464  do ikir=1,Nk_irr 
01465   ik=numMK(ikir)
01466   MAT_VXC(:,:)=0.0d0 
01467   call calc_MAT_VXC(Nk_irr,Mb,NTG,NGI(1),KGI(1,1,1),vxcr(1,1,1),CIR(1,1+Ns(ik),ikir),ikir,&
01468   nrx2,nry2,nrz2,nfft1,nfft2,Nl123,wfunc(1),fftwk(1),fs,MAT_VXC(1,1))
01469   VXCirr(:,:,ikir)=MAT_VXC(:,:)
01470   write(6,*)'ikir=',ikir 
01471  enddo!ikir 
01472  deallocate(MAT_VXC,vxcr) 
01473  !
01474  allocate(VXC(Mb,Mb,NTK));VXC(:,:,:)=0.0d0 
01475  do ik=1,NTK 
01476   !
01477   ikir=numirr(ik) 
01478   !
01479   if(trs(ik)==1)then 
01480    VXC(:,:,ik)=VXCirr(:,:,ikir) 
01481   elseif(trs(ik)==-1)then 
01482    do ib=1,Nb(ik)
01483     do jb=1,Nb(ik) 
01484      VXC(ib,jb,ik)=VXCirr(jb,ib,ikir) 
01485     enddo!jb 
01486    enddo!ib 
01487   endif 
01488  enddo!ik 
01489  write(6,*)'#finish make VXC'
01490  !
01491  allocate(pf(-Na1:Na1,-Na2:Na2,-Na3:Na3,NTK));pf=0.0d0     
01492  do ik=1,NTK 
01493   do ia3=-Na3,Na3 
01494    do ia2=-Na2,Na2 
01495     do ia1=-Na1,Na1 
01496      phase=tpi*(SK0(1,ik)*dble(ia1)+SK0(2,ik)*dble(ia2)+SK0(3,ik)*dble(ia3)) 
01497      pf(ia1,ia2,ia3,ik)=exp(-ci*phase) 
01498     enddo  
01499    enddo  
01500   enddo  
01501  enddo  
01502  write(6,*)'#finish make pf'
01503  !
01504  allocate(MAT_VXC_R(NWF,NWF,-Na1:Na1,-Na2:Na2,-Na3:Na3));MAT_VXC_R(:,:,:,:,:)=0.0d0 
01505  do ia3=-Na3,Na3
01506   do ia2=-Na2,Na2
01507    do ia1=-Na1,Na1
01508     do iw=1,NWF
01509      do jw=1,NWF 
01510      SUM_CMPX=0.0D0 
01511      if(.true.)then!.true.=off-diag, .false.=diag
01512       do ik=1,NTK 
01513        do ib=1,Nb(ik) 
01514         do jb=1,Nb(ik) 
01515          SUM_CMPX=SUM_CMPX+CONJG(UNT(ib,iw,ik))*VXC(ib,jb,ik)*UNT(jb,jw,ik)*pf(ia1,ia2,ia3,ik) 
01516         enddo 
01517        enddo 
01518       enddo 
01519      else 
01520       do ik=1,NTK 
01521        do ib=1,Nb(ik) 
01522         SUM_CMPX=SUM_CMPX+CONJG(UNT(ib,iw,ik))*VXC(ib,ib,ik)*UNT(ib,jw,ik)*pf(ia1,ia2,ia3,ik) 
01523        enddo 
01524       enddo 
01525      endif 
01526      MAT_VXC_R(iw,jw,ia1,ia2,ia3)=SUM_CMPX/DBLE(NTK)
01527      enddo 
01528     enddo 
01529    enddo 
01530   enddo 
01531  enddo 
01532  deallocate(VXC,pf) 
01533  !
01534  write(6,*) 
01535  write(6,*)'================='
01536  write(6,*)' MAT_VXC_R in eV '
01537  write(6,*)'================='
01538  write(6,*) 
01539  do ia1=-Na1,Na1 
01540   do ia2=-Na2,Na2 
01541    do ia3=-Na3,Na3 
01542     write(6,*) ia1,ia2,ia3           
01543     do iw=1,NWF
01544      write(6,'(300F15.8)')(dble(MAT_VXC_R(iw,jw,ia1,ia2,ia3)*au),jw=1,NWF)
01545     enddo!iw  
01546     write(6,*)  
01547    enddo!ia3       
01548   enddo!ia2 
01549  enddo!ia1 
01550  write(6,*)'finish MAT_VXC_R' 
01551  !
01552  !write(5002) MAT_VXC_R 
01553  !
01554 endif!myrank.eq.0
01555 !
01556 !##############
01557 !  OUTPUT VXC 
01558 !##############
01559 !
01560 if(myrank.eq.0)then 
01561  ierr=CHDIR("./dir-gw") 
01562  call system('pwd') 
01563  !
01564  !OPEN(152,W,FILE='energy_vs_vxc') 
01565  !
01566  OPEN(152,FILE='./dat.energy_vs_vxc') 
01567  rewind(152) 
01568  do ikir=1,Nk_irr 
01569   ik=numMK(ikir)
01570   do jb=1,Nb(ik)!Mb 
01571    en=(E_EIGI(jb+Ns(ik),ikir)-FermiEnergy)*au 
01572    write(152,'(3f20.10)') en,VXCirr(jb,jb,ikir)*au 
01573   enddo!jb 
01574  enddo!ikir 
01575  close(152)
01576  !
01577  deallocate(VXCirr) 
01578  !
01579  !OPEN(153,W,FILE='vxc_mat_r') 
01580  !
01581  OPEN(153,FILE='./dat.vxc_mat_r') 
01582  rewind(153)
01583  write(153,'(a)')'#vxc_mat_r'
01584  write(153,'(a)')'#1:R1, 2:R2, 3:R3 (lattice vector)'
01585  write(153,'(a)')'#1:i, 2:j, 3:Re(vxc_ij) [eV], 4:Im(vxc_ij) [eV]' 
01586  do ia1=-Na1,Na1 
01587   do ia2=-Na2,Na2 
01588    do ia3=-Na3,Na3 
01589     write(153,*) ia1,ia2,ia3           
01590     do iw=1,NWF
01591      do jw=1,NWF 
01592       write(153,'(i5,i5,2f20.10)') iw,jw,MAT_VXC_R(iw,jw,ia1,ia2,ia3)*au 
01593      enddo!jw  
01594     enddo!iw  
01595     write(153,*)  
01596    enddo!ia3       
01597   enddo!ia2 
01598  enddo!ia1  
01599  close(153) 
01600  ierr=CHDIR("..") 
01601  call system('pwd') 
01602 endif!myrank.eq.0
01603 !
01604 !######################
01605 !  GW-DOS CALCULATION
01606 !######################
01607 !
01608 if(myrank.eq.0)then 
01609   write(6,'(a)')'## gw-dos calc start ##'
01610   allocate(ksdos(nsgm)); ksdos=0.0d0 
01611   allocate(gwdos(nsgm)); gwdos=0.0d0  
01612   allocate(gw_sigma_dos(nsgm)); gw_sigma_dos=0.0d0  
01613   ! 
01614   call calc_ksdos(NWF,NTK,nsgm,Na1,Na2,Na3,nkb1,nkb2,nkb3,idlt,dmna,dmnr,FermiEnergy,a1(1),a2(1),a3(1),&
01615   b1(1),b2(1),b3(1),sgmw(1),SK0(1,1),H_MAT_R(1,1,-Na1,-Na2,-Na3),ksdos(1))
01616   !
01617   call calc_gwdos(NWF,NTK,nsgm,Na1,Na2,Na3,nkb1,nkb2,nkb3,idlt,dmna,dmnr,FermiEnergy,a1(1),a2(1),a3(1),&
01618   b1(1),b2(1),b3(1),sgmw(1),SK0(1,1),H_MAT_R(1,1,-Na1,-Na2,-Na3),MAT_VXC_R(1,1,-Na1,-Na2,-Na3),&
01619   MAT_SX_R(1,1,-Na1,-Na2,-Na3),MAT_SC_R(1,1,-Na1,-Na2,-Na3,1),shift_value,gwdos(1),gw_sigma_dos(1))
01620   !
01621   write(6,'(a)')'finish gw-dos calc' 
01622 endif!myrank.eq.0 
01623 !
01624 !##############
01625 !  OUTPUT DOS 
01626 !##############
01627 !
01628 if(myrank.eq.0)then 
01629  ierr=CHDIR("./dir-gw") 
01630  call system('pwd') 
01631  !
01632  !OPEN(159,W,FILE='dat.gwdos') 
01633  !
01634  OPEN(159,FILE='./dat.gwdos') 
01635  rewind(159) 
01636  do ie=1,nsgm 
01637   write(159,'(2f20.10)') sgmw(ie)*au,gwdos(ie) 
01638  enddo 
01639  close(159) 
01640  !
01641  !OPEN(160,W,FILE='dat.ksdos') 
01642  !
01643  OPEN(160,FILE='./dat.ksdos') 
01644  rewind(160)
01645  do ie=1,nsgm
01646   write(160,'(2f20.10)') sgmw(ie)*au,ksdos(ie)
01647  enddo 
01648  close(160) 
01649  !
01650  !OPEN(170,W,FILE='dat.gw_sigma_dos') 
01651  !
01652  OPEN(170,FILE='./dat.gw_sigma_dos') 
01653  rewind(170)
01654  do ie=1,nsgm
01655   write(170,'(3f20.10)') sgmw(ie)*au,gw_sigma_dos(ie)
01656  enddo 
01657  close(170) 
01658  !
01659  ierr=CHDIR("..") 
01660  call system('pwd') 
01661 endif!myrank.eq.0 
01662 !
01663 !######################
01664 !  GW-AKW CALCULATION
01665 !######################
01666 !
01667 if(myrank.eq.0)then
01668   write(6,'(a)')'## gw-akw calc start ##'
01669   NSK_BAND_DISP=Ndiv*(N_sym_points-1)+1
01670   allocate(SK_BAND_DISP(3,NSK_BAND_DISP)); SK_BAND_DISP(:,:)=0.0d0 
01671   call makekpts(Ndiv,N_sym_points,NSK_BAND_DISP,SK_sym_pts(1,1),SK_BAND_DISP(1,1))
01672   !
01673   !write(5001) H_MAT_R 
01674   !H_MAT_R=H_MAT_R-MAT_VXC_R-MAT_SX_R  
01675   !
01676   allocate(kdata(NSK_BAND_DISP)); kdata=0.0d0 
01677   allocate(EKS(NWF,NSK_BAND_DISP)); EKS=0.0d0 
01678   call calc_band_disp(Ndiv,N_sym_points,NTK,NSK_BAND_DISP,Na1,Na2,Na3,NWF,SK_BAND_DISP(1,1),&
01679   H_MAT_R(1,1,-Na1,-Na2,-Na3),nkb1,nkb2,nkb3,a1(1),a2(1),a3(1),b1(1),b2(1),b3(1),kdata,EKS(1,1))  
01680   !
01681   allocate(gwakw(NSK_BAND_DISP,nsgm)); gwakw=0.0d0
01682   allocate(gw_sigma_kw(NSK_BAND_DISP,nsgm)); gw_sigma_kw=0.0d0
01683   call calc_gwakw(NWF,NTK,nsgm,Na1,Na2,Na3,nkb1,nkb2,nkb3,NSK_BAND_DISP,idlt,shift_value,&
01684   SK_BAND_DISP(1,1),a1(1),a2(1),a3(1),sgmw(1),H_MAT_R(1,1,-Na1,-Na2,-Na3),MAT_VXC_R(1,1,-Na1,-Na2,-Na3),&
01685   MAT_SX_R(1,1,-Na1,-Na2,-Na3),MAT_SC_R(1,1,-Na1,-Na2,-Na3,1),gwakw(1,1),gw_sigma_kw(1,1))            
01686   !
01687   write(6,'(a)')'finish gw-akw calc' 
01688 endif!myrank.eq.0
01689 !
01690 !##############
01691 !  OUTPUT AKW 
01692 !##############
01693 !
01694 if(myrank.eq.0)then 
01695  ierr=CHDIR("./dir-gw") 
01696  call system('pwd') 
01697  !
01698  !OPEN(114,W,FILE='dat.iband') 
01699  !
01700  OPEN(114,FILE='./dat.iband') 
01701  write(114,'(a)')'#Wannier interpolaed band'
01702  write(114,'(a)')'#1:k, 2:Energy [eV]' 
01703  REVERSE=.TRUE.        
01704  do ib=1,NWF 
01705   if(REVERSE)then 
01706    do ik=1,NSK_BAND_DISP                     
01707     write(114,'(2f20.10)') kdata(ik),EKS(ib,ik)*au
01708    enddo!ik        
01709    REVERSE=.FALSE.        
01710   else         
01711    do ik=NSK_BAND_DISP,1,-1          
01712     write(114,'(2f20.10)') kdata(ik),EKS(ib,ik)*au
01713    enddo!ik        
01714    REVERSE=.TRUE.        
01715   endif!REVERSE                   
01716  enddo!ib 
01717  close(114) 
01718  ! 
01719  !OPEN(161,W,FILE='dat.gwakw') 
01720  !
01721  OPEN(161,FILE='./dat.gwakw') 
01722  rewind(161) 
01723  do ie=1,nsgm 
01724   do ik=1,NSK_BAND_DISP 
01725    if(gwakw(ik,ie)>=500.0d0) gwakw(ik,ie)=500.0d0 
01726    write(161,'(3f20.10)') kdata(ik),sgmw(ie)*au,gwakw(ik,ie) 
01727   enddo 
01728   write(161,*) 
01729  enddo 
01730  close(161) 
01731  ! 
01732  !OPEN(171,W,FILE='dat.gw_sigma_kw') 
01733  !
01734  OPEN(171,FILE='./dat.gw_sigma_kw') 
01735  rewind(171) 
01736  do ie=1,nsgm 
01737   do ik=1,NSK_BAND_DISP 
01738    write(171,'(3f20.10)') kdata(ik),sgmw(ie)*au,gw_sigma_kw(ik,ie) 
01739   enddo 
01740   write(171,*) 
01741  enddo 
01742  close(171) 
01743  ! 
01744  ierr=CHDIR("..") 
01745  call system('pwd') 
01746 endif!myrank.eq.0 
01747 !
01748 !##########
01749 !  FINISH 
01750 !##########
01751 !
01752 call MPI_BARRIER(comm,ierr)
01753 end_time=MPI_Wtime()
01754 diff_time=end_time-start_time 
01755 if(myrank.eq.0) then 
01756  call system('pwd') 
01757  write(6,*)'#TOTAL TIME=',diff_time 
01758 endif 
01759 !
01760 call MPI_FINALIZE(ierr)
01761 STOP 
01762 END