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