1subroutine mk2cse(iband,ikk,lossfn, &
2& vol,pi,nwpt,wmax,nbcore,nbocc,ncband,ngkpt,natom,xred,projwf,nlmn, &
3& pwmatel,tpwmatel, &
4& kg,kgq,enrgy,enrgyq,cg,cgq,npwt,npwtq,bantot,bantotq,ncg,ncgq, &
5& indxkpw,indxkpwq,indxkbnd,indxkbndq,indxkcg,indxkcgq,npwarr,npwarrq, &
6& kpt,kptq,nkpt,nkptq,nqpt,nsymk,nsymkq,symk,symkq,nsym,nsymq,symrel,syminv, &
7& ihlf,ihlfq,lvtrans,lvtransq,bmet,blat,ipaw, &
8& ipwx,ipwndx,npwndx,ntpwndx,napwndx, &
9& npwc,npwx,invpw2ndx,pwsymndx,iqsymndx, &
10& igmx,igmn,igndx,igndxq,ikndx,ikndxq,iqndx,isymndx,isymndxq,npw,npwq, &
11& nband,nbandq,nsppol,shiftk,shiftkq,zz,cse,xse)
12implicit none
13integer :: iband,ikk(3),nwpt,nbcore,nbocc,ncband,ngkpt(3),natom,nlmn
14integer :: igmn(3),igmx(3)
15integer :: npwt,npwtq,bantot,bantotq,ncg,ncgq,nkpt,nkptq,nqpt,nsym,nsymq,nsppol
16integer :: npw,npwc,npwx,npwq,ipw1,npwndx,ntpwndx,napwndx,ipaw
17double precision :: vol,pi,wmax,xred(3,natom)
18double complex :: lossfn(nwpt,nqpt,ntpwndx)
19double complex :: projwf(natom,nlmn,nkpt,ncband)
20integer :: kg(3,npwt),kgq(3,npwtq)
21double precision :: enrgy(bantot),enrgyq(bantotq)
22double complex :: cg(ncg),cgq(ncgq)
23integer :: indxkpw(nkpt),indxkpwq(nkptq),indxkbnd(nkpt),indxkbndq(nkptq)
24integer :: indxkcg(nkpt),indxkcgq(nkptq),npwarr(nkpt),npwarrq(nkptq)
25double precision :: kpt(3,nkpt),kptq(3,nkptq),shiftk(3),shiftkq(3)
26integer :: nsymk(nkpt),nsymkq(nkptq),symk(nkpt,nsym*2),symkq(nkptq,nsymq*2)
27integer :: symrel(3,3,nsym),syminv(3,3,nsymq)
28integer :: lvtrans(3,ngkpt(1),ngkpt(2),ngkpt(3))
29integer :: lvtransq(3,ngkpt(1),ngkpt(2),ngkpt(3))
30integer :: ihlf(nkpt),ihlfq(nkptq)
31double precision :: bmet(3,3),blat(3,3)
32integer :: ipwx(3,npwx),ipwndx(2,napwndx)
33integer :: igndx(nkpt,igmn(1):igmx(1),igmn(2):igmx(2),igmn(3):igmx(3))
34integer :: igndxq(nkptq,igmn(1):igmx(1),igmn(2):igmx(2),igmn(3):igmx(3))
35integer :: ikndx(ngkpt(1),ngkpt(2),ngkpt(3)),ikndxq(ngkpt(1),ngkpt(2),ngkpt(3))
36integer :: iqndx(ngkpt(1),ngkpt(2),ngkpt(3))
37integer :: isymndx(ngkpt(1),ngkpt(2),ngkpt(3)),isymndxq(ngkpt(1),ngkpt(2),ngkpt(3))
38integer :: nband(nkpt*nsppol),nbandq(nkptq*nsppol)
39double complex :: sei,seipw,seiold,seibc,csepw,cseold,csebc
40double complex, allocatable :: ssi(:)
41double complex :: ssc(-nwpt:nwpt),cse,dcse,cse1(-nwpt:nwpt),zz,ctest,cse2(-nwpt:nwpt)
42double precision :: xse,ssx,ssx1,ssx2,xse2
43integer :: ii,jj,kk,ix,iy,iz,iskip,isign
44integer :: iqq(3),iqqp(3),jka(3),jkb(3),jkk(3),iqv(3),ictr(3)
45integer :: ikpt,ikptq,iks(3),ikslf1(3),ikslf2(3)
46integer :: igg(3),igh(3),igg0(3),isym,isymq,ibp,iqpt,iqsym,iw,iww,ie,je,ies
47double precision :: xck(3),xckq(3),qadj(6)
48double complex :: cmatel,cmatel2,amatel,amatel2,vqmat2(ngkpt(1),ngkpt(2),ngkpt(3))
49double complex :: vfactor(ngkpt(1),ngkpt(2),ngkpt(3),nwpt),vfv(8,nwpt),vfx(8)
50double precision :: vq2(8),vqvtx0(4),vqvtx(4),qkcvt(3,3)
51double precision :: omega(ngkpt(1),ngkpt(2),ngkpt(3))
52double precision :: whi,wlo,wwhi,wwlo,ww,www,enval(8),dw,eshift,evtx0(4),evtx(4)
53double precision :: rrpyr(3,4),vpyr0(4,nwpt),kvtx(3,4),xk(3,4),rg(3,3),tvol,vpyr(4),xpyr0(4),xpyr(4)
54double precision :: de21,de31,de32,de41,de42,de43,thresh
55double precision :: fbx(4),fb(4),cmx(3,4),cm(3,4),xkt(3,3)
56double precision :: aa0(nwpt),av(3,nwpt),xme,xv(3)
57double precision :: avec(3),bgrad(3),xmult
58double precision :: abr,rlr
59integer :: iwh,iwl,ibmin,ibmax,iwhi,iwlo,iehi,ielo
60integer :: indxe(4),iwwhi,iwwlo
61double precision :: vq(ngkpt(1),ngkpt(2),ngkpt(3)),qq(3),qp(3),qq2,qp2,qs(3),ek,ekmq
62double precision :: stvec(ngkpt(3)),stvec2(ngkpt(3)),temparray(ngkpt(1),ngkpt(2),ngkpt(3))
63integer :: iipw,jjpw
64integer :: iqsymndx(ngkpt(1),ngkpt(2),ngkpt(3)),invpw2ndx(npwx,npwx)
65integer :: pwsymndx(npwc,2*nsym)
66integer :: ipw2,jw
67double precision :: eps1,eps2,eps
68double precision :: eval,brd,esprd
69double precision :: gfo,gamma(3),pln(3),dist(3)
70double precision :: sint1,sint1a,sint1b,svec(3),sveca(3),svecb(3)
71double complex :: wint,wint0,w2int,w2int0,gterm,vcentr
72double precision :: wcentr(-nwpt:nwpt),wgrid(nwpt)
73double complex :: pwmatel(nlmn,nlmn,npwx,ngkpt(1),ngkpt(2),ngkpt(3)), &
74&                tpwmatel(nlmn,nlmn,npwx,ngkpt(1),ngkpt(2),ngkpt(3))
75logical :: lqsing,lqcentr
76integer :: idum
77logical :: lx,lc
78double precision :: xdum,vdum(3)
79double precision :: linterp
80double precision :: rr(3,8)
81integer :: ivndx(4,6),itet,iv
82character*4 :: label
83external linterp
84data rr(1:3,1) /0.0d0,0.0d0,0.0d0/
85data rr(1:3,2) /1.0d0,0.0d0,0.0d0/
86data rr(1:3,3) /0.0d0,1.0d0,0.0d0/
87data rr(1:3,4) /1.0d0,1.0d0,0.0d0/
88data rr(1:3,5) /0.0d0,0.0d0,1.0d0/
89data rr(1:3,6) /1.0d0,0.0d0,1.0d0/
90data rr(1:3,7) /0.0d0,1.0d0,1.0d0/
91data rr(1:3,8) /1.0d0,1.0d0,1.0d0/
92data ivndx(1:4,1) /1,2,3,5/
93data ivndx(1:4,2) /3,5,6,7/
94data ivndx(1:4,3) /2,3,5,6/
95data ivndx(1:4,4) /2,3,4,6/
96data ivndx(1:4,5) /3,4,6,7/
97data ivndx(1:4,6) /4,6,7,8/
98
99abr=1.d-6
100rlr=1.d-6
101sei=0.d0
102xse=(0.d0,0.d0)
103ctest=(0.d0,0.d0)
104cse=(0.d0,0.d0)
105zz=(0.d0,0.d0)
106igg0=(/0,0,0/)
107do ie=-nwpt,nwpt
108  cse1(ie)=(0.d0,0.d0)
109enddo
110xse=0.d0
111dw=wmax/dble(nwpt)
112do iw=1,nwpt
113  wgrid(iw)=iw*dw
114enddo
115ikpt=ikndx(ikk(1),ikk(2),ikk(3))
116isym=isymndx(ikk(1),ikk(2),ikk(3))
117ek=enrgy(indxkbnd(ikpt)+iband)
118do ii=1,3
119do jj=1,3
120  qkcvt(ii,jj)=blat(ii,jj)/dble(ngkpt(ii))
121enddo
122enddo
123
124do ix=1,ngkpt(1)
125do iy=1,ngkpt(2)
126do iz=1,ngkpt(3)
127  temparray(ix,iy,iz)=(0.d0,0.d0)
128enddo
129enddo
130enddo
131
132gamma=(blat(1,:)+blat(2,:)+blat(3,:))/2.d0
133do ii=1,3
134  jj=mod(ii,3)+1
135  kk=mod(ii+1,3)+1
136  pln(1)=blat(jj,2)*blat(kk,3)-blat(jj,3)*blat(kk,2)
137  pln(2)=-blat(jj,1)*blat(kk,3)+blat(jj,3)*blat(kk,1)
138  pln(3)=blat(jj,1)*blat(kk,2)-blat(jj,2)*blat(kk,1)
139  dist(ii)=dot_product(gamma,pln)/sqrt(dot_product(pln,pln))
140enddo
141gfo=minval(dist)
142
143!ipaw=0
144!write(6,'(a)') '      finding contibutions from:'
145!write(6,'(a)') 'band  plane waves    correlation                exchange'
146do ibp=nbcore+1,ncband
147!do ibp=1,1
148!do ibp=16,16
149  seibc=sei
150  cse2=cse1
151  xse2=xse
152  if (ibp.gt.nbocc) then
153    isign=-1
154  else
155    isign=1
156  endif
157  do iipw=1,napwndx
158!  do iipw=1,1
159!  do iipw=15,15
160!    write(6,'(a,i4)') 'iipw ',iipw
161    seipw=sei
162    do ie=-nwpt,nwpt
163      ssc(ie)=(0.d0,0.d0)
164    enddo
165    ipw1=ipwndx(1,iipw)
166    ipw2=ipwndx(2,iipw)
167!    if (ipw1.ne.ipw2) cycle
168!    write(6,'(38x,i3,14x,2i3)') ibp,ipw1,ipw2
169    lx=ipw1.eq.ipw2.and.ibp.le.nbocc
170    lc=iipw.le.ntpwndx
171    if ((.not.lx).and.(.not.lc)) cycle
172
173    do ix=1,ngkpt(1)
174    do iy=1,ngkpt(2)
175    do iz=1,ngkpt(3)
176!    do ix=2,2
177!    do iy=2,2
178!    do iz=4,4
179      iqq=(/ix,iy,iz/)
180!write(6,*) '>>>>> iqq = ',iqq
181!write(6,'(a,3f10.5)') 'kpt = ',kpt(:,ikpt)
182      iks=iqq-ngkpt/2
183      jka=ikk-iks
184      jkk=mod(jka,ngkpt(ii))
185      do ii=1,3
186        if (jkk(ii).le.0) jkk(ii)=jkk(ii)+ngkpt(ii)
187      enddo
188      if (iks(1).eq.0.and.iks(2).eq.0.and.iks(3).eq.0) then
189!      if (.true.) then
190        ikptq=ikndxq(jkk(1),jkk(2),jkk(3))
191        jkb=nint(kptq(:,ikptq)*dble(ngkpt))+ngkpt/2
192        qq=kpt(:,ikpt)-kptq(:,ikptq)
193!write(6,'(a,3f10.5)') 'kptq = ',kptq(:,ikptq)
194      else
195        ikptq=ikndx(jkk(1),jkk(2),jkk(3))
196        jkb=nint(kpt(:,ikptq)*dble(ngkpt))+ngkpt/2
197        qq=kpt(:,ikpt)-kpt(:,ikptq)
198!write(6,'(a,3f10.5)') 'kptq = ',kptq(:,ikptq)
199      endif
200      igg=nint(dble(jkb-jka)/dble(ngkpt))
201      qp=qq+igg+ipwx(:,ipw2)
202      qq=qq+igg+ipwx(:,ipw1)
203      qq2=0.d0
204      qp2=0.d0
205      do ii=1,3
206      do jj=1,3
207        qq2=qq2+qq(ii)*bmet(ii,jj)*qq(jj)
208        qp2=qp2+qp(ii)*bmet(ii,jj)*qp(jj)
209      enddo
210      enddo
211      vq(ix,iy,iz)=4.d0*pi/sqrt(qq2*qp2)
212!write(6,'(a,3f10.5)') 'qq = ',qq
213!write(6,'(a,3f10.5)') 'qp = ',qp
214!write(6,'(a,3f10.5)') 'qq2 = ',qq2
215!write(6,'(a,3f10.5)') 'qp2 = ',qp2
216!write(6,*) 'vq = ',vq(ix,iy,iz)
217!write(6,*) 'jka = ',jka
218!write(6,*) 'jkk = ',jkk
219!write(6,*) 'jkb = ',jkb
220!write(6,*) 'igg = ',igg
221!write(6,*) 'ikpt = ',ikpt
222!write(6,*) 'ikptq = ',ikptq
223      if (iks(1).eq.0.and.iks(2).eq.0.and.iks(3).eq.0) then
224        call mkmatelX(ibp,iband,ikpt,ikptq,ipwx(:,ipw1),igg, &
225&         ncg,ncgq,nkpt,nkptq,npwt,npwtq,igmx,igmn,igndx,igndxq, &
226&         isym,isymq,symrel,syminv,nsym,nsymq,ihlf,ihlfq,kpt,kptq, &
227&         lvtrans(1:3,ikk(1),ikk(2),ikk(3)),lvtransq(1:3,ikk(1),ikk(2),ikk(3)), &
228&         cg,cgq,indxkcg,indxkcgq,indxkpw,indxkpwq,npwarr,npwarrq,kg,kgq, &
229&         cmatel)
230        if (ipaw.ne.0) then
231!          call mkmatelP(pi,xred,natom,ibp,iband,ikpt,ikptq,qq,igg0,ngkpt, &
232!&               pwmatel(:,:,ipw1,ix,iy,iz),tpwmatel(:,:,ipw1,ix,iy,iz), &
233!&               projwf,nlmn,nkpt,nkptq,ncband,amatel)
234        else
235          amatel=(0.d0,0.d0)
236        endif
237        if (ipw1.eq.ipw2) then
238          cmatel2=cmatel
239          amatel2=amatel
240        else
241          call mkmatelX(ibp,iband,ikpt,ikptq,ipwx(:,ipw2),igh, &
242&         ncg,ncgq,nkpt,nkptq,npwt,npwtq,igmx,igmn,igndx,igndxq, &
243&         isym,isymq,symrel,syminv,nsym,nsymq,ihlf,ihlfq,kpt,kptq, &
244&         lvtrans(1:3,ikk(1),ikk(2),ikk(3)),lvtransq(1:3,ikk(1),ikk(2),ikk(3)), &
245&         cg,cgq,indxkcg,indxkcgq,indxkpw,indxkpwq,npwarr,npwarrq,kg,kgq, &
246&         cmatel2)
247          if (ipaw.ne.0) then
248!            call mkmatelP(pi,xred,natom,ibp,iband,ikpt,ikptq,qp,igg0,ngkpt, &
249!&                 pwmatel(:,:,ipw2,ix,iy,iz),tpwmatel(:,:,ipw2,ix,iy,iz), &
250!&                 projwf,nlmn,nkpt,nkptq,ncband,amatel2)
251          else
252            amatel2=(0.d0,0.d0)
253          endif
254        endif
255        omega(ix,iy,iz)=enrgyq(indxkbndq(ikptq)+ibp)-ek+dw/2
256      else
257        call mkmatelX(ibp,iband,ikpt,ikptq,ipwx(:,ipw1),igg, &
258&         ncg,ncg,nkpt,nkpt,npwt,npwt,igmx,igmn,igndx,igndx, &
259&         isym,isymq,symrel,syminv,nsym,nsym,ihlf,ihlf,kpt,kpt, &
260&         lvtrans(1:3,ikk(1),ikk(2),ikk(3)),lvtrans(1:3,ikk(1),ikk(2),ikk(3)), &
261&         cg,cg,indxkcg,indxkcg,indxkpw,indxkpw,npwarr,npwarr,kg,kg, &
262&         cmatel)
263        if (ipaw.ne.0) then
264!          call mkmatelP(pi,xred,natom,ibp,iband,ikpt,ikptq,qq,igg0,ngkpt, &
265!&               pwmatel(:,:,ipw1,ix,iy,iz),tpwmatel(:,:,ipw1,ix,iy,iz), &
266!&               projwf,nlmn,nkpt,nkpt,ncband,amatel)
267        else
268          amatel=(0.d0,0.d0)
269        endif
270        if (ipw1.eq.ipw2) then
271          cmatel2=cmatel
272          amatel2=amatel
273        else
274          call mkmatelX(ibp,iband,ikpt,ikptq,ipwx(:,ipw2),igh, &
275&         ncg,ncg,nkpt,nkpt,npwt,npwt,igmx,igmn,igndx,igndx, &
276&         isym,isymq,symrel,syminv,nsym,nsym,ihlf,ihlf,kpt,kpt, &
277&         lvtrans(1:3,ikk(1),ikk(2),ikk(3)),lvtrans(1:3,ikk(1),ikk(2),ikk(3)), &
278&         cg,cg,indxkcg,indxkcg,indxkpw,indxkpw,npwarr,npwarr,kg,kg, &
279&         cmatel2)
280          if (ipaw.ne.0) then
281!            call mkmatelP(pi,xred,natom,ibp,iband,ikpt,ikptq,qp,igg0,ngkpt, &
282!&                 pwmatel(:,:,ipw2,ix,iy,iz),tpwmatel(:,:,ipw2,ix,iy,iz), &
283!&                 projwf,nlmn,nkpt,nkpt,ncband,amatel2)
284          else
285            amatel2=(0.d0,0.d0)
286          endif
287        endif
288        omega(ix,iy,iz)=enrgy(indxkbnd(ikptq)+ibp)-ek+dw/2
289      endif
290      vqmat2(ix,iy,iz)=(cmatel+amatel)*conjg(cmatel2+amatel2)*vq(ix,iy,iz)
291!write(6,*) 'cmatel  = ',cmatel
292!write(6,*) 'amatel  = ',amatel
293!write(6,*) 'cmatel+amatel  = ',cmatel+amatel
294!write(6,*) '|cmatel+amatel|^2  = ',dble((cmatel+amatel)*conjg(cmatel+amatel))
295!write(6,*) 'cmatel2 = ',cmatel2
296!write(6,*) 'vqmat2 = ',vqmat2(iqq(1),iqq(2),iqq(3))
297!write(56,*) '>>>>>',iqq,vqmat2(iqq(1),iqq(2),iqq(3))
298!      stvec(iqq(3))=vq(ix,iy,iz)
299!      stvec(iqq(3))=sqrt(qq2)
300!      stvec(iqq(3))=dble(vqmat2(iqq(1),iqq(2),iqq(3)))
301!      stvec2(iqq(3))=dimag(vqmat2(iqq(1),iqq(2),iqq(3)))
302!      stvec(iqq(3))=dble(cmatel)
303!      stvec2(iqq(3))=dimag(cmatel)
304!      stvec(iqq(3))=dble(cmatel2)
305!      stvec2(iqq(3))=dimag(cmatel2)
306!      stvec(iqq(3))=dble(cmatel*conjg(cmatel))
307!      stvec(iqq(3))=dble(amatel*conjg(amatel))
308!      stvec(iqq(3))=dble((cmatel+amatel)*conjg(cmatel+amatel))
309!      temparray(ix,iy,iz)=temparray(ix,iy,iz)+stvec(iqq(3))
310    enddo
311!    write(76,'(10f8.2)') stvec(:)
312!    write(76,'(10es8.1)') stvec(:)
313!    write(76,'(10es8.1)') stvec2(:)
314!    write(76,'(10f8.5)') 10.d0**(log10(stvec(:))-int(log10(stvec(:))))
315!    write(76,*)
316    enddo
317!    write(76,*) iqq(1)+1,'--------------------------------------------------',ibp
318    enddo
319!stop
320    do ix=1,ngkpt(1)
321    do iy=1,ngkpt(2)
322    do iz=1,ngkpt(3)
323      iqq=(/ix,iy,iz/)
324      call fhilo(omega,iqq,ngkpt,whi,wlo)
325      if (ix.eq.1.and.iy.eq.1.and.iz.eq.1) then
326        wwhi=whi
327        wwlo=wlo
328      else
329        wwhi=max(wwhi,whi)
330        wwlo=min(wwlo,wlo)
331      endif
332    enddo
333    enddo
334    enddo
335
336    iwlo=nint(wwlo*dble(nwpt)/wmax)
337    iwhi=nint(wwhi*dble(nwpt)/wmax)
338    if (ibp.gt.nbocc) then
339      ielo=iwlo+1
340      iehi=iwhi+nwpt
341    else
342      ielo=iwlo-nwpt
343      iehi=iwhi-1
344    endif
345    allocate (ssi(ielo:iehi))
346    do ie=ielo,iehi
347      ssi(ie)=0.d0
348    enddo
349    ssx=0.d0
350
351    if (lc) then
352      do iw=1,nwpt
353        do ix=1,ngkpt(1)
354        do iy=1,ngkpt(2)
355        do iz=1,ngkpt(3)
356          iqq=(/ix,iy,iz/)
357          iqpt=iqndx(iqq(1),iqq(2),iqq(3))
358          call locateelement(iqq,ipw1,ipw2,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpw)
359          if (jjpw.ne.0) then
360            vfactor(iqq(1),iqq(2),iqq(3),iw)=vqmat2(iqq(1),iqq(2),iqq(3))*dimag(lossfn(iw,iqpt,jjpw))
361          else
362            vfactor(iqq(1),iqq(2),iqq(3),iw)=(0.d0,0.d0)
363          endif
364        enddo
365        enddo
366        enddo
367      enddo
368    endif
369
370    if (ipw1.eq.1.or.ipw2.eq.1) then
371      ix=ngkpt(1)/2
372      iy=ngkpt(2)/2
373      iz=ngkpt(3)/2
374      qadj(1)=abs(vqmat2(ix,iy,iz)/vqmat2(ix+1,iy,iz))
375      qadj(2)=abs(vqmat2(ix,iy,iz)/vqmat2(ix-1,iy,iz))
376      qadj(3)=abs(vqmat2(ix,iy,iz)/vqmat2(ix,iy+1,iz))
377      qadj(4)=abs(vqmat2(ix,iy,iz)/vqmat2(ix,iy-1,iz))
378      qadj(5)=abs(vqmat2(ix,iy,iz)/vqmat2(ix,iy,iz+1))
379      qadj(6)=abs(vqmat2(ix,iy,iz)/vqmat2(ix,iy,iz-1))
380      if (qadj(1).gt.1.d3.and.qadj(2).gt.1.d3.and.qadj(3).gt.1.d3.and.  &
381&         qadj(4).gt.1.d3.and.qadj(5).gt.1.d3.and.qadj(6).gt.1.d3) then
382        lqsing=.true.
383      else
384        lqsing=.false.
385      endif
386    else
387      lqsing=.false.
388    endif
389
390    ictr=(/ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2/)
391    do ix=1,ngkpt(1)
392    do iy=1,ngkpt(2)
393    do iz=1,ngkpt(3)
394!    do ix=6,6
395!    do iy=5,5
396!    do iz=5,5
397!    do ix=1,1
398!    do iy=1,1
399!    do iz=1,1
400      iqq=(/ix,iy,iz/)
401      ssx1=ssx
402      call fval(omega,iqq,iqqp,ngkpt,enval)
403      if (lqsing) call fval(vq,iqq,iqqp,ngkpt,vq2)
404      if (lc) then
405        do iw=1,nwpt
406          call fpol(vfactor(:,:,:,iw),ngkpt,iqq,iqqp,vfv(:,iw))
407        enddo
408      endif
409      if (lx) call fpol(vqmat2(:,:,:),ngkpt,iqq,iqqp,vfx)
410      do itet=1,6
411!      do itet=1,1
412        ssx2=ssx
413        lqcentr=.false.
414        do iv=1,4
415          evtx0(iv)=enval(ivndx(iv,itet))
416          rrpyr(1:3,iv)=rr(1:3,ivndx(iv,itet))
417          if (lc) then
418            do iw=1,nwpt
419              vpyr0(iv,iw)=vfv(ivndx(iv,itet),iw)
420            enddo
421          endif
422          if (lx) xpyr0(iv)=vfx(ivndx(iv,itet))
423          if (lqsing.and..not.lqcentr) then
424            iqv=iqq+nint(rrpyr(:,iv))
425            if (iqv(1).eq.ictr(1).and.iqv(2).eq.ictr(2).and.iqv(3).eq.ictr(3)) then
426              lqcentr=.true.
427            endif
428          endif
429        enddo
430        call indxhpsort(4,4,evtx0,indxe)
431        evtx=evtx0(indxe)    ! order vertices by energy
432!write(6,'(4es15.5)') evtx
433!        kvtx(1:3,1:4)=rrpyr(1:3,indxe)
434        do ii=1,4
435          jj=indxe(ii)
436          do kk=1,3
437            kvtx(kk,ii)=dot_product(iqq+rrpyr(:,jj)-ictr,qkcvt(:,kk))
438          enddo
439        enddo
440        xk(1:3,1)=(/0.d0,0.d0,0.d0/)
441        do ii=2,4
442          xk(1:3,ii)=kvtx(1:3,ii)-kvtx(1:3,1)
443        enddo
444        call cross(xk(:,3),xk(:,4),vdum)
445        tvol=dot_product(vdum,xk(:,2))
446        do jj=1,3  ! make contragradient
447          call cross(xk(1:3,mod(jj,3)+2),xk(1:3,mod(jj+1,3)+2),rg(1:3,jj))
448        enddo
449        rg=rg/tvol
450        tvol=abs(tvol)
451        iwwlo=nint(evtx(1)/dw)
452        iwwhi=nint(evtx(4)/dw)
453        de21=evtx(2)-evtx(1)
454        de31=evtx(3)-evtx(1)
455        de32=evtx(3)-evtx(2)
456        de41=evtx(4)-evtx(1)
457        de42=evtx(4)-evtx(2)
458        de43=evtx(4)-evtx(3)
459        thresh=1.d-5*de41
460        if (lqcentr) then            ! integral around coulomb singularity
461!        if (.true.) then
462          do iv=1,4
463            vqvtx0(iv)=vq2(ivndx(iv,itet))
464          enddo
465          bgrad=(/0.d0,0.d0,0.d0/)   ! energy gradient
466          do jj=1,3
467            bgrad=bgrad+(evtx(jj+1)-evtx(1))*rg(:,jj)
468          enddo
469          xmult=4.d0*pi/sqrt(dot_product(bgrad,bgrad))
470          if (lc) then
471            do iw=1,nwpt
472              vpyr(:)=vpyr0(:,iw)/vqvtx0
473              aa0(iw)=vpyr(indxe(1))
474              av(1:3,iw)=(/0.d0,0.d0,0.d0/)
475              do jj=1,3
476                av(1:3,iw)=av(1:3,iw) &
477&                         +(vpyr(indxe(jj+1))-vpyr(indxe(1)))*rg(1:3,jj)
478              enddo
479            enddo
480          endif
481          if (lx) then
482            xpyr=xpyr0/vqvtx0
483            xme=xpyr(indxe(1))
484            xv=(/0.d0,0.d0,0.d0/)
485            do jj=1,3
486              xv=xv+(xpyr(indxe(jj+1))-xpyr(indxe(1)))*rg(:,jj)
487            enddo
488          endif
489          do iww=iwwlo,iwwhi
490            www=dw*iww
491            if (www.lt.evtx(1)) then
492              cycle
493            elseif (www.lt.evtx(2)) then
494              call singint(1,www,xk,kvtx(:,1),evtx,sint1,svec)
495            elseif (www.lt.evtx(3)) then
496              if (de21.lt.thresh.and.de43.lt.thresh) then
497                call singint2(www,xk,kvtx(:,1),evtx,sint1b,svecb)
498                sint1a=0.d0
499                sveca=(/0.d0,0.d0,0.d0/)
500              elseif (de21.ge.de43) then
501                call singint(2,www,xk,kvtx(:,1),evtx,sint1a,sveca)
502                call singint(1,www,xk,kvtx(:,1),evtx,sint1b,svecb)
503              else
504                call singint(3,www,xk,kvtx(:,1),evtx,sint1a,sveca)
505                call singint(4,www,xk,kvtx(:,1),evtx,sint1b,svecb)
506              endif
507              sint1=sint1b-sint1a
508              svec=svecb-sveca
509            elseif (www.lt.evtx(4)) then
510              call singint(4,www,xk,kvtx(:,1),evtx,sint1,svec)
511            else
512              cycle
513            endif
514            sint1=sint1*xmult
515            svec=svec*xmult
516            if (lc) then
517              do iw=1,nwpt
518                ie=iww-isign*iw
519                ssi(ie)=ssi(ie)+(aa0(iw)*sint1+dot_product(av(:,iw),svec))*dw
520              enddo
521            endif
522            if (lx) ssx=ssx+(xme*sint1+dot_product(xv,svec))*dw
523          enddo
524        else                         ! non-singular, tetrahedron method
525          fbx(1)=tvol/(2.d0*de21*de31*de41)
526          if (de21.ge.de43) then
527            fbx(2)=tvol/(2.d0*de21*de32*de42)
528          else
529            fbx(3)=tvol/(2.d0*de31*de32*de43)
530          endif
531          fbx(4)=tvol/(2.d0*de41*de42*de43)
532          do ii=1,4
533            cmx(1:3,ii)=(/0.d0,0.d0,0.d0/)
534            do jj=1,4
535              if (jj.ne.ii) cmx(1:3,ii)=cmx(1:3,ii)+(xk(1:3,jj)-xk(1:3,ii)) &
536&                                                  /(3*(evtx(jj)-evtx(ii)))
537            enddo
538          enddo
539          if (lc) then
540            do iw=1,nwpt
541              aa0(iw)=vpyr0(indxe(1),iw)   ! base value of function
542              av(1:3,iw)=(/0.d0,0.d0,0.d0/) ! gradient of function
543              do jj=1,3
544                av(1:3,iw)=av(1:3,iw) &
545&                         +(vpyr0(indxe(jj+1),iw)-vpyr0(indxe(1),iw))*rg(1:3,jj)
546              enddo
547            enddo
548          endif
549          if (lx) then
550            xme=xpyr0(indxe(1))
551            xv=(/0.d0,0.d0,0.d0/)
552            do jj=1,3
553              xv=xv+(xpyr0(indxe(jj+1))-xpyr0(indxe(1)))*rg(:,jj)
554            enddo
555          endif
556          do iww=iwwlo,iwwhi
557            www=dw*iww
558            if (www.lt.evtx(1)) then
559              cycle
560            elseif (www.lt.evtx(2)) then
561label='1  F'
562              sint1=fbx(1)*(www-evtx(1))**2
563              svec=(xk(1:3,1)+(www-evtx(1))*cmx(1:3,1))*sint1
564            elseif (www.lt.evtx(3)) then
565              if (de21.lt.thresh.and.de43.lt.thresh) then
566                xkt(:,1)=xk(:,3)*(www-evtx(1))/de31
567                xkt(:,2)=xk(:,4)*(www-evtx(1))/de41
568                xkt(:,3)=(xk(:,4)-xk(:,2))*(www-evtx(2))/de42+xk(:,2)
569                sint1=tvol*(2*(www-evtx(1))/(de31*de41) &
570&                     -(www-evtx(1))**2*(de31+de41)/(de31*de41)**2)/2
571                svec=((xkt(:,1)+xkt(:,3))/2.d0)*sint1
572              elseif (de21.ge.de43) then
573                fb(1)=fbx(1)*(www-evtx(1))**2
574                fb(2)=fbx(2)*(www-evtx(2))**2
575                sint1=fb(1)-fb(2)
576                cm(1:3,1)=xk(1:3,1)+(www-evtx(1))*cmx(1:3,1)
577                cm(1:3,2)=xk(1:3,2)+(www-evtx(2))*cmx(1:3,2)
578                svec=cm(:,1)*fb(1)-cm(:,2)*fb(2)
579              else
580                fb(3)=fbx(3)*(www-evtx(3))**2
581                fb(4)=fbx(4)*(www-evtx(4))**2
582                sint1=fb(4)-fb(3)
583                cm(1:3,3)=xk(1:3,3)+(www-evtx(3))*cmx(1:3,3)
584                cm(1:3,4)=xk(1:3,4)+(www-evtx(4))*cmx(1:3,4)
585                svec=cm(:,4)*fb(4)-cm(:,3)*fb(3)
586              endif
587            elseif (www.lt.evtx(4)) then
588              sint1=fbx(4)*(www-evtx(4))**2
589              svec=(xk(1:3,4)+(www-evtx(4))*cmx(1:3,4))*sint1
590            else
591              cycle
592            endif
593            if (lc) then
594              do iw=1,nwpt
595                ie=iww-isign*iw
596                ssi(ie)=ssi(ie)+(aa0(iw)*sint1+dot_product(av(:,iw),svec))*dw
597              enddo
598            endif
599            if (lx) ssx=ssx+(xme*sint1+dot_product(xv,svec))*dw
600          enddo
601        endif
602!if (lqcentr) then
603!write(47,'(3i3,3x,i3,f16.8,3x,a)') ix,iy,iz,itet,ssx-ssx2,'T'
604!else
605!write(47,'(3i3,3x,i3,f16.8,3x,a)') ix,iy,iz,itet,ssx-ssx2,'F'
606!endif
607      enddo
608!      temparray(ix,iy,iz)=(ssx-ssx1)/((2*pi)**3)
609    enddo
610    enddo
611    enddo
612!!!!!!!!!!    ssi=ssi/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol*pi)
613    if (lc) ssi=-isign*ssi/((2*pi)**3*pi)
614    if (lx) ssx=-ssx/((2*pi)**3)
615
616    if (lc) then
617      do ie=-nwpt,nwpt
618        eps1=dble(ie)*dw-dw/2
619!        ssc(ie)=(0.d0,0.d0)
620        do je=ielo,iehi
621          if (ie.eq.je) then
622            ssc(ie)=ssc(ie)+(0.d0,1.d0)*pi*ssi(je)
623          else
624            eps2=dble(je)*dw-dw/2
625            ssc(ie)=ssc(ie)+isign*ssi(je)*dw/(eps2-eps1)
626          endif
627!write(6,'(2i4,2f10.3,es12.3,2(2x,2es11.3))') ie,je,eps1*27.2114,eps2*27.2114,1/(eps2-eps1),ssc(ie)*27.2114,dble(ssi(je)*27.2114)
628        enddo
629      enddo
630    endif
631    deallocate (ssi)
632
633    if (lc) cse1=cse1+ssc
634    if (lx) xse=xse+ssx
635!    if (lx) then
636!      write(6,'(i4,4x,i5,3x,3i3,5x,"(",f10.6,",",f10.6,")",5x,f10.6)')  &
637!&          ibp,ipw1,ipwx(:,ipw1),((ssc(1)+ssc(0))/2)*27.2114,dble(ssx)*27.2114
638!    else
639!      write(6,'(i4,4x,i5,3x,3i3,5x,"(",f10.6,",",f10.6,")")') ibp,ipw1,ipwx(:,ipw1), &
640!&          ((ssc(1)+ssc(0))/2)*27.2114
641!    endif
642
643  enddo
644!  if (lx) then
645!    write(6,'(i4,5x,"(",f10.6,",",f10.6,")",5x,f10.6)') ibp, &
646!& (cse1(1)+cse1(0)-cse2(1)-cse2(0))*27.2114/2.d0, &
647!& dble(xse-xse2)*27.2114
648!  else
649!    write(6,'(i4,5x,"(",f10.6,",",f10.6,")",f10.6)') ibp, &
650!& (cse1(1)+cse1(0)-cse2(1)-cse2(0))*27.2114/2.d0
651!  endif
652enddo
653
654do ie=-nwpt,nwpt
655  eps1=dble(ie)*dw-dw/2
656!  write(12,'(i4,f10.3,2(3x,2es12.3))') ie,eps1*27.2114,cse1(ie)*27.2114
657enddo
658cse=(cse1(0)+cse1(1))/2
659dcse=(cse1(1)-cse1(0))/dw
660zz=1.d0/(1.d0-dcse)
661!write(6,*) cse*27.2114
662!write(6,*) xse*27.2114
663!write(6,*) zz
664!write(6,*) dble(zz)*dimag(cse)*27.2114
665
666!temparray=log10(temparray)
667!temparray=10**(temparray-int(temparray))
668!do ix=1,ngkpt(1)
669!  write(76,*) ix,'--------------------------------------------------'
670!  do iy=1,ngkpt(2)
671!!    write(76,'(10es8.1)') dble(temparray(ix,iy,:))
672!    write(76,'(10f8.5)') dble(temparray(ix,iy,:))
673!  enddo
674!enddo
675
676return
677end subroutine mk2cse
678
679!***************************************************************************
680!
681! Over a surface S defined by the triangle with vertices xkvtx+xkp(1:3,1)
682!                                                        xkvtx+xkp(1:3,2)
683!                                                        xkvtx+xkp(1:3,3)
684! where xkp defined to give S of constant energy
685! (see G. Lehmann and M. Taut, Phys. stat. sol. (b) 54 469 (1972))
686! and vector k measured from tip of xkvtx
687! find      sint1 = integral dS 1/q^2
688! and       svec  = integral dS k/q^2
689! where q=k+xkvtx
690! On S, define coordinates X and Y such that k=xkp(1:3,3)+Y*xq1+X*xq2
691! normal=cross product(xq1,xq2)
692! dot_product(k+xkvtx,k+xkvtx) = aa*Y^2+Y*(bb+cc*X)+dd*X+ee*X^2+FF
693! integral dS -> |normal|*integral^1_0 dX integral^{1-X}_0 dY
694! integral^{1-X}_0 dY 1/(dot_product(k+xkvtx,k+xkvtx)) = xfn1
695! integral^{1-X}_0 dY X/(dot_product(k+xkvtx,k+xkvtx)) = X*xfn1 = xfn3
696! integral^{1-X}_0 dY Y/(dot_product(k+xkvtx,k+xkvtx)) = (xfn2-bb*xfn1-cc*xfn3)/(2*aa)
697! numerically integrate xfn1,xfn2,xfn3
698! sint1 = integral^1_0 dX xfn1
699! sint2 = integral^1_0 dX xfn2
700! sint3 = integral^1_0 dX xfn3
701! svec = xkp(:,3)*sint1+xq1*(sint2-bb*sint1-cc*sint3)/(2*aa)+xq2*sint3
702
703subroutine singint(ive,ww,xk,xkvtx,evtx,sint1,svec)
704implicit none
705integer :: ive
706double precision :: xk(3,4),xkvtx(3),evtx(4),ww
707double precision :: xkp(3,3),xq1(3),xq2(3),xkv0(3),xcvt
708double precision :: aa,bb,cc,dd,ee,ff,qq
709double precision :: xmin,xmax,abr,rlr,xsing(20),error
710double precision :: xfn1,xfn2,xfn3,grater
711double precision :: sint1,sint2,sint3,svec(3),root1,root2
712double precision :: xnorm(3),area
713double precision :: xdum
714integer :: ii,jj,ll
715integer :: nsing,numcal,maxns,nroots
716integer :: iter  !######
717common /fn/ aa,bb,cc,dd,ee,ff,iter !######
718external xfn1,xfn2,xfn3,grater
719
720!write(6,*) ive
721!write(6,'(4es15.5)') evtx
722!write(6,'(3f10.6)') xk(:,1)
723!write(6,*)
724!write(6,*) 'boundary vectors of tetrahedron'
725!write(6,'(3f10.6)') xk(:,2)
726!write(6,'(3f10.6)') xk(:,3)
727!write(6,'(3f10.6)') xk(:,4)
728  if (ive.eq.1) then
729    xkp(:,1)=(ww-evtx(1))*xk(:,2)/(evtx(2)-evtx(1))
730    xkp(:,2)=(ww-evtx(1))*xk(:,3)/(evtx(3)-evtx(1))
731    xkp(:,3)=(ww-evtx(1))*xk(:,4)/(evtx(4)-evtx(1))
732  elseif (ive.eq.2) then
733    xkp(:,1)=(ww-evtx(1))*xk(:,2)/(evtx(2)-evtx(1))
734    xkp(:,2)=(ww-evtx(2))*(xk(:,3)-xk(:,2))/(evtx(3)-evtx(2)) + xk(:,2)
735    xkp(:,3)=(ww-evtx(2))*(xk(:,4)-xk(:,2))/(evtx(4)-evtx(2)) + xk(:,2)
736  elseif (ive.eq.3) then
737    xkp(:,1)=(ww-evtx(1))*xk(:,3)/(evtx(3)-evtx(1))
738    xkp(:,2)=(ww-evtx(2))*(xk(:,3)-xk(:,2))/(evtx(3)-evtx(2)) + xk(:,2)
739    xkp(:,3)=(ww-evtx(3))*(xk(:,4)-xk(:,3))/(evtx(4)-evtx(3)) + xk(:,3)
740  else
741    xkp(:,1)=(ww-evtx(2))*(xk(:,4)-xk(:,2))/(evtx(4)-evtx(2)) + xk(:,2)
742    xkp(:,2)=(ww-evtx(3))*(xk(:,4)-xk(:,3))/(evtx(4)-evtx(3)) + xk(:,3)
743    xkp(:,3)=(ww-evtx(1))*xk(:,4)/(evtx(4)-evtx(1))
744  endif
745!write(6,*)
746!write(6,*) 'vertices of S'
747!write(6,'(3f10.6)') xkp(:,1)
748!write(6,'(3f10.6)') xkp(:,2)
749!write(6,'(3f10.6)') xkp(:,3)
750  xq1=xkp(:,1)-xkp(:,3)
751  xq2=xkp(:,2)-xkp(:,3)
752!write(6,*)
753!write(6,*) 'boundary vectors xq'
754!write(6,'(3f10.6)') xq1
755!write(6,'(3f10.6)') xq2
756  xkv0=xkp(:,3)+xkvtx
757!write(6,'(3f10.6)') xkvtx
758!write(6,'(3f10.6)') xkv0
759  call cross(xq1,xq2,xnorm)
760  area=sqrt(dot_product(xnorm,xnorm))
761!write(6,*)
762!write(6,'(a,es15.5)') 'area  ',area
763  aa=dot_product(xq1,xq1)
764  bb=2*dot_product(xq1,xkv0)
765  cc=2*dot_product(xq1,xq2)
766  dd=2*dot_product(xq2,xkv0)
767  ee=dot_product(xq2,xq2)
768  ff=dot_product(xkv0,xkv0)
769!write(6,'("coef",6f10.6)') aa,bb,cc,dd,ee,ff
770  nsing=0
771  call rquadroots(ee,dd,ff,nroots,root1,root2)
772!write(6,'(i4,2es15.5)') nroots,root1,root2
773  if (nroots.gt.0.and.root1.gt.0.d0.and.root1.lt.1.d0) then
774    nsing=nsing+1
775    xsing(nsing)=root1
776  endif
777  if (nroots.gt.1.and.root2.gt.0.d0.and.root2.lt.1.d0) then
778    nsing=nsing+1
779    xsing(nsing)=root2
780  endif
781  call rquadroots(1.d0+cc/2+ee,bb/2+dd,ff,nroots,root1,root2)
782!write(6,'(i4,2es15.5)') nroots,root1,root2
783  if (nroots.gt.0.and.root1.gt.0.d0.and.root1.lt.1.d0) then
784    nsing=nsing+1
785    xsing(nsing)=root1
786  endif
787  if (nroots.gt.1.and.root2.gt.0.d0.and.root2.lt.1.d0) then
788    nsing=nsing+1
789    xsing(nsing)=root2
790  endif
791!write(6,'(i4,4es15.5)') nsing,xsing(1:nsing)
792  call hpsort(nsing,nsing,xsing(1:nsing))
793!write(6,'("sings",i4,4es15.5)') nsing,xsing(1:nsing)
794  xmin=0.d0
795  xmax=1.d0
796  abr=1.d-12
797  rlr=1.d-6
798  sint1=grater(xfn1,xmin,xmax,abr,rlr,nsing,xsing,error,numcal,maxns)*area
799  sint2=grater(xfn2,xmin,xmax,abr,rlr,nsing,xsing,error,numcal,maxns)*area
800  sint3=grater(xfn3,xmin,xmax,abr,rlr,nsing,xsing,error,numcal,maxns)*area
801!write(6,*)
802!write(6,'("sints",6f10.6)') sint1,sint2,sint3
803  svec=xkp(:,3)*sint1+xq1*(sint2-bb*sint1-cc*sint3)/(2*aa)+xq2*sint3
804!write(6,'("svec",6f10.6)') svec
805!tint=(aa0*sint1+dot_product(av,svec))/bmag
806
807end subroutine singint
808
809!***************************************************************************
810
811double precision function xfn1(xx)
812implicit none
813double precision :: xx,xx2
814double precision :: aa,bb,cc,dd,ee,ff
815double precision :: xt1,xt2,xl1,qq,dx
816double precision :: sqq,xterm,df,d2f
817integer :: iter
818common /fn/ aa,bb,cc,dd,ee,ff,iter
819!integer :: ict
820!data ict /0/
821!save ict
822
823  xx2=xx*xx
824  xt1=bb+cc*xx
825  xt2=dd*xx+ee*xx2+ff
826  dx=2.d0*aa*(1.d0-xx)
827  xl1=dx+xt1
828  qq=4.d0*aa*xt2-xt1**2
829  if (qq.lt.0.d0) then
830    sqq=sqrt(-qq)
831!    xfn1=(log(abs((sqq+xl1)/(sqq-xl1)))-log(abs((sqq+xt1)/(sqq-xt1))))/sqq
832    xfn1=(log(abs(((sqq+xl1)*(sqq-xt1))/((sqq-xl1)*(sqq+xt1)))))/sqq
833!write(iter,'(" B",es15.5,4x,2es15.5)') xx,xfn1
834  else
835!    df=2.d0*dx/(qq+xt1**2)
836!    d2f=df**2*xt1/2.d0
837!    if (d2f/df.lt.1.d-6) then
838!      xfn1=df-d2f
839!    else
840      sqq=sqrt(qq)
841      if (abs(dx/(qq+xt1**2)).gt.abs(1.d8*(dx**2*xt1/(qq+xt1**2)**2))) then
842! atan(z+d)->atan(z)+d/(1+z^2)-d^2z/(1+z^2)^2+..., d->0 (taylor expansion)
843! for small dx, xfn1->dx/(qq+xt1^2)
844        xfn1=2*(dx/(qq+xt1**2)-dx**2*xt1/(qq+xt1**2)**2)
845      elseif (abs(1.d0/xl1-1.d0/xt1).gt.abs(1.d8*qq*(1.d0/xl1**3-1.d0/xt1**3))) then
846! atan(z)->pi/2-1/z+1/(3z^3)-1/(5z^5)+..., z->infty
847! => at small sqq or large xt1, xl1, xfn1->-1/xl1+1/xt1
848        xfn1=2*(-1.d0/xl1+1.d0/xt1+(qq/3.d0)*(1.d0/xl1**3-1.d0/xt1**3))
849      else
850        xfn1=2*(atan(xl1/sqq)-atan(xt1/sqq))/sqq
851      endif
852!write(iter,'(" A",es15.5,3x,4es15.5)') xx,xfn1,xt1,bb,cc
853!    endif
854  endif
855!write(6,'(es15.5,4es15.5)') xx,qq,xl1,xt1,xfn1
856!write(6,'(es15.5,4es15.5)') xx-1.d0,2/(qq+xt1**2),df,df-d2f,xfn1
857!write(6,'(es15.5,4es15.5)') xx-1.d0,d2f/df,df,df-d2f,xfn1
858!write(47,'(es15.7,4es15.7)') xx-1.d0,2/(qq+xt1**2),df,df-d2f,xfn1
859!write(6,'(es15.7,4es15.7)') xx,2/(qq+xt1**2),qq,xt1**2,4.d0*aa*xt2
860!write(6,'(es15.7,4es15.7)') xx,qq,xl1,xt1,xfn1
861!write(47,'(es15.7,4es15.7)') xx,qq,xl1,xt1,xfn1
862!write(6,'(es15.5,4es15.5)') xx,dx,qq,xt1**2,xfn1
863!write(iter,'(es15.5,4x,2es15.5)') xx,xfn1
864!ict=ict+1
865!if (ict.gt.3000) stop
866
867end function xfn1
868
869!***************************************************************************
870
871double precision function xfn2(xx)
872implicit none
873double precision :: xx,xx2,omx
874double precision :: aa,bb,cc,dd,ee,ff
875double precision :: xt1,xt2,xl1,qq
876double precision :: xterm,sqq
877integer :: iter  !######
878common /fn/ aa,bb,cc,dd,ee,ff,iter !######
879
880  xx2=xx*xx
881  omx=1.d0-xx
882  xterm=(aa*omx*omx+omx*(bb+cc*xx)+dd*xx+ee*xx2+ff)/(dd*xx+ee*xx2+ff)
883  xfn2=log(abs(xterm))
884!  xt1=bb+cc*xx
885!  xt2=dd*xx+ee*xx2+ff
886!  xl1=2.d0*aa*omx+xt1
887!  qq=4.d0*aa*xt2-xt1**2
888!  if (qq.lt.0.d0) then
889!    sqq=sqrt(-qq)
890!    xfn2=log(abs(xterm))+xt1*(log(abs((sqq+xl1)/(sqq-xl1)))-log(abs((sqq+xt1)/(sqq-xt1))))/sqq
891!  else
892!    sqq=sqrt(qq)
893!    xfn2=log(abs(xterm))+2*xt1*(atan((xl1)/sqq)-atan(xt1/sqq))/sqq
894!  endif
895!!write(47,'(4es15.5)') xx,xfn2,(aa+cc+ee)*xx2+(bb+dd)*xx+ff,dd*xx+ee*xx2+ff
896
897end function xfn2
898
899!***************************************************************************
900
901double precision function xfn3(xx)
902implicit none
903double precision :: xfn1,xx
904external xfn1
905
906  xfn3=xx*xfn1(xx)
907
908end function xfn3
909
910!***************************************************************************
911!
912! Over a surface S defined by the parallelogram with vertices
913!           xkvtx+xkp(1:3,1)
914!           xkvtx+xkp(1:3,2)
915!           xkvtx+xkp(1:3,3)
916!           xkvtx+xkp(1:3,1)+(xkp(1:3,2)-xkp(1:3,3))
917! where xkp defined to give S of constant energy
918! (see G. Lehmann and M. Taut, Phys. stat. sol. (b) 54 469 (1972))
919! and vector k measured from tip of xkvtx
920! find      sint1 = integral dS 1/q^2
921! and       svec  = integral dS k/q^2
922! where q=k+xkvtx
923! On S, define coordinates X and Y such that k=xkp(1:3,3)+Y*xq1+X*xq2
924! normal=cross product(xq1,xq2)
925! dot_product(k+xkvtx,k+xkvtx) = aa*Y^2+Y*(bb+cc*X)+dd*X+ee*X^2+FF
926! integral dS -> |normal|*integral^1_0 dX integral^1_0 dY
927! integral^1_0 dY 1/(dot_product(k+xkvtx,k+xkvtx)) = xfn4
928! integral^1_0 dY X/(dot_product(k+xkvtx,k+xkvtx)) = X*xfn4 = xfn6
929! integral^1_0 dY Y/(dot_product(k+xkvtx,k+xkvtx)) = (xfn5-bb*xfn4-cc*xfn6)/(2*aa)
930! numerically integrate xfn4,xfn5,xfn6
931! sint1 = integral^1_0 dX xfn4
932! sint2 = integral^1_0 dX xfn5
933! sint3 = integral^1_0 dX xfn6
934! svec = xkp(:,3)*sint1+xq1*(sint2-bb*sint1-cc*sint3)/(2*aa)+xq2*sint3
935
936subroutine singint2(ww,xk,xkvtx,evtx,sint1,svec)
937implicit none
938double precision :: xk(3,4),xkvtx(3),evtx(4),ww
939double precision :: xkp(3,3),xq1(3),xq2(3),xkv0(3),xcvt
940double precision :: aa,bb,cc,dd,ee,ff,qq
941double precision :: xmin,xmax,abr,rlr,xsing(20),error
942double precision :: xfn4,xfn5,xfn6,grater
943double precision :: sint1,sint2,sint3,svec(3),root1,root2
944double precision :: xnorm(3),area
945double precision :: xdum
946integer :: ii,jj,ll
947integer :: nsing,numcal,maxns,nroots
948integer :: iter  !######
949common /fn/ aa,bb,cc,dd,ee,ff,iter !######
950external xfn4,xfn5,xfn6,grater
951
952!write(6,*) ive
953!write(6,'(4es15.5)') evtx
954!write(6,'(3f10.6)') xk(:,1)
955!write(6,*)
956!write(6,*) 'boundary vectors of tetrahedron'
957!write(6,'(3f10.6)') xk(:,2)
958!write(6,'(3f10.6)') xk(:,3)
959!write(6,'(3f10.6)') xk(:,4)
960  xkp(:,1)=(ww-evtx(1))*xk(:,3)/(evtx(3)-evtx(1))
961  xkp(:,2)=(ww-evtx(1))*xk(:,4)/(evtx(4)-evtx(1))
962  xkp(:,3)=(ww-evtx(2))*(xk(:,4)-xk(:,2))/(evtx(4)-evtx(2)) + xk(:,2)
963!write(6,*)
964!write(6,*) 'vertices of S'
965!write(6,'(3f10.6)') xkp(:,1)
966!write(6,'(3f10.6)') xkp(:,2)
967!write(6,'(3f10.6)') xkp(:,3)
968  xq1=xkp(:,1)-xkp(:,3)
969  xq2=xkp(:,2)-xkp(:,3)
970!write(6,*)
971!write(6,*) 'boundary vectors xq'
972!write(6,'(3f10.6)') xq1
973!write(6,'(3f10.6)') xq2
974  xkv0=xkp(:,3)+xkvtx
975!write(6,'(3f10.6)') xkvtx
976!write(6,'(3f10.6)') xkv0
977  call cross(xq1,xq2,xnorm)
978  area=sqrt(dot_product(xnorm,xnorm))
979!write(6,*)
980!write(6,'(a,es15.5)') 'area  ',area
981  aa=dot_product(xq1,xq1)
982  bb=2*dot_product(xq1,xkv0)
983  cc=2*dot_product(xq1,xq2)
984  dd=2*dot_product(xq2,xkv0)
985  ee=dot_product(xq2,xq2)
986  ff=dot_product(xkv0,xkv0)
987!write(6,'("coef",6f10.6)') aa,bb,cc,dd,ee,ff
988  nsing=0
989  call rquadroots(ee,dd,ff,nroots,root1,root2)
990!write(6,'(i4,2es15.5)') nroots,root1,root2
991  if (nroots.gt.0.and.root1.gt.0.d0.and.root1.lt.1.d0) then
992    nsing=nsing+1
993    xsing(nsing)=root1
994  endif
995  if (nroots.gt.1.and.root2.gt.0.d0.and.root2.lt.1.d0) then
996    nsing=nsing+1
997    xsing(nsing)=root2
998  endif
999  call rquadroots(1.d0+cc/2+ee,bb/2+dd,ff,nroots,root1,root2)
1000!write(6,'(i4,2es15.5)') nroots,root1,root2
1001  if (nroots.gt.0.and.root1.gt.0.d0.and.root1.lt.1.d0) then
1002    nsing=nsing+1
1003    xsing(nsing)=root1
1004  endif
1005  if (nroots.gt.1.and.root2.gt.0.d0.and.root2.lt.1.d0) then
1006    nsing=nsing+1
1007    xsing(nsing)=root2
1008  endif
1009!write(6,'(i4,4es15.5)') nsing,xsing(1:nsing)
1010  call hpsort(nsing,nsing,xsing(1:nsing))
1011!write(6,'("sings",i4,4es15.5)') nsing,xsing(1:nsing)
1012  xmin=0.d0
1013  xmax=1.d0
1014  abr=1.d-12
1015  rlr=1.d-6
1016  sint1=grater(xfn4,xmin,xmax,abr,rlr,nsing,xsing,error,numcal,maxns)*area
1017  sint2=grater(xfn5,xmin,xmax,abr,rlr,nsing,xsing,error,numcal,maxns)*area
1018  sint3=grater(xfn6,xmin,xmax,abr,rlr,nsing,xsing,error,numcal,maxns)*area
1019!write(6,*)
1020!write(6,'("sints",6f10.6)') sint1,sint2,sint3
1021  svec=xkp(:,3)*sint1+xq1*(sint2-bb*sint1-cc*sint3)/(2*aa)+xq2*sint3
1022!write(6,'("svec",6f10.6)') svec
1023!tint=(aa0*sint1+dot_product(av,svec))/bmag
1024
1025end subroutine singint2
1026
1027!***************************************************************************
1028
1029double precision function xfn4(xx)
1030implicit none
1031double precision :: xx,xx2
1032double precision :: aa,bb,cc,dd,ee,ff
1033double precision :: xt1,xt2,xl1,qq
1034double precision :: sqq,xterm
1035integer :: iter  !######
1036common /fn/ aa,bb,cc,dd,ee,ff,iter !######
1037
1038  xx2=xx*xx
1039  xt1=bb+cc*xx
1040  xt2=dd*xx+ee*xx2+ff
1041  xl1=2.d0*aa+xt1
1042  qq=4.d0*aa*xt2-xt1**2
1043  if (qq.lt.0.d0) then
1044    sqq=sqrt(-qq)
1045    xfn4=(log(abs((sqq+xl1)/(sqq-xl1)))-log(abs((sqq+xt1)/(sqq-xt1))))/sqq
1046  else
1047    sqq=sqrt(qq)
1048    xfn4=2*(atan((xl1)/sqq)-atan(xt1/sqq))/sqq
1049  endif
1050!write(47,'(es15.5,4x,2es15.5)') xx,xfn4
1051
1052end function xfn4
1053
1054!***************************************************************************
1055
1056double precision function xfn5(xx)
1057implicit none
1058double precision :: xx,xx2
1059double precision :: aa,bb,cc,dd,ee,ff
1060double precision :: xt1,xt2,xl1,qq
1061double precision :: xterm,sqq
1062integer :: iter  !######
1063common /fn/ aa,bb,cc,dd,ee,ff,iter !######
1064
1065  xx2=xx*xx
1066  xterm=(aa+bb+(cc+dd)*xx+ee*xx2+ff)/(dd*xx+ee*xx2+ff)
1067  xfn5=log(abs(xterm))
1068!  xt1=bb+cc*xx
1069!  xt2=dd*xx+ee*xx2+ff
1070!  xl1=2.d0*aa+xt1
1071!  qq=4.d0*aa*xt2-xt1**2
1072!  if (qq.lt.0.d0) then
1073!    sqq=sqrt(-qq)
1074!    xfn5=log(abs(xterm))+xt1*(log(abs((sqq+xl1)/(sqq-xl1)))-log(abs((sqq+xt1)/(sqq-xt1))))/sqq
1075!  else
1076!    sqq=sqrt(qq)
1077!    xfn5=log(abs(xterm))+2*xt1*(atan((xl1)/sqq)-atan(xt1/sqq))/sqq
1078!  endif
1079
1080end function xfn5
1081
1082!***************************************************************************
1083
1084double precision function xfn6(xx)
1085implicit none
1086double precision :: xfn4,xx
1087external xfn4
1088
1089  xfn6=xx*xfn4(xx)
1090
1091end function xfn6
1092
1093