1subroutine mkppse(iband,jband,ikk,omegap2,xkf,lossfn,fsum, &
2& vol,pi,nwpt,wmax,nbcore,nbocc,ncband,ngkpt,natom,xred,projwf,nlmn, &
3& pwmatel,tpwmatel, &
4& kg,kgq,enrgy,enrgyq,cg,cgq,npwt,npwtq,bantot,bantotq,ncg,ncgq, &
5& indxkpw,indxkpwq,indxkbnd,indxkbndq,indxkcg,indxkcgq,npwarr,npwarrq, &
6& kpt,kptq,nkpt,nkptq,nqpt,nsymk,nsymkq,symk,symkq,nsym,nsymq,symrel,syminv, &
7& ihlf,ihlfq,lvtrans,lvtransq,bmet,blat,ipaw, &
8& ipwndx,npwndx,ntpwndx,napwndx, &
9& npwc,npwx,invpw2ndx,pwsymndx,iqsymndx,invpwndx,pwcmin,pwcmax, &
10& igmx,igmn,igndx,igndxq,ikndx,ikndxq,iqndx,isymndx,isymndxq,npw,npwq, &
11& nband,nbandq,nsppol,shiftk,shiftkq,zz,zzpp,zzcorr,cse,csepp,csecorr,xse)
12implicit none
13integer :: iband,jband,ikk(3),nwpt,nbcore,nbocc,ncband,ngkpt(3),natom,nlmn
14integer :: igmn(3),igmx(3)
15integer :: npwt,npwtq,bantot,bantotq,ncg,ncgq,nkpt,nkptq,nqpt,nsym,nsymq,nsppol
16integer :: npw,npwq,ipw1,npwndx,ntpwndx,napwndx,ipwa,ipaw
17double precision :: vol,pi,wmax,xred(3,natom)
18double precision :: omegap2,xkf,ppfct
19double complex :: lossfn(nwpt,nqpt,ntpwndx)
20double complex :: projwf(natom,nlmn,nkpt,ncband)
21integer :: kg(3,npwt),kgq(3,npwtq)
22double precision :: enrgy(bantot),enrgyq(bantotq)
23double complex :: cg(ncg),cgq(ncgq)
24integer :: indxkpw(nkpt),indxkpwq(nkptq),indxkbnd(nkpt),indxkbndq(nkptq)
25integer :: indxkcg(nkpt),indxkcgq(nkptq),npwarr(nkpt),npwarrq(nkptq)
26double precision :: kpt(3,nkpt),kptq(3,nkptq),shiftk(3),shiftkq(3)
27integer :: nsymk(nkpt),nsymkq(nkptq),symk(nkpt,nsym*2),symkq(nkptq,nsymq*2)
28integer :: symrel(3,3,nsym),syminv(3,3,nsymq)
29integer :: lvtrans(3,ngkpt(1),ngkpt(2),ngkpt(3))
30integer :: lvtransq(3,ngkpt(1),ngkpt(2),ngkpt(3))
31integer :: ihlf(nkpt),ihlfq(nkptq)
32double precision :: bmet(3,3),blat(3,3)
33integer :: ipwndx(2,napwndx)
34integer :: igndx(nkpt,igmn(1):igmx(1),igmn(2):igmx(2),igmn(3):igmx(3))
35integer :: igndxq(nkptq,igmn(1):igmx(1),igmn(2):igmx(2),igmn(3):igmx(3))
36integer :: ikndx(ngkpt(1),ngkpt(2),ngkpt(3)),ikndxq(ngkpt(1),ngkpt(2),ngkpt(3))
37integer :: iqndx(ngkpt(1),ngkpt(2),ngkpt(3))
38integer :: isymndx(ngkpt(1),ngkpt(2),ngkpt(3)),isymndxq(ngkpt(1),ngkpt(2),ngkpt(3))
39integer :: nband(nkpt*nsppol),nbandq(nkptq*nsppol)
40double complex :: sei,seipw,seiold,seibc,cseold,csebc
41double complex :: xse1,xse2
42double complex :: cse(nbcore+1:ncband),cse1(nbcore+1:ncband),cse2(nbcore+1:ncband)
43double complex :: csepp(nbcore+1:ncband),cse1pp(nbcore+1:ncband)
44double complex :: csecorr(nbcore+1:ncband),cse1corr(nbcore+1:ncband)
45double complex :: zz(nbcore+1:ncband),zz1(nbcore+1:ncband),zz2(nbcore+1:ncband)
46double complex :: zzpp(nbcore+1:ncband),zz1pp(nbcore+1:ncband)
47double complex :: zzcorr(nbcore+1:ncband),zz1corr(nbcore+1:ncband)
48double precision :: xse
49integer :: ii,jj,kk,ix,iy,iz,iqq(3),iqqp(3),jka(3),jkb(3),jkk(3),iskip,isign,iloss,insd
50integer :: ikpt,ikptq,iks(3),ikslf1(3),ikslf2(3)
51integer :: pwcmin(3),pwcmax(3)
52integer :: igg(3),igpw(3),igh(3),igg0(3),icenter,isym,isymq,ibp,iqpt,iqsym,iw,ie,je,ies
53double precision :: xck(3),xckq(3)
54double complex :: cmatel,cmatel2,amatel,amatel2,vqmat2(ngkpt(1),ngkpt(2),ngkpt(3))
55double complex :: vfactor(ngkpt(1),ngkpt(2),ngkpt(3)),vfv(8),lossv(8)
56double precision :: vq2(8)
57double precision :: omega(ngkpt(1),ngkpt(2),ngkpt(3))
58double precision :: whi,wlo,wwhi,wwlo,ww,enval(8),dw,eshift,wwq
59double precision :: abr,rlr
60double complex :: tseincr,tseincrx,dse,dcse
61integer :: iwh,iwl,ibmin,ibmax
62double precision :: vq(ngkpt(1),ngkpt(2),ngkpt(3)),qq(3),qp(3),qq2,qp2,qs(3),xk(3),xkmq(3),ek,ekmq,ekpw,xpw(3)
63double precision :: stvec(ngkpt(3)),stvec2(ngkpt(3))
64integer :: npwc,npwx,iipw,jjpw
65integer :: iqsymndx(ngkpt(1),ngkpt(2),ngkpt(3)),invpw2ndx(npwx,npwx), &
66& invpwndx(pwcmin(1):pwcmax(1),pwcmin(2):pwcmax(2),pwcmin(3):pwcmax(3))
67integer :: pwsymndx(npwc,2*nsym)
68integer :: ipw2,jw
69double precision :: fsum(npwc,nqpt)
70double precision :: eps1(nbcore+1:ncband),eps2,eps
71double precision :: eval,brd,esprd
72double complex :: cenrgy
73double precision :: gfo,gamma(3),pln(3),dist(3)
74double complex :: wint(nbcore+1:ncband),wint0(nbcore+1:ncband)
75double complex :: w2int(nbcore+1:ncband),w2int0(nbcore+1:ncband)
76double complex :: wppint(nbcore+1:ncband),wppint0(nbcore+1:ncband)
77double complex :: wpp2int(nbcore+1:ncband),wpp2int0(nbcore+1:ncband)
78double complex :: cmgamma2,gterm
79double complex :: pwmatel(nlmn,nlmn,npwx,ngkpt(1),ngkpt(2),ngkpt(3)), &
80&                tpwmatel(nlmn,nlmn,npwx,ngkpt(1),ngkpt(2),ngkpt(3))
81integer :: iqsing
82integer :: idum
83logical :: lx,lc,lcorr
84double precision :: xdum
85double precision :: linterp
86external linterp
87
88abr=1.d-6
89rlr=1.d-6
90sei=0.d0
91xse=0.d0
92cse=(0.d0,0.d0)
93csecorr=(0.d0,0.d0)
94zz=(0.d0,0.d0)
95igg0=(/0,0,0/)
96dw=wmax/dble(nwpt)
97ikpt=ikndx(ikk(1),ikk(2),ikk(3))
98isym=isymndx(ikk(1),ikk(2),ikk(3))
99ek=enrgy(indxkbnd(ikpt)+iband)
100
101gamma=(blat(1,:)+blat(2,:)+blat(3,:))/2.d0
102do ii=1,3
103  jj=mod(ii,3)+1
104  kk=mod(ii+1,3)+1
105  pln(1)=blat(jj,2)*blat(kk,3)-blat(jj,3)*blat(kk,2)
106  pln(2)=-blat(jj,1)*blat(kk,3)+blat(jj,3)*blat(kk,1)
107  pln(3)=blat(jj,1)*blat(kk,2)-blat(jj,2)*blat(kk,1)
108  dist(ii)=dot_product(gamma,pln)/sqrt(dot_product(pln,pln))
109enddo
110gfo=minval(dist)
111
112!write(6,'(a)') '      finding contibutions from band:      plane waves:'
113do ibp=nbcore+1,ncband
114!do ibp=1,1
115  xse2=xse
116  cse2=cse
117  if (ibp.gt.nbocc) then
118    isign=-1
119  else
120    isign=1
121  endif
122  if (ibp.ge.ibmin.and.ibp.le.ibmax) then
123    iloss=1
124  else
125    iloss=0
126  endif
127  do ipw1=1,npwarr(ikpt)
128!  do ipw1=1,1
129    ipw2=ipw1
130    igpw=kg(:,indxkpw(ikpt)+ipw1)
131    if (igpw(1).ge.pwcmin(1).and.igpw(1).le.pwcmax(1).and.  &
132&       igpw(2).ge.pwcmin(2).and.igpw(2).le.pwcmax(2).and.  &
133&       igpw(3).ge.pwcmin(3).and.igpw(3).le.pwcmax(3)) then
134      ipwa=invpwndx(igpw(1),igpw(2),igpw(3))
135    else
136      ipwa=0
137    endif
138    xpw=dble(igpw)+dble(ikk)/dble(ngkpt)
139    ekpw=0.d0
140    do ii=1,3
141    do jj=1,3
142      ekpw=ekpw+xpw(ii)*bmet(ii,jj)*xpw(jj)
143    enddo
144    enddo
145    if (ekpw.gt.8.d0) cycle
146    do ix=1,ngkpt(1)
147    do iy=1,ngkpt(2)
148    do iz=1,ngkpt(3)
149!    do ix=2,2
150!    do iy=2,2
151!    do iz=4,4
152      iqq=(/ix,iy,iz/)
153!write(6,*) '>>>>> iqq = ',iqq
154!write(6,'(a,3f10.5)') 'kpt = ',kpt(:,ikpt)
155      iks=iqq-ngkpt/2
156      jka=ikk-iks
157      jkk=mod(jka,ngkpt(ii))
158      do ii=1,3
159        if (jkk(ii).le.0) jkk(ii)=jkk(ii)+ngkpt(ii)
160      enddo
161      if (iks(1).eq.0.and.iks(2).eq.0.and.iks(3).eq.0) then
162!      if (.true.) then
163        ikptq=ikndxq(jkk(1),jkk(2),jkk(3))
164        jkb=nint(kptq(:,ikptq)*dble(ngkpt))+ngkpt/2
165        qq=kpt(:,ikpt)-kptq(:,ikptq)
166!write(6,'(a,3f10.5)') 'kptq = ',kptq(:,ikptq)
167      else
168        ikptq=ikndx(jkk(1),jkk(2),jkk(3))
169        jkb=nint(kpt(:,ikptq)*dble(ngkpt))+ngkpt/2
170        qq=kpt(:,ikpt)-kpt(:,ikptq)
171!write(6,'(a,3f10.5)') 'kptq = ',kptq(:,ikptq)
172      endif
173      igg=nint(dble(jkb-jka)/dble(ngkpt))
174      qq=qq+igg+igpw
175      qq2=0.d0
176      do ii=1,3
177      do jj=1,3
178        qq2=qq2+qq(ii)*bmet(ii,jj)*qq(jj)
179      enddo
180      enddo
181      vq(ix,iy,iz)=4.d0*pi/qq2
182      wwq=sqrt(omegap2+(xkf**2/3.d0)*qq2+qq2**2/4.d0)
183!write(6,'(a,3f10.5)') 'qq = ',qq
184!write(6,'(a,3f10.5)') 'qq2 = ',qq2
185!write(6,*) 'vq = ',vq(ix,iy,iz)
186!write(6,*) 'jka = ',jka
187!write(6,*) 'jkk = ',jkk
188!write(6,*) 'jkb = ',jkb
189!write(6,*) 'igg = ',igg
190!write(6,*) 'ikpt = ',ikpt
191!write(6,*) 'ikptq = ',ikptq
192      if (iks(1).eq.0.and.iks(2).eq.0.and.iks(3).eq.0) then
193!      if (.true.) then
194        call mkmatelX(ibp,iband,ikpt,ikptq,igpw,igg, &
195&         ncg,ncgq,nkpt,nkptq,npwt,npwtq,igmx,igmn,igndx,igndxq, &
196&         isym,isymq,symrel,syminv,nsym,nsymq,ihlf,ihlfq,kpt,kptq, &
197&         lvtrans(1:3,ikk(1),ikk(2),ikk(3)),lvtransq(1:3,ikk(1),ikk(2),ikk(3)), &
198&         cg,cgq,indxkcg,indxkcgq,indxkpw,indxkpwq,npwarr,npwarrq,kg,kgq, &
199&         cmatel)
200        if (ipaw.ne.0) then
201!*** need to calculate pwmatel for all plane waves
202!*** ipw1 indexes full plane wave set, not the local field plane wave set
203!*** pwmatel is calculated for
204!          call mkmatelP(pi,xred,natom,ibp,iband,ikpt,ikptq,qq,igg0,ngkpt, &
205!&               pwmatel(:,:,ipw1,ix,iy,iz),tpwmatel(:,:,ipw1,ix,iy,iz), &
206!&               projwf,nlmn,nkpt,nkptq,ncband,amatel)
207        else
208          amatel=0.d0
209        endif
210        if (ibp.gt.nbocc) then
211          omega(ix,iy,iz)=-enrgyq(indxkbndq(ikptq)+ibp)
212        else
213          omega(ix,iy,iz)=enrgyq(indxkbndq(ikptq)+ibp)
214        endif
215      else
216        call mkmatelX(ibp,iband,ikpt,ikptq,igpw,igg, &
217&         ncg,ncg,nkpt,nkpt,npwt,npwt,igmx,igmn,igndx,igndx, &
218&         isym,isymq,symrel,syminv,nsym,nsym,ihlf,ihlf,kpt,kpt, &
219&         lvtrans(1:3,ikk(1),ikk(2),ikk(3)),lvtrans(1:3,ikk(1),ikk(2),ikk(3)), &
220&         cg,cg,indxkcg,indxkcg,indxkpw,indxkpw,npwarr,npwarr,kg,kg, &
221&         cmatel)
222        if (ipaw.ne.0) then
223!          call mkmatelP(pi,xred,natom,ibp,iband,ikpt,ikptq,qq,igg0,ngkpt, &
224!&               pwmatel(:,:,ipw1,ix,iy,iz),tpwmatel(:,:,ipw1,ix,iy,iz), &
225!&               projwf,nlmn,nkpt,nkpt,ncband,amatel)
226        else
227          amatel=0.d0
228        endif
229        if (ibp.gt.nbocc) then
230          omega(ix,iy,iz)=-enrgy(indxkbnd(ikptq)+ibp)
231        else
232          omega(ix,iy,iz)=enrgy(indxkbnd(ikptq)+ibp)
233        endif
234      endif
235      if (iband.eq.jband) then
236        cmatel2=cmatel
237        amatel2=amatel
238      elseif (iks(1).eq.0.and.iks(2).eq.0.and.iks(3).eq.0) then
239!      if (.true.) then
240        call mkmatelX(ibp,jband,ikpt,ikptq,igpw,igg, &
241&         ncg,ncgq,nkpt,nkptq,npwt,npwtq,igmx,igmn,igndx,igndxq, &
242&         isym,isymq,symrel,syminv,nsym,nsymq,ihlf,ihlfq,kpt,kptq, &
243&         lvtrans(1:3,ikk(1),ikk(2),ikk(3)),lvtransq(1:3,ikk(1),ikk(2),ikk(3)), &
244&         cg,cgq,indxkcg,indxkcgq,indxkpw,indxkpwq,npwarr,npwarrq,kg,kgq, &
245&         cmatel2)
246        if (ipaw.ne.0) then
247!*** need to calculate pwmatel for all plane waves
248!*** ipw1 indexes full plane wave set, not the local field plane wave set
249!*** pwmatel is calculated for
250          call mkmatelP(pi,xred,natom,ibp,jband,ikpt,ikptq,qq,igg0,ngkpt, &
251&               pwmatel(:,:,ipw1,ix,iy,iz),tpwmatel(:,:,ipw1,ix,iy,iz), &
252&               projwf,nlmn,nkpt,nkptq,ncband,amatel2)
253        else
254          amatel2=0.d0
255        endif
256      else
257        call mkmatelX(ibp,jband,ikpt,ikptq,igpw,igg, &
258&         ncg,ncg,nkpt,nkpt,npwt,npwt,igmx,igmn,igndx,igndx, &
259&         isym,isymq,symrel,syminv,nsym,nsym,ihlf,ihlf,kpt,kpt, &
260&         lvtrans(1:3,ikk(1),ikk(2),ikk(3)),lvtrans(1:3,ikk(1),ikk(2),ikk(3)), &
261&         cg,cg,indxkcg,indxkcg,indxkpw,indxkpw,npwarr,npwarr,kg,kg, &
262&         cmatel2)
263        if (ipaw.ne.0) then
264          call mkmatelP(pi,xred,natom,ibp,jband,ikpt,ikptq,qq,igg0,ngkpt, &
265&               pwmatel(:,:,ipw1,ix,iy,iz),tpwmatel(:,:,ipw1,ix,iy,iz), &
266&               projwf,nlmn,nkpt,nkpt,ncband,amatel2)
267        else
268          amatel2=0.d0
269        endif
270      endif
271      vqmat2(ix,iy,iz)=(cmatel+amatel)*conjg(cmatel2+amatel2)*vq(ix,iy,iz)
272      if (ix.eq.ngkpt(1)/2.and.iy.eq.ngkpt(2)/2.and.iz.eq.ngkpt(3)/2) then
273        cmgamma2=(cmatel+amatel)*conjg(cmatel+amatel)
274      endif
275!write(6,*) 'cmatel  = ',cmatel
276!write(6,*) 'vqmat2 = ',vqmat2(iqq(1),iqq(2),iqq(3))
277!write(56,*) '>>>>>',iqq,vqmat2(iqq(1),iqq(2),iqq(3))
278!      stvec(iqq(3))=vq(ix,iy,iz)
279!      stvec(iqq(3))=sqrt(qq2)
280!      stvec(iqq(3))=dble(vqmat2(iqq(1),iqq(2),iqq(3)))
281!      stvec2(iqq(3))=dimag(vqmat2(iqq(1),iqq(2),iqq(3)))
282!      stvec(iqq(3))=dble(cmatel)
283!      stvec2(iqq(3))=dimag(cmatel)
284!      stvec(iqq(3))=dble(cmatel*conjg(cmatel))
285    enddo
286!    write(76,'(10f8.2)') stvec(:)
287!    write(76,'(10es8.1)') stvec(:)
288!    write(76,'(10es8.1)') stvec2(:)
289!    write(76,*)
290    enddo
291!    write(76,*) iqq(1)+1,'--------------------------------------------------',ibp
292    enddo
293!stop
294!return
295! Cut material for testing energy dependence of self energy at bottom of file
296
297    lx=ibp.le.nbocc
298    lcorr=ipwa.ne.0
299    if (lx) xse1=(0.d0,0.d0)
300    do ie=nbcore+1,ncband
301      cse1(ie)=(0.d0,0.d0)
302      if (lcorr) cse1corr(ie)=(0.d0,0.d0)
303      zz1(ie)=(0.d0,0.d0)
304      eps=enrgy(indxkbnd(ikpt)+ie)
305      if (ibp.gt.nbocc) then
306        eps1(ie)=eps
307      else
308        eps1(ie)=-eps
309      endif
310    enddo
311    brd=dw
312
313
314    if (abs(vqmat2(ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2)/vqmat2(ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2-1)).gt.1.d3 &
315&   .and.ipw1.eq.1) then
316!    if (.true.) then
317      iqsing=1
318      iqq=(/ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2/)
319      iqpt=iqndx(iqq(1),iqq(2),iqq(3))
320!      call fespread(iqq,ngkpt,omega,esprd)
321!      brd=max(esprd,dw)
322      qq2=4.d0*pi/vq(ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2)
323      wwq=sqrt(omegap2+(xkf**2/3.d0)*qq2+qq2**2/4.d0)
324      ppfct=pi*omegap2/(2.d0*wwq)
325      if (lcorr) then
326        call locateelement(iqq,ipwa,ipwa,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpw)
327      else
328        jjpw=0
329      endif
330      do ie=nbcore+1,ncband
331        eval=omega(iqq(1),iqq(2),iqq(3))+eps1(ie)-wwq
332        cenrgy=eval*(1.d0,0.d0)+brd*(0.d0,1.d0)
333        wppint0(ie)=dble(ppfct/cenrgy)
334        wpp2int0(ie)=-dble(ppfct/cenrgy**2)
335        wint0(ie)=(0.d0,0.d0)
336        w2int0(ie)=(0.d0,0.d0)
337        if (jjpw.ne.0) then
338          do iw=1,nwpt
339            ww=iw*dw
340            eval=omega(iqq(1),iqq(2),iqq(3))+eps1(ie)-ww
341            cenrgy=eval*(1.d0,0.d0)+brd*(0.d0,1.d0)
342            wint0(ie)=wint0(ie)+dw*dimag(lossfn(iw,iqpt,jjpw))/cenrgy
343            w2int0(ie)=w2int0(ie)-dw*dimag(lossfn(iw,iqpt,jjpw))/cenrgy**2
344          enddo
345        endif
346      enddo
347    else
348      iqsing=0
349    endif
350    do ix=1,ngkpt(1)
351    do iy=1,ngkpt(2)
352    do iz=1,ngkpt(3)
353!    do ix=1,1
354!    do iy=1,1
355!    do iz=1,1
356      iqq=(/ix,iy,iz/)
357      iqpt=iqndx(ix,iy,iz)
358!      call fespread(iqq,ngkpt,omega,esprd)
359!      brd=3*max(esprd,dw)
360!write(6,'(3i3,es20.10)') ix,iy,iz,esprd*27.2114
361      qq2=4.d0*pi/vq(ix,iy,iz)
362      wwq=sqrt(omegap2+(xkf**2/3.d0)*qq2+qq2**2/4.d0)
363      ppfct=pi*omegap2/(2.d0*wwq)
364      if (lcorr) then
365        call locateelement(iqq,ipwa,ipwa,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpw)
366      else
367        jjpw=0
368      endif
369      do ie=nbcore+1,ncband
370        eval=omega(ix,iy,iz)+eps1(ie)-wwq
371        cenrgy=eval*(1.d0,0.d0)+brd*(0.d0,1.d0)
372        wppint(ie)=ppfct/cenrgy
373        wpp2int(ie)=-ppfct/cenrgy**2
374        wint(ie)=(0.d0,0.d0)
375        w2int(ie)=(0.d0,0.d0)
376        if (jjpw.ne.0) then
377          do iw=1,nwpt
378            ww=iw*dw
379            eval=omega(ix,iy,iz)+eps1(ie)-ww
380            cenrgy=eval*(1.d0,0.d0)+brd*(0.d0,1.d0)
381            wint(ie)=wint(ie)+dw*dimag(lossfn(iw,iqpt,jjpw))/cenrgy
382            w2int(ie)=w2int(ie)-dw*dimag(lossfn(iw,iqpt,jjpw))/cenrgy**2
383          enddo
384        endif
385      enddo
386      if (iqsing.eq.0) then
387        if (lx) xse1=xse1+vqmat2(ix,iy,iz)
388        cse1pp=cse1pp+wppint*vqmat2(ix,iy,iz)
389        zz1pp=zz1pp+wpp2int*vqmat2(ix,iy,iz)
390        if (lcorr) then
391          cse1=cse1+wint*vqmat2(ix,iy,iz)
392          zz1=zz1+w2int*vqmat2(ix,iy,iz)
393          cse1corr=cse1corr &
394&                +wppint*vqmat2(ix,iy,iz)*(1.d0-fsum(ipwa,iqpt))
395          zz1corr=zz1corr &
396&                +wpp2int*vqmat2(ix,iy,iz)*(1.d0-fsum(ipwa,iqpt))
397        endif
398      else
399! Method one of handling singularity
400        if (sqrt(qq2).lt.gfo) then
401          if (ix.eq.ngkpt(1)/2.and.iy.eq.ngkpt(2)/2.and.iz.eq.ngkpt(3)/2) cycle
402          gterm=cmgamma2*(1+cos(pi*sqrt(qq2)/gfo))/(2.d0)
403          if (lx) xse1=xse1+vqmat2(ix,iy,iz)-gterm*vq(ix,iy,iz)
404          cse1pp=cse1pp+wppint*vqmat2(ix,iy,iz)-wppint0*gterm*vq(ix,iy,iz)
405          zz1pp=zz1pp+wpp2int*vqmat2(ix,iy,iz)-wpp2int0*gterm*vq(ix,iy,iz)
406          if (lcorr) then
407            cse1=cse1+wint*vqmat2(ix,iy,iz)-wint0*gterm*vq(ix,iy,iz)
408            zz1=zz1+w2int*vqmat2(ix,iy,iz)-w2int0*gterm*vq(ix,iy,iz)
409            cse1corr=cse1corr &
410&                 +wppint*vqmat2(ix,iy,iz)*(1.d0-fsum(ipwa,iqpt)) &
411&                 -wppint0*gterm*vq(ix,iy,iz)*(1.d0-fsum(ipwa,1))
412            zz1corr=zz1corr &
413&                 +wpp2int*vqmat2(ix,iy,iz)*(1.d0-fsum(ipwa,iqpt)) &
414&                 -wpp2int0*gterm*vq(ix,iy,iz)*(1.d0-fsum(ipwa,1))
415          endif
416! Method two of handling singularity
417!        if (ix.eq.ngkpt(1)/2.and.iy.eq.ngkpt(2)/2.and.iz.eq.ngkpt(3)/2) then
418!          cycle
419! End singularity selection
420        else
421          if (lx) xse1=xse1+vqmat2(ix,iy,iz)
422          cse1pp=cse1pp+wppint*vqmat2(ix,iy,iz)
423          zz1pp=zz1pp+wpp2int*vqmat2(ix,iy,iz)
424          if (lcorr) then
425            cse1=cse1+wint*vqmat2(ix,iy,iz)
426            zz1=zz1+w2int*vqmat2(ix,iy,iz)
427            cse1corr=cse1corr &
428&                 +wppint*vqmat2(ix,iy,iz)*(1.d0-fsum(ipwa,iqpt))
429            zz1corr=zz1corr &
430&                 +wpp2int*vqmat2(ix,iy,iz)*(1.d0-fsum(ipwa,iqpt))
431          endif
432        endif
433      endif
434    enddo
435    enddo
436    enddo
437    if (lx) xse1=xse1/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol)
438    cse1pp=cse1pp/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol*pi)
439    zz1pp=zz1pp/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol*pi)
440    if (lcorr) then
441      cse1=cse1/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol*pi)
442      zz1=zz1/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol*pi)
443      cse1corr=cse1corr/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol*pi)
444      zz1pp=zz1pp/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol*pi)
445    endif
446    if (iqsing.eq.1) then
447! Method one of handling singularity
448! 1/pi Integral (d^3q/(2 pi)^3) wint0 gterm vq
449      if (lx) xse1=xse1+cmgamma2*8.d0*pi*pi*gfo/((2.d0*pi)**3)
450      cse1pp=cse1pp+wppint0*cmgamma2*8.d0*pi*pi*gfo/(pi*(2.d0*pi)**3)
451      zz1pp=zz1pp+wpp2int0*cmgamma2*8.d0*pi*pi*gfo/(pi*(2.d0*pi)**3)
452      if (lcorr) then
453        cse1=cse1+wint0*cmgamma2*8.d0*pi*pi*gfo/(pi*(2.d0*pi)**3)
454        zz1=zz1+w2int0*cmgamma2*8.d0*pi*pi*gfo/(pi*(2.d0*pi)**3)
455        cse1corr=cse1corr+ &
456&         wppint0*cmgamma2*8.d0*pi*pi*gfo/(pi*(2.d0*pi)**3)*(1.d0-fsum(ipwa,1))
457        zz1corr=zz1corr+ &
458&         wpp2int0*cmgamma2*8.d0*pi*pi*gfo/(pi*(2.d0*pi)**3)*(1.d0-fsum(ipwa,1))
459      endif
460! Method two of handling singularity
461!      xdum=7.44*(((2.d0*pi)**3/vol)/ngkpt(1)*ngkpt(2)*ngkpt(3))**(-2.d0/3.d0) &
462!&           /(vol*ngkpt(1)*ngkpt(2)*ngkpt(3))
463!      if (lx) xse1=xse1+cmgamma2*xdum
464!      cse1pp=cse1pp+wppint0*cmgamma2*xdum/pi
465!      zz1pp=zz1pp+wpp2int0*cmgamma2*xdum/pi
466!      if (lcorr) then
467!        cse1=cse1+wint0*cmgamma2*xdum/pi
468!        zz1=zz1+w2int0*cmgamma2*xdum/pi
469!        cse1corr=cse1corr+ &
470!&         wppint0*cmgamma2*xdum*(1.d0-fsum(ipwa,1))/pi
471!        zz1corr=zz1corr+ &
472!&         wpp2int0*cmgamma2*xdum*(1.d0-fsum(ipwa,1))/pi
473!      endif
474! End singularity selection
475    endif
476    if (lx) xse1=-xse1
477    cse1pp=-cse1pp*isign
478    if (lcorr) then
479      cse1=-cse1*isign
480      zz1=-zz1*isign
481      cse1corr=-cse1corr*isign
482      zz1corr=-zz1corr*isign
483    endif
484    if (lx) xse=xse+dble(xse1)
485    csepp=csepp+cse1pp
486    if (lcorr) then
487      cse=cse+cse1
488      zz=zz+zz1
489      csecorr=csecorr+cse1corr
490      zzcorr=zzcorr+zz1corr
491    endif
492!if (lx.and.lcorr) then
493!  write(6,'(i4,4x,i5,3x,3i3,5x,"(",f10.6,",",f10.6,")",5x,f10.6)') &
494!&        ibp,ipw1,kg(:,indxkpw(ikpt)+ipw1),cse1(iband)*27.2114,dble(xse1)*27.2114
495!elseif (lcorr) then
496!  write(6,'(i4,4x,i5,3x,3i3,5x,"(",f10.6,",",f10.6,")")') &
497!&        ibp,ipw1,kg(:,indxkpw(ikpt)+ipw1),cse1(iband)*27.2114
498!endif
499  enddo
500!  if (lx) then
501!    write(6,'(i4,5x,"(",f10.6,",",f10.6,")",5x,f10.6)') ibp,(cse(iband)-cse2(iband))*27.2114,dble(xse-xse2)*27.2114
502!  else
503!    write(6,'(i4,5x,"(",f10.6,",",f10.6,")",f10.6)') ibp,(cse(iband)-cse2(iband))*27.2114
504!  endif
505enddo
506
507!zz=1.d0/(1.d0-zz)
508!zzpp=1.d0/(1.d0-zzpp)
509!zzcorr=1.d0/(1.d0-zz-zzcorr)
510!write(6,*) cse*27.2114
511!write(6,*) csecorr*27.2114
512!write(6,*) xse*27.2114
513!write(6,*) zz
514
515
516return
517end subroutine mkppse
518
519!**************************************************************************
520
521subroutine tebound(iqq,ngkpt,omega,insd)
522implicit none
523integer :: iqq(3),ngkpt(3),insd
524double precision :: omega(ngkpt(1),ngkpt(2),ngkpt(3))
525integer :: ix,iy,iz,iqp(3)
526double precision :: e1,e0
527
528insd=0
529e0=omega(iqq(1),iqq(2),iqq(3))
530do ix=0,1
531do iy=0,1
532do iz=0,1
533  iqp=mod(iqq+(/ix,iy,iz/)-1+ngkpt,ngkpt)+1
534  e1=omega(iqp(1),iqp(2),iqp(3))
535  if (e0*e1.lt.0.d0) insd=1
536enddo
537enddo
538enddo
539
540return
541end subroutine tebound
542
543!**************************************************************************
544! cut material from mkcse for testing energy dependence of self energy
545
546!    do ix=1,ngkpt(1)
547!    do iy=1,ngkpt(2)
548!    do iz=1,ngkpt(3)
549!      iqq=(/ix,iy,iz/)
550!      call fhilo(omega,iqq,ngkpt,whi,wlo)
551!      if (ix.eq.1.and.iy.eq.1.and.iz.eq.1) then
552!        wwhi=whi
553!        wwlo=wlo
554!      else
555!        wwhi=max(wwhi,whi)
556!        wwlo=min(wwlo,wlo)
557!      endif
558!    enddo
559!    enddo
560!    enddo
561!    if (ibp.gt.nbocc) then
562!      ioff=nint(-wwlo*dble(nwpt)/wmax)
563!    else
564!      ioff=nint(wwhi*dble(nwpt)/wmax)-nwpt
565!    endif
566!    ioff2=nint(ek*dble(nwpt)/wmax)
567!    do ie=1,nwpt
568!      eps=(ie+ioff-ioff2)*dw+ek+dw/2.d0
569!      ssi(ie)=(0.d0,0.d0)
570!!IF (.TRUE.) THEN
571!IF (.FALSE.) THEN
572!      if (ibp.gt.nbocc) then
573!        iwh=nint((eps+wwhi)*dble(nwpt)/wmax)
574!        iwl=nint((eps+wwlo)*dble(nwpt)/wmax)
575!      else
576!        iwh=nint((wwhi-eps)*dble(nwpt)/wmax)
577!        iwl=nint((wwlo-eps)*dble(nwpt)/wmax)
578!      endif
579!      do iw=max(iwl,1),min(iwh,nwpt)
580!        if (ibp.gt.nbocc) then
581!          ww=iw*dw-eps
582!        else
583!          ww=iw*dw+eps
584!        endif
585!!        write(6,*) ie,iw,ww*27.2114
586!        do ix=1,ngkpt(1)
587!        do iy=1,ngkpt(2)
588!        do iz=1,ngkpt(3)
589!          iqq=(/ix,iy,iz/)
590!          iqpt=iqndx(iqq(1),iqq(2),iqq(3))
591!          call locateelement(iqq,ipw1,ipw2,ngkpt,iqsymndx,npwc,nsym,pwsymndx,invpw2ndx,jjpw)
592!          if (jjpw.ne.0) then
593!            vfactor(iqq(1),iqq(2),iqq(3))=vqmat2(iqq(1),iqq(2),iqq(3))*dimag(lossfn(iw,iqpt,jjpw))
594!          else
595!            vfactor(iqq(1),iqq(2),iqq(3))=(0.d0,0.d0)
596!          endif
597!!          stvec(iqq(3))=dble(vfactor(iqq(1),iqq(2),iqq(3)))
598!!          stvec2(iqq(3))=dimag(vfactor(iqq(1),iqq(2),iqq(3)))
599!!          stvec(iqq(3))=dimag(lossfn(iw,iqpt,jjpw))
600!!          stvec(iqq(3))=dble(iqpt)
601!!          stvec2(iqq(3))=dble(jjpw)
602!        enddo
603!!        write(77,'(10es8.1)') stvec(:)
604!!        write(77,'(10es8.1)') stvec2(:)
605!!        write(77,'(10i8)') nint(stvec(:))
606!!        write(77,'(10i8)') nint(stvec2(:))
607!!        write(77,*)
608!        enddo
609!!        write(77,*) iqq(1)+1,'--------------------------------------------------',ibp
610!        enddo
611!        seiold=sei
612!        do ix=1,ngkpt(1)
613!        do iy=1,ngkpt(2)
614!        do iz=1,ngkpt(3)
615!          iqq=(/ix,iy,iz/)
616!          call fval(omega,iqq,iqqp,ngkpt,enval)
617!          call fpol(vfactor,ngkpt,iqq,iqqp,vfv)
618!!          write(42,*) '>>>>',iqq,ww
619!!          write(42,*)
620!!          write(42,'(2es12.3,2(4x,2es12.3))') enval(1:2),dble(vfv(1:2)),dimag(vfv(1:2))
621!!          write(42,'(2es12.3,2(4x,2es12.3))') enval(3:4),dble(vfv(3:4)),dimag(vfv(3:4))
622!!          write(42,*)
623!!          write(42,'(2es12.3,2(4x,2es12.3))') enval(5:6),dble(vfv(5:6)),dimag(vfv(5:6))
624!!          write(42,'(2es12.3,2(4x,2es12.3))') enval(7:8),dble(vfv(7:8)),dimag(vfv(7:8))
625!!          write(42,*)
626!          call cubeint(enval,ww,vfv,tseincr)
627!!          write(42,*) tseincr
628!          call fcenter(iqq,ngkpt,icenter)
629!          if (abs(tseincr).ne.0.d0.and.icenter.ne.0 &
630!&   .and.(ipw1.eq.1.or.ipw2.eq.1)) then
631!            rlr=1.d-2
632!            abr=1.d-2
633!!            write(42,*) '>>>>',iqq,ww
634!!            write(6,*) tseincr
635!            call fvq2(iqq,shiftk,shiftkq,bmet,ngkpt,ipw1,ipw2,ipwlf,npwlf,vq2)
636!!            write(42,*)
637!!            write(42,'(2es12.3,5x,2es12.3)') vq2(1:2),dble(vfv(1:2))
638!!            write(42,'(2es12.3,5x,2es12.3)') vq2(3:4),dble(vfv(3:4))
639!!            write(42,*)
640!!            write(42,'(2es12.3,5x,2es12.3)') vq2(5:6),dble(vfv(5:6))
641!!            write(42,'(2es12.3,5x,2es12.3)') vq2(7:8),dble(vfv(7:8))
642!!            write(42,*)
643!            if (ipw1.ge.ipw2) then
644!              lossv=vfv*vq2/(4.d0*pi)
645!              call subint2(iqq,shiftk,shiftkq,lossv,enval,ww,iw,bmet,iqndx,ngkpt,nwpt,nqpt,ikndxq,vol,ipw1,ipw2,ipwlf,npwlf,pi,rlr,abr,tseincr)
646!            else
647!              lossv=-conjg(vfv)*vq2/(4.d0*pi)
648!              tseincrx=-conjg(tseincr)
649!              call subint2(iqq,shiftk,shiftkq,lossv,enval,ww,iw,bmet,iqndx,ngkpt,nwpt,nqpt,ikndxq,vol,ipw1,ipw2,ipwlf,npwlf,pi,rlr,abr,tseincrx)
650!              tseincr=-conjg(tseincrx)
651!            endif
652!!            write(6,*) tseincr
653!          endif
654!!          write(42,*) tseincr
655!          dse=tseincr/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol)
656!          ssi(ie)=ssi(ie)+dse*dw*isign/pi
657!!          write(66,'(3i4,5x,es10.3)') iqq,dse
658!!          write(66,'(3i4,5x,f20.6)') iqq,dble(dse)*27.2114
659!!          write(66,'(3i4,a)') iqq,'-----------------------------------------------'
660!!          write(66,'(es10.3,5x,2es10.3)') ww,dse
661!!          write(66,*)
662!!          write(66,'(2es10.3,5x,2es10.3)') dble(vfv(1)),dble(vfv(2)),enval(1:2)-ww
663!!          write(66,'(2es10.3,5x,2es10.3)') dble(vfv(3)),dble(vfv(4)),enval(3:4)-ww
664!!          write(66,*)
665!!          write(66,'(2es10.3,5x,2es10.3)') dble(vfv(5)),dble(vfv(6)),enval(5:6)-ww
666!!          write(66,'(2es10.3,5x,2es10.3)') dble(vfv(7)),dble(vfv(8)),enval(7:8)-ww
667!!          write(42,*) ssi(iw)
668!!          stvec(iqq(3))=dble(tseincr)
669!!          stvec2(iqq(3))=dimag(tseincr)
670!!          stvec(iqq(3))=dble(icenter)
671!        enddo
672!!        write(78,'(10es8.1)') stvec(:)
673!!        write(78,'(10es8.1)') stvec2(:)
674!!        write(78,*)
675!!        write(78,'(10i4)') nint(stvec(:))
676!        enddo
677!!        write(78,*) iqq(1)+1,'--------------------------------------------------',ibp
678!        enddo
679!      enddo
680!ELSE
681!      brd=dw
682!      if (ibp.gt.nbocc) then
683!        eps1=eps
684!      else
685!        eps1=-eps
686!      endif
687!      if (abs(vqmat2(ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2)/vqmat2(ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2-1)).gt.1.d3 &
688!&     .and.(ipw1.eq.1.or.ipw2.eq.1)) then
689!        iqsing=1
690!        iqpt=iqndx(ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2)
691!        call locateelement(iqq,ipw1,ipw2,ngkpt,iqsymndx,npwc,nsym,pwsymndx,invpw2ndx,jjpw)
692!        wint0=(0.d0,0.d0)
693!        if (jjpw.ne.0) then
694!          do iw=1,nwpt
695!            ww=iw*dw
696!            eval=omega(ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2)+eps1-ww
697!            wint0=wint0+dw*dimag(lossfn(iw,iqpt,jjpw))/(eval,brd)
698!          enddo
699!        endif
700!!        ssi(ie)=ssi(ie)+wint0*cmgamma2*gfo/pi
701!!        ssi(ie)=ssi(ie)+wint0*cmgamma2*4.d0*pi*gfo**3*(1.d0/6.d0-1.d0/(pi**2))
702!      else
703!        iqsing=0
704!      endif
705!      do ix=1,ngkpt(1)
706!      do iy=1,ngkpt(2)
707!      do iz=1,ngkpt(3)
708!        iqq=(/ix,iy,iz/)
709!        iqpt=iqndx(ix,iy,iz)
710!        wint=(0.d0,0.d0)
711!        call locateelement(iqq,ipw1,ipw2,ngkpt,iqsymndx,npwc,nsym,pwsymndx,invpw2ndx,jjpw)
712!        if (jjpw.eq.0) cycle
713!        do iw=1,nwpt
714!          ww=iw*dw
715!          eval=omega(ix,iy,iz)+eps1-ww
716!          wint=wint+dw*dimag(lossfn(iw,iqpt,jjpw))/(eval,brd)
717!        enddo
718!        if (iqsing.eq.0) then
719!          ssi(ie)=ssi(ie)+wint*vqmat2(ix,iy,iz)
720!        else
721!          qq2=4.d0*pi/vq(ix,iy,iz)
722!          if (sqrt(qq2).lt.gfo) then
723!            gterm=cmgamma2*(1+cos(pi*sqrt(qq2)/gfo))/(2.d0)
724!            ssi(ie)=ssi(ie)+(wint*vqmat2(ix,iy,iz)/vq(ix,iy,iz)-wint0*gterm)*vq(ix,iy,iz)
725!          else
726!            gterm=0.d0
727!            ssi(ie)=ssi(ie)+wint*vqmat2(ix,iy,iz)
728!          endif
729!        endif
730!      enddo
731!      enddo
732!      enddo
733!      ssi(ie)=ssi(ie)/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol*pi)
734!      if (iqsing.eq.1) then
735!! 1/pi Integral (d^3q/(2 pi)^3) wint0 gterm vq
736!        ssi(ie)=ssi(ie)+wint0*cmgamma2*8.d0*pi*pi*gfo/(pi*(2.d0*pi)**3)
737!      endif
738!ENDIF
739!    enddo
740!
741!    do ie=1,nwpt
742!!      eps1=ek+dble(ie-nwpt/2)*dw+dw/2.d0
743!!      ssc(ie)=(0.d0,0.d0)
744!!      do je=1,nwpt
745!!        if (ie.eq.(je+ioff-ioff2+nwpt/2)) then
746!!          ssc(ie)=ssc(ie)-(0.d0,pi)*ssi(je)
747!!        else
748!!          eps2=dble(je+ioff-ioff2)*dw+ek+dw/2.d0
749!!          ssc(ie)=ssc(ie)+ssi(je)*dw/(eps2-eps1)
750!!        endif
751!!write(6,'(2i4,2f10.3,es12.3,2(2x,2es11.3))') ie,je,eps1*27.2114,eps2*27.2114,1/(eps2-eps1),ssc(ie)*27.2114,dble(ssi(je)*27.2114)
752!!      enddo
753!!      write(69,'(i4,f10.3,2(3x,2es12.3))') ie,eps1*27.2114,ssc(ie)*27.2114,ssi(ie)*27.2114
754!      sec(ie)=sec(ie)+ssi(ie)
755!    enddo
756