1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief deal with the Fermi distribution, compute it, fix mu, get derivs 8!> \author Joost VandeVondele 9!> \date 09.2008 10! ************************************************************************************************** 11MODULE fermi_utils 12 13 USE kahan_sum, ONLY: accurate_sum 14 USE kinds, ONLY: dp 15#include "base/base_uses.f90" 16 17 IMPLICIT NONE 18 19 PRIVATE 20 21 PUBLIC :: Fermi, FermiFixed, FermiFixedDeriv, Fermikp, Fermikp2 22 23 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fermi_utils' 24 INTEGER, PARAMETER, PRIVATE :: BISECT_MAX_ITER = 400 25 26CONTAINS 27! ************************************************************************************************** 28!> \brief returns occupations according to Fermi-Dirac statistics 29!> for a given set of energies and fermi level. 30!> Note that singly occupied orbitals are assumed 31!> \param f occupations 32!> \param N total number of electrons (output) 33!> \param kTS ... 34!> \param e eigenvalues 35!> \param mu Fermi level (input) 36!> \param T electronic temperature 37!> \param maxocc maximum occupation of an orbital 38!> \param estate excited state in core level spectroscopy 39!> \param festate occupation of the excited state in core level spectroscopy 40!> \date 09.2008 41!> \par History 42!> - Made estate and festate optional (LT, 2014/02/26) 43!> \author Joost VandeVondele 44! ************************************************************************************************** 45 SUBROUTINE Fermi(f, N, kTS, e, mu, T, maxocc, estate, festate) 46 47 REAL(KIND=dp), INTENT(out) :: f(:), N, kTS 48 REAL(KIND=dp), INTENT(IN) :: e(:), mu, T, maxocc 49 INTEGER, INTENT(IN), OPTIONAL :: estate 50 REAL(KIND=dp), INTENT(IN), OPTIONAL :: festate 51 52 CHARACTER(len=*), PARAMETER :: routineN = 'Fermi', routineP = moduleN//':'//routineN 53 54 INTEGER :: I, Nstate 55 REAL(KIND=dp) :: arg, occupation, term1, term2, tmp, & 56 tmp2, tmp3, tmp4, tmplog 57 58 Nstate = SIZE(e) 59 kTS = 0.0_dp 60 ! kTS is the entropic contribution to the energy i.e. -TS 61 ! kTS= kT*[f ln f + (1-f) ln (1-f)] 62 63 DO i = 1, Nstate 64 IF (PRESENT(estate) .AND. PRESENT(festate)) THEN 65 IF (i == estate) THEN 66 occupation = festate 67 ELSE 68 occupation = maxocc 69 END IF 70 ELSE 71 occupation = maxocc 72 END IF 73 ! have the result of exp go to zero instead of overflowing 74 IF (e(i) > mu) THEN 75 arg = -(e(i) - mu)/T 76 ! tmp is smaller than 1 77 tmp = EXP(arg) 78 tmp4 = tmp + 1.0_dp 79 tmp2 = tmp/tmp4 80 tmp3 = 1.0_dp/tmp4 81 ! log(1+eps), might need to be written more accurately 82 tmplog = -LOG(tmp4) 83 term1 = tmp2*(arg + tmplog) 84 term2 = tmp3*tmplog 85 ELSE 86 arg = (e(i) - mu)/T 87 ! tmp is smaller than 1 88 tmp = EXP(arg) 89 tmp4 = tmp + 1.0_dp 90 tmp2 = 1.0_dp/tmp4 91 tmp3 = tmp/tmp4 92 tmplog = -LOG(tmp4) 93 term1 = tmp2*tmplog 94 term2 = tmp3*(arg + tmplog) 95 END IF 96 97 f(i) = occupation*tmp2 98 kTS = kTS + T*occupation*(term1 + term2) 99 END DO 100 101 N = accurate_sum(f) 102 103 END SUBROUTINE Fermi 104 105! ************************************************************************************************** 106!> \brief Fermi function for KP cases 107!> \param f Occupation numbers 108!> \param nel Number of electrons (total) 109!> \param kTS Entropic energy contribution 110!> \param e orbital (band) energies 111!> \param mu chemical potential 112!> \param wk kpoint weights 113!> \param t Temperature 114!> \param maxocc Maximum occupation 115! ************************************************************************************************** 116 SUBROUTINE Fermi2(f, nel, kTS, e, mu, wk, t, maxocc) 117 118 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: f 119 REAL(KIND=dp), INTENT(OUT) :: nel, kTS 120 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: e 121 REAL(KIND=dp), INTENT(IN) :: mu 122 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wk 123 REAL(KIND=dp), INTENT(IN) :: t, maxocc 124 125 CHARACTER(len=*), PARAMETER :: routineN = 'Fermi2', routineP = moduleN//':'//routineN 126 127 INTEGER :: ik, is, nkp, nmo 128 REAL(KIND=dp) :: arg, beta, term1, term2, tmp, tmp2, & 129 tmp3, tmp4, tmplog 130 131 nmo = SIZE(e, 1) 132 nkp = SIZE(e, 2) 133 kTS = 0.0_dp 134 ! kTS is the entropic contribution to the energy i.e. -TS 135 ! kTS= kT*[f ln f + (1-f) ln (1-f)] 136 IF (t > 1.0e-14_dp) THEN 137 beta = 1.0_dp/t 138 DO ik = 1, nkp 139 DO is = 1, nmo 140 IF (e(is, ik) > mu) THEN 141 arg = -(e(is, ik) - mu)*beta 142 tmp = EXP(arg) 143 tmp4 = tmp + 1.0_dp 144 tmp2 = tmp/tmp4 145 tmp3 = 1.0_dp/tmp4 146 tmplog = -LOG(tmp4) 147 term1 = tmp2*(arg + tmplog) 148 term2 = tmp3*tmplog 149 ELSE 150 arg = (e(is, ik) - mu)*beta 151 tmp = EXP(arg) 152 tmp4 = tmp + 1.0_dp 153 tmp2 = 1.0_dp/tmp4 154 tmp3 = tmp/tmp4 155 tmplog = -LOG(tmp4) 156 term1 = tmp2*tmplog 157 term2 = tmp3*(arg + tmplog) 158 END IF 159 160 f(is, ik) = maxocc*tmp2 161 kTS = kTS + t*maxocc*(term1 + term2)*wk(ik) 162 END DO 163 END DO 164 ELSE 165 DO ik = 1, nkp 166 DO is = 1, nmo 167 IF (e(is, ik) <= mu) THEN 168 f(is, ik) = maxocc 169 ELSE 170 f(is, ik) = 0.0_dp 171 END IF 172 END DO 173 END DO 174 END IF 175 176 nel = 0.0_dp 177 DO ik = 1, nkp 178 nel = nel + accurate_sum(f(1:nmo, ik))*wk(ik) 179 END DO 180 181 END SUBROUTINE Fermi2 182 183! ************************************************************************************************** 184!> \brief returns occupations according to Fermi-Dirac statistics 185!> for a given set of energies and number of electrons. 186!> Note that singly occupied orbitals are assumed. 187!> could fail if the fermi level lies out of the range of eigenvalues 188!> (to be fixed) 189!> \param f occupations 190!> \param mu Fermi level (output) 191!> \param kTS ... 192!> \param e eigenvalues 193!> \param N total number of electrons (input) 194!> \param T electronic temperature 195!> \param maxocc maximum occupation of an orbital 196!> \param estate excited state in core level spectroscopy 197!> \param festate occupation of the excited state in core level spectroscopy 198!> \date 09.2008 199!> \par History 200!> - Made estate and festate optional (LT, 2014/02/26) 201!> \author Joost VandeVondele 202! ************************************************************************************************** 203 SUBROUTINE FermiFixed(f, mu, kTS, e, N, T, maxocc, estate, festate) 204 REAL(KIND=dp), INTENT(OUT) :: f(:), mu, kTS 205 REAL(KIND=dp), INTENT(IN) :: e(:), N, T, maxocc 206 INTEGER, INTENT(IN), OPTIONAL :: estate 207 REAL(KIND=dp), INTENT(IN), OPTIONAL :: festate 208 209 CHARACTER(len=*), PARAMETER :: routineN = 'FermiFixed', routineP = moduleN//':'//routineN 210 211 INTEGER :: iter, my_estate 212 REAL(KIND=dp) :: mu_max, mu_min, mu_now, my_festate, & 213 N_max, N_min, N_now 214 215 IF (PRESENT(estate) .AND. PRESENT(festate)) THEN 216 my_estate = estate 217 my_festate = festate 218 ELSE 219 my_estate = NINT(maxocc) 220 my_festate = my_estate 221 END IF 222 223! bisection search to find N 224! first bracket 225 226 mu_min = MINVAL(e) 227 iter = 0 228 DO 229 iter = iter + 1 230 CALL Fermi(f, N_min, kTS, e, mu_min, T, maxocc, my_estate, my_festate) 231 IF (N_min > N .OR. iter > 20) THEN 232 mu_min = mu_min - T 233 ELSE 234 EXIT 235 ENDIF 236 ENDDO 237 238 mu_max = MAXVAL(e) 239 iter = 0 240 DO 241 iter = iter + 1 242 CALL Fermi(f, N_max, kTS, e, mu_max, T, maxocc, my_estate, my_festate) 243 IF (N_max < N .OR. iter > 20) THEN 244 mu_max = mu_max + T 245 ELSE 246 EXIT 247 ENDIF 248 ENDDO 249 250 ! now bisect 251 iter = 0 252 DO WHILE (mu_max - mu_min > EPSILON(mu)*MAX(1.0_dp, ABS(mu_max), ABS(mu_min))) 253 iter = iter + 1 254 mu_now = (mu_max + mu_min)/2.0_dp 255 CALL Fermi(f, N_now, kTS, e, mu_now, T, maxocc, my_estate, my_festate) 256 iter = iter + 1 257 258 IF (N_now <= N) THEN 259 mu_min = mu_now 260 ELSE 261 mu_max = mu_now 262 ENDIF 263 264 IF (iter > BISECT_MAX_ITER) THEN 265 CPWARN("Maximum number of iterations reached while finding the Fermi energy") 266 EXIT 267 ENDIF 268 ENDDO 269 270 mu = (mu_max + mu_min)/2.0_dp 271 CALL Fermi(f, N_now, kTS, e, mu, T, maxocc, my_estate, my_festate) 272 273 END SUBROUTINE FermiFixed 274 275! ************************************************************************************************** 276!> \brief Bisection search to find mu for a given nel (kpoint case) 277!> \param f Occupation numbers 278!> \param mu chemical potential 279!> \param kTS Entropic energy contribution 280!> \param e orbital (band) energies 281!> \param nel Number of electrons (total) 282!> \param wk kpoint weights 283!> \param t Temperature 284!> \param maxocc Maximum occupation 285! ************************************************************************************************** 286 SUBROUTINE Fermikp(f, mu, kTS, e, nel, wk, t, maxocc) 287 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: f 288 REAL(KIND=dp), INTENT(OUT) :: mu, kTS 289 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: e 290 REAL(KIND=dp), INTENT(IN) :: nel 291 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wk 292 REAL(KIND=dp), INTENT(IN) :: t, maxocc 293 294 CHARACTER(len=*), PARAMETER :: routineN = 'Fermikp', routineP = moduleN//':'//routineN 295 REAL(KIND=dp), PARAMETER :: epsocc = 1.0e-12_dp 296 297 INTEGER :: iter 298 REAL(KIND=dp) :: de, mu_max, mu_min, N_now 299 300 ! bisection search to find mu for a given nel 301 de = t*LOG((1.0_dp - epsocc)/epsocc) 302 de = MAX(de, 0.5_dp) 303 mu_min = MINVAL(e) - de 304 mu_max = MAXVAL(e) + de 305 iter = 0 306 DO WHILE (mu_max - mu_min > EPSILON(mu)*MAX(1.0_dp, ABS(mu_max), ABS(mu_min))) 307 iter = iter + 1 308 mu = (mu_max + mu_min)/2.0_dp 309 CALL Fermi2(f, N_now, kTS, e, mu, wk, t, maxocc) 310 311 IF (ABS(N_now - nel) < nel*epsocc) EXIT 312 313 IF (N_now <= nel) THEN 314 mu_min = mu 315 ELSE 316 mu_max = mu 317 ENDIF 318 319 IF (iter > BISECT_MAX_ITER) THEN 320 CPWARN("Maximum number of iterations reached while finding the Fermi energy") 321 EXIT 322 ENDIF 323 ENDDO 324 325 mu = (mu_max + mu_min)/2.0_dp 326 CALL Fermi2(f, N_now, kTS, e, mu, wk, t, maxocc) 327 328 END SUBROUTINE Fermikp 329 330! ************************************************************************************************** 331!> \brief Bisection search to find mu for a given nel (kpoint case) 332!> \param f Occupation numbers 333!> \param mu chemical potential 334!> \param kTS Entropic energy contribution 335!> \param e orbital (band) energies 336!> \param nel Number of electrons (total) 337!> \param wk kpoint weights 338!> \param t Temperature 339! ************************************************************************************************** 340 SUBROUTINE Fermikp2(f, mu, kTS, e, nel, wk, t) 341 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: f 342 REAL(KIND=dp), INTENT(OUT) :: mu, kTS 343 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: e 344 REAL(KIND=dp), INTENT(IN) :: nel 345 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wk 346 REAL(KIND=dp), INTENT(IN) :: t 347 348 CHARACTER(len=*), PARAMETER :: routineN = 'Fermikp2', routineP = moduleN//':'//routineN 349 REAL(KIND=dp), PARAMETER :: epsocc = 1.0e-12_dp 350 351 INTEGER :: iter 352 REAL(KIND=dp) :: de, kTSa, kTSb, mu_max, mu_min, N_now, & 353 na, nb 354 355 ! only do spin polarized case 356 CPASSERT(SIZE(f, 3) == 2 .AND. SIZE(e, 3) == 2) 357 358 ! bisection search to find mu for a given nel 359 de = t*LOG((1.0_dp - epsocc)/epsocc) 360 de = MAX(de, 0.5_dp) 361 mu_min = MINVAL(e) - de 362 mu_max = MAXVAL(e) + de 363 iter = 0 364 DO WHILE (mu_max - mu_min > EPSILON(mu)*MAX(1.0_dp, ABS(mu_max), ABS(mu_min))) 365 iter = iter + 1 366 mu = (mu_max + mu_min)/2.0_dp 367 CALL Fermi2(f(:, :, 1), na, kTSa, e(:, :, 1), mu, wk, t, 1.0_dp) 368 CALL Fermi2(f(:, :, 2), nb, kTSb, e(:, :, 2), mu, wk, t, 1.0_dp) 369 N_now = na + nb 370 371 IF (ABS(N_now - nel) < nel*epsocc) EXIT 372 373 IF (N_now <= nel) THEN 374 mu_min = mu 375 ELSE 376 mu_max = mu 377 ENDIF 378 379 IF (iter > BISECT_MAX_ITER) THEN 380 CPWARN("Maximum number of iterations reached while finding the Fermi energy") 381 EXIT 382 ENDIF 383 ENDDO 384 385 mu = (mu_max + mu_min)/2.0_dp 386 CALL Fermi2(f(:, :, 1), na, kTSa, e(:, :, 1), mu, wk, t, 1.0_dp) 387 CALL Fermi2(f(:, :, 2), nb, kTSb, e(:, :, 2), mu, wk, t, 1.0_dp) 388 kTS = kTSa + kTSb 389 390 END SUBROUTINE Fermikp2 391 392! ************************************************************************************************** 393!> \brief returns f and dfde for a given set of energies and number of electrons 394!> it is a numerical derivative, trying to use a reasonable step length 395!> it ought to yield an accuracy of approximately EPSILON()^(2/3) (~10^-11) 396!> l ~ 10*T yields best accuracy 397!> Note that singly occupied orbitals are assumed. 398!> To be fixed: this could be parallellized for better efficiency 399!> \param dfde derivatives of the occupation numbers with respect to the eigenvalues 400!> the ith column is the derivative of f wrt to e_i 401!> \param f occupations 402!> \param mu Fermi level (input) 403!> \param kTS ... 404!> \param e eigenvalues 405!> \param N total number of electrons (output) 406!> \param T electronic temperature 407!> \param maxocc ... 408!> \param l typical lenght scale (~ 10 * T) 409!> \param estate ... 410!> \param festate ... 411!> \date 09.2008 412!> \par History 413!> - Made estate and festate optional (LT, 2014/02/26) 414!> - Changed order of input, so l is before the two optional varaibles 415!> (LT, 2014/02/26) 416!> \author Joost VandeVondele 417! ************************************************************************************************** 418 SUBROUTINE FermiFixedDeriv(dfde, f, mu, kTS, e, N, T, maxocc, l, estate, festate) 419 REAL(KIND=dp), INTENT(OUT) :: dfde(:, :), f(:), mu, kTS 420 REAL(KIND=dp), INTENT(IN) :: e(:), N, T, maxocc, l 421 INTEGER, INTENT(IN), OPTIONAL :: estate 422 REAL(KIND=dp), INTENT(IN), OPTIONAL :: festate 423 424 CHARACTER(len=*), PARAMETER :: routineN = 'FermiFixedDeriv', & 425 routineP = moduleN//':'//routineN 426 427 INTEGER :: handle, I, my_estate, Nstate 428 REAL(KIND=dp) :: h, mux, my_festate 429 REAL(KIND=dp), ALLOCATABLE :: ex(:), fx(:) 430 431 CALL timeset(routineN, handle) 432 433 Nstate = SIZE(e) 434 ALLOCATE (ex(Nstate), fx(Nstate)) 435 436 IF (PRESENT(estate) .AND. PRESENT(festate)) THEN 437 my_estate = estate 438 my_festate = festate 439 ELSE 440 my_estate = NINT(maxocc) 441 my_festate = my_estate 442 END IF 443 444 DO I = 1, Nstate 445 ! NR 5.7.8 446 ! the problem here is that each f_i 'seems to have' a different lenght scale 447 ! and it would be to expensive to compute each single df_i/de_i using a finite difference 448 h = (EPSILON(h)**(1.0_dp/3.0_dp))*l 449 ! get an exact machine representable number close to this h 450 h = 2.0_dp**EXPONENT(h) 451 ! this should write three times the same number 452 ! write(*,*) h,(e(i)+h)-e(i),(e(i)-h)-e(i) 453 ! and the symmetric finite difference 454 ex(:) = e 455 ex(i) = e(i) + h 456 CALL FermiFixed(fx, mux, kTS, ex, N, T, maxocc, my_estate, my_festate) 457 dfde(:, I) = fx 458 ex(i) = e(i) - h 459 CALL FermiFixed(fx, mux, kTS, ex, N, T, maxocc, my_estate, my_festate) 460 dfde(:, I) = (dfde(:, I) - fx)/(2.0_dp*h) 461 ENDDO 462 DEALLOCATE (ex, fx) 463 464 CALL FermiFixed(f, mu, kTS, e, N, T, maxocc, my_estate, my_festate) 465 466 CALL timestop(handle) 467 468 END SUBROUTINE FermiFixedDeriv 469 470END MODULE fermi_utils 471