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