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_chem_pot 13 14 use precision, only : dp 15 16 use units, only : Pi, eV, Kelvin 17 18 use m_ts_io_ctype, only : C_N_NAME_LEN 19 20 implicit none 21 22 integer, public, parameter :: NAME_MU_LEN = 32 23 24 real(dp), public, parameter :: mu_same = 1.e-8_dp 25 integer, private, parameter :: def_poles = 8 26 27 type :: ts_mu 28 ! name of the chemical potential 29 character(len=NAME_MU_LEN) :: name = ' ' 30 ! ID 31 integer :: ID = 0 32 ! Number of poles for the equilibrium contour 33 integer :: N_poles = 0 34 ! the chemical potential 35 real(dp) :: mu = 0._dp 36 character(len=NAME_MU_LEN) :: cmu 37 ! the temperature associated with this chemical potential 38 real(dp) :: kT = 0._dp 39 character(len=NAME_MU_LEN) :: ckT 40 ! number of electrodes having this chemical potential 41 integer :: N_El = 0 42 ! array of electrode indices (conforming with the Elecs-array) 43 integer, pointer :: el(:) => null() 44 ! We must have a container which determines the contour segments 45 ! that gets attached to the chemical potential 46 character(len=C_N_NAME_LEN), allocatable :: Eq_seg(:) 47 end type ts_mu 48 public :: ts_mu 49 50 interface hasC 51 module procedure hasCio 52 module procedure hasCeq 53 end interface hasC 54 public :: hasC 55 56 interface hasEl 57 module procedure hasEl_i 58 end interface hasEl 59 public :: hasEl 60 61 interface name 62 module procedure name_ 63 end interface name 64 public :: name 65 66 interface copy 67 module procedure copy_ 68 end interface copy 69 public :: copy 70 71 public :: Eq_segs 72 73 public :: chem_pot_add_Elec 74 75 public :: fdf_nmu, fdffake_mu, fdf_mu 76 77 public :: print_mus_block 78 79 private 80 81contains 82 83 function fdf_nmu(prefix,kT,this_n) result(n) 84 use fdf 85 86 character(len=*), intent(in) :: prefix 87 real(dp), intent(in) :: kT 88 type(ts_mu), intent(inout), allocatable :: this_n(:) 89 integer :: n 90 91 ! prepare to read in the data... 92 type(block_fdf) :: bfdf 93 type(parsed_line), pointer :: pline => null() 94 integer :: i 95 96 real(dp) :: E_pole 97 logical :: found 98 99 n = 0 100 found = fdf_block(trim(prefix)//'.ChemPots',bfdf) 101 if ( .not. found ) return 102 103 ! first count the number of electrodes 104 n = 0 105 do while ( fdf_bline(bfdf,pline) ) 106 if ( fdf_bnnames(pline) == 0 ) cycle 107 n = n + 1 108 end do 109 110 allocate(this_n(n)) 111 this_n(:)%N_poles = fdf_get('TS.Contours.Eq.Pole.N',def_poles) 112 E_pole = fdf_get('TS.Contours.Eq.Pole',1.5_dp*eV,'Ry') 113 if ( E_pole > 0._dp ) then 114 call E2Npoles(E_pole,kT,i) 115 this_n(:)%N_poles = i 116 end if 117 118 ! rewind to read again 119 call fdf_brewind(bfdf) 120 121 n = 0 122 do while ( fdf_bline(bfdf,pline) ) 123 if ( fdf_bnnames(pline) == 0 ) cycle 124 n = n + 1 125 ! Attach the name 126 this_n(n)%Name = trim(fdf_bnames(pline,1)) 127 ! Save the ID of the chemical potential 128 this_n(n)%ID = n 129 if ( n > 1 ) then 130 ! Check that no name is the same 131 do i = 1 , n - 1 132 if ( leqi(name(this_n(i)),name(this_n(n))) ) then 133 call die('Chemical potential names must not be the same') 134 end if 135 end do 136 end if 137 end do 138 139 end function fdf_nmu 140 141 function fdffake_mu(this_n,kT,Volt) result(n) 142 use fdf 143 144 type(ts_mu), intent(inout), allocatable :: this_n(:) 145 ! SIESTA temperature 146 real(dp), intent(in) :: kT, Volt 147 real(dp) :: E_pole 148 integer :: n, n_pole 149 ! In case the user whishes to utilise the standard 150 ! transiesta setup we fake the chemical potentials 151 n = 2 152 allocate(this_n(n)) 153 154 ! Read in number of poles 155 this_n(:)%N_poles = fdf_get('TS.Contours.Eq.Pole.N',def_poles) 156 E_pole = fdf_get('TS.Contours.Eq.Pole',1.5_dp*eV,'Ry') 157 ! If the energy is larger than zero, the user requests 158 ! number of poles 159 if ( E_pole > 0._dp ) then 160 call E2Npoles(E_pole,kT,n_pole) 161 this_n(:)%N_poles = n_pole 162 end if 163 164 ! Set the temperature 165 this_n(:)%kT = kT 166 this_n(:)%ckT = ' ' 167 168 ! We star-mark the defaultet contours... 169 ! this will let us construct them readily... 170 171 ! *** NOTE this is hard-coded together with the ts_contour_eq 172 173 ! Create the left chemical potential... 174 this_n(1)%name = 'Left' 175 this_n(1)%mu = Volt * 0.5_dp 176 this_n(1)%cmu = 'V/2' 177 allocate(this_n(1)%Eq_seg(3)) ! one fake for poles 178 this_n(1)%Eq_seg(1) = '*c-left' 179 this_n(1)%Eq_seg(2) = '*t-left' 180 this_n(1)%ID = 1 181 182 this_n(2)%name = 'Right' 183 this_n(2)%mu = - Volt * 0.5_dp 184 this_n(2)%cmu = '-V/2' 185 allocate(this_n(2)%Eq_seg(3)) ! one fake for poles 186 this_n(2)%Eq_seg(1) = '*c-right' 187 this_n(2)%Eq_seg(2) = '*t-right' 188 this_n(2)%ID = 2 189 190 end function fdffake_mu 191 192 function fdf_mu(prefix,this, kT, Volt) result(found) 193 use fdf 194 use m_ts_io_ctype, only: pline_E_parse 195 196 character(len=*), intent(in) :: prefix 197 type(ts_mu), intent(inout) :: this 198 ! The SIESTA Electronic temperature, NOT the chemical potential 199 real(dp), intent(in) :: kT 200 real(dp), intent(in) :: Volt 201 logical :: found 202 203 ! prepare to read in the data... 204 type(block_fdf) :: bfdf 205 type(parsed_line), pointer :: pline => null() 206 logical :: info(2), bool_pole(2) 207 real(dp) :: E_pole 208 character(len=256) :: ln 209 210 found = fdf_block(trim(prefix)//'.ChemPot.'//trim(Name(this)),bfdf) 211 if ( .not. found ) return 212 213 info(:) = .false. 214#ifdef TBTRANS 215 ! TBtrans does not need the equilbrium contour information 216 info(2) = .true. 217#ifdef TBT_PHONON 218 ! PHtrans does not need to specify the chemical potential 219 info(1) = .true. 220#endif 221#endif 222 223 bool_pole = .false. 224 225 ! Initialize the temperature for this chemical potential 226 this%kT = kT 227 this%ckT = ' ' 228 229 do while ( fdf_bline(bfdf,pline) ) 230 if ( fdf_bnnames(pline) == 0 ) cycle 231 232 ln = trim(fdf_bnames(pline,1)) 233 234 ! We select the input 235 if ( leqi(ln,'chemical-shift') .or. & 236 leqi(ln,'mu') ) then 237 if ( fdf_bnvalues(pline) < 1 .and. & 238 fdf_bnnames(pline) < 2 ) call die('Chemical-shift not supplied') 239 240 ! utilize the io_ctype E-parser to grap the chemical potential 241 ! we force the energies to be before the last assignment 242 call pline_E_parse(pline,1,ln, & 243 val = this%mu, V = Volt, kT = kT, before=3) 244 this%cmu = trim(ln) 245 246 info(1) = .true. 247 248 else if ( leqi(ln,'contour.eq') ) then 249 ! we automatically make room for one pole contour 250 call read_contour_names('Equilibrium',this%Eq_seg,fakes=1) 251 info(2) = .true. 252 253 else if ( leqi(ln,'temp') .or. leqi(ln,'kT') .or. & 254 leqi(ln,'Electronic.Temperature') .or. & 255 leqi(ln,'ElectronicTemperature') ) then 256 257 call pline_E_parse(pline,1,ln, & 258 val = this%kT, kT = kT, before=3) 259 this%ckT = trim(ln) 260 261 else if ( leqi(ln,'contour.eq.pole.n') ) then 262 263 if ( fdf_bnintegers(pline) < 1 ) & 264 call die('You have not specified a number for & 265 &number of poles.') 266 267 this%N_poles = fdf_bintegers(pline,1) 268 bool_pole(1) = .true. 269 270 else if ( leqi(ln,'contour.eq.pole') ) then 271 272 call pline_E_parse(pline,1,ln, & 273 val = E_pole, kT = kT, before=3) 274 275 bool_pole(2) = .true. 276 277 else 278 279 call die('Unrecognized option "'//trim(ln)//'" & 280 &for chemical potential: '//trim(this%name)) 281 282 end if 283 284 end do 285 286 ! Check whether the Eq.Pole has been set, 287 if ( bool_pole(2) ) then 288 ! Update the number of poles 289 call E2Npoles(E_pole, this%kT, this%N_poles) 290 else if ( .not. bool_pole(1) ) then 291 E_pole = fdf_get('TS.Contours.Eq.Pole',1.5_dp*eV,'Ry') 292 if ( E_pole > 0._dp ) then 293 call E2Npoles(E_pole,this%kT,this%N_poles) 294 end if 295 end if 296 297 if ( .not. info(2) ) then 298 ! The contour.eq hasn't been defined 299 ! we assume the user want to use the continued 300 ! fraction method 301 302 ! we allow this 303 allocate(this%Eq_seg(1)) 304 ! The equilibrium here 305 this%Eq_seg(1) = 'cont-frac-'//trim(this%name) 306 307 if ( bool_pole(2) ) then 308 ! Update the number of poles (only up to 70% as 309 ! they become sparse 310 this%N_poles = int(E_pole * 0.7_dp / Pi / this%kT) 311 else if ( .not. bool_pole(1) ) then 312 ! Another default is 60 * Pi * kT ~ 60 poles 313 ! scale for the sparse top 314 E_pole = Pi * 60._dp * this%kT * 0.7_dp 315 E_pole = fdf_get('TS.Contours.Eq.Pole',E_pole,'Ry') 316 if ( E_pole > 0._dp ) then 317 this%N_poles = int(E_pole / Pi / this%kT) 318 end if 319 end if 320 321#ifndef TRANSIESTA_OVERRIDE 322 if ( this%N_poles < 20 ) then 323 call die('The continued fraction method requires at least 20 poles.') 324 end if 325#endif 326 info(2) = .true. 327 end if 328 329 if ( .not. all(info) ) then 330 write(*,*)'You need to supply at least:' 331 write(*,*)' - chemical-shift' 332 write(*,*)' - contour.eq' 333 call die('You have not supplied all chemical potential information') 334 end if 335 336 337 if ( this%N_poles < 1 ) then 338 print '(a)','Electrode: '//trim(this%name) 339 print '(a,i0)',' Number of poles: ',this%N_poles 340 print '(a,f10.2,a)',' Temperature: ', this%kT/Kelvin,' K' 341 call die('Number of poles must be larger than or equal to 1') 342 end if 343 344 contains 345 346 subroutine read_contour_names(name,con,fakes) 347 character(len=*), intent(in) :: name 348 character(len=C_N_NAME_LEN), intent(inout), allocatable :: con(:) 349 integer, intent(in), optional :: fakes 350 integer :: i, empty 351 352 character(len=256) :: ln 353 354 if ( allocated(con) ) deallocate(con) 355 356 ! we need to read in the equilibrium contour 357 ! skip to "begin" 358 if ( .not. fdf_bline(bfdf,pline) ) & 359 call die("Chemical potential block ended prematurely.") 360 361 ! read in the begin ... end block 362 ln = fdf_bnames(pline,1) 363 if ( .not. leqi(ln,"begin") ) & 364 call die(trim(name)//" contour errorneously formatted. & 365 &First line *must* be begin!") 366 367 ! Count lines 368 i = 0 369 empty = 0 370 do 371 if ( .not. fdf_bline(bfdf,pline) ) & 372 call die("Chemical potential block ended prematurely.") 373 if ( fdf_bnnames(pline) < 1 ) then 374 empty = empty + 1 375 else 376 ln = fdf_bnames(pline,1) 377 if ( leqi(ln,'end') ) exit 378 i = i + 1 379 end if 380 end do 381 382 ! allocate names 383 if ( present(fakes) ) then 384 allocate(con(i+fakes)) 385 empty = empty - fakes 386 else 387 allocate(con(i)) 388 end if 389 con = ' ' 390 do i = 0 , size(con) + empty 391 if ( .not. fdf_bbackspace(bfdf) ) & 392 call die("Backspacing too much...") 393 end do 394 i = 0 395 do 396 if ( .not. fdf_bline(bfdf,pline) ) & 397 call die("Chemical potential block ended prematurely.") 398 if ( fdf_bnnames(pline) < 1 ) cycle 399 ln = fdf_bnames(pline,1) 400 if ( leqi(ln,'end') ) exit 401 i = i + 1 402 if ( len_trim(ln) > C_N_NAME_LEN ) then 403 call die('Contour name: '//trim(ln)//' is too long, please use a & 404 &shorter name.') 405 end if 406 con(i) = trim(ln) 407 end do 408 409 end subroutine read_contour_names 410 411 end function fdf_mu 412 413 elemental function Eq_segs(this) result(count) 414 type(ts_mu), intent(in) :: this 415 integer :: count 416 count = size(this%Eq_seg) 417 end function Eq_segs 418 419 subroutine chem_pot_add_Elec(this,iEl) 420 type(ts_mu), intent(inout) :: this 421 integer, intent(in) :: iEl 422 integer, pointer :: tmp(:), tmp2(:) 423 424 nullify(tmp) 425 allocate(tmp(this%N_El+1)) 426 427 if ( this%N_El == 0 ) then 428 this%N_El = 1 429 this%el => tmp 430 this%el(1) = iEl 431 else if ( all(this%el /= iEl) ) then 432 ! copy over 433 tmp(1:this%N_El) = this%el(:) 434 this%N_El = this%N_El + 1 435 tmp(this%N_El) = iEl 436 ! clean up 437 tmp2 => this%el 438 this%el => tmp 439 deallocate(tmp2) 440 end if 441 end subroutine chem_pot_add_Elec 442 443 elemental function hasEl_i(this,iEl) result(has) 444 type(ts_mu), intent(in) :: this 445 integer, intent(in) :: iEl 446 logical :: has 447 has = this%N_El > 0 448 if ( has ) has = any(this%el == iEl) 449 end function hasEl_i 450 451 elemental function hasCio(this,c_io) result(has) 452 use m_ts_io_ctype, only : ts_c_io 453 type(ts_mu), intent(in) :: this 454 type(ts_c_io), intent(in) :: c_io 455 logical :: has 456 integer :: i 457 do i = 1 , Eq_segs(this) 458 has = this%Eq_seg(i) .eq. c_io%name 459 if ( has ) return 460 end do 461 end function hasCio 462 463 elemental function hasCeq(this,c) result(has) 464 use m_ts_cctype, only : ts_cw 465 type(ts_mu), intent(in) :: this 466 type(ts_cw), intent(in) :: c 467 logical :: has 468 has = hasCio(this,c%c_io) 469 end function hasCeq 470 471 elemental function name_(this) result(name) 472 type(ts_mu), intent(in) :: this 473 character(len=NAME_MU_LEN) :: name 474 name = this%name 475 end function name_ 476 477 subroutine print_mus_block(prefix,N_mu,mus) 478 use parallel, only : IONode 479 character(len=*), intent(in) :: prefix 480 integer, intent(in) :: N_mu 481 type(ts_mu), intent(in) :: mus(N_mu) 482 integer :: i 483 484 if (IONode) then 485 write(*,'(2a)') '%block ',trim(prefix)//'.ChemPots' 486 do i = 1 , N_mu 487 write(*,'(t3,a)') trim(mus(i)%name) 488 end do 489 write(*,'(2a,/)') '%endblock ',trim(prefix)//'.ChemPots' 490 end if 491 492 do i = 1 , N_mu 493 call print_mu_block(prefix,mus(i)) 494 end do 495 496 end subroutine print_mus_block 497 498 subroutine print_mu_block(prefix,this) 499 use parallel, only : IONode 500 use fdf 501 character(len=*), intent(in) :: prefix 502 type(ts_mu), intent(in) :: this 503 real(dp) :: E_pole 504 integer :: i, def_pole 505 character(len=64) :: ln 506 507 if ( .not. IONode ) return 508 509 def_pole = fdf_get('TS.Contours.Eq.Pole.N',def_poles) 510 E_pole = fdf_get('TS.Contours.Eq.Pole',1.5_dp*eV,'Ry') 511 if ( E_pole > 0._dp ) then 512 call E2Npoles(E_pole,this%kT, def_pole) 513 end if 514 515 ! Start by writing out the block beginning 516 write(*,'(a,a)') '%block '//trim(prefix)//'.ChemPot.',trim(this%name) 517 write(*,'(t3,a,tr2,a)') 'mu',trim(this%cmu) 518 write(*,'(t3,a)') 'contour.eq' 519 write(*,'(t4,a)') 'begin' 520 do i = 1 , Eq_segs(this) - 1 ! no pole-allowed 521 ln = this%Eq_seg(i) 522 if ( ln(1:1) == '*' ) then 523 write(*,'(t5,a)') trim(ln(2:)) 524 else 525 write(*,'(t5,a)') trim(ln) 526 end if 527 end do 528 write(*,'(t4,a)') 'end' 529 if ( def_pole /= this%N_poles ) then 530 write(*,'(t3,a,tr2,i0)') 'contour.eq.pole.n',this%N_poles 531 end if 532 if ( len_trim(this%ckT) > 0 ) then 533 write(*,'(t3,2a)') 'Temp ',this%ckT 534 end if 535 write(*,'(a,a)') '%endblock '//trim(prefix)//'.ChemPot.',trim(this%name) 536 537 end subroutine print_mu_block 538 539 subroutine E2Npoles(E_pole,kT,n_pole) 540 real(dp), intent(in) :: kT, E_pole 541 integer, intent(out) :: n_pole 542 real(dp) :: tmp 543 544 ! Calculate number of \pi kT 545 tmp = E_pole / ( Pi * kT ) 546 ! Using ceiling should also force it to at least 547 ! 1 pole. 548 n_pole = ceiling( tmp / 2._dp ) 549 550 end subroutine E2Npoles 551 552 subroutine copy_(this,copy) 553 type(ts_mu), intent(inout) :: this, copy 554 copy%name = this%name 555 copy%ID = this%ID 556 copy%N_poles = this%N_poles 557 copy%mu = this%mu 558 copy%cmu = this%cmu 559 copy%kT = this%kT 560 copy%ckT = this%ckT 561 copy%N_El = this%N_El 562 allocate(copy%Eq_seg(size(this%Eq_seg))) 563 copy%Eq_seg(:) = this%Eq_seg(:) 564 end subroutine copy_ 565 566end module m_ts_chem_pot 567 568