1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8! This code segment has been fully created by: 9! Nick Papior Andersen, 2013, nickpapior@gmail.com 10! 11 12module m_ts_contour_eq 13 14 use precision, only : dp 15 16! Use the type associated with the contour 17! Maybe they should be collected to this module. 18! However, I like this partition. 19 use m_ts_electype 20 21 use m_ts_cctype 22 use m_ts_chem_pot 23 use m_ts_io_contour 24 use m_ts_io_ctype 25 use m_ts_aux 26 27 implicit none 28 29 ! equilibrium contour IO-segments 30 integer, save, public :: N_Eq 31 type(ts_c_io), pointer, save, public :: Eq_io(:) => null() 32 type(ts_cw) , pointer, save, public :: Eq_c(:) => null() 33 34 ! The contours for the equilibrium density are attributed a fruitful discussion with 35 ! Hans Skriver. Previously the routine names reflected his contribution. 36 ! However, for programming clarity we have employed a different naming scheme. 37 ! In order to retain the contributions it is encouraged to keep this sentiment for his 38 ! contributions. 39 40 ! For further attributions see the original paper by Brandbyge, et. al, 2002: DOI: 10.1103/PhysRevB.65.165401 41 42 ! this is heavily linked with CONTOUR_NEQ from m_ts_contour_neq 43 integer, parameter, public :: CONTOUR_EQ = 1 44 45 public :: read_contour_eq_options 46 public :: print_contour_eq_options 47 public :: print_contour_eq_block 48 public :: io_contour_eq 49 public :: N_Eq_E 50 public :: get_c_io_index 51 public :: Eq_linear_index 52 public :: Eq_E 53 public :: c2energy 54 public :: c2weight_eq 55 public :: ID2idx 56 57 interface ID2idx 58 module procedure ID2idx_cw 59 module procedure ID2idx_cidx 60 end interface 61 62 private 63 64contains 65 66 subroutine read_contour_eq_options(N_mu, mus, Volt) 67 68 use units, only : Pi, Kelvin, eV 69 use fdf 70 71 integer, intent(in) :: N_mu 72 type(ts_mu), intent(inout) :: mus(N_mu) 73 real(dp), intent(in) :: Volt 74 75 integer :: i, j, k, c_mu 76 integer :: N 77 character(len=C_N_NAME_LEN), allocatable :: tmp(:), nContours(:) 78 character(len=C_N_NAME_LEN) :: tmp_one 79 integer :: cur 80 integer, allocatable :: idx(:) 81 logical :: isStart, isTail 82 real(dp) :: r1, r2 83 84 call fdf_obsolete('TS.ComplexContour.NPoles') 85 86 87 ! We only allow the user to either use the old input format, or the new 88 ! per-electrode input 89 do i = 1 , N_mu 90 if ( Eq_segs(mus(i)) == 1 ) then 91 mus(i)%Eq_seg(1) = 'cont-frac-'//trim(mus(i)%name) 92 else 93 write(mus(i)%Eq_seg(Eq_segs(mus(i))),'(a,i0)') 'pole',i 94 end if 95 end do 96 97 ! Count the number of contour segments 98 N = minval(Eq_segs(mus)) 99 if ( N == 0 ) then 100 call die('You have not assigned any contours for the chemical & 101 &potentials.') 102 end if 103 N = sum(Eq_segs(mus)) 104 105 ! collect all equilibrium names 106 allocate(tmp(N)) 107 tmp(:) = ' ' 108 N = 0 109 do i = 1 , N_mu 110 do j = 1 , Eq_segs(mus(i)) - 1 111 N = N + 1 112 tmp(N) = mus(i)%Eq_seg(j) 113 end do 114 end do 115 ! we need to have the pole contours as the last ones 116 ! to be able to easily differ them from the read in values 117 do i = 1 , N_mu 118 j = Eq_segs(mus(i)) 119 N = N + 1 120 tmp(N) = mus(i)%Eq_seg(j) 121 end do 122 123 ! find all unique equilibrium names 124 N = 1 125 uniq_names: do i = 2 , size(tmp) 126 j = 0 127 do 128 j = j + 1 129 if ( i <= j ) exit 130 if ( tmp(i) .eq. tmp(j) ) cycle uniq_names 131 end do 132 N = N + 1 133 end do uniq_names 134 135 allocate(nContours(N)) ! unique names 136 nContours = ' ' 137 138 ! populate unique equilibrium names 139 N = 0 140 uniq_names_pop: do i = 1 , size(tmp) 141 do j = 1 , N 142 if ( tmp(i) .eq. nContours(j) ) cycle uniq_names_pop 143 end do 144 N = N + 1 145 if ( N > size(nContours) ) call die('Unique contour setup went wrong, contact devs') 146 nContours(N) = tmp(i) 147 end do uniq_names_pop 148 if ( N /= size(nContours) ) then 149 call die('ERROR: We have not populated enough contours') 150 end if 151 deallocate(tmp) 152 153 ! Allocate space for the equilibrium contours 154 nullify(Eq_io,Eq_c) 155 N_Eq = N 156 allocate(Eq_io(N_Eq)) ! the equilibrium contour io container 157 allocate(Eq_c(N_Eq)) ! the equilibrium contour container 158 Eq_io(:)%type = 'eq' ! set the equilibrium tag 159 160 ! Attach all the contours 161 do i = 1 , N 162 Eq_c(i)%c_io => Eq_io(i) 163 end do 164 165 ! read in the equilibrium contours (the last will be the poles) 166 do i = 1 , N - N_mu 167 168 ! Save the name of the contour (to be able to find the 169 ! associated chemical potential) 170 Eq_io(i)%name = nContours(i) 171 172 ! Get index of associated chemical potential 173 c_mu = 0 174 do j = 1 , N_mu 175 if ( hasC(mus(j),Eq_io(i)) ) then 176 c_mu = j 177 exit 178 end if 179 end do 180 if ( c_mu == 0 ) then 181 print *,Eq_io(i)%name 182 call die('Error in determining associated & 183 &chemical potential') 184 end if 185 ! Remove the name, it is needed to be ' ' 186 Eq_io(i)%name = ' ' 187 188 tmp_one = nContours(i) 189 190 ! read in the contour 191 if ( tmp_one(1:1) == '*' .and. & 192 .not. ts_exists_contour_block('TS','',tmp_one(2:)//' ') ) then 193 194 if ( N_mu > 2 ) then 195 ! here we can have a zero bias case where 196 ! N_mu == 1 :) 197 call die('When using other than 2 electrodes and & 198 &chemical potentials you are forced to setup & 199 &the contour input yourself.') 200 end if 201 202 ! this is a fake-contour, read it in 203 ! if it exists, else create it... 204 ! *** NOTE this is hard coded against the chem_pot code 205 206 Eq_io(i)%name = tmp_one 207 if ( tmp_one(2:2) == 'c' ) then 208 ! the circle contour 209 Eq_io(i)%part = 'circle' 210 Eq_io(i)%N = 25 211 write(Eq_io(i)%cN,'(i0)') Eq_io(i)%N 212 if ( tmp_one(4:4) == 'l' ) then ! left 213 Eq_io(i)%ca = '-40. eV + V/2' 214 Eq_io(i)%a = - 40._dp * eV + Volt * .5_dp 215 Eq_io(i)%cb = '-10 kT + V/2' 216 Eq_io(i)%b = -10._dp * mus(c_mu)%kT + Volt * .5_dp 217 else ! must be right 218 Eq_io(i)%ca = '-40. eV - V/2' 219 Eq_io(i)%a = - 40._dp * eV - Volt * .5_dp 220 Eq_io(i)%cb = '-10 kT - V/2' 221 Eq_io(i)%b = -10._dp * mus(c_mu)%kT - Volt * .5_dp 222 end if 223 Eq_io(i)%method = 'g-legendre' 224 call c_io_add_opt(Eq_io(i),'right','right') 225 226 else 227 Eq_io(i)%part = 'tail' 228 Eq_io(i)%N = 10 229 write(Eq_io(i)%cN,'(i0)') Eq_io(i)%N 230 Eq_io(i)%ca = 'prev' 231 if ( tmp_one(4:4) == 'l' ) then ! left 232 Eq_io(i)%a = -10._dp * mus(c_mu)%kT + Volt * .5_dp 233 else ! must be right 234 Eq_io(i)%a = -10._dp * mus(c_mu)%kT - Volt * .5_dp 235 end if 236 Eq_io(i)%cb = 'inf' 237 Eq_io(i)%b = huge(1._dp) 238 Eq_io(i)%method = 'g-fermi' 239 end if 240 else 241 call ts_read_contour_block('TS','',nContours(i),Eq_io(i), & 242 mus(c_mu)%kT, Volt) 243 end if 244 245 end do 246 247 ! We here create the "fake" pole contours 248 j = 1 249 do i = N - N_mu + 1, N 250 if ( Eq_segs(mus(j)) == 1 ) then 251 Eq_io(i)%part = 'cont-frac' 252 Eq_io(i)%method = 'continued-fraction' 253 ! Currently this is just a very high number 254 Eq_io(i)%a = 1.e10_dp * eV ! the continued fraction infinity point 255 else 256 Eq_io(i)%part = 'pole' 257 Eq_io(i)%method = 'residual' 258 Eq_io(i)%a = mus(j)%mu 259 end if 260 ! assign name to the Eq_io 261 Eq_io(i)%name = nContours(i) 262 Eq_io(i)%N = mus(j)%N_poles 263 Eq_io(i)%b = mus(j)%kT ! Save kT in the contour 264 Eq_io(i)%d = mus(j)%mu ! save the chemical potential here 265 j = j + 1 266 267 end do 268 269 ! fix the read-in constants 270 do i = 1 , N_mu 271 272 cur = get_c_io_index(mus(i)%Eq_seg(1)) 273 if ( cur == 0 ) then 274 call die('A terrible error has occured, please inform the & 275 &developers') 276 end if 277 ! Two cases can happen, 278 ! 1) A continued fraction method 279 ! 2) The regular circle contour 280 if ( leqi(Eq_io(cur)%part,'cont-frac') ) then 281 282 ! The segment _has_ to be alone 283 if ( Eq_segs(mus(i)) > 1 ) then 284 call die('Continued fraction contours & 285 &are individual and cannot be connected.') 286 end if 287 288 else 289 290 ! first fix the contours (not the poles) 291 k = Eq_segs(mus(i)) - 1 292 293 allocate(idx(k)) 294 295 ! Create index-table 296 do j = 1 , k 297 idx(j) = get_c_io_index(mus(i)%Eq_seg(j)) 298 end do 299 300 ! Fix contours 301 call ts_fix_contours(N_Eq, Eq_io, k, idx) 302 303 deallocate(idx) 304 305 ! Check equilibrium settings 306 isStart = leqi(Eq_io(cur)%part,'circle') .or. & 307 leqi(Eq_io(cur)%part,'square') 308 isTail = leqi(Eq_io(cur)%part,'tail') 309 310 ! we should not check the pole 311 do j = 1 , k 312 cur = get_c_io_index(mus(i)%Eq_seg(j)) 313 call consecutive_types(Eq_io(cur),isStart,isTail, & 314 mus(i)%mu, mus(i)%kT) 315 end do 316 317 ! check that the last contour actually lies on the RHS 318 ! of the chemical potential, at least 10 kT across! 319 if ( Eq_io(cur)%b < mus(i)%mu + 10._dp * mus(i)%kT ) then 320 print *,'Contour name: ',trim(Eq_io(cur)%name) 321 write(*,*) 'Energies are too close: ',Eq_io(cur)%b,mus(i)%mu + 10._dp * mus(i)%kT 322 call eq_die('The last contour of the chemical potential: & 323 &'//trim(Name(mus(i)))//' lies too close to the & 324 &chemical potential. It must be at least 10 kT from mu.') 325 end if 326 327 end if 328 329 end do 330 331 ! Allocate sizes of the contours 332 do i = 1 , N_Eq 333 ! number of different weights for each energy-point 334 k = count(hasC(mus,Eq_c(i))) 335 if ( k == 0 ) then 336 call eq_die('No chemical potentials has been assigned this contour: '// & 337 trim(Eq_c(i)%c_io%name)) 338 end if 339 340 ! the id's referring to the chemical potential in the 341 ! mus array. 342 allocate(Eq_c(i)%ID(k)) 343 k = 0 344 do j = 1 , N_mu 345 if ( hasC(mus(j),Eq_c(i)) ) then 346 k = k + 1 347 Eq_c(i)%ID(k) = j 348 349 ! We need to make sure that the number 350 ! of poles is the same if they overlap 351 if ( k > 1 ) then 352 ! We need to ensure that the energy of the poles coincide 353 ! This will happen very rarely!!! 354 r1 = Pi * mus(j)%kT * 2._dp * mus(j)%N_poles 355 c_mu = Eq_c(i)%ID(k-1) 356 r2 = Pi * mus(c_mu)%kT * 2._dp * mus(c_mu)%N_poles 357 358 ! Less than 5 K in difference 359 if ( abs(r1 - r2) > Kelvin ) then 360 call die('If two equilibrium contours & 361 &should overlap, they should have the same & 362 &energy where the line crosses the Fermi energy!') 363 end if 364 end if 365 end if 366 end do 367 368 if ( Eq_c(i)%c_io%N < 1 ) then 369 write(*,*) 'Contour: '//trim(Eq_c(i)%c_io%Name)//' has 0 points.' 370 write(*,*) 'Please ensure at least 1 point in each contour...' 371 call die('Contour number of points, invalid') 372 end if 373 374 ! Allocate the different weights and initialize 375 allocate(Eq_c(i)%c(Eq_c(i)%c_io%N)) 376 Eq_c(i)%c = 0._dp 377 allocate(Eq_c(i)%w(k,Eq_c(i)%c_io%N)) 378 Eq_c(i)%w = 0._dp 379 380 end do 381 382 do i = 1 , N_mu 383 384 call setup_Eq_contour(mus(i)) 385 386 end do 387 388 ! TODO correct empty cycles i.e. when two contours share an energy-point 389 390 contains 391 392 subroutine consecutive_types(c_io,isStart,isTail,mu,kT) 393 type(ts_c_io), intent(in) :: c_io 394 logical, intent(inout) :: isStart, isTail 395 real(dp), intent(in) :: mu, kT 396 397 if ( (leqi(c_io%part,'circle') .or. leqi(c_io%part,'square')) & 398 .and. .not. isStart ) then 399 call eq_die('The circle/square contour must be the first contour type & 400 &as well as connected to other circle/square contours.') 401 end if 402 isStart = leqi(c_io%part,'circle') .or. leqi(c_io%part,'square') 403 404 if ( isStart ) then 405 ! We need to check whether the contour lies well below the 406 ! chemical potential 407 if ( c_io%b > mu - 3._dp * kT ) then 408 call eq_die('The circle/square contour lies too close or across the & 409 &chemical potential. This is not allowed!') 410 end if 411 end if 412 413 if ( .not. leqi(c_io%part,'tail') .and. isTail ) then 414 call eq_die('The tail contour must be the last contour & 415 &as well as connected to other tail contours.') 416 end if 417 isTail = leqi(c_io%part,'tail') 418 419 end subroutine consecutive_types 420 421 end subroutine read_contour_eq_options 422 423 424 ! This routine assures that we have setup all the 425 ! equilibrium contours for the passed electrode 426 subroutine setup_Eq_contour(mu) 427 use fdf, only: leqi 428 type(ts_mu), intent(in) :: mu 429 430 ! Local variables 431 integer :: i, j, l, idx 432 real(dp) :: a, b, R, cR, lift 433 integer :: sq_prev, sq_next 434 435 if ( Eq_segs(mu) == 0 ) & 436 call die('Chemical potential: '//trim(Name(mu))//' has & 437 &no equilibrium contours.') 438 439 ! retrieve circle bounds for the electrode 440 call mu_circle_bounds(mu,a,b) 441 442 ! Calculate the circle entries 443 call calc_Circle(a,b,mu%N_poles,mu%kT,R,cR,lift) 444 445 do i = 1 , Eq_segs(mu) 446 447 ! ensures we progress the contour from start to end 448 ! and not in "random" order 449 idx = get_c_io_index(mu%Eq_seg(i)) 450 451 if ( leqi(Eq_c(idx)%c_io%part,'user') ) then 452 453 call contour_file(Eq_c(idx),mu,lift) 454 455 else if ( leqi(Eq_c(idx)%c_io%part,'circle') ) then 456 457 call contour_Circle(Eq_c(idx),mu,R,cR) 458 459 else if ( leqi(Eq_c(idx)%c_io%part,'square') ) then 460 461 ! find the previous and following square 462 sq_prev = idx 463 sq_next = idx 464 do j = 1 , Eq_segs(mu) 465 l = get_c_io_index(mu%Eq_seg(j)) 466 if ( leqi(Eq_c(l)%c_io%part,'square') ) then 467 if ( l == idx - 1 ) sq_prev = l 468 if ( l == idx + 1 ) sq_next = l 469 end if 470 end do 471 472 call contour_Square(Eq_c(idx),mu, lift, i==1, & 473 Eq_c(sq_prev), Eq_c(sq_next)) 474 475 else if ( leqi(Eq_c(idx)%c_io%part,'line') ) then 476 477 call contour_line(Eq_c(idx),mu,lift) 478 479 else if ( leqi(Eq_c(idx)%c_io%part,'tail') ) then 480 481 call contour_tail(Eq_c(idx),mu,lift) 482 483 else if ( leqi(Eq_c(idx)%c_io%part,'pole') ) then 484 485 ! the poles all have the same weight (Pi*kT*2) 486 call contour_poles(Eq_c(idx),Eq_c(idx)%c_io%d,mu%kT) 487 488 else if ( leqi(Eq_c(idx)%c_io%part,'cont-frac') ) then 489 490 ! the poles all have the same weight (Pi*kT*2) 491 call contour_continued_fraction(Eq_c(idx),mu,lift) 492 493 else 494 495 call die('Unrecognized contour type for the & 496 &equilibrium part.') 497 498 end if 499 500 end do 501 502 contains 503 504 subroutine calc_Circle(a,b,N_poles,kT,R,cR,lift) 505 use units, only : Pi 506 real(dp), intent(in) :: a,b ! the circle bounds 507 integer, intent(in) :: N_poles ! number of poles the 'b' point is lifted 508 real(dp), intent(in) :: kT 509 real(dp), intent(out) :: R, cR, lift ! the radius, center(real) 510 511 ! local variables 512 real(dp) :: alpha 513 514 ! The poles are @ \pi kT, 3\pi kT, 5\pi kT, ... 515 ! Middle point are @ 2\pi kT, 4\pi kT, 6\pi kT 516 lift = Pi * kT * 2._dp * N_poles 517 ! this means that we place the line contour right in the middle of two poles! 518 519 cR = b - a 520 ! the angle between the axis and the line from the start 521 ! of the circle to the end of the circle contour 522 alpha = datan( lift / cR ) 523 524 ! the radius can be calculated using two triangles in the circle 525 ! there is no need to use the cosine-relations 526 R = 0.5_dp * cR / cos(alpha) ** 2 527 528 ! the real-axis center 529 cR = a + R 530 531 end subroutine calc_Circle 532 533 ! retrieve the bounds of the circle contour 534 ! this allows to split the circle integral into as many parts as needed 535 subroutine mu_circle_bounds(mu,a,b) 536 type(ts_mu), intent(in) :: mu 537 real(dp), intent(out) :: a,b 538 if ( Eq_segs(mu) < 1 ) call die('Error Eq_seg CB') 539 idx = get_c_io_index(mu%Eq_seg(1)) 540 a = Eq_c(idx)%c_io%a 541 b = Eq_c(idx)%c_io%b 542 do i = 1 , Eq_segs(mu) 543 idx = get_c_io_index(mu%Eq_seg(i)) 544 if ( Eq_c(idx)%c_io%type == 'eq' .and. & 545 leqi(Eq_c(idx)%c_io%part, 'circle') ) then 546 b = Eq_c(idx)%c_io%b 547 end if 548 end do 549 end subroutine Mu_circle_bounds 550 551 end subroutine setup_Eq_contour 552 553 subroutine ID2idx_cw(c,ID,idx) 554 type(ts_cw), intent(in) :: c 555 integer, intent(in) :: ID 556 integer, intent(out) :: idx 557 do idx = 1 , size(c%ID) 558 if ( c%ID(idx) == ID ) return 559 end do 560 idx = -1 561 end subroutine ID2idx_cw 562 563 subroutine ID2idx_cidx(c,ID,idx) 564 type(ts_c_idx), intent(in) :: c 565 integer, intent(in) :: ID 566 integer, intent(out) :: idx 567 if ( c%idx(1) /= CONTOUR_EQ ) call die('Could not locate ID') 568 call ID2idx_cw(Eq_c(c%idx(2)),ID,idx) 569 end subroutine ID2idx_cidx 570 571 pure subroutine c2weight_eq(c,idx,k,W,ZW) 572 type(ts_c_idx), intent(in) :: c 573 integer, intent(in) :: idx 574 real(dp), intent(in) :: k 575 complex(dp), intent(out) :: W,ZW 576 577 if ( c%idx(1) /= CONTOUR_EQ ) then 578 W = 0._dp 579 ZW = 0._dp 580 return 581 end if 582 583 W = k * Eq_c(c%idx(2))%w(idx,c%idx(3)) 584 ZW = W * Eq_c(c%idx(2))%c(c%idx(3)) 585 586 end subroutine c2weight_eq 587 588 pure function c2energy(c) result(Z) 589 type(ts_c_idx), intent(in) :: c 590 complex(dp) :: Z 591 592 Z = Eq_c(c%idx(2))%c(c%idx(3)) 593 594 end function c2energy 595 596 subroutine contour_Circle(c,mu,R,cR) 597 use fdf, only: leqi 598 use m_integrate 599 use m_gauss_quad 600 type(ts_cw), intent(inout) :: c 601 type(ts_mu), intent(in) :: mu 602 real(dp), intent(in) :: R, cR 603 604 ! local variables 605 character(len=c_N) :: tmpC 606 logical :: set_c 607 complex(dp) :: ztmp 608 integer :: i, idx 609 real(dp) :: a,b, tmp 610 real(dp), allocatable :: ce(:), cw(:) 611 612 if ( .not. leqi(c%c_io%part,'circle') ) & 613 call die('Contour is not a circle') 614 615 ! notice the switch (the circle has increasing contour 616 ! from right to left) 617 call calc_angle(c%c_io%a,R,cR,a) 618 call calc_angle(c%c_io%b,R,cR,b) 619 if ( b > a ) then 620 call die('Contour equilibrium circle went wrong') 621 end if 622 623 allocate(ce(c%c_io%N)) 624 allocate(cw(c%c_io%N)) 625 626 select case ( method(c%c_io) ) 627 case ( CC_G_LEGENDRE ) 628 629 if ( c_io_has_opt(c%c_io,'right') ) then 630 631 deallocate(ce,cw) 632 allocate(ce(c%c_io%N*2)) 633 allocate(cw(c%c_io%N*2)) 634 635 call Gauss_Legendre_Rec(2*c%c_io%N,0,2._dp*a-b,b,ce,cw) 636 637 do i = 1 ,c%c_io%N 638 ce(i) = ce(i+c%c_io%N) 639 cw(i) = cw(i+c%c_io%N) 640 end do 641 642 else if ( c_io_has_opt(c%c_io,'left') ) then 643 644 deallocate(ce,cw) 645 allocate(ce(c%c_io%N*2)) 646 allocate(cw(c%c_io%N*2)) 647 648 call Gauss_Legendre_Rec(2*c%c_io%N,0,a,2._dp*b-a,ce,cw) 649 650 else 651 652 call Gauss_Legendre_Rec(c%c_io%N,0,a,b,ce,cw) 653 654 end if 655 656 case ( CC_TANH_SINH ) 657 658 ! we should also gain an option for this 659 if ( c_io_has_opt(c%c_io,'precision') ) then 660 tmpC = c_io_get_opt(c%c_io,'precision') 661 read(tmpC,'(g20.10)') tmp 662 else 663 tmp = 2.e-2_dp * abs(b-a) / real(c%c_io%N,dp) 664 write(tmpC,'(g20.10)') tmp 665 call c_io_add_opt(c%c_io,'precision',tmpC) 666 end if 667 668 if ( c_io_has_opt(c%c_io,'right') ) then 669 deallocate(ce,cw) 670 allocate(ce(c%c_io%N*2)) 671 allocate(cw(c%c_io%N*2)) 672 673 call TanhSinh_Exact(c%c_io%N*2,ce,cw,2._dp*a-b,b, p=tmp) 674 675 do i = 1 ,c%c_io%N 676 ce(i) = ce(i+c%c_io%N) 677 cw(i) = cw(i+c%c_io%N) 678 end do 679 680 else if ( c_io_has_opt(c%c_io,'left') ) then 681 682 deallocate(ce,cw) 683 allocate(ce(c%c_io%N*2)) 684 allocate(cw(c%c_io%N*2)) 685 686 call TanhSinh_Exact(c%c_io%N*2,ce,cw,a,2._dp*b-a, p=tmp) 687 688 else 689 690 call TanhSinh_Exact(c%c_io%N,ce,cw,a,b, p=tmp) 691 692 end if 693 694 case ( CC_BOOLE_MIX ) 695 696 call Booles_Simpson_38_3_rule(c%c_io%N, ce, cw, a, b) 697 698 case ( CC_SIMP_MIX ) 699 700 call Simpson_38_3_rule(c%c_io%N,ce,cw,a,b) 701 702 case ( CC_MID ) 703 704 call Mid_Rule(c%c_io%N,ce,cw,a,b) 705 706 case default 707 call die('Unknown method for the circle integral, please correct') 708 end select 709 710 ! I know this is "bad" practice, however, zero is a well-defined quantity in FPU 711 set_c = sum(abs(c%c(:))) == 0._dp 712 713 ! get the index in the ID array (same index in w-array) 714 call ID2idx(c,mu%ID,idx) 715 716 do i = 1 , c%c_io%N 717 ztmp = R * cdexp(cmplx(0._dp,ce(i), dp)) 718 719 if ( set_c ) then 720 c%c(i) = cmplx(cR,0._dp, dp) + ztmp 721 else 722 if ( abs(c%c(i) - (cmplx(cR,0._dp, dp) + ztmp) ) > 1.e-10_dp ) then 723 print *,c%c(i),cmplx(cR,0._dp,dp) + ztmp 724 call die('contours does not match') 725 end if 726 end if 727 728 ! Factor i, comes from Ed\theta=dE=iR e^{i\theta} 729 c%w(idx,i) = cw(i) * nf((cmplx(cR,0._dp,dp)+ztmp-mu%mu)/mu%kT) & 730 * cmplx(0._dp,1._dp,dp) * ztmp 731 732 end do 733 734 deallocate(ce,cw) 735 736 contains 737 738 subroutine calc_angle(a,R,cR,angle) 739 use units, only : Pi 740 real(dp), intent(in) :: a,R,cR 741 real(dp), intent(out) :: angle 742 real(dp) :: E 743 744 E = abs(cR - a) 745 746 ! the integration angles 747 if ( abs(E - R) < 1.e-8_dp ) then 748 ! we are at the radius of the circle 749 angle = Pi 750 else 751 angle = asin(E/R) 752 ! correct for left/right side of radius 753 if ( a < cR ) then 754 angle = Pi * 0.5_dp + angle 755 else if ( a > cR ) then 756 angle = Pi * 0.5_dp - angle 757 end if 758 end if 759 760 end subroutine calc_angle 761 762 end subroutine contour_Circle 763 764 subroutine contour_Square(c,mu,lift, first, prev, next) 765 use fdf, only: leqi 766 use m_integrate 767 use m_gauss_quad 768 type(ts_cw), intent(inout) :: c 769 type(ts_mu), intent(in) :: mu 770 real(dp), intent(in) :: lift 771 ! Whether this is the first square contour (starts from 0) 772 logical, intent(in) :: first 773 ! the previous and following square contour 774 type(ts_cw), intent(inout) :: prev, next 775 776 ! local variables 777 character(len=c_N) :: tmpC 778 logical :: set_c 779 complex(dp) :: ztmp(2) 780 ! Segments 781 integer :: N_seg, Ni(3) 782 real(dp) :: li(3) 783 integer :: i, j, idx 784 ! The previous, current and next imaginary part 785 real(dp) :: im(3) 786 real(dp) :: a, b 787 real(dp), allocatable :: ce(:), cw(:) 788 789 if ( .not. leqi(c%c_io%part,'square') ) & 790 call die('Contour is not a square') 791 792 if ( first .and. .not. (prev%c_io .eq. c%c_io) ) then 793 call die('Error in setting up square contour.') 794 end if 795 796 ! We need to split the square quadrature up into different parts 797 ! Specifically this is a requirement as the Gaussian quadratures 798 ! does not have an equi spaced distance between abscissas. 799 ! Thus when climbing/descending on the complex contour we find different 800 ! prefactors which is sub-optimal. 801 802 ! Thus we figure out the number of segments and divide it into those 803 ! partitions 804 805 ! Create parametrization lengths 806 if ( first ) then 807 im(1) = 0._dp 808 else 809 im(1) = read_eta(prev, lift) 810 end if 811 im(2) = read_eta(c, lift) 812 if ( next%c_io .eq. c%c_io ) then 813 im(3) = lift 814 N_seg = 3 815 else 816 ! this is because the following square 817 ! contour steps up! 818 im(3) = im(2) 819 N_seg = 2 820 end if 821 822 if ( im(2) < 0._dp ) then 823 print *,'final imaginary part: ',im(2) 824 call die("Your square contour crosses the real axis. & 825 &This is not allowed") 826 end if 827 828 ! Allocate total weights 829 allocate(ce(c%c_io%N)) 830 allocate(cw(c%c_io%N)) 831 ! I know this is "bad" practice, however, zero is a well-defined quantity in FPU 832 set_c = sum(abs(c%c(:))) == 0._dp 833 call ID2idx(c,mu%ID,idx) 834 835 ! the real axis start/end 836 a = c%c_io%a 837 b = c%c_io%b 838 839 ! Now 'N_seg' contains number of segments 840 ! dL1 = abs(im(2) - im(1)) 841 li(1) = im(2) - im(1) 842 ! dL2 = b - a 843 li(2) = b - a 844 ! dL3 = abs(im(3) - im(2)) ! which may be zero... 845 li(3) = im(3) - im(2) 846 847 ! Ensure enough points for remaining lines 848 Ni = 2 849 call n_distribute(c%c_io%N, N_seg, li, Ni) 850 if ( any(Ni < 2) ) then 851 call die('Could not divide contour points appropriately. & 852 &Choose a different contour method or increase contour points.') 853 end if 854 855 ! Get first segment 856 call get_abscissas(Ni(1),0._dp,li(1), ce(1), cw(1)) 857 ! translate to correct contours 858 do i = 1 , Ni(1) 859 ztmp(1) = cmplx( a, im(1) + ce(i), dp) 860 ztmp(2) = cmplx( 0._dp , cw(i), dp) 861 call set_abscissas(c%c(i),c%w(idx,i), ztmp) 862 end do 863 864 ! Get first segment 865 call get_abscissas(Ni(2),0._dp,li(2), ce(1), cw(1)) 866 ! translate to correct contours 867 j = Ni(1) 868 do i = 1 , Ni(2) 869 ! Convert to parameterisation 870 ztmp(1) = cmplx( a + ce(i) , im(2), dp) 871 ztmp(2) = cmplx( cw(i) , 0._dp, dp) 872 call set_abscissas(c%c(j+i),c%w(idx,j+i), ztmp) 873 end do 874 875 if ( N_seg == 3 ) then 876 ! Get first segment 877 call get_abscissas(Ni(3),0._dp,li(3), ce(1), cw(1)) 878 ! translate to correct contours 879 j = Ni(1) + Ni(2) 880 do i = 1 , Ni(3) 881 ztmp(1) = cmplx( b, im(2) + ce(i), dp) 882 ztmp(2) = cmplx( 0._dp , cw(i), dp) 883 call set_abscissas(c%c(j+i),c%w(idx,j+i), ztmp) 884 end do 885 end if 886 887 deallocate(ce,cw) 888 889 contains 890 891 subroutine get_abscissas(N,a,b,ce,cw) 892 integer, intent(in) :: N 893 real(dp), intent(in) :: a, b 894 real(dp), intent(inout) :: ce(N), cw(N) 895 real(dp), allocatable :: tce(:), tcw(:) 896 real(dp) :: tmp 897 898 select case ( method(c%c_io) ) 899 case ( CC_G_LEGENDRE ) 900 901 if ( c_io_has_opt(c%c_io,'right') ) then 902 903 allocate(tce(N*2),tcw(N*2)) 904 905 call Gauss_Legendre_Rec(2*N,0,2._dp*a-b,b,tce,tcw) 906 907 do i = 1 ,N 908 ce(i) = tce(i+N) 909 cw(i) = tcw(i+N) 910 end do 911 912 deallocate(tce,tcw) 913 914 else if ( c_io_has_opt(c%c_io,'left') ) then 915 916 allocate(tce(N*2),tcw(N*2)) 917 918 call Gauss_Legendre_Rec(2*N,0,a,2._dp*b-a,tce,tcw) 919 920 do i = 1 ,N 921 ce(i) = tce(i) 922 cw(i) = tcw(i) 923 end do 924 925 deallocate(tce,tcw) 926 927 else 928 929 call Gauss_Legendre_Rec(N,0,a,b,ce,cw) 930 931 end if 932 933 case ( CC_TANH_SINH ) 934 935 ! we should also gain an option for this 936 if ( c_io_has_opt(c%c_io,'precision') ) then 937 tmpC = c_io_get_opt(c%c_io,'precision') 938 read(tmpC,'(g20.10)') tmp 939 else 940 tmp = 2.e-2_dp * abs(b-a) / real(N,dp) 941 write(tmpC,'(g20.10)') tmp 942 call c_io_add_opt(c%c_io,'precision',tmpC) 943 end if 944 945 if ( c_io_has_opt(c%c_io,'right') ) then 946 947 allocate(tce(N*2),tcw(N*2)) 948 949 call TanhSinh_Exact(N*2,tce,tcw,2._dp*a-b,b, p=tmp) 950 951 do i = 1 ,N 952 ce(i) = tce(i+N) 953 cw(i) = tcw(i+N) 954 end do 955 956 else if ( c_io_has_opt(c%c_io,'left') ) then 957 958 allocate(tce(N*2),tcw(N*2)) 959 960 call TanhSinh_Exact(N*2,tce,tcw,a,2._dp*b-a, p=tmp) 961 962 do i = 1 ,N 963 ce(i) = tce(i) 964 cw(i) = tcw(i) 965 end do 966 967 deallocate(tce,tcw) 968 969 else 970 971 call TanhSinh_Exact(N,ce,cw,a,b, p=tmp) 972 973 end if 974 975 case ( CC_BOOLE_MIX ) 976 977 call Booles_Simpson_38_3_rule(N, ce, cw, a, b) 978 979 case ( CC_SIMP_MIX ) 980 981 call Simpson_38_3_rule(N,ce,cw,a,b) 982 983 case ( CC_MID ) 984 985 call Mid_Rule(N,ce,cw,a,b) 986 987 case default 988 call die('Unknown method for the square integral, please correct') 989 end select 990 991 end subroutine get_abscissas 992 993 subroutine set_abscissas(e,w,ztmp) 994 complex(dp), intent(inout) :: e, w 995 complex(dp), intent(inout) :: ztmp(2) 996 real(dp) :: tmp 997 998 if ( set_c ) then 999 e = ztmp(1) 1000 else 1001 if ( abs(e - ztmp(1) ) > 1.e-10_dp ) then 1002 print *,e,ztmp(1) 1003 call die('contours does not match') 1004 end if 1005 end if 1006 1007 ! the parameterisation does not alter the fermi function 1008 ! importantly we define the band-bottom to have a vanishing 1009 ! fermi function quantity. 1010 tmp = (real(ztmp(1),dp) - mu%mu) / mu%kT 1011 ztmp(1) = (ztmp(1) - mu%mu) / mu%kT 1012 w = nf(ztmp(1)) * ztmp(2) 1013 1014 end subroutine set_abscissas 1015 1016 function read_eta(c, lift) result(eta) 1017 use parse, only: parsed_line, digest, destroy 1018 type(ts_cw), intent(in) :: c 1019 real(dp), intent(in) :: lift 1020 real(dp) :: eta 1021 1022 type(parsed_line), pointer :: pline 1023 character(len=c_N) :: tmpC 1024 1025 if ( c_io_has_opt(c%c_io,'eta-add') ) then 1026 tmpC = c_io_get_opt(c%c_io,'eta-add') 1027 pline => digest(tmpC) 1028 eta = lift + ts_c_bphysical(pline,0,'Ry') 1029 call destroy(pline) 1030 else if ( c_io_has_opt(c%c_io,'eta') ) then 1031 tmpC = c_io_get_opt(c%c_io,'eta') 1032 pline => digest(tmpC) 1033 eta = ts_c_bphysical(pline,0,'Ry') 1034 call destroy(pline) 1035 else 1036 ! We default to add 1.5 Ry to the square 1037 ! to keep it well of the axis 1038 eta = lift + 1.5_dp 1039 end if 1040 1041 end function read_eta 1042 1043 end subroutine contour_Square 1044 1045 subroutine contour_line(c,mu,Eta) 1046 use fdf, only: leqi 1047 use m_integrate 1048 use m_gauss_quad 1049 1050 type(ts_cw), intent(inout) :: c 1051 type(ts_mu), intent(in) :: mu 1052 ! The lifting into the complex plane 1053 real(dp), intent(in) :: Eta 1054 1055 ! local variables 1056 character(len=c_N) :: tmpC 1057 logical :: set_c 1058 integer :: i, idx 1059 real(dp) :: a,b, tmp 1060 real(dp), allocatable :: ce(:), cw(:) 1061 1062 if ( .not. leqi(c%c_io%part,'line') ) & 1063 call die('Contour is not a line') 1064 1065 ! get bounds 1066 a = c%c_io%a 1067 b = c%c_io%b 1068 1069 allocate(ce(c%c_io%N)) 1070 allocate(cw(c%c_io%N)) 1071 1072 select case ( method(c%c_io) ) 1073 case ( CC_MID ) 1074 1075 call Mid_Rule(c%c_io%N,ce,cw,a,b) 1076 1077 case ( CC_SIMP_MIX ) 1078 1079 call Simpson_38_3_rule(c%c_io%N,ce,cw,a,b) 1080 1081 case ( CC_BOOLE_MIX ) 1082 1083 call Booles_Simpson_38_3_rule(c%c_io%N,ce,cw,a,b) 1084 1085 case ( CC_G_LEGENDRE ) 1086 1087 if ( c_io_has_opt(c%c_io,'right') ) then 1088 1089 deallocate(ce,cw) 1090 allocate(ce(c%c_io%N*2)) 1091 allocate(cw(c%c_io%N*2)) 1092 1093 call Gauss_Legendre_Rec(c%c_io%N*2,0,2._dp*a-b,b,ce,cw) 1094 1095 do i = 1 ,c%c_io%N 1096 ce(i) = ce(i+c%c_io%N) 1097 cw(i) = cw(i+c%c_io%N) 1098 end do 1099 1100 else if ( c_io_has_opt(c%c_io,'left') ) then 1101 1102 deallocate(ce,cw) 1103 allocate(ce(c%c_io%N*2)) 1104 allocate(cw(c%c_io%N*2)) 1105 1106 call Gauss_Legendre_Rec(c%c_io%N*2,0,a,2._dp*b-a,ce,cw) 1107 1108 else 1109 1110 call Gauss_Legendre_Rec(c%c_io%N,0,a,b,ce,cw) 1111 1112 end if 1113 1114 case ( CC_TANH_SINH ) 1115 1116 ! we should also gain an option for this 1117 if ( c_io_has_opt(c%c_io,'precision') ) then 1118 tmpC = c_io_get_opt(c%c_io,'precision') 1119 read(tmpC,'(g20.10)') tmp 1120 else 1121 tmp = 2.e-2_dp * abs(b-a) / real(c%c_io%N,dp) 1122 write(tmpC,'(g20.10)') tmp 1123 call c_io_add_opt(c%c_io,'precision',tmpC) 1124 end if 1125 1126 if ( c_io_has_opt(c%c_io,'right') ) then 1127 deallocate(ce,cw) 1128 allocate(ce(c%c_io%N*2)) 1129 allocate(cw(c%c_io%N*2)) 1130 1131 call TanhSinh_Exact(c%c_io%N*2,ce,cw,2._dp*a-b,b, p=tmp) 1132 1133 do i = 1 ,c%c_io%N 1134 ce(i) = ce(i+c%c_io%N) 1135 cw(i) = cw(i+c%c_io%N) 1136 end do 1137 1138 else if ( c_io_has_opt(c%c_io,'left') ) then 1139 1140 deallocate(ce,cw) 1141 allocate(ce(c%c_io%N*2)) 1142 allocate(cw(c%c_io%N*2)) 1143 1144 call TanhSinh_Exact(c%c_io%N*2,ce,cw,a,2._dp*b-a, p=tmp) 1145 1146 else 1147 1148 call TanhSinh_Exact(c%c_io%N,ce,cw,a,b, p=tmp) 1149 1150 end if 1151 1152 case ( CC_USER ) 1153 1154 ! Read the file information 1155 call contour_file(c,mu,Eta) 1156 1157 case default 1158 write(*,*) 'Method for contour ',trim(c%c_io%name), & 1159 ' could not be deciphered: ', c%c_io%method 1160 call die('Could not determine the line-integral') 1161 end select 1162 1163 ! get the index in the ID array (same index in w-array) 1164 call ID2idx(c,mu%ID,idx) 1165 1166 if ( method(c%c_io) == CC_USER ) then 1167 1168 do i = 1 , c%c_io%N 1169 c%w(idx,i) = c%w(idx,i) * nf((real(c%c(i), dp) - mu%mu) / mu%kT) 1170 end do 1171 1172 else 1173 1174 ! I know this is "bad" practice, however, zero is a well-defined quantity. 1175 set_c = sum(abs(c%c(:))) == 0._dp 1176 1177 do i = 1 , c%c_io%N 1178 if ( set_c ) then 1179 c%c(i) = cmplx(ce(i),Eta, dp) 1180 else 1181 if ( abs(c%c(i) - cmplx(ce(i),Eta, dp)) > 1.e-10_dp ) then 1182 call die('contour_line: Error on contour match') 1183 end if 1184 end if 1185 1186 c%w(idx,i) = cw(i) * nf((ce(i) - mu%mu) / mu%kT) 1187 1188 end do 1189 1190 end if 1191 1192 deallocate(ce,cw) 1193 1194 end subroutine contour_line 1195 1196 subroutine contour_tail(c,mu,Eta) 1197 use fdf, only: leqi 1198 use m_gauss_fermi_inf 1199 use m_gauss_fermi_30 1200 use m_gauss_fermi_28 1201 use m_gauss_fermi_26 1202 use m_gauss_fermi_24 1203 use m_gauss_fermi_22 1204 use m_gauss_fermi_20 1205 use m_gauss_fermi_19 1206 use m_gauss_fermi_18 1207 use m_gauss_fermi_17 1208 1209 type(ts_cw), intent(inout) :: c 1210 type(ts_mu), intent(in) :: mu 1211 ! This describes the lifting of the tail integral into the complex plane 1212 real(dp), intent(in) :: Eta 1213 1214 ! local variables 1215 integer :: idx, offset, infinity 1216 real(dp) :: a,b 1217 real(dp), allocatable :: ce(:), cw(:) 1218 1219 if ( .not. leqi(c%c_io%part,'tail') ) & 1220 call die('Contour is not a tail contour') 1221 1222 ! get bounds 1223 a = c%c_io%a 1224 b = c%c_io%b 1225 1226 allocate(ce(c%c_io%N)) 1227 allocate(cw(c%c_io%N)) 1228 1229 select case ( method(c%c_io) ) 1230 case ( CC_G_NF_MIN:CC_G_NF_MAX ) 1231 1232 offset = nint((c%c_io%a-mu%mu)/mu%kT) 1233 1234 if ( abs(offset * mu%kT - (c%c_io%a-mu%mu)) > 1.e-7_dp ) then 1235 call die('The integer value of the kT offset for the & 1236 &Gauss-Fermi integral is not valid, please check input') 1237 end if 1238 if ( c%c_io%b < 0._dp ) then 1239 if ( c%c_io%b < -30.5_dp * mu%kT ) then 1240 infinity = huge(1) 1241 else 1242 infinity = nint(abs(c%c_io%b-mu%mu)/mu%kT) 1243 end if 1244 else 1245 if ( c%c_io%b > 30.5_dp * mu%kT ) then 1246 infinity = huge(1) 1247 else 1248 infinity = nint(abs(c%c_io%b-mu%mu)/mu%kT) 1249 end if 1250 end if 1251 if ( infinity > 30 ) infinity = huge(1) 1252 1253 ! calculate the offset from the energy chemical potential tail 1254 select case ( infinity ) 1255 case ( huge(1) ) 1256 call GaussFermi_inf(offset,c%c_io%N,ce,cw) 1257 case ( 30 ) 1258 call GaussFermi_30(offset,c%c_io%N,ce,cw) 1259 case ( 28 ) 1260 call GaussFermi_28(offset,c%c_io%N,ce,cw) 1261 case ( 26 ) 1262 call GaussFermi_26(offset,c%c_io%N,ce,cw) 1263 case ( 24 ) 1264 call GaussFermi_24(offset,c%c_io%N,ce,cw) 1265 case ( 22 ) 1266 call GaussFermi_22(offset,c%c_io%N,ce,cw) 1267 case ( 20 ) 1268 call GaussFermi_20(offset,c%c_io%N,ce,cw) 1269 case ( 19 ) 1270 call GaussFermi_19(offset,c%c_io%N,ce,cw) 1271 case ( 18 ) 1272 call GaussFermi_18(offset,c%c_io%N,ce,cw) 1273 case ( 17 ) 1274 call GaussFermi_17(offset,c%c_io%N,ce,cw) 1275 case default 1276 call die('Unknown tail integral ending') 1277 end select 1278 1279 ! we might as well correct the method 1280 1281 ! the Gauss-Fermi quadrature is wrt. E'->E/kT 1282 ce = ce * mu%kT + mu%mu 1283 cw = cw * mu%kT 1284 1285 call ID2idx(c,mu%ID,idx) 1286 1287 ! move over the weights and the contour values 1288 c%c(:) = cmplx(ce,Eta, dp) 1289 c%w(idx,:) = cw 1290 1291 case default 1292 1293 ! we revert so that we can actually use the line-integral 1294 ! The tail and line are equivalent in the sense that the 1295 ! fermi functions are applied to the weights 1296 c%c_io%part = 'line' 1297 1298 call contour_line(c,mu,Eta) 1299 1300 end select 1301 1302 deallocate(ce,cw) 1303 1304 end subroutine contour_tail 1305 1306 1307 ! The residuals of the fermi-function at a real-energy 1308 subroutine contour_poles(c, E, kT) 1309 use units, only : Pi 1310 1311! *********************** 1312! * OUTPUT variables * 1313! *********************** 1314 type(ts_cw), intent(inout) :: c 1315 1316! *********************** 1317! * INPUT variables * 1318! *********************** 1319 real(dp), intent(in) :: E ! at energy 1320 real(dp), intent(in) :: kT ! temperature in Ry 1321 1322! *********************** 1323! * LOCAL variables * 1324! *********************** 1325 integer :: i 1326 1327 ! all pole-weights have the same weight (negative due to contour method) 1328 c%w(:,:) = - cmplx(0._dp, Pi * kT * 2._dp, dp) 1329 ! Residuals 1330 do i = 1 , c%c_io%N 1331 c%c(i) = cmplx(E , Pi * kT * (2._dp*(i-1)+1._dp), dp) 1332 end do 1333 1334 end subroutine contour_poles 1335 1336 1337 subroutine contour_continued_fraction(c,mu,Eta) 1338 use fdf, only: leqi 1339 use units, only: Pi 1340 1341 type(ts_cw), intent(inout) :: c 1342 type(ts_mu), intent(in) :: mu 1343 ! The lifting into the complex plane 1344 real(dp), intent(in) :: Eta 1345 1346 ! local variables 1347 logical :: set_c 1348 integer :: i, idx 1349 complex(dp) :: cc 1350 real(dp), allocatable :: ce(:), cw(:) 1351 1352 if ( .not. leqi(c%c_io%part,'cont-frac') ) & 1353 call die('Contour is not a continued fraction') 1354 1355 allocate(ce(c%c_io%N-1)) 1356 allocate(cw(c%c_io%N-1)) 1357 1358 select case ( method(c%c_io) ) 1359 case ( CC_CONTINUED_FRAC ) 1360 1361 call Ozaki_residue(c%c_io%N-1,ce,cw) 1362 1363 case default 1364 write(*,*) 'Method for contour ',trim(c%c_io%name), & 1365 ' could not be deciphered: ', c%c_io%method 1366 call die('Could not determine the pole-integral') 1367 end select 1368 1369 set_c = sum(abs(c%c(:))) == 0._dp 1370 1371 ! get the index in the ID array (same index in w-array) 1372 call ID2idx(c,mu%ID,idx) 1373 1374 do i = 1 , c%c_io%N - 1 1375 1376 ! Calculate current contour point 1377 cc = cmplx(mu%mu, ce(i) * mu%kT, dp) 1378 1379 if ( set_c ) then 1380 c%c(i) = cc 1381 else 1382 if ( abs(c%c(i) - cc) > 1.e-10_dp ) then 1383 call die('contour_cont_frac: Error on contour match') 1384 end if 1385 end if 1386 1387 ! Extra minus in implementation and Im[] 1388 ! We also divide the weight by Pi in the loop (and it should 1389 ! not exist in the continued fraction scheme) 1390 c%w(idx,i) = cmplx( 0._dp , 2._dp * cw(i) * mu%kT * Pi, dp) 1391 1392 end do 1393 1394 ! The zero'th moment lies infinitely far and is from -inf -- inf 1395 cc = cmplx( mu%mu , c%c_io%a, dp) 1396 1397 if ( set_c ) then 1398 ! The last pole is set 1399 c%c(c%c_io%N) = cc 1400 else 1401 if ( abs(c%c(c%c_io%N) - cc) > 1.e-10_dp ) then 1402 call die('contour_cont_frac: Error on contour match') 1403 end if 1404 end if 1405 1406 ! The zeroth moment (extra minus in implementation and Im[]) 1407 ! w = iR, but from -\Im we get w = R 1408 ! And remove the loop division by Pi 1409 c%w(idx,c%c_io%N) = cmplx( 0.5_dp * c%c_io%a * Pi, 0._dp, dp) 1410 1411 deallocate(ce,cw) 1412 1413 contains 1414 1415 subroutine Ozaki_residue(N,c,w) 1416 ! The number of poles 1417 integer, intent(in) :: N 1418 real(dp), intent(out) :: c(N), w(N) 1419 1420 ! Diagonalization matrices 1421 real(dp), allocatable :: A(:,:), B(:,:), we(:), work(:) 1422 integer :: i, j, nn 1423 1424 real(dp), external :: ddot 1425 1426 ! The last weight of the residue problem is the R->infinity 1427 ! point. Hence we remove one point from the poles. 1428 nn = 2*N 1429 allocate(A(nn,nn), B(nn,nn)) 1430 allocate(we(nn),work(3*nn)) 1431 1432 A = 0._dp 1433 B = 0._dp 1434 do i = 1 , nn - 1 1435 B(i,i) = 2._dp * i - 1._dp 1436 A(i,i+1) = -0.5_dp 1437 A(i+1,i) = -0.5_dp 1438 end do 1439 B(nn,nn) = 2._dp * nn - 1._dp 1440 1441 ! Matrices have been initialized, diagonalize 1442 ! to find residues from eigenvalues 1443 call dsygv(1,'V','U',nn,A,nn,B,nn,we,work,3*nn,i) 1444 if ( i /= 0 ) then 1445 write(*,'(a)')'error in Ozaki diagonalization of weights.' 1446 call die('Error in diagonalization of weights') 1447 end if 1448 1449 1450 ! Transfer weights and eigenvalues 1451 j = nn 1452 do i = 1 , N 1453 1454 ! Eigenvalue 1455 c(i) = 1._dp / we(j) 1456 1457 ! Residual 1458 w(i) = - (A(1,j)*c(i))**2 * 0.25_dp 1459 1460 j = j - 1 1461 1462 end do 1463 1464 ! Checked with Ozaki and got the same 1465 !print '(1000(/,2(tr2,f24.12)))',( (/c(i),w(i)/) ,i=1,N) 1466 1467 deallocate(A,B,we,work) 1468 1469 end subroutine Ozaki_residue 1470 1471 end subroutine contour_continued_fraction 1472 1473 ! This routine will read the contour points from a file 1474 subroutine contour_file(c,mu,Eta) 1475 use m_io, only: io_assign, io_close 1476 use fdf, only: fdf_convfac 1477 1478 type(ts_cw), intent(inout) :: c 1479 type(ts_mu), intent(in) :: mu 1480 ! The lifting into the complex plane 1481 real(dp), intent(in) :: Eta 1482 1483 integer :: iu, iostat, ne, idx 1484 logical :: exist 1485 complex(dp) :: E , W 1486 real(dp) :: rE, iE, rW, iW, conv 1487 character(len=512) :: file, line 1488 character(len=16) :: unit 1489 1490 ! The contour type contains the file name in: 1491 ! c%c_io%cN (weirdly enough) 1492 file = c%c_io%cN 1493 call ID2idx(c,mu%ID,idx) 1494 1495 call io_assign(iu) 1496 1497 ! Open the contour file 1498 inquire(file=trim(file), exist=exist) 1499 if ( .not. exist ) then 1500 call die('The file: '//trim(file)//' could not be found to read contour points!') 1501 end if 1502 1503 ! Open the file 1504 open(iu, file=trim(file), form='formatted', status='old') 1505 1506 ne = 0 1507 ! The default unit is eV. 1508 ! On every line an optional unit-specificer may be used to specify the 1509 ! subsequent lines units (until another unit is specified) 1510 conv = fdf_convfac('eV', 'Ry') 1511 do 1512 ! Now we have the line 1513 read(iu, '(a)', iostat=iostat) line 1514 if ( iostat /= 0 ) exit 1515 if ( len_trim(line) == 0 ) cycle 1516 line = trim(adjustl(line)) 1517 if ( line(1:1) == '#' ) cycle 1518 1519 ! We have a line with energy and weight 1520 ne = ne + 1 1521 ! There are three optional ways of reading this 1522 ! 1. ReE ImE, ReW ImW [unit] 1523 read(line, *, iostat=iostat) rE, iE, rW, iW, unit 1524 if ( iostat == 0 ) then 1525 conv = fdf_convfac(unit, 'Ry') 1526 else 1527 read(line, *, iostat=iostat) rE, iE, rW, iW 1528 end if 1529 if ( iostat == 0 ) then 1530 c%c(ne) = cmplx(rE,iE, dp) * conv 1531 c%w(idx,ne) = cmplx(rW,iW, dp) * conv 1532 cycle 1533 end if 1534 1535 ! 2. ReE ImE, ReW [unit] 1536 iW = 0._dp 1537 read(line, *, iostat=iostat) rE, iE, rW, unit 1538 if ( iostat == 0 ) then 1539 conv = fdf_convfac(unit, 'Ry') 1540 else 1541 read(line, *, iostat=iostat) rE, iE, rW 1542 end if 1543 if ( iostat == 0 ) then 1544 c%c(ne) = cmplx(rE,iE,dp) * conv 1545 c%w(idx,ne) = cmplx(rW,iW,dp) * conv 1546 cycle 1547 end if 1548 1549 ! 3. ReE , ReW [unit] 1550 iE = Eta 1551 iW = 0._dp 1552 read(line, *, iostat=iostat) rE, rW, unit 1553 if ( iostat == 0 ) then 1554 conv = fdf_convfac(unit, 'Ry') 1555 else 1556 read(line, *, iostat=iostat) rE, rW 1557 end if 1558 if ( iostat == 0 ) then 1559 c%c(ne) = cmplx(rE * conv,iE,dp) 1560 c%w(idx,ne) = cmplx(rW,iW,dp) * conv 1561 cycle 1562 end if 1563 1564 call die('Contour file: '//trim(file)//' is not formatted correctly. & 1565 &Please read the documentation!') 1566 1567 end do 1568 1569 call io_close(iu) 1570 1571 if ( c%c_io%N /= ne ) then 1572 call die('Error in reading the contour points from file: '//trim(file)) 1573 end if 1574 1575 end subroutine contour_file 1576 1577 function Eq_E(id,step) result(c) 1578 integer, intent(in) :: id 1579 integer, intent(in), optional :: step 1580 type(ts_c_idx) :: c ! the configuration of the energy-segment 1581 integer :: lstep, i, PN 1582 lstep = 1 1583 if ( present(step) ) lstep = step 1584 PN = N_Eq_E() 1585 if ( id <= PN ) then 1586 c = get_c(id) 1587 return 1588 end if 1589 c = get_c(-1) 1590 i = MOD(PN,lstep) 1591 if ( i /= 0 ) PN = PN + lstep - i 1592 if ( id <= PN ) then 1593 c%exist = .true. 1594 c%fake = .true. 1595 end if 1596 end function Eq_E 1597 1598 function get_c(id) result(c) 1599 integer, intent(in) :: id 1600 type(ts_c_idx) :: c 1601 integer :: i,j,iE 1602 c%exist = .false. 1603 c%fake = .false. 1604 c%e = cmplx(0._dp,0._dp, dp) 1605 c%idx = 0 1606 if ( id < 1 ) return 1607 1608 iE = 0 1609 do j = 1 , N_Eq ! number of contours 1610 if ( iE + Eq_c(j)%c_io%N < id ) then 1611 iE = iE + Eq_c(j)%c_io%N 1612 cycle 1613 end if 1614 i = id - iE 1615 if ( i <= Eq_c(j)%c_io%N ) then 1616 c%exist = .true. 1617 c%e = Eq_c(j)%c(i) 1618 c%idx(1) = CONTOUR_EQ ! designates the equilibrium contours 1619 c%idx(2) = j ! designates the index of the equilibrium contour 1620 c%idx(3) = i ! is the index of the equilibrium contour 1621 return 1622 end if 1623 end do 1624 1625 end function get_c 1626 1627 function N_Eq_E() result(N) 1628 integer :: N, i 1629 N = 0 1630 do i = 1 , N_Eq 1631 N = N + Eq_c(i)%c_io%N 1632 end do 1633 end function N_Eq_E 1634 1635 !< Calculate the linear (total) index of a given index+local contour index 1636 function Eq_linear_index(idx, ie) result(lidx) 1637 !< The index of the equilibrium contour 1638 integer, intent(in) :: idx 1639 !< The index of the energy on the `idx` contour 1640 integer, intent(in) :: ie 1641 !< The global index of the `idx,ie` contour point 1642 integer :: lidx 1643 1644 integer :: i 1645 1646 lidx = ie 1647 do i = 1, idx - 1 1648 lidx = lidx + Eq_c(i)%c_io%N 1649 end do 1650 1651 end function Eq_linear_index 1652 1653 1654 subroutine print_contour_eq_block(prefix) 1655 use fdf, only: leqi 1656 character(len=*), intent(in) :: prefix 1657 integer :: i 1658 1659 do i = 1 , N_Eq 1660 if ( leqi(Eq_io(i)%part,'pole') ) cycle 1661 if ( leqi(Eq_io(i)%part,'cont-frac') ) cycle 1662 call ts_print_contour_block(trim(prefix)//'.Contour.',Eq_io(i)) 1663 end do 1664 end subroutine print_contour_eq_block 1665 1666 subroutine print_contour_eq_options(prefix) 1667 1668 use fdf, only: leqi 1669 use parallel, only : IONode 1670 use units, only : Pi 1671 1672 use m_ts_io_contour 1673 1674 character(len=*), intent(in) :: prefix 1675 character(len=200) :: chars 1676 real(dp) :: tmp 1677 integer :: i, N 1678 type(ts_c_opt_ll), pointer :: opt 1679 1680 if ( .not. IONode ) return 1681 1682 write(*,opt_n) ' ----------------- Contour ----------------- ' 1683 1684 N = 0 1685 write(*,opt_n) ' >> Residual contour << ' 1686 do i = 1 , N_eq 1687 chars = trim(eq_io(i)%part) 1688 if ( .not. (leqi(chars,'pole') .or. & 1689 leqi(chars,'cont-frac')) ) cycle 1690 N = N + 1 1691 if ( leqi(chars,'cont-frac') ) then 1692 call write_e('Continued fraction chemical potential',eq_io(i)%d) 1693 else 1694 call write_e('Pole chemical potential',eq_io(i)%d) 1695 end if 1696 call write_e(' Chemical potential temperature',eq_io(i)%b, & 1697 unit = 'K') 1698 write(*,opt_int) ' Number of poles',eq_io(i)%N 1699 if ( .not. leqi(chars,'cont-frac') ) then 1700 ! Calculate energy of middle pole-point 1701 tmp = Pi * eq_io(i)%b * 2._dp * eq_io(i)%N 1702 call write_e(' Top energy point',tmp) 1703 end if 1704 end do 1705 1706 if ( N < N_Eq ) & 1707 write(*,opt_n) ' >> Equilibrium contour << ' 1708 do i = 1 , N_eq 1709 if ( leqi(eq_io(i)%part,'pole') ) cycle 1710 if ( leqi(eq_io(i)%part,'cont-frac') ) cycle 1711 chars = ' '//trim(eq_io(i)%part) 1712 ! the starred contours are "fakes" 1713 if ( eq_io(i)%name(1:1) == '*' ) then 1714 write(*,opt_c) 'Contour name',trim(prefix)//'.Contour.'//trim(eq_io(i)%name(2:)) 1715 else 1716 write(*,opt_c) 'Contour name',trim(prefix)//'.Contour.'//trim(eq_io(i)%name) 1717 end if 1718 1719 call write_e(trim(chars)//' contour E_min',eq_io(i)%a) 1720 call write_e(trim(chars)//' contour E_max',eq_io(i)%b) 1721 write(*,opt_int) trim(chars)//' contour points',eq_io(i)%N 1722 write(*,opt_c) trim(chars)//' contour method', & 1723 trim(longmethod2str(eq_io(i))) 1724 opt => eq_io(i)%opt 1725 do while ( associated(opt) ) 1726 if ( len_trim(opt%val) > 0 ) then 1727 write(*,opt_cc) ' Option for contour method',trim(opt%opt),trim(opt%val) 1728 else 1729 write(*,opt_c) ' Option for contour method',trim(opt%opt) 1730 end if 1731 opt => opt%next 1732 end do 1733 end do 1734 1735 end subroutine print_contour_eq_options 1736 1737 function get_c_io_index(Name) result(idx) 1738 character(len=*), intent(in) :: Name 1739 integer :: idx 1740 do idx = 1 , N_Eq 1741 if ( trim(name) == trim(Eq_io(idx)%name) ) return 1742 end do 1743 idx = 0 1744 end function get_c_io_index 1745 1746 subroutine io_contour_eq(mus,slabel,suffix) 1747 type(ts_mu), intent(in) :: mus(:) 1748 character(len=*), intent(in) :: slabel 1749 character(len=*), intent(in), optional :: suffix 1750 1751 integer :: i 1752 1753 do i = 1 , size(mus) 1754 call io_contour_eq_mu(mus(i),slabel,suffix) 1755 end do 1756 1757 end subroutine io_contour_eq 1758 1759 subroutine io_contour_eq_mu(mu,slabel,suffix) 1760 use parallel, only : IONode 1761 use fdf, only : leqi 1762 use units, only: eV, Kelvin 1763 type(ts_mu), intent(in) :: mu 1764 character(len=*), intent(in) :: slabel 1765 character(len=*), intent(in), optional :: suffix 1766 1767! ********************* 1768! * LOCAL variables * 1769! ********************* 1770 character(len=256) :: fname 1771 integer :: i, unit, idx 1772 type(ts_c_idx) :: cidx 1773 1774 if ( .not. IONode ) return 1775 1776 if ( present(suffix) ) then 1777 fname = trim(slabel)//trim(suffix) 1778 else 1779 fname = trim(slabel)//'.TSCCEQ-'//trim(Name(mu)) 1780 end if 1781 1782 call io_assign( unit ) 1783 open( unit, file=fname, status='unknown' ) 1784 write(unit,'(a)') '# Contour path for the equilibrium contour segment.' 1785 write(unit,'(a)') '# This segment belongs to the chemical potential: '//trim(Name(mu)) 1786 write(unit,'(a)') '# Chemical potential:' 1787 if ( mu%mu < 0._dp ) then 1788 write(unit, '(a,g10.4,a)')'# - ', -mu%mu / eV, ' eV' 1789 else 1790 write(unit, '(a,g10.4,a)')'# + ', mu%mu / eV, ' eV' 1791 end if 1792 write(unit,'(a)') '# Electronic temperature:' 1793 write(unit, '(a,g10.4,a)')'# ', mu%kT / Kelvin, ' K' 1794 write(unit,'(a,a24,3(tr1,a25))') '#','Re(c) [eV]','Im(c) [eV]','Re(w) [eV]','Im(w) [eV]' 1795 1796 cidx%idx(1) = CONTOUR_EQ 1797 do i = 1 , Eq_segs(mu) 1798 1799 cidx%idx(2) = get_c_io_index(mu%Eq_seg(i)) 1800 if ( cidx%idx(2) < 1 ) call die('io_contour_eq_mu: Error in code, C-ID') 1801 call ID2idx(Eq_c(cidx%idx(2)),mu%ID,idx) 1802 if ( idx < 1 ) call die('io_contour_eq_mu: Error in code') 1803 1804 !print *,trim(Name(mu)),mu%Eq_seg(i) 1805 1806 ! we have now retrieved the electrode index in the contour part 1807 ! write it out 1808 call io_contour_c(unit,cidx,idx) 1809 1810 end do 1811 1812 call io_close( unit ) 1813 1814 end subroutine io_contour_eq_mu 1815 1816! Write out the contour to a contour file 1817 subroutine io_contour_c(unit,cidx,idx) 1818 use parallel, only : IONode 1819 use units, only : Pi, eV 1820 use fdf, only: leqi 1821 integer, intent(in) :: unit 1822 type(ts_c_idx), intent(inout) :: cidx 1823 integer, intent(in) :: idx 1824 1825! ********************* 1826! * LOCAL variables * 1827! ********************* 1828 integer :: i 1829 logical :: is_cont_frac 1830 type(ts_cw), pointer :: c 1831 complex(dp) :: W, ZW 1832 1833 if ( .not. IONode ) return 1834 c => Eq_c(cidx%idx(2)) 1835 1836 is_cont_frac = leqi(c%c_io%part,'cont-frac') 1837 1838 do i = 1 , size(c%c) 1839 cidx%e = c%c(i) 1840 cidx%idx(3) = i 1841 call c2weight_eq(cidx,idx,1._dp,W,ZW) 1842 if ( is_cont_frac ) then 1843 write(unit,'(4(e25.17,tr1))') c%c(i)/eV, W/(eV*Pi) 1844 else 1845 write(unit,'(4(e25.17,tr1))') c%c(i)/eV, W/eV 1846 end if 1847 end do 1848 1849 end subroutine io_contour_c 1850 1851 subroutine eq_die(msg) 1852 character(len=*), intent(in) :: msg 1853 write(*,*) 'Killing... printing out so-far gathered information' 1854 call print_contour_eq_options('TS') 1855 call print_contour_eq_block('TS') 1856 call die(msg) 1857 end subroutine eq_die 1858 1859end module m_ts_contour_eq 1860