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, debug_pop_sub, debug_push_sub
  use global_oct_m, only: &
    m_one, 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
  use profiling_oct_m, only: prof_vars, profiling_memory, profiling_memory_allocate, profiling_memory_deallocate
  implicit none
 
  private
 
  public :: &
    simplex_t, &
    simplex_init, &
    simplex_end, &
    simplex_weights, &
    simplex_dos
 
  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
 
  interface simplex_weights
    module procedure simplex_weights_single, simplex_weights_array
  end interface simplex_weights
 
  interface simplex_dos
    module procedure simplex_dos_single, simplex_dos_array
  end interface simplex_dos
 
contains
 
  pure subroutine simplex_sort_2(values, idx)
    real(real64), intent(inout) :: values(2)
    integer, intent(inout) :: idx(2)
 
    call simplex_compare_swap(values(1), values(2), idx(1), idx(2))
  end subroutine simplex_sort_2
 
  pure subroutine simplex_sort_3(values, idx)
    real(real64), intent(inout) :: values(3)
    integer, intent(inout) :: idx(3)
 
    call simplex_compare_swap(values(1), values(2), idx(1), idx(2))
    call simplex_compare_swap(values(2), values(3), idx(2), idx(3))
    call simplex_compare_swap(values(1), values(2), idx(1), idx(2))
  end subroutine simplex_sort_3
 
  pure subroutine simplex_sort_4(values, idx)
    real(real64), intent(inout) :: values(4)
    integer, intent(inout) :: idx(4)
 
    call simplex_compare_swap(values(1), values(2), idx(1), idx(2))
    call simplex_compare_swap(values(3), values(4), idx(3), idx(4))
    call simplex_compare_swap(values(1), values(3), idx(1), idx(3))
    call simplex_compare_swap(values(2), values(4), idx(2), idx(4))
    call simplex_compare_swap(values(2), values(3), idx(2), idx(3))
  end subroutine simplex_sort_4
 
  pure subroutine simplex_compare_swap(a, b, ia, ib)
    real(real64), intent(inout) :: a, b
    integer, intent(inout) :: ia, ib
 
    real(real64) :: tmp_a
    integer :: tmp_i
 
    if (a > b) then
      tmp_a = a
      a = b
      b = tmp_a
 
      tmp_i = ia
      ia = ib
      ib = tmp_i
    end if
  end subroutine simplex_compare_swap
 
  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_single(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(:)
 
    real(real64) :: weights_array(size(weights), 1), dos_array(size(dos), 1)
 
    ! no PUSH_SUB, called too often
 
    call simplex_weights_array(rdim, esimplex, [ef], weights_array, dos_array)
    weights(:) = weights_array(:, 1)
    dos(:) = dos_array(:, 1)
  end subroutine simplex_weights_single
 
  subroutine simplex_weights_array(rdim, esimplex, eFs, weights, dos)
    integer, intent(in) :: rdim
    real(real64), intent(in) :: esimplex(:)
    real(real64), intent(in) :: efs(:)
    real(real64), intent(out) :: weights(:,:)
    real(real64), intent(out) :: dos(:,:)
 
    ! no PUSH_SUB, called too often
 
    assert(size(weights, 1) == rdim + 1)
    assert(size(dos, 1) == rdim + 1)
    assert(size(weights, 2) == size(efs))
    assert(size(dos, 2) == size(efs))
 
    select case (rdim)
    case (1)
      call simplex_weights_1d(esimplex, efs, weights, dos)
    case (2)
      call simplex_weights_2d(esimplex, efs, weights, dos)
    case (3)
      call simplex_weights_3d(esimplex, efs, weights, dos)
    case default
      assert(.false.)
    end select
  end subroutine simplex_weights_array
 
  subroutine simplex_dos_single(rdim, esimplex, eF, dos)
    integer, intent(in) :: rdim
    real(real64), intent(in) :: esimplex(:)
    real(real64), intent(in) :: ef
    real(real64), intent(out) :: dos(:)
 
    real(real64) :: dos_array(size(dos), 1)
 
    ! no PUSH_SUB, called too often
 
    call simplex_dos_array(rdim, esimplex, [ef], dos_array)
    dos(:) = dos_array(:, 1)
  end subroutine simplex_dos_single
 
  subroutine simplex_dos_array(rdim, esimplex, eFs, dos)
    integer, intent(in) :: rdim
    real(real64), intent(in) :: esimplex(:)
    real(real64), intent(in) :: efs(:)
    real(real64), intent(out) :: dos(:,:)
 
    ! no PUSH_SUB, called too often
 
    assert(size(dos, 1) == rdim + 1)
    assert(size(dos, 2) == size(efs))
 
    select case (rdim)
    case (1)
      call simplex_dos_1d(esimplex, efs, dos)
    case (2)
      call simplex_dos_2d(esimplex, efs, dos)
    case (3)
      call simplex_dos_3d(esimplex, efs, dos)
    case default
      assert(.false.)
    end select
  end subroutine simplex_dos_array
 
  subroutine simplex_weights_1d(esegment, eFs, weights, dos)
    real(real64), intent(in) :: esegment(:)
    real(real64), intent(in) :: eFs(:)
    real(real64), intent(out) :: weights(:,:)
    real(real64), intent(out) :: dos(:,:)
 
    real(real64) :: E(2), E1, E2, eF
    real(real64) :: w(2), d(2), sumE, bloechl_corr(2)
    integer :: idx(2), ie, ne
    logical :: apply_bloechl
 
    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 simplex_sort_2(e, idx)
    e1 = e(1)
    e2 = e(2)
    ne = size(efs)
    apply_bloechl = (size(esegment) == 2)
    if (apply_bloechl) then
      sume = sum(e)
      bloechl_corr(:) = (sume - 2.0_real64 * e) / 12.0_real64
    end if
 
    do ie = 1, ne
      ef = efs(ie)
 
      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, ie) = d
      weights(idx, ie) = w
      if (apply_bloechl) weights(idx, ie) = weights(idx, ie) + sum(d) * bloechl_corr
    end do
  end subroutine simplex_weights_1d
 
  subroutine simplex_dos_1d(esegment, eFs, dos)
    real(real64), intent(in) :: esegment(:)
    real(real64), intent(in) :: eFs(:)
    real(real64), intent(out) :: dos(:,:)
 
    real(real64) :: E(2), E1, E2, eF
    real(real64) :: d(2)
    integer :: idx(2), ie, ne
 
    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 simplex_sort_2(e, idx)
    e1 = e(1)
    e2 = e(2)
    ne = size(efs)
 
    do ie = 1, ne
      ef = efs(ie)
 
      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, ie) = d
    end do
  end subroutine simplex_dos_1d
 
  subroutine simplex_weights_2d(etriangle, eFs, weights, dos)
    real(real64), intent(in) :: etriangle(:)
    real(real64), intent(in) :: eFs(:)
    real(real64), intent(out) :: weights(:,:)
    real(real64), intent(out) :: dos(:,:)
 
    real(real64) :: E(3), E1, E2, E3, eF
    real(real64) :: w(3), d(3), sumE, bloechl_corr(3)
    integer :: idx(3), ie, ne
    logical :: apply_bloechl
 
    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 simplex_sort_3(e, idx)
    e1 = e(1)
    e2 = e(2)
    e3 = e(3)
    ne = size(efs)
    apply_bloechl = (size(etriangle) == 3)
    if (apply_bloechl) then
      sume = sum(e)
      bloechl_corr(:) = (sume - 3.0_real64 * e) / 24.0_real64
    end if
 
    do ie = 1, ne
      ef = efs(ie)
 
      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, ie) = d
      weights(idx, ie) = w
      if (apply_bloechl) weights(idx, ie) = weights(idx, ie) + sum(d) * bloechl_corr
    end do
  end subroutine simplex_weights_2d
 
  subroutine simplex_dos_2d(etriangle, eFs, dos)
    real(real64), intent(in) :: etriangle(:)
    real(real64), intent(in) :: eFs(:)
    real(real64), intent(out) :: dos(:,:)
 
    real(real64) :: E(3), E1, E2, E3, eF
    real(real64) :: d(3)
    integer :: idx(3), ie, ne
 
    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 simplex_sort_3(e, idx)
    e1 = e(1)
    e2 = e(2)
    e3 = e(3)
    ne = size(efs)
 
    do ie = 1, ne
      ef = efs(ie)
 
      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, ie) = d
    end do
  end subroutine simplex_dos_2d
 
  subroutine simplex_weights_3d(etetra, eFs, weights, dos)
    real(real64), intent(in) :: etetra(:)
    real(real64), intent(in) :: eFs(:)
    real(real64), intent(out) :: weights(:,:)
    real(real64), intent(out) :: dos(:,:)
 
    real(real64) :: E(4), E1, E2, E3, E4, eF
    real(real64) :: w(4), d(4), sumE, bloechl_corr(4)
    integer :: idx(4), ie, ne
    logical :: apply_bloechl
 
    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 simplex_sort_4(e, idx)
    e1 = e(1)
    e2 = e(2)
    e3 = e(3)
    e4 = e(4)
    ne = size(efs)
    apply_bloechl = (size(etetra) == 4)
    if (apply_bloechl) then
      sume = sum(e)
      bloechl_corr(:) = (sume - 4.0_real64 * e) / 40.0_real64
    end if
 
    do ie = 1, ne
      ef = efs(ie)
 
      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, ie) = d
      weights(idx, ie) = w
      if (apply_bloechl) weights(idx, ie) = weights(idx, ie) + sum(d) * bloechl_corr
    end do
  end subroutine simplex_weights_3d
 
  subroutine simplex_dos_3d(etetra, eFs, dos)
    real(real64), intent(in) :: etetra(:)
    real(real64), intent(in) :: eFs(:)
    real(real64), intent(out) :: dos(:,:)
 
    real(real64) :: E(4), E1, E2, E3, E4, eF
    real(real64) :: d(4)
    integer :: idx(4), ie, ne
 
    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 simplex_sort_4(e, idx)
    e1 = e(1)
    e2 = e(2)
    e3 = e(3)
    e4 = e(4)
    ne = size(efs)
 
    do ie = 1, ne
      ef = efs(ie)
 
      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, ie) = d
    end do
  end subroutine simplex_dos_3d
 
end module simplex_oct_m