kpoints.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.
!!
 
#include "global.h"
 
module kpoints_oct_m
  use debug_oct_m
  use distributed_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lattice_vectors_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use simplex_oct_m, only: simplex_t
  use sort_oct_m
  use symm_op_oct_m
  use symmetries_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use utils_oct_m
 
  implicit none
 
  private
 
  public ::                       &
    kpoints_grid_t,               &
    kpoints_t,                    &
    kpoints_init,                 &
    kpoints_allocate_and_init,    &
    kpoints_end,                  &
    kpoints_copy,                 &
    kpoints_number,               &
    kpoints_point_is_gamma,       &
    kpoints_get_symmetry_ops,     &
    kpoints_get_num_symmetry_ops, &
    kpoints_get_full_symmetry_op_index, &
    kpoints_kweight_denominator,  &
    kpoints_grid_generate,        &
    kpoints_to_absolute,          &
    kpoints_get_path_coord,       &
    kpoints_grid_init,            &
    kpoints_path_generate,        &
    kpoints_to_reduced,           &
    kpoints_fold_to_1bz,          &
    kpoints_is_compatible_downsampling, &
    kpoints_is_valid_symmetry,    &
    kpoints_grid_end,             &
    kpoints_nkpt_in_path,         &
    kpoints_lattice_vectors_update, &
    kpoints_read_mp,              &
    kpoints_read_path,            &
    kpoints_read_user_kpoints,    &
    kpoints_deterministic_sort
 
  type kpoints_grid_t
    ! Components are public by default
    real(real64), allocatable :: point(:, :)
    real(real64), allocatable :: point1BZ(:, :)
    real(real64), allocatable :: red_point(:, :)
    real(real64), allocatable :: weight(:)
    type(simplex_t),  pointer :: simplex => null()
    integer            :: nshifts = 1
    real(real64), allocatable :: shifts(:, :)
    integer            :: npoints = 0
    integer            :: dim = 0
  end type kpoints_grid_t
 
  type kpoints_t
    ! Components are public by default
 
    type(lattice_vectors_t) :: latt
 
    type(kpoints_grid_t) :: full
    type(kpoints_grid_t) :: reduced
 
    integer              :: method = 0
 
    logical              :: use_symmetries = .false.
    logical              :: use_time_reversal = .false.
    integer              :: nik_skip = 0
 
    integer, allocatable          :: nik_axis(:)
    integer, allocatable          :: niq_axis(:)
    integer, allocatable          :: lk123(:, :)
    integer, allocatable, private :: symmetry_ops(:, :)
    integer, allocatable, private :: num_symmetry_ops(:)
    integer, allocatable, private :: map_full_to_reduced(:)
    integer, allocatable, private :: map_full_to_symidx(:)
 
    real(real64), allocatable, private   :: coord_along_path(:)
 
    integer, allocatable          :: downsampling(:)
 
    type(symmetries_t), pointer :: symm => null()
 
  contains
    procedure :: have_zero_weight_path => kpoints_have_zero_weight_path
    procedure :: get_kpoint_method => kpoints_get_kpoint_method
    procedure :: gamma_only => kpoints_gamma_only
    procedure :: get_weight => kpoints_get_weight
    procedure :: get_equiv => kpoints_get_equiv
    procedure :: get_full_symmetry_op_index => kpoints_get_full_symmetry_op_index
    procedure :: get_point =>  kpoints_get_point
    procedure :: write_info =>  kpoints_write_info
    procedure :: nkpt_in_path => kpoints_nkpt_in_path
  end type kpoints_t
 
  integer, public, parameter ::        &
    KPOINTS_GAMMA       =  2,          &
    kpoints_monkh_pack  =  4,          &
    kpoints_user        =  8,          &
    kpoints_path        =  16
 
contains
 
  ! ---------------------------------------------------------
  subroutine kpoints_grid_init(dim, this, npoints, nshifts)
    integer,              intent(in)  :: dim
    type(kpoints_grid_t), intent(out) :: this
    integer,              intent(in)  :: npoints
    integer,              intent(in)  :: nshifts
 
    push_sub(kpoints_grid_init)
 
    this%dim = dim
    this%npoints = npoints
    this%nshifts = nshifts
    safe_allocate(this%red_point(1:dim, 1:npoints))
    safe_allocate(this%point(1:dim, 1:npoints))
    safe_allocate(this%point1bz(1:dim,1:npoints))
    safe_allocate(this%weight(1:npoints))
    safe_allocate(this%shifts(1:dim,1:nshifts))
 
    pop_sub(kpoints_grid_init)
  end subroutine kpoints_grid_init
 
  ! ---------------------------------------------------------
  subroutine kpoints_grid_end(this)
    type(kpoints_grid_t), intent(inout) :: this
 
    push_sub(kpoints_grid_end)
 
    safe_deallocate_a(this%red_point)
    safe_deallocate_a(this%point)
    safe_deallocate_a(this%point1BZ)
    safe_deallocate_a(this%weight)
    safe_deallocate_a(this%shifts)
    safe_deallocate_p(this%simplex)
 
    pop_sub(kpoints_grid_end)
  end subroutine kpoints_grid_end
 
  ! ---------------------------------------------------------
  subroutine kpoints_grid_copy(bb, aa)
    type(kpoints_grid_t), intent(in)    :: bb
    type(kpoints_grid_t), intent(inout) :: aa
 
    push_sub(kpoints_grid_copy)
 
    call kpoints_grid_end(aa)
    call kpoints_grid_init(bb%dim, aa, bb%npoints, bb%nshifts)
    aa%weight(1:bb%npoints) = bb%weight(1:bb%npoints)
    aa%point(1:bb%dim, 1:bb%npoints)  = bb%point(1:bb%dim, 1:bb%npoints)
    aa%point1BZ(1:bb%dim, 1:bb%npoints) = bb%point1BZ(1:bb%dim, 1:bb%npoints)
    aa%red_point(1:bb%dim, 1:bb%npoints) = bb%red_point(1:bb%dim, 1:bb%npoints)
    aa%shifts(1:bb%dim, 1:bb%nshifts) = bb%shifts(1:bb%dim, 1:bb%nshifts)
 
    pop_sub(kpoints_grid_copy)
  end subroutine kpoints_grid_copy
 
  ! ---------------------------------------------------------
  subroutine kpoints_grid_addto(this, that)
    type(kpoints_grid_t), intent(inout) :: this
    type(kpoints_grid_t), intent(in)    :: that
 
    type(kpoints_grid_t) :: old_grid
 
    push_sub(kpoints_grid_addto)
 
    if (.not. allocated(that%point)) then
      pop_sub(kpoints_grid_addto)
      return
    end if
 
    if (.not. allocated(this%point)) then
      call kpoints_grid_copy(that, this)
      pop_sub(kpoints_grid_addto)
      return
    end if
 
    call kpoints_grid_copy(this, old_grid)
 
    call kpoints_grid_end(this)
    call kpoints_grid_init(old_grid%dim, this, that%npoints + old_grid%npoints, old_grid%nshifts)
 
    this%red_point = m_zero
    this%point = m_zero
    this%weight = m_zero
    this%shifts = m_zero
 
    ! Fill the the result with values form this
    this%red_point(1:old_grid%dim, 1:old_grid%npoints)= old_grid%red_point(1:old_grid%dim, 1:old_grid%npoints)
    this%point(1:old_grid%dim, 1:old_grid%npoints)    = old_grid%point(1:old_grid%dim, 1:old_grid%npoints)
    this%point1BZ(1:old_grid%dim, 1:old_grid%npoints) = old_grid%point1BZ(1:old_grid%dim, 1:old_grid%npoints)
    this%weight(1:old_grid%npoints)                   = old_grid%weight(1:old_grid%npoints)
    this%shifts(1:old_grid%dim, 1:old_grid%nshifts)   = old_grid%shifts(1:old_grid%dim, 1:old_grid%nshifts)
 
    ! Fill the result with that
    this%red_point(1:old_grid%dim, old_grid%npoints+1:this%npoints)= that%red_point(1:that%dim, 1:that%npoints)
    this%point(1:old_grid%dim, old_grid%npoints+1:this%npoints)    = that%point(1:that%dim, 1:that%npoints)
    this%point1BZ(1:old_grid%dim, old_grid%npoints+1:this%npoints) = that%point1BZ(1:that%dim, 1:that%npoints)
    this%weight(old_grid%npoints+1:this%npoints)                   = that%weight(1:that%npoints)
 
    call kpoints_grid_end(old_grid)
 
 
    pop_sub(kpoints_grid_addto)
  end subroutine kpoints_grid_addto
 
 
  function kpoints_set_method(namespace, use_symmetries, only_gamma) result(method)
    type(namespace_t),  intent(in)   :: namespace
    logical,            intent(in)   :: use_symmetries
    logical,            intent(in)   :: only_gamma
 
    integer :: method
 
    push_sub(kpoints_set_method)
 
    method = 0
 
    if (only_gamma) then
      method = kpoints_gamma
      pop_sub(kpoints_set_method)
      return
    end if
 
    !Monkhorst Pack grid
    if (parse_is_defined(namespace, 'KPointsGrid')) then
      method = ior(method, kpoints_monkh_pack)
    end if
 
    !User-defined kpoints
    if (parse_is_defined(namespace, 'KPointsReduced').or. parse_is_defined(namespace, 'KPoints')) then
      method = ior(method, kpoints_user)
      if (use_symmetries) then
        write(message(1), '(a)') "User-defined k-points are not compatible with KPointsUseSymmetries=yes."
        call messages_warning(1, namespace=namespace)
      end if
    end if
 
    !User-defined k-points path
    if (parse_is_defined(namespace, 'KPointsPath')) then
      method = ior(method, kpoints_path)
      if (use_symmetries) then
        write(message(1), '(a)') "KPointsPath is not compatible with KPointsUseSymmetries=yes."
        call messages_warning(1, namespace=namespace)
      end if
    end if
 
    if (method == 0) then
      write(message(1), '(a)') "Unable to determine the method for defining k-points."
      write(message(2), '(a)') "Octopus will continue assuming a Monkhorst Pack grid."
      call messages_warning(2, namespace=namespace)
      method = kpoints_monkh_pack
    end if
 
    pop_sub(kpoints_set_method)
 
  end function kpoints_set_method
 
  subroutine kpoints_allocate_and_init(this, symm, dim, latt)
    type(kpoints_t),            intent(inout) :: this
    type(symmetries_t), target, intent(in)    :: symm
    integer,                    intent(in)    :: dim
    type(lattice_vectors_t),    intent(in)    :: latt
 
    push_sub(kpoints_allocate_and_init)
 
    this%latt = latt
    this%symm => symm
 
    safe_allocate(this%nik_axis(1:dim))
    safe_allocate(this%niq_axis(1:dim))
    safe_allocate(this%downsampling(1:dim))
 
    this%nik_axis = 1
    this%niq_axis = 1
    this%downsampling = 1
    this%nik_skip = 0
    safe_deallocate_a(this%map_full_to_reduced)
    safe_deallocate_a(this%map_full_to_symidx)
 
    pop_sub(kpoints_allocate_and_init)
  end subroutine kpoints_allocate_and_init
 
  ! ---------------------------------------------------------
  subroutine kpoints_init(this, namespace, symm, dim, periodic_dim, latt)
    type(kpoints_t),            intent(out) :: this
    type(namespace_t),          intent(in)  :: namespace
    type(symmetries_t), target, intent(in)  :: symm
    integer,                    intent(in)  :: dim
    integer,                    intent(in)  :: periodic_dim
    type(lattice_vectors_t),    intent(in)  :: latt
 
    logical :: only_gamma, default_timereversal
 
    push_sub(kpoints_init)
 
    !%Variable KPointsUseSymmetries
    !%Type logical
    !%Default no
    !%Section Mesh::KPoints
    !%Description
    !% This variable defines whether symmetries are taken into account
    !% or not for the choice of <i>k</i>-points. If it is set to no, the <i>k</i>-point
    !% sampling will range over the full Brillouin zone.
    !%
    !% When a perturbation is applied to the system, the full
    !% symmetries of the system cannot be used. In this case you must
    !% not use symmetries or use the <tt>SymmetryBreakDir</tt> to tell
    !% Octopus the direction of the perturbation (for the moment this
    !% has to be done by hand by the user, in the future it will be
    !% automatic).
    !%
    !%End
    call parse_variable(namespace, 'KPointsUseSymmetries', .false., this%use_symmetries)
 
    !%Variable KPointsUseTimeReversal
    !%Type logical
    !%Section Mesh::KPoints
    !%Description
    !% If symmetries are used to reduce the number of <i>k</i>-points,
    !% this variable defines whether time-reversal symmetry is taken
    !% into account or not. If it is set to no, the <i>k</i>-point
    !% sampling will not be reduced according to time-reversal
    !% symmetry.
    !%
    !% The default is yes, unless symmetries are broken in one
    !% direction by the <tt>SymmetryBreakDir</tt> block.
    !%
    !% Warning: For time propagation runs with an external field,
    !% time-reversal symmetry should not be used.
    !%
    !%End
    default_timereversal = this%use_symmetries .and. .not. symmetries_have_break_dir(symm)
    call parse_variable(namespace, 'KPointsUseTimeReversal', default_timereversal, this%use_time_reversal)
 
    only_gamma = (periodic_dim == 0)
    this%method = kpoints_set_method(namespace, this%use_symmetries, only_gamma)
    call kpoints_allocate_and_init(this, symm, dim, latt)
 
    if (only_gamma) then
      call kpoints_read_mp(this, namespace, symm, dim, gamma_only=.true.)
      pop_sub(kpoints_init)
      return
    end if
 
    if (bitand(this%method, kpoints_monkh_pack) /= 0) then
      call kpoints_read_mp(this, namespace, symm, dim, gamma_only=.false.)
    end if
 
    if (bitand(this%method, kpoints_user) /= 0) then
      call kpoints_read_user_kpoints(this, namespace, dim)
    end if
 
    if (bitand(this%method, kpoints_path) /= 0) then
      call kpoints_read_path(this, namespace, dim)
    end if
 
    call kpoints_print_info(this, namespace)
 
    pop_sub(kpoints_init)
  end subroutine kpoints_init
 
  subroutine kpoints_read_mp(this, namespace, symm, dim, gamma_only)
    type(kpoints_t),    intent(inout) :: this
    type(namespace_t),  intent(in)    :: namespace
    type(symmetries_t), intent(in)    :: symm
    integer,            intent(in)    :: dim
    logical,            intent(in)    :: gamma_only
 
    logical       :: gamma_only_
    integer       :: ii, ik, is, ncols, nshifts
    type(block_t) :: blk
    integer,      allocatable :: symm_ops(:, :), num_symm_ops(:)
    real(real64), allocatable :: shifts(:, :)
 
    push_sub(kpoints_read_mp)
 
    call messages_obsolete_variable(namespace, 'KPointsMonkhorstPack', 'KPointsGrid')
 
    !%Variable KPointsGrid
    !%Type block
    !%Default <math>\Gamma</math>-point only
    !%Section Mesh::KPoints
    !%Description
    !% When this block is given (and the <tt>KPoints</tt> block is not present),
    !% <i>k</i>-points are distributed in a uniform grid, according to a modified
    !% version of the Monkhorst-Pack scheme. For the original MP scheme, see
    !% James D. Pack and Hendrik J. Monkhorst,
    !% <i>Phys. Rev. B</i> <b>13</b>, 5188 (1976) and <i>Phys. Rev. B</i> <b>16</b>, 1748 (1977).
    !%
    !% The number of columns should be equal to <tt>Dimensions</tt>,
    !% but the grid and shift numbers should be 1 and zero in finite directions.
    !%
    !% The first row of the block is a set of integers defining
    !% the number of <i>k</i>-points to be used along each direction
    !% in reciprocal space. The numbers refer to the whole Brillouin
    !% zone, and the actual number of <i>k</i>-points is usually
    !% reduced exploiting the symmetries of the system.  By default
    !% the grid will always include the <math>\Gamma</math>-point.
    !%
    !% Optional rows can be added to specify multiple shifts in the <i>k</i>-points (between 0.0 and 1.0),
    !% in units of the Brillouin zone divided by the number in the first row.
    !% Please note that without specifying any shift, an implicit shift of -0.5 is added for odd number
    !% of k-points, such that the Gamma point is always included.
    !%
    !% For example, the following input samples the BZ with 100 points in the
    !% <i>xy</i>-plane of reciprocal space:
    !%
    !% <tt>%KPointsGrid
    !% <br>&nbsp;&nbsp;10 | 10 | 1
    !% <br>%</tt>
    !%
    !% is equivalent to
    !%
    !% <tt>%KPointsGrid
    !% <br> 10 | 10 | 1
    !% <br> 0 | 0 | -0.5
    !% <br>%</tt>
    !%
    !%End
 
    gamma_only_ = gamma_only
    if (.not. gamma_only_) then
      gamma_only_ = (parse_block(namespace, 'KPointsGrid', blk) /= 0)
    end if
 
    this%nik_axis = 1
 
    if (.not. gamma_only_) then
      nshifts = max(parse_block_n(blk)-1,1)
    else
      nshifts = 1
    end if
    safe_allocate(shifts(1:dim,1:nshifts))
    shifts = m_zero
 
    if (.not. gamma_only_) then
      ncols = parse_block_cols(blk, 0)
      if (ncols /= dim) then
        write(message(1),'(a,i3,a,i3)') 'KPointsGrid first row has ', ncols, ' columns but should have ', dim
        if (ncols < dim) then
          call messages_fatal(1, namespace=namespace)
        else
          write(message(2),'(a)') 'Continuing, but ignoring the additional values.'
          call messages_warning(2, namespace=namespace)
        end if
      end if
      do ii = 1, dim
        call parse_block_integer(blk, 0, ii - 1, this%nik_axis(ii))
      end do
 
      if (any(this%nik_axis < 1)) then
        message(1) = 'Input: KPointsGrid is not valid.'
        call messages_fatal(1, namespace=namespace)
      end if
 
      if (parse_block_n(blk) > 1) then ! we have a shift, or even more
        ncols = parse_block_cols(blk, 1)
        if (ncols /= dim) then
          write(message(1),'(a,i3,a,i3)') 'KPointsGrid shift has ', ncols, ' columns but must have ', dim
          call messages_fatal(1, namespace=namespace)
        end if
        do is = 1, nshifts
          do ii = 1, dim
            call parse_block_float(blk, is, ii - 1, shifts(ii,is))
          end do
        end do
      else
        !We include a default 0.5 shift for odd number of k-points
        do ii = 1, dim
          if (mod(this%nik_axis(ii), 2) /= 0) then
            shifts(ii,1) = m_half
          end if
        end do
      end if
 
      call parse_block_end(blk)
    else
      shifts(:, 1) = -m_half
    end if
 
    !%Variable QPointsGrid
    !%Type block
    !%Default KPointsGrid
    !%Section Mesh::KPoints
    !%Description
    !% This block allows to define a q-point grid used for the calculation of the Fock operator
    !% with k-points. The <i>q</i>-points are distributed in a uniform grid, as done for the
    !% <tt>KPointsGrid</tt> variable.
    !% See J. Chem Phys. 124, 154709 (2006) for details
    !%
    !% For each dimension, the number of q point must be a divider of the number of  k point
    !%
    !% <tt>%QPointsGrid
    !% <br>&nbsp;&nbsp;2 | 2 | 1
    !% <br>%</tt>
    !%
    !% At the moment, this is not compatible with k-point symmetries.
    !%
    !%End
    this%niq_axis(:) = this%nik_axis(:)
    this%downsampling(:) = 1
    if (parse_is_defined(namespace, 'QPointsGrid')) then
      if (parse_block(namespace, 'QPointsGrid', blk) == 0) then
        ncols = parse_block_cols(blk, 0)
        if (ncols /= dim) then
          write(message(1),'(a,i3,a,i3)') 'QPointsGrid first row has ', ncols, ' columns but must have ', dim
          call messages_fatal(1, namespace=namespace)
        end if
        do ii = 1, dim
          call parse_block_integer(blk, 0, ii - 1, this%niq_axis(ii))
        end do
 
        if (any(this%nik_axis/this%niq_axis /= nint(this%nik_axis/real(this%niq_axis)))) then
          message(1) = 'Input: QPointsGrid is not compatible with the KPointsGrid.'
          call messages_fatal(1, namespace=namespace)
        end if
 
        this%downsampling = this%nik_axis/this%niq_axis
 
        if (any(this%downsampling /= 1) .and. this%use_symmetries) then
          call messages_not_implemented('QPointsGrid together with k-point symmetries', namespace=namespace)
        end if
 
        call parse_block_end(blk)
      end if
    end if
 
    call kpoints_grid_init(dim, this%full, product(this%nik_axis)*nshifts, nshifts)
 
    !We move the k-points into this%shifts
    this%full%shifts(:, :) = shifts(:, :)
    safe_deallocate_a(shifts)
 
    ! The grid labels are kept, as they give the position of a k-point on the
    ! Monkhorst-Pack grid without having to recover it from its coordinates.
    safe_allocate(this%lk123(1:this%full%npoints, 1:dim))
    call kpoints_grid_generate(dim, this%nik_axis, this%full%nshifts, this%full%shifts, this%full%red_point, &
      lk123 = this%lk123)
 
    do ik = 1, this%full%npoints
      call kpoints_to_absolute(this%latt, this%full%red_point(:, ik), this%full%point(:, ik))
    end do
 
    this%full%weight = m_one / this%full%npoints
 
    if (this%use_symmetries) then
      message(1) = "Checking if the generated full k-point grid is symmetric"
      call messages_info(1, namespace=namespace)
      call kpoints_check_symmetries(this%full, symm, dim, this%use_time_reversal, namespace)
    end if
 
    call kpoints_grid_copy(this%full, this%reduced)
 
    if (this%use_symmetries) then
      safe_allocate(num_symm_ops(1:this%full%npoints))
 
      if (this%use_time_reversal) then
        safe_allocate(symm_ops(1:this%full%npoints, 1:2*symmetries_number(symm)))
      else
        safe_allocate(symm_ops(1:this%full%npoints, 1:symmetries_number(symm)))
      end if
 
      safe_allocate(this%map_full_to_reduced(1:this%full%npoints))
      safe_allocate(this%map_full_to_symidx(1:this%full%npoints))
      call kpoints_grid_reduce(symm, this%use_time_reversal, &
        this%reduced%npoints, dim, this%reduced%red_point, this%reduced%weight, this%map_full_to_reduced, &
        this%map_full_to_symidx, symm_ops, num_symm_ops)
 
      assert(maxval(num_symm_ops) >= 1)
      if (this%use_time_reversal) then
        assert(maxval(num_symm_ops) <= 2*symmetries_number(symm))
      else
        assert(maxval(num_symm_ops) <= symmetries_number(symm))
      end if
 
      ! the total number of symmetry operations in the list has to be equal to the number of k-points
      assert(sum(num_symm_ops(1:this%reduced%npoints)) == this%full%npoints)
 
      do ik = 1, this%reduced%npoints
        assert(all(symm_ops(ik, 1:num_symm_ops(ik)) <= symm%nops))
      end do
 
      safe_allocate(this%symmetry_ops(1:this%reduced%npoints, 1:maxval(num_symm_ops)))
      safe_allocate(this%num_symmetry_ops(1:this%reduced%npoints))
 
      this%num_symmetry_ops(1:this%reduced%npoints) = num_symm_ops(1:this%reduced%npoints)
      this%symmetry_ops(1:this%reduced%npoints, 1:maxval(num_symm_ops)) = &
        symm_ops(1:this%reduced%npoints, 1:maxval(num_symm_ops))
 
      safe_deallocate_a(num_symm_ops)
      safe_deallocate_a(symm_ops)
    else
      safe_allocate(this%num_symmetry_ops(1:this%reduced%npoints))
      safe_allocate(this%symmetry_ops(1:this%reduced%npoints, 1))
      this%num_symmetry_ops(1:this%reduced%npoints) = 1
      this%symmetry_ops(1:this%reduced%npoints, 1) = 1
      safe_allocate(this%map_full_to_reduced(1:this%full%npoints))
      do ik = 1, this%full%npoints
        this%map_full_to_reduced(ik) = ik
      end do
      safe_allocate(this%map_full_to_symidx(1:this%full%npoints))
      this%map_full_to_symidx(:) = 1
    end if
 
    do ik = 1, this%reduced%npoints
      call kpoints_to_absolute(this%latt, this%reduced%red_point(:, ik), this%reduced%point(:, ik))
    end do
 
    call kpoints_fold_to_1bz(this%full, this%latt)
    call kpoints_fold_to_1bz(this%reduced, this%latt)
 
    block
      integer(int64) :: dos_method, smearing_function
      logical :: use_simplex, use_simplex_opt
 
      call parse_variable(namespace, 'DOSMethod', option__dosmethod__smear, dos_method)
      call parse_variable(namespace, 'SmearingFunction', option__smearingfunction__semiconducting, smearing_function)
 
      use_simplex = dos_method == option__dosmethod__tetrahedra .or. &
        smearing_function == option__smearingfunction__tetrahedra
      use_simplex_opt = dos_method == option__dosmethod__tetrahedra_opt .or. &
        smearing_function == option__smearingfunction__tetrahedra_opt
 
      if (use_simplex .and. use_simplex_opt) then
        call messages_not_implemented("DOSMethod and SmearingFunction with different tetrahedron methods", namespace=namespace)
      end if
 
      if (use_simplex .or. use_simplex_opt) then
        safe_deallocate_p(this%full%simplex)
        safe_deallocate_p(this%reduced%simplex)
      end if
 
      if (use_simplex) then
        this%full%simplex => simplex_t(dim, this%nik_axis, this%full%nshifts, this%full%shifts, this%full%red_point, &
          opt = .false.)
        this%reduced%simplex => simplex_t(dim, this%nik_axis, this%full%nshifts, this%full%shifts, this%full%red_point, &
          equiv = this%map_full_to_reduced, opt = .false.)
      elseif (use_simplex_opt) then
        this%full%simplex => simplex_t(dim, this%nik_axis, this%full%nshifts, this%full%shifts, this%full%red_point, &
          opt = .true.)
        this%reduced%simplex => simplex_t(dim, this%nik_axis, this%full%nshifts, this%full%shifts, this%full%red_point, &
          equiv = this%map_full_to_reduced, opt = .true.)
      end if
    end block
 
    pop_sub(kpoints_read_mp)
  end subroutine kpoints_read_mp
 
 
  subroutine kpoints_read_path(this, namespace, dim)
    type(kpoints_t),   intent(inout) :: this
    type(namespace_t), intent(in)    :: namespace
    integer,           intent(in)    :: dim
 
    type(block_t) :: blk
    integer :: nshifts, nkpoints, nhighsympoints, nsegments
    integer :: icol, ik, idir, ncols, npoints
    integer, allocatable :: resolution(:)
    real(real64), allocatable   :: highsympoints(:, :)
    type(kpoints_grid_t) :: path_kpoints_grid
    integer, allocatable :: symm_ops(:, :), num_symm_ops(:)
 
    push_sub(kpoints_read_path)
 
    !%Variable KPointsPath
    !%Type block
    !%Section Mesh::KPoints
    !%Description
    !% When this block is given, <i>k</i>-points are generated along a path
    !% defined by the points of the list.
    !% The points must be given in reduced coordinates.
    !%
    !% The first row of the block is a set of integers defining
    !% the number of <i>k</i>-points for each segments of the path.
    !% The number of columns should be equal to <tt>Dimensions</tt>,
    !% and the k-points coordinate should be zero in finite directions.
    !%
    !% For example, the following input samples the BZ with 15 points:
    !%
    !% <tt>%KPointsPath
    !% <br>&nbsp;&nbsp;10 | 5
    !% <br>&nbsp;&nbsp; 0 | 0 | 0
    !% <br>&nbsp;&nbsp; 0.5 | 0 | 0
    !% <br>&nbsp;&nbsp; 0.5 | 0.5 | 0.5
    !% <br>%</tt>
    !%
    !%End
 
    if (parse_block(namespace, 'KPointsPath', blk) /= 0) then
      write(message(1),'(a)') 'Internal error while reading KPointsPath.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    ! There is one high symmetry k-point per line
    nsegments = parse_block_cols(blk, 0)
    nhighsympoints = parse_block_n(blk) - 1
    if (nhighsympoints /= nsegments + 1) then
      write(message(1),'(a,i3,a,i3)') 'The first row of KPointsPath is not compatible with the number of specified k-points.'
      call messages_fatal(1, namespace=namespace)
    end if
 
    safe_allocate(resolution(1:nsegments))
    do icol = 1, nsegments
      call parse_block_integer(blk, 0, icol-1, resolution(icol))
    end do
    !Total number of points in the segment
    nkpoints = sum(resolution) + 1
 
    safe_allocate(highsympoints(1:dim, 1:nhighsympoints))
    do ik = 1, nhighsympoints
      ncols = parse_block_cols(blk, ik)
      if (ncols /= dim) then
        write(message(1),'(a,i8,a,i3,a,i3)') 'KPointsPath row ', ik, ' has ', ncols, ' columns but must have ', dim
        call messages_fatal(1, namespace=namespace)
      end if
 
      do idir = 1, dim
        call parse_block_float(blk, ik, idir-1, highsympoints(idir, ik))
      end do
    end do
 
    call parse_block_end(blk)
 
    !We do not have shifts
    nshifts = 1
    call kpoints_grid_init(dim, path_kpoints_grid, nkpoints, nshifts)
 
    ! For the output of band-structures
    safe_allocate(this%coord_along_path(1:nkpoints))
 
    call kpoints_path_generate(dim, this%latt, nkpoints, nsegments, resolution, highsympoints, path_kpoints_grid%point, &
      this%coord_along_path)
 
    safe_deallocate_a(resolution)
    safe_deallocate_a(highsympoints)
 
    !Use zero weight for the path so that it can be used with any kind of calculation mode
    !without interfering with the physical BZ integral.
    if (this%method == kpoints_path) then
      path_kpoints_grid%weight = m_one/path_kpoints_grid%npoints
    else
      path_kpoints_grid%weight = m_zero
      this%nik_skip = this%nik_skip + path_kpoints_grid%npoints
    end if
 
    !The points have been generated in absolute coordinates
    do ik = 1, path_kpoints_grid%npoints
      call kpoints_to_reduced(this%latt, path_kpoints_grid%point(:, ik), path_kpoints_grid%red_point(:, ik))
    end do
 
    call kpoints_fold_to_1bz(path_kpoints_grid, this%latt)
 
    !We need to copy the arrays containing the information on the symmetries
    !Before calling kpoints_grid_addto
    if (allocated(this%symmetry_ops)) then
      npoints = this%reduced%npoints
      safe_allocate(num_symm_ops(1:npoints))
      safe_allocate(symm_ops(1:npoints, 1:maxval(this%num_symmetry_ops)))
 
      num_symm_ops(1:npoints) = this%num_symmetry_ops(1:npoints)
      symm_ops(1:npoints, 1:maxval(num_symm_ops)) = &
        this%symmetry_ops(1:npoints, 1:maxval(num_symm_ops))
    end if
 
    call kpoints_grid_addto(this%full   , path_kpoints_grid)
    call kpoints_grid_addto(this%reduced, path_kpoints_grid)
 
    if (allocated(this%symmetry_ops)) then
      safe_deallocate_a(this%num_symmetry_ops)
      safe_deallocate_a(this%symmetry_ops)
      safe_allocate(this%num_symmetry_ops(1:this%reduced%npoints))
      safe_allocate(this%symmetry_ops(1:this%reduced%npoints,1:maxval(num_symm_ops)))
 
      this%num_symmetry_ops(1:npoints) = num_symm_ops(1:npoints)
      this%symmetry_ops(1:npoints, 1:maxval(num_symm_ops)) = &
        symm_ops(1:npoints, 1:maxval(num_symm_ops))
      this%num_symmetry_ops(npoints+1:this%reduced%npoints) = 1
      this%symmetry_ops(npoints+1:this%reduced%npoints, 1) = 1
 
      safe_deallocate_a(num_symm_ops)
      safe_deallocate_a(symm_ops)
    else if(this%use_symmetries) then ! If symmetries are activated but no KPointsGrid specified
      safe_allocate(this%num_symmetry_ops(1:this%reduced%npoints))
      safe_allocate(this%symmetry_ops(1:this%reduced%npoints,1:1))
 
      this%num_symmetry_ops(1:this%reduced%npoints) = 1
      this%symmetry_ops(1:this%reduced%npoints, 1) = 1
    end if
 
    call kpoints_grid_end(path_kpoints_grid)
 
    pop_sub(kpoints_read_path)
  end subroutine kpoints_read_path
 
 
  subroutine kpoints_read_user_kpoints(this, namespace, dim)
    type(kpoints_t),   intent(inout) :: this
    type(namespace_t), intent(in)    :: namespace
    integer,           intent(in)    :: dim
 
    type(block_t) :: blk
    logical :: reduced
    integer :: ik, idir
    real(real64) :: weight_sum
    type(kpoints_grid_t) :: user_kpoints_grid
 
    push_sub(kpoints_read_user_kpoints)
 
    !%Variable KPoints
    !%Type block
    !%Section Mesh::KPoints
    !%Description
    !% This block defines an explicit set of <i>k</i>-points and their weights for
    !% a periodic-system calculation. The first column is the weight
    !% of each <i>k</i>-point and the following are the components of the <i>k</i>-point
    !% vector. You only need to specify the components for the
    !% periodic directions. Note that the <i>k</i>-points should be given in
    !% Cartesian coordinates (not in reduced coordinates), in the units of inverse length.
    !% The weights will be renormalized so they sum to 1 (and must be rational numbers).
    !%
    !% For example, if you want to include only the Gamma point, you can
    !% use:
    !%
    !% <tt>%KPoints
    !% <br>&nbsp;&nbsp;1.0 | 0 | 0 | 0
    !% <br>%</tt>
    !%
    !%End
 
    !%Variable KPointsReduced
    !%Type block
    !%Section Mesh::KPoints
    !%Description
    !% Same as the block <tt>KPoints</tt> but this time the input is given in reduced
    !% coordinates, <i>i.e.</i>
    !% what <tt>Octopus</tt> writes in a line in the ground-state standard output as
    !%
    !% <tt>#k =   1, k = (    0.154000,    0.154000,    0.154000)</tt>.
    !%End
 
    reduced = .false.
    if (parse_block(namespace, 'KPoints', blk) /= 0) then
      if (parse_block(namespace, 'KPointsReduced', blk) == 0) then
        reduced = .true.
      else
        ! This case should really never happen
        write(message(1),'(a)') 'Internal error loading user-defined k-point list.'
        call messages_fatal(1, namespace=namespace)
      end if
    end if
 
    call kpoints_grid_init(dim, user_kpoints_grid, parse_block_n(blk), 1)
 
    user_kpoints_grid%red_point = m_zero
    user_kpoints_grid%point = m_zero
    user_kpoints_grid%weight = m_zero
    user_kpoints_grid%shifts = m_zero
 
    if (reduced) then
      do ik = 1, user_kpoints_grid%npoints
        call parse_block_float(blk, ik - 1, 0, user_kpoints_grid%weight(ik))
        do idir = 1, dim
          call parse_block_float(blk, ik - 1, idir, user_kpoints_grid%red_point(idir, ik))
        end do
        ! generate also the absolute coordinates
        call kpoints_to_absolute(this%latt, user_kpoints_grid%red_point(:, ik), user_kpoints_grid%point(:, ik))
      end do
    else
      do ik = 1, user_kpoints_grid%npoints
        call parse_block_float(blk, ik - 1, 0, user_kpoints_grid%weight(ik))
        do idir = 1, dim
          call parse_block_float(blk, ik - 1, idir, user_kpoints_grid%point(idir, ik), unit_one/units_inp%length)
        end do
        ! generate also the reduced coordinates
        call kpoints_to_reduced(this%latt, user_kpoints_grid%point(:, ik), user_kpoints_grid%red_point(:, ik))
      end do
    end if
    call parse_block_end(blk)
 
    this%nik_skip = 0
    if (any(user_kpoints_grid%weight < m_epsilon)) then
      call messages_experimental('K-points with zero weight')
      message(1) = "Found k-points with zero weight. They are excluded from density calculation"
      call messages_warning(1, namespace=namespace)
      ! count k-points with zero weight and  make sure the points are given in
      ! a block after all regular k-points. This is for convenience, so they can be skipped
      ! easily and not a big restraint for the user who has to provide the k-points
      ! explicitly anyway.
      do ik=1,user_kpoints_grid%npoints
        if (user_kpoints_grid%weight(ik) < m_epsilon) then
          ! check there are no points with positive weight following a zero weighted one
          if (ik < user_kpoints_grid%npoints) then
            if (user_kpoints_grid%weight(ik+1) > m_epsilon) then
              message(1) = "K-points with zero weight must follow all regular k-points in a block"
              call messages_fatal(1, namespace=namespace)
            end if
          end if
          this%nik_skip = this%nik_skip + 1
          user_kpoints_grid%weight(ik) = m_zero
        end if
      end do
    end if
    ! renormalize weights
    weight_sum = sum(user_kpoints_grid%weight)
    if (weight_sum < m_epsilon) then
      message(1) = "k-point weights must sum to a positive number."
      call messages_fatal(1, namespace=namespace)
    end if
    user_kpoints_grid%weight = user_kpoints_grid%weight / weight_sum
 
    ! for the moment we do not apply symmetries to user kpoints
    ! call kpoints_grid_copy(this%full, this%reduced)
 
    call kpoints_fold_to_1bz(user_kpoints_grid, this%latt)
 
    call kpoints_grid_addto(this%full   , user_kpoints_grid)
    call kpoints_grid_addto(this%reduced, user_kpoints_grid)
 
    write(message(1), '(a,i4,a)') 'Input: ', user_kpoints_grid%npoints, ' k-points were read from the input file'
    call messages_info(1, namespace=namespace)
 
    call kpoints_grid_end(user_kpoints_grid)
 
    pop_sub(kpoints_read_user_kpoints)
  end subroutine kpoints_read_user_kpoints
 
  ! ---------------------------------------------------------
  subroutine kpoints_end(this)
    type(kpoints_t), intent(inout) :: this
 
    push_sub(kpoints_end)
 
    call kpoints_grid_end(this%full)
    call kpoints_grid_end(this%reduced)
 
    safe_deallocate_a(this%nik_axis)
    safe_deallocate_a(this%niq_axis)
    safe_deallocate_a(this%lk123)
    safe_deallocate_a(this%downsampling)
    safe_deallocate_a(this%symmetry_ops)
    safe_deallocate_a(this%num_symmetry_ops)
    safe_deallocate_a(this%map_full_to_reduced)
    safe_deallocate_a(this%map_full_to_symidx)
    safe_deallocate_a(this%coord_along_path)
 
    pop_sub(kpoints_end)
  end subroutine kpoints_end
 
 
  subroutine kpoints_print_info(this, namespace)
    type(kpoints_t),            intent(in)  :: this
    type(namespace_t),          intent(in)  :: namespace
 
    integer :: idir, is, ik, dim
    character(len=100) :: str_tmp
 
    push_sub(kpoints_print_info)
 
    dim = this%reduced%dim
 
    if (bitand(this%method, kpoints_monkh_pack) /= 0) then
      write(message(1),'(a)') ' '
      write(message(2),'(1x,i5,a)') this%reduced%npoints, ' k-points generated from parameters :'
      write(message(3),'(1x,a)') '---------------------------------------------------'
      write(message(4),'(4x,a)') 'n ='
      do idir = 1, dim
        write(str_tmp,'(i5)') this%nik_axis(idir)
        message(4) = trim(message(4)) // trim(str_tmp)
      end do
      call messages_info(4, namespace=namespace)
 
      do is = 1, this%reduced%nshifts
        write(message(1),'(a)') ' '
        write(message(2),'(4x,a,i1,a)') 's', is, '  ='
        do idir = 1, dim
          write(str_tmp,'(f6.2)') this%reduced%shifts(idir,is)
          message(2) = trim(message(2)) // trim(str_tmp)
        end do
        call messages_info(2, namespace=namespace)
      end do
    end if
 
    write(message(1),'(a)') ' '
    write(message(2),'(a)') '   index |      weight  |             coordinates              |'
    call messages_info(2, namespace=namespace)
 
    do ik = 1, this%reduced%npoints
      write(str_tmp,'(i8,a,f12.6,a)') ik, " | ", this%reduced%weight(ik), " |"
      message(1) =  str_tmp
      do idir = 1, dim
        write(str_tmp,'(f12.6)') this%reduced%red_point(idir, ik)
        message(1) = trim(message(1)) // trim(str_tmp)
      end do
      write(str_tmp,'(a)') "  |"
      message(1) = trim(message(1)) // trim(str_tmp)
      call messages_info(1, namespace=namespace)
    end do
 
    write(message(1),'(a)') ' '
    call messages_info(1, namespace=namespace)
 
    pop_sub(kpoints_print_info)
 
  end subroutine kpoints_print_info
 
  ! ---------------------------------------------------------
  subroutine kpoints_to_absolute(latt, kin, kout)
    type(lattice_vectors_t), intent(in)  :: latt
    real(real64),   intent(in)  :: kin(:)
    real(real64),   intent(out) :: kout(:)
 
    push_sub(kpoints_to_absolute)
 
    kout = matmul(latt%klattice, kin)
 
    pop_sub(kpoints_to_absolute)
  end subroutine kpoints_to_absolute
 
  ! ---------------------------------------------------------
  subroutine kpoints_to_reduced(latt, kin, kout)
    type(lattice_vectors_t), intent(in)  :: latt
    real(real64),            intent(in)  :: kin(:)
    real(real64),            intent(out) :: kout(:)
 
    push_sub(kpoints_to_reduced)
 
    kout = matmul(kin, latt%rlattice) / (m_two*m_pi)
 
    pop_sub(kpoints_to_reduced)
  end subroutine kpoints_to_reduced
 
 
  ! ---------------------------------------------------------
  subroutine kpoints_copy(kin, kout)
    type(kpoints_t), intent(in)    :: kin
    type(kpoints_t), intent(inout) :: kout
 
    push_sub(kpoints_copy)
 
    call kpoints_end(kout)
 
    kout%method = kin%method
 
    call kpoints_grid_copy(kin%full, kout%full)
    call kpoints_grid_copy(kin%reduced, kout%reduced)
 
    kout%use_symmetries = kin%use_symmetries
    kout%use_time_reversal = kin%use_time_reversal
 
    safe_allocate(kout%nik_axis(1:kin%full%dim))
    safe_allocate(kout%niq_axis(1:kin%full%dim))
    safe_allocate(kout%downsampling(1:kin%full%dim))
    kout%nik_axis(:) = kin%nik_axis(:)
    kout%niq_axis(:) = kin%niq_axis(:)
    kout%downsampling(:) = kin%downsampling(:)
 
    if (allocated(kin%lk123)) then
      safe_allocate(kout%lk123(1:kin%full%npoints, 1:kin%full%dim))
      kout%lk123(:, :) = kin%lk123(:, :)
    end if
 
    if (allocated(kin%coord_along_path)) then
      safe_allocate(kout%coord_along_path(1:kin%full%npoints))
      kout%coord_along_path(1:kin%full%npoints) = kin%coord_along_path(1:kin%full%npoints)
    end if
 
    if (allocated(kin%symmetry_ops)) then
      safe_allocate(kout%num_symmetry_ops(1:kin%reduced%npoints))
      safe_allocate(kout%symmetry_ops(1:kin%reduced%npoints, 1:maxval(kin%num_symmetry_ops)))
      kout%num_symmetry_ops(1:kin%reduced%npoints) = kin%num_symmetry_ops(1:kin%reduced%npoints)
      kout%symmetry_ops(1:kin%reduced%npoints, 1:maxval(kin%num_symmetry_ops)) = &
        kin%symmetry_ops(1:kin%reduced%npoints, 1:maxval(kin%num_symmetry_ops))
    end if
    if (allocated(kin%map_full_to_symidx)) then
      safe_allocate(kout%map_full_to_symidx(1:kin%full%npoints))
      kout%map_full_to_symidx(1:kin%full%npoints) = kin%map_full_to_symidx(1:kin%full%npoints)
    end if
    if (allocated(kin%map_full_to_reduced)) then
      safe_allocate(kout%map_full_to_reduced(1:kin%full%npoints))
      kout%map_full_to_reduced(1:kin%full%npoints) = kin%map_full_to_reduced(1:kin%full%npoints)
    end if
 
    pop_sub(kpoints_copy)
  end subroutine kpoints_copy
 
 
  ! ----------------------------------------------------------
  integer pure function kpoints_number(this) result(number)
    type(kpoints_t), intent(in) :: this
 
    number = this%reduced%npoints
 
  end function kpoints_number
 
 
  ! ----------------------------------------------------------
  pure function kpoints_get_point(this, ik, absolute_coordinates) result(point)
    class(kpoints_t),  intent(in) :: this
    integer,           intent(in) :: ik
    logical, optional, intent(in) :: absolute_coordinates
    real(real64)                  :: point(1:this%full%dim)
 
    if (optional_default(absolute_coordinates, .true.)) then
      point(1:this%full%dim) = this%reduced%point(1:this%full%dim, ik)
    else
      point(1:this%full%dim) = this%reduced%red_point(1:this%full%dim, ik)
    end if
 
  end function kpoints_get_point
 
 
  ! ----------------------------------------------------------
  real(real64) pure function kpoints_get_weight(this, ik) result(weight)
    class(kpoints_t), intent(in) :: this
    integer,          intent(in) :: ik
 
    weight = this%reduced%weight(ik)
 
  end function kpoints_get_weight
 
 
  ! ----------------------------------------------------------
  integer function kpoints_get_equiv(this, ik) result(equiv)
    class(kpoints_t), intent(in) :: this
    integer,          intent(in) :: ik
 
    assert(allocated(this%map_full_to_reduced))
    equiv = this%map_full_to_reduced(ik)
 
  end function kpoints_get_equiv
 
 
  ! ----------------------------------------------------------
  subroutine kpoints_grid_generate(dim, naxis, nshifts, shift, kpoints, lk123)
    integer,           intent(in)  :: dim
    integer,           intent(in)  :: naxis(1:dim)
    integer,           intent(in)  :: nshifts
    real(real64),      intent(in)  :: shift(:, :)
    real(real64),      intent(out) :: kpoints(:, :)
    integer, optional, intent(out) :: lk123(:, :)
    !                                                 !< running from 0 to naxis(1:3).
 
    real(real64)         :: dx(dim)
    integer              :: ii, jj, divisor, ik, idir, npoints, is, ix(dim)
    integer, allocatable :: lk123_(:, :), idx(:)
 
    push_sub(kpoints_grid_generate)
 
    dx = m_one/(m_two*naxis)
 
    npoints = product(naxis)
 
    safe_allocate(idx(1:npoints*nshifts))
    if (present(lk123)) then
      safe_allocate(lk123_(1:npoints*nshifts,1:dim))
    end if
 
    do is = 1, nshifts
      do ii = 0, npoints - 1
        ik = npoints*is - ii
        jj = ii
        divisor = npoints
 
        do idir = 1, dim
          divisor = divisor / naxis(idir)
          ix(idir) = jj / divisor + 1
          jj = mod(jj, divisor)
 
          kpoints(idir, ik) = (m_two*ix(idir) - m_one*naxis(idir) + m_two*shift(idir,is))*dx(idir)
          !A default shift of +0.5 is including in case if (mod(naxis(idir), 2) /= 0)
 
          !bring back point to first Brillouin zone, except for points at 1/2
          if (abs(kpoints(idir, ik) - m_half) > 1e-14_real64) then
            kpoints(idir, ik) = mod(kpoints(idir, ik) + m_half, m_one) - m_half
          end if
 
        end do
        if (present(lk123)) then
          lk123_(ik, :) = ix
          idx(ik) = ik
        end if
      end do
    end do
 
    call kpoints_deterministic_sort(naxis, kpoints, nshifts=nshifts, indices=idx)
 
    if (present(lk123)) then
      do ik = 1, npoints*nshifts
        lk123(ik,:) = lk123_(idx(ik),:)
      end do
      safe_deallocate_a(lk123_)
    end if
 
    safe_deallocate_a(idx)
 
    pop_sub(kpoints_grid_generate)
 
  end subroutine kpoints_grid_generate
 
 
  subroutine kpoints_deterministic_sort(nik_axis, kpoints, nshifts, indices)
    integer,      intent(in   ) :: nik_axis(:)
    real(real64), intent(inout) :: kpoints(:, :)
    !                                              In: Unsorted
    !                                              Out: Sorted
    integer, optional,  intent(in) :: nshifts
    integer, optional, intent(out) :: indices(:)
    ! original column chosen for sorted position ik_sorted
 
    integer                   :: ndim, nk, ik, nshifts_
    real(real64)              :: dk(size(kpoints, 1))
    integer,      allocatable :: idx(:)
    real(real64), allocatable :: shell(:), kpoints_unsorted(:, :)
 
    push_sub(kpoints_deterministic_sort)
 
    ndim = size(kpoints, 1)
    nshifts_ = optional_default(nshifts, 1)
    nk = product(nik_axis) * nshifts_
 
    ! The MP generator is defined as k = (2*ix - N + 2*shift) * dk
    ! Therefore, for a single step where ix = 1,
    ! Delta k = 2 * dk, so to be consistent dk = 1 / 2N, such that Delta k = 1/N
    dk = m_one / (m_two * nik_axis)
 
    safe_allocate(shell(1:nk))
    safe_allocate(idx(1:nk))
    safe_allocate_source(kpoints_unsorted(1:ndim, 1:nk), kpoints)
 
    ! Squared length of the index-space coordinate
    do ik = 1, nk
      shell(ik) = sum((kpoints_unsorted(:, ik) / dk)**2)
    enddo
 
    call robust_sort_by_abs(shell, kpoints_unsorted, idx)
    safe_deallocate_a(shell)
 
    do ik = 1, nk
      kpoints(:, ik) = kpoints_unsorted(:, idx(ik))
    end do
 
    if (present(indices)) indices = idx
 
    safe_deallocate_a(idx)
    safe_deallocate_a(kpoints_unsorted)
 
    pop_sub(kpoints_deterministic_sort)
 
  end subroutine kpoints_deterministic_sort
 
 
  subroutine kpoints_path_generate(dim, latt, nkpoints, nsegments, resolution, highsympoints, kpoints, coord)
    integer,                 intent(in)  :: dim
    type(lattice_vectors_t), intent(in)  :: latt
    integer,                 intent(in)  :: nkpoints
    integer,                 intent(in)  :: nsegments
    integer,                 intent(in)  :: resolution(:)
    real(real64),            intent(in)  :: highsympoints(:,:)
    real(real64),            intent(out) :: kpoints(1:dim, 1:nkpoints)
    real(real64),            intent(out) :: coord(1:nkpoints)
 
    integer :: is, ik, kpt_ind
    real(real64)   :: length, total_length, accumulated_length
    real(real64)   :: kpt1(dim), kpt2(dim), vec(dim)
 
    push_sub(kpoints_path_generate)
 
    total_length = m_zero
    assert(ubound(highsympoints, dim=2) >= nsegments+1)
    !We first compute the total length of the k-point path
    do is = 1, nsegments
      ! We need to work in abolute coordinates to get the correct path length
      call kpoints_to_absolute(latt, highsympoints(:,is), kpt1)
      call kpoints_to_absolute(latt, highsympoints(:,is+1), kpt2)
 
      vec = kpt2 - kpt1
      length = norm2(vec)
      if (resolution(is) > 0) total_length = total_length + length
    end do
 
    accumulated_length = m_zero
    kpt_ind = 0
    !Now we generate the points
    do is = 1, nsegments
      ! We need to work in abolute coordinates to get the correct path length
      call kpoints_to_absolute(latt, highsympoints(:, is), kpt1)
      call kpoints_to_absolute(latt, highsympoints(:, is+1), kpt2)
 
      vec = kpt2 - kpt1
      length = norm2(vec)
      if (length > m_epsilon) vec = vec/length
 
      do ik = 1, resolution(is)
        kpt_ind = kpt_ind +1
        coord(kpt_ind) = accumulated_length + (ik-1)*length/resolution(is)
        kpoints(:, kpt_ind) = kpt1 + (ik-1)*length/resolution(is)*vec
      end do
      if (resolution(is) > 0) accumulated_length = accumulated_length + length
    end do
    !We add the last point
    kpt_ind = kpt_ind +1
    call kpoints_to_absolute(latt, highsympoints(:,nsegments+1), kpt1)
    coord(kpt_ind) = accumulated_length
    kpoints(:, kpt_ind) =  kpt1
 
    !The length of the total path is arbitrarily put to 1
    coord = coord/total_length
 
    pop_sub(kpoints_path_generate)
  end subroutine kpoints_path_generate
 
 
  ! --------------------------------------------------------------------------------------------
  subroutine kpoints_grid_reduce(symm, time_reversal, nkpoints, dim, kpoints, weights, &
    map_full_to_reduced, map_full_to_symidx, symm_ops, num_symm_ops)
    type(symmetries_t), intent(in)    :: symm
    logical,            intent(in)    :: time_reversal
    integer,            intent(inout) :: nkpoints
    integer,            intent(in)    :: dim
    real(real64),       intent(inout) :: kpoints(1:dim,1:nkpoints)
    real(real64),       intent(out)   :: weights(1:nkpoints)
    integer,            intent(out)   :: map_full_to_reduced(1:nkpoints)
    integer,            intent(out)   :: map_full_to_symidx(1:nkpoints)
    integer,            intent(out)   :: symm_ops(:, :)
    integer,            intent(out)   :: num_symm_ops(:)
 
    integer :: nreduced
    real(real64), allocatable :: reduced(:, :)
 
    integer ik, iop, ik2
    real(real64) :: tran(dim), diff(dim)
    real(real64), allocatable :: kweight(:)
 
    real(real64) :: prec
 
    push_sub(kpoints_grid_reduce)
 
    ! In case of really dense k-point grid, 1/nkpoints is might be smaller
    ! the symprec, causing problems
    ! Therefore we use PREC in the following
    prec = min(symprec, m_one/(nkpoints*100))
 
    ! reduce to irreducible zone
 
    ! kmap is used to mark reducible k-points and also to
    ! map reducible to irreducible k-points
 
    safe_allocate(kweight(1:nkpoints))
    safe_allocate(reduced(1:dim, 1:nkpoints))
 
    kweight = m_one / nkpoints
 
    nreduced = 0
    map_full_to_reduced(:) = -1 ! for safety
    map_full_to_symidx(:) = -1 ! for safety
 
    num_symm_ops = 1
    symm_ops(:, 1) = symmetries_identity_index(symm)
 
    do ik = 1, nkpoints
      if (kweight(ik) < prec) cycle
 
      ! new irreducible point
      ! has reduced non-zero weight
      nreduced = nreduced + 1
      reduced(:, nreduced) = kpoints(:, ik)
      map_full_to_reduced(ik) = nreduced
      map_full_to_symidx(ik) = 1
 
      !No need to check Gamma
      if (maxval(abs(kpoints(:, ik))) < m_epsilon) cycle
 
      if (ik == nkpoints) cycle
 
      ! operate with the symmetry operations
      do iop = 1, symmetries_number(symm)
        if (iop == symmetries_identity_index(symm) .and. &
          .not. time_reversal) cycle ! no need to check for the identity
 
        call symmetries_apply_kpoint_red(symm, iop, reduced(1:dim, nreduced), tran)
        !We remove potential umklapp
        tran = tran - anint(tran + m_half*prec)
 
        ! remove (mark) k-points related to irreducible reduced by symmetry
        do ik2 = ik + 1, nkpoints
          if (kweight(ik2) < prec) cycle
 
          if (.not. iop == symmetries_identity_index(symm)) then ! no need to check for the identity
            diff = tran - kpoints(:, ik2)
            diff = diff - anint(diff)
 
            ! both the transformed rk ...
            if (sum(abs(diff)) < prec) then
              kweight(ik) = kweight(ik) + kweight(ik2)
              kweight(ik2) = m_zero
              map_full_to_reduced(ik2) = nreduced
              map_full_to_symidx(ik2) = num_symm_ops(nreduced) + 1
              weights(nreduced) = kweight(ik)
              num_symm_ops(nreduced) = num_symm_ops(nreduced) + 1
              symm_ops(nreduced, num_symm_ops(nreduced)) = iop
              cycle
            end if
          end if
 
          if (time_reversal) then
            diff = tran + kpoints(:, ik2)
            diff = diff - anint(diff)
 
            ! and its inverse
            if (sum(abs(diff)) < prec) then
              kweight(ik) = kweight(ik) + kweight(ik2)
              kweight(ik2) = m_zero
              map_full_to_reduced(ik2) = nreduced
              map_full_to_symidx(ik2) = num_symm_ops(nreduced) + 1
              weights(nreduced) = kweight(ik)
              num_symm_ops(nreduced) = num_symm_ops(nreduced) + 1
              !We mark the symmetry+time-reversal operation as negative
              symm_ops(nreduced, num_symm_ops(nreduced)) = -iop
            end if
          end if
        end do
      end do
    end do
 
    assert(sum(weights(1:nreduced)) - m_one < prec)
    assert(all(1 <= map_full_to_reduced(:)) .and. all(map_full_to_reduced(:) <= nreduced))
    assert(all(1 <= map_full_to_symidx(:)) .and. all(map_full_to_symidx(:) <= size(symm_ops, 2)))
 
    nkpoints = nreduced
    do ik = 1, nreduced
      kpoints(:, ik) = reduced(:, ik)
    end do
 
    safe_deallocate_a(kweight)
    safe_deallocate_a(reduced)
 
    pop_sub(kpoints_grid_reduce)
  end subroutine kpoints_grid_reduce
 
 
  ! --------------------------------------------------------------------------------------------
  subroutine kpoints_fold_to_1bz(grid, latt)
    type(kpoints_grid_t),    intent(inout) :: grid
    type(lattice_vectors_t), intent(in)    :: latt
 
    integer :: ig1, ig2, ik
    real(real64) :: Gvec(grid%dim, 3**grid%dim), Gvec_cart(grid%dim, 3**grid%dim)
    real(real64) :: vec(grid%dim), kpt(grid%dim)
    real(real64) :: d, dmin
 
    push_sub(kpoints_fold_to_1bz)
 
    !We only need to compute the first G-vectors
    do ig1 = 1, grid%dim
      do ig2 = 1, 3**grid%dim
        gvec(ig1, ig2) = modulo((ig2 - 1)/grid%dim**(grid%dim - ig1), 3) - 1
      end do
    end do
 
    do ig1 = 1, 3**grid%dim
      call kpoints_to_absolute(latt, gvec(:,ig1), gvec_cart(:,ig1))
    end do
 
    do ik = 1, grid%npoints
 
      dmin = 1e10_real64
      do ig1 = 1, 3**grid%dim
        vec = gvec_cart(:,ig1) - grid%point(:,ik)
        d = real(sum(vec**2), real32) !Conversion to simple precision
        !                        To avoid numerical error problems
        if (d < dmin) then
          dmin = d
          ig2 = ig1
        end if
      end do
      kpt = grid%red_point(:,ik) - gvec(:,ig2)
      call kpoints_to_absolute(latt, kpt, grid%point1BZ(:, ik))
    end do
 
    pop_sub(kpoints_fold_to_1bz)
  end subroutine kpoints_fold_to_1bz
 
 
  ! ---------------------------------------------------------
  subroutine kpoints_write_info(this, iunit, namespace, absolute_coordinates)
    class(kpoints_t),            intent(in) :: this
    integer,           optional, intent(in) :: iunit
    type(namespace_t), optional, intent(in) :: namespace
    logical,           optional, intent(in) :: absolute_coordinates
 
    integer :: ik, idir
    character(len=100) :: str_tmp
    character :: index
 
    push_sub(kpoints_write_info)
 
    call messages_print_with_emphasis('Brillouin zone sampling', iunit, namespace)
 
    if (this%method == kpoints_monkh_pack) then
 
      call messages_write('Dimensions of the k-point grid      =')
      do idir = 1, this%full%dim
        call messages_write(this%nik_axis(idir), fmt = '(i3,1x)')
      end do
      call messages_new_line()
 
      if (.not. all(this%downsampling(1:this%full%dim) == 1)) then
        call messages_write('Dimensions of the q-point grid      =')
        do idir = 1, this%full%dim
          call messages_write(this%niq_axis(idir), fmt = '(i3,1x)')
        end do
        call messages_new_line()
      end if
 
      call messages_write('Total number of k-points            =')
      call messages_write(this%full%npoints)
      call messages_new_line()
 
      call messages_write('Number of symmetry-reduced k-points =')
      call messages_write(this%reduced%npoints)
 
      call messages_info(iunit=iunit, namespace=namespace)
 
    else
 
      call messages_write('Total number of k-points            =')
      call messages_write(this%full%npoints)
      call messages_new_line()
      call messages_info(iunit=iunit, namespace=namespace)
 
    end if
 
    call messages_new_line()
    call messages_write('List of k-points:')
    call messages_info(iunit=iunit, namespace=namespace)
 
    write(message(1), '(6x,a)') 'ik'
    do idir = 1, this%full%dim
      index = index2axis(idir)
      write(str_tmp, '(9x,2a)') 'k_', index
      message(1) = trim(message(1)) // trim(str_tmp)
    end do
    write(str_tmp, '(6x,a)') 'Weight'
    message(1) = trim(message(1)) // trim(str_tmp)
    message(2) = '---------------------------------------------------------'
    call messages_info(2, iunit)
 
    do ik = 1, kpoints_number(this)
      write(message(1),'(i8,1x)') ik
      do idir = 1, this%full%dim
        if (optional_default(absolute_coordinates, .false.)) then
          write(str_tmp,'(f12.4)') this%reduced%point(idir, ik)
        else
          write(str_tmp,'(f12.4)') this%reduced%red_point(idir, ik)
        end if
        message(1) = trim(message(1)) // trim(str_tmp)
      end do
      write(str_tmp,'(f12.4)') kpoints_get_weight(this, ik)
      message(1) = trim(message(1)) // trim(str_tmp)
      call messages_info(1, iunit, namespace=namespace)
    end do
 
    call messages_info(iunit=iunit, namespace=namespace)
 
    call messages_print_with_emphasis(iunit=iunit, namespace=namespace)
 
    pop_sub(kpoints_write_info)
  end subroutine kpoints_write_info
 
 
  ! ---------------------------------------------------------
  logical pure function kpoints_point_is_gamma(this, ik) result(is_gamma)
    class(kpoints_t),    intent(in) :: this
    integer,            intent(in) :: ik
 
    is_gamma = (maxval(abs(kpoints_get_point(this, ik))) < m_epsilon)
 
  end function kpoints_point_is_gamma
 
  !--------------------------------------------------------
 
  integer pure function kpoints_get_num_symmetry_ops(this, ik) result(num)
    type(kpoints_t),    intent(in) :: this
    integer,            intent(in) :: ik
 
    if (this%use_symmetries) then
      num = this%num_symmetry_ops(ik)
    else
      num = 1
    end if
 
  end function kpoints_get_num_symmetry_ops
  !--------------------------------------------------------
 
  integer pure function kpoints_get_symmetry_ops(this, ik, index) result(iop)
    type(kpoints_t),    intent(in) :: this
    integer,            intent(in) :: ik
    integer,            intent(in) :: index
 
    if (this%use_symmetries) then
      iop = this%symmetry_ops(ik, index)
    else
      iop = 1
    end if
 
  end function kpoints_get_symmetry_ops
 
  !--------------------------------------------------------
  integer pure function kpoints_get_full_symmetry_op_index(this, ifull) result(iiop)
    class(kpoints_t), intent(in) :: this
    integer,          intent(in) :: ifull
 
    if (this%use_symmetries) then
      iiop = this%map_full_to_symidx(ifull)
    else
      iiop = 1
    end if
 
  end function kpoints_get_full_symmetry_op_index
 
  !--------------------------------------------------------
  logical pure function kpoints_is_valid_symmetry(this, ik, index) result(valid)
    type(kpoints_t),    intent(in) :: this
    integer,            intent(in) :: ik
    integer,            intent(in) :: index
 
    integer :: iop
 
    valid = .false.
    if (this%use_symmetries) then
      do iop = 1, this%num_symmetry_ops(ik)
        if (this%symmetry_ops(ik, iop) == index) then
          valid = .true.
          return
        end if
      end do
    else
      valid = .true.
    end if
 
  end function kpoints_is_valid_symmetry
 
 
  !--------------------------------------------------------
  integer function kpoints_kweight_denominator(this)
    type(kpoints_t),    intent(in) :: this
 
    integer :: denom, nik, nik_skip
 
    push_sub(kpoints_kweight_denominator)
 
    nik = this%full%npoints
    nik_skip = this%nik_skip
 
    if (this%method == kpoints_monkh_pack) then
      kpoints_kweight_denominator = this%full%npoints
    else
      kpoints_kweight_denominator = 0
      ! NB largest reasonable value is: # k-points x 48. from space-group symmetries
      do denom = 1, 100000
        if (all(abs(int(this%full%weight(1:nik-nik_skip)*denom + 10_real64*m_epsilon) - &
          this%full%weight(1:nik-nik_skip)*denom) < 100_real64*m_epsilon)) then
          kpoints_kweight_denominator = denom
          exit
        end if
      end do
    end if
 
    pop_sub(kpoints_kweight_denominator)
  end function kpoints_kweight_denominator
 
  !--------------------------------------------------------
  logical  pure function kpoints_have_zero_weight_path(this) result(have_zerow)
    class(kpoints_t),    intent(in) :: this
 
    if (this%nik_skip > 0) then
      have_zerow = .true.
    else
      have_zerow = .false.
    end if
 
  end function kpoints_have_zero_weight_path
 
  !--------------------------------------------------------
  integer pure function kpoints_get_kpoint_method(this)
    class(kpoints_t),    intent(in) :: this
 
    kpoints_get_kpoint_method = this%method
  end function kpoints_get_kpoint_method
 
  !--------------------------------------------------------
  real(real64) pure function kpoints_get_path_coord(this, ind) result(coord)
    type(kpoints_t),    intent(in) :: this
    integer,            intent(in) :: ind
 
    coord = this%coord_along_path(ind)
  end function
 
 
 
  !--------------------------------------------------------
  subroutine kpoints_check_symmetries(grid, symm, dim, time_reversal, namespace)
    type(kpoints_grid_t), intent(in) :: grid
    type(symmetries_t),   intent(in) :: symm
    integer,              intent(in) :: dim
    logical,              intent(in) :: time_reversal
    type(namespace_t),    intent(in) :: namespace
 
    integer, allocatable :: kmap(:)
    real(real64) :: kpt(dim), diff(dim)
    integer :: nk, ik, ik2, iop
    type(distributed_t) :: kpt_dist
 
    push_sub(kpoints_check_symmetries)
 
    nk = grid%npoints
 
    !We distribute the k-points here for this routine, independently of the rest of the code
    call distributed_init(kpt_dist, nk, mpi_world%comm, "kpt_check")
 
    !A simple map to tell if the k-point as a matching symmetric point or not
    safe_allocate(kmap(kpt_dist%start:kpt_dist%end))
 
    do iop = 1, symmetries_number(symm)
      if (iop == symmetries_identity_index(symm) .and. &
        .not. time_reversal) cycle
 
      do ik = kpt_dist%start, kpt_dist%end
        kmap(ik) = ik
      end do
 
      do ik = kpt_dist%start, kpt_dist%end
        !We apply the symmetry
        call symmetries_apply_kpoint_red(symm, iop, grid%red_point(:, ik), kpt)
        !We remove potential umklapp
        kpt = kpt - anint(kpt + m_half*symprec)
 
        ! remove (mark) k-points which already have a symmetric point
        do ik2 = 1, nk
 
          if (iop /= symmetries_identity_index(symm)) then
            diff = kpt - grid%red_point(:, ik2)
            diff = diff - anint(diff)
            !We found point corresponding to the symmetric kpoint
            if (sum(abs(diff)) < symprec) then
              kmap(ik) = -ik2
              exit
            end if
          end if
 
          if (time_reversal) then
            diff = kpt + grid%red_point(:, ik2)
            diff = diff - anint(diff)
            !We found point corresponding to the symmetric kpoint
            if (sum(abs(diff)) < symprec) then
              kmap(ik) = -ik2
              exit
            end if
          end if
 
        end do
        !In case we have not found a symnetric k-point...
        if (kmap(ik) == ik) then
          write(message(1),'(a,i5,a2,3(f7.3,a2),a)') "The reduced k-point ", ik, " (", &
            grid%red_point(1, ik), ", ", grid%red_point(2, ik), ", ", grid%red_point(3, ik),  &
            ") ", "has no symmetric in the k-point grid for the following symmetry"
          write(message(2),'(i5,1x,a,2x,3(3i4,2x))') iop, ':', transpose(symm_op_rotation_matrix_red(symm%ops(iop)))
          message(3) = "Change your k-point grid or use KPointsUseSymmetries=no."
          call messages_fatal(3, namespace=namespace)
        end if
      end do
    end do
 
    safe_deallocate_a(kmap)
 
    call distributed_end(kpt_dist)
 
    pop_sub(kpoints_check_symmetries)
  end subroutine kpoints_check_symmetries
 
  !--------------------------------------------------------
  logical function kpoints_is_compatible_downsampling(kpt, ik, iq) result(compatible)
    type(kpoints_t), intent(in) :: kpt
    integer,              intent(in) :: ik
    integer,              intent(in) :: iq
 
    integer :: dim, idim
    real(real64) :: diff(kpt%reduced%dim), red(kpt%reduced%dim)
 
    push_sub(kpoints_is_compatible_downsampling)
 
    compatible = .true.
    dim = kpt%reduced%dim
 
    !No downsampling. We use all k-points
    if (all(kpt%downsampling == 1)) then
      pop_sub(kpoints_is_compatible_downsampling)
      return
    end if
 
    assert(kpt%method == kpoints_monkh_pack)
 
    diff = kpt%reduced%red_point(:, ik) - kpt%reduced%red_point(:, iq)
    do idim = 1, dim
      !We remove potential umklapp
      diff(idim) = diff(idim) - anint(diff(idim) + m_half*symprec)
      red(idim) = diff(idim)*kpt%nik_axis(idim)/real(kpt%downsampling(idim))
      if (abs(red(idim) - anint(red(idim))) > m_epsilon) then
        compatible = .false.
        pop_sub(kpoints_is_compatible_downsampling)
        return
      end if
    end do
 
    pop_sub(kpoints_is_compatible_downsampling)
  end function kpoints_is_compatible_downsampling
 
  ! ---------------------------------------------------------
  integer function kpoints_nkpt_in_path(this) result(npath)
    class(kpoints_t), intent(in) :: this
 
    npath = 0
    if (allocated(this%coord_along_path)) then
      npath = SIZE(this%coord_along_path)
    end if
 
  end function kpoints_nkpt_in_path
 
  ! ---------------------------------------------------------
  logical function kpoints_gamma_only(this) result(gamma_only)
    class(kpoints_t), intent(in) :: this
 
    push_sub(kpoints_gamma_only)
 
    gamma_only = kpoints_number(this) == 1 .and. kpoints_point_is_gamma(this, 1)
 
    pop_sub(kpoints_gamma_only)
  end function kpoints_gamma_only
 
  ! ---------------------------------------------------------
  subroutine kpoints_lattice_vectors_update(this, new_latt)
    class(kpoints_t),        intent(inout) :: this
    type(lattice_vectors_t), intent(in)    :: new_latt
 
    integer :: ik
 
    push_sub(kpoints_lattice_vectors_update)
 
    this%latt = new_latt
 
    do ik = 1, this%reduced%npoints
      call kpoints_to_absolute(this%latt, this%reduced%red_point(:, ik), this%reduced%point(:, ik))
    end do
 
    do ik = 1, this%full%npoints
      call kpoints_to_absolute(this%latt, this%full%red_point(:, ik), this%full%point(:, ik))
    end do
 
    pop_sub(kpoints_lattice_vectors_update)
  end subroutine kpoints_lattice_vectors_update
 
 
end module kpoints_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: