1subroutine mkploss(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,total_se,ploss)
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
43double complex, allocatable :: ssi(:)
44double precision :: ploss
45double precision :: ploss1(-nwpt:nwpt)
46double precision :: total_se
47double complex :: ctest
48integer :: ii,jj,kk,ll,ix,iy,iz,iskip,ocsign,iiq
49integer :: iqq(3),jka(3),jkb(3),jkk(3),iqv(3),ictr(3),iqr(3),ixr(3)
50integer :: ixx(3),ixxp(3),ixp2(3),ixp1(3),ixm1(3),ixm2(3)
51integer :: ikpt,ikptq
52integer :: igg(3),igh(3),igg0(3),isym,isymq,ibp,iqpt,iqsym,iw,iww,ie,je,ies,iqctr
53double precision :: xck(3),xckq(3),qadj(6)
54double complex :: cmatel,cmatel2,amatel,amatel2
55double complex :: jmatel(3),jmatel2(3),vmatel(3),vmatel2(3)
56double complex :: vqmat2(ngkpt(1),ngkpt(2),ngkpt(3))
57double complex :: xmat2(ngkpt(1),ngkpt(2),ngkpt(3))
58double complex :: jmat2(3,3), rjmat2(3,3), rjmatdum
59double precision :: smat2
60double complex :: vfactor(ngkpt(1),ngkpt(2),ngkpt(3),nwpt),vjfactor(3,3,nwpt),vsfactor(3,3,nwpt)
61double complex :: vfv(8,nwpt),vfx(8),vsm1(3,3,nwpt),vsm2(3,3,nwpt)
62double precision :: vq2(8),vqvtx0(4),vqvtx(4),qkcvt(3,3)
63double precision :: omega(ngkpt(1),ngkpt(2),ngkpt(3)), omega_t
64double precision :: whi,wlo,wwhi,wwlo,ww,www,wwwlo,wwwhi,enval(8),dw,rw,eshift,evtx0(4),evtx(4)
65double precision :: rrpyr(3,4),kvtx(3,4),xk(3,4),ckvtx(3,4),cxk(3,4),rg(3,3),tvol,xpyr0(4),xpyr(4)
66double precision :: xkc(3,4),rgc(3,3),tvolc,kvtxc(3,4)
67double complex :: vpyr0(4,nwpt),vpyr(4)
68double precision :: de21,de31,de32,de41,de42,de43,thresh
69double precision :: fbx(4),fb(4),cmx(3,4),cm(3,4),xkt(3,3)
70double complex :: aa0(nwpt),av(3,nwpt)
71double precision :: xme,xv(3)
72double precision :: avec(3),bgrad(3),xmult
73double precision :: abr,rlr
74integer :: iwh,iwl,ibmin,ibmax,iwhi,iwlo,iehi,ielo
75integer :: indxe(4),iwwlim(4)
76double precision :: vq(ngkpt(1),ngkpt(2),ngkpt(3)),qq(3),qp(3),qq2,qp2,qs(3),ek,ekmq
77double precision :: stvec(ngkpt(3)),stvec2(ngkpt(3)),temparray(ngkpt(1),ngkpt(2),ngkpt(3))
78integer :: iipw,jjpw,jjpwt
79integer :: iqsymndx(ngkpt(1),ngkpt(2),ngkpt(3)),invpw2ndx(npwx,npwx)
80integer :: pwsymndx(npwc,2*nsym)
81integer :: ipw2,jpw1,jpw2,jw,iqcentr,jqcentr
82double precision :: eps1,eps2,eps
83double precision :: eval,brd,esprd
84double precision :: gfo,gamma(3),pln(3),dist(3)
85double precision :: sint1,sint1a,sint1b,svec(3),sveca(3),svecb(3),stens(3,3)
86double complex :: cint
87double complex :: wint,wint0,w2int,w2int0,gterm,vcentr
88double precision :: wcentr(-nwpt:nwpt),wgrid(nwpt)
89integer :: ntypepaw
90double complex :: pwmatel(ntypepaw,nlmn,nlmn,npwx,ngkpt(1),ngkpt(2),ngkpt(3)), &
91&                tpwmatel(ntypepaw,nlmn,nlmn,npwx,ngkpt(1),ngkpt(2),ngkpt(3)), &
92&                 pwjmatel(3,ntypepaw,nlmn,nlmn),               &
93&                tpwjmatel(3,ntypepaw,nlmn,nlmn)
94logical :: lqsing,lqcentr
95logical :: lc
96double precision :: linterp
97double precision :: rr(3,8)
98integer :: ivndx(4,6),itet,iv
99integer :: nsing
100double precision :: wsing(20)
101character*4 :: label
102external linterp
103data rr(1:3,1) /0.0d0,0.0d0,0.0d0/
104data rr(1:3,2) /1.0d0,0.0d0,0.0d0/
105data rr(1:3,3) /0.0d0,1.0d0,0.0d0/
106data rr(1:3,4) /1.0d0,1.0d0,0.0d0/
107data rr(1:3,5) /0.0d0,0.0d0,1.0d0/
108data rr(1:3,6) /1.0d0,0.0d0,1.0d0/
109data rr(1:3,7) /0.0d0,1.0d0,1.0d0/
110data rr(1:3,8) /1.0d0,1.0d0,1.0d0/
111data ivndx(1:4,1) /1,2,3,5/
112data ivndx(1:4,2) /3,5,6,7/
113data ivndx(1:4,3) /2,3,5,6/
114data ivndx(1:4,4) /2,3,4,6/
115data ivndx(1:4,5) /3,4,6,7/
116data ivndx(1:4,6) /4,6,7,8/
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)
130
131!write(6,*) "qopf ",total_se,ploss
132i_prt_DEBUG = 0
133j_prt_DEBUG = 0
134k_prt_DEBUG = 0
135i_ct_DEBUG = 0
136j_ct_DEBUG = 0
137k_ct_DEBUG = 0
138l_ct_DEBUG = 0
139i_ref_DEBUG = -1
140k_ref_DEBUG = -1
141abr=1.d-4
142rlr=1.d-4
143sei=0.d0
144ctest=(0.d0,0.d0)
145igg0=(/0,0,0/)
146ploss=0.d0
147do ie=-nwpt,nwpt
148  ploss1(ie)=0.d0
149enddo
150dw=wmax/dble(nwpt)
151do iw=1,nwpt
152  wgrid(iw)=iw*dw
153enddo
154volelmnt=(vol*ngkpt(1)*ngkpt(2)*ngkpt(3))/((2.d0*pi)**3)
155volred = 1.d0
156ikpt=ikndx(ikk(1),ikk(2),ikk(3))
157isym=isymndx(ikk(1),ikk(2),ikk(3))
158igsymk=isymg(:,ikpt,isym)
159ek=enrgy(indxkbnd(ikpt)+iband)
160do ii=1,3
161do jj=1,3
162  qkcvt(ii,jj)=blat(ii,jj)/dble(ngkpt(ii))
163enddo
164enddo
165ictr=(/ngkpt(1)/2,ngkpt(2)/2,ngkpt(3)/2/)
166iqctr=iqndx(ictr(1),ictr(2),ictr(3))
167nsing=0
168
169do ii=1,3
170do kk=1,3
171  bmet2(ii,kk) = 0.d0                   ! bmet2 = bmet*bmet
172  do jj=1,3
173    bmet2(ii,kk) = bmet2(ii,kk) + bmet(ii,jj)*bmet(jj,kk)
174  enddo
175enddo
176enddo
177
178do ix=1,ngkpt(1)
179do iy=1,ngkpt(2)
180do iz=1,ngkpt(3)
181  temparray(ix,iy,iz)=(0.d0,0.d0)
182enddo
183enddo
184enddo
185
186gamma=(blat(1,:)+blat(2,:)+blat(3,:))/2.d0
187do ii=1,3
188  jj=mod(ii,3)+1
189  kk=mod(ii+1,3)+1
190  pln(1)=blat(jj,2)*blat(kk,3)-blat(jj,3)*blat(kk,2)
191  pln(2)=-blat(jj,1)*blat(kk,3)+blat(jj,3)*blat(kk,1)
192  pln(3)=blat(jj,1)*blat(kk,2)-blat(jj,2)*blat(kk,1)
193  dist(ii)=dot_product(gamma,pln)/sqrt(dot_product(pln,pln))
194enddo
195gfo=minval(dist)
196
197!ipaw=0
198!write(6,'(a)') '      finding contibutions from:'
199!write(6,'(a)') 'band  plane waves    correlation                exchange'
200do ibp=test_bands_se(1),test_bands_se(2)
201!do ibp=2,2
202!do ibp=iband,iband
203  write(6,'(3(a,i4))') 'quasiparticle band ',iband,', polarized band ',ibp,' out of ',test_bands_se(2)
204  seibc=sei
205  if (ibp.gt.nbocc) then
206    ocsign=-1
207  else
208    ocsign=1
209  endif
210  do iipw=1,napwndx
211!  do iipw=napwndx,napwndx
212!  do iipw=1,1
213!  do iipw=15,15
214!    write(6,'(a,i4)') 'iipw ',iipw
215    seipw=sei
216    ipw1=ipwndx(1,iipw)
217    ipw2=ipwndx(2,iipw)
218!    if (ipw1.ne.ipw2) cycle
219!    write(6,'(1x,i3,4x,2i4)') ibp,ipw1,ipw2
220    lc=iipw.le.ntpwndx  ! do correlation
221    if (.not.lc) cycle
222
223! Step 1: Compute matrix elements
224    do ix=1,ngkpt(1)
225    do iy=1,ngkpt(2)
226    do iz=1,ngkpt(3)
227!    do ix=1,1
228!    do iy=10,10
229!    do iz=10,10
230      vqmat2(ix,iy,iz)=(0.d0,0.d0)
231      ixx=(/ix,iy,iz/)
232      iqq=ixx-ictr
233      qq = dble(iqq)/dble(ngkpt)
234      jka=ikk-iqq
235      jkk=mod(jka,ngkpt)
236      do ii=1,3
237        if (jkk(ii).le.0) jkk(ii)=jkk(ii)+ngkpt(ii)
238      enddo
239      ikptq=ikndx(jkk(1),jkk(2),jkk(3))
240      isymq=isymndx(jkk(1),jkk(2),jkk(3))
241      igsymq=isymg(:,ikptq,isymq)
242      jkb=nint(kpt(:,ikptq)*dble(ngkpt))+ictr
243      igg=nint(dble(jkk-jka)/dble(ngkpt))
244!write(6,*) '>>>>> iqq = ',iqq
245!write(6,*) '>>>>> ikk = ',ikk
246!write(6,*) '>>>>> jkk = ',jkk
247!write(6,*) '>>>>> jka = ',jka
248!write(6,*) '>>>>> igg = ',igg
249!write(6,'(a,3f10.5)') 'kpt = ',kpt(:,ikpt)
250!write(6,'(a,3f10.5)') 'kptq = ',kpt(:,ikptq)
251!read(*,*)
252      qp=qq+ipwx(:,ipw2)
253      qq=qq+ipwx(:,ipw1)
254      qq2=0.d0
255      qp2=0.d0
256      do ii=1,3
257      do jj=1,3
258        qq2=qq2+qq(ii)*bmet(ii,jj)*qq(jj)
259        qp2=qp2+qp(ii)*bmet(ii,jj)*qp(jj)
260      enddo
261      enddo
262      vq(ix,iy,iz)=4.d0*pi/sqrt(qq2*qp2)
263!write(6,'(a,3f10.5)') 'qq = ',qq
264!write(6,'(a,3f10.5)') 'qp = ',qp
265!write(6,'(a,3f10.5)') 'qq2 = ',qq2
266!write(6,'(a,3f10.5)') 'qp2 = ',qp2
267!write(6,*) 'vq = ',vq(ix,iy,iz)
268!write(6,*) 'jka = ',jka
269!write(6,*) 'jkk = ',jkk
270!write(6,*) 'jkb = ',jkb
271!write(6,*) 'igg = ',igg
272!write(6,*) 'ikpt = ',ikpt
273!write(6,*) 'ikptq = ',ikptq
274      omega_t=enrgy(indxkbnd(ikptq)+ibp)-ek
275      omega(ix,iy,iz)=enrgy(indxkbnd(ikptq)+ibp)-ek+dw/2
276      if (iqq(1).eq.0.and.iqq(2).eq.0.and.iqq(3).eq.0) then
277        if (ipw1.eq.1.and.ipw2.eq.1.and.iband.eq.ibp) then
278! Find singular factor
279!          call mkmatelX1(ibp,iband,ikpt,ikptq,ipwx(:,ipw1),igg, &
280!&           ncg,nkpt,npwt,igmx,igmn,igndx, &
281!&           isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, &
282!&           lvtrans(1:3,ikk(1),ikk(2),ikk(3)), &
283!&           cg,indxkcg,indxkpw,npwarr,kg, &
284!&           cmatel)
285!          if (ipaw.ne.0) then
286!            call mkPAWmatelX1(pi,ibp,iband,ikk,jkk,qq,igg0, &
287!&                 pwmatel,tpwmatel,nbcore,ncband,amatel)
288!          else
289!            amatel2=(0.d0,0.d0)
290!          endif
291!          smat2=4.d0*pi*dble((cmatel+amatel)*conjg(cmatel+amatel))
292!          if (smat2.ne.0.d0) then
293!            lqsing=.true.
294!          else
295!            lqsing=.false.
296!          endif
297          smat2=4.d0*pi*omega_t
298          lqsing=.true.
299        else
300          smat2=0.d0
301          lqsing=.false.
302        endif
303! If not singular, find tensor response at q=0
304        if (.not.lqsing) then
305          if (ipwx(1,ipw1).eq.0.and.ipwx(2,ipw1).eq.0.and.ipwx(3,ipw1).eq.0) then
306            call mkmatelJ1(ibp,iband,ikpt,ikptq,ipwx(:,ipw1),igg, &
307&             ncg,nkpt,npwt,igmx,igmn,igndx, &
308&             isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, &
309&             lvtrans(1:3,ikk(1),ikk(2),ikk(3)), &
310&             cg,indxkcg,indxkpw,npwarr,kg, &
311&             jmatel)
312            if (ipaw.ne.0) then
313              call mkPAWmatelJ1(pi,ibp,iband,ikk,jkk, &
314&                           pwjmatel,tpwjmatel,nbcore,ncband,vmatel)
315              jmatel(:) = jmatel(:) + vmatel(:)
316            endif
317            do ii=1,3
318              jmatel(ii) = jmatel(ii)/(enrgy(indxkbnd(ikptq)+ibp)-ek)
319            enddo
320          else
321            call mkmatelX1(ibp,iband,ikpt,ikptq,ipwx(:,ipw1),igg, &
322&             ncg,nkpt,npwt,igmx,igmn,igndx, &
323&             isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, &
324&             lvtrans(1:3,ikk(1),ikk(2),ikk(3)), &
325&             cg,indxkcg,indxkpw,npwarr,kg, &
326&             cmatel)
327            if (ipaw.ne.0) then
328              call mkPAWmatelX1(pi,ibp,iband,ikk,jkk,qq,igg0, &
329&                   pwmatel,tpwmatel,nbcore,ncband,amatel)
330            else
331              amatel=(0.d0,0.d0)
332            endif
333            jmatel = (cmatel+amatel)*qq/qq2
334          endif
335          if (ipw1.eq.ipw2) then
336            jmatel2=jmatel
337          else
338            if (ipwx(1,ipw2).eq.0.and.ipwx(2,ipw2).eq.0.and.ipwx(3,ipw2).eq.0) then
339              call mkmatelJ1(ibp,iband,ikpt,ikptq,ipwx(:,ipw2),igh, &
340&               ncg,nkpt,npwt,igmx,igmn,igndx, &
341&               isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, &
342&               lvtrans(1:3,ikk(1),ikk(2),ikk(3)), &
343&               cg,indxkcg,indxkpw,npwarr,kg, &
344&               jmatel2)
345              if (ipaw.ne.0) then
346                call mkPAWmatelJ1(pi,ibp,iband,ikk,jkk, &
347&                             pwjmatel,tpwjmatel,nbcore,ncband,vmatel)
348                jmatel2(:) = jmatel2(:) + vmatel(:)
349              endif
350              do ii=1,3
351                jmatel2(ii) = jmatel2(ii)/(enrgy(indxkbnd(ikptq)+ibp)-ek)
352              enddo
353            else
354              call mkmatelX1(ibp,iband,ikpt,ikptq,ipwx(:,ipw2),igh, &
355&             ncg,nkpt,npwt,igmx,igmn,igndx, &
356&             isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, &
357&             lvtrans(1:3,ikk(1),ikk(2),ikk(3)), &
358&             cg,indxkcg,indxkpw,npwarr,kg, &
359&             cmatel2)
360              if (ipaw.ne.0) then
361                call mkPAWmatelX1(pi,ibp,iband,ikk,jkk,qp,igg0, &
362&                     pwmatel,tpwmatel,nbcore,ncband,amatel2)
363              else
364                amatel2=(0.d0,0.d0)
365              endif
366              jmatel2 = (cmatel2+amatel2)*qp/qp2
367            endif
368          endif
369!! DEBUG
370!write(6,*) "jmatel"
371!do ii=1,3
372!write(6,*) jmatel(ii)
373!enddo
374!write(6,*) "jmatel2"
375!do ii=1,3
376!write(6,*) jmatel2(ii)
377!enddo
378!read(*,*)
379          do ii=1,3
380          do jj=1,3
381            jmat2(ii,jj)=4.d0*pi*jmatel(ii)*conjg(jmatel2(jj))*omega_t
382          enddo
383          enddo
384        else
385          do ii=1,3
386          do jj=1,3
387            jmat2(ii,jj)=(0.d0,0.d0)
388          enddo
389          enddo
390        endif
391      else
392        call mkmatelX1(ibp,iband,ikpt,ikptq,ipwx(:,ipw1),igg, &
393&         ncg,nkpt,npwt,igmx,igmn,igndx, &
394&         isym,isymq,symrel,syminv,nsym,ihlf,igsymk,igsymq,kpt, &
395&         lvtrans(1:3,ikk(1),ikk(2),ikk(3)), &
396&         cg,indxkcg,indxkpw,npwarr,kg, &
397&         cmatel)
398        if (ipaw.ne.0) then
399          call mkPAWmatelX1(pi,ibp,iband,ikk,jkk,qq,igg0, &
400&               pwmatel,tpwmatel,nbcore,ncband,amatel)
401        else
402          amatel=(0.d0,0.d0)
403        endif
404        if (ipw1.eq.ipw2) then
405          cmatel2=cmatel
406          amatel2=amatel
407        else
408          call mkmatelX1(ibp,iband,ikpt,ikptq,ipwx(:,ipw2),igh, &
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&         cmatel2)
414          if (ipaw.ne.0) then
415            call mkPAWmatelX1(pi,ibp,iband,ikk,jkk,qp,igg0, &
416&                 pwmatel,tpwmatel,nbcore,ncband,amatel2)
417          else
418            amatel2=(0.d0,0.d0)
419          endif
420        endif
421        xmat2(ix,iy,iz)=(cmatel+amatel)*conjg(cmatel2+amatel2)*omega_t
422        vqmat2(ix,iy,iz)=xmat2(ix,iy,iz)*vq(ix,iy,iz)
423      endif
424    enddo
425    enddo
426    enddo
427
428! Alternate way of getting jmat2
429! Interpolation of density matrix elements / q^2 to q=0
430! For equally spaced samples A(-2) = -2 q, A(-1) = -1 q, A(1) = 1 q, A(2) = 2 q
431! average of two quadratic interpolations gives (see Numerical Recipes 3.1)
432! A(0) = (2/3) A(-1) + (2/3) A(1) - (1/6) A(-2) - (1/6) A(2)
433! Expect vqmat2 to converge to well defined value at q = 0, interpolate this
434! sum(ii,jj) qq(ii) * (bmet * jmat2 * bmet)(ii,jj) * qq(jj) / qq^2 = vqmat
435! define rjmat2 = bmet * jmat2 * bmet
436! find rjmat2 as q->0, then
437! jmat2(q->0) = bmet^(-1)*rjmat2(q->0)*bmet^(-1)
438!
439! 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)
440! so longitudinal response only depends on sum of symmetric off-diagonal terms,
441! independent of difference.
442! Thus, we can safely set rjmat2(ii,jj) = (rjmat2(ii,jj)+rjmat2(jj,ii))/2
443    if (.true.) then
444! first do diagonal parts
445! e.g., along direction 1
446! qq = qq(1) * blat(1,:)
447! qq^2 = qq(1)*qq(1) * dot_product(blat(1,:),blat(1,:))
448! sum(ii,jj) qq(ii) * rjmat2(ii,jj) * qq(jj) / qq^2
449! = (1/qq^2) * rjmat2(1,1) * qq(1)*qq(1)
450! = rjmat2(1,1) / dot_product(blat(1,:),blat(1,:))
451! -> rjmat2(1,1) = dot_product(blat(1,:),blat(1,:)) * [response in direction 1]
452      do ii=1,3
453        ixx = (/0,0,0/)
454        ixx(ii) = 1
455        ixp2 = ictr+2*ixx
456        ixp1 = ictr+ixx
457        ixm1 = ictr-ixx
458        ixm2 = ictr-2*ixx
459        qq2 = dot_product(blat(ii,:),blat(ii,:))
460        rjmatdum = 2.d0/3.d0 * vqmat2(ixp1(1),ixp1(2),ixp1(3)) &
461&                + 2.d0/3.d0 * vqmat2(ixm1(1),ixm1(2),ixm1(3)) &
462&                - 1.d0/6.d0 * vqmat2(ixp2(1),ixp2(2),ixp2(3)) &
463&                - 1.d0/6.d0 * vqmat2(ixm2(1),ixm2(2),ixm2(3))
464        rjmat2(ii,ii) =  qq2 * rjmatdum
465      enddo
466! find off-diagonals using qq along directions (1,1,0), (1,0,1), and (0,1,1)
467! e.g., along direction (1,1,0)
468! qq = qq(1) * blat(1,:) + qq(2) * blat(2,:)
469! qq^2 =    qq(1)*qq(1) * dot_product(blat(1,:),blat(1,:))
470!       +   qq(2)*qq(2) * dot_product(blat(2,:),blat(2,:))
471!       + 2*qq(1)*qq(2) * dot_product(blat(1,:),blat(2,:))
472!      =    dot_product(blat(1,:),blat(1,:)) + dot_product(blat(2,:),blat(2,:))
473!       + 2*dot_product(blat(1,:),blat(2,:))
474! sum(ii,jj) qq(ii) * rjmat2(ii,jj) * qq(jj) / qq^2
475! = (1/qq^2) * ( rjmat2(1,1) * qq(1)*qq(1) + rjmat2(2,2) * qq(2)*qq(2)
476!               + qq(1)*qq(2) * (rjmat2(1,2) + rjmat2(2,1)) )
477! = (1/qq^2) * ( rjmat2(1,1) * qq(1)*qq(1) + rjmat2(2,2) * qq(2)*qq(2)
478!               + 2 * rjmat2(1,2) qq(1)*qq(2))
479! = (1/qq^2) * ( rjmat2(1,1) + rjmat2(2,2) + 2 * rjmat2(1,2) )
480! -> rjmat2(1,2) = (qq^2 * [response in direction 1] - (rjmat2(1,1) + rjmat2(2,2)))/2
481      do ii=1,3
482        ixx = (/1,1,1/)
483        ixx(ii) = 0
484        ixp2 = ictr+2*ixx
485        ixp1 = ictr+ixx
486        ixm1 = ictr-ixx
487        ixm2 = ictr-2*ixx
488        qq2 = dot_product(blat(jj,:),blat(jj,:))   &
489&           + dot_product(blat(kk,:),blat(kk,:))   &
490&           + dot_product(blat(jj,:),blat(kk,:))*2
491        rjmatdum = 2.d0/3.d0 * vqmat2(ixp1(1),ixp1(2),ixp1(3)) &
492&                + 2.d0/3.d0 * vqmat2(ixm1(1),ixm1(2),ixm1(3)) &
493&                - 1.d0/6.d0 * vqmat2(ixp2(1),ixp2(2),ixp2(3)) &
494&                - 1.d0/6.d0 * vqmat2(ixm2(1),ixm2(2),ixm2(3))
495        jj = mod(ii,3)+1   ! note ii = mod(ii-1,3)+1
496        kk = mod(ii+1,3)+1
497        rjmat2(jj,kk) = (qq2*rjmatdum - (rjmat2(jj,jj)+rjmat2(kk,kk)))/2.d0
498        rjmat2(kk,jj) = rjmat2(jj,kk)
499!! DEBUG
500!qq = ixx/dble(ngkpt)
501!qq2=0.d0
502!do kk=1,3
503!do jj=1,3
504!qq2=qq2+qq(kk)*bmet(kk,jj)*qq(jj)
505!enddo
506!enddo
507!cvecdum = bmet(:,1)*qq(1)+bmet(:,2)*qq(2)+bmet(:,3)*qq(3)
508!cvecdum2 = jmat2(:,1)*cvecdum(1)+jmat2(:,2)*cvecdum(2)+jmat2(:,3)*cvecdum(3)
509!cvecdum = bmet(:,1)*cvecdum2(1)+bmet(:,2)*cvecdum2(2)+bmet(:,3)*cvecdum2(3)
510!write(6,*) "qq: ",qq
511!write(6,*) (cvecdum(1)*qq(1)+cvecdum(2)*qq(2)+cvecdum(3)*qq(3))/qq2
512!write(6,*) ixx
513!write(6,*) vqmat2(ixp2(1),ixp2(2),ixp2(3))
514!write(6,*) vqmat2(ixp1(1),ixp1(2),ixp1(3))
515!write(6,*) rjmatdum
516!write(6,*) vqmat2(ixm1(1),ixm1(2),ixm1(3))
517!write(6,*) vqmat2(ixm2(1),ixm2(2),ixm2(3))
518      enddo
519! DEBUG
520!do ii=1,3
521!do jj=1,3
522!if (ii.eq.jj) then
523!rjmat2(ii,jj) = 2.9357115896d0
524!else
525!rjmat2(ii,jj) = 0.d0;
526!endif
527!enddo
528!enddo
529      bmet_t=bmet
530      call inverse(bmet_t,bmetinv,3,3)
531      do ii=1,3
532      do jj=1,3
533        jmat2(ii,jj) = (0.d0,0.d0)
534        do kk=1,3
535        do ll=1,3
536          jmat2(ii,jj) = jmat2(ii,jj) + bmetinv(ii,kk)*rjmat2(kk,ll)*bmetinv(ll,jj)
537        enddo
538        enddo
539      enddo
540      enddo
541    endif
542!! DEBUG
543!write(6,*) "rjmat2"
544!do ii=1,3
545!write(6,*) dble(rjmat2(:,ii))
546!enddo
547!write(6,*)
548!do ii=1,3
549!write(6,*) dimag(rjmat2(:,ii))
550!enddo
551!write(6,*)
552!write(6,*) "jmat2"
553!do ii=1,3
554!write(6,*) dble(jmat2(:,ii))
555!enddo
556!write(6,*)
557!do ii=1,3
558!write(6,*) dimag(jmat2(:,ii))
559!enddo
560!write(6,*)
561!write(6,*) "in Cartesian coordinates"
562!do ii=1,3
563!do jj=1,3
564!ctensdum(ii,jj) = (0.d0,0.d0)
565!do kk=1,3
566!do ll=1,3
567!ctensdum(ii,jj) = ctensdum(ii,jj) + blat(kk,ii)*jmat2(kk,ll)*blat(ll,jj)
568!enddo
569!enddo
570!enddo
571!enddo
572!do ii=1,3
573!write(6,*) dble(ctensdum(:,ii))
574!enddo
575!write(6,*)
576!do ii=1,3
577!write(6,*) dimag(ctensdum(:,ii))
578!enddo
579!stop
580! DEBUG
581!write(6,*) "Matrix elements computed"
582
583    if (itetrahedron.eq.0) then
584      write(6,*) "Sum over points not yet implemented in self energy calculation"
585      write(6,*) "set itetrahedron=1 and run again"
586    else
587! tetrahedron integration
588
589! Don't waste time integrating over energies with no contribution
590      do ix=1,ngkpt(1)
591      do iy=1,ngkpt(2)
592      do iz=1,ngkpt(3)
593        ixx=(/ix,iy,iz/)
594        call fhilo(omega,ixx,ngkpt,whi,wlo)
595        if (ix.eq.1.and.iy.eq.1.and.iz.eq.1) then
596          wwhi=whi
597          wwlo=wlo
598        else
599          wwhi=max(wwhi,whi)
600          wwlo=min(wwlo,wlo)
601        endif
602      enddo
603      enddo
604      enddo
605
606      iwlo=nint(wwlo*dble(nwpt)/wmax)
607      iwhi=nint(wwhi*dble(nwpt)/wmax)
608      if (ibp.gt.nbocc) then
609        ielo=iwlo+1
610        iehi=iwhi+nwpt
611      else
612        ielo=iwlo-nwpt
613        iehi=iwhi-1
614      endif
615      allocate (ssi(ielo:iehi))
616      do ie=ielo,iehi
617        ssi(ie)=0.d0
618      enddo
619
620      if (lc) then
621        do iw=1,nwpt
622          do ix=1,ngkpt(1)
623          do iy=1,ngkpt(2)
624          do iz=1,ngkpt(3)
625            ixx=(/ix,iy,iz/)
626            iqpt=iqndx(ixx(1),ixx(2),ixx(3))
627            call locateelement(ixx,ipw1,ipw2,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpw)
628            call locateelement(ixx,ipw2,ipw1,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpwt) ! transpose matrix coordiantes, to find anti-Hermitian part of lossfn
629            if (iqpt.eq.iqctr) then
630              vfactor(ixx(1),ixx(2),ixx(3),iw)=(0.d0,0.d0)
631            else if (jjpw.ne.0.and.jjpwt.ne.0) then
632              vfactor(ixx(1),ixx(2),ixx(3),iw)= &
633&                 vqmat2(ixx(1),ixx(2),ixx(3))* &
634&                 (lossfn(iw,iqpt,jjpw)-conjg(lossfn(iw,iqpt,jjpwt)))*0.5d0
635            else
636              vfactor(ixx(1),ixx(2),ixx(3),iw)=(0.d0,0.d0)
637            endif
638          enddo
639          enddo
640          enddo
641!write(6,*) "iw = ", iw
642          call locateelement(ictr,ipw1,ipw2,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpw)
643          call locateelement(ictr,ipw2,ipw1,ngkpt,iqsymndx,npwc,npwx,nsym,pwsymndx,invpw2ndx,jjpwt)
644          do ii=1,3
645          do jj=1,3
646            if (jjpw.ne.0.and.jjpwt.ne.0) then
647              vjfactor(ii,jj,iw)=(0.d0,0.d0)
648              do kk=1,3
649              do ll=1,3
650                iiq=ii + 3*(kk-1)
651                vjfactor(ii,jj,iw)= &
652&                     vjfactor(ii,jj,iw) + &
653&                     0.5d0*(lossfn(iw,nqpt+iiq,jjpw)-conjg(lossfn(iw,nqpt+iiq,jjpwt))) &
654&                     *bmet(kk,ll)*jmat2(ll,jj)
655              enddo
656              enddo
657              iiq=ii + 3*(jj-1)
658              vsfactor(ii,jj,iw)=smat2*(lossfn(iw,nqpt+iiq,jjpw)-conjg(lossfn(iw,nqpt+iiq,jjpwt)))*0.5d0
659            else
660              vjfactor(ii,jj,iw)=(0.d0,0.d0)
661              vsfactor(ii,jj,iw)=(0.d0,0.d0)
662            endif
663          enddo
664          enddo
665!do ii=1,3
666!write(6,*) dble(vjfactor(ii,:,iw))
667!enddo
668!write(6,*)
669!do ii=1,3
670!write(6,*) dimag(vjfactor(ii,:,iw))
671!enddo
672!write(6,*)
673!do ii=1,3
674!write(6,*) dble(vsfactor(ii,:,iw))
675!enddo
676!write(6,*)
677!do ii=1,3
678!write(6,*) dimag(vsfactor(ii,:,iw))
679!enddo
680!read(*,*)
681          do ii=1,3
682          do jj=1,3
683            vsm1(ii,jj,iw)=(0.d0,0.d0)
684            do kk=1,3
685              vsm1(ii,jj,iw)=vsm1(ii,jj,iw)+bmet(ii,kk)*vsfactor(kk,jj,iw)
686            enddo
687          enddo
688          enddo
689          do ii=1,3
690          do jj=1,3
691            vsm2(ii,jj,iw)=(0.d0,0.d0)
692            do kk=1,3
693              vsm2(ii,jj,iw)=vsm2(ii,jj,iw)+vsm1(ii,kk,iw)*bmet(kk,jj)
694            enddo
695          enddo
696          enddo
697        enddo
698      endif
699
700! DEBUG
701!write(6,*) "Beginning integration"
702
703! Do the integration
704      do ix=1,ngkpt(1)
705      do iy=1,ngkpt(2)
706      do iz=1,ngkpt(3)
707        ixx=(/ix,iy,iz/)
708        iqq=ixx-ictr
709        call fval(omega,ixx,ixxp,ngkpt,enval)
710        if (lqsing) call fval(vq,ixx,ixxp,ngkpt,vq2)
711        if (lc) then
712          do iw=1,nwpt
713            call fpol(vfactor(:,:,:,iw),ngkpt,ixx,ixxp,vfv(:,iw))
714          enddo
715        endif
716        do itet=1,6
717          centerint = .false.
718          do iv=1,4
719            evtx0(iv)=enval(ivndx(iv,itet))
720            rrpyr(1:3,iv)=rr(1:3,ivndx(iv,itet))
721            if (lc) then
722              do iw=1,nwpt
723                vpyr0(iv,iw)=vfv(ivndx(iv,itet),iw)
724              enddo
725            endif
726            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
727              centerint = .true.  ! does tetrahedron contain center vertex?
728              vcenterint = iv     ! index of center vertex
729            endif
730          enddo
731!write(6,*) ixx
732!write(6,*) itet,centerint
733!if (centerint) then
734!write(6,*) itet
735!do iv=1,4
736!write(6,*) ixx+rrpyr(:,iv)
737!enddo
738!!read(*,*)
739!endif
740          call indxhpsort(4,4,evtx0,indxe)
741          evtx=evtx0(indxe)
742          do ii=1,4
743            jj=indxe(ii)
744            kvtx(:,ii) = (iqq(:) + rrpyr(:,jj))
745            do kk=1,3
746              kvtxc(kk,ii)=dot_product(iqq+rrpyr(:,jj),qkcvt(:,kk))
747            enddo
748          enddo
749          xk(1:3,1)=(/0.d0,0.d0,0.d0/)
750          xkc(1:3,1)=(/0.d0,0.d0,0.d0/)
751          do ii=1,3
752            xk(:,ii+1)=rrpyr(:,indxe(ii+1))-rrpyr(:,indxe(1))
753          enddo
754          do ii=2,4
755            xkc(1:3,ii)=kvtxc(1:3,ii)-kvtxc(1:3,1)
756          enddo
757          do jj=1,3  ! make contragredient, reduced coordinates
758            call cross(xk(:,mod(jj,3)+2),xk(:,mod(jj+1,3)+2),rg(:,jj))
759          enddo
760          do jj=1,3  ! make contragredient, Cartesian coordinates
761            call cross(xkc(:,mod(jj,3)+2),xkc(:,mod(jj+1,3)+2),rgc(:,jj))
762          enddo
763          tvol=dot_product(xk(1:3,2),rg(1:3,1))
764          rg=rg/tvol
765          tvol=abs(tvol)
766          tvolc=dot_product(xkc(1:3,2),rgc(1:3,1))
767          rgc=rgc/tvolc
768          tvolc=abs(tvolc)
769          iwwlim=nint(evtx/dw)
770          if (lqsing .and. centerint) then
771! singular part
772! DEBUG
773!xdum = 0.d0
774!write(6,*)
775!write(6,*) "singular integration"
776!write(6,*) 'ixx: ',ixx
777!write(6,*) 'itet: ',itet,'  centerint: ',centerint
778!write(6,*) 'iwwlim: ',iwwlim
779            do iv=1,4
780! subtraction of integrated singular function from values at vertices
781              ixr=ixx+rrpyr(:,iv)
782              iqr = ixr - ictr
783              if (ixr(1).eq.ictr(1).and.ixr(2).eq.ictr(2).and.ixr(3).eq.ictr(3)) then
784                if (lc) then
785                  do iw=1,nwpt
786                    vpyr0(iv,iw) = 0.d0
787                  enddo
788                endif
789              else
790                qq = dble(iqr)/dble(ngkpt)
791                qq2 = 0.d0
792                do ii=1,3
793                do jj=1,3
794                  qq2 = qq2 + qq(ii)*bmet(ii,jj)*qq(jj)
795                enddo
796                enddo
797                if (lc) then
798                  do iw=1,nwpt
799                    do ii=1,3
800                    do jj=1,3
801                      vpyr0(iv,iw)=vpyr0(iv,iw)-qq(ii)*vsm2(ii,jj,iw)*qq(jj)/(qq2**2)
802                    enddo
803                    enddo
804                  enddo
805                endif
806              endif
807            enddo
808            bgrad=(/0.d0,0.d0,0.d0/)   ! energy gradient
809            do jj=1,3
810              bgrad=bgrad+(evtx(jj+1)-evtx(1))*rgc(:,jj)
811            enddo
812            xmult=1.d0/sqrt(dot_product(bgrad,bgrad))
813! DEBUG
814!cdum = 0.d0
815! DEBUG
816!write(6,*) "first set"
817            do iww=iwwlim(1),iwwlim(2) ! energies between evtx(1) & evtx(2)
818! DEBUG
819!write(6,*) iww
820              wwwlo=max(dw*(iww-0.5d0),evtx(1))
821              wwwhi=min(dw*(iww+0.5d0),evtx(2))
822              if (abs(wwwhi-wwwlo)*1.d14.lt.max(abs(wwwhi),abs(wwwlo),dw)) cycle
823              do ii=1,3
824              do jj=1,3
825! DEBUG
826k_ct_DEBUG = k_ct_DEBUG + 1
827!write(6,*) "k_ct_DEBUG: ",k_ct_DEBUG,ii,jj
828!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "E_LIM:      ",wwwlo,wwwhi
829!if (k_ct_DEBUG.le.9) write(77,*) "# ",ii,jj
830                call tsinggrater(ii,jj,wwwlo,wwwhi,1,xk,kvtx(:,1),ngkpt,evtx,abr,rlr,nsing,wsing,stens(ii,jj))
831              enddo
832              enddo
833              stens = stens * xmult
834              if (lc) then
835                do iw=1,nwpt
836                  ie=iww-ocsign*iw
837!! DEBUG
838!cdum = 0.d0
839!do ii=1,3
840!do jj=1,3
841!ctensdum2(ii,jj) = (0.d0,0.d0)
842!do kk=1,3
843!do ll=1,3
844!ctensdum2(ii,jj) = ctensdum2(ii,jj) + blat(kk,ii)*vsfactor(kk,ll,iw)*blat(ll,jj)
845!enddo
846!enddo
847!enddo
848!enddo
849                  do ii=1,3
850                  do jj=1,3
851                    ssi(ie)=ssi(ie)+vsm2(ii,jj,iw)*stens(jj,ii)
852!! DEBUG
853!if (ie.eq.0.or.ie.eq.1) then
854!cdum = cdum+vsm2(ii,jj,iw)*stens(jj,ii)
855!endif
856!if (iw.eq.158) then
857!write(6,*) ii,jj,dimag(ctensdum2(ii,jj)),dimag(vsfactor(ii,jj,iw)),dimag(vsm2(ii,jj,iw)),stens(jj,ii)
858!endif
859                  enddo
860                  enddo
861! DEBUG
862!if (iw.eq.158) then
863!write(6,'(i4,23e14.3)') iw,cdum
864!stop
865!endif
866!if (cdum.ne.(0.d0,0.d0)) then
867!write(6,'(i4,23e14.3)') iw,cdum
868!read(*,*)
869!endif
870                enddo
871              endif
872! DEBUG
873!write(6,'(5i3,2x,2f16.10)') ix,iy,iz,itet,iww,cdum
874            enddo
875! DEBUG
876!write(6,*) 'Khst7',xdum
877!read(*,*)
878!write(6,*) "second set"
879            do iww=iwwlim(2),iwwlim(3) ! energies between evtx(2) & evtx(3)
880! DEBUG
881!write(6,*) iww
882              wwwlo=max(dw*(iww-0.5d0),evtx(2))
883              wwwhi=min(dw*(iww+0.5d0),evtx(3))
884              if (abs(wwwhi-wwwlo)*1.d14.lt.max(abs(wwwhi),abs(wwwlo),dw)) cycle
885              do ii=1,3
886              do jj=1,3
887! DEBUG
888k_ct_DEBUG = k_ct_DEBUG + 1
889!write(6,*) "k_ct_DEBUG: ",k_ct_DEBUG,ii,jj
890!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "E_LIM:      ",wwwlo,wwwhi
891                call tsinggrater(ii,jj,wwwlo,wwwhi,2,xk,kvtx(:,1),ngkpt,evtx,abr,rlr,nsing,wsing,sint1a)
892! DEBUG
893k_ct_DEBUG = k_ct_DEBUG + 1
894!write(6,*) "k_ct_DEBUG: ",k_ct_DEBUG,ii,jj
895!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "E_LIM:      ",wwwlo,wwwhi
896                call tsinggrater(ii,jj,wwwlo,wwwhi,3,xk,kvtx(:,1),ngkpt,evtx,abr,rlr,nsing,wsing,sint1b)
897                stens(ii,jj) = sint1a+sint1b
898              enddo
899              enddo
900              stens = stens * xmult
901! DEBUG
902!cdum = 0.d0
903              if (lc) then
904                do iw=1,nwpt
905                  ie=iww-ocsign*iw
906                  do ii=1,3
907                  do jj=1,3
908                      ssi(ie)=ssi(ie)+vsm2(ii,jj,iw)*stens(jj,ii)
909! DEBUG
910!if (ie.eq.0.or.ie.eq.1) then
911!cdum = cdum+vsm2(ii,jj,iw)*stens(jj,ii)
912!endif
913                  enddo
914                  enddo
915                enddo
916              endif
917! DEBUG
918!write(6,'(5i3,2x,2f16.10)') ix,iy,iz,itet,iww,cdum
919            enddo
920! DEBUG
921!write(6,*) 'Khst7',xdum
922!write(6,*) "third set"
923            do iww=iwwlim(3),iwwlim(4) ! energies between evtx(3) & evtx(4)
924! DEBUG
925!write(6,*) iww
926              wwwlo=max(dw*(iww-0.5d0),evtx(3))
927              wwwhi=min(dw*(iww+0.5d0),evtx(4))
928              if (abs(wwwhi-wwwlo)*1.d14.lt.max(abs(wwwhi),abs(wwwlo),dw)) cycle
929              do ii=1,3
930              do jj=1,3
931! DEBUG
932k_ct_DEBUG = k_ct_DEBUG + 1
933!write(6,*) "k_ct_DEBUG: ",k_ct_DEBUG,ii,jj
934!if (k_ct_DEBUG.eq.k_ref_DEBUG) write(6,*) "E_LIM:      ",wwwlo,wwwhi
935                call tsinggrater(ii,jj,wwwlo,wwwhi,4,xk,kvtx(:,1),ngkpt,evtx,abr,rlr,nsing,wsing,stens(ii,jj))
936              enddo
937              enddo
938              stens = stens * xmult
939! DEBUG
940!cdum = 0.d0
941              if (lc) then
942                do iw=1,nwpt
943                  ie=iww-ocsign*iw
944                  do ii=1,3
945                  do jj=1,3
946                      ssi(ie)=ssi(ie)+vsm2(ii,jj,iw)*stens(jj,ii)
947! DEBUG
948!if (ie.eq.0.or.ie.eq.1) then
949!cdum = cdum+vsm2(ii,jj,iw)*stens(jj,ii)
950!endif
951                  enddo
952                  enddo
953                enddo
954              endif
955! DEBUG
956!write(6,'(5i3,2x,2f16.10)') ix,iy,iz,itet,iww,cdum
957            enddo
958! DEBUG
959!write(6,*) 'Khst7',xdum
960!write(32,'(4i3,2x,f14.10,2x,f14.10,"  sing")') ixx,itet,(cdum/2.d0)
961!write(6,*) "finished singular integration"
962          elseif (centerint) then
963! anisotropic part
964! DEBUG
965!cdum = 0.d0
966            if (lc) then
967              do iww=iwwlim(1)-1,iwwlim(4)+1
968                wwwlo=dw*(iww-0.5d0)
969!                wwwhi=dw*(iww+0.5d0)
970                do iw=1,nwpt
971                  call vnitetint(rrpyr,evtx0,wwwlo,dw,vjfactor(:,:,iw),bmet,vcenterint,cint)
972                  ie=iww-ocsign*iw
973                  ssi(ie)=ssi(ie)+cint/volelmnt
974! DEBUG
975!if (ie.eq.0.or.ie.eq.1) cdum = cdum+cint/volelmnt
976                enddo
977              enddo
978            endif
979! DEBUG
980!write(32,'(4i3,2x,f14.10,2x,f14.10,"  cent")') ixx,itet,(cdum/2.d0)
981          endif
982! can use basic tetrahedron integration for everything else
983! DEBUG
984!cdum = 0.d0
985          if (lc) then
986            do iw=1,nwpt
987              aa0(iw)=vpyr0(indxe(1),iw)   ! base value of function
988              av(1:3,iw)=(/0.d0,0.d0,0.d0/) ! gradient of function
989              do jj=1,3
990                av(1:3,iw)=av(1:3,iw) &
991&                      +(vpyr0(indxe(jj+1),iw)-vpyr0(indxe(1),iw))*rg(1:3,jj)
992              enddo
993            enddo
994            do iww=iwwlim(1)-1,iwwlim(4)+1
995              wwwlo=dw*(iww-0.5d0)
996!              wwwhi=dw*(iww+0.5d0)
997              call mkvsint(rrpyr(1:3,indxe),xk,volred,evtx,wwwlo,dw,sint1,svec)
998              do iw=1,nwpt
999                ie=iww-ocsign*iw
1000                ssi(ie)=ssi(ie)+(aa0(iw)*sint1+dot_product(svec,av(:,iw)))/volelmnt
1001! DEBUG
1002!if (ie.eq.0.or.ie.eq.1) cdum = cdum+(aa0(iw)*sint1+dot_product(svec,av(:,iw)))/volelmnt
1003              enddo
1004            enddo
1005          endif
1006! DEBUG
1007!write(32,'(4i3,2x,f14.10,2x,f14.10)') ixx,itet,(cdum/2.d0)
1008        enddo
1009      enddo
1010      enddo
1011      enddo
1012!    if (lc) ssi=-ssi/(dble(ngkpt(1)*ngkpt(2)*ngkpt(3))*vol*pi)
1013      if (lc) ssi=-ocsign*ssi/((2*pi)**3*pi)
1014
1015! DEBUG
1016!write(6,*) "Integration complete"
1017
1018      if (lc) then
1019        do ie=max(-nwpt,ielo),min(nwpt,iehi)
1020          ploss1(ie)=ploss1(ie)+abs(dimag(ssi(ie)))*pi
1021        enddo
1022      endif
1023      deallocate (ssi)
1024    endif
1025
1026  enddo
1027enddo
1028
1029do ie=-nwpt,nwpt
1030  eps1=dble(ie)*dw-dw/2
1031!  write(14,'(i4,f10.3,2(3x,2es12.3))') ie,eps1*27.2114,ploss1(ie)
1032enddo
1033write(14,*)
1034if (itetrahedron.eq.1) then
1035  iw = int(total_se/dw)
1036  rw = total_se - iw*dw
1037!  write(14,*) iw,rw
1038  ploss=ploss1(iw)*((dw-rw)/dw)+ploss1(iw+1)*(rw/dw)
1039!  ploss=(ploss1(0)+ploss1(1))/2.
1040!  write(14,*) ploss
1041endif
1042
1043!write(6,*) "kwls ",total_se,ploss
1044
1045!temparray=log10(temparray)
1046!temparray=10**(temparray-int(temparray))
1047!do ix=1,ngkpt(1)
1048!  write(76,*) ix,'--------------------------------------------------'
1049!  do iy=1,ngkpt(2)
1050!!    write(76,'(10es8.1)') dble(temparray(ix,iy,:))
1051!    write(76,'(10f8.5)') dble(temparray(ix,iy,:))
1052!  enddo
1053!enddo
1054
1055return
1056end subroutine mkploss
1057
1058