utility_wgauss Function

  1. utility.F90
  2. w90_utility
  3. utility_wgauss

public function utility_wgauss(x, n)

Uses

  • proc~~utility_wgauss~~UsesGraph proc~utility_wgauss utility_wgauss module~w90_constants w90_constants proc~utility_wgauss->module~w90_constants

this function computes the approximate theta function for the given order n, at the point x.

(n>=0) : Methfessel-Paxton case. See PRB 40, 3616 (1989).

(n=-1 ): Cold smearing (Marzari-Vanderbilt). See PRL 82, 3296 (1999) 1/2erf(x-1/sqrt(2)) + 1/sqrt(2pi)exp(-(x-1/sqrt(2))*2) + 1/2

(n=-99): Fermi-Dirac case: 1.0/(1.0+exp(-x)).

Arguments

Type IntentOptional AttributesName
real(kind=dp) :: x
integer :: n

Return Value real(kind=dp)


Calls

proc~~utility_wgauss~~CallsGraph proc~utility_wgauss utility_wgauss proc~qe_erf qe_erf proc~utility_wgauss->proc~qe_erf proc~gauss_freq gauss_freq proc~utility_wgauss->proc~gauss_freq proc~qe_erfc qe_erfc proc~qe_erf->proc~qe_erfc proc~gauss_freq->proc~qe_erfc proc~qe_erfc->proc~qe_erf

Contents

Source Code

utility_wgauss

Source Code

  function utility_wgauss(x, n)
    !-----------------------------------------------------------------------
    !
    !! this function computes the approximate theta function for the
    !! given order n, at the point x.
    !!
    !! (n>=0) : Methfessel-Paxton case. See PRB 40, 3616 (1989).
    !!
    !! (n=-1 ): Cold smearing (Marzari-Vanderbilt). See PRL 82, 3296 (1999)
    !!       1/2*erf(x-1/sqrt(2)) + 1/sqrt(2*pi)*exp(-(x-1/sqrt(2))**2) + 1/2
    !!
    !! (n=-99): Fermi-Dirac case: 1.0/(1.0+exp(-x)).
    !
    use w90_constants, only: dp, pi

    implicit none
    real(kind=dp) :: utility_wgauss, x
    !! output: the value of the function
    !! input: the argument of the function
    integer :: n
    !! input: the order of the function
    !
    !    the local variables
    !

    real(kind=dp) :: a, hp, arg, hd, xp
    ! the coefficient a_n
    ! the hermitean function
    ! the argument of the exponential
    ! the hermitean function
    ! auxiliary variable (cold smearing)
    integer :: i, ni
    ! counter on the n indices
    ! counter on 2n
!    real(kind=dp), external :: gauss_freq, qe_erf
    real(kind=dp), parameter :: maxarg = 200.0_dp
    ! maximum value for the argument of the exponential

    ! Fermi-Dirac smearing
    if (n .eq. -99) then
      if (x .lt. -maxarg) then
        utility_wgauss = 0.0_dp
      elseif (x .gt. maxarg) then
        utility_wgauss = 1.0_dp
      else
        utility_wgauss = 1.00_dp/(1.00_dp + exp(-x))
      endif
      return

    endif
    ! Cold smearing
    if (n .eq. -1) then
      xp = x - 1.00_dp/sqrt(2.00_dp)
      arg = min(maxarg, xp**2)
      utility_wgauss = 0.50_dp*qe_erf(xp) + 1.00_dp/sqrt(2.00_dp*pi)*exp(- &
                                                                         arg) + 0.50_dp
      return

    endif
    ! Methfessel-Paxton
    utility_wgauss = gauss_freq(x*sqrt(2.00_dp))
    if (n .eq. 0) return
    hd = 0.0_dp
    arg = min(maxarg, x**2)
    hp = exp(-arg)
    ni = 0
    a = 1.0_dp/sqrt(pi)
    do i = 1, n
      hd = 2.00_dp*x*hp - 2.00_dp*DBLE(ni)*hd
      ni = ni + 1
      a = -a/(DBLE(i)*4.00_dp)
      utility_wgauss = utility_wgauss - a*hd
      hp = 2.00_dp*x*hd - 2.00_dp*DBLE(ni)*hp
      ni = ni + 1
    enddo
    return
  end function utility_wgauss