smear.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
! Some pieces glued from PW/w(0,1)gauss.f90 from PWSCF
 
#include "global.h"
 
module smear_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use kpoints_oct_m
  use loct_math_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use simplex_oct_m, only: simplex_t, simplex_weights
  use sort_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
 
  implicit none
 
  private
  public ::                           &
    smear_t,                          &
    smear_init,                       &
    smear_copy,                       &
    smear_find_fermi_energy,          &
    smear_fill_occupations,           &
    smear_calc_entropy,               &
    smear_delta_function,             &
    smear_step_function,              &
    smear_entropy_function,           &
    smear_is_semiconducting,          &
    smear_write_info
 
  type smear_t
    private
    integer, public :: method
    real(real64),   public :: dsmear
    real(real64)    :: dsmear_cond
    real(real64),   public :: e_fermi
    real(real64)    :: e_fermi_cond
 
    integer, public :: el_per_state
    real(real64),   public :: ef_occ
    real(real64)    :: ef_occ_cond
    logical, public :: integral_occs
    integer, public :: MP_n
    integer         :: fermi_count
    integer         :: nik_factor
    integer         :: nspins
 
    type(simplex_t), pointer :: tetra_simplex => null()
 
    ! PhotonDoping
    logical, public:: photodop
    real(real64)   :: nephotodop
    integer        :: photodop_bandmin
 
  end type smear_t
 
  integer, parameter, public ::       &
    SMEAR_SEMICONDUCTOR     = 1,      &
    smear_fermi_dirac       = 2,      &
    smear_cold              = 3,      &
    smear_methfessel_paxton = 4,      &
    smear_spline            = 5,      &
    smear_fixed_occ         = 6,      &
    smear_gaussian          = option__smearingfunction__gaussian_smearing, &
    smear_lorentzian        = option__smearingfunction__lorentzian_smearing, &
    smear_tetrahedra        = option__smearingfunction__tetrahedra, &
    smear_tetrahedra_opt    = option__smearingfunction__tetrahedra_opt
 
  real(real64), parameter :: TOL_SMEAR = 1e-6_real64
 
contains
 
  !--------------------------------------------------
  subroutine smear_init(this, namespace, ispin, fixed_occ, integral_occs, kpoints)
    type(smear_t),             intent(out) :: this
    type(namespace_t),         intent(in)  :: namespace
    integer,                   intent(in)  :: ispin
    logical,                   intent(in)  :: fixed_occ
    logical,                   intent(in)  :: integral_occs
    type(kpoints_t),           intent(in)  :: kpoints
 
    push_sub(smear_init)
 
    this%integral_occs = integral_occs
 
    !%Variable SmearingFunction
    !%Type integer
    !%Default semiconducting
    !%Section States
    !%Description
    !% This is the function used to smear the electronic occupations.
    !% It is ignored if the <tt>Occupations</tt> block is set.
    !%Option semiconducting 1
    !% Semiconducting occupations, <i>i.e.</i> the lowest lying states are occupied
    !% until no more electrons are left.
    !%Option fermi_dirac 2
    !% Simple Fermi-Dirac distribution. In this case, <tt>Smearing</tt> has
    !% the meaning of an electronic temperature. DN Mermin, <i>Phys. Rev.</i> <b>137</b>, A1441 (1965).
    !%Option cold_smearing 3
    !% N Marzari, D Vanderbilt, A De Vita, and MC Payne, <i>Phys. Rev. Lett.</i> <b>82</b>, 3296 (1999).
    !%Option methfessel_paxton 4
    !% M Methfessel and AT Paxton, <i>Phys. Rev. B</i> <b>40</b>, 3616 (1989).
    !% In this case, the variable <tt>SmearingMPOrder</tt> sets the order of the smearing.
    !% Occupations may be negative.
    !%Option spline_smearing 5
    !% Nearly identical to Gaussian smearing.
    !% JM Holender, MJ Gillan, MC Payne, and AD Simpson, <i>Phys. Rev. B</i> <b>52</b>, 967 (1995).
    !%Option gaussian_smearing
    !% Gaussian smearing.
    !%Option lorentzian_smearing
    !% Lorentzian smearing.
    !%Option tetrahedra
    !% Linear tetrahedron method as in P. E. Bloechl, et al., <i>Phys. Rev. B</i> <b>49</b>, 16223 (1994).
    !% Requires a regular Monkhorst-Pack generated by Octopus.
    !%Option tetrahedra_opt
    !% Improved tetrahedron method as in M. Kawamura, et al., <i>Phys. Rev. B</i> <b>89</b>, 094515 (2014).
    !% Requires a regular Monkhorst-Pack generated by Octopus.
    !%End
    if (fixed_occ) then
      this%method = smear_fixed_occ
    else
      call parse_variable(namespace, 'SmearingFunction', smear_semiconductor, this%method)
      if (.not. varinfo_valid_option('SmearingFunction', this%method)) call messages_input_error(namespace, 'SmearingFunction')
      call messages_print_var_option('SmearingFunction', this%method, namespace=namespace)
    end if
 
    !%Variable Smearing
    !%Type float
    !%Default 0.1 eV
    !%Section States
    !%Description
    !% If <tt>Occupations</tt> is not set, <tt>Smearing</tt> is the
    !% smearing width used in the <tt>SmearingFunction</tt> to distribute the electrons
    !% among the existing states.
    !%End
    this%dsmear = 1e-14_real64
    if (this%method /= smear_semiconductor .and. this%method /= smear_fixed_occ) then
      call parse_variable(namespace, 'Smearing', 0.1_real64 / (m_two * p_ry), this%dsmear, units_inp%energy)
    end if
 
    !%Variable PhotoDopingSmearing
    !%Type float
    !%Default 0.1 eV
    !%Section States
    !%Description
    !% If <tt>Occupations</tt> is not set, <tt>PhotoDopingSmearing</tt> is the
    !% smearing width used in the <tt>SmearingFunction</tt> to distribute the electrons in conduction bands
    !% among the existing states.
    !%End
    this%dsmear_cond = this%dsmear
    if (this%method /= smear_semiconductor .and. this%method /= smear_fixed_occ) then
      call parse_variable(namespace, 'PhotoDopingSmearing', 0.1_real64 / (m_two * p_ry), this%dsmear_cond, units_inp%energy)
    end if
 
    !%Variable PhotoDopingNumElectrons
    !%Type float
    !%Default 0
    !%Section States
    !%Description
    !% This variable allows to set the number of valence electrons taken out to put
    !% into the conduction bands. The method follows Marini and Calandra, PRB 104, 144103 (2021).
    !%End
    call parse_variable(namespace, 'PhotoDopingNumElectrons', m_zero, this%nephotodop)
 
    if (this%nephotodop > m_min_occ) then
      this%photodop=.true.
    else
      this%photodop=.false.
    endif
 
 
    !%Variable PhotoDopingBand
    !%Type integer
    !%Default 1
    !%Section States
    !%Description
    !% This variable specifies where the minium (starting) conduction band index for the photodoping electrons.
    !%End
    call parse_variable(namespace, 'PhotoDopingBand', 1, this%photodop_bandmin)
 
 
    call messages_obsolete_variable(namespace, 'ElectronicTemperature', 'Smearing')
 
    if (ispin == 1) then ! unpolarized
      this%el_per_state = 2
    else
      this%el_per_state = 1
    end if
 
    if (ispin == 2) then
      this%nspins = 2
    else
      this%nspins = 1
    end if
 
    if (this%method == smear_semiconductor) then
      this%nik_factor = kpoints_kweight_denominator(kpoints)
 
      if (this%nik_factor == 0) then
        message(1) = "k-point weights in KPoints or KPointsReduced blocks must be rational numbers for semiconducting smearing."
        call messages_fatal(1, namespace=namespace)
      end if
    else if (this%method == smear_tetrahedra .or. this%method == smear_tetrahedra_opt) then
      if (.not. associated(kpoints%reduced%simplex)) then
        message(1) = "Tetrahedron smearing requires a simplex mesh for the k-point grid."
        call messages_fatal(1, namespace=namespace)
      end if
 
      this%tetra_simplex => kpoints%reduced%simplex
    end if
 
    this%MP_n = 0
    if (this%method == smear_methfessel_paxton) then
      !%Variable SmearingMPOrder
      !%Type integer
      !%Default 1
      !%Section States
      !%Description
      !% Sets the order of the Methfessel-Paxton smearing function.
      !%End
      call parse_variable(namespace, 'SmearingMPOrder', 1, this%MP_n)
    end if
 
    pop_sub(smear_init)
  end subroutine smear_init
 
 
  !--------------------------------------------------
  subroutine smear_copy(to, from)
    type(smear_t), intent(out) :: to
    type(smear_t), intent(in)  :: from
 
    push_sub(smear_copy)
 
    to%method       = from%method
    to%dsmear       = from%dsmear
    to%e_fermi      = from%e_fermi
    to%el_per_state = from%el_per_state
    to%ef_occ       = from%ef_occ
    to%MP_n         = from%MP_n
    to%fermi_count  = from%fermi_count
    to%nik_factor   = from%nik_factor
    to%tetra_simplex => from%tetra_simplex
 
    to%dsmear_cond  = from%dsmear_cond
    to%e_fermi_cond = from%e_fermi_cond
    to%ef_occ_cond  = from%ef_occ_cond
    to%integral_occs = from%integral_occs
    to%nspins       = from%nspins
    to%photodop     = from%photodop
    to%nephotodop   = from%nephotodop
    to%photodop_bandmin = from%photodop_bandmin
 
    pop_sub(smear_copy)
  end subroutine smear_copy
 
 
  !--------------------------------------------------
  subroutine smear_find_fermi_energy(this, namespace, eigenvalues, occupations, &
    qtot, nik, nst, kweights)
    type(smear_t), intent(inout) :: this
    type(namespace_t), intent(in) :: namespace
    real(real64), contiguous, intent(in) :: eigenvalues(:,:), occupations(:,:)
    real(real64), intent(in) :: qtot
    real(real64), contiguous, intent(in) :: kweights(:)
    integer, intent(in) :: nik, nst
 
    real(real64), parameter   :: tol = 1.0e-10_real64
    integer            :: ist, ik, iter, maxq, weight, sumq_int, sum_weight
    real(real64)       :: sumq_frac
    logical            :: conv
    real(real64),   allocatable :: eigenval_list(:)
    integer, allocatable :: k_list(:), reorder(:)
    integer            :: fermi_count_up, fermi_count_down
 
    push_sub(smear_find_fermi_energy)
 
    maxq = this%el_per_state * nst * this%nspins
    if (maxq - qtot <= -tol) then ! not enough states
      message(1) = 'Not enough states'
      write(message(2),'(6x,a,f12.6,a,i10)')'(total charge = ', qtot, &
        ' max charge = ', maxq
      call messages_fatal(2, namespace=namespace)
    end if
 
    conv = .true.
    if (this%method == smear_fixed_occ) then ! Fermi energy: last of occupied states
      ist_cycle: do ist = nst, 1, -1
        if (any(occupations(ist, :) > m_min_occ)) then
          this%e_fermi =  eigenvalues(ist, 1)
          this%ef_occ   = occupations(ist, 1) / this%el_per_state
          do ik = 2, nik
            if (eigenvalues(ist, ik) > this%e_fermi .and. occupations(ist, ik) > m_min_occ) then
              this%e_fermi  = eigenvalues(ist, ik)
              this%ef_occ   = occupations(ist, ik) / this%el_per_state
            end if
          end do
          exit ist_cycle
        end if
      end do ist_cycle
 
    else if (this%method == smear_semiconductor) then
      ! first we sort the eigenvalues
      safe_allocate(eigenval_list(1:nst * nik))
      safe_allocate(       k_list(1:nst * nik))
      safe_allocate(      reorder(1:nst * nik))
 
      iter = 1
      do ist = 1, nst
        do ik = 1, nik
          eigenval_list(iter) = eigenvalues(ist, ik)
          k_list(iter) = ik
          reorder(iter) = iter
          iter = iter + 1
        end do
      end do
 
      call sort(eigenval_list, reorder)
 
      sumq_int = floor(qtot) * this%nik_factor
      sumq_frac = qtot * this%nik_factor - sumq_int
      if ( sumq_frac + tol > this%nik_factor ) then
        sumq_int = sumq_int + this%nik_factor
        sumq_frac = sumq_frac - this%nik_factor
      end if
      if (sumq_frac < tol) sumq_frac = m_zero
 
      do iter = 1, nst * nik
        weight = int(kweights(k_list(reorder(iter))) * this%nik_factor + m_half)
        if (weight <= 0) cycle
        this%e_fermi = eigenval_list(iter)
        this%ef_occ  = (sumq_int + sumq_frac) / (weight * this%el_per_state)
 
        if (sumq_int - weight * this%el_per_state <= 0) then
          ! count how many occupied states are at the fermi level,
          ! this is required later to fill the states
          this%fermi_count = 1
          fermi_count_up = 1
          fermi_count_down = 1
          sum_weight = weight
          do
            if (iter - fermi_count_down < 1) exit
            if (abs(this%e_fermi - eigenval_list(iter - fermi_count_down)) > tol_smear) exit
            weight = int(kweights(k_list(reorder(iter - fermi_count_down))) * this%nik_factor + m_half)
            if (weight > 0) then
              sumq_int = sumq_int + weight * this%el_per_state
              sum_weight = sum_weight + weight
            end if
            fermi_count_down = fermi_count_down + 1
          end do
          do
            if (iter + fermi_count_up > nst*nik) exit
            if (abs(this%e_fermi - eigenval_list(iter + fermi_count_up)) > tol_smear) exit
            weight = int(kweights(k_list(reorder(iter + fermi_count_up))) * this%nik_factor + m_half)
            if (weight > 0) then
              sum_weight = sum_weight + weight
            end if
            fermi_count_up = fermi_count_up + 1
          end do
          this%fermi_count = fermi_count_up + fermi_count_down - 1
          this%ef_occ  = (sumq_int + sumq_frac) / (sum_weight * this%el_per_state)
          exit
        end if
 
        sumq_int = sumq_int - weight * this%el_per_state
      end do
      assert(this%ef_occ < m_one + 10.0_real64*m_epsilon)
 
      safe_deallocate_a(eigenval_list)
      safe_deallocate_a(k_list)
      safe_deallocate_a(reorder)
 
    else if (this%method == smear_tetrahedra .or. this%method == smear_tetrahedra_opt) then
      if (this%photodop) then
        call messages_not_implemented('Photodoping for SmearingFunction = tetrahedra/tetrahedra_opt')
      end if
      call bisection_find_fermi_energy_tetra(this, namespace, tol, eigenvalues, nik, nst, qtot, this%e_fermi)
    else ! bisection
      if (this%photodop) then
        ! find fermi energy for valence electrons
        call bisection_find_fermi_energy(this, namespace, this%dsmear, tol, &
          eigenvalues, kweights, nik, qtot-this%nephotodop, 1, this%photodop_bandmin-1, &
          this%e_fermi, this%ef_occ)
 
        ! find fermi energy for conduction electrons
        call bisection_find_fermi_energy(this, namespace, this%dsmear_cond, tol, &
          eigenvalues, kweights, nik, this%nephotodop, this%photodop_bandmin, nst, &
          this%e_fermi_cond, this%ef_occ_cond)
      else
        call bisection_find_fermi_energy(this, namespace, this%dsmear, tol, &
          eigenvalues, kweights, nik, qtot, 1, nst, this%e_fermi, this%ef_occ)
      end if
    end if
 
    pop_sub(smear_find_fermi_energy)
  end subroutine smear_find_fermi_energy
 
  ! ---------------------------------------------------------
  subroutine bisection_find_fermi_energy(this, namespace, dsmear_in, tol, &
    eigenvalues, kweights, nik, q_in, start_band, end_band, e_fermi, ef_occ)
    type(smear_t), intent(inout) :: this
    type(namespace_t), intent(in) :: namespace
    real(real64), intent(in) :: dsmear_in
    real(real64), intent(in) :: tol
    integer, intent(in) :: nik
    real(real64), contiguous, intent(in) :: eigenvalues(:,:)
    real(real64), intent(in) :: q_in
    real(real64), contiguous, intent(in) :: kweights(:)
    integer, intent(in) :: start_band, end_band
    real(real64), intent(out) :: e_fermi, ef_occ
 
    integer, parameter :: nitmax = 200
    integer            :: ist, ik, iter
    real(real64)       :: drange, dsmear, emin, emax, xx, sumq
    logical            :: conv
 
    push_sub(bisection_find_fermi_energy)
 
    dsmear = max(1e-14_real64, dsmear_in)
    drange = dsmear * sqrt(-log(tol * 0.01_real64))
 
    emin = minval(eigenvalues) - drange
    emax = maxval(eigenvalues) + drange
 
    do iter = 1, nitmax
      e_fermi = m_half * (emin + emax)
      sumq         = m_zero
 
      do ik = 1, nik
        do ist = start_band, end_band
          xx   = (e_fermi - eigenvalues(ist, ik))/dsmear
          sumq = sumq + kweights(ik) * this%el_per_state * &
            smear_step_function(this, xx)
        end do
      end do
 
      conv = (abs(sumq - q_in) <= tol)
      if (conv) exit
 
      if (sumq <= q_in) emin = e_fermi
      if (sumq >= q_in) emax = e_fermi
 
      ef_occ = smear_step_function(this, m_zero)
    end do
 
    if (.not. conv) then
      message(1) = 'Fermi: did not converge.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    pop_sub(bisection_find_fermi_energy)
  end subroutine bisection_find_fermi_energy
 
  ! ---------------------------------------------------------
  subroutine bisection_find_fermi_energy_tetra(this, namespace, tol, eigenvalues, nik, nst, q_in, e_fermi)
    type(smear_t), intent(inout) :: this
    type(namespace_t), intent(in) :: namespace
    real(real64), intent(in) :: tol
    integer, intent(in) :: nik, nst
    real(real64), contiguous, intent(in) :: eigenvalues(:,:)
    real(real64), intent(in) :: q_in
    real(real64), intent(out) :: e_fermi
 
    integer, parameter :: nitmax = 200
    integer :: iter
    real(real64) :: emin, emax, sumq
    logical :: conv
 
    integer :: ist, ii, ll, is, ik, ns, nik_red, nvertices
    real(real64) :: e_simplex(20), w_simplex(4), dos_simplex(4)
 
    push_sub_with_profile(bisection_find_fermi_energy_tetra)
 
    assert(associated(this%tetra_simplex))
    assert(this%tetra_simplex%n_simplices > 0)
    assert(this%tetra_simplex%sdim > 0)
 
    ns = this%nspins
    assert(ns > 0)
    assert(mod(nik, ns) == 0)
 
    nik_red = nik / ns
    assert(maxval(this%tetra_simplex%simplices) <= nik_red)
 
    nvertices = this%tetra_simplex%rdim + 1
    assert(nvertices > 0)
    assert(nvertices <= this%tetra_simplex%sdim)
 
    emin = minval(eigenvalues)
    emax = maxval(eigenvalues)
    if (.not. (emin < emax)) then
      write(message(1), '(a)') 'Tetra bisection: invalid energy bracket from eigenvalues.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    do iter = 1, nitmax
      e_fermi = m_half * (emin + emax)
      sumq = m_zero
      !$omp parallel do collapse(2) private(ist, ii, is, ll, ik, e_simplex, w_simplex, dos_simplex) reduction(+:sumq)
      do ist = 1, nst
        do ii = 1, this%tetra_simplex%n_simplices
          do is = 0, ns - 1
            do ll = 1, this%tetra_simplex%sdim
              ik = ns * (this%tetra_simplex%simplices(ii, ll) - 1) + 1
              e_simplex(ll) = eigenvalues(ist, ik + is)
            end do
            w_simplex(:) = m_zero
 
            call simplex_weights(this%tetra_simplex%rdim, e_simplex(1:this%tetra_simplex%sdim), e_fermi, &
              w_simplex(1:nvertices), dos_simplex(1:nvertices))
 
            do ll = 1, nvertices
              sumq = sumq + this%el_per_state * w_simplex(ll) / this%tetra_simplex%n_points
            end do
          end do
        end do
      end do
      !$omp end parallel do
 
      conv = (abs(sumq - q_in) <= tol)
      if (conv) exit
 
      if (sumq <= q_in) emin = e_fermi
      if (sumq >= q_in) emax = e_fermi
    end do
 
    if (.not. conv) then
      message(1) = 'Fermi (tetrahedra): did not converge.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    pop_sub_with_profile(bisection_find_fermi_energy_tetra)
  end subroutine bisection_find_fermi_energy_tetra
 
  ! ---------------------------------------------------------
  subroutine smear_fill_occupations(this, eigenvalues, occupations, kweights, nik, nst)
    type(smear_t), intent(in) :: this
    real(real64), contiguous, intent(in) :: eigenvalues(:,:)
    real(real64), contiguous, intent(inout) :: occupations(:,:)
    real(real64), contiguous, intent(in) :: kweights(:)
    integer, intent(in) :: nik, nst
 
    integer :: ik, ist, ifermi
    real(real64)   :: dsmear, xx, dsmear_cond
 
    push_sub(smear_fill_occupations)
 
    if (this%method == smear_fixed_occ) then
      ! do nothing
    else if (this%method == smear_semiconductor) then
      assert(this%fermi_count > 0 .and. this%fermi_count <= nik*nst)
 
      ifermi = 0
      do ik = 1, nik
        do ist = 1, nst
          xx = eigenvalues(ist, ik) - this%e_fermi
          if (xx < -tol_smear) then
            occupations(ist, ik) = this%el_per_state
          else if (abs(xx) <= tol_smear .and. ifermi < this%fermi_count) then
            occupations(ist, ik) = this%ef_occ * this%el_per_state
            ifermi = ifermi + 1
          else
            occupations(ist, ik) = m_zero
          end if
 
        end do
      end do
    else if (this%method == smear_tetrahedra .or. this%method == smear_tetrahedra_opt) then
      if (this%photodop) then
        call messages_not_implemented('Photodoping for SmearingFunction = tetrahedra/tetrahedra_opt')
      end if
      block
        integer :: ii, ll, is, ns, nik_red, nvertices
        real(real64) :: e_simplex(20), w_simplex(4), dos_simplex(4)
 
        assert(associated(this%tetra_simplex))
        assert(this%tetra_simplex%n_simplices > 0)
        assert(this%tetra_simplex%sdim > 0)
 
        ns = this%nspins
        assert(ns > 0)
        assert(mod(nik, ns) == 0)
 
        nik_red = nik / ns
        assert(maxval(this%tetra_simplex%simplices) <= nik_red)
 
        nvertices = this%tetra_simplex%rdim + 1
        assert(nvertices > 0)
        assert(nvertices <= this%tetra_simplex%sdim)
 
        occupations(:, :) = m_zero
        !$omp parallel do collapse(2) private(ist, ii, is, ll, ik, e_simplex, w_simplex, dos_simplex) reduction(+:occupations)
        do ist = 1, nst
          do ii = 1, this%tetra_simplex%n_simplices
            do is = 0, ns - 1
              do ll = 1, this%tetra_simplex%sdim
                ik = ns * (this%tetra_simplex%simplices(ii, ll) - 1) + 1
                e_simplex(ll) = eigenvalues(ist, ik + is)
              end do
              w_simplex(:) = m_zero
 
              call simplex_weights(this%tetra_simplex%rdim, e_simplex(1:this%tetra_simplex%sdim), this%e_fermi, &
                w_simplex(1:nvertices), dos_simplex(1:nvertices))
 
              do ll = 1, nvertices
                ik = ns * (this%tetra_simplex%simplices(ii, ll) - 1) + 1
                assert(kweights(ik + is) > m_epsilon)
                occupations(ist, ik + is) = occupations(ist, ik + is) + this%el_per_state * &
                  w_simplex(ll) / (this%tetra_simplex%n_points * kweights(ik + is))
              end do
            end do
          end do
        end do
        !$omp end parallel do
      end block
    else
      dsmear = max(1e-14_real64, this%dsmear)
      if (this%photodop) then
        dsmear_cond = max(1e-14_real64, this%dsmear_cond)
        do ik = 1, nik
          do ist = 1, nst
            if (ist < this%photodop_bandmin) then
              ! valence electrons
              xx = (this%e_fermi - eigenvalues(ist, ik))/dsmear
            else
              ! conduction electrons
              xx = (this%e_fermi_cond - eigenvalues(ist, ik))/dsmear_cond
            end if
            occupations(ist, ik) = smear_step_function(this, xx) * this%el_per_state
          end do
        end do
      else
        do ik = 1, nik
          do ist = 1, nst
            xx = (this%e_fermi - eigenvalues(ist, ik))/dsmear
            occupations(ist, ik) = smear_step_function(this, xx) * this%el_per_state
          end do
        end do
      end if
    end if
 
    pop_sub(smear_fill_occupations)
  end subroutine smear_fill_occupations
 
  !--------------------------------------------------
  real(real64) function smear_calc_entropy(this, eigenvalues, &
    nik, nst, kweights, occ) result(entropy)
    type(smear_t), intent(in) :: this
    real(real64), contiguous, intent(in) :: eigenvalues(:,:)
    real(real64), contiguous, intent(in) :: kweights(:)
    integer, intent(in) :: nik, nst
    real(real64), contiguous, intent(in) :: occ(:, :)
 
    integer :: ist, ik
    real(real64) :: dsmear, xx, term, ff
 
    push_sub(smear_calc_entropy)
 
    entropy = m_zero
 
    if (this%method == smear_fixed_occ .or. this%method == smear_semiconductor) then
      ! Fermi-Dirac entropy, not quite the same as will be obtained with true smearing
      ! RM Wentzcovitch, JL Martins, and PB Allen, Phys. Rev. B 45, 11372 (1992) eqn (5)
      ! also N Marzari PhD thesis p 117, http://quasiamore.mit.edu/phd/Marzari_PhD.pdf
      do ik = 1, nik
        do ist = 1, nst
          ff = occ(ist, ik) / this%el_per_state
          if (ff > m_zero .and. ff  <  m_one) then
            term = ff * log(ff) + (1 - ff) * log(1 - ff)
          else ! we have semiconducting smearing, or perverse occupations as in Methfessel-Paxton
            term = m_zero
          end if
          entropy = entropy - kweights(ik) * this%el_per_state * term
        end do
      end do
    else
      dsmear = max(1e-14_real64, this%dsmear)
 
      do ik = 1, nik
        do ist = 1, nst
          if (eigenvalues(ist, ik) < huge(m_one)) then
            xx = (this%e_fermi - eigenvalues(ist, ik)) / dsmear
            entropy = entropy - kweights(ik) * this%el_per_state *  &
              smear_entropy_function(this, xx)
          end if
        end do
      end do
    end if
 
    pop_sub(smear_calc_entropy)
  end function smear_calc_entropy
 
 
  ! ---------------------------------------------------------
  real(real64) function smear_delta_function(this, xx) result(deltaf)
    type(smear_t), intent(in) :: this
    real(real64),  intent(in) ::  xx
 
    real(real64), parameter :: maxarg = 200.0_real64
    real(real64) :: xp, arg, hd, hp, aa
    integer :: ii, ni
 
    ! no PUSH_SUB, called too often
 
    ! smear_delta_function is not defined for SMEAR_FIXED_OCC
    assert(this%method /= smear_fixed_occ)
 
    deltaf = m_zero
    select case (this%method)
    case (smear_semiconductor)
      if (abs(xx) <= m_epsilon) then
        deltaf = this%ef_occ
      end if
 
    case (smear_fermi_dirac)
      if (abs(xx) <= 36.0_real64) then
        deltaf = m_one / (m_two + exp(-xx) + exp(xx))
      end if
 
    case (smear_cold)
      xp  = xx - m_one / sqrt(m_two)
      arg = min(maxarg, xp**2)
 
      deltaf = exp(-arg) / sqrt(m_pi) * (m_two - sqrt(m_two) * xx)
 
    case (smear_methfessel_paxton)
      arg    = min(maxarg, xx**2)
      deltaf = exp(-arg) / sqrt(m_pi)
 
      if (this%MP_n > 0) then ! recursion
        hd = m_zero
        hp = exp(-arg)
        ni = 0
        aa = m_one / sqrt(m_pi)
        do ii = 1, this%MP_n
          hd = m_two * xx * hp - m_two * ni * hd
          ni = ni + 1
          aa = -aa / (m_four * ii)
          hp = m_two * xx * hd - m_two * ni * hp
          ni = ni + 1
          deltaf = deltaf + aa * hp
        end do
      end if
 
    case (smear_spline)
      xp     = abs(xx) + m_one / sqrt(m_two)
      deltaf = sqrt(m_e) * xp * exp(-xp * xp)
 
    case (smear_gaussian)
      deltaf = exp(-xx * xx) / sqrt(m_pi)
 
    case (smear_lorentzian)
      deltaf = m_one / (m_pi * (xx * xx + m_one))
 
    end select
 
  end function smear_delta_function
 
 
  ! ---------------------------------------------------------
  real(real64) function smear_step_function(this, xx) result(stepf)
    type(smear_t), intent(in) :: this
    real(real64),  intent(in) ::  xx
 
    real(real64), parameter :: maxarg = 200.0_real64
    real(real64) :: xp, arg, hd, hp, aa
    integer :: ii, ni
 
    push_sub(smear_step_function)
 
    ! smear_step_function is not defined for SMEAR_FIXED_OCC
    assert(this%method /= smear_fixed_occ)
 
    stepf = m_zero
    select case (this%method)
    case (smear_semiconductor)
      if (xx > m_zero) then
        stepf = m_one
      else if (abs(xx) <= m_epsilon) then
        stepf = this%ef_occ
      end if
 
    case (smear_fermi_dirac)
      if (xx > maxarg) then
        stepf = m_one
      else if (xx > -maxarg) then
        stepf = m_one / (m_one + exp(-xx))
      end if
 
    case (smear_cold)
      xp  = xx - m_one / sqrt(m_two)
      arg = min(maxarg, xp**2)
 
      stepf = m_half * erf(xp) + &
        m_one / sqrt(m_two * m_pi) * exp(-arg) + m_half
 
    case (smear_methfessel_paxton)
      stepf = m_half * erfc(-xx)
 
      if (this%MP_n > 0) then ! recursion
        hd = m_zero
        arg = min(maxarg, xx**2)
        hp = exp(-arg)
        ni = 0
        aa = m_one / sqrt(m_pi)
        do ii = 1, this%MP_n
          hd = m_two * xx * hp - m_two * ni * hd
          ni = ni + 1
          aa = -aa / (m_four * ii)
          stepf = stepf - aa * hd
          hp = m_two * xx * hd - m_two * ni * hp
          ni = ni + 1
        end do
      end if
 
    case (smear_spline)
      if (xx <= m_zero) then
        xp = xx - m_one / sqrt(m_two)
        stepf = m_half * sqrt(m_e) * exp(-xp * xp)
      else
        xp = xx + m_one / sqrt(m_two)
        stepf = m_one - m_half * sqrt(m_e) * exp(-xp * xp)
      end if
 
    case (smear_gaussian)
      stepf = m_half * erfc(-xx)
 
    case (smear_lorentzian)
      stepf = atan(xx) / m_pi + m_half
 
    end select
 
    pop_sub(smear_step_function)
  end function smear_step_function
 
 
  ! ---------------------------------------------------------
  real(real64) function smear_entropy_function(this, xx) result(entropyf)
    type(smear_t), intent(in) :: this
    real(real64),  intent(in) ::  xx
 
    real(real64), parameter :: maxarg = 200.0_real64
    real(real64) :: xp, arg, hd, hp, hpm1, aa
    integer :: ii, ni
 
    push_sub(smear_entropy_function)
 
    ! smear_entropy_function is not defined for SMEAR_FIXED_OCC
    assert(this%method /= smear_fixed_occ)
 
    entropyf = m_zero
    select case (this%method)
    case (smear_semiconductor)
 
    case (smear_fermi_dirac)
      if (abs(xx) <= 36.0_real64) then
        xp = m_one / (m_one + exp(-xx))
        entropyf = xp * log(xp) + (m_one - xp) * log(m_one - xp)
      end if
 
    case (smear_cold)
      xp  = xx - m_one / sqrt(m_two)
      arg = min(maxarg, xp**2)
 
      entropyf =  m_one / sqrt(m_two * m_pi) * xp * exp(-arg)
 
    case (smear_methfessel_paxton)
      arg = min(maxarg, xx**2)
      entropyf = -m_half * exp(-arg) / sqrt(m_pi)
 
      if (this%MP_n > 0) then ! recursion
        hd = m_zero
        hp = exp(-arg)
        ni = 0
        aa = m_one / sqrt(m_pi)
        do ii = 1, this%MP_n
          hd = m_two * xx * hp - m_two * ni * hd
          ni = ni + 1
          hpm1 = hp
          hp = m_two * xx * hd - m_two * ni * hp
          ni = ni + 1
          aa = -aa / (m_four * ii)
          entropyf = entropyf - aa * (m_half * hp + hpm1 * ni)
        end do
      end if
 
    case (smear_spline)
      xp = abs(xx) + m_one / sqrt(m_two)
      entropyf = -sqrt(m_e) * (abs(xx) * exp(-xp * xp) / m_two + sqrt(m_pi) / m_four * erfc(xp))
 
    case (smear_gaussian)
      entropyf = -exp(-xx * xx) / m_two / sqrt(m_pi)
 
    case (smear_lorentzian)
      entropyf = m_ninf
 
    end select
 
    pop_sub(smear_entropy_function)
  end function smear_entropy_function
 
 
  ! ---------------------------------------------------------
  logical pure function smear_is_semiconducting(this) result(answer)
    type(smear_t), intent(in) :: this
 
    answer = this%method == smear_semiconductor
 
  end function smear_is_semiconducting
 
  subroutine smear_write_info(this, namespace, iunit)
    type(smear_t),     intent(in) :: this
    type(namespace_t), intent(in) :: namespace
    integer, optional, intent(in) :: iunit
 
 
    if (this%method /= smear_semiconductor .and. this%method /= smear_fixed_occ) then
      if (this%photodop) then
        write(message(1), '(a,f12.6,1x,a)') "Fermi energy (valence   ) = ", &
          units_from_atomic(units_out%energy, this%e_fermi), units_abbrev(units_out%energy)
        write(message(2), '(a,f12.6,1x,a)') "Fermi energy (conduction) = ", &
          units_from_atomic(units_out%energy, this%e_fermi_cond), units_abbrev(units_out%energy)
        call messages_info(2, iunit=iunit, namespace=namespace)
      else
        write(message(1), '(a,f12.6,1x,a)') "Fermi energy = ", &
          units_from_atomic(units_out%energy, this%e_fermi), units_abbrev(units_out%energy)
        call messages_info(1, iunit=iunit, namespace=namespace)
      end if
    end if
 
  end subroutine smear_write_info
 
 
end module smear_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: