1subroutine mkse(iband,ikk,lossfn, & 2& vol,pi,nwpt,wmax,nbcore,nbocc,ncband,ngkpt,natom,xred,projwf,nlmn, & 3& test_bands_se, & 4& ntypepaw, & 5& pwmatel,tpwmatel, pwjmatel, tpwjmatel, & 6& kg,enrgy,cg,npwt,bantot,ncg, & 7& indxkpw,indxkbnd,indxkcg,npwarr, & 8& kpt,nkpt,nqpt,nsymk,symk,nsym,symrel,syminv, & 9& ihlf,lvtrans,bmet,blat,ipaw,itetrahedron, & 10& ipwx,ipwndx,npwndx,ntpwndx,napwndx, & 11& npwc,npwx,invpw2ndx,pwsymndx,iqsymndx, & 12& igmx,igmn,igndx,ikndx,iqndx,isymndx,isymg,npw, & 13& nband,nsppol,shiftk,zz,cse,xse) 14implicit none 15integer :: iband,ikk(3),nwpt,nbcore,nbocc,ncband,ngkpt(3),natom,nlmn 16integer :: igmn(3),igmx(3) 17integer :: npwt,bantot,ncg,nkpt,nqpt,nsym,nsppol 18integer :: npw,npwc,npwx,ipw1,npwndx,ntpwndx,napwndx,ipaw,itetrahedron 19double precision :: vol,pi,wmax,xred(3,natom) 20double complex :: lossfn(nwpt,nqpt+9,ntpwndx) 21integer :: test_bands_se(2) 22double complex :: projwf(natom,nlmn,nkpt,ncband) 23integer :: kg(3,npwt) 24double precision :: enrgy(bantot) 25double complex :: cg(ncg) 26integer :: indxkpw(nkpt),indxkbnd(nkpt) 27integer :: indxkcg(nkpt),npwarr(nkpt) 28double precision :: kpt(3,nkpt),shiftk(3) 29integer :: nsymk(nkpt),symk(nkpt,nsym*2) 30integer :: symrel(3,3,nsym),syminv(3,3,nsym) 31integer :: lvtrans(3,ngkpt(1),ngkpt(2),ngkpt(3)) 32integer :: ihlf(nkpt) 33double precision :: bmet(3,3),blat(3,3),bmetinv(3,3),bmet_t(3,3),bmet2(3,3) 34double precision :: volelmnt, volred 35integer :: ipwx(3,npwx),ipwndx(2,napwndx) 36integer :: igndx(nkpt,igmn(1):igmx(1),igmn(2):igmx(2),igmn(3):igmx(3)) 37integer :: ikndx(ngkpt(1),ngkpt(2),ngkpt(3)) 38integer :: iqndx(ngkpt(1),ngkpt(2),ngkpt(3)) 39integer :: isymndx(ngkpt(1),ngkpt(2),ngkpt(3)) 40integer :: isymg(3,nkpt,nsym),igsymk(3),igsymq(3) 41integer :: nband(nkpt*nsppol) 42double complex :: sei,seipw,seiold,seibc,csepw,cseold,csebc 43double complex, allocatable :: ssi(:) 44double complex :: ssc(-nwpt:nwpt),cse,dcse,cse0,dcse0,cse1(-nwpt:nwpt),zz,ctest,cse2(-nwpt:nwpt) 45double precision :: xse,xse0,ssx,ssx1,ssx2,xse2,ssx3 46integer :: ii,jj,kk,ll,ix,iy,iz,iskip,ocsign,iiq 47integer :: iqq(3),jka(3),jkb(3),jkk(3),iqv(3),ictr(3),iqr(3),ixr(3) 48double precision :: qr(3,4),qqr(3,4) 49integer :: ixx(3),ixxp(3),ixp2(3),ixp1(3),ixm1(3),ixm2(3) 50integer :: ikpt,ikptq 51integer :: igg(3),igh(3),igg0(3),isym,isymq,ibp,iqpt,iqsym,iw,iww,ie,je,ies,iqctr 52double precision :: xck(3),xckq(3),qadj(6) 53double complex :: cmatel,cmatel2,amatel,amatel2 54double complex :: jmatel(3),jmatel2(3),vmatel(3),vmatel2(3) 55double complex :: vqmat2(ngkpt(1),ngkpt(2),ngkpt(3)) 56double complex :: xmat2(ngkpt(1),ngkpt(2),ngkpt(3)) 57double complex :: jmat2(3,3), rjmat2(3,3), rjmatdum 58double precision :: smat2 59double complex :: vfactor(ngkpt(1),ngkpt(2),ngkpt(3),nwpt),vjfactor(3,3,nwpt),vsfactor(3,3,nwpt) 60double complex :: vfv(8,nwpt),vfx(8),vsm1(3,3,nwpt),vsm2(3,3,nwpt) 61double precision :: vq2(8),vqvtx0(4),vqvtx(4),qkcvt(3,3) 62double precision :: omega(ngkpt(1),ngkpt(2),ngkpt(3)) 63double precision :: whi,wlo,wwhi,wwlo,ww,www,wwwlo,wwwhi,enval(8),dw,eshift,evtx0(4),evtx(4) 64double precision :: rrpyr(3,4),kvtx(3,4),xk(3,4),ckvtx(3,4),cxk(3,4),rg(3,3),tvol,xpyr0(4),xpyr(4) 65double precision :: xkc(3,4),rgc(3,3),tvolc,kvtxc(3,4) 66double complex :: vpyr0(4,nwpt),vpyr(4) 67double precision :: de21,de31,de32,de41,de42,de43,thresh 68double precision :: fbx(4),fb(4),cmx(3,4),cm(3,4),xkt(3,3) 69double complex :: aa0(nwpt),av(3,nwpt) 70double precision :: xme,xv(3) 71double precision :: avec(3),bgrad(3),xmult 72double precision :: abr,rlr,rlr0 73integer :: iwh,iwl,ibmin,ibmax,iwhi,iwlo,iehi,ielo 74integer :: indxe(4),iwwlim(4) 75double precision :: vq(ngkpt(1),ngkpt(2),ngkpt(3)),qq(3),qp(3),qq2,qp2,qs(3),ek,ekmq 76double precision :: stvec(ngkpt(3)),stvec2(ngkpt(3)),temparray(ngkpt(1),ngkpt(2),ngkpt(3)) 77integer :: iipw,jjpw,jjpwt 78integer :: iqsymndx(ngkpt(1),ngkpt(2),ngkpt(3)),invpw2ndx(npwx,npwx) 79integer :: pwsymndx(npwc,2*nsym) 80integer :: ipw2,jpw1,jpw2,jw,iqcentr,jqcentr 81double precision :: eps1,eps2,eps 82double precision :: eval,brd,esprd 83double precision :: gfo,gamma(3),pln(3),dist(3) 84double precision :: sint1,sint1a,sint1b,svec(3),sveca(3),svecb(3),stens(3,3) 85double complex :: cint 86double complex :: wint,wint0,w2int,w2int0,gterm,vcentr 87double precision :: wcentr(-nwpt:nwpt),wgrid(nwpt) 88integer :: ntypepaw 89double complex :: pwmatel(ntypepaw,nlmn,nlmn,npwx,ngkpt(1),ngkpt(2),ngkpt(3)), & 90& tpwmatel(ntypepaw,nlmn,nlmn,npwx,ngkpt(1),ngkpt(2),ngkpt(3)), & 91& pwjmatel(3,ntypepaw,nlmn,nlmn), & 92& tpwjmatel(3,ntypepaw,nlmn,nlmn) 93logical :: lqsing,lqcentr 94logical :: lx,lc 95double precision :: linterp 96double precision :: rr(3,8) 97integer :: ivndx(4,6),itet,iv 98integer :: nsing 99double precision :: wsing(20) 100character*4 :: label 101external linterp 102data rr(1:3,1) /0.0d0,0.0d0,0.0d0/ 103data rr(1:3,2) /1.0d0,0.0d0,0.0d0/ 104data rr(1:3,3) /0.0d0,1.0d0,0.0d0/ 105data rr(1:3,4) /1.0d0,1.0d0,0.0d0/ 106data rr(1:3,5) /0.0d0,0.0d0,1.0d0/ 107data rr(1:3,6) /1.0d0,0.0d0,1.0d0/ 108data rr(1:3,7) /0.0d0,1.0d0,1.0d0/ 109data rr(1:3,8) /1.0d0,1.0d0,1.0d0/ 110data ivndx(1:4,1) /1,2,3,5/ 111data ivndx(1:4,2) /3,5,6,7/ 112data ivndx(1:4,3) /2,3,5,6/ 113data ivndx(1:4,4) /2,3,4,6/ 114data ivndx(1:4,5) /3,4,6,7/ 115data ivndx(1:4,6) /4,6,7,8/ 116double precision :: aa,bb,cc,dd,ee,ff 117integer :: iter 118data iter /30/ 119logical centerint 120integer vcenterint 121integer :: i_prt_DEBUG, j_prt_DEBUG, k_prt_DEBUG, & 122& i_ct_DEBUG, j_ct_DEBUG, k_ct_DEBUG, l_ct_DEBUG, & 123& k_ref_DEBUG, i_ref_DEBUG 124common /tsingint_DEBUG/ i_prt_DEBUG, j_prt_DEBUG, k_prt_DEBUG, & 125& i_ct_DEBUG, j_ct_DEBUG, k_ct_DEBUG, l_ct_DEBUG, & 126& k_ref_DEBUG, i_ref_DEBUG 127integer :: idum 128double precision :: xdum,xdum2,xdum3,vdum(3),mtxdum(3,3),mtxdum2(3,3) 129double complex :: cdum,cdum2,cdum3,cvecdum(3),cvecdum2(3),ctensdum(3,3),ctensdum2(3,3) 130double precision :: xunit(3,3) 131 132do ii=1,3 133do jj=1,3 134 if (ii.eq.jj) then 135 xunit(ii,jj)=1.d0 136 else 137 xunit(ii,jj)=0.d0 138 endif 139enddo 140enddo 141 142i_prt_DEBUG = 0 143j_prt_DEBUG = 0 144k_prt_DEBUG = 0 145i_ct_DEBUG = 0 146j_ct_DEBUG = 0 147k_ct_DEBUG = 0 148l_ct_DEBUG = 0 149i_ref_DEBUG = -1 150k_ref_DEBUG = -1 151abr=1.d-100 152rlr0=3.d-2 153rlr=1.d-3 154sei=0.d0 155ctest=(0.d0,0.d0) 156cse=(0.d0,0.d0) 157dcse=(0.d0,0.d0) 158zz=(0.d0,0.d0) 159igg0=(/0,0,0/) 160do ie=-nwpt,nwpt 161 cse1(ie)=(0.d0,0.d0) 162enddo 163xse=0.d0 164dw=wmax/dble(nwpt) 165do iw=1,nwpt 166 wgrid(iw)=iw*dw 167enddo 168volelmnt=(vol*ngkpt(1)*ngkpt(2)*ngkpt(3))/((2.d0*pi)**3) 169volred = 1.d0 170ikpt=ikndx(ikk(1),ikk(2),ikk(3)) 171isym=isymndx(ikk(1),ikk(2),ikk(3)) 172igsymk=isymg(:,ikpt,isym) 173ek=enrgy(indxkbnd(ikpt)+iband) 174do ii=1,3 175do jj=1,3 176 qkcvt(ii,jj)=blat(ii,jj)/dble(ngkpt(ii)) 177enddo 178enddo 179ictr=(/ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2/) 180iqctr=iqndx(ictr(1),ictr(2),ictr(3)) 181nsing=0 182 183do ii=1,3 184do kk=1,3 185 bmet2(ii,kk) = 0.d0 ! bmet2 = bmet*bmet 186 do jj=1,3 187 bmet2(ii,kk) = bmet2(ii,kk) + bmet(ii,jj)*bmet(jj,kk) 188 enddo 189enddo 190enddo 191 192do ix=1,ngkpt(1) 193do iy=1,ngkpt(2) 194do iz=1,ngkpt(3) 195 temparray(ix,iy,iz)=(0.d0,0.d0) 196enddo 197enddo 198enddo 199 200gamma=(blat(1,:)+blat(2,:)+blat(3,:))/2.d0 201do ii=1,3 202 jj=mod(ii,3)+1 203 kk=mod(ii+1,3)+1 204 pln(1)=blat(jj,2)*blat(kk,3)-blat(jj,3)*blat(kk,2) 205 pln(2)=-blat(jj,1)*blat(kk,3)+blat(jj,3)*blat(kk,1) 206 pln(3)=blat(jj,1)*blat(kk,2)-blat(jj,2)*blat(kk,1) 207 dist(ii)=dot_product(gamma,pln)/sqrt(dot_product(pln,pln)) 208enddo 209gfo=minval(dist) 210 211!ipaw=0 212!write(6,'(a)') ' finding contibutions from:' 213!write(6,'(a)') 'band plane waves correlation exchange' 214do ibp=test_bands_se(1),test_bands_se(2) 215!do ibp=1,1 216 write(6,'(3(a,i4))') 'quasiparticle band ',iband,', polarized band ',ibp,' out of ',test_bands_se(2) 217 seibc=sei 218 cse2=cse1 219 xse2=xse 220 if (ibp.gt.nbocc) then 221 ocsign=-1 222 else 223 ocsign=1 224 endif 225 do iipw=1,napwndx 226! do iipw=1,1 227! write(6,'(a,i4)') 'iipw ',iipw 228 seipw=sei 229 do ie=-nwpt,nwpt 230 ssc(ie)=(0.d0,0.d0) 231 enddo 232 ipw1=ipwndx(1,iipw) 233 ipw2=ipwndx(2,iipw) 234! if (ipw1.ne.ipw2) cycle 235! write(6,'(1x,i3,4x,2i4)') ibp,ipw1,ipw2 236 lx=ipw1.eq.ipw2.and.ibp.le.nbocc ! do exchange 237 lc=iipw.le.ntpwndx ! do correlation 238 if ((.not.lx).and.(.not.lc)) cycle 239 240! Step 1: Compute matrix elements 241 do ix=1,ngkpt(1) 242 do iy=1,ngkpt(2) 243 do iz=1,ngkpt(3) 244! do ix=1,1 245! do iy=10,10 246! do iz=10,10 247 vqmat2(ix,iy,iz)=(0.d0,0.d0) 248 ixx=(/ix,iy,iz/) 249 iqq=ixx-ictr 250 qq = dble(iqq)/dble(ngkpt) 251 jka=ikk-iqq 252 jkk=mod(jka,ngkpt) 253 do ii=1,3 254 if (jkk(ii).le.0) jkk(ii)=jkk(ii)+ngkpt(ii) 255 enddo 256 ikptq=ikndx(jkk(1),jkk(2),jkk(3)) 257 isymq=isymndx(jkk(1),jkk(2),jkk(3)) 258 igsymq=isymg(:,ikptq,isymq) 259 jkb=nint(kpt(:,ikptq)*dble(ngkpt))+ictr 260 igg=nint(dble(jkk-jka)/dble(ngkpt)) 261!write(6,*) '>>>>> iqq = ',iqq 262!write(6,*) '>>>>> ikk = ',ikk 263!write(6,*) '>>>>> jkk = ',jkk 264!write(6,*) '>>>>> jka = ',jka 265!write(6,*) '>>>>> igg = ',igg 266!write(6,'(a,3f10.5)') 'kpt = ',kpt(:,ikpt) 267!write(6,'(a,3f10.5)') 'kptq = ',kpt(:,ikptq) 268!read(*,*) 269 qp=qq+ipwx(:,ipw2) 270 qq=qq+ipwx(:,ipw1) 271 qq2=0.d0 272 qp2=0.d0 273 do ii=1,3 274 do jj=1,3 275 qq2=qq2+qq(ii)*bmet(ii,jj)*qq(jj) 276 qp2=qp2+qp(ii)*bmet(ii,jj)*qp(jj) 277 enddo 278 enddo 279 vq(ix,iy,iz)=4.d0*pi/sqrt(qq2*qp2) 280!write(6,'(a,3f10.5)') 'qq = ',qq 281!write(6,'(a,3f10.5)') 'qp = ',qp 282!write(6,'(a,3f10.5)') 'qq2 = ',qq2 283!write(6,'(a,3f10.5)') 'qp2 = ',qp2 284!write(6,*) 'vq = ',vq(ix,iy,iz) 285!write(6,*) 'jka = ',jka 286!write(6,*) 'jkk = ',jkk 287!write(6,*) 'jkb = ',jkb 288!write(6,*) 'igg = ',igg 289!write(6,*) 'ikpt = ',ikpt 290!write(6,*) 'ikptq = ',ikptq 291 omega(ix,iy,iz)=enrgy(indxkbnd(ikptq)+ibp)-ek+dw/2 292 if (iqq(1).eq.0.and.iqq(2).eq.0.and.iqq(3).eq.0) then 293 if (ipw1.eq.1.and.ipw2.eq.1.and.iband.eq.ibp) then 294! Find singular factor 295! call mkmatelX1(ibp,iband,ikpt,ikptq,ipwx(:,ipw1),igg, & 296!& ncg,nkpt,npwt,igmx,igmn,igndx, & 297!& isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, & 298!& lvtrans(1:3,ikk(1),ikk(2),ikk(3)), & 299!& cg,indxkcg,indxkpw,npwarr,kg, & 300!& cmatel) 301! if (ipaw.ne.0) then 302! call mkPAWmatelX1(pi,ibp,iband,ikk,jkk,qq,igg0, & 303!& pwmatel,tpwmatel,nbcore,ncband,amatel) 304! else 305! amatel2=(0.d0,0.d0) 306! endif 307! smat2=4.d0*pi*dble((cmatel+amatel)*conjg(cmatel+amatel)) 308! if (smat2.ne.0.d0) then 309! lqsing=.true. 310! else 311! lqsing=.false. 312! endif 313 smat2=4.d0*pi 314 lqsing=.true. 315 else 316 smat2=0.d0 317 lqsing=.false. 318 endif 319! If not singular, find tensor response at q=0 320 if (.not.lqsing) then 321 if (ipwx(1,ipw1).eq.0.and.ipwx(2,ipw1).eq.0.and.ipwx(3,ipw1).eq.0) then 322 call mkmatelJ1(ibp,iband,ikpt,ikptq,ipwx(:,ipw1),igg, & 323& ncg,nkpt,npwt,igmx,igmn,igndx, & 324& isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, & 325& lvtrans(1:3,ikk(1),ikk(2),ikk(3)), & 326& cg,indxkcg,indxkpw,npwarr,kg, & 327& jmatel) 328 if (ipaw.ne.0) then 329 call mkPAWmatelJ1(pi,ibp,iband,ikk,jkk, & 330& pwjmatel,tpwjmatel,nbcore,ncband,vmatel) 331 jmatel(:) = jmatel(:) + vmatel(:) 332 endif 333 do ii=1,3 334 jmatel(ii) = jmatel(ii)/(enrgy(indxkbnd(ikptq)+ibp)-ek) 335 enddo 336 else 337 call mkmatelX1(ibp,iband,ikpt,ikptq,ipwx(:,ipw1),igg, & 338& ncg,nkpt,npwt,igmx,igmn,igndx, & 339& isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, & 340& lvtrans(1:3,ikk(1),ikk(2),ikk(3)), & 341& cg,indxkcg,indxkpw,npwarr,kg, & 342& cmatel) 343 if (ipaw.ne.0) then 344 call mkPAWmatelX1(pi,ibp,iband,ikk,jkk,qq,igg0, & 345& pwmatel,tpwmatel,nbcore,ncband,amatel) 346 else 347 amatel=(0.d0,0.d0) 348 endif 349 jmatel = (cmatel+amatel)*qq/qq2 350 endif 351 if (ipw1.eq.ipw2) then 352 jmatel2=jmatel 353 else 354 if (ipwx(1,ipw2).eq.0.and.ipwx(2,ipw2).eq.0.and.ipwx(3,ipw2).eq.0) then 355 call mkmatelJ1(ibp,iband,ikpt,ikptq,ipwx(:,ipw2),igh, & 356& ncg,nkpt,npwt,igmx,igmn,igndx, & 357& isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, & 358& lvtrans(1:3,ikk(1),ikk(2),ikk(3)), & 359& cg,indxkcg,indxkpw,npwarr,kg, & 360& jmatel2) 361 if (ipaw.ne.0) then 362 call mkPAWmatelJ1(pi,ibp,iband,ikk,jkk, & 363& pwjmatel,tpwjmatel,nbcore,ncband,vmatel) 364 jmatel2(:) = jmatel2(:) + vmatel(:) 365 endif 366 do ii=1,3 367 jmatel2(ii) = jmatel2(ii)/(enrgy(indxkbnd(ikptq)+ibp)-ek) 368 enddo 369 else 370 call mkmatelX1(ibp,iband,ikpt,ikptq,ipwx(:,ipw2),igh, & 371& ncg,nkpt,npwt,igmx,igmn,igndx, & 372& isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, & 373& lvtrans(1:3,ikk(1),ikk(2),ikk(3)), & 374& cg,indxkcg,indxkpw,npwarr,kg, & 375& cmatel2) 376 if (ipaw.ne.0) then 377 call mkPAWmatelX1(pi,ibp,iband,ikk,jkk,qp,igg0, & 378& pwmatel,tpwmatel,nbcore,ncband,amatel2) 379 else 380 amatel2=(0.d0,0.d0) 381 endif 382 jmatel2 = (cmatel2+amatel2)*qp/qp2 383 endif 384 endif 385!! DEBUG 386!write(6,*) "jmatel" 387!do ii=1,3 388!write(6,*) jmatel(ii) 389!enddo 390!write(6,*) "jmatel2" 391!do ii=1,3 392!write(6,*) jmatel2(ii) 393!enddo 394!read(*,*) 395 do ii=1,3 396 do jj=1,3 397 jmat2(ii,jj)=4.d0*pi*jmatel(ii)*conjg(jmatel2(jj)) 398 enddo 399 enddo 400 else 401 do ii=1,3 402 do jj=1,3 403 jmat2(ii,jj)=(0.d0,0.d0) 404 enddo 405 enddo 406 endif 407 else 408 call mkmatelX1(ibp,iband,ikpt,ikptq,ipwx(:,ipw1),igg, & 409& ncg,nkpt,npwt,igmx,igmn,igndx, & 410& isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, & 411& lvtrans(1:3,ikk(1),ikk(2),ikk(3)), & 412& cg,indxkcg,indxkpw,npwarr,kg, & 413& cmatel) 414 if (ipaw.ne.0) then 415 call mkPAWmatelX1(pi,ibp,iband,ikk,jkk,qq,igg0, & 416& pwmatel,tpwmatel,nbcore,ncband,amatel) 417 else 418 amatel=(0.d0,0.d0) 419 endif 420 if (ipw1.eq.ipw2) then 421 cmatel2=cmatel 422 amatel2=amatel 423 else 424 call mkmatelX1(ibp,iband,ikpt,ikptq,ipwx(:,ipw2),igh, & 425& ncg,nkpt,npwt,igmx,igmn,igndx, & 426& isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, & 427& lvtrans(1:3,ikk(1),ikk(2),ikk(3)), & 428& cg,indxkcg,indxkpw,npwarr,kg, & 429& cmatel2) 430 if (ipaw.ne.0) then 431 call mkPAWmatelX1(pi,ibp,iband,ikk,jkk,qp,igg0, & 432& pwmatel,tpwmatel,nbcore,ncband,amatel2) 433 else 434 amatel2=(0.d0,0.d0) 435 endif 436 endif 437 xmat2(ix,iy,iz)=(cmatel+amatel)*conjg(cmatel2+amatel2) 438 vqmat2(ix,iy,iz)=xmat2(ix,iy,iz)*vq(ix,iy,iz) 439 endif 440 enddo 441 enddo 442 enddo 443 444 445! Alternate way of getting jmat2 446! Interpolation of density matrix elements / q^2 to q=0 447! For equally spaced samples A(-2) = -2 q, A(-1) = -1 q, A(1) = 1 q, A(2) = 2 q 448! average of two quadratic interpolations gives (see Numerical Recipes 3.1) 449! A(0) = (2/3) A(-1) + (2/3) A(1) - (1/6) A(-2) - (1/6) A(2) 450! Expect vqmat2 to converge to well defined value at q = 0, interpolate this 451! sum(ii,jj) qq(ii) * (bmet * jmat2 * bmet)(ii,jj) * qq(jj) / qq^2 = vqmat 452! define rjmat2 = bmet * jmat2 * bmet 453! find rjmat2 as q->0, then 454! jmat2(q->0) = bmet^(-1)*rjmat2(q->0)*bmet^(-1) 455! 456! note qq * rjmat2 * qq = rjmat2(1,1)*qq(1)^2 + rjmat2(2,2)*qq(2)^2 + rjmat2(3,3)*qq(3)^2 + (rjmat2(1,2)+rjmat2(2,1))*qq(1)*qq(2) + (rjmat2(1,3)+rjmat2(3,1))*qq(1)*qq(3) + (rjmat2(2,3)+rjmat2(3,2))*qq(2)*qq(3) 457! so longitudinal response only depends on sum of symmetric off-diagonal terms, 458! independent of difference. 459! Thus, we can safely set rjmat2(ii,jj) = (rjmat2(ii,jj)+rjmat2(jj,ii))/2 460 if (.true.) then 461! first do diagonal parts 462! e.g., along direction 1 463! qq = qq(1) * blat(1,:) 464! qq^2 = qq(1)*qq(1) * dot_product(blat(1,:),blat(1,:)) 465! sum(ii,jj) qq(ii) * rjmat2(ii,jj) * qq(jj) / qq^2 466! = (1/qq^2) * rjmat2(1,1) * qq(1)*qq(1) 467! = rjmat2(1,1) / dot_product(blat(1,:),blat(1,:)) 468! -> rjmat2(1,1) = dot_product(blat(1,:),blat(1,:)) * [response in direction 1] 469 do ii=1,3 470 ixx = (/0,0,0/) 471 ixx(ii) = 1 472 ixp2 = ictr+2*ixx 473 ixp1 = ictr+ixx 474 ixm1 = ictr-ixx 475 ixm2 = ictr-2*ixx 476 qq2 = dot_product(blat(ii,:),blat(ii,:)) 477 rjmatdum = 2.d0/3.d0 * vqmat2(ixp1(1),ixp1(2),ixp1(3)) & 478& + 2.d0/3.d0 * vqmat2(ixm1(1),ixm1(2),ixm1(3)) & 479& - 1.d0/6.d0 * vqmat2(ixp2(1),ixp2(2),ixp2(3)) & 480& - 1.d0/6.d0 * vqmat2(ixm2(1),ixm2(2),ixm2(3)) 481 rjmat2(ii,ii) = qq2 * rjmatdum 482 enddo 483! find off-diagonals using qq along directions (1,1,0), (1,0,1), and (0,1,1) 484! e.g., along direction (1,1,0) 485! qq = qq(1) * blat(1,:) + qq(2) * blat(2,:) 486! qq^2 = qq(1)*qq(1) * dot_product(blat(1,:),blat(1,:)) 487! + qq(2)*qq(2) * dot_product(blat(2,:),blat(2,:)) 488! + 2*qq(1)*qq(2) * dot_product(blat(1,:),blat(2,:)) 489! = dot_product(blat(1,:),blat(1,:)) + dot_product(blat(2,:),blat(2,:)) 490! + 2*dot_product(blat(1,:),blat(2,:)) 491! sum(ii,jj) qq(ii) * rjmat2(ii,jj) * qq(jj) / qq^2 492! = (1/qq^2) * ( rjmat2(1,1) * qq(1)*qq(1) + rjmat2(2,2) * qq(2)*qq(2) 493! + qq(1)*qq(2) * (rjmat2(1,2) + rjmat2(2,1)) ) 494! = (1/qq^2) * ( rjmat2(1,1) * qq(1)*qq(1) + rjmat2(2,2) * qq(2)*qq(2) 495! + 2 * rjmat2(1,2) qq(1)*qq(2)) 496! = (1/qq^2) * ( rjmat2(1,1) + rjmat2(2,2) + 2 * rjmat2(1,2) ) 497! -> rjmat2(1,2) = (qq^2 * [response in direction 1] - (rjmat2(1,1) + rjmat2(2,2)))/2 498 do ii=1,3 499 ixx = (/1,1,1/) 500 ixx(ii) = 0 501 ixp2 = ictr+2*ixx 502 ixp1 = ictr+ixx 503 ixm1 = ictr-ixx 504 ixm2 = ictr-2*ixx 505 qq2 = dot_product(blat(jj,:),blat(jj,:)) & 506& + dot_product(blat(kk,:),blat(kk,:)) & 507& + dot_product(blat(jj,:),blat(kk,:))*2 508 rjmatdum = 2.d0/3.d0 * vqmat2(ixp1(1),ixp1(2),ixp1(3)) & 509& + 2.d0/3.d0 * vqmat2(ixm1(1),ixm1(2),ixm1(3)) & 510& - 1.d0/6.d0 * vqmat2(ixp2(1),ixp2(2),ixp2(3)) & 511& - 1.d0/6.d0 * vqmat2(ixm2(1),ixm2(2),ixm2(3)) 512 jj = mod(ii,3)+1 ! note ii = mod(ii-1,3)+1 513 kk = mod(ii+1,3)+1 514 rjmat2(jj,kk) = (qq2*rjmatdum - (rjmat2(jj,jj)+rjmat2(kk,kk)))/2.d0 515 rjmat2(kk,jj) = rjmat2(jj,kk) 516!! DEBUG 517!qq = ixx/dble(ngkpt) 518!qq2=0.d0 519!do kk=1,3 520!do jj=1,3 521!qq2=qq2+qq(kk)*bmet(kk,jj)*qq(jj) 522!enddo 523!enddo 524!cvecdum = bmet(:,1)*qq(1)+bmet(:,2)*qq(2)+bmet(:,3)*qq(3) 525!cvecdum2 = jmat2(:,1)*cvecdum(1)+jmat2(:,2)*cvecdum(2)+jmat2(:,3)*cvecdum(3) 526!cvecdum = bmet(:,1)*cvecdum2(1)+bmet(:,2)*cvecdum2(2)+bmet(:,3)*cvecdum2(3) 527!write(6,*) "qq: ",qq 528!write(6,*) (cvecdum(1)*qq(1)+cvecdum(2)*qq(2)+cvecdum(3)*qq(3))/qq2 529!write(6,*) ixx 530!write(6,*) vqmat2(ixp2(1),ixp2(2),ixp2(3)) 531!write(6,*) vqmat2(ixp1(1),ixp1(2),ixp1(3)) 532!write(6,*) rjmatdum 533!write(6,*) vqmat2(ixm1(1),ixm1(2),ixm1(3)) 534!write(6,*) vqmat2(ixm2(1),ixm2(2),ixm2(3)) 535 enddo 536! DEBUG 537!do ii=1,3 538!do jj=1,3 539!if (ii.eq.jj) then 540!rjmat2(ii,jj) = 2.9357115896d0 541!else 542!rjmat2(ii,jj) = 0.d0; 543!endif 544!enddo 545!enddo 546 bmet_t=bmet 547 call inverse(bmet_t,bmetinv,3,3) 548 do ii=1,3 549 do jj=1,3 550 jmat2(ii,jj) = (0.d0,0.d0) 551 do kk=1,3 552 do ll=1,3 553 jmat2(ii,jj) = jmat2(ii,jj) + bmetinv(ii,kk)*rjmat2(kk,ll)*bmetinv(ll,jj) 554 enddo 555 enddo 556 enddo 557 enddo 558 endif 559!! DEBUG 560!write(6,*) "rjmat2" 561!do ii=1,3 562!write(6,*) dble(rjmat2(:,ii)) 563!enddo 564!write(6,*) 565!do ii=1,3 566!write(6,*) dimag(rjmat2(:,ii)) 567!enddo 568!write(6,*) 569!write(6,*) "jmat2" 570!do ii=1,3 571!write(6,*) dble(jmat2(:,ii)) 572!enddo 573!write(6,*) 574!do ii=1,3 575!write(6,*) dimag(jmat2(:,ii)) 576!enddo 577!write(6,*) 578!write(6,*) "in Cartesian coordinates" 579!do ii=1,3 580!do jj=1,3 581!ctensdum(ii,jj) = (0.d0,0.d0) 582!do kk=1,3 583!do ll=1,3 584!ctensdum(ii,jj) = ctensdum(ii,jj) + blat(kk,ii)*jmat2(kk,ll)*blat(ll,jj) 585!enddo 586!enddo 587!enddo 588!enddo 589!do ii=1,3 590!write(6,*) dble(ctensdum(:,ii)) 591!enddo 592!write(6,*) 593!do ii=1,3 594!write(6,*) dimag(ctensdum(:,ii)) 595!enddo 596!stop 597! DEBUG 598!write(6,*) "Matrix elements computed" 599 600 if (itetrahedron.eq.0) then 601 write(6,*) "Sum over points not yet implemented in self energy calculation" 602 write(6,*) "set itetrahedron=1 and run again" 603 else 604! tetrahedron integration 605 606! Don't waste time integrating over energies with no contribution 607 do ix=1,ngkpt(1) 608 do iy=1,ngkpt(2) 609 do iz=1,ngkpt(3) 610 ixx=(/ix,iy,iz/) 611 call fhilo(omega,ixx,ngkpt,whi,wlo) 612 if (ix.eq.1.and.iy.eq.1.and.iz.eq.1) then 613 wwhi=whi 614 wwlo=wlo 615 else 616 wwhi=max(wwhi,whi) 617 wwlo=min(wwlo,wlo) 618 endif 619 enddo 620 enddo 621 enddo 622 623 iwlo=nint(wwlo*dble(nwpt)/wmax) 624 iwhi=nint(wwhi*dble(nwpt)/wmax) 625 if (ibp.gt.nbocc) then 626 ielo=iwlo+1 627 iehi=iwhi+nwpt 628 else 629 ielo=iwlo-nwpt 630 iehi=iwhi-1 631 endif 632 allocate (ssi(ielo:iehi)) 633 do ie=ielo,iehi 634 ssi(ie)=0.d0 635 enddo 636 ssx=0.d0 637 638 if (lc) then 639 do iw=1,nwpt 640 do ix=1,ngkpt(1) 641 do iy=1,ngkpt(2) 642 do iz=1,ngkpt(3) 643 ixx=(/ix,iy,iz/) 644 iqpt=iqndx(ixx(1),ixx(2),ixx(3)) 645 call locateelement(ixx,ipw1,ipw2,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpw) 646 call locateelement(ixx,ipw2,ipw1,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpwt) ! transpose matrix coordiantes, to find anti-Hermitian part of lossfn 647 if (iqpt.eq.iqctr) then 648 vfactor(ixx(1),ixx(2),ixx(3),iw)=(0.d0,0.d0) 649 else if (jjpw.ne.0.and.jjpwt.ne.0) then 650 vfactor(ixx(1),ixx(2),ixx(3),iw)= & 651& vqmat2(ixx(1),ixx(2),ixx(3))* & 652& (lossfn(iw,iqpt,jjpw)-conjg(lossfn(iw,iqpt,jjpwt)))*0.5d0 653 else 654 vfactor(ixx(1),ixx(2),ixx(3),iw)=(0.d0,0.d0) 655 endif 656 enddo 657 enddo 658 enddo 659!write(6,*) "iw = ", iw 660 call locateelement(ictr,ipw1,ipw2,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpw) 661 call locateelement(ictr,ipw2,ipw1,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpwt) 662 do ii=1,3 663 do jj=1,3 664 if (jjpw.ne.0.and.jjpwt.ne.0) then 665 vjfactor(ii,jj,iw)=(0.d0,0.d0) 666 do kk=1,3 667 do ll=1,3 668 iiq=ii + 3*(kk-1) 669 vjfactor(ii,jj,iw)= & 670& vjfactor(ii,jj,iw) + & 671& 0.5d0*(lossfn(iw,nqpt+iiq,jjpw)-conjg(lossfn(iw,nqpt+iiq,jjpwt))) & 672& *bmet(kk,ll)*jmat2(ll,jj) 673 enddo 674 enddo 675 iiq=ii + 3*(jj-1) 676 vsfactor(ii,jj,iw)=smat2*(lossfn(iw,nqpt+iiq,jjpw)-conjg(lossfn(iw,nqpt+iiq,jjpwt)))*0.5d0 677 else 678 vjfactor(ii,jj,iw)=(0.d0,0.d0) 679 vsfactor(ii,jj,iw)=(0.d0,0.d0) 680 endif 681 enddo 682 enddo 683!do ii=1,3 684!write(6,*) dble(vjfactor(ii,:,iw)) 685!enddo 686!write(6,*) 687!do ii=1,3 688!write(6,*) dimag(vjfactor(ii,:,iw)) 689!enddo 690!write(6,*) 691!do ii=1,3 692!write(6,*) dble(vsfactor(ii,:,iw)) 693!enddo 694!write(6,*) 695!do ii=1,3 696!write(6,*) dimag(vsfactor(ii,:,iw)) 697!enddo 698!read(*,*) 699 do ii=1,3 700 do jj=1,3 701 vsm1(ii,jj,iw)=(0.d0,0.d0) 702 do kk=1,3 703 vsm1(ii,jj,iw)=vsm1(ii,jj,iw)+bmet(ii,kk)*vsfactor(kk,jj,iw) 704 enddo 705 enddo 706 enddo 707 do ii=1,3 708 do jj=1,3 709 vsm2(ii,jj,iw)=(0.d0,0.d0) 710 do kk=1,3 711 vsm2(ii,jj,iw)=vsm2(ii,jj,iw)+vsm1(ii,kk,iw)*bmet(kk,jj) 712 enddo 713 enddo 714 enddo 715 enddo 716 endif 717 718! DEBUG 719!write(6,*) "Beginning integration" 720 721! Do the integration 722 do ix=1,ngkpt(1) 723 do iy=1,ngkpt(2) 724 do iz=1,ngkpt(3) 725 ixx=(/ix,iy,iz/) 726 iqq=ixx-ictr 727 ssx1=ssx 728 call fval(omega,ixx,ixxp,ngkpt,enval) 729 if (lqsing) call fval(vq,ixx,ixxp,ngkpt,vq2) 730 if (lc) then 731 do iw=1,nwpt 732 call fpol(vfactor(:,:,:,iw),ngkpt,ixx,ixxp,vfv(:,iw)) 733 enddo 734 endif 735 if (lx) call fpol(vqmat2(:,:,:),ngkpt,ixx,ixxp,vfx) 736 do itet=1,6 737 centerint = .false. 738 do iv=1,4 739 evtx0(iv)=enval(ivndx(iv,itet)) 740 rrpyr(1:3,iv)=rr(1:3,ivndx(iv,itet)) 741 if (lc) then 742 do iw=1,nwpt 743 vpyr0(iv,iw)=vfv(ivndx(iv,itet),iw) 744 enddo 745 endif 746 if (lx) xpyr0(iv)=vfx(ivndx(iv,itet)) 747 if ((ixx(1)+rrpyr(1,iv).eq.ictr(1)).and.(ixx(2)+rrpyr(2,iv).eq.ictr(2)).and.(ixx(3)+rrpyr(3,iv).eq.ictr(3))) then 748 centerint = .true. ! does tetrahedron contain center vertex? 749 vcenterint = iv ! index of center vertex 750 endif 751 enddo 752!write(6,*) ixx 753!write(6,*) itet,centerint 754!if (centerint) then 755!write(6,*) itet 756!do iv=1,4 757!write(6,*) ixx+rrpyr(:,iv) 758!enddo 759!!read(*,*) 760!endif 761 call indxhpsort(4,4,evtx0,indxe) 762 evtx=evtx0(indxe) 763 do ii=1,4 764 jj=indxe(ii) 765 kvtx(:,ii) = (iqq(:) + rrpyr(:,jj)) 766 do kk=1,3 767 kvtxc(kk,ii)=dot_product(iqq+rrpyr(:,jj),qkcvt(:,kk)) 768 enddo 769 enddo 770 xk(1:3,1)=(/0.d0,0.d0,0.d0/) 771 xkc(1:3,1)=(/0.d0,0.d0,0.d0/) 772 do ii=1,3 773 xk(:,ii+1)=rrpyr(:,indxe(ii+1))-rrpyr(:,indxe(1)) 774 enddo 775 do ii=2,4 776 xkc(1:3,ii)=kvtxc(1:3,ii)-kvtxc(1:3,1) 777 enddo 778 do jj=1,3 ! make contragredient, reduced coordinates 779 call cross(xk(:,mod(jj,3)+2),xk(:,mod(jj+1,3)+2),rg(:,jj)) 780 enddo 781 do jj=1,3 ! make contragredient, Cartesian coordinates 782 call cross(xkc(:,mod(jj,3)+2),xkc(:,mod(jj+1,3)+2),rgc(:,jj)) 783 enddo 784 tvol=dot_product(xk(1:3,2),rg(1:3,1)) 785 rg=rg/tvol 786 tvol=abs(tvol) 787 tvolc=dot_product(xkc(1:3,2),rgc(1:3,1)) 788 rgc=rgc/tvolc 789 tvolc=abs(tvolc) 790 iwwlim=nint(evtx/dw) 791 if (lqsing .and. centerint) then 792! singular part 793! DEBUG 794!xdum = 0.d0 795!write(6,*) 796!write(6,*) "singular integration" 797!write(6,*) 'ixx: ',ixx 798!write(6,*) 'itet: ',itet,' centerint: ',centerint 799!write(6,*) 'iwwlim: ',iwwlim 800 do iv=1,4 801! subtraction of integrated singular function from values at vertices 802 ixr=ixx+rrpyr(:,iv) 803 iqr = ixr - ictr 804 qr(:,iv) = dble(iqr) 805 if (ixr(1).eq.ictr(1).and.ixr(2).eq.ictr(2).and.ixr(3).eq.ictr(3)) then 806 if (lc) then 807 do iw=1,nwpt 808 vpyr0(iv,iw) = 0.d0 809 enddo 810 endif 811 if (lx) then 812 xpyr0(iv) = 0.d0 813 endif 814 else 815 qq = qr(:,iv)/dble(ngkpt) 816 qqr(:,iv)=qq 817 qq2 = 0.d0 818 do ii=1,3 819 do jj=1,3 820 qq2 = qq2 + qq(ii)*bmet(ii,jj)*qq(jj) 821 enddo 822 enddo 823 if (lc) then 824 do iw=1,nwpt 825 do ii=1,3 826 do jj=1,3 827 vpyr0(iv,iw)=vpyr0(iv,iw)-qq(ii)*vsm2(ii,jj,iw)*qq(jj)/(qq2**2) 828 enddo 829 enddo 830 enddo 831 endif 832 if (lx) then 833 xpyr0(iv)=xpyr0(iv)-smat2/qq2 834 endif 835 endif 836 enddo 837 bgrad=(/0.d0,0.d0,0.d0/) ! energy gradient 838 do jj=1,3 839 bgrad=bgrad+(evtx(jj+1)-evtx(1))*rgc(:,jj) 840 enddo 841 xmult=1.d0/sqrt(dot_product(bgrad,bgrad)) 842 if (lc) then 843! DEBUG 844!cdum = 0.d0 845 do iww=iwwlim(1)-1,iwwlim(4)+1 846 wwwlo=dw*(iww-0.5d0) 847 call singular_adaptive_tetint(qqr,evtx,wwwlo,dw,vcenterint,bmet,abr,rlr,rlr0,stens) 848 do iw=1,nwpt 849 ie=iww-ocsign*iw 850 if (ie.ge.ielo.and.ie.le.iehi) then 851 do ii=1,3 852 do jj=1,3 853 ssi(ie)=ssi(ie)+vsm2(ii,jj,iw)*stens(jj,ii) 854! DEBUG 855!if (ie.eq.0) cdum = cdum+vsm2(ii,jj,iw)*stens(jj,ii) 856 enddo 857 enddo 858 endif 859 enddo 860 enddo 861!write(6,*) cdum 862 endif 863 if (lx) then 864 wwwlo = min(evtx(1),evtx(2),evtx(3),evtx(4)) 865 dw = max(evtx(1),evtx(2),evtx(3),evtx(4))-wwwlo 866!xdum = time() 867 call singular_adaptive_tetint(qqr,evtx,wwwlo,dw,vcenterint,bmet,abr,rlr,rlr0,stens) 868!! DEBUG 869!cdum = 0.d0 870 do ii=1,3 871 do jj=1,3 872 ssx=ssx+smat2*bmet(ii,jj)*stens(jj,ii) 873!! DEBUG 874!cdum = cdum+smat2*bmet(ii,jj)*stens(jj,ii) 875 enddo 876 enddo 877!write(6,*) cdum 878 endif 879!write(6,*) 'Khst7',xdum 880!write(32,'(4i3,2x,f14.10,2x,f14.10," sing")') ixx,itet,(cdum/2.d0) 881!write(6,*) "finished singular integration" 882 elseif (centerint) then 883! anisotropic part 884! DEBUG 885!cdum = 0.d0 886 if (lc) then 887 do iww=iwwlim(1)-1,iwwlim(4)+1 888 wwwlo=dw*(iww-0.5d0) 889! wwwhi=dw*(iww+0.5d0) 890 do iw=1,nwpt 891 call vnitetint(rrpyr,evtx0,wwwlo,dw,vjfactor(:,:,iw),bmet,vcenterint,cint) 892 ie=iww-ocsign*iw 893 if (ie.ge.ielo.and.ie.le.iehi) ssi(ie)=ssi(ie)+cint/volelmnt 894! DEBUG 895!if (ie.eq.0.or.ie.eq.1) cdum = cdum+cint/volelmnt 896 enddo 897 enddo 898 endif 899 if (lx) then 900 call vnitetint(rrpyr,evtx0,evtx0(indxe(1)),evtx0(indxe(4))-evtx0(indxe(1)),jmat2,bmet,vcenterint,cint) 901 ssx=ssx+cint/volelmnt 902 endif 903! DEBUG 904!write(32,'(4i3,2x,f14.10,2x,f14.10," cent")') ixx,itet,(cdum/2.d0) 905 endif 906! can use basic tetrahedron integration for everything else 907! DEBUG 908!cdum = 0.d0 909 if (lc) then 910 do iw=1,nwpt 911 aa0(iw)=vpyr0(indxe(1),iw) ! base value of function 912 av(1:3,iw)=(/0.d0,0.d0,0.d0/) ! gradient of function 913 do jj=1,3 914 av(1:3,iw)=av(1:3,iw) & 915& +(vpyr0(indxe(jj+1),iw)-vpyr0(indxe(1),iw))*rg(1:3,jj) 916 enddo 917 enddo 918 do iww=iwwlim(1)-1,iwwlim(4)+1 919 wwwlo=dw*(iww-0.5d0) 920! wwwhi=dw*(iww+0.5d0) 921 call mkvsint(rrpyr(1:3,indxe),xk,volred,evtx,wwwlo,dw,sint1,svec) 922 do iw=1,nwpt 923 ie=iww-ocsign*iw 924 if (ie.ge.ielo.and.ie.le.iehi) ssi(ie)=ssi(ie)+(aa0(iw)*sint1+dot_product(svec,av(:,iw)))/volelmnt 925! DEBUG 926!if (ie.eq.0.or.ie.eq.1) cdum = cdum+(aa0(iw)*sint1+dot_product(svec,av(:,iw)))/volelmnt 927 enddo 928 enddo 929 endif 930! DEBUG 931!write(32,'(4i3,2x,f14.10,2x,f14.10)') ixx,itet,(cdum/2.d0) 932 if (lx) then 933 xme=xpyr0(indxe(1)) 934 xv=(/0.d0,0.d0,0.d0/) 935 do jj=1,3 936 xv=xv+(xpyr0(indxe(jj+1))-xpyr0(indxe(1)))*rg(:,jj) 937 enddo 938 call mkvsint(rrpyr(1:3,indxe),xk,volred,evtx,evtx(1),evtx(4)-evtx(1),sint1,svec) 939 ssx=ssx+(xme*sint1+dot_product(svec,xv))/volelmnt 940! DEBUG 941!write(32,'(4i3,2x,f16.10)') ixx,itet,(xme*sint1+dot_product(svec,xv))/volelmnt 942!write(32,'(4i3,2x,f14.10,2x,3f14.10)') ixx,itet, & 943!& (xme*sint1+dot_product(svec,xv))/volelmnt, & 944!& (xpyr0(indxe(1))+xpyr0(indxe(2))+xpyr0(indxe(3))+xpyr0(indxe(4))) & 945!& /(volelmnt*6.d0*4.d0) 946 endif 947 enddo 948! DEBUG 949!write(32,'(3i3,2x,4f16.10)') ixx,ssi(0),ssi(1) 950 enddo 951 enddo 952 enddo 953! if (lc) ssi=-ssi/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol*pi) 954! if (lx) ssx=-ssx/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol) 955 if (lc) ssi=-ocsign*ssi/((2*pi)**3*pi) 956 if (lx) ssx=-ssx/((2*pi)**3) 957 958! DEBUG 959!write(6,*) "Integration complete" 960 961 if (lc) then 962 do ie=-nwpt,nwpt 963 eps1=dble(ie)*dw-dw/2 964! ssc(ie)=(0.d0,0.d0) 965 do je=ielo,iehi 966 if (ie.eq.je) then 967 ssc(ie)=ssc(ie)+pi*ssi(je) 968 else 969 eps2=dble(je)*dw-dw/2 970 ssc(ie)=ssc(ie)-(0.d0,1.d0)*ocsign*ssi(je)*dw/(eps2-eps1) 971 endif 972!if (ie.eq.0) write(16,'(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) 973!if (ie.eq.0) write(16,'(2x,i4,2x,f10.3,4x,f10.6,2x,f10.6)') je,eps2*27.2114,ssi(je)*27.2114 974 enddo 975 enddo 976 endif 977!do je=ielo,iehi 978!write(6,*) je,ssi(je) 979!enddo 980 deallocate (ssi) 981!write(6,'("lstv",2x,f14.10,4x,"(",f14.10,",",f14.10,")")') ssx,(ssc(0)+ssc(1))/2. 982 983 if (lc) cse1=cse1+ssc 984 if (lx) xse=xse+ssx 985! if (lx.and.lc) then 986! write(6,'(i4,4x,i5,2i3,3x,3i3,5x,"(",f10.6,",",f10.6,")",5x,f10.6)') & 987!& ibp,iipw,ipw1,ipw2,ipwx(:,ipw1),((ssc(1)+ssc(0))/2)*27.2114,dble(ssx)*27.2114 988! elseif (lc) then 989! write(6,'(i4,4x,i5,2i3,3x,3i3,5x,"(",f10.6,",",f10.6,")")') ibp,iipw,ipw1,ipw2,ipwx(:,ipw1), & 990!& ((ssc(1)+ssc(0))/2)*27.2114 991! else 992! write(6,'(i4,4x,i5,2i3,3x,3i3,28x,5x,f10.6)') & 993!& ibp,iipw,ipw1,ipw2,ipwx(:,ipw1),dble(ssx)*27.2114 994! endif 995 endif 996 997 enddo 998! if (lx) then 999! write(6,'(i4,5x,"(",f10.6,",",f10.6,")",5x,f10.6)') ibp, & 1000!& (cse1(1)+cse1(0)-cse2(1)-cse2(0))*27.2114/2.d0, & 1001!& dble(xse-xse2)*27.2114 1002! else 1003! write(6,'(i4,5x,"(",f10.6,",",f10.6,")",f10.6)') ibp, & 1004!& (cse1(1)+cse1(0)-cse2(1)-cse2(0))*27.2114/2.d0 1005! endif 1006enddo 1007 1008!do ie=-nwpt,nwpt 1009! eps1=dble(ie)*dw-dw/2 1010! write(12,'(i4,f10.3,2(3x,2es12.3))') ie,eps1*27.2114,cse1(ie)*27.2114 1011!enddo 1012if (itetrahedron.eq.1) then 1013 cse=(cse1(0)+cse1(1))/2 1014 dcse=(cse1(1)-cse1(0))/dw 1015endif 1016zz=1.d0/(1.d0-dcse) 1017!write(6,*) cse*27.2114 1018!write(6,*) xse*27.2114 1019!write(6,*) zz 1020!write(6,*) dble(zz)*dimag(cse)*27.2114 1021 1022return 1023end subroutine mkse 1024 1025!**************************************************************************** 1026! evaluate Integral_wmin^wmax dww tsingint 1027 1028subroutine tsinggrater(iv,jv,wmin,wmax,ive,xk,xkvtx,ngkpt,evtx,abr,rlr,nsing,wsing,sint) 1029implicit none 1030integer :: iv,jv,numcal,maxns,mx,nstack,ii,jj,kk,icount,ive,nsing 1031double precision :: wsing(20) 1032parameter (mx=1500) 1033integer :: ngkpt(3) 1034double precision :: xk(3,4),xkvtx(3),evtx(4),sint 1035double precision :: value,valu,fval(3,mx) 1036double precision :: wmin,wmax,del,del1,dif,frac 1037double precision :: abr,rlr,error 1038double precision :: wleft(mx),dw(3),wt(3),ww 1039double precision :: wt9(9) 1040logical :: atsing 1041save dw,wt,wt9 1042data dw /0.1127016653792583,0.5,0.8872983346207417/ 1043data wt /0.277777777777777778,0.4444444444444444444,0.2777777777777777778/ 1044data wt9 /0.0616938806304841571,0.108384229110206161, & 1045& 0.0398463603260281088,0.175209035316976464, & 1046& 0.229732989232610220,0.175209035316976464, & 1047& 0.0398463603260281088,0.108384229110206161, & 1048& 0.0616938806304841571 / 1049integer :: i_prt_DEBUG, j_prt_DEBUG, k_prt_DEBUG, & 1050& i_ct_DEBUG, j_ct_DEBUG, k_ct_DEBUG, l_ct_DEBUG, & 1051& k_ref_DEBUG, i_ref_DEBUG 1052common /tsingint_DEBUG/ i_prt_DEBUG, j_prt_DEBUG, k_prt_DEBUG, & 1053& i_ct_DEBUG, j_ct_DEBUG, k_ct_DEBUG, l_ct_DEBUG, & 1054& k_ref_DEBUG, i_ref_DEBUG 1055 1056! nstack is the number of different intervals into which the 1057! integration region is currently divided. The number of regions can 1058! grow if more accuracy is needed by dividing the right-most region 1059! into three regions. The number of regions shrinks when the integral 1060! over the right-most region is accurate enough, in which case that 1061! integral is added to the total (stored in grater) and the region 1062! is removed from consideration (and a new region is the right-most) 1063 nstack=nsing+1 1064 maxns=nstack 1065 error=0.d0 1066 sint=0.d0 1067! The array xleft stores the boundary points of the regions. 1068! The singular points just govern the initial placement of the regions. 1069 wleft(1)=wmin 1070 wleft(nsing+2)=wmax 1071 if(nsing.gt.0) then 1072 do jj=1,nsing 1073 wleft(jj+1)=wsing(jj) 1074 enddo 1075 endif 1076! For each region, calculate the function and store at three selected points. 1077 do ii=1,nstack 1078 del=wleft(ii+1)-wleft(ii) 1079 do jj=1,3 1080 ww=wleft(ii)+del*dw(jj) 1081if (k_ct_DEBUG.eq.k_ref_DEBUG) i_ct_DEBUG = i_ct_DEBUG+1 1082!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "%%%%%AGS ",ww,i_ct_DEBUG 1083 call tsingint(iv,jv,ive,ww,xk,xkvtx,ngkpt,evtx,fval(jj,ii)) 1084 enddo 1085 enddo 1086 numcal=nstack*3 1087 do 1088 if(nstack+3.ge.mx) then 1089 write(6,*) 'tsinggrater: too many regions' 1090 stop 1091 endif 1092! Divide the rightmost region into three subregions. 1093 del=wleft(nstack+1)-wleft(nstack) 1094 wleft(nstack+3)=wleft(nstack+1) 1095 wleft(nstack+1)=wleft(nstack)+del*dw(1)*2.d0 1096 wleft(nstack+2)=wleft(nstack+3)-del*dw(1)*2.d0 1097! The three data points already found for the region become the 1098! middle data points (number 2 in first index of fval) for each region. 1099 fval(2,nstack+2)=fval(3,nstack) 1100 fval(2,nstack+1)=fval(2,nstack) 1101 fval(2,nstack)=fval(1,nstack) 1102! Now do the integral over the right-most region in two different ways- 1103! a three point integral (valu) over each of the three subregions 1104! and a more accurate nine-point integral (value) over whole region. 1105! valu is used only for the error estimate. 1106 icount=0 1107 value=0.d0 1108 valu=0.d0 1109 do jj=nstack,nstack+2 1110 del1=wleft(jj+1)-wleft(jj) 1111 ww=wleft(jj)+dw(1)*del1 1112if (k_ct_DEBUG.eq.k_ref_DEBUG) i_ct_DEBUG = i_ct_DEBUG+1 1113!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "%%%%%AGS ",ww,i_ct_DEBUG 1114 call tsingint(iv,jv,ive,ww,xk,xkvtx,ngkpt,evtx,fval(1,jj)) 1115 ww=wleft(jj)+dw(3)*del1 1116if (k_ct_DEBUG.eq.k_ref_DEBUG) i_ct_DEBUG = i_ct_DEBUG+1 1117!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "%%%%%AGS ",ww,i_ct_DEBUG 1118 call tsingint(iv,jv,ive,ww,xk,xkvtx,ngkpt,evtx,fval(3,jj)) 1119 numcal=numcal+2 1120!if (numcal.gt.1000) stop 1121 do kk=1,3 1122 icount=icount+1 1123 value=value+wt9(icount)*fval(kk,jj)*del 1124 valu=valu+fval(kk,jj)*wt(kk)*del1 1125 enddo 1126 enddo 1127 dif=abs(value-valu) 1128! If the following condition is true, add in this integral to the total, 1129! and reduce the number of regions under consideration. 1130 frac=del/(wmax-wmin) 1131 atsing=.false. 1132 if (frac.le.1.0d-8) atsing=.true. 1133 if (dif.le.abr*frac.or.dif.le.rlr*abs(value).or. & 1134& (atsing.and.(frac.le.1.0d-15.or.dif.le.abr*0.1))) then 1135 sint=sint+value 1136 error=error+abs(dif) 1137 nstack=nstack-1 1138! If no more regions, we are done. 1139 if(nstack.le.0) return 1140 else 1141! If the integration is insufficiently accurate, make each of the 1142! three subregions of the right-most region into regions. 1143! On next pass the right-most of these is the new current region. 1144 nstack=nstack+2 1145 maxns = max(maxns,nstack) 1146 endif 1147 enddo 1148end subroutine tsinggrater 1149 1150!*************************************************************************** 1151! 1152! Over a surface S defined by the triangle with vertices xkvtx+xkp(1:3,1) 1153! xkvtx+xkp(1:3,2) 1154! xkvtx+xkp(1:3,3) 1155! where xkp defined to give S of constant energy 1156! (see G. Lehmann and M. Taut, Phys. stat. sol. (b) 54 469 (1972)) 1157! and vector k measured from tip of xkvtx 1158! find sint = integral dS q(iv)q(jv)/q^4 1159! where q=k+xkvtx 1160! On S, define coordinates X and Y such that k=xkp(1:3,3)+X*xq1+Y*xq2 1161! normal=cross product(xq1,xq2) 1162! integral dS -> |normal|*integral^1_0 dX integral^{1-X}_0 dY 1163! integral^{1-X}_0 dY ((k(iv)+xkvtx(iv))*(k(jv)+xkvtx(jv)))/(dot_product(k+xkvtx,k+xkvtx))**2 = tfn1 1164! numerically integrate tfn1 1165! sint = integral^1_0 dX tfn1 1166 1167subroutine tsingint(iv,jv,ive,ww,xk,xkvtx,ngkpt,evtx,sint) 1168use geometry 1169implicit none 1170integer :: iv,jv,ive,ivv,jvv 1171integer :: ngkpt(3) 1172double precision :: xk(3,4),xkvtx(3),evtx(4),ww 1173double precision :: xkp(3,3),xq1(3),xq2(3),xkv0(3),xcvt,xkvmag(3) 1174double precision :: xkv1(3),xkv2(3) 1175double precision :: vt1(3),vt2(3),vt3(3) 1176double precision :: aa0,aa1,aa2,bb0,bb1,cc,dd0,dd1,dd2,ee0,ee1,ff,qq,delta,brd 1177double precision :: xmin,xmax,abr,rlr,xsing(30),error 1178double precision :: tfn1,xfn2,xfn3,grater 1179double precision :: sint,sint1,sint2,sint3 1180double precision :: root1,root2,root1a,root2a,xroot1,xroot2,xroot1a,xroot2a 1181double precision :: xnorm(3),area 1182double precision :: xdum,xmx,rlxk 1183double precision :: vtemp(3),utemp(3),wtemp(3),vmag2,umag2,wmag2 1184double precision polyterm(0:4) 1185integer :: ii,jj,kk,ll 1186integer :: nsing,numcal,maxns,nroots,nrootsa 1187integer :: iter 1188integer :: indxq(3) ! sorts vertices in ascending magnitude of q-vector 1189double precision :: xd01,xd02 ! xx value where delta is 0 1190double precision :: xdt11,xdt12,xdt21,xdt22 ! for xx = xd0 + dx, linear and quadratic term in dx for delta 1191double precision :: xx_dd1, xx_dd2 ! roots of dd = 0 1192double precision :: xx_def1, xx_def2 ! roots of dd + ee*zz + ff*zz^2 = 0 1193logical :: usexd1,usexd2,usexx_dd1,usexx_dd2 ! are roots in integration region 1194double precision :: denom,snangle,snangle1,snangle2,snangle3 1195double precision :: aat1,aat2,bbt1,bbt2,ddt1,ddt2,eet1,eet2,zzt1,zzt2,deltat1,deltat2,aa2ddt1,aa2ddt2 1196double precision :: RRA1,RRA2,RRB1,RRB2,RRC1,RRC2,d0e0f,d1e1e02f,d2e1f 1197double precision :: qqmag(2), q02k1, q22qdiff, qqdiff(3), qqdiffmag 1198double precision :: xi1,xi2,zeta12,eta0,eta2,nu2,kappa,lambda 1199double precision :: omxi12,omxi22,ometa22,omzeta122 1200double precision :: deltaprefact, d_angle 1201double precision :: rrootd,irootd,rrootdelta,irootdelta,rrootdezfzz,irootdezfzz 1202double precision :: gammapfroot1,gammapfroot2,gammapfpf 1203integer :: gammapfnroot 1204double precision :: denomroot1,denomroot2,denomroot3,denomroot4,denompf 1205integer :: denomnroot 1206common /tfn/ aa0,aa1,aa2,bb0,bb1,cc,dd0,dd1,dd2,ee0,ee1,ff,brd, & 1207& xkv0,xq1,ivv,jvv, & 1208& xd01,xd02,xdt11,xdt12,xdt21,xdt22,xx_dd1,xx_dd2,xx_def1,xx_def2, & 1209& usexd1,usexd2,usexx_dd1,usexx_dd2, & 1210& aat1,aat2,bbt1,bbt2,ddt1,ddt2,eet1,eet2,zzt1,zzt2, & 1211& deltat1,deltat2,aa2ddt1,aa2ddt2, qqmag, qqdiffmag, & 1212& q02k1, q22qdiff, deltaprefact, umag2, & 1213& rrootd,irootd,rrootdelta,irootdelta,rrootdezfzz,irootdezfzz, & 1214& gammapfroot1,gammapfroot2,denomroot1,denomroot2,denomroot3,denomroot4, & 1215& gammapfnroot,denomnroot, & 1216& gammapfpf,denompf, & 1217& RRA1,RRA2,RRB1,RRB2,RRC1,RRC2,d0e0f, & 1218& xi1,xi2,zeta12,eta0,eta2,kappa,lambda,omxi12,omxi22,ometa22,omzeta122 1219integer :: i_prt_DEBUG, j_prt_DEBUG, k_prt_DEBUG, & 1220& i_ct_DEBUG, j_ct_DEBUG, k_ct_DEBUG, l_ct_DEBUG, & 1221& k_ref_DEBUG, i_ref_DEBUG 1222common /tsingint_DEBUG/ i_prt_DEBUG, j_prt_DEBUG, k_prt_DEBUG, & 1223& i_ct_DEBUG, j_ct_DEBUG, k_ct_DEBUG, l_ct_DEBUG, & 1224& k_ref_DEBUG, i_ref_DEBUG 1225double precision :: sqrtm1 1226external tfn1,xfn2,xfn3,grater,sqrtm1 1227double precision :: xdummy, ydummy, zdummy,xxdummy,yydummy,zzdummy 1228integer :: idummy, jdummy, kdummy 1229double complex :: cdummy1,cdummy2 1230 1231! broadening: 1/|q|^2 -> 1/(|q|^2+delta^2) 1232! delta=min(dot_product(xk(:,2),xk(:,2)),dot_product(xk(:,3),xk(:,3)),dot_product(xk(:,4),xk(:,4)))*1.d-4 1233 ivv=iv 1234 jvv=jv 1235 if (ive.eq.1) then 1236 xkp(:,1)=(ww-evtx(1))*xk(:,2)/(evtx(2)-evtx(1)) 1237 xkp(:,2)=(ww-evtx(1))*xk(:,3)/(evtx(3)-evtx(1)) 1238 xkp(:,3)=(ww-evtx(1))*xk(:,4)/(evtx(4)-evtx(1)) 1239 elseif (ive.eq.2) then 1240 xkp(:,1)=(ww-evtx(2))*(xk(:,4)-xk(:,2))/(evtx(4)-evtx(2)) + xk(:,2) 1241 xkp(:,2)=(ww-evtx(2))*(xk(:,3)-xk(:,2))/(evtx(3)-evtx(2)) + xk(:,2) 1242 xkp(:,3)=(ww-evtx(1))*(xk(:,4)-xk(:,1))/(evtx(4)-evtx(1)) 1243 elseif (ive.eq.3) then 1244 xkp(:,1)=(ww-evtx(1))*xk(:,3)/(evtx(3)-evtx(1)) 1245 xkp(:,2)=(ww-evtx(2))*(xk(:,3)-xk(:,2))/(evtx(3)-evtx(2)) + xk(:,2) 1246 xkp(:,3)=(ww-evtx(1))*(xk(:,4)-xk(:,1))/(evtx(4)-evtx(1)) 1247 else 1248 xkp(:,1)=(ww-evtx(2))*(xk(:,4)-xk(:,2))/(evtx(4)-evtx(2)) + xk(:,2) 1249 xkp(:,2)=(ww-evtx(3))*(xk(:,4)-xk(:,3))/(evtx(4)-evtx(3)) + xk(:,3) 1250 xkp(:,3)=(ww-evtx(1))*xk(:,4)/(evtx(4)-evtx(1)) 1251 endif 1252 do kk=1,3 1253 xkvmag(kk) = 0.d0 1254 do ii=1,3 1255 do jj=1,3 1256 xkvmag(kk) = xkvmag(kk) + ((xkp(ii,kk)+xkvtx(ii))/ngkpt(ii))*bmet(ii,jj)*((xkp(jj,kk)+xkvtx(jj))/ngkpt(jj)) 1257 enddo 1258 enddo 1259 enddo 1260 call indxhpsort(3,3,xkvmag,indxq) 1261 xq1=(xkp(:,1)-xkp(:,3))/ngkpt 1262 xq2=(xkp(:,2)-xkp(:,3))/ngkpt 1263 xkv0=(xkp(:,3)+xkvtx)/ngkpt 1264 xkv1=(xkp(:,1)+xkvtx)/ngkpt 1265 xkv2=(xkp(:,2)+xkvtx)/ngkpt 1266 qqdiff=xkv1-xkv2 1267 vt1=xkv0+xq2 1268 vt2=xq1-xq2 1269 call crossreduced(xq1,xq2,xnorm) 1270 area=sqrt(dot_product(xnorm,xnorm)) 1271 aa0=xkv0(iv)*xkv0(jv) 1272 aa1=xq1(iv)*xkv0(jv)+xq1(jv)*xkv0(iv) 1273 aa2=xq1(iv)*xq1(jv) 1274 bb0=xq2(iv)*xkv0(jv)+xq2(jv)*xkv0(iv) 1275 bb1=xq1(iv)*xq2(jv)+xq1(jv)*xq2(iv) 1276 cc=xq2(iv)*xq2(jv) 1277 dd0 = 0.d0 1278 dd1 = 0.d0 1279 dd2 = 0.d0 1280 ee0 = 0.d0 1281 ee1 = 0.d0 1282 ff = 0.d0 1283 d0e0f = 0.d0 1284 d1e1e02f = 0.d0 1285 d2e1f = 0.d0 1286 eta0 = 0.d0 1287 eta2 = 0.d0 1288 nu2 = 0.d0 1289 qqmag(1) = 0.d0 1290 qqmag(2) = 0.d0 1291 qqdiffmag = 0.d0 1292 do ii=1,3 1293 do jj=1,3 1294 dd0 = dd0 + xkv0(ii)*bmet(ii,jj)*xkv0(jj) 1295 dd1 = dd1 + xq1(ii)*bmet(ii,jj)*xkv0(jj)*2 1296 dd2 = dd2 + xq1(ii)*bmet(ii,jj)*xq1(jj) 1297 ee0 = ee0 + xq2(ii)*bmet(ii,jj)*xkv0(jj)*2 1298 ee1 = ee1 + xq1(ii)*bmet(ii,jj)*xq2(jj)*2 1299 ff = ff + xq2(ii)*bmet(ii,jj)*xq2(jj) 1300 d0e0f = d0e0f + vt1(ii)*bmet(ii,jj)*vt1(jj) 1301 d1e1e02f = d1e1e02f + vt1(ii)*bmet(ii,jj)*vt2(jj)*2 1302 d2e1f = d2e1f + vt2(ii)*bmet(ii,jj)*vt2(jj) 1303 eta0 = eta0 + qqdiff(ii)*bmet(ii,jj)*xkv0(jj) 1304 eta2 = eta2 + qqdiff(ii)*bmet(ii,jj)*xkv2(jj) 1305 nu2 = nu2 + qqdiff(ii)*bmet(ii,jj)*xq2(jj) 1306 qqmag(1) = qqmag(1) + xkv1(ii)*bmet(ii,jj)*xkv1(jj) 1307 qqmag(2) = qqmag(2) + xkv2(ii)*bmet(ii,jj)*xkv2(jj) 1308 qqdiffmag = qqdiffmag + qqdiff(ii)*bmet(ii,jj)*qqdiff(jj) 1309 enddo 1310 enddo 1311 xi1 = dd1/(2*sqrt(dd0*dd2)) 1312 xi2 = ee0/(2*sqrt(dd0*ff)) 1313 zeta12 = ee1/(2*sqrt(dd2*ff)) 1314 eta0 = eta0/sqrt(qqdiffmag*dd0) 1315 eta2 = eta2/sqrt(qqdiffmag*qqmag(2)) 1316 nu2 = nu2/sqrt(qqdiffmag*ff) 1317 call sinreducedangle(xkv0,xq1,omxi12) 1318 call sinreducedangle(xkv0,xq2,omxi22) 1319 call sinreducedangle(xq1,xq2,omzeta122) 1320 call sinreducedangle(xkv2,qqdiff,ometa22) 1321 q02k1 = sqrt(dd0/dd2) 1322 q22qdiff = sqrt(qqmag(2)/qqdiffmag) 1323 rrootd = -q02k1*xi1 1324 irootd = q02k1*omxi12 != q02k1*sqrt(1-xi1**2) 1325 call crossreduced(xq2,xkv0,utemp) 1326 call crossreduced(xq2,xq1,vtemp) 1327 utemp = 2*utemp 1328 vtemp = 2*vtemp 1329 call cross(utemp,vtemp,wtemp) 1330 vmag2 = dot_product(vtemp,vtemp) 1331 umag2 = dot_product(utemp,utemp) 1332 wmag2 = dot_product(wtemp,wtemp) 1333 deltaprefact = vmag2 ! = 4*dd2*ff-ee1*ee1 1334 if (vmag2.eq.0.d0) then 1335 rrootdelta = 0.d0 1336 irootdelta = 0.d0 1337 else 1338 rrootdelta = -dot_product(utemp,vtemp)/vmag2 1339 irootdelta = sqrt(wmag2)/vmag2 1340 endif 1341! rrootdelta = -q02k1*kappa 1342! irootdelta = q02k1*lambda 1343 rrootdezfzz = -q22qdiff*eta2 1344 irootdezfzz = q22qdiff*ometa22 != q22qdiff*sqrt(1-eta2**2) 1345!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "q1*q2: ",dd0+0.5*(dd1+ee0+ee1) 1346!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "eta2: ",eta2 1347!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "xi1: ",xi1 1348!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "xi2: ",xi2 1349!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "zeta12: ",zeta12 1350!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "q02k1: ",q02k1 1351!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "dd2: ",dd2 1352!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "ff: ",ff 1353!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "1-zeta12**2:",1-zeta12**2 1354!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "deltaprefact: ",deltaprefact 1355!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "umag2: ",umag2 1356!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "vmag2: ",vmag2 1357!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "wmag2: ",wmag2 1358!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "rroot (d): ",rrootd 1359!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "iroot (d): ",irootd 1360!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "rroot (delta): ",rrootdelta, & 1361!& -q02k1*(xi1-zeta12*xi2)/(1-zeta12**2) 1362!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "iroot (delta): ",irootdelta, & 1363!& q02k1*sqrt(1.d0 + 2*xi1*xi2*zeta12 - xi1*xi1 - xi2*xi2 - zeta12*zeta12)/(1-zeta12**2) 1364!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "rroot (d+ez+fz^2): ",rrootdezfzz 1365!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "iroot (d+ez+fz^2): ",irootdezfzz 1366!write(77,*) xx, aa, aa2*((xx+q02k1*xi1)**2 + (1-xi1**2)*(q02k1)**2) 1367! brd=max(abs(dd0),abs(dd1),abs(dd2),abs(ee0),abs(ee1),abs(ff))*1.e-10 1368 nsing=0 1369 root1 = rrootd-100*irootd 1370 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1371 nsing=nsing+1 1372 xsing(nsing)=root1 1373!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "D: ",xsing(nsing) 1374 endif 1375 root1 = rrootd-10*irootd 1376 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1377 nsing=nsing+1 1378 xsing(nsing)=root1 1379!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "D: ",xsing(nsing) 1380 endif 1381 root1 = rrootd 1382 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1383 nsing=nsing+1 1384 xsing(nsing)=root1 1385!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "D: ",xsing(nsing) 1386 endif 1387 root1 = rrootd+10*irootd 1388 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1389 nsing=nsing+1 1390 xsing(nsing)=root1 1391!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "D: ",xsing(nsing) 1392 endif 1393 root1 = rrootd+100*irootd 1394 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1395 nsing=nsing+1 1396 xsing(nsing)=root1 1397!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "D: ",xsing(nsing) 1398 endif 1399 root1 = rrootdelta-100*irootdelta 1400 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1401 nsing=nsing+1 1402 xsing(nsing)=root1 1403!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "Delta: ",xsing(nsing) 1404 endif 1405 root1 = rrootdelta-10*irootdelta 1406 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1407 nsing=nsing+1 1408 xsing(nsing)=root1 1409!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "Delta: ",xsing(nsing) 1410 endif 1411 root1 = rrootdelta 1412 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1413 nsing=nsing+1 1414 xsing(nsing)=root1 1415!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "Delta: ",xsing(nsing) 1416 endif 1417 root1 = rrootdelta+10*irootdelta 1418 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1419 nsing=nsing+1 1420 xsing(nsing)=root1 1421!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "Delta: ",xsing(nsing) 1422 endif 1423 root1 = rrootdelta+100*irootdelta 1424 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1425 nsing=nsing+1 1426 xsing(nsing)=root1 1427!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "Delta: ",xsing(nsing) 1428 endif 1429 root1 = rrootdezfzz-100*irootdezfzz 1430 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1431 nsing=nsing+1 1432 xsing(nsing)=root1 1433!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "D+EZ+FZ^2: ",xsing(nsing) 1434 endif 1435 root1 = rrootdezfzz-10*irootdezfzz 1436 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1437 nsing=nsing+1 1438 xsing(nsing)=root1 1439!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "D+EZ+FZ^2: ",xsing(nsing) 1440 endif 1441 root1 = rrootdezfzz 1442 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1443 nsing=nsing+1 1444 xsing(nsing)=root1 1445!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "D+EZ+FZ^2: ",xsing(nsing) 1446 endif 1447 root1 = rrootdezfzz+10*irootdezfzz 1448 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1449 nsing=nsing+1 1450 xsing(nsing)=root1 1451!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "D+EZ+FZ^2: ",xsing(nsing) 1452 endif 1453 root1 = rrootdezfzz+100*irootdezfzz 1454 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1455 nsing=nsing+1 1456 xsing(nsing)=root1 1457!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "D+EZ+FZ^2: ",xsing(nsing) 1458 endif 1459 call hpsort(nsing,nsing,xsing(1:nsing)) 1460 xmin=0.d0 1461 xmax=1.d0 1462 ii = 1 1463 if (nsing.gt.0) then 1464 do 1465 root1 = xsing(nsing)*100 1466 if (root1.ge.xmax) then 1467 exit 1468 else 1469 nsing = nsing+1 1470 xsing(nsing)=root1 1471 endif 1472 enddo 1473 endif 1474 root1 = -ee0/ee1 1475 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1476 nsing=nsing+1 1477 xsing(nsing)=root1 1478!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "E: ",xsing(nsing) 1479 endif 1480 root1 = -(ee0+2*ff)/(ee1-2*ff) 1481 if (root1.gt.0.d0.and.root1.lt.1.d0) then 1482 nsing=nsing+1 1483 xsing(nsing)=root1 1484!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "E+2FZ: ",xsing(nsing) 1485 endif 1486 do 1487 if (ii.ge.nsing) exit 1488 if (xsing(ii).eq.xsing(ii+1)) then 1489 nsing = nsing-1 1490 do jj=ii+1,nsing 1491 xsing(jj) = xsing(jj+1) 1492 enddo 1493 endif 1494 ii = ii + 1 1495 enddo 1496 call hpsort(nsing,nsing,xsing(1:nsing)) 1497 1498 gammapfpf = 4*(cc*dd2+ff*aa2)-2*bb1*ee1 1499 call rquadroots(4*(cc*dd2+ff*aa2)-2*bb1*ee1, & 1500& 4*(cc*dd1+ff*aa1)-2*(bb0*ee1+bb1*ee0), & 1501& 4*(cc*dd0+ff*aa0)-2*bb0*ee0, & 1502& gammapfnroot,gammapfroot1,gammapfroot2) 1503 1504 polyterm(4) = 2*ff*aa2*dd2 - aa2*ee1**2 + ff*aa2*ee1 + dd2*bb1*ee1 & 1505& - 2*ff*dd2*bb1 + cc*dd2*ee1 - 2*cc*dd2*dd2 1506 polyterm(3) = 2*ff*(aa1*dd2+aa2*dd1) & 1507& - (aa1*ee1**2 + 2*aa2*ee0*ee1) & 1508& - ff*(aa2*(ee1-ee0)-aa1*ee1) & 1509& + (dd2*(bb0*ee1+bb1*ee0)+dd1*bb1*ee1) & 1510& + 2*ff*(dd2*(bb1-bb0)-dd1*bb1) & 1511& - cc*(dd2*(ee1-ee0)-dd1*ee1) & 1512& - 4*cc*dd1*dd2 1513 polyterm(2) = 2*ff*(aa2*dd0+aa1*dd1+aa0*dd2) & 1514& - (aa0*ee1**2+2*aa1*ee0*ee1+aa2*ee0**2) & 1515& - ff*(-aa0*ee1+aa1*(ee1-ee0)+aa2*ee0) & 1516& + (dd0*bb1*ee1+dd1*(bb0*ee1+bb1*ee0)+dd2*bb0*ee0) & 1517& + 2*ff*(-dd0*bb1+dd1*(bb1-bb0)+dd2*bb0) & 1518& - cc*(-dd0*ee1+dd1*(ee1-ee0)+dd2*ee0) & 1519& - 2*cc*(2*dd0*dd2+dd1**2) 1520 polyterm(1) = 2*ff*(aa0*dd1+aa1*dd0) & 1521& - (2*aa0*ee0*ee1+aa1*ee0**2) & 1522& - ff*(aa0*(ee1-ee0)+aa1*ee0) & 1523& + (dd0*(bb0*ee1+bb1*ee0)+dd1*bb0*ee0) & 1524& + 2*ff*(dd0*(bb1-bb0)+dd1*bb0) & 1525& - cc*(dd0*(ee1-ee0)+dd1*ee0) & 1526& - 4*cc*(dd0*dd1) 1527 polyterm(0) = 2*ff*aa0*dd0 - aa0*ee0**2 - ff*aa0*ee0 + dd0*bb0*ee0 & 1528& + 2*ff*dd0*bb0 - cc*dd0*ee0 - 2*cc*dd0**2 1529 denompf = polyterm(4) 1530!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "AAA" 1531 call rquarticroots(polyterm(4),polyterm(3),polyterm(2),polyterm(1),polyterm(0), & 1532& denomnroot,denomroot1,denomroot2,denomroot3,denomroot4) 1533!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "BBB" 1534 1535 abr=1.d-12 1536 rlr=1.d-6 1537if (k_ct_DEBUG.eq.k_ref_DEBUG) then 1538write(6,*) ww,i_ct_DEBUG 1539if (ww.ge.0.002125.and.ww.le.0.0021271) then 1540write(6,*) "nsing: ",nsing 1541do ii=1,nsing 1542write(6,*) xsing(ii) 1543enddo 1544write(6,*) "|q0|: ",sqrt(dd0) 1545write(6,*) "d prefact: ",dd2 1546write(6,*) "rroot (d): ",rrootd 1547write(6,*) "iroot (d): ",irootd 1548write(6,*) "delta prefact: ",deltaprefact 1549write(6,*) "rroot (delta): ",rrootdelta 1550write(6,*) "iroot (delta): ",irootdelta 1551write(6,*) "(d+ez+fz^2) prfct: ",qqdiffmag 1552write(6,*) "rroot (d+ez+fz^2): ",rrootdezfzz 1553write(6,*) "iroot (d+ez+fz^2): ",irootdezfzz 1554write(6,*) "x^4 term: ",polyterm(4),denompf 1555write(6,*) "x^3 term: ",polyterm(3), & 1556& -denompf*(denomroot1+denomroot2+2*denomroot3) 1557write(6,*) "x^2 term: ",polyterm(2), & 1558& denompf*(denomroot1*denomroot2+2*denomroot1*denomroot3 & 1559& +2*denomroot2*denomroot3+denomroot3*denomroot3+denomroot4*denomroot4) 1560write(6,*) "x^1 term: ",polyterm(1), & 1561& denompf*((denomroot1+denomroot2)*(denomroot3*denomroot3+denomroot4*denomroot4) & 1562& +2*denomroot1*denomroot2*denomroot3) 1563write(6,*) "x^0 term: ",polyterm(0), denompf*(denomroot1*denomroot2*(denomroot3*denomroot3+denomroot4*denomroot4)) 1564write(6,*) "gammapfpf: ",gammapfpf 1565write(6,*) "gammapfnroot: ",gammapfnroot 1566write(6,*) "roots: " 1567write(6,*) gammapfroot1 1568write(6,*) gammapfroot2 1569write(6,*) "denompf: ",denompf 1570write(6,*) "denomnroot: ",denomnroot 1571write(6,*) "roots: " 1572write(6,*) denomroot1 1573write(6,*) denomroot2 1574write(6,*) denomroot3 1575write(6,*) denomroot4 1576write(6,*) 1577l_ct_debug = l_ct_DEBUG + 1 1578i_prt_DEBUG = 1 1579else 1580i_prt_DEBUG = 0 1581endif 1582if (l_ct_DEBUG.ge.100) stop 1583!write(6,*) 1584!write(6,*) "q0^2: ",xkvmag(indxq(1)) 1585!write(6,*) "q1^2: ",xkvmag(indxq(2)) 1586!write(6,*) "q2^2: ",xkvmag(indxq(3)) 1587!write(6,*) "|q0|: ",sqrt(dd0) 1588!write(6,*) "|k1|: ",sqrt(dd2) 1589!write(6,*) "|k2|: ",sqrt(ff) 1590!write(6,*) "q0: ",xkv0," (reduced)" 1591!write(6,*) "q1: ",(xkp(:,indxq(2))+xkvtx)/ngkpt," (reduced)" 1592!write(6,*) "q2: ",(xkp(:,indxq(3))+xkvtx)/ngkpt," (reduced)" 1593!write(6,*) "k1: ",xq1," (reduced)" 1594!write(6,*) "k2: ",xq2," (reduced)" 1595!write(6,*) "d: ",dd0,dd1,dd2 1596!write(6,*) "e: ",ee0,ee1 1597!write(6,*) "f: ",ff 1598!write(6,*) "d0+e0+f: ",dd0+ee0+ff,d0e0f 1599!write(6,*) "d1+e1-e0-2f: ",dd1+ee1-ee0-2*ff,d1e1e02f 1600!write(6,*) "d2-e1+f: ",dd2-ee1+ff,d2e1f 1601endif 1602 sint=grater(tfn1,xmin,xmax,abr,rlr,nsing,xsing,error,numcal,maxns)*area 1603!!write(6,*) "%%%%%%%% ",ww,sint,i_ct_DEBUG 1604if (k_ct_DEBUG.eq.k_ref_DEBUG) then 1605write(78,*) ww,sint,i_ct_DEBUG 1606if (i_ct_DEBUG.gt.1000) stop 1607if (j_ct_DEBUG+1.gt.1) stop 1608endif 1609!if (k_ct_DEBUG.gt.k_ref_DEBUG) stop 1610end subroutine tsingint 1611 1612!*************************************************************************** 1613 1614double precision function tfn1(xx) 1615implicit none 1616double precision :: xx 1617double precision :: aa0,aa1,aa2,bb0,bb1,cc,dd0,dd1,dd2,ee0,ee1,ff,brd 1618double precision :: xkv0(3),xq1(3) 1619double precision :: aa,bb,dd,ee 1620double precision :: xx2,zz,delta,sqdelta,xterm0,xterm1,xterm2,dterm1,dterm2, & 1621& gammat,dgammaD,dgammaD_old 1622double precision :: xx0d,dxxd 1623double precision :: xterm1inv 1624double precision :: logterm1,logterm2,logterm3,logterm4,atf,darg,dargoe 1625double precision :: arg1,arg2,sum1,sum2,prod1,temp 1626double precision :: test1,test2,test3,test4,testlo1,testlo2,testlo3 1627double precision :: dtarg 1628integer :: xi 1629integer :: ii,jj,kk,ll 1630integer :: ivv,jvv 1631logical :: lowdarg 1632double precision :: xdiff1,xdiff2,dquad 1633double precision :: h1,h2,h3,h4,h5,h6,h7,h8,term1,term2,term3,term4 1634double precision :: aa2dd 1635double precision :: hiorderexp,quadterm1,quadterm2,quadterm3 1636double precision :: tt0,tt1,tt2,ta0,ta1,ta2,tx0,tx1,tx2,tc1,tc2 1637double precision :: xd01,xd02,xd0,dx ! xd0 -> xx value where delta is 0 1638double precision :: xdt11,xdt12,xdt21,xdt22,xdt1,xdt2 ! for xx = xd0 + dx, linear and quadratic term in dx for delta 1639double precision :: xx_dd1, xx_dd2 ! roots of dd = 0 1640double precision :: xx_def1, xx_def2 ! roots of dd + ee*zz + ff*zz^2 = 0 1641logical :: usexd1,usexd2,usexx_dd1,usexx_dd2 ! are roots in integration region 1642double precision :: aat1,aat2,bbt1,bbt2,ddt1,ddt2,eet1,eet2,zzt1,zzt2,deltat1,deltat2,aa2ddt1,aa2ddt2 1643double precision :: RRA1,RRA2,RRB1,RRB2,RRC1,RRC2,d0e0f,d1e1e02f,d2e1f 1644double precision :: qqmag(2), q02k1, q22qdiff, qqdiffmag 1645double precision :: xi1,xi2,zeta12,eta0,eta2,nu2,kappa,lambda 1646double precision :: omxi12,omxi22,ometa22,omzeta122 1647double precision :: deltaprefact, umag2 1648double precision :: rrootd,irootd,rrootdelta,irootdelta,rrootdezfzz,irootdezfzz 1649double precision :: gammapfroot1,gammapfroot2,gammapfpf 1650integer :: gammapfnroot 1651double precision :: denomroot1,denomroot2,denomroot3,denomroot4,denompf 1652integer :: denomnroot 1653common /tfn/ aa0,aa1,aa2,bb0,bb1,cc,dd0,dd1,dd2,ee0,ee1,ff,brd, & 1654& xkv0,xq1,ivv,jvv, & 1655& xd01,xd02,xdt11,xdt12,xdt21,xdt22,xx_dd1,xx_dd2,xx_def1,xx_def2, & 1656& usexd1,usexd2,usexx_dd1,usexx_dd2, & 1657& aat1,aat2,bbt1,bbt2,ddt1,ddt2,eet1,eet2,zzt1,zzt2, & 1658& deltat1,deltat2,aa2ddt1,aa2ddt2, qqmag, qqdiffmag, & 1659& q02k1, q22qdiff, deltaprefact,umag2, & 1660& rrootd,irootd,rrootdelta,irootdelta,rrootdezfzz,irootdezfzz, & 1661& gammapfroot1,gammapfroot2,denomroot1,denomroot2,denomroot3,denomroot4, & 1662& gammapfnroot,denomnroot, & 1663& gammapfpf,denompf, & 1664& RRA1,RRA2,RRB1,RRB2,RRC1,RRC2,d0e0f, & 1665& xi1,xi2,zeta12,eta0,eta2,kappa,lambda,omxi12,omxi22,ometa22,omzeta122 1666double precision :: xdummy, ydummy, zdummy,xxdummy,yydummy,zzdummy 1667integer :: idummy,jdummy,kdummy,ldummy 1668integer :: i_prt_DEBUG, j_prt_DEBUG, k_prt_DEBUG, & 1669& i_ct_DEBUG, j_ct_DEBUG, k_ct_DEBUG, l_ct_DEBUG, & 1670& k_ref_DEBUG, i_ref_DEBUG 1671common /tsingint_DEBUG/ i_prt_DEBUG, j_prt_DEBUG, k_prt_DEBUG, & 1672& i_ct_DEBUG, j_ct_DEBUG, k_ct_DEBUG, l_ct_DEBUG, & 1673& k_ref_DEBUG, i_ref_DEBUG 1674 1675xdummy = 0.d0 1676ydummy = 0.d0 1677xxdummy = 0.d0 1678yydummy = 0.d0 1679 xx2 = xx*xx 1680 zz = 1.d0-xx 1681 xx0d = -q02k1*xi1 1682 dxxd = xx - xx0d 1683 aa = (xkv0(ivv)+xx*xq1(ivv))*(xkv0(jvv)+xx*xq1(jvv)) ! =aa0 + aa1*xx + aa2*xx2 1684 bb = bb0 + bb1*xx 1685 dd = dd2*((xx-rrootd)**2 + (irootd)**2) ! =dd0 + dd1*xx + dd2*xx2 1686 ee = ee0 + ee1*xx 1687 aa2dd = aa/dd 1688 if (dd.eq.0.d0) aa2dd = (aa1+(2*xx0d+dxxd)*aa2)/(dd1+(2*xx0d+dxxd)*dd2) 1689 xterm0 = zz*(ff*zz + ee) 1690 xterm1 = qqdiffmag*((xx-rrootdezfzz)**2 + (irootdezfzz)**2) ! = dd + xterm0 1691 xterm2 = xterm0/xterm1 1692 darg = 2*ff*zz 1693 dterm1 = ee+darg 1694 lowdarg = (abs(darg/ee).lt.1.d-4) 1695 if (deltaprefact.eq.0.d0) then 1696 delta = umag2 1697 else 1698 delta = deltaprefact*((xx-rrootdelta)**2 + (irootdelta)**2) ! =4*dd*ff-ee*ee 1699 endif 1700xdummy = 4*dd*ff-ee*ee 1701ydummy = dd0 + dd1*xx + dd2*xx2 1702zdummy = ydummy + xterm0 1703xxdummy = (atan((ee+darg)/sqrt(xdummy)) - atan(ee/sqrt(xdummy)))/sqrt(xdummy) 1704yydummy = aa0 + aa1*xx + aa2*xx2 1705 if (lowdarg) then 1706 dargoe = darg/ee 1707 test1 = (dargoe+dargoe**2+(dargoe**3)/3.d0)/((ee+darg)**3) 1708 test2 = delta*(dargoe+2*dargoe**2+2*dargoe**3+dargoe**4+(dargoe**5)/5.d0)/((ee+darg)**5) 1709 else 1710 test1 = (1.d0/ee**3 - 1.d0/(ee+darg)**3)/3.d0 1711 test2 = delta*(1.d0/ee**5 - 1.d0/(ee+darg)**5)/5.d0 1712 endif 1713 if (abs(test2/test1).lt.1.d-4) then 1714! need series expansion to avoid numerical instability due to differences in large numbers 1715! dgammaD = (gamma - (1/ee - 1/(ee+darg)))/delta = (gamma - (1/ee)(darg/(ee+darg)))/delta 1716idummy=0 1717 dgammaD = test2-test1 1718 ii = 2 1719 xi = -1 1720 do 1721 jj = 2*ii+3 1722 if (lowdarg) then 1723 ll = 1 ! binomial coefficient 1724 sum1 = 0.d0 1725 prod1 = 1.d0 1726 do kk=1,jj ! binomial series expansion for ((ee+darg)**jj - ee**jj)/ee**jj 1727 ll = (ll*(jj-kk+1))/(kk) 1728 prod1 = prod1*dargoe 1729 sum1 = sum1 + ll*prod1 1730 enddo 1731 test3 = delta**ii * (sum1/(dble(jj)*(ee+darg)**jj))/dble(jj) 1732 else 1733 test3 = delta**ii * (1.d0/ee**jj - 1.d0/(ee+darg)**jj)/dble(jj) 1734 endif 1735 dgammaD_old = dgammaD 1736 dgammaD = dgammaD + xi*test3 1737 if (dgammaD.eq.dgammaD_old) exit 1738 xi = -1*xi 1739 ii = ii+1 1740 enddo 1741 temp = zz/(ee*dterm1*xterm1) 1742 tt0 = (aa/dd)*((2*dd*ff-dterm1*(ff*zz+ee))*temp + 4*dd*ff*dgammaD) 1743 tt1 = -bb*(ee*temp + 2*ee*dgammaD) 1744 tt2 = cc*((2*dd + dterm1*zz)*temp + 4*dd*ff*dgammaD) 1745 tfn1 = tt0 + tt1 + tt2 1746 else 1747 if (delta.lt.0.d0) then 1748idummy=1 1749 write(6,*) "ERROR in function tfn1: negative delta is unphysical" 1750 stop 1751 else 1752idummy=2 1753 sqdelta = sqrt(delta) 1754 dtarg = darg/sqdelta 1755 arg1 = ee/sqdelta 1756 arg2 = arg1 + dtarg 1757 testlo1 = zz/(2*dd) ! = darg/(delta+ee*ee) 1758 testlo2 = -zz*zz*ee/(4*dd*dd) ! = darg**2*ee/(delta+ee*ee)**2 1759 testlo3 = zz**3 * (ee*ee - delta/3.d0) / (8*dd**3) 1760 test3 = darg/(ee*(ee+darg)) 1761 test4 = ((ee*ee*darg+ee*darg*darg+darg**3/3.d0)/(ee*(ee+darg))**3)*delta 1762 ta0 = 2*dd*ff-ee*ee-ee*ff*zz 1763 ta1 = ee+2*ff*zz 1764 ta2 = ee*zz+2*dd 1765 if ((abs(testlo1).gt.abs(1.d3*testlo2)).and.(abs(testlo1).gt.abs(1.d6*testlo3))) then 1766! if difference in atan arguments small, may lead to numerical error 1767! atan(z+d)->atan(z)+d/(1+z^2)-d^2z/(1+z^2)^2+d^3(z^2-1/3)/(1+z^2)^3+..., d->0 (taylor expansion) 1768 call atan_expansion(darg/sqdelta,arg1,hiorderexp,3) 1769 hiorderexp = hiorderexp/sqdelta 1770 gammat = testlo1+testlo2+testlo3 + hiorderexp 1771jdummy=1 1772 tt0 = aa2dd*(zz*ta0/xterm1 + 4*ff*dd*gammat)/delta 1773 tt1 = bb*zz*zz*(dd*delta-ee*zz*(2*dd*ff-ee*ee-ee*ff*zz)) & 1774& /(2*dd*dd*delta*xterm1) & 1775& - 2*ee*bb*hiorderexp/delta 1776 tt2 = cc*(zz**3*(2*ff*dd - ff*zz*ee - ee*ee))/(dd*delta*xterm1) & 1777& + 4*dd*hiorderexp*cc/delta 1778 tfn1 = tt0 + tt1 + tt2 1779! atf = gammat*sqdelta 1780 elseif (abs(test3).gt.abs(1.d6*test4)) then 1781! if both atan arguments large, may lead to numerical error 1782! atan(z)->pi/2-1/z+1/(3z^3)-1/(5z^5)+..., z->infty 1783! atan(z)->-pi/2-1/z+1/(3z^3)-1/(5z^5)+..., z->-infty 1784 gammat = test3-test4 1785! atf = gammat*sqdelta 1786 tt0 = aa2dd*(zz*ta0/xterm1 + 4*ff*dd*gammat)/delta 1787 tt1 = -bb*(-zz*ta1/xterm1 + 2*ee*gammat)/delta 1788 tt2 = cc*(-zz*ta2/xterm1 + 4*dd*gammat)/delta 1789 tfn1 = tt0 + tt1 + tt2 1790jdummy=2 1791 else 1792 atf = atan(arg2) - atan(arg1) 1793 gammat = atf/sqdelta 1794 if (gammapfnroot.gt.0) then 1795 quadterm3 = gammapfpf*(xx-gammapfroot1)*(xx-gammapfroot2) 1796 else 1797 quadterm3 = gammapfpf*((xx-gammapfroot1)**2 + gammapfroot2**2) 1798 endif 1799 if (denomnroot.ge.4) then 1800 quadterm1 = (xx-denomroot3)*(xx-denomroot4) 1801 else 1802 quadterm1 = (xx-denomroot3)**2+denomroot4**2 1803 endif 1804 if (denomnroot.ge.2) then 1805 quadterm2 = ((xx-denomroot1)*(xx-denomroot2))*denompf 1806 else 1807 quadterm2 = ((xx-denomroot1)**2+denomroot2**2)*denompf 1808 endif 1809 tc1 = zz*quadterm1*quadterm2/(dd*xterm1)/delta 1810 tc2 = quadterm3*gammat/delta 1811 test1 = tc1+tc2 1812 tt0 = aa2dd*(zz*ta0/xterm1 + 4*ff*dd*gammat)/delta 1813 tt1 = -bb*(-zz*ta1/xterm1 + 2*ee*gammat)/delta 1814 tt2 = cc*(-zz*ta2/xterm1 + 4*dd*gammat)/delta 1815 test2 = tt0 + tt1 + tt2 1816! test1 and test2 are mathematically equivalent, choose the representation that minimizes error due to difference of similar large numbers 1817 if (max(abs(tc1),abs(tc2)).lt.max(abs(tt0),abs(tt1),abs(tt2))) then 1818 tfn1 = test1 1819 else 1820 tfn1 = test2 1821 endif 1822jdummy=0 1823 endif 1824 endif 1825 endif 1826 1827!if (i_prt_DEBUG.gt.0) then 1828!write(79+l_ct_DEBUG,*) xx,gammat,atf,atan(arg2),atan(arg1),sqdelta 1829!endif 1830if (i_ct_DEBUG.eq.i_ref_DEBUG) then 1831j_ct_DEBUG = j_ct_DEBUG+1 1832!write(77,*) xx,tfn1,idummy,jdummy,i_ct_DEBUG 1833!write(77,*) xx,tfn1,tt0,tt1,tt2 1834!write(77,*) xx, & 1835!& (aa2dd/delta)*(zz*ta0/xterm1 + 4*ff*dd*gammat), & 1836!& (-bb/delta)*(-zz*ta1/xterm1 + 2*ee*gammat), & 1837!& (cc/delta)*(-zz*ta2/xterm1 + 4*dd*gammat), & 1838!& (yydummy/(ydummy*xdummy))*(zz*(2*ydummy*ff-ee*ee-ee*ff*zz)/zdummy+4*ydummy*ff*xxdummy), & 1839!& -(bb/delta)*(-zz*(ee+2*ff*zz)/zdummy+2*ee*xxdummy), & 1840!& (cc/delta)*(-zz*(ee*zz+2*ydummy)/zdummy+4*ydummy*xxdummy) 1841!write(77,*) xx,tfn1,idummy,jdummy,(-zz*ta2/xterm1 + 4*dd*gammat),(zz**3*(2*ff*dd - ff*zz*ee - ee*ee))/(dd*xterm1) 1842!write(77,*) xx,tfn1,tt2,-zz*ta2/xterm1,4*dd*gammat 1843!write(77,*) xx,tfn1,delta,4*ff*(sqrt(dd0)*omxi22 + xx*sqrt(dd2)*omzeta122)**2 1844!if (jdummy.eq.1.and.idummy.eq.2) write(77,*) xx,tfn1,test1,test2,darg,ee,delta 1845! 1846!if (j_ct_DEBUG+1.gt.1000) return 1847if (j_ct_DEBUG+1.gt.1000) then 1848write(6,*) "Programmed DEBUG stop" 1849if (idummy.eq.0) then 1850 write(6,*) "Gamma from series expansion in Delta, idummy = ", idummy 1851else if (idummy.eq.1) then 1852 write(6,*) "Gamma from natural logarithm function, idummy = ", idummy 1853else 1854 write(6,*) "Gamma from arctangent function, idummy = ", idummy 1855endif 1856stop 1857endif 1858endif 1859 1860if (.not.(tfn1.eq.tfn1)) then 1861write(6,*) "NAN error in tfn1" 1862write(6,*) "xx = ",xx 1863write(6,*) "zz = ",zz 1864if (idummy.eq.0) then 1865 write(6,*) "Gamma from series expansion in Delta" 1866else if (idummy.eq.1) then 1867 write(6,*) "Gamma from natural logarithm function" 1868else 1869 write(6,*) "Gamma from arctangent function" 1870endif 1871write(6,*) idummy,jdummy,kdummy 1872write(6,*) xdummy,ydummy,zdummy 1873write(6,*) xxdummy,yydummy,zzdummy 1874write(6,*) "ta0 = ",ta0 1875write(6,*) "ta1 = ",ta1 1876write(6,*) "ta2 = ",ta2 1877write(6,*) "tt0 = ",tt0 1878write(6,*) "tt1 = ",tt1 1879write(6,*) "tt2 = ",tt2 1880write(6,*) "aa = ",aa 1881write(6,*) "bb = ",bb 1882write(6,*) "cc = ",cc 1883write(6,*) "dd = ",dd 1884write(6,*) "ee = ",ee 1885write(6,*) "ff = ",ff 1886write(6,*) "dd0 = ",dd0 1887write(6,*) "dd1 = ",dd1 1888write(6,*) "dd2 = ",dd2 1889write(6,*) "ee0 = ",ee0 1890write(6,*) "ee1 = ",ee1 1891write(6,*) "rrootd = ",rrootd 1892write(6,*) "irootd = ",irootd 1893write(6,*) "delta = ",delta 1894write(6,*) "deltaprefact = ",deltaprefact 1895write(6,*) "q02k1 = ",q02k1 1896write(6,*) "kappa = ",kappa 1897write(6,*) "lambda = ",lambda 1898write(6,*) "darg = ",darg 1899write(6,*) "xterm1 = ",xterm1 1900write(6,*) "gammat = ",gammat 1901write(6,*) "atf = ",atf 1902write(6,*) "sqdelta = ",sqdelta 1903stop 1904endif 1905 1906end function tfn1 1907 1908!*************************************************************************** 1909! 1910! expand (1+h1*x)(1+h2*x)(1+h3*x)(1+h4*x) = 1 + x*term1 + x^2*term2 + x^3*term3 + x^4*term4 1911 1912subroutine expand4(h1,h2,h3,h4,term1,term2,term3,term4) 1913double precision :: h1,h2,h3,h4,term1,term2,term3,term4 1914term1 = h1+h2+h3+h4 1915term2 = h1*h2 + h1*h3 + h1*h4 + h2*h3 + h2*h4 + h3*h4 1916term3 = h2*h3*h4 + h1*h3*h4 + h1*h2*h4 + h1*h2*h3 1917term4 = h1*h2*h3*h4 1918return 1919end subroutine expand4 1920 1921!*************************************************************************** 1922! 1923! sine of angle between two vectors expressed in reduced coordinates 1924 1925subroutine sinreducedangle(v1,v2,ss) 1926use geometry 1927double precision :: v1(3),v2(3),ss 1928double precision :: w1(3),w2(3),w3(3),w1mag,w2mag,w3mag 1929integer :: ii 1930do ii=1,3 1931 w1(ii) = dot_product(v1,blat(:,ii)) 1932 w2(ii) = dot_product(v2,blat(:,ii)) 1933enddo 1934call cross(w1,w2,w3) 1935w1mag = sqrt(dot_product(w1,w1)) 1936w2mag = sqrt(dot_product(w2,w2)) 1937w3mag = sqrt(dot_product(w3,w3)) 1938ss = w3mag/(w1mag*w2mag) 1939return 1940end subroutine sinreducedangle 1941 1942!*************************************************************************** 1943! 1944! cross product (returned in cartesian coordinates) of two vectors expressed in reduced coordinates 1945 1946subroutine crossreduced(v1,v2,w3) 1947use geometry 1948double precision :: v1(3),v2(3),ss 1949double precision :: w1(3),w2(3),w3(3) 1950integer :: ii 1951do ii=1,3 1952 w1(ii) = dot_product(v1,blat(:,ii)) 1953 w2(ii) = dot_product(v2,blat(:,ii)) 1954enddo 1955call cross(w1,w2,w3) 1956return 1957end subroutine crossreduced 1958 1959!*************************************************************************** 1960! atan(zz+dz) - atan(zz) 1961! using taylor expansion of atan(zz+dz) in powers of dz 1962! startterm > 1 ignores the first startterm-1 terms in the series 1963! only recommended for dz << zz 1964 1965subroutine atan_expansion(dz,zz,atex,startterm) 1966 implicit none 1967 double precision :: dz,zz,atex 1968 integer :: startterm 1969 integer :: nterm 1970 parameter (nterm = 50) 1971 double precision :: termprefact(0:nterm) 1972 integer :: ii,jj,kk,ll,mm 1973 integer :: iparity,jparity 1974 double precision :: sum1, sum1_old, sum2, term 1975 double precision :: zzpow, dzpow, denompow, denomterm 1976 1977! first term is 1/(1+zz*zz) 1978! build up subsequent terms by induction 1979! (d/d(zz)) zz^(ii-mm)/(1+zz*zz)^ii = (-(ii+mm)*zz^(ii-mm+1) + (ii-mm)*zz^(ii-mm-1))/(1+zz*zz)^(ii+1) 1980! allows power series in zz of numerator. 1981! Coefficients of zz^jj stored in termprefact(jj) 1982! Because any term has only odd or even powers of zz, can interlace coefficients from current and next terms in one array (and can forget previous terms) 1983! Might as well also include factorial part of Taylor expansion in termprefact (division by ii+1) 1984 termprefact(0) = 1.d0; 1985 sum1 = 0.d0 1986 dzpow = 1.d0 1987 denomterm = 1 + zz*zz 1988 denompow = 1.d0 1989 do ii=1,nterm 1990 iparity = mod(ii-1,2)+1 1991 jparity = mod(ii,2)+1 1992!write(6,*) ii,iparity,jparity 1993 dzpow = dzpow * dz 1994 denompow = denompow * denomterm 1995 if (iparity.eq.1) then 1996 zzpow = 1.d0 1997 else 1998 zzpow = zz 1999 endif 2000 sum2 = 0.d0 2001 do jj = jparity-1,ii,2 2002 termprefact(jj) = 0.d0 2003 enddo 2004 do jj = iparity-1,ii-1,2 2005 mm = ii-jj 2006 termprefact(jj+1) = termprefact(jj+1) + termprefact(jj)*(-(ii+mm))/dble(ii+1) 2007 if (jj.gt.0) then 2008 termprefact(jj-1) = termprefact(jj-1) + termprefact(jj)*(ii-mm)/dble(ii+1) 2009 endif 2010!write(6,'(1x,3i3,4x,2i3,f10.3,i3,f10.3,f10.3)') ii,jj,mm,ii+1,jj+1, & 2011!& termprefact(jj),-(ii+mm),(-(ii+mm))/dble(ii+1), & 2012!& termprefact(jj+1) 2013!if (jj.gt.0) then 2014!write(6,'(14x,2i3,f10.3,i3,f10.3,f10.3)') ii+1,jj-1, & 2015!& termprefact(jj),(ii-mm),(ii-mm)/dble(ii+1), & 2016!& termprefact(jj-1) 2017!else 2018!write(6,*) 2019!endif 2020 sum2 = sum2 + termprefact(jj)*zzpow 2021 zzpow = zzpow * zz * zz 2022 enddo 2023 if (ii.ge.startterm) then 2024 term = dzpow * sum2 / (denompow) 2025 sum1_old = sum1 2026 sum1 = sum1 + term 2027!write(6,*) "**** ",ii," ",term,sum1 2028!read(*,*) 2029 if (sum1_old.eq.sum1) exit 2030 endif 2031 enddo 2032 atex = sum1 2033 2034end subroutine atan_expansion 2035 2036 2037 2038 2039 2040 2041 2042 2043 2044 2045