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 module m_fermid 9 10 use precision, only : dp 11 use parallel, only : IOnode 12 use fdf, only : fdf_boolean, fdf_integer, fdf_string, leqi 13 use m_errorf, only: derfc 14 use sys, only: die 15 16 implicit none 17 18 public :: fermid, fermispin, stepf 19 private :: enpy, whg, hp 20 21 private 22 23 CONTAINS 24 25 subroutine fermid(nspinor, maxspn, NK, WK, maxe, NE, E, 26 . temp, qtot, WKE, EF, entropy ) 27 28C ********************************************************************* 29C Finds the Fermi energy and the occupation weights of states. 30C Written by J.M.Soler. August'96. 31C Simple single excitation introduced by E. Artacho August 2002. 32C Alternative occupation functions introduced by P. Ordejon, May'03. 33C ********** INPUT **************************************************** 34C INTEGER nspinor : Number of spinors (1 or 2) 35! nspinor = 1 for spin-less calculations 36! nspinor = 2 for collinear or non-collinear spin 37C INTEGER maxspn : Number of separate spin blocks (1 or 2) 38! in the E and WKE matrices (that is, the length of 39! the second dimension of E and WKE). 40! For NC and SOC this should be 1, since the states 41! are not classified by spin. 42! For collinear spin, maxspn=2 43! For the spinless case, maxspn=1 44! 45! Note that in the NC/SOC case, the actual arguments E and WKE are of 46! the form E(2*neigwanted,nk) or E(2*no_u,nk) 47! and in this routine they are unpacked to the E(*,1,nk) form 48! 49C INTEGER NK : Number of K-points 50C REAL*8 WK(NK) : Sampling weights of k-points (must sum 1) 51C INTEGER maxe : First dimension of E and WKE 52C For NC/SOC this should be large enough to hold twice the number of 53C wanted eigenstates (2*neigwanted). 54C INTEGER NE : Number of bands (per spin in the collinear case) 55C REAL*8 E(maxe,maxspn,NK) : State eigenvalues 56C REAL*8 temp : Temperature (in the same units of E) 57C REAL*8 qtot : Total valence charge (number of electrons) 58C ********** OUTPUT *************************************************** 59C REAL*8 WKE(maxe,maxspn,NK) : Occupations multiplied by k-point weights 60C (sum qtot) 61C REAL*8 EF : Fermi energy 62C REAL*8 entropy : Entropy contribution to the electronic 63C Free Energy 64C ********************************************************************* 65 66C Passed variables 67 integer, intent(in) :: nspinor, maxspn 68 integer :: maxe 69 integer :: ne 70 integer, intent(in) :: nk 71 72 real(dp), intent(in) :: E(maxe,maxspn,nk) 73 real(dp) :: entropy 74 real(dp) :: qtot 75 real(dp) :: temp 76 real(dp) :: wke(maxe,maxspn,nk) 77 real(dp), intent(in) :: wk(nk) 78 real(dp) :: ef 79 80C Local variables 81 integer :: ie 82 integer :: ief 83 integer :: ik 84 integer :: ispin 85 integer :: iter 86 integer, save :: nitmax = 150 87 logical, save :: blread = .false. 88 logical, save :: excitd = .false. 89 real(dp), save :: tol = 1.0d-10 90 real(dp) :: sumq, emin, emax 91 real(dp) :: t, tinv, drange, wkebuf, W, eik 92 93#ifdef DEBUG 94 call write_debug( ' PRE fermid' ) 95#endif 96 97C Reading whether to do an excited state 98 if ( .not. blread) then 99 excitd = fdf_boolean('SingleExcitation', .false.) 100 if ( ionode ) then 101 if ( excitd ) write(6,'(/a)') 102 . 'fermid: Calculating for lowest-exciton excited state' 103 endif 104 blread = .true. 105 end if 106 107C Zero occupancies, including those not explicitly 108C calculated here if ne < maxe 109 wke(1:maxe,1:maxspn,1:nk) = 0.0d0 110 111C Determine Fermi level 112 sumq = 0.0d0 113 emin = e(1,1,1) 114 emax = e(1,1,1) 115 if ( nspinor == 1 ) then ! also maxspn == 1, but this is implicit 116 ! Each spinless state can have two electrons, so we double the occupancies 117 ! This will be reflected further down as well 118 do ik = 1,nk 119 sumq = sumq + wk(ik) * ne * 2._dp 120 do ie = 1,ne 121 wke(ie,1,ik) = wk(ik) * 2._dp 122 emin = min(emin,e(ie,1,ik)) 123 emax = max(emax,e(ie,1,ik)) 124 end do 125 end do 126 else 127 do ik = 1,nk 128 sumq = sumq + wk(ik) * maxspn * ne 129 do ispin = 1, maxspn 130 do ie = 1, ne 131 wke(ie,ispin,ik) = wk(ik) 132 emin = min(emin,e(ie,ispin,ik)) 133 emax = max(emax,e(ie,ispin,ik)) 134 end do 135 end do 136 end do 137 end if 138 139 ef = emax 140 141 if (abs(sumq-qtot).lt.tol) then 142 if (excitd) then 143 if (ionode) then 144 write (6,'(/a)') 145 . 'Fermid: Bands full, no excitation possible' 146 endif 147 call die() 148 else 149#ifdef DEBUG 150 call write_debug( ' POS fermid' ) 151#endif 152 return 153 endif 154 endif 155 if (sumq.lt.qtot) then 156 if (ionode) then 157 write(6,*) 'Fermid: Not enough states' 158 write(6,*) 'Fermid: qtot,sumq=',qtot,sumq 159 endif 160 call die() 161 endif 162 T = max(temp,1.d-6) 163 Tinv = 1._dp / T 164 drange = T*sqrt(-log(tol*0.01d0)) 165 emin = emin - drange 166 emax = emax + drange 167 do iter = 1,nitmax 168 ef = 0.5d0*(emin + emax) 169 170 sumq = 0.0d0 171 if ( nspinor == 1 ) then ! also maxspn == 1, but this is implicit 172 do ik = 1,nk 173 do ie = 1,ne 174 wke(ie,1,ik) = wk(ik)* 2._dp * 175 . stepf((e(ie,1,ik)-ef)*Tinv) 176 sumq = sumq + wke(ie,1,ik) 177 end do 178 end do 179 else 180 do ik = 1,nk 181 do ispin = 1, maxspn 182 do ie = 1,ne 183 wke(ie,ispin,ik) = wk(ik) * 184 . stepf((e(ie,ispin,ik)-ef)*Tinv) 185 sumq = sumq + wke(ie,ispin,ik) 186 end do 187 end do 188 end do 189 end if 190 191C If the Fermi level was found..................... 192 if (abs(sumq-qtot).lt.tol) then 193 194C If excited state is to be calculated, find the level above Ef for 195C k=1 and spin=1 196 if (excitd) then 197 loop: do ie = 1, ne 198 if ( e(ie,1,1) .gt. ef ) then 199 ief = ie ! LUMO index 200 exit loop 201 endif 202 enddo loop 203 204C and swap populations (meaningful only for T close to 0): 205C if nspinor=1 populations of homo and lumo are just swapped for is=1 206C if nspinor=2 populations of homo and lumo become equal. 207 if ( nspinor == 1 ) then 208 wkebuf = ( wke(ief-1,1,1) - wke(ief,1,1) ) * 0.5_dp 209 else 210 wkebuf = wke(ief-1,1,1) - wke(ief,1,1) 211 end if 212 wke(ief,1,1) = wke(ief,1,1) + wkebuf 213 wke(ief-1,1,1) = wke(ief-1,1,1) - wkebuf 214 endif 215 216C Obtain the electronic entropy 217 entropy = 0.0_dp 218 if ( nspinor == 1 ) then ! also maxspn == 1, but this is implicit 219 do ik = 1 , nk 220 do ie = 1 , ne 221 W = 0.5_dp * wke(ie,1,ik) / wk(ik) 222 eik = (e(ie,1,ik)-ef) * Tinv 223 entropy = entropy + 2.0d0 * wk(ik) * enpy(eik,W) 224 end do 225 end do 226 else 227 do ik = 1 , nk 228 do ispin = 1 , maxspn 229 do ie = 1 , ne 230 W = wke(ie,ispin,ik) / wk(ik) 231 eik = (e(ie,ispin,ik)-ef) * Tinv 232 entropy = entropy + wk(ik) * enpy(eik,W) 233 end do 234 end do 235 end do 236 end if 237 238#ifdef DEBUG 239 call write_debug( ' POS fermid' ) 240#endif 241 return 242 243 endif 244 245 if (sumq.le.qtot) emin = ef 246 if (sumq.ge.qtot) emax = ef 247 enddo 248 249 if (ionode) then 250 write(6,*) 'Fermid: Iteration has not converged.' 251 write(6,*) 'Fermid: qtot,sumq=',qtot,sumq 252 endif 253 call die('Fermid: Iteration has not converged.') 254 255 end subroutine fermid 256 257!--------------------------------------------------------- 258 259 260 subroutine fermispin( nspin, maxspn, NK, WK, maxe, NE, E, 261 . temp, qtot, WKE, EF, entropy ) 262 263C ********************************************************************* 264C Finds the Fermi energy and the occupation weights of states, 265C for the case where the total spin of the calculation is fixed 266C to a given value. 267C Written by J.M.Soler. August'96. 268C Version modified for fixed spin configurations: P. Ordejon'03-04 269C ********** INPUT **************************************************** 270C INTEGER nspin : Number of different spin polarizations (1 or 2) 271C INTEGER maxspn : Maximum number of different spin polarizations (1 or 2) 272C for E and WKE matrices dimensions 273C INTEGER NK : Number of K-points 274C REAL*8 WK(NK) : Sampling weights of k-points (must sum 1) 275C INTEGER maxe : First dimension of E and WKE 276C INTEGER NE : Number of bands 277C REAL*8 E(maxe,maxspn,NK) : State eigenvalues 278C REAL*8 temp : Temperature (in the same units of E) 279C REAL*8 qtot(maxspn) : Total valence charge (number of electrons) 280C for each spin component 281C ********** OUTPUT *************************************************** 282C REAL*8 WKE(maxe,maxspn,NK) : Occupations multiplied by k-point weights 283C (sum qtot) 284C REAL*8 EF(nspin) : Fermi energy (for each spin, if qtot 285C is different for each spin component. 286C REAL*8 entropy : Entropy contribution to the electronic 287C Free Energy 288C ********************************************************************* 289 290 integer :: maxe, maxspn, nk, nspin, nitmax, ispin 291 integer :: ik, ie, ne, iter 292 real(dp) :: entropy 293 294 real(dp) :: E(maxe,maxspn,NK),emin(4),emax(4), 295 . EF(nspin),qtot(maxspn),sumq(4),temp, 296 . WKE(maxe,maxspn,NK),WK(NK) 297 real(dp) :: t, drange, w, eik, tol 298 logical :: conv 299 parameter (tol=1.0D-10,nitmax=150) 300 301 conv = .false. 302 303 do ispin = 1,nspin 304 sumq(ispin) = 0.0d0 305 enddo 306 do ispin = 1,nspin 307 emin(ispin) = E(1,ispin,1) 308 emax(ispin) = E(1,ispin,1) 309 enddo 310 311C Zero occupancies, including those not explicitly 312C calculated here if ne < maxe 313 wke(1:maxe,1:nspin,1:nk) = 0.0d0 314 315 do ik = 1,nk 316 do ispin = 1,nspin 317 do ie = 1,ne 318 wke(ie,ispin,IK) = wk(ik)*2.0d0/dble(nspin) 319 sumq(ispin) = sumq(ispin) + wke(ie,ispin,ik) 320 emin(ispin) = min(emin(ispin),e(ie,ispin,ik)) 321 emax(ispin) = max(emax(ispin),e(ie,ispin,ik)) 322 enddo 323 enddo 324 enddo 325 do ispin = 1,nspin 326 ef(ispin) = emax(ispin) 327 enddo 328 conv = .true. 329 do ispin = 1,nspin 330 if (abs(sumq(ispin)-qtot(ispin)).gt.tol) conv = .false. 331 enddo 332 333 if (conv) goto 100 334 335 do ispin = 1,nspin 336 if (sumq(ispin).lt.qtot(ispin)) then 337 if (ionode) then 338 write(6,*) 'Fermid: Not enough states' 339 write(6,*) 'Fermid: ispin,qtot,sumq=', 340 . ispin,qtot(ispin),sumq(ispin) 341 endif 342 call die() 343 ENDIF 344 enddo 345 T = max(temp,1.0d-6) 346 drange = T*sqrt(-log(tol*0.01d0)) 347 do ispin = 1,nspin 348 emin(ispin) = emin(ispin) - drange 349 emax(ispin) = emax(ispin) + drange 350 enddo 351 do iter = 1,nitmax 352 do ispin = 1,nspin 353 ef(ispin) = 0.5d0*(emin(ispin)+emax(ispin)) 354 sumq(ispin) = 0.0d0 355 enddo 356 do ik = 1,nk 357 do ispin = 1,nspin 358 do ie = 1,ne 359 wke(ie,ispin,ik) = wk(ik)* 360 . stepf((E(ie,ispin,ik)-ef(ispin))/T)*2.0d0/dble(nspin) 361 sumq(ispin)=sumq(ispin)+wke(ie,ispin,ik) 362 enddo 363 enddo 364 enddo 365 conv = .true. 366 do ispin = 1,nspin 367 if (abs(sumq(ispin)-qtot(ispin)).gt.tol) conv = .false. 368 enddo 369 if (conv) goto 100 370 do ispin = 1,nspin 371 if (sumq(ispin).le.qtot(ispin)) emin(ispin) = ef(ispin) 372 if (sumq(ispin).ge.qtot(ispin)) emax(ispin) = ef(ispin) 373 enddo 374 enddo 375 376 if (ionode) then 377 write (6,*) 'Fermid: Iteration has not converged.' 378 do ispin = 1,nspin 379 write (6,*) 'Fermid: ispin,qtot,sumq=', 380 . ispin,qtot(ispin),sumq(ispin) 381 enddo 382 endif 383 call die('Fermid: Iteration has not converged.') 384 385100 continue 386 entropy = 0.0d0 387 388 do ik = 1,nk 389 do ie = 1,ne 390 do ispin = 1,nspin 391 392 W = (nspin / 2.0d0) * wke(ie,ispin,ik) / wk(ik) 393 eik = (E(ie,ispin,ik)-ef(ispin)) / T 394 395 entropy = entropy + ( 2.0d0 * wk(ik) / nspin ) * 396 . enpy(eik,W) 397 398 enddo 399 enddo 400 enddo 401 402 return 403 404 end subroutine fermispin 405 406 407 real(dp) function stepf(X) 408 real(dp), intent(in) :: x 409 410 integer :: i, j 411 real(dp) :: pi, a, gauss 412 413C Local variables 414 character(len=22), save :: ocf = 'FD' 415 integer, save :: nh = 1 416 integer, save :: ocupfnct 417 logical, save :: ocfread = .false. 418 419 parameter (PI = 3.14159265358979D0) 420 421C Reading which electronic occupation function to use ------- 422 423 if ( .not. ocfread) then 424 425 ocf = fdf_string('OccupationFunction','FD') 426 427 if (leqi(ocf,'FD')) then 428 ocupfnct = 1 429 if ( ionode ) then 430 write(6,'(/a)') 'stepf: Fermi-Dirac step function' 431 endif 432 else if (leqi(ocf,'MP')) then 433 ocupfnct = 2 434 nh = fdf_integer('OccupationMPOrder',1) 435 if ( ionode ) then 436 write(6,'(/a)') 'stepf: Methfessel-Paxton step function' 437 write(6,'(a,i2)') 438 . ' Using Hermite-Gauss polynomials of order ',nh 439 endif 440 else 441 call die('fermid: Error: Allowed values ' 442 $ // 'for OccupationFunction are FD and MP') 443 endif 444 ocfread = .true. 445 endif 446c-------------------------------------- 447 448C Complementary error function. Ref: Fu & Ho, PRB 28, 5480 (1983) 449* STEPF=DERFC(X) - not available 450 451 if (ocupfnct .eq. 1) then 452 453C Fermi-Dirac distribution 454 if (x.gt.100.D0) then 455 stepf = 0.D0 456 elseif (x.lt.-100.d0) then 457 stepf = 1.d0 458 else 459 stepf = 1.d0 / ( 1.d0 + exp(x) ) 460 endif 461 462 463 else if (ocupfnct .eq. 2) then 464 465C Improved step function. Ref: Methfessel & Paxton PRB40 (15/Aug/89) 466C NH is the order of the Hemite polynomial expansion. 467 468 stepf = 0.5d0 * derfc(x) 469 a = 1.0d0/sqrt(pi) 470 471 do i = 1,nh 472 473C Get coefficients in Hermite-Gauss expansion 474 A = -A / (I * 4.0d0) 475 476C Get contribution to step function at order I 477 gauss = dexp(-x*x) 478 J = 2*I -1 479 if (gauss .gt. 1.d-20) stepf = stepf + a * hp(x,j) * gauss 480 481 enddO 482 483 else 484 485 call die( 'Stepf: Incorrect step function') 486 487 endif 488 489 end function stepf 490 491 492 real(dp) function enpy(E,W) 493 real(dp), intent(in) :: e, w 494 495C Computes the contribution of a given state with energy E 496C (refered to the Fermi energy, in units of the smearing 497C temperature) and occupation W to the electronic entropy 498C P. Ordejon, June 2003 499 500C Local variables 501 character(len=22), save :: ocf = 'FD' 502 integer, save :: nh = 1 503 integer, save :: ocupfnct 504 logical, save :: ocfread = .false. 505 real(dp) :: tiny, wo, we 506 507 data tiny /1.d-15/ 508 509c Reading which electronic occupation function to use ------- 510 511 if ( .not. ocfread) then 512 513 ocf = fdf_string('OccupationFunction', 'FD') 514 515 if (leqi(ocf,'FD')) then 516 ocupfnct = 1 517 else if (leqi(ocf,'MP')) then 518 ocupfnct = 2 519 nh = fdf_integer('OccupationMPOrder',1) 520 else 521 call die('fermid: Error: Allowed values ' 522 $ // 'for OccupationFunction are FD and MP') 523 endif 524 ocfread = .true. 525 endif 526c-------------------------------------- 527 528 529 if (ocupfnct .eq. 1) then 530 531C Mermin entropy for the Fermi-Dirac distribution 532 wo = max( w, tiny ) 533 we = 1.0d0 - wo 534 we = max( we, tiny ) 535 536 enpy = - 1.0 * ( wo*log(wo) + we*log(we) ) 537 . 538 539 else if (ocupfnct .eq. 2) then 540 541C Entropy for the Improved step function. 542C Ref: Methfessel & Paxton PRB40 (15/Aug/89) 543 544 enpy = whg(e,nh) 545 546 else 547 548 call die( 'Stepf: Incorrect step function') 549 550 endif 551 552 END function enpy 553 554 555 real(dp) function whg(x,n) 556C 557C Computes the factors to get the entropy term 558C for the Methfessel-Paxton smearing with Hermite 559C polynomials of order N 560C 561C P. Ordejon, June '03 562C 563 564C Passed variables 565 integer :: n 566 real(dp) :: x 567 568C Local variables 569 integer :: i 570 real(dp) :: a 571 real(dp) :: gauss 572 real(dp), save :: pi = 3.14159265358979d0 573 real(dp) :: x2 574 575 x2 = x**2.0d0 576 577C Get coefficients 578 579 a = 1.0d0/sqrt(pi) 580 do i = 1,n 581 a = - a / (dble(I) * 4.0d0) 582 enddo 583 584 gauss = dexp(-x2) 585 whg = 0.0d0 586 if (gauss .gt. 1.0d-20) whg = 0.5d0*a*hp(x,2*n)*gauss 587 588 return 589 end function whg 590 591 592 real(dp) function hp(x,n) 593C 594C Returns the value of the Hermite polynomial of degree N 595C evaluated at X. 596C 597C H_0 (x) = 1 598C H_1 (x) = 2x 599C ... 600C H_n+1(x) = 2 x H_n(x) - 2 n H_n-1(x) 601C 602C P. Ordejon, June 2003 603C 604 605C Passed variables 606 integer :: n 607 real(dp) :: x 608 609C Local variables 610 integer :: i 611 real(dp) :: hm1 612 real(dp) :: hm2 613 614 if (n .gt. 1000) 615 . call die('Fermid: Order of Hermite polynomial too large') 616 617 618 hp = 1.0d0 619 hm2 = 0.0d0 620 hm1 = 1.0d0 621 622 do i = 1,n 623 624 hp = 2.0_dp * (x * hm1 - dble(i-1) * hm2) 625 hm2 = hm1 626 hm1 = hp 627 628 enddo 629 630 end function hp 631 632 end module m_fermid 633