1!==========================================================================
2!
3! lin_denominator_m   Originally By Liang Z. Tan         2/8/2013
4!
5!    computes the integral of
6!
7!      Theta(Ef-Ev(k))Theta(Ec(k+q)-Ef)/(Ev(k)-Ec(k+q)+w+i0+)
8!
9!    where the range of integration over k is the mini-BZ and
10!    Ec(k+q) and Ev(k) are approximated as linear functions of k
11!    (i.e. planes).
12!    Currently implemented for 2D.
13!    z-direction is the out-of-plane direction.
14!
15!    input: k, energies and velocities at k, k+q, frequency, Ef.
16!    output: (complex) value of the integral.
17!
18!==========================================================================
19
20
21
22#include "f_defs.h"
23
24module lin_denominator_m
25
26  use global_m
27  use sort_m
28  implicit none
29
30  !> define boundaries as a useful type
31  !! x = slope * y + intercept
32  type boundary
33    real(DP) :: slope, intercept
34    real(DP) :: sgn
35  end type boundary
36
37  private
38  public :: integrate_mbz
39
40contains
41
42  !> this integrates log|ax+b| over x in [0,1].
43  real(DP) function integratelog(a,b)
44    real(DP), intent(in) :: a,b
45
46    PUSH_SUB(integratelog)
47
48    if(abs(a) > TOL_Zero .and. abs(b)>TOL_Zero) then
49      if(abs(a+b) > TOL_Zero) then
50        integratelog = -1.0d0 -(b/a)*log(abs(b))+((a+b)/a)*log(abs(a+b))
51      else
52        integratelog = -1.0d0 -(b/a)*log(abs(b))
53      endif
54    elseif(abs(a) > TOL_Zero .and. abs(b)<TOL_Zero) then
55      if(abs(a+b) > TOL_Zero) then
56        integratelog = -1.0d0 +((a+b)/a)*log(abs(a+b))
57      else
58        integratelog = -1.0d0
59      endif
60    else
61      integratelog = log(abs(b))
62    endif
63    !write(6,*) "integratelog: a= ",a," b= ",b," ans= ",integratelog
64
65    POP_SUB(integratelog)
66
67  end function integratelog
68
69  !> evaluate a linear function at a point
70  real(DP) function evalb(b,y)
71    type(boundary), intent(in) :: b
72    real(DP), intent(in) :: y
73
74    PUSH_SUB(evalb)
75
76    evalb = b%slope * y + b%intercept
77
78    POP_SUB(evalb)
79  end function evalb
80
81  !> intersection point of two lines
82  real(DP) function critfn(b1,b2)
83    type(boundary), intent(in) :: b1, b2
84
85    PUSH_SUB(critfn)
86    if(b1%slope /= b2%slope) then
87      critfn = (b2%intercept - b1%intercept)/(b1%slope - b2%slope)
88    else
89      critfn = 100
90    end if
91
92    POP_SUB(critfn)
93  end function critfn
94
95
96  !> imaginary part of the integral 1/(ux*kx+uy*ky+1+J*eta)
97  !! over a trapezoid in 2D
98  real(DP) function trapzimagpart(ux,uy,y0,y1,bdl,bdu)
99    real(DP),intent(in) :: ux,uy,y0,y1
100    type(boundary),intent(in) :: bdl, bdu
101
102    real(DP) :: yl,yu,x1,x0,x0l,x0u,x1l,x1u,pi,yc
103    type(boundary) :: bdsing
104
105    PUSH_SUB(trapzimagpart)
106
107    pi = 4.d0*datan(1.d0)
108
109    if(evalb(bdl,(y0+y1)/2) > evalb(bdu,(y0+y1)/2)) then
110      trapzimagpart = 0d0
111    elseif(abs(ux)<TOL_Zero .and. abs(uy)<TOL_Zero) then
112      trapzimagpart = 0d0
113    elseif (abs(ux)<TOL_Zero) then
114      yc = -1.0d0/uy
115      if( yc<y1 .and. yc>y0) then
116        trapzimagpart = -pi*(evalb(bdu,yc)-evalb(bdl,yc))/abs(uy)
117      else
118        trapzimagpart = 0d0
119      endif
120    else
121      bdsing%slope = -uy/ux
122      bdsing%intercept = -1.0d0/ux
123      yu = critfn(bdsing,bdu)
124      yl = critfn(bdsing,bdl)
125      x1 = evalb(bdsing,y1)
126      x0 = evalb(bdsing,y0)
127      x0l= evalb(bdl,y0)
128      x0u= evalb(bdu,y0)
129      x1l= evalb(bdl,y1)
130      x1u= evalb(bdu,y1)
131      if ((x0<x0l .and. x1<x1l) .or. (x0>x0u .and. x1>x1u)) then
132        trapzimagpart=0d0
133      elseif(x0>x0l .and. x0<x0u .and. x1>x1l .and. x1<x1u) then
134        trapzimagpart= -pi*abs(y1-y0)/abs(ux)
135      elseif((x0<x0l .and. x1>x1u) .or. (x0>x0u .and. x1<x1l))then
136        trapzimagpart = -pi*abs(yl-yu)/abs(ux)
137      elseif(x0<x0l .and. x1>x1l .and. x1<x1u) then
138        trapzimagpart = -pi*abs(y1-yl)/abs(ux)
139      elseif(x0>x0u .and. x1>x1l .and. x1<x1u) then
140        trapzimagpart = -pi*abs(y1-yu)/abs(ux)
141      elseif(x1<x1l .and. x0>x0l .and. x0<x0u) then
142        trapzimagpart = -pi*abs(yl-y0)/abs(ux)
143      !elseif(x1>x1u .and. x0>x0l .and. x0<x0u) then
144      else
145        trapzimagpart = -pi*abs(yu-y0)/abs(ux)
146      endif
147    endif
148
149    POP_SUB(trapzimagpart)
150  end function trapzimagpart
151
152  !> real part of the integral 1/(ux*kx+uy*ky+1+J*eta)
153  !! over a trapezoid in 2D
154  real(DP) function trapzrealpart(ux,uy,y0,y1,bdl,bdu)
155    real(DP),intent(in) :: ux,uy,y0,y1
156    type(boundary),intent(in) :: bdl, bdu
157
158    real(DP) :: a,b,part1, part2
159
160    PUSH_SUB(trapzrealpart)
161
162    if(evalb(bdl,(y0+y1)/2) > evalb(bdu,(y0+y1)/2)) then
163      trapzrealpart = 0d0
164    elseif( abs(ux)<TOL_Small .and. abs(uy)<TOL_Small) then
165      trapzrealpart = 0.5*(y1-y0)&
166      *(evalb(bdu,y0)-evalb(bdl,y0)+evalb(bdu,y1)-evalb(bdl,y1))
167    elseif(abs(ux)<TOL_Small) then
168      a = (bdu%slope-bdl%slope)/uy
169      b = bdu%intercept-bdl%intercept - a
170      part1 = (y1-y0)*a
171      part2 = (b/uy)*(log(abs(uy*y1+1.0d0))-log(abs(uy*y0+1.0d0)))
172      trapzrealpart = part1 + part2
173    else
174      a = (ux*bdu%slope+uy)*(y1-y0)
175      b = (ux*bdu%slope+uy)*y0+ux*bdu%intercept+1.0d0
176      part1 = (y1-y0)/ux * integratelog(a,b)
177      a = (ux*bdl%slope+uy)*(y1-y0)
178      b = (ux*bdl%slope+uy)*y0+ux*bdl%intercept+1.0d0
179      part2 = (y1-y0)/ux * integratelog(a,b)
180      trapzrealpart = part1 - part2
181    endif
182
183    POP_SUB(trapzrealpart)
184  end function trapzrealpart
185
186  !> picks the relevant boundaries from a list
187  subroutine actb(bs,y,b0,b1)
188    type(boundary), intent(in) :: bs(4)
189    type(boundary), intent(out) :: b0,b1
190
191    type(boundary) :: b,bs0(4), bs1(4)
192    integer :: i,nbs0,nbs1, xidx
193    real(DP), intent(in) :: y
194    real(DP) :: x,xtry
195
196    PUSH_SUB(actb)
197
198    nbs0=0
199    nbs1=0
200    do i = 1,4
201      b = bs(i)
202      if(b%sgn>0) then
203        nbs0 = nbs0 +1
204        bs0(nbs0) = b
205      else
206        nbs1 = nbs1 +1
207        bs1(nbs1) = b
208      endif
209    enddo
210
211    x = -100
212    do i = 1,nbs0
213      xtry = evalb(bs0(i),y)
214      if(xtry> x) then
215        x = xtry
216        xidx = i
217      endif
218    enddo
219    b0 = bs0(xidx)
220
221    x = 100
222    do i = 1,nbs1
223      xtry = evalb(bs1(i),y)
224      if(xtry<x) then
225        x = xtry
226        xidx = i
227      endif
228    enddo
229    b1 = bs1(xidx)
230    POP_SUB(actb)
231  end subroutine actb
232
233  !> computes the interval (ymin,ymax) from a set of n linear inequalities
234  !! that y satisfies.
235  !! Each inequality is modelled as a boundary a*y + c > 0 or a*y + c < 0
236  !! where a = slope, c = intercept, > : sgn=1.0, < : sgn=-1.0
237  !! Implicit bounds (0,1) assumed.
238  !! If there is no solution, a zero-length interval is returned.
239  subroutine intervals(n,ineqs,ymin,ymax)
240    integer,intent(in) :: n
241    type(boundary),intent(in) :: ineqs(:)
242    real(DP),intent(out) :: ymin,ymax
243
244    integer :: ii
245    real(DP) :: a,c,s,y2
246
247    PUSH_SUB(intervals)
248
249    ymin = 0.0d0
250    ymax = 1.0d0
251
252    do ii = 1,n
253      a=ineqs(ii)%slope
254      c=ineqs(ii)%intercept
255      s=ineqs(ii)%sgn
256      if(abs(a)<TOL_Small)then
257        if(s*c<0.0d0) then
258          ymax = ymin
259          exit
260        endif
261      else
262        if(s*a>0.0) then
263          y2 = -c/a
264          if(y2>ymin)then
265            ymin = y2
266          endif
267          if(ymin>ymax) then
268            ymin = ymax
269            exit
270          endif
271        else
272          y2 = -c/a
273          if(y2<ymax)then
274            ymax = y2
275          endif
276          if(ymax<ymin) then
277            ymax = ymin
278            exit
279          endif
280        endif
281      endif
282    enddo
283
284    POP_SUB(intervals)
285  end subroutine intervals
286
287  !> exact value of the integral Theta(a.r+c)Theta(b.r+c)/(u.r+1)
288  !! over the unit square in 2D.
289  complex(DP) function exactmodel1(ax,ay,c,bx,by,d,ux,uy)
290    real(DP), intent(in) :: ax,ay,c,bx,by,d,ux,uy
291
292    type(boundary) :: bd0, bd1, bd2, bd3, bds(4)
293    type(boundary) :: bl,bu,bls(7),bus(7),intv1,intv2
294    real(DP) :: critpts0(5), y(7), yval, accu_real, accu_imag,ymin,ymax
295    complex(DP) :: accu
296    integer ::  critpts0_idx(5)
297    integer :: num_crit
298    integer :: i
299
300    PUSH_SUB(exactmodel1)
301
302    !make boundaries
303    bd0%slope = 0d0
304    bd0%intercept = 0d0
305    bd0%sgn = 1.0d0
306    bd1%slope = 0d0
307    bd1%intercept = 1.0d0
308    bd1%sgn = -1.0d0
309
310    if(abs(ax)>TOL_Small .and. abs(bx)>TOL_Small) then
311      bd2%slope = -ay/ax
312      bd2%intercept = -c/ax
313      bd2%sgn = sign(1.0d0,ax)
314      bd3%slope = -by/bx
315      bd3%intercept = -d/bx
316      bd3%sgn = sign(1.0d0,bx)
317
318      ymin = 0.0d0
319      ymax = 1.0d0
320    elseif(abs(ax)<TOL_Small .and. abs(bx)>TOL_Small) then
321      bd2%slope = 1.0d0 !< dummy boundary
322      bd2%intercept = 100.0 !<dummy boundary
323      bd2%sgn = -1.0d0 !<dummy boundary
324
325      bd3%slope = -by/bx
326      bd3%intercept = -d/bx
327      bd3%sgn = sign(1.0d0,bx)
328
329      intv1 = boundary(ay,c,1.0d0)
330      call intervals(1,(/intv1/),ymin,ymax)
331      !write(6,*) "ymin= ",ymin,"ymax= ",ymax
332
333
334    elseif(abs(ax)>TOL_Small .and. abs(bx)<TOL_Small) then
335      bd2%slope = -ay/ax
336      bd2%intercept = -c/ax
337      bd2%sgn = sign(1.0d0,ax)
338      bd3%slope = 1.0d0 !< dummy boundary
339      bd3%intercept = 100.0d0 !<dummy boundary
340      bd3%sgn = -1.0d0 !<dummy boundary
341
342      intv1 = boundary(by,d,1.0d0)
343      call intervals(1,(/intv1/),ymin,ymax)
344
345    else !both ax and bx small
346      bd2%slope = 0.9d0 !< dummy boundary
347      bd2%intercept = 101.0d0 !<dummy boundary
348      bd2%sgn = -1.0d0 !<dummy boundary
349      bd3%slope = 1.0d0 !< dummy boundary
350      bd3%intercept = 100.0d0 !<dummy boundary
351      bd3%sgn = -1.0d0 !<dummy boundary
352
353      intv1 = boundary(ay,c,1.0d0)
354      intv2 = boundary(by,d,1.0d0)
355      call intervals(2,(/intv1,intv2/),ymin,ymax)
356      !write(6,*) "ymin= ",ymin,"ymax= ",ymax
357    endif
358
359    bds(1) = bd0
360    bds(2) = bd1
361    bds(3) = bd2
362    bds(4) = bd3
363
364    !first split integration region into trapezoids
365    !start by finding where the boundaries intersect
366    critpts0(1) = critfn(bd0,bd2)
367    critpts0(2) = critfn(bd0,bd3)
368    critpts0(3) = critfn(bd1,bd2)
369    critpts0(4) = critfn(bd1,bd3)
370    critpts0(5) = critfn(bd2,bd3)
371
372    call sortrx(5,critpts0,critpts0_idx)
373    num_crit = 2
374    y(1) = ymin
375    do i = 1,5
376      yval = critpts0(critpts0_idx(i))
377      if(yval>ymin .and. yval<ymax) then
378        y(num_crit) = yval
379        num_crit = num_crit + 1
380      endif
381    enddo
382    y(num_crit) = ymax
383
384
385    !make the list of lower and upper boundaries
386    do i = 1,num_crit-1
387      call actb(bds,(y(i)+y(i+1))/2,bl,bu)
388      bls(i) = bl
389      bus(i) = bu
390    enddo
391
392    !sum over all the trapezoids
393    accu = (0,0)
394    do i = 1,num_crit-1
395      !write(6,*) "ux= ",ux,"uy= ",uy,"yi= ",y(i),"yi1= ",y(i+1),&
396      !"bls= ",bls(i),"bus= ",bus(i)
397      accu_real = trapzrealpart(ux,uy,y(i),y(i+1),bls(i),bus(i))
398      accu_imag = trapzimagpart(ux,uy,y(i),y(i+1),bls(i),bus(i))
399      accu = accu + CMPLX(accu_real,accu_imag)
400    enddo
401
402    exactmodel1 = accu
403
404    !write(6,*) "exactmodel1= ",exactmodel1
405
406    POP_SUB(exactmodel1)
407  end function exactmodel1
408
409  !> compute integral in miniBZ from energies and velocities
410  !! NOTE: currently for 2D only. z-direction is the out-of-plane direction.
411  complex(DP) function integrate_mbz(kpts,k,ek,ekq,vk,vkq,w,ef)
412    type(kpoints), intent(in) :: kpts
413    real(DP), intent(in) :: k(3)
414    real(DP), intent(in):: ek, ekq, vk(2), vkq(2), w,ef
415
416    real(DP) :: kx,ky,dkx,dky,kxmin,kymin,vkx,vky,vkqx,vkqy, &
417                cc,ux,uy,ax,ay,c,bx,by,d
418    complex(DP) :: accu
419    integer :: ii, ix, iy
420
421    PUSH_SUB(integrate_mbz)
422
423    ! need to know which is the out of plane direction
424    ix = 1
425    iy = 2
426    do ii = 1,3
427      if(kpts%kgrid(ii) == 1) then
428        if(ii <= ix) then
429          ix = ix +1
430        endif
431        if(ii <= iy) then
432          iy = iy +1
433        endif
434      endif
435    enddo
436
437    dkx = 1.0d0/real(kpts%kgrid(ix))
438    dky = 1.0d0/real(kpts%kgrid(iy))
439    kx = k(ix)
440    ky = k(iy)
441    vkx = vk(1)
442    vky = vk(2)
443    vkqx = vkq(1)
444    vkqy = vkq(2)
445    kxmin = kx-dkx/2.0
446    kymin = ky-dky/2.0
447
448    cc = (vkx-vkqx)*kxmin + (vky-vkqy)*kymin + &
449         (ek - ekq +(vkqx-vkx)*kx +(vkqy-vky)*ky + w)
450    ux = (vkx-vkqx)*dkx/cc
451    uy = (vky-vkqy)*dky/cc
452
453
454    accu = (0,0)
455    !k occ, kq unocc
456    ax = vkx*dkx
457    ay = vky*dky
458    c = vkx*kxmin+vky*kymin-vkx*kx-vky*ky+ek-ef
459    bx = -vkqx*dkx
460    by = -vkqy*dky
461    d = -vkqx*kxmin-vkqy*kymin+vkqx*kx+vkqy*ky+(ef-ekq)
462    accu = accu -  exactmodel1(ax,ay,c,bx,by,d,ux,uy)/cc
463
464    !if(dble(accu)> 0.0) then
465    !  write(6,*) "accu1= ",accu,"vkx= ",vkx,"vky= ",vky,&
466    !  "vkqx= ",vkqx,"vkqy= ",vkqy,&
467    !  "ax",ax,"ay",ay,"bx",bx,"by",by,&
468    !  "c= ",c,"d= ",d,"cc= ",cc
469    !endif
470
471    ! k occ, kq unocc
472    ax = -vkx*dkx
473    ay = -vky*dky
474    c = -1.0d0*(vkx*kxmin+vky*kymin-vkx*kx-vky*ky+ek-ef)
475    bx = vkqx*dkx
476    by = vkqy*dky
477    d = -1.0d0*(-vkqx*kxmin-vkqy*kymin+vkqx*kx+vkqy*ky+(ef-ekq))
478    accu = accu +  exactmodel1(ax,ay,c,bx,by,d,ux,uy)/cc
479
480    !if(dble(accu)> 0.0) then
481    !  write(6,*) "accu2= ",accu,"vkx= ",vkx,"vky= ",vky,&
482    !  "vkqx= ",vkqx,"vkqy= ",vkqy,&
483    !  "ax",ax,"ay",ay,"bx",bx,"by",by,&
484    !  "c= ",c,"d= ",d,"cc= ",cc
485    !endif
486
487    integrate_mbz = accu
488    POP_SUB(integrate_mbz)
489  end function integrate_mbz
490
491end module lin_denominator_m
492