simplex.F90

Go to the documentation of this file.

!! Copyright (C) 2025 Octopus Developers
!!
!! 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.
!!
 
#include "global.h"
 
module simplex_oct_m
  use, intrinsic :: iso_fortran_env, only: real64
  use debug_oct_m, only: debug
#ifndef NDEBUG
  use debug_oct_m, only: debug_pop_sub, debug_push_sub
#endif
  use global_oct_m, only: &
    m_one, m_two, m_zero, global_alloc_err, global_alloc_errmsg, global_sizeof, not_in_openmp, operator(+), stderr
  use messages_oct_m, only: alloc_error, dealloc_error, message, messages_fatal, messages_not_implemented
  use profiling_oct_m, only: &
    prof_vars, profiling_in, profiling_memory, profiling_memory_allocate, profiling_memory_deallocate, profiling_out
  use sort_oct_m, only: sort
 
  implicit none
 
  private
 
  public :: &
    simplex_t, &
    simplex_init, &
    simplex_end, &
    simplex_weights, &
    simplex_dos, &
    simplex_weights_1d, &
    simplex_dos_1d, &
    simplex_weights_2d, &
    simplex_dos_2d, &
    simplex_weights_3d, &
    simplex_dos_3d
 
  type simplex_t
    integer :: rdim
    integer :: sdim
    integer :: n_simplices
    integer :: n_points
    integer, allocatable :: simplices(:,:)
  contains
    final :: simplex_end
  end type simplex_t
 
  interface simplex_t
    module procedure simplex_init
  end interface simplex_t
 
contains
 
  function simplex_init(dim, naxis, nshifts, shift, kpoints, equiv, opt) result(this)
    integer, intent(in) :: dim
    integer, intent(in) :: naxis(1:dim)
    integer, intent(in) :: nshifts
    real(real64), intent(in) :: shift(:,:)
    real(real64), intent(in) :: kpoints(:,:)
    integer, intent(in), optional :: equiv(:)
    logical, intent(in) :: opt
    type(simplex_t), pointer :: this
 
    real(real64) :: kmin(dim)
    integer :: ik, npoints
 
    integer :: ix(dim)
    integer, allocatable :: kl123(:,:,:)
 
    integer :: rdim, raxis(3)
 
    push_sub(simplex_init)
 
    if (nshifts /= 1) then
      message(1) = "The linear tetrahedron method only works for automatic k-point grids with a single shift"
      call messages_fatal(1)
    end if
 
    safe_allocate(this)
    safe_allocate_source(kl123(1:naxis(1), 1:naxis(2), 1:naxis(3)), -1)
 
    npoints = product(naxis)
    kmin = minval(kpoints, 2)
 
    do ik = 1, npoints
      ix(:) = nint((kpoints(:,ik) - kmin) * naxis + 1)
      assert(kl123(ix(1), ix(2), ix(3)) == -1)
      if (present(equiv)) then
        kl123(ix(1), ix(2), ix(3)) = equiv(ik)
      else
        kl123(ix(1), ix(2), ix(3)) = ik
      end if
    end do
 
    rdim = sum(merge(1, 0, naxis > 1))
    raxis(1:rdim) = pack(naxis, naxis > 1)
 
    if (any(raxis(1:rdim) /= naxis(1:rdim))) then
      message(1) = "The periodic dimensions must be consecutive"
      call messages_fatal(1)
    end if
 
    select case (rdim)
    case default
      assert(.false.)
    case (1)
      block
        integer, parameter :: submesh_segments(1,2) = reshape([ &
          1, 2 ], shape(submesh_segments), order=[2, 1])
        ! coordinates of corners and neighbors in barycentric coordinates
        integer, parameter :: b(4,2) = reshape([ &
          1 , 0 , & ! k1
          0 , 1 , & ! k2
          2 , -1, & ! 2 * k1 - k2
          -1, 2   & ! 2 * k2 - k1
          ], shape(b), order=[2, 1])
 
        integer :: i, ip1, it, n
        integer :: corners(2,1), v(2,1), c(1)
        integer :: this_segment(2), this_corner
 
        this%n_points = npoints
        this%n_simplices = npoints
        this%rdim = rdim
        this%sdim = merge(4, 2, opt)
        safe_allocate(this%simplices(this%n_simplices, this%sdim))
 
        do i = 1, raxis(1)
          ip1 = modulo(i, raxis(1)) + 1
          corners(:,:) = reshape([ i , ip1 ], shape(corners), order=[2, 1])
 
          do it = 1, size(submesh_segments, 1)
            n = (it - 1) + 1 * (i - 1) + 1
            this_segment(:) = submesh_segments(it, :)
            v(1,:) = corners(this_segment(1), :)
            v(2,:) = corners(this_segment(2), :)
            do ik = 1, this%sdim
              c(:) = b(ik,1) * v(1,:) + b(ik,2) * v(2,:)
              c(:) = modulo(c(:) - 1, raxis(1:rdim)) + 1
              this_corner = kl123(c(1), 1, 1)
              this%simplices(n,ik) = this_corner
            end do
          end do
        end do
      end block
    case (2)
      block
        integer, parameter :: submesh_triangles(2,3) = reshape([ &
          1, 2, 3, &
          1, 4, 3], shape(submesh_triangles), order=[2, 1])
        ! coordinates of corners and neighbors in barycentric coordinates
        integer, parameter :: b(10,3) = reshape([ &
          1 , 0 , 0 , & ! k1
          0 , 1 , 0 , & ! k2
          0 , 0 , 1 , & ! k3
          2 , -1, 0 , & ! 2 * k1 - k2
          0 , 2 , -1, & ! 2 * k2 - k3
          2 , 0 , -1, & ! 2 * k1 - k3
          -1, 0 , 2 , & ! 2 * k3 - k1
          -1, 2 , 0 , & ! 2 * k2 - k1
          0 , -1, 2 , & ! 2 * k3 - k2
          1 , -1, 1   & ! k1 - k2 + k3
          ], shape(b), order=[2, 1])
 
        integer :: i, j, ip1, jp1, it, n
        integer :: corners(4,2), v(3,2), c(2)
        integer :: this_triangle(3), this_corner
 
        this%n_points = npoints
        this%n_simplices = 2 * npoints
        this%rdim = rdim
        this%sdim = merge(10, 3, opt)
        safe_allocate(this%simplices(this%n_simplices, this%sdim))
 
        do i = 1, raxis(1)
          do j = 1, raxis(2)
            ip1 = modulo(i, raxis(1)) + 1
            jp1 = modulo(j, raxis(2)) + 1
            corners(:,:) = reshape([ &
              i    , j    , &
              ip1  , j    , &
              ip1  , jp1  , &
              i    , jp1  ], shape(corners), order=[2, 1])
 
            do it = 1, size(submesh_triangles, 1)
              n = (it - 1) + 2 * ((j - 1) + raxis(2) * (i - 1)) + 1
              this_triangle(:) = submesh_triangles(it, :)
              v(1,:) = corners(this_triangle(1), :)
              v(2,:) = corners(this_triangle(2), :)
              v(3,:) = corners(this_triangle(3), :)
              do ik = 1, this%sdim
                c(:) = b(ik,1) * v(1,:) + b(ik,2) * v(2,:) + b(ik,3) * v(3,:)
                c(:) = modulo(c(:) - 1, raxis(1:rdim)) + 1
                this_corner = kl123(c(1), c(2), 1)
                this%simplices(n,ik) = this_corner
              end do
            end do
          end do
        end do
      end block
    case(3)
      block
        integer, parameter :: submesh_tetras(6,4) = reshape([ &
          1, 2, 3, 6, &
          1, 3, 5, 6, &
          3, 5, 6, 7, &
          3, 6, 7, 8, &
          3, 4, 6, 8, &
          2, 3, 4, 6], shape(submesh_tetras), order=[2, 1])
        ! coordinates of corners and neighbors in barycentric coordinates
        integer, parameter :: b(20,4) = reshape([ &
          1 , 0 , 0 , 0 , & ! k1
          0 , 1 , 0 , 0 , & ! k2
          0 , 0 , 1 , 0 , & ! k3
          0 , 0 , 0 , 1 , & ! k4
          2 , -1, 0 , 0 , & ! 2 * k1 - k2
          0 , 2 , -1, 0 , & ! 2 * k2 - k3
          0 , 0 , 2 , -1, & ! 2 * k3 - k4
          -1, 0 , 0 , 2 , & ! 2 * k4 - k1
          2 , 0 , -1, 0 , & ! 2 * k1 - k3
          0 , 2 , 0 , -1, & ! 2 * k2 - k4
          -1, 0 , 2 , 0 , & ! 2 * k3 - k1
          0 , -1, 0 , 2 , & ! 2 * k4 - k2
          2 , 0 , 0 , -1, & ! 2 * k1 - k4
          -1, 2 , 0 , 0 , & ! 2 * k2 - k1
          0 , -1, 2 , 0 , & ! 2 * k3 - k2
          0 , 0 , -1, 2 , & ! 2 * k4 - k3
          -1, 1 , 0 , 1 , & ! k4 - k1 + k2
          1 , -1, 1 , 0 , & ! k1 - k2 + k3
          0 , 1 , -1, 1 , & ! k2 - k3 + k4
          1 , 0 , 1 , -1  & ! k3 - k4 + k1
          ], shape(b), order=[2, 1])
 
        integer :: i, j, k, ip1, jp1, kp1, it, n
        integer :: corners(8,3), v(4,3), c(3)
        integer :: this_tetra(4), this_corner
 
        this%n_points = npoints
        this%n_simplices = 6 * npoints
        this%rdim = rdim
        this%sdim = merge(20, 4, opt)
        safe_allocate(this%simplices(this%n_simplices, this%sdim))
 
        do i = 1, raxis(1)
          do j = 1, raxis(2)
            do k = 1, raxis(3)
              ip1 = modulo(i, raxis(1)) + 1
              jp1 = modulo(j, raxis(2)) + 1
              kp1 = modulo(k, raxis(3)) + 1
              corners(:,:) = reshape([ &
                i    , j    , k    , &
                ip1  , j    , k    , &
                i    , jp1  , k    , &
                ip1  , jp1  , k    , &
                i    , j    , kp1  , &
                ip1  , j    , kp1  , &
                i    , jp1  , kp1  , &
                ip1  , jp1  , kp1  ], shape(corners), order=[2, 1])
 
              do it = 1, size(submesh_tetras, 1)
                n = (it - 1) + 6 * ((k - 1) + raxis(3) * ((j - 1) + raxis(2) * (i - 1))) + 1
                this_tetra(:) = submesh_tetras(it, :)
                v(1,:) = corners(this_tetra(1), :)
                v(2,:) = corners(this_tetra(2), :)
                v(3,:) = corners(this_tetra(3), :)
                v(4,:) = corners(this_tetra(4), :)
                do ik = 1, this%sdim
                  c(:) = b(ik,1) * v(1,:) + b(ik,2) * v(2,:) + b(ik,3) * v(3,:) + b(ik,4) * v(4,:)
                  c(:) = modulo(c(:) - 1, raxis(1:rdim)) + 1
                  this_corner = kl123(c(1), c(2), c(3))
                  this%simplices(n,ik) = this_corner
                end do
              end do
            end do
          end do
        end do
      end block
    end select
 
    safe_deallocate_a(kl123)
 
    pop_sub(simplex_init)
  end function simplex_init
 
  subroutine simplex_end(this)
    type(simplex_t), intent(inout) :: this
    push_sub(simplex_end)
    safe_deallocate_a(this%simplices)
    pop_sub(simplex_end)
  end subroutine simplex_end
 
  subroutine simplex_weights(rdim, esimplex, eF, weights, dos)
    integer, intent(in) :: rdim
    real(real64), intent(in) :: esimplex(:)
    real(real64), intent(in) :: ef
    real(real64), intent(out) :: weights(:)
    real(real64), intent(out) :: dos(:)
 
    ! no PUSH_SUB, called too often
 
    select case (rdim)
    case (1)
      call simplex_weights_1d(esimplex, ef, weights, dos)
    case (2)
      call simplex_weights_2d(esimplex, ef, weights, dos)
    case (3)
      call simplex_weights_3d(esimplex, ef, weights, dos)
    case default
      assert(.false.)
    end select
  end subroutine simplex_weights
 
  subroutine simplex_dos(rdim, esimplex, eF, dos)
    integer, intent(in) :: rdim
    real(real64), intent(in) :: esimplex(:)
    real(real64), intent(in) :: ef
    real(real64), intent(out) :: dos(:)
 
    ! no PUSH_SUB, called too often
 
    select case (rdim)
    case (1)
      call simplex_dos_1d(esimplex, ef, dos)
    case (2)
      call simplex_dos_2d(esimplex, ef, dos)
    case (3)
      call simplex_dos_3d(esimplex, ef, dos)
    case default
      assert(.false.)
    end select
  end subroutine simplex_dos
 
  subroutine simplex_weights_1d(esegment, eF, weights, dos)
    real(real64), intent(in) :: esegment(:)
    real(real64), intent(in) :: ef
    real(real64), intent(out) :: weights(:)
    real(real64), intent(out) :: dos(:)
 
    real(real64) :: e(2), e1, e2
    real(real64) :: w(2), d(2)
    integer :: idx(2)
 
    real(real64), parameter :: vt_vg = 1.0_real64 
    real(real64), parameter :: vt_2vg = vt_vg / 2.0_real64 
 
    real(real64), parameter :: p(2,4) = 1.0_real64 / 60.0_real64 * reshape([ &
      64 , 1  , -3 , -2 , &
      1  , 64 , -2 , -3 ], shape(p), order=[2, 1]) 
 
    ! no PUSH_SUB, called too often
 
    select case (size(esegment))
    case (2)
      e(:) = esegment(:)
    case (4)
      e(:) = m_zero
      block
        integer :: i
        do i = 1, size(esegment)
          e(:) = e(:) + p(:,i) * esegment(i)
        end do
      end block
    case default
      assert(.false.)
    end select
 
    idx = [1,2]
    call sort(e, idx)
    e1 = e(1)
    e2 = e(2)
 
    if (ef <= e1) then
      w(:) = m_zero
      d(:) = m_zero
    elseif (e2 < ef) then
      w(:) = vt_2vg
      d(:) = m_zero
    elseif (e1 < ef .and. ef <= e2) then
      block
        real(real64) :: e21, c
        e21 = e2 - e1
        c = vt_2vg * (ef - e1) / e21
 
        w(:) = c * [ &
          2.0_real64 - (ef - e1) / e21, &
          (ef - e1) / e21]
 
        d(:) = vt_vg / e21 * [ &
          m_one - (ef - e1) / e21, &
          (ef - e1) / e21]
      end block
    else
      assert(.false.)
    end if
 
    dos(idx) = d
    select case (size(esegment))
    case (2)
      weights(idx) = w + sum(d) * (sum(e) - 2.0_real64 * e) / 12.0_real64
    case (4)
      weights(idx) = w ! no Bloechl correction
    case default
      assert(.false.)
    end select
  end subroutine simplex_weights_1d
 
  subroutine simplex_dos_1d(esegment, eF, dos)
    real(real64), intent(in) :: esegment(:)
    real(real64), intent(in) :: ef
    real(real64), intent(out) :: dos(:)
 
    real(real64) :: e(2), e1, e2
    real(real64) :: d(2)
    integer :: idx(2)
 
    real(real64), parameter :: vt_vg = 1.0_real64 
 
    real(real64), parameter :: p(2,4) = 1.0_real64 / 60.0_real64 * reshape([ &
      64 , 1  , -3 , -2 , &
      1  , 64 , -2 , -3 ], shape(p), order=[2, 1]) 
 
    ! no PUSH_SUB, called too often
 
    select case (size(esegment))
    case (2)
      e(:) = esegment(:)
    case (4)
      e(:) = m_zero
      block
        integer :: i
        do i = 1, size(esegment)
          e(:) = e(:) + p(:,i) * esegment(i)
        end do
      end block
    case default
      assert(.false.)
    end select
 
    idx = [1, 2]
    call sort(e, idx)
    e1 = e(1)
    e2 = e(2)
 
    if (ef <= e1 .or. e2 < ef) then
      d(:) = m_zero
    elseif (e1 < ef .and. ef <= e2) then
      block
        real(real64) :: e21
        e21 = e2 - e1
 
        d(:) = vt_vg / e21 * [ &
          m_one - (ef - e1) / e21, &
          (ef - e1) / e21]
      end block
    else
      assert(.false.)
    end if
 
    dos(idx) = d
  end subroutine simplex_dos_1d
 
  subroutine simplex_weights_2d(etriangle, eF, weights, dos)
    real(real64), intent(in) :: etriangle(:)
    real(real64), intent(in) :: ef
    real(real64), intent(out) :: weights(:)
    real(real64), intent(out) :: dos(:)
 
    real(real64) :: e(3), e1, e2, e3
    real(real64) :: w(3), d(3)
    integer :: idx(3)
 
    real(real64), parameter :: vt_vg = 1.0_real64 / 2.0_real64 
    real(real64), parameter :: vt_3vg = vt_vg / 3.0_real64 
 
    real(real64), parameter :: p(3,10) = 1.0_real64 / 360.0_real64 * reshape([ &
      402 , 0   , 6   , -13 , 5   , -17 , -13 , -11 , 7   , -6  , &
      6   , 396 , 6   , -9  , -15 , 3   , 3   , -15 , -9  , -6  , &
      6   , 0   , 402 , 7   , -11 , -13 , -17 , 5   , -13 , -6    &
      ], shape(p), order=[2, 1]) 
 
    ! no PUSH_SUB, called too often
 
    select case (size(etriangle))
    case (3)
      e(:) = etriangle(:)
    case (10)
      e(:) = m_zero
      block
        integer :: i
        do i = 1, size(etriangle)
          e(:) = e(:) + p(:,i) * etriangle(i)
        end do
      end block
    case default
      assert(.false.)
    end select
 
    idx = [1,2,3]
    call sort(e, idx)
    e1 = e(1)
    e2 = e(2)
    e3 = e(3)
 
    if (ef <= e1) then
      w(:) = m_zero
      d(:) = m_zero
    elseif (e3 < ef) then
      w(:) = vt_3vg
      d(:) = m_zero
    elseif (e1 < ef .and. ef <= e2) then
      block
        real(real64) :: e21, e31, c
        e21 = e2 - e1
        e31 = e3 - e1
        c = vt_3vg * (ef - e1) ** 2 / (e21 * e31)
 
        w(:) = c * [ &
          3.0_real64 - (ef - e1) * (m_one / e21 + m_one / e31), &
          (ef - e1) / e21, &
          (ef - e1) / e31]
 
        d(:) = vt_vg * (ef - e1) / (e21 * e31) * [&
          2.0_real64 - (ef - e1) * (m_one / e31 + m_one / e21), &
          (ef - e1) / e21, &
          (ef - e1) / e31]
      end block
    elseif (e2 < ef .and. ef <= e3) then
      block
        real(real64) :: e23, e31, c1, c2
        e23 = e2 - e3
        e31 = e3 - e1
        c1 = vt_3vg
        c2 = vt_3vg * (ef - e3) ** 2 / (e23 * e31)
 
        w(:) = [ &
          c1 - c2 * (ef - e3) / e31, &
          c1 + c2 * (ef - e3) / e23, &
          c1 + c2 * (3.0_real64 - (ef - e3) * (m_one / e23 - m_one / e31))]
 
        d(:) = vt_vg * (ef - e3) / (e23 * e31) * [ &
          - (ef - e3) / e31, &
          (ef - e3) / e23, &
          2.0_real64 - (ef - e3) * (m_one / e23 - m_one / e31)]
      end block
    else
      assert(.false.)
    end if
 
    dos(idx) = d
    select case (size(etriangle))
    case (3)
      weights(idx) = w + sum(d) * (sum(e) - 3.0_real64 * e) / 24.0_real64
    case (10)
      weights(idx) = w ! no Bloechl correction
    case default
      assert(.false.)
    end select
  end subroutine simplex_weights_2d
 
  subroutine simplex_dos_2d(etriangle, eF, dos)
    real(real64), intent(in) :: etriangle(:)
    real(real64), intent(in) :: ef
    real(real64), intent(out) :: dos(:)
 
    real(real64) :: e(3), e1, e2, e3
    real(real64) :: d(3)
    integer :: idx(3)
 
    real(real64), parameter :: vt_vg = 1.0_real64 / 2.0_real64 
 
    real(real64), parameter :: p(3,10) = 1.0_real64 / 360.0_real64 * reshape([ &
      402 , 0   , 6   , -13 , 5   , -17 , -13 , -11 , 7   , -6  , &
      6   , 396 , 6   , -9  , -15 , 3   , 3   , -15 , -9  , -6  , &
      6   , 0   , 402 , 7   , -11 , -13 , -17 , 5   , -13 , -6    &
      ], shape(p), order=[2, 1]) 
 
    ! no PUSH_SUB, called too often
 
    select case (size(etriangle))
    case (3)
      e(:) = etriangle(:)
    case (10)
      e(:) = m_zero
      block
        integer :: i
        do i = 1, size(etriangle)
          e(:) = e(:) + p(:,i) * etriangle(i)
        end do
      end block
    case default
      assert(.false.)
    end select
 
    idx = [1, 2, 3]
    call sort(e, idx)
    e1 = e(1)
    e2 = e(2)
    e3 = e(3)
 
    if (ef <= e1 .or. e3 < ef) then
      d(:) = m_zero
    elseif (e1 < ef .and. ef <= e2) then
      block
        real(real64) :: e21, e31
        e21 = e2 - e1
        e31 = e3 - e1
 
        d(:) = vt_vg * (ef - e1) / (e21 * e31) * [&
          (2.0_real64 - (ef - e1) * (m_one / e31 + m_one / e21)), &
          (ef - e1) / e21, &
          (ef - e1) / e31]
      end block
    elseif (e2 < ef .and. ef <= e3) then
      block
        real(real64) :: e23, e31
        e23 = e2 - e3
        e31 = e3 - e1
 
        d(:) = vt_vg * (ef - e3) / (e23 * e31) * [ &
          - (ef - e3) / e31, &
          (ef - e3) / e23, &
          2.0_real64 - (ef - e3) * (m_one / e23 - m_one / e31)]
      end block
    else
      assert(.false.)
    end if
 
    dos(idx) = d
  end subroutine simplex_dos_2d
 
  subroutine simplex_weights_3d(etetra, eF, weights, dos)
    real(real64), intent(in) :: etetra(:)
    real(real64), intent(in) :: ef
    real(real64), intent(out) :: weights(:)
    real(real64), intent(out) :: dos(:)
 
    real(real64) :: e(4), e1, e2, e3, e4
    real(real64) :: w(4), d(4)
    integer :: idx(4)
 
    real(real64), parameter :: vt_vg = 1.0_real64 / 6.0_real64 
    real(real64), parameter :: vt_4vg = vt_vg / 4.0_real64 
 
    real(real64), parameter :: p(4,20) = 1.0_real64 / 1260.0_real64 * reshape([ &
      1440, 0   , 30  , 0   , -38 , 7   , 17  , -28 , -56 , 9   , -46 , 9   , -38 , -28 , 17  , 7   , -18 , -18 , 12  , -18 , &
      0   , 1440, 0   , 30  , -28 , -38 , 7   , 17  , 9   , -56 , 9   , -46 , 7   , -38 , -28 , 17  , -18 , -18 , -18 , 12  , &
      30  , 0   , 1440, 0   , 17  , -28 , -38 , 7   , -46 , 9   , -56 , 9   , 17  , 7   , -38 , -28 , 12  , -18 , -18 , -18 , &
      0   , 30  , 0   , 1440, 7   , 17  , -28 , -38 , 9   , -46 , 9   , -56 , -28 , 17  , 7   , -38 , -18 , 12  , -18 , -18   &
      ], shape(p), order=[2, 1]) 
 
    ! no PUSH_SUB, called too often
 
    select case (size(etetra))
    case (4)
      e(:) = etetra(:)
    case (20)
      e(:) = m_zero
      block
        integer :: i
        do i = 1, size(etetra)
          e(:) = e(:) + p(:,i) * etetra(i)
        end do
      end block
    case default
      assert(.false.)
    end select
 
    idx = [1,2,3,4]
    call sort(e, idx)
    e1 = e(1)
    e2 = e(2)
    e3 = e(3)
    e4 = e(4)
 
    if (e1 >= ef) then
      w(:) = m_zero
      d(:) = m_zero
    elseif (e4 < ef) then
      w(:) = vt_4vg
      d(:) = m_zero
    elseif (e1 < ef .and. ef <= e2) then
      block
        real(real64) :: e21, e31, e41, c
        e21 = e2 - e1
        e31 = e3 - e1
        e41 = e4 - e1
        c = vt_4vg * (ef - e1) ** 3 / (e21 * e31 * e41)
 
        w(:) = c * [ &
          4.0_real64 - (ef - e1) * (m_one / e21 + m_one / e31 + m_one / e41), &
          (ef - e1) / e21, &
          (ef - e1) / e31, &
          (ef - e1) / e41]
      end block
      block
        real(real64) :: f12, f13, f14, f21, f31, f41, g
        f21 = (ef - e1) / (e2 - e1)
        f31 = (ef - e1) / (e3 - e1)
        f41 = (ef - e1) / (e4 - e1)
        f12 = m_one - f21
        f13 = m_one - f31
        f14 = m_one - f41
        g = f31 * f41 / (e2 - e1)
        d(:) = vt_vg * g * [&
          f12 + f13 + f14, &
          f21, &
          f31, &
          f41]
      end block
    elseif (e2 < ef .and. ef <= e3) then
      block
        real(real64) :: e21, e31, e32, e41, e42, c1, c2, c3
        e21 = e2 - e1
        e31 = e3 - e1
        e32 = e3 - e2
        e41 = e4 - e1
        e42 = e4 - e2
        c1 = vt_4vg * (ef - e1) ** 2 / (e41 * e31)
        c2 = vt_4vg * (ef - e1) * (ef - e2) * (e3 - ef) / (e41 * e32 * e31)
        c3 = vt_4vg * (ef - e2) ** 2 * (e4 - ef) / (e42 * e32 * e41)
 
        w(:) = [ &
          c1 + (c1 + c2) * (e3 - ef) / e31 + (c1 + c2 + c3) * (e4 - ef) / e41, &
          c1 + c2 + c3 + (c2 + c3) * (e3 - ef) / e32 + c3 * (e4 - ef) / e42, &
          (c1 + c2) * (ef - e1) / e31 + (c2 + c3) * (ef - e2) / e32, &
          (c1 + c2 + c3) * (ef - e1) / e41 + c3 * (ef - e2) / e42]
      end block
      block
        real(real64) :: f13, f14, f23, f24, f31, f32, f41, f42, g, delta
        delta = e4 - e1
        f31 = (ef - e1) / (e3 - e1)
        f41 = (ef - e1) / (e4 - e1)
        f32 = (ef - e2) / (e3 - e2)
        f42 = (ef - e2) / (e4 - e2)
        f13 = m_one - f31
        f14 = m_one - f41
        f23 = m_one - f32
        f24 = m_one - f42
        g = 3.0_real64 / delta * (f23 * f31 + f32 * f24)
        d(:) = vt_vg * [&
          g * f14 / 3.0_real64 + f13 * f31 * f23 / delta, &
          g * f23 / 3.0_real64 + f24 * f24 * f32 / delta, &
          g * f32 / 3.0_real64 + f31 * f31 * f23 / delta, &
          g * f41 / 3.0_real64 + f42 * f24 * f32 / delta]
      end block
    elseif (e3 < ef .and. ef <= e4) then
      block
        real(real64) :: e41, e42, e43, c
        e41 = e4 - e1
        e42 = e4 - e2
        e43 = e4 - e3
        c = vt_4vg * (e4 - ef) ** 3 / (e41 * e42 * e43)
 
        w(:) = vt_4vg - c * [ &
          (e4 - ef) / e41, &
          (e4 - ef) / e42, &
          (e4 - ef) / e43, &
          4.0_real64 - (e4 - ef) * (m_one / e41 + m_one / e42 + m_one / e43)]
      end block
      block
        real(real64) :: f14, f24, f34, f41, f42, f43, g
        f14 = (ef - e4) / (e1 - e4)
        f24 = (ef - e4) / (e2 - e4)
        f34 = (ef - e4) / (e3 - e4)
        f41 = m_one - f14
        f42 = m_one - f24
        f43 = m_one - f34
        g = f14 * f24 / (e4 - e3)
        d(:) = vt_vg * g *[ &
          f14, &
          f24, &
          f34, &
          f41 + f42 + f43]
      end block
    else
      assert(.false.)
    end if
 
    dos(idx) = d
    select case (size(etetra))
    case (4)
      weights(idx) = w + sum(d) * (sum(e) - 4.0_real64 * e) / 40.0_real64
    case (20)
      weights(idx) = w ! no Bloechl correction
    case default
      assert(.false.)
    end select
  end subroutine simplex_weights_3d
 
  subroutine simplex_dos_3d(etetra, eF, dos)
    real(real64), intent(in) :: etetra(:)
    real(real64), intent(in) :: ef
    real(real64), intent(out) :: dos(:)
 
    real(real64) :: e(4), e1, e2, e3, e4
    real(real64) :: d(4)
    integer :: idx(4)
 
    real(real64), parameter :: vt_vg = 1.0_real64 / 6.0_real64 
 
    real(real64), parameter :: p(4,20) = 1.0_real64 / 1260.0_real64 * reshape([ &
      1440, 0   , 30  , 0   , -38 , 7   , 17  , -28 , -56 , 9   , -46 , 9   , -38 , -28 , 17  , 7   , -18 , -18 , 12  , -18 , &
      0   , 1440, 0   , 30  , -28 , -38 , 7   , 17  , 9   , -56 , 9   , -46 , 7   , -38 , -28 , 17  , -18 , -18 , -18 , 12  , &
      30  , 0   , 1440, 0   , 17  , -28 , -38 , 7   , -46 , 9   , -56 , 9   , 17  , 7   , -38 , -28 , 12  , -18 , -18 , -18 , &
      0   , 30  , 0   , 1440, 7   , 17  , -28 , -38 , 9   , -46 , 9   , -56 , -28 , 17  , 7   , -38 , -18 , 12  , -18 , -18   &
      ], shape(p), order=[2, 1]) 
 
    ! no PUSH_SUB, called too often
 
    select case (size(etetra))
    case (4)
      e(:) = etetra(:)
    case (20)
      e(:) = m_zero
      block
        integer :: i
        do i = 1, size(etetra)
          e(:) = e(:) + p(:,i) * etetra(i)
        end do
      end block
    case default
      assert(.false.)
    end select
 
    idx = [1, 2, 3, 4]
    call sort(e, idx)
    e1 = e(1)
    e2 = e(2)
    e3 = e(3)
    e4 = e(4)
 
    if (e1 >= ef .or. e4 < ef) then
      d(:) = m_zero
    elseif (e1 < ef .and. ef <= e2) then
      block
        real(real64) :: f12, f13, f14, f21, f31, f41, g
        f21 = (ef - e1) / (e2 - e1)
        f31 = (ef - e1) / (e3 - e1)
        f41 = (ef - e1) / (e4 - e1)
        f12 = m_one - f21
        f13 = m_one - f31
        f14 = m_one - f41
        g = f31 * f41 / (e2 - e1)
        d(:) = vt_vg * g * [&
          f12 + f13 + f14, &
          f21, &
          f31, &
          f41]
      end block
    elseif (e2 < ef .and. ef <= e3) then
      block
        real(real64) :: f13, f14, f23, f24, f31, f32, f41, f42, g, delta
        delta = e4 - e1
        f31 = (ef - e1) / (e3 - e1)
        f41 = (ef - e1) / (e4 - e1)
        f32 = (ef - e2) / (e3 - e2)
        f42 = (ef - e2) / (e4 - e2)
        f13 = m_one - f31
        f14 = m_one - f41
        f23 = m_one - f32
        f24 = m_one - f42
        g = 3.0_real64 / delta * (f23 * f31 + f32 * f24)
        d(:) = vt_vg * [&
          g * f14 / 3.0_real64 + f13 * f31 * f23 / delta, &
          g * f23 / 3.0_real64 + f24 * f24 * f32 / delta, &
          g * f32 / 3.0_real64 + f31 * f31 * f23 / delta, &
          g * f41 / 3.0_real64 + f42 * f24 * f32 / delta]
      end block
    elseif (e3 < ef .and. ef <= e4) then
      block
        real(real64) :: f14, f24, f34, f41, f42, f43, g
        f14 = (ef - e4) / (e1 - e4)
        f24 = (ef - e4) / (e2 - e4)
        f34 = (ef - e4) / (e3 - e4)
        f41 = m_one - f14
        f42 = m_one - f24
        f43 = m_one - f34
        g = f14 * f24 / (e4 - e3)
        d(:) = vt_vg * g *[ &
          f14, &
          f24, &
          f34, &
          f41 + f42 + f43]
      end block
    else
      assert(.false.)
    end if
 
    dos(idx) = d
  end subroutine simplex_dos_3d
 
end module simplex_oct_m