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