1subroutine setetint(ebkq,qq2,vme2,W,omega,vsign,lqsing,wmax, &
2& ipwa,ipwb,ngkpt,iqndx,iqsymndx,npwup, &
3& nsym,pwsymndx,invpw2ndx,nwpt,nqpt,ntpwndx, &
4& cse1)
5! Tetrahedron integration for self energy
6use geometry
7implicit none
8integer :: ipwa,ipwb,ngkpt(3),npwup,nsym,nwpt,nqpt,ntpwndx
9double precision :: ebkq(ngkpt(1),ngkpt(2),ngkpt(3))
10double precision :: qq2(ngkpt(1),ngkpt(2),ngkpt(3))
11double complex :: vme2(ngkpt(1),ngkpt(2),ngkpt(3))
12double complex :: W(nwpt,nqpt,ntpwndx)
13double precision :: omega(ngkpt(1),ngkpt(2),ngkpt(3))
14double complex :: cse1(-nwpt:nwpt)
15double precision :: wmax
16integer :: vsign
17logical :: lqsing
18integer :: iqsymndx(ngkpt(1),ngkpt(2),ngkpt(3)),invpw2ndx(npwup,npwup)
19integer :: pwsymndx(npwup,2*nsym)
20integer :: iqndx(ngkpt(1),ngkpt(2),ngkpt(3))
21
22double precision :: pi,fourpi
23double precision :: dw,eps1,eps2,wwlo,wwhi,whi,wlo
24integer :: iwwlo,iwwhi,iwlo,iwhi,ielo,iehi
25double precision :: de21,de31,de32,de41,de42,de43
26double precision :: fbx(4),fb(4),cmx(3,4),cm(3,4),xkt(3,3)
27double precision :: thresh
28integer :: iqq(3),iqqp(3),iqpt,iqv(3),ictr(3)
29double complex :: vfactor(ngkpt(1),ngkpt(2),ngkpt(3),nwpt)
30double precision :: qkcvt(3,3)
31double precision :: enval(8),vqa(8),qqa(8)
32double complex :: vfv(8,nwpt)
33double precision :: evtx0(4),evtx(4),rrpyr(3,4),vpyr0(4,nwpt),kvtx(3,4), &
34& xk(3,4),rg(3,3),tvol,vdum(3),vpyr(4),xpyr0(4),xpyr(4),vqvtx0(4),vqvtx(4)
35double precision :: aa0(nwpt),av(3,nwpt)
36integer :: indxe(4)
37logical :: lqcentr
38double precision :: bgrad(3),xmult
39double precision :: sint1,sint1a,sint1b,svec(3),sveca(3),svecb(3)
40double complex, allocatable :: ssi(:)
41!double complex :: ssi(-20*nwpt:20*nwpt)
42double complex :: ssc(-nwpt:nwpt)
43double precision :: rr(3,8)
44integer :: ivndx(4,6)
45data rr(1:3,1) /0.0d0,0.0d0,0.0d0/
46data rr(1:3,2) /1.0d0,0.0d0,0.0d0/
47data rr(1:3,3) /0.0d0,1.0d0,0.0d0/
48data rr(1:3,4) /1.0d0,1.0d0,0.0d0/
49data rr(1:3,5) /0.0d0,0.0d0,1.0d0/
50data rr(1:3,6) /1.0d0,0.0d0,1.0d0/
51data rr(1:3,7) /0.0d0,1.0d0,1.0d0/
52data rr(1:3,8) /1.0d0,1.0d0,1.0d0/
53data ivndx(1:4,1) /1,2,3,5/
54data ivndx(1:4,2) /3,5,6,7/
55data ivndx(1:4,3) /2,3,5,6/
56data ivndx(1:4,4) /2,3,4,6/
57data ivndx(1:4,5) /3,4,6,7/
58data ivndx(1:4,6) /4,6,7,8/
59parameter (pi=3.1415926535897932384626433832795, &
60& fourpi=12.566370614359172953850573533118d0)
61integer :: ii,jj,kk,iv,iw,ix,iy,iz,itet,iww,ie,je
62double precision :: www
63
64  dw=wmax/dble(nwpt)
65  do ii=1,3
66  do jj=1,3
67    qkcvt(ii,jj)=blat(ii,jj)/dble(ngkpt(ii))
68  enddo
69  enddo
70  ictr=ngkpt/2
71
72  do ie=-nwpt,nwpt
73    ssc(ie)=(0.d0,0.d0)
74  enddo
75  do ix=1,ngkpt(1)
76  do iy=1,ngkpt(2)
77  do iz=1,ngkpt(3)
78    iqq=(/ix,iy,iz/)
79    call fhilo(omega,iqq,ngkpt,whi,wlo)
80    if (ix.eq.1.and.iy.eq.1.and.iz.eq.1) then
81      wwhi=whi
82      wwlo=wlo
83    else
84      wwhi=max(wwhi,whi)
85      wwlo=min(wwlo,wlo)
86    endif
87  enddo
88  enddo
89  enddo
90
91  iwlo=nint(wwlo*dble(nwpt)/wmax)
92  iwhi=nint(wwhi*dble(nwpt)/wmax)
93  if (vsign.eq.1) then
94    ielo=iwlo+1
95    iehi=iwhi+nwpt
96  else
97    ielo=iwlo-nwpt
98    iehi=iwhi-1
99  endif
100  allocate (ssi(ielo:iehi))
101  do ie=ielo,iehi
102    ssi(ie)=0.d0
103  enddo
104  call mkvfactor(ipwa,ipwb,ngkpt,iqndx,iqsymndx,npwup,nsym,pwsymndx, &
105& invpw2ndx,nwpt,nqpt,ntpwndx,vme2,W,vfactor)
106
107  do ix=1,ngkpt(1)
108  do iy=1,ngkpt(2)
109  do iz=1,ngkpt(3)
110    iqq=(/ix,iy,iz/)
111    call fval(omega,iqq,iqqp,ngkpt,enval)
112    if (lqsing) then
113      call fval(qq2,iqq,iqqp,ngkpt,qqa)
114      vqa=fourpi/qqa
115    endif
116    do iw=1,nwpt
117      call fpol(vfactor(:,:,:,iw),ngkpt,iqq,iqqp,vfv(:,iw))
118    enddo
119    do itet=1,6
120      lqcentr=.false.
121      do iv=1,4
122        evtx0(iv)=enval(ivndx(iv,itet))
123        rrpyr(1:3,iv)=rr(1:3,ivndx(iv,itet))
124        do iw=1,nwpt
125          vpyr0(iv,iw)=vfv(ivndx(iv,itet),iw)
126        enddo
127        if (lqsing.and..not.lqcentr) then
128          iqv=iqq+nint(rrpyr(:,iv))
129          if (iqv(1).eq.ictr(1).and.iqv(2).eq.ictr(2).and. &
130& iqv(3).eq.ictr(3)) lqcentr=.true.
131        endif
132      enddo
133      call indxhpsort(4,4,evtx0,indxe)
134      evtx=evtx0(indxe)
135      do ii=1,4
136        jj=indxe(ii)
137        do kk=1,3
138          kvtx(kk,ii)=dot_product(iqq+rrpyr(:,jj)-ictr,qkcvt(:,kk))
139        enddo
140      enddo
141      xk(1:3,1)=(/0.d0,0.d0,0.d0/)
142      do ii=2,4
143        xk(1:3,ii)=kvtx(1:3,ii)-kvtx(1:3,1)
144      enddo
145      call cross(xk(:,3),xk(:,4),vdum)
146      tvol=dot_product(vdum,xk(:,2))
147      do jj=1,3  ! make contragradient
148        call cross(xk(1:3,mod(jj,3)+2),xk(1:3,mod(jj+1,3)+2),rg(1:3,jj))
149      enddo
150      rg=rg/tvol
151      tvol=abs(tvol)
152      iwwlo=nint(evtx(1)/dw)
153      iwwhi=nint(evtx(4)/dw)
154      de21=evtx(2)-evtx(1)
155      de31=evtx(3)-evtx(1)
156      de32=evtx(3)-evtx(2)
157      de41=evtx(4)-evtx(1)
158      de42=evtx(4)-evtx(2)
159      de43=evtx(4)-evtx(3)
160      thresh=1.d-5*de41
161      if (lqcentr) then            ! integral around coulomb singularity
162        do iv=1,4
163          vqvtx0(iv)=vqa(ivndx(iv,itet))
164        enddo
165        bgrad=(/0.d0,0.d0,0.d0/)   ! energy gradient
166        do jj=1,3
167          bgrad=bgrad+(evtx(jj+1)-evtx(1))*rg(:,jj)
168        enddo
169        xmult=4.d0*pi/sqrt(dot_product(bgrad,bgrad))
170        do iw=1,nwpt
171          vpyr(:)=vpyr0(:,iw)/vqvtx0
172          aa0(iw)=vpyr(indxe(1))
173          av(1:3,iw)=(/0.d0,0.d0,0.d0/)
174          do jj=1,3
175            av(1:3,iw)=av(1:3,iw) &
176&                     +(vpyr(indxe(jj+1))-vpyr(indxe(1)))*rg(1:3,jj)
177          enddo
178        enddo
179        do iww=iwwlo,iwwhi
180          www=dw*iww
181          if (www.lt.evtx(1)) then
182            cycle
183          elseif (www.lt.evtx(2)) then
184            call singint(1,www,xk,kvtx(:,1),evtx,sint1,svec)
185          elseif (www.lt.evtx(3)) then
186            if (de21.lt.thresh.and.de43.lt.thresh) then
187              call singint2(www,xk,kvtx(:,1),evtx,sint1a,sveca)
188            elseif (de21.ge.de43) then
189              call singint(2,www,xk,kvtx(:,1),evtx,sint1a,sveca)
190              call singint(1,www,xk,kvtx(:,1),evtx,sint1b,svecb)
191            else
192              call singint(3,www,xk,kvtx(:,1),evtx,sint1a,sveca)
193              call singint(4,www,xk,kvtx(:,1),evtx,sint1b,svecb)
194            endif
195            sint1=sint1b-sint1a
196            svec=svecb-sveca
197          elseif (www.lt.evtx(4)) then
198            call singint(4,www,xk,kvtx(:,1),evtx,sint1,svec)
199          else
200            cycle
201          endif
202          sint1=sint1*xmult
203          svec=svec*xmult
204          do iw=1,nwpt
205            ie=iww-vsign*iw
206            ssi(ie)=ssi(ie)+(aa0(iw)*sint1+dot_product(av(:,iw),svec))*dw
207          enddo
208        enddo
209      else
210        fbx(1)=tvol/(2.d0*de21*de31*de41)
211        if (de21.ge.de43) then
212          fbx(2)=tvol/(2.d0*de21*de32*de42)
213        else
214          fbx(3)=tvol/(2.d0*de31*de32*de43)
215        endif
216        fbx(4)=tvol/(2.d0*de41*de42*de43)
217        do ii=1,4
218          cmx(1:3,ii)=(/0.d0,0.d0,0.d0/)
219          do jj=1,4
220            if (jj.ne.ii) cmx(1:3,ii)=cmx(1:3,ii)+(xk(1:3,jj)-xk(1:3,ii)) &
221&                                                /(3*(evtx(jj)-evtx(ii)))
222          enddo
223        enddo
224        do iw=1,nwpt
225          aa0(iw)=vpyr0(indxe(1),iw)   ! base value of function
226          av(1:3,iw)=(/0.d0,0.d0,0.d0/) ! gradient of function
227          do jj=1,3
228            av(1:3,iw)=av(1:3,iw) &
229&                     +(vpyr0(indxe(jj+1),iw)-vpyr0(indxe(1),iw))*rg(1:3,jj)
230          enddo
231        enddo
232        do iww=iwwlo,iwwhi
233          www=dw*iww
234          if (www.lt.evtx(1)) then
235            cycle
236          elseif (www.lt.evtx(2)) then
237            sint1=fbx(1)*(www-evtx(1))**2
238            svec=(xk(1:3,1)+(www-evtx(1))*cmx(1:3,1))*sint1
239          elseif (www.lt.evtx(3)) then
240            if (de21.lt.thresh.and.de43.lt.thresh) then
241              xkt(:,1)=xk(:,3)*(www-evtx(1))/de31
242              xkt(:,2)=xk(:,4)*(www-evtx(1))/de41
243              xkt(:,3)=(xk(:,4)-xk(:,2))*(www-evtx(2))/de42+xk(:,2)
244              sint1=tvol*(2*(www-evtx(1))/(de31*de41) &
245&                   -(www-evtx(1))**2*(de31+de41)/(de31*de41)**2)/2
246              svec=((xkt(:,1)+xkt(:,3))/2.d0)*sint1
247            elseif (de21.ge.de43) then
248              fb(1)=fbx(1)*(www-evtx(1))**2
249              fb(2)=fbx(2)*(www-evtx(2))**2
250              sint1=fb(1)-fb(2)
251              cm(1:3,1)=xk(1:3,1)+(www-evtx(1))*cmx(1:3,1)
252              cm(1:3,2)=xk(1:3,2)+(www-evtx(2))*cmx(1:3,2)
253              svec=cm(:,1)*fb(1)-cm(:,2)*fb(2)
254            else
255              fb(3)=fbx(3)*(www-evtx(3))**2
256              fb(4)=fbx(4)*(www-evtx(4))**2
257              sint1=fb(4)-fb(3)
258              cm(1:3,3)=xk(1:3,3)+(www-evtx(3))*cmx(1:3,3)
259              cm(1:3,4)=xk(1:3,4)+(www-evtx(4))*cmx(1:3,4)
260              svec=cm(:,4)*fb(4)-cm(:,3)*fb(3)
261            endif
262          elseif (www.lt.evtx(4)) then
263            sint1=fbx(4)*(www-evtx(4))**2
264            svec=(xk(1:3,4)+(www-evtx(4))*cmx(1:3,4))*sint1
265          else
266            cycle
267          endif
268          do iw=1,nwpt
269            ie=iww-vsign*iw
270            ssi(ie)=ssi(ie)+(aa0(iw)*sint1+dot_product(av(:,iw),svec))*dw
271          enddo
272        enddo
273      endif
274    enddo
275  enddo
276  enddo
277  enddo
278  ssi=-vsign*ssi/((2*pi)**3*pi)
279
280  do ie=-nwpt,nwpt
281    eps1=dble(ie)*dw-dw/2
282    do je=ielo,iehi
283      if (ie.eq.je) then
284        ssc(ie)=ssc(ie)+(0.d0,pi)*ssi(je)
285      else
286        eps2=dble(je)*dw-dw/2
287        ssc(ie)=ssc(ie)+vsign*ssi(je)*dw/(eps2-eps1)
288      endif
289    enddo
290  enddo
291  deallocate (ssi)
292
293  cse1=cse1+ssc
294
295!  cse=(cse1(0)+cse1(1))/2
296!  dcse=(cse1(1)-cse1(0))/dw
297
298return
299end subroutine setetint
300
301!*****************************************************************************
302
303subroutine mkvfactor(ipwa,ipwb,ngkpt,iqndx,iqsymndx,npwup,nsym,pwsymndx,invpw2ndx,nwpt,nqpt,ntpwndx,vme2,W,vfactor)
304implicit none
305integer :: ipwa,ipwb,ngkpt(3),npwup,nsym,nwpt,nqpt,ntpwndx
306integer :: iqsymndx(ngkpt(1),ngkpt(2),ngkpt(3)),invpw2ndx(npwup,npwup)
307integer :: pwsymndx(npwup,2*nsym)
308integer :: iqndx(ngkpt(1),ngkpt(2),ngkpt(3))
309double complex :: vme2(ngkpt(1),ngkpt(2),ngkpt(3))
310double complex :: W(nwpt,nqpt,ntpwndx)
311double complex :: vfactor(ngkpt(1),ngkpt(2),ngkpt(3),nwpt)
312integer iw,ix,iy,iz,iqq(3),iqpt,jjpw
313
314  do iw=1,nwpt
315    do ix=1,ngkpt(1)
316    do iy=1,ngkpt(2)
317    do iz=1,ngkpt(3)
318      iqq=(/ix,iy,iz/)
319      iqpt=iqndx(iqq(1),iqq(2),iqq(3))
320      call locateelement(iqq,ipwa,ipwb,ngkpt,iqsymndx,npwup,nsym,pwsymndx,invpw2ndx,jjpw)
321      if (jjpw.ne.0) then
322        vfactor(ix,iy,iz,iw)=vme2(ix,iy,iz)*dimag(W(iw,iqpt,jjpw))
323      else
324        vfactor(ix,iy,iz,iw)=(0.d0,0.d0)
325      endif
326    enddo
327    enddo
328    enddo
329  enddo
330
331return
332end subroutine mkvfactor
333
334!*****************************************************************************
335
336subroutine seqptsum(ikpt,lqsing,lx,lc,ipwa,ipwb,ngkpt,iqndx,omegap2,xkf,qq2,W, &
337& npwup,nsym,iqsymndx,pwsymndx,nbcore,nbocc,ncband,nwpt,nqpt,ntpwndx, &
338& wmax,vme2,ame2,falloff,vsign,ebkq,indxkbnd,nkpt,eigen,bantot,invpw2ndx, &
339& exch,ppcor,cor)
340use geometry
341implicit none
342logical :: lqsing,lx,lc
343integer :: ikpt,ipwa,ipwb,nbcore,nbocc,ncband,nwpt,nqpt,ntpwndx,vsign
344integer :: nsym,npwup,nkpt,bantot
345integer :: ngkpt(3)
346integer :: iqndx(ngkpt(1),ngkpt(2),ngkpt(3))
347integer :: indxkbnd(nkpt)
348double precision :: qq2(ngkpt(1),ngkpt(2),ngkpt(3))
349double precision :: omegap2,xkf,wmax,falloff
350integer :: iqsymndx(ngkpt(1),ngkpt(2),ngkpt(3)),invpw2ndx(npwup,npwup)
351integer :: pwsymndx(npwup,2*nsym)
352double complex :: W(nwpt,nqpt,ntpwndx)
353double complex :: vme2(ngkpt(1),ngkpt(2),ngkpt(3))
354double complex :: ame2(ngkpt(1),ngkpt(2),ngkpt(3))
355double precision :: exch
356double complex :: ppcor(nbcore+1:ncband),cor(nbcore+1:ncband)
357double complex :: cse(-nwpt:nwpt)
358double precision :: ebkq(ngkpt(1),ngkpt(2),ngkpt(3))
359double precision :: eigen(bantot)
360
361double precision :: exch1
362double complex :: ppcor1(nbcore+1:ncband),cor1(nbcore+1:ncband)
363double complex :: cse1(-nwpt:nwpt)
364
365double precision :: pi,fourpi
366double precision :: ww,wwq,ppfct,omegat,brd,dw,volelmnt
367double complex :: cedenom,cterm
368double complex :: wppint0(nbcore+1:ncband),wpp2int0(nbcore+1:ncband)
369double complex :: wint0(nbcore+1:ncband),w2int0(nbcore+1:ncband)
370double complex :: wppint(nbcore+1:ncband),wpp2int(nbcore+1:ncband)
371double complex :: wint(nbcore+1:ncband),w2int(nbcore+1:ncband)
372integer :: ictr(3),iqpt,jjpw,iqq(3)
373integer :: ie,iw,ix,iy,iz,ii,jj,kk
374double precision :: singterm
375parameter (pi=3.1415926535897932384626433832795, &
376& fourpi=12.566370614359172953850573533118d0)
377
378  volelmnt=vol*ngkpt(1)*ngkpt(2)*ngkpt(3)
379  dw=wmax/dble(nwpt)
380  brd=dw
381  ictr=ngkpt/2
382  if (lqsing) then
383    iqpt=iqndx(ictr(1),ictr(2),ictr(3))
384    wwq=sqrt(omegap2+(xkf**2/3.d0)*qq2(ictr(1),ictr(2),ictr(3))  &
385&                 +qq2(ictr(1),ictr(2),ictr(3))**2/4.d0)
386    ppfct=pi*omegap2/(2.d0*wwq)
387    if (lc) then
388      call locateelement(ictr,ipwa,ipwb,ngkpt,iqsymndx,npwup,nsym,pwsymndx,invpw2ndx,jjpw)
389    else
390      jjpw=0
391    endif
392    do ie=nbcore+1,ncband
393      omegat=eigen(indxkbnd(ikpt)+ie)-ebkq(ictr(1),ictr(2),ictr(3))
394      cedenom=(omegat-vsign*wwq)*(1.d0,0.d0)+vsign*brd*(0.d0,1.d0)
395      wppint0(ie)=(ppfct/cedenom)/pi
396      wpp2int0(ie)=-(ppfct/cedenom**2)/pi
397      wint0(ie)=(0.d0,0.d0)
398      w2int0(ie)=(0.d0,0.d0)
399      if (jjpw.ne.0) then
400        do iw=1,nwpt
401          ww=iw*dw
402          cedenom=(omegat-vsign*ww)*(1.d0,0.d0)+vsign*brd*(0.d0,1.d0)
403          wint0(ie)=wint0(ie)+dw*(dimag(W(iw,iqpt,jjpw))/cedenom)/pi
404          w2int0(ie)=w2int0(ie)-dw*(dimag(W(iw,iqpt,jjpw))/cedenom**2)/pi
405        enddo
406      endif
407    enddo
408  endif
409
410  exch1=0.d0
411  do ie=nbcore+1,ncband
412    ppcor1(ie)=(0.d0,0.d0)
413    cor1(ie)=(0.d0,0.d0)
414  enddo
415
416  do ix=1,ngkpt(1)
417  do iy=1,ngkpt(2)
418  do iz=1,ngkpt(3)
419    iqq=(/ix,iy,iz/)
420    iqpt=iqndx(iqq(1),iqq(2),iqq(3))
421    wwq=sqrt(omegap2+(xkf**2/3.d0)*qq2(ix,iy,iz)  &
422&           +qq2(ix,iy,iz)**2/4.d0)
423    ppfct=pi*omegap2/(2.d0*wwq)
424    if (lc) then
425      call locateelement(iqq,ipwa,ipwb,ngkpt,iqsymndx,npwup,nsym,pwsymndx,invpw2ndx,jjpw)
426    else
427      jjpw=0
428    endif
429    do ie=nbcore+1,ncband
430      omegat=eigen(indxkbnd(ikpt)+ie)-ebkq(ix,iy,iz)
431      cedenom=(omegat-vsign*wwq)*(1.d0,0.d0)+vsign*brd*(0.d0,1.d0)
432      wppint(ie)=(ppfct/cedenom)/pi
433      wpp2int(ie)=-(ppfct/cedenom**2)/pi
434      wint(ie)=(0.d0,0.d0)
435      w2int(ie)=(0.d0,0.d0)
436      if (jjpw.ne.0) then
437        do iw=1,nwpt
438          ww=iw*dw
439          cedenom=(omegat-vsign*ww)*(1.d0,0.d0)+vsign*brd*(0.d0,1.d0)
440          wint(ie)=wint(ie)+dw*(dimag(W(iw,iqpt,jjpw))/cedenom)/pi
441          w2int(ie)=w2int(ie)-dw*(dimag(W(iw,iqpt,jjpw))/cedenom**2)/pi
442        enddo
443      endif
444    enddo
445    if (.not.lqsing) then
446      if (lx) exch1=exch1-dble(vme2(ix,iy,iz))
447      ppcor1=ppcor1+vme2(ix,iy,iz)*wppint
448      if (lc) cor1=cor1+vme2(ix,iy,iz)*wint
449    else
450!*** method 1 for dealing with singularity ***
451!      if (qq2(ix,iy,iz).lt.falloff**2) then
452!        if (ix.eq.ictr(1).and.iy.eq.ictr(2).and.iz.eq.ictr(3)) cycle
453!        cterm=ame2(ictr(1),ictr(2),ictr(3))*(1+cos(pi*sqrt(qq2(ix,iy,iz))/falloff))/2.d0
454!        if (lx) exch1=exch1-dble(vme2(ix,iy,iz)-cterm*fourpi/qq2(ix,iy,iz))
455!        ppcor1=ppcor1+vme2(ix,iy,iz)*wppint-cterm*wppint0*fourpi/qq2(ix,iy,iz)
456!        if (lc) cor1=cor1+vme2(ix,iy,iz)*wint-cterm*wint0*fourpi/qq2(ix,iy,iz)
457!*** method 2 for dealing with singularity ***
458      if (ix.eq.ictr(1).and.iy.eq.ictr(2).and.iz.eq.ictr(3)) then
459        cycle
460!*** end alternate singularity method
461      else
462        if (lx) exch1=exch1-dble(vme2(ix,iy,iz))
463        ppcor1=ppcor1+vme2(ix,iy,iz)*wppint
464        if (lc) cor1=cor1+vme2(ix,iy,iz)*wint
465      endif
466    endif
467  enddo
468  enddo
469  enddo
470  if (lx) exch1=exch1/volelmnt
471  ppcor1=ppcor1/volelmnt
472  if (lc) cor1=cor1/volelmnt
473
474  if (lqsing) then
475!*** method 1 for dealing with singularity ***
476!    if (lx) exch1=exch1-ame2(ictr(1),ictr(2),ictr(3))*falloff/pi
477!    ppcor1=ppcor1+wppint0*ame2(ictr(1),ictr(2),ictr(3))*falloff/pi
478!    if (lc) cor1=cor1+wint0*ame2(ictr(1),ictr(2),ictr(3))*falloff/pi
479!*** method 2 for dealing with singularity ***
480    singterm=7.44d0*(vbz/(ngkpt(1)*ngkpt(2)*ngkpt(3)))**(-2.d0/3.d0)/volelmnt
481    if (lx) exch1=exch1-ame2(ictr(1),ictr(2),ictr(3))*singterm
482    ppcor1=ppcor1+wppint0*ame2(ictr(1),ictr(2),ictr(3))*singterm
483    if (lc) cor1=cor1+wint0*ame2(ictr(1),ictr(2),ictr(3))*singterm
484!*** end alternate singularity method
485  endif
486  exch=exch+exch1
487  ppcor=ppcor+ppcor1
488  cor=cor+cor1
489
490return
491end subroutine seqptsum
492