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