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