fermi_utils.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief deal with the Fermi distribution, compute it, fix mu, get derivs
!> \author Joost VandeVondele
!> \date 09.2008
! **************************************************************************************************
MODULE fermi_utils
 
   USE kahan_sum,                       ONLY: accurate_sum
   USE kinds,                           ONLY: dp
#include "base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   PUBLIC :: fermi, fermifixed, fermifixedderiv, fermikp, fermikp2
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fermi_utils'
   INTEGER, PARAMETER, PRIVATE          :: BISECT_MAX_ITER = 400
 
CONTAINS
! **************************************************************************************************
!> \brief   returns occupations according to Fermi-Dirac statistics
!>          for a given set of energies and fermi level.
!>          Note that singly occupied orbitals are assumed
!> \param   f occupations
!> \param   N total number of electrons (output)
!> \param kTS ...
!> \param   e eigenvalues
!> \param   mu Fermi level (input)
!> \param   T  electronic temperature
!> \param   maxocc maximum occupation of an orbital
!> \param   estate excited state in core level spectroscopy
!> \param   festate occupation of the excited state in core level spectroscopy
!> \date    09.2008
!> \par History
!>          - Made estate and festate optional (LT, 2014/02/26)
!> \author  Joost VandeVondele
! **************************************************************************************************
   SUBROUTINE fermi(f, N, kTS, e, mu, T, maxocc, estate, festate)
 
      REAL(kind=dp), INTENT(out)                         :: f(:), n, kts
      REAL(kind=dp), INTENT(IN)                          :: e(:), mu, t, maxocc
      INTEGER, INTENT(IN), OPTIONAL                      :: estate
      REAL(kind=dp), INTENT(IN), OPTIONAL                :: festate
 
      INTEGER                                            :: i, nstate
      REAL(kind=dp)                                      :: arg, occupation, term1, term2, tmp, &
                                                            tmp2, tmp3, tmp4, tmplog
 
      nstate = SIZE(e)
      kts = 0.0_dp
      ! kTS is the entropic contribution to the energy i.e. -TS
      ! kTS= kT*[f ln f + (1-f) ln (1-f)]
 
      DO i = 1, nstate
         IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
            IF (i == estate) THEN
               occupation = festate
            ELSE
               occupation = maxocc
            END IF
         ELSE
            occupation = maxocc
         END IF
         ! have the result of exp go to zero instead of overflowing
         IF (e(i) > mu) THEN
            arg = -(e(i) - mu)/t
            ! tmp is smaller than 1
            tmp = exp(arg)
            tmp4 = tmp + 1.0_dp
            tmp2 = tmp/tmp4
            tmp3 = 1.0_dp/tmp4
            ! log(1+eps), might need to be written more accurately
            tmplog = -log(tmp4)
            term1 = tmp2*(arg + tmplog)
            term2 = tmp3*tmplog
         ELSE
            arg = (e(i) - mu)/t
            ! tmp is smaller than 1
            tmp = exp(arg)
            tmp4 = tmp + 1.0_dp
            tmp2 = 1.0_dp/tmp4
            tmp3 = tmp/tmp4
            tmplog = -log(tmp4)
            term1 = tmp2*tmplog
            term2 = tmp3*(arg + tmplog)
         END IF
 
         f(i) = occupation*tmp2
         kts = kts + t*occupation*(term1 + term2)
      END DO
 
      n = accurate_sum(f)
 
   END SUBROUTINE fermi
 
! **************************************************************************************************
!> \brief Fermi function for KP cases
!> \param f   Occupation numbers
!> \param nel Number of electrons (total)
!> \param kTS Entropic energy contribution
!> \param e   orbital (band) energies
!> \param mu  chemical potential
!> \param wk  kpoint weights
!> \param t   Temperature
!> \param maxocc Maximum occupation
! **************************************************************************************************
   SUBROUTINE fermi2(f, nel, kTS, e, mu, wk, t, maxocc)
 
      REAL(kind=dp), DIMENSION(:, :), INTENT(OUT)        :: f
      REAL(kind=dp), INTENT(OUT)                         :: nel, kts
      REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: e
      REAL(kind=dp), INTENT(IN)                          :: mu
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: wk
      REAL(kind=dp), INTENT(IN)                          :: t, maxocc
 
      INTEGER                                            :: ik, is, nkp, nmo
      REAL(kind=dp)                                      :: arg, beta, term1, term2, tmp, tmp2, &
                                                            tmp3, tmp4, tmplog
 
      nmo = SIZE(e, 1)
      nkp = SIZE(e, 2)
      kts = 0.0_dp
      ! kTS is the entropic contribution to the energy i.e. -TS
      ! kTS= kT*[f ln f + (1-f) ln (1-f)]
      IF (t > 1.0e-14_dp) THEN
         beta = 1.0_dp/t
         DO ik = 1, nkp
            DO is = 1, nmo
               IF (e(is, ik) > mu) THEN
                  arg = -(e(is, ik) - mu)*beta
                  tmp = exp(arg)
                  tmp4 = tmp + 1.0_dp
                  tmp2 = tmp/tmp4
                  tmp3 = 1.0_dp/tmp4
                  tmplog = -log(tmp4)
                  term1 = tmp2*(arg + tmplog)
                  term2 = tmp3*tmplog
               ELSE
                  arg = (e(is, ik) - mu)*beta
                  tmp = exp(arg)
                  tmp4 = tmp + 1.0_dp
                  tmp2 = 1.0_dp/tmp4
                  tmp3 = tmp/tmp4
                  tmplog = -log(tmp4)
                  term1 = tmp2*tmplog
                  term2 = tmp3*(arg + tmplog)
               END IF
 
               f(is, ik) = maxocc*tmp2
               kts = kts + t*maxocc*(term1 + term2)*wk(ik)
            END DO
         END DO
      ELSE
         DO ik = 1, nkp
            DO is = 1, nmo
               IF (e(is, ik) <= mu) THEN
                  f(is, ik) = maxocc
               ELSE
                  f(is, ik) = 0.0_dp
               END IF
            END DO
         END DO
      END IF
 
      nel = 0.0_dp
      DO ik = 1, nkp
         nel = nel + accurate_sum(f(1:nmo, ik))*wk(ik)
      END DO
 
   END SUBROUTINE fermi2
 
! **************************************************************************************************
!> \brief   returns occupations according to Fermi-Dirac statistics
!>          for a given set of energies and number of electrons.
!>          Note that singly occupied orbitals are assumed.
!>          could fail if the fermi level lies out of the range of eigenvalues
!>          (to be fixed)
!> \param   f occupations
!> \param   mu Fermi level (output)
!> \param kTS ...
!> \param   e eigenvalues
!> \param   N total number of electrons (input)
!> \param   T  electronic temperature
!> \param   maxocc maximum occupation of an orbital
!> \param   estate excited state in core level spectroscopy
!> \param   festate occupation of the excited state in core level spectroscopy
!> \date    09.2008
!> \par History
!>          - Made estate and festate optional (LT, 2014/02/26)
!> \author  Joost VandeVondele
! **************************************************************************************************
   SUBROUTINE fermifixed(f, mu, kTS, e, N, T, maxocc, estate, festate)
      REAL(kind=dp), INTENT(OUT)                         :: f(:), mu, kts
      REAL(kind=dp), INTENT(IN)                          :: e(:), n, t, maxocc
      INTEGER, INTENT(IN), OPTIONAL                      :: estate
      REAL(kind=dp), INTENT(IN), OPTIONAL                :: festate
 
      INTEGER                                            :: iter, my_estate
      REAL(kind=dp)                                      :: mu_max, mu_min, mu_now, my_festate, &
                                                            n_max, n_min, n_now
 
      IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
         my_estate = estate
         my_festate = festate
      ELSE
         my_estate = nint(maxocc)
         my_festate = my_estate
      END IF
 
! bisection search to find N
! first bracket
 
      mu_min = minval(e)
      iter = 0
      DO
         iter = iter + 1
         CALL fermi(f, n_min, kts, e, mu_min, t, maxocc, my_estate, my_festate)
         IF (n_min > n .OR. iter > 20) THEN
            mu_min = mu_min - t
         ELSE
            EXIT
         END IF
      END DO
 
      mu_max = maxval(e)
      iter = 0
      DO
         iter = iter + 1
         CALL fermi(f, n_max, kts, e, mu_max, t, maxocc, my_estate, my_festate)
         IF (n_max < n .OR. iter > 20) THEN
            mu_max = mu_max + t
         ELSE
            EXIT
         END IF
      END DO
 
      ! now bisect
      iter = 0
      DO WHILE (mu_max - mu_min > epsilon(mu)*max(1.0_dp, abs(mu_max), abs(mu_min)))
         iter = iter + 1
         mu_now = (mu_max + mu_min)/2.0_dp
         CALL fermi(f, n_now, kts, e, mu_now, t, maxocc, my_estate, my_festate)
         iter = iter + 1
 
         IF (n_now <= n) THEN
            mu_min = mu_now
         ELSE
            mu_max = mu_now
         END IF
 
         IF (iter > bisect_max_iter) THEN
            cpwarn("Maximum number of iterations reached while finding the Fermi energy")
            EXIT
         END IF
      END DO
 
      mu = (mu_max + mu_min)/2.0_dp
      CALL fermi(f, n_now, kts, e, mu, t, maxocc, my_estate, my_festate)
 
   END SUBROUTINE fermifixed
 
! **************************************************************************************************
!> \brief Bisection search to find mu for a given nel (kpoint case)
!> \param f   Occupation numbers
!> \param mu  chemical potential
!> \param kTS Entropic energy contribution
!> \param e   orbital (band) energies
!> \param nel Number of electrons (total)
!> \param wk  kpoint weights
!> \param t   Temperature
!> \param maxocc Maximum occupation
! **************************************************************************************************
   SUBROUTINE fermikp(f, mu, kTS, e, nel, wk, t, maxocc)
      REAL(kind=dp), DIMENSION(:, :), INTENT(OUT)        :: f
      REAL(kind=dp), INTENT(OUT)                         :: mu, kts
      REAL(kind=dp), DIMENSION(:, :), INTENT(IN)         :: e
      REAL(kind=dp), INTENT(IN)                          :: nel
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: wk
      REAL(kind=dp), INTENT(IN)                          :: t, maxocc
 
      REAL(kind=dp), PARAMETER                           :: epsocc = 1.0e-12_dp
 
      INTEGER                                            :: iter
      REAL(kind=dp)                                      :: de, mu_max, mu_min, n_now
 
      ! bisection search to find mu for a given nel
      de = t*log((1.0_dp - epsocc)/epsocc)
      de = max(de, 0.5_dp)
      mu_min = minval(e) - de
      mu_max = maxval(e) + de
      iter = 0
      DO WHILE (mu_max - mu_min > epsilon(mu)*max(1.0_dp, abs(mu_max), abs(mu_min)))
         iter = iter + 1
         mu = (mu_max + mu_min)/2.0_dp
         CALL fermi2(f, n_now, kts, e, mu, wk, t, maxocc)
 
         IF (abs(n_now - nel) < nel*epsocc) EXIT
 
         IF (n_now <= nel) THEN
            mu_min = mu
         ELSE
            mu_max = mu
         END IF
 
         IF (iter > bisect_max_iter) THEN
            cpwarn("Maximum number of iterations reached while finding the Fermi energy")
            EXIT
         END IF
      END DO
 
      mu = (mu_max + mu_min)/2.0_dp
      CALL fermi2(f, n_now, kts, e, mu, wk, t, maxocc)
 
   END SUBROUTINE fermikp
 
! **************************************************************************************************
!> \brief Bisection search to find mu for a given nel (kpoint case)
!> \param f   Occupation numbers
!> \param mu  chemical potential
!> \param kTS Entropic energy contribution
!> \param e   orbital (band) energies
!> \param nel Number of electrons (total)
!> \param wk  kpoint weights
!> \param t   Temperature
! **************************************************************************************************
   SUBROUTINE fermikp2(f, mu, kTS, e, nel, wk, t)
      REAL(kind=dp), DIMENSION(:, :, :), INTENT(OUT)     :: f
      REAL(kind=dp), INTENT(OUT)                         :: mu, kts
      REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN)      :: e
      REAL(kind=dp), INTENT(IN)                          :: nel
      REAL(kind=dp), DIMENSION(:), INTENT(IN)            :: wk
      REAL(kind=dp), INTENT(IN)                          :: t
 
      REAL(kind=dp), PARAMETER                           :: epsocc = 1.0e-12_dp
 
      INTEGER                                            :: iter
      REAL(kind=dp)                                      :: de, ktsa, ktsb, mu_max, mu_min, n_now, &
                                                            na, nb
 
      ! only do spin polarized case
      cpassert(SIZE(f, 3) == 2 .AND. SIZE(e, 3) == 2)
 
      ! bisection search to find mu for a given nel
      de = t*log((1.0_dp - epsocc)/epsocc)
      de = max(de, 0.5_dp)
      mu_min = minval(e) - de
      mu_max = maxval(e) + de
      iter = 0
      DO WHILE (mu_max - mu_min > epsilon(mu)*max(1.0_dp, abs(mu_max), abs(mu_min)))
         iter = iter + 1
         mu = (mu_max + mu_min)/2.0_dp
         CALL fermi2(f(:, :, 1), na, ktsa, e(:, :, 1), mu, wk, t, 1.0_dp)
         CALL fermi2(f(:, :, 2), nb, ktsb, e(:, :, 2), mu, wk, t, 1.0_dp)
         n_now = na + nb
 
         IF (abs(n_now - nel) < nel*epsocc) EXIT
 
         IF (n_now <= nel) THEN
            mu_min = mu
         ELSE
            mu_max = mu
         END IF
 
         IF (iter > bisect_max_iter) THEN
            cpwarn("Maximum number of iterations reached while finding the Fermi energy")
            EXIT
         END IF
      END DO
 
      mu = (mu_max + mu_min)/2.0_dp
      CALL fermi2(f(:, :, 1), na, ktsa, e(:, :, 1), mu, wk, t, 1.0_dp)
      CALL fermi2(f(:, :, 2), nb, ktsb, e(:, :, 2), mu, wk, t, 1.0_dp)
      kts = ktsa + ktsb
 
   END SUBROUTINE fermikp2
 
! **************************************************************************************************
!> \brief   returns f and dfde for a given set of energies and number of electrons
!>          it is a numerical derivative, trying to use a reasonable step length
!>          it ought to yield an accuracy of approximately EPSILON()^(2/3) (~10^-11)
!>          l ~ 10*T yields best accuracy
!>          Note that singly occupied orbitals are assumed.
!>          To be fixed: this could be parallellized for better efficiency
!> \param   dfde derivatives of the occupation numbers with respect to the eigenvalues
!>               the ith column is the derivative of f wrt to e_i
!> \param   f occupations
!> \param   mu Fermi level (input)
!> \param kTS ...
!> \param   e eigenvalues
!> \param   N total number of electrons (output)
!> \param   T  electronic temperature
!> \param maxocc ...
!> \param   l  typical length scale (~ 10 * T)
!> \param estate ...
!> \param festate ...
!> \date    09.2008
!> \par   History
!>          - Made estate and festate optional (LT, 2014/02/26)
!>          - Changed order of input, so l is before the two optional variables
!>            (LT, 2014/02/26)
!> \author  Joost VandeVondele
! **************************************************************************************************
   SUBROUTINE fermifixedderiv(dfde, f, mu, kTS, e, N, T, maxocc, l, estate, festate)
      REAL(kind=dp), INTENT(OUT)                         :: dfde(:, :), f(:), mu, kts
      REAL(kind=dp), INTENT(IN)                          :: e(:), n, t, maxocc, l
      INTEGER, INTENT(IN), OPTIONAL                      :: estate
      REAL(kind=dp), INTENT(IN), OPTIONAL                :: festate
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'FermiFixedDeriv'
 
      INTEGER                                            :: handle, i, my_estate, nstate
      REAL(kind=dp)                                      :: h, mux, my_festate
      REAL(kind=dp), ALLOCATABLE                         :: ex(:), fx(:)
 
      CALL timeset(routinen, handle)
 
      nstate = SIZE(e)
      ALLOCATE (ex(nstate), fx(nstate))
 
      IF (PRESENT(estate) .AND. PRESENT(festate)) THEN
         my_estate = estate
         my_festate = festate
      ELSE
         my_estate = nint(maxocc)
         my_festate = my_estate
      END IF
 
      DO i = 1, nstate
         ! NR 5.7.8
         ! the problem here is that each f_i 'seems to have' a different length scale
         ! and it would be to expensive to compute each single df_i/de_i using a finite difference
         h = (epsilon(h)**(1.0_dp/3.0_dp))*l
         ! get an exact machine representable number close to this h
         h = 2.0_dp**exponent(h)
         ! this should write three times the same number
         ! write(*,*) h,(e(i)+h)-e(i),(e(i)-h)-e(i)
         ! and the symmetric finite difference
         ex(:) = e
         ex(i) = e(i) + h
         CALL fermifixed(fx, mux, kts, ex, n, t, maxocc, my_estate, my_festate)
         dfde(:, i) = fx
         ex(i) = e(i) - h
         CALL fermifixed(fx, mux, kts, ex, n, t, maxocc, my_estate, my_festate)
         dfde(:, i) = (dfde(:, i) - fx)/(2.0_dp*h)
      END DO
      DEALLOCATE (ex, fx)
 
      CALL fermifixed(f, mu, kts, e, n, t, maxocc, my_estate, my_festate)
 
      CALL timestop(handle)
 
   END SUBROUTINE fermifixedderiv
 
END MODULE fermi_utils