kmesh_get Subroutine

  subroutine kmesh_get(kmesh_input, kmesh_info, print_output, kpt_latt, real_lattice, num_kpts, &
                       gamma_only, seedname, stdout, timer, error, comm)
    !================================================
    !
    !! Main routine to calculate the b-vectors
    !
    !================================================

    use w90_utility, only: utility_compar, utility_recip_lattice, utility_frac_to_cart, &
                           utility_cart_to_frac, utility_inverse_mat
    use w90_io, only: io_stopwatch_start, io_stopwatch_stop
    use w90_types, only: kmesh_info_type, kmesh_input_type, print_output_type, timer_list_type

    implicit none

    ! arguments
    type(print_output_type), intent(in) :: print_output
    type(kmesh_info_type), intent(inout) :: kmesh_info
    type(kmesh_input_type), intent(inout) :: kmesh_input
    type(timer_list_type), intent(inout) :: timer
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm

    character(len=*), intent(in)  :: seedname
    integer, intent(in) :: num_kpts
    integer, intent(in) :: stdout
    real(kind=dp), intent(in) :: real_lattice(3, 3)
    real(kind=dp), intent(in) :: kpt_latt(:, :)
    logical, intent(in) :: gamma_only

    ! local variables
    real(kind=dp), allocatable :: bvec_tmp(:, :), bvec_inp(:, :, :) ! bvec_inp is allocated in kmesh_shell_from_file()
    real(kind=dp), allocatable :: kpt_cart(:, :)
    real(kind=dp), allocatable :: bk_local(:, :, :)
    real(kind=dp), parameter :: eta = 99999999.0_dp    ! eta = very large
    real(kind=dp) :: dist, dnn0, dnn1, bb1, bbn, ddelta
    real(kind=dp) :: dnn(max(kmesh_input%search_shells, 6*kmesh_input%higher_order_n))
    real(kind=dp) :: recip_lattice(3, 3), volume
    real(kind=dp) :: vkpp(3), vkpp2(3), kpbvec(3)
    integer :: bnum, nbvec ! context of kmesh_input%kmesh_shell_from_file

    ! higher-order finite-difference
    integer, allocatable :: lmn(:, :) ! Order in which to search the cells (ordered in dist from origin)
    integer, allocatable :: nnlist_tmp(:, :), nncell_tmp(:, :, :) ![ysl]
    integer, allocatable :: nnshell(:, :)
    integer :: ifound, counter, na, nap, loop_s, loop_b, shell !, nbvec, bnum
    integer :: ifpos, ifneg, ierr, multi(max(kmesh_input%search_shells, 6*kmesh_input%higher_order_n))
    integer :: nlist, nkp, nkp2, l, m, n, ndnn, ndnnx, ndnntot
    integer :: nnsh, nn, nnx, loop, i, j
    integer :: num_first_shells, ndnn2, nnx2, multi_cumulative, lmn_temp(3)
    integer :: num_x((1 + kmesh_input%higher_order_n)*(1 + 2*kmesh_input%higher_order_n))
    integer :: num_y((1 + kmesh_input%higher_order_n)*(1 + 2*kmesh_input%higher_order_n))
    integer :: num_z((1 + kmesh_input%higher_order_n)*(1 + 2*kmesh_input%higher_order_n))
    real(kind=dp) :: bk_latt(3), inv_lattice(3, 3)
    real(kind=dp) :: bweight(kmesh_input%max_shells_h)
    real(kind=dp) :: wb_local(kmesh_input%num_nnmax_h)

    if (print_output%timing_level > 0) call io_stopwatch_start('kmesh: get', timer)

    allocate (bk_local(3, kmesh_input%num_nnmax_h, num_kpts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating bk_local in kmesh_get', comm)
      return
    end if

    allocate (lmn(3, (2*kmesh_input%search_supcell_size + 1)**3), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating lmn in kmesh_get', comm)
      return
    end if

    allocate (nnshell(num_kpts, max(kmesh_input%search_shells, 6*kmesh_input%higher_order_n)), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating nnshell in kmesh_get', comm)
      return
    end if

    call utility_recip_lattice(real_lattice, recip_lattice, volume, error, comm)
    if (allocated(error)) return
    call utility_inverse_mat(recip_lattice, inv_lattice)
    if (print_output%iprint > 0) write (stdout, '(/1x,a)') &
      '*---------------------------------- K-MESH ----------------------------------*'

    ! Sort the cell neighbours so we loop in order of distance from the home shell
    call kmesh_supercell_sort(print_output, recip_lattice, lmn, timer)

    allocate (kpt_cart(3, num_kpts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error allocating kpt_cart in kmesh_get', comm)
      return
    end if
    do nkp = 1, num_kpts
      call utility_frac_to_cart(kpt_latt(:, nkp), kpt_cart(:, nkp), recip_lattice)
    end do

    ! find the distance between k-point 1 and its nearest-neighbour shells
    ! if we have only one k-point, the n-neighbours are its periodic images

    dnn0 = 0.0_dp
    dnn1 = eta
    ndnntot = 0
    do nlist = 1, kmesh_input%search_shells
      do nkp = 1, num_kpts
        do loop = 1, (2*kmesh_input%search_supcell_size + 1)**3
          l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop)
          !
          vkpp = kpt_cart(:, nkp) + matmul(lmn(:, loop), recip_lattice)
          dist = sqrt((kpt_cart(1, 1) - vkpp(1))**2 &
                      + (kpt_cart(2, 1) - vkpp(2))**2 + (kpt_cart(3, 1) - vkpp(3))**2)
          if ((dist .gt. kmesh_input%tol) .and. (dist .gt. dnn0 + kmesh_input%tol)) then
            if (dist .lt. dnn1 - kmesh_input%tol) then
              dnn1 = dist  ! found a closer shell
              counter = 0
            end if
            if (dist .gt. (dnn1 - kmesh_input%tol) .and. dist .lt. (dnn1 + kmesh_input%tol)) then
              counter = counter + 1 ! count the multiplicity of the shell
            end if
          end if
        end do
      end do
      if (dnn1 .lt. eta - kmesh_input%tol) ndnntot = ndnntot + 1
      dnn(nlist) = dnn1
      multi(nlist) = counter
      dnn0 = dnn1
      dnn1 = eta
    end do

    if (print_output%iprint > 0) then
      write (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
      write (stdout, '(1x,a)') '|                    Distance to Nearest-Neighbour Shells                    |'
      write (stdout, '(1x,a)') '|                    ------------------------------------                    |'
      if (trim(print_output%length_unit) == 'Ang') then
        write (stdout, '(1x,a)') '|          Shell             Distance (Ang^-1)          Multiplicity         |'
        write (stdout, '(1x,a)') '|          -----             -----------------          ------------         |'
      else
        write (stdout, '(1x,a)') '|          Shell             Distance (Bohr^-1)         Multiplicity         |'
        write (stdout, '(1x,a)') '|          -----             ------------------         ------------         |'
      end if
      do ndnn = 1, ndnntot
        write (stdout, '(1x,a,11x,i3,17x,f10.6,19x,i4,12x,a)') '|', ndnn, &
          dnn(ndnn)/print_output%lenconfac, multi(ndnn), '|'
      end do
      write (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
    end if

    if (print_output%iprint >= 4) then
      ! Write out all the bvectors
      if (print_output%iprint > 0) then
        write (stdout, '(1x,"|",76(" "),"|")')
        write (stdout, '(1x,a)') '|         Complete list of b-vectors and their lengths                       |'
        write (stdout, '(1x,"|",76(" "),"|")')
        write (stdout, '(1x,"+",76("-"),"+")')
      end if

      allocate (bvec_tmp(3, maxval(multi)), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating bvec_tmp in kmesh_get', comm)
        return
      end if
      bvec_tmp = 0.0_dp
      counter = 0
      do shell = 1, kmesh_input%search_shells
        call kmesh_get_bvectors(kmesh_input, print_output, bvec_tmp(:, 1:multi(shell)), kpt_cart, &
                                recip_lattice, dnn(shell), lmn, 1, multi(shell), num_kpts, &
                                timer, error, comm)
        if (allocated(error)) return

        do loop = 1, multi(shell)
          counter = counter + 1
          if (print_output%iprint > 0) write (stdout, '(a,I4,1x,a,2x,3f12.6,2x,a,2x,f12.6,a)') ' | b-vector  ', counter, ': (', &
            bvec_tmp(:, loop)/print_output%lenconfac, ')', dnn(shell)/print_output%lenconfac, '  |'
        end do
      end do
      deallocate (bvec_tmp, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating bvec_tmp in kmesh_get', comm)
        return
      end if
      if (print_output%iprint > 0) write (stdout, '(1x,"|",76(" "),"|")')
      if (print_output%iprint > 0) write (stdout, '(1x,"+",76("-"),"+")')
    end if

    ! Get the shell weights to satisfy the B1 condition
    if (kmesh_input%kmesh_shell_from_file) then
      ! note: B1 condition is not tested in kmesh_shell_from_file(); test occurs below
      call kmesh_shell_from_file(kmesh_input, print_output, bvec_inp, bweight, dnn, kpt_cart, &
                                 recip_lattice, lmn, multi, num_kpts, seedname, stdout, timer, &
                                 error, comm)
      if (allocated(error)) return

    elseif (kmesh_input%num_shells == 0) then
      call kmesh_shell_automatic(kmesh_input, print_output, bweight, dnn, kpt_cart, recip_lattice, &
                                 lmn, multi, num_kpts, stdout, timer, error, comm)
      if (allocated(error)) return

    elseif (kmesh_input%num_shells > 0) then
      call kmesh_shell_fixed(kmesh_input, print_output, bweight, dnn, kpt_cart, recip_lattice, &
                             lmn, multi, num_kpts, stdout, timer, error, comm)
      if (allocated(error)) return
    end if

    num_first_shells = kmesh_input%num_shells ! for convenience in printout error

    if (print_output%iprint > 0) then
      if (kmesh_input%higher_order_nearest_shells) then
        write (stdout, '(1x,a)', advance='no') '| The following shells and their multiples are used: '
      else
        write (stdout, '(1x,a)', advance='no') '| The following shells are used:                     '
      end if
      do ndnn = 1, kmesh_input%num_shells
        if (ndnn .eq. kmesh_input%num_shells) then
          write (stdout, '(i3,1x)', advance='no') kmesh_input%shell_list(ndnn)
        else
          write (stdout, '(i3,",")', advance='no') kmesh_input%shell_list(ndnn)
        end if
      end do
      do l = 1, 6 - kmesh_input%num_shells
        write (stdout, '(4x)', advance='no')
      end do
      if (kmesh_input%higher_order_nearest_shells) then
        write (stdout, '(20x)', advance='no')
      end if
      write (stdout, '("|")')
    end if

    kmesh_info%nntot = 0
    do loop_s = 1, kmesh_input%num_shells
      kmesh_info%nntot = kmesh_info%nntot + multi(kmesh_input%shell_list(loop_s))
    end do

    if (kmesh_info%nntot > kmesh_input%num_nnmax_h) then
      if (print_output%iprint > 0) then
        write (stdout, '(a,i2,a)') ' **WARNING: kmesh has found >', kmesh_input%num_nnmax_h, ' nearest neighbours**'
        write (stdout, '(a)') ' '
        write (stdout, '(a)') ' This is probably caused by an error in your unit cell specification'
        write (stdout, '(a)') ' '
        write (stdout, '(a)') ' If you think this is not the problem; please send your *.win file to the '
        write (stdout, '(a)') ' wannier90 developers'
        write (stdout, '(a)') ' '
        write (stdout, '(a)') ' The problem may be caused by having accidentally degenerate shells of '
        write (stdout, '(a)') ' kpoints. The solution is then to rerun wannier90 specifying the b-vectors '
        write (stdout, '(a)') ' in each shell.  Give kmesh_shell_from_file=T in the *.win file'
        write (stdout, '(a)') ' and create a *.kshell file:'
        write (stdout, '(a)') ' '
        write (stdout, '(a)') ' $>   cat hexagonal.kshell'
        write (stdout, '(a)') ' $>   1 2'
        write (stdout, '(a)') ' $>   5 6 7 8'
        write (stdout, '(a)') ' '
        write (stdout, '(a)') ' Where each line is a new shell (so num_shells in total)'
        write (stdout, '(a)') ' The elements are the bvectors labelled according to the following '
        write (stdout, '(a)') ' list (last column is distance)'
        write (stdout, '(a)') ' '
      end if

      allocate (bvec_tmp(3, maxval(multi)), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error allocating bvec_tmp in kmesh_get', comm)
        return
      end if
      bvec_tmp = 0.0_dp
      counter = 0
      do shell = 1, kmesh_input%search_shells
        call kmesh_get_bvectors(kmesh_input, print_output, bvec_tmp(:, 1:multi(shell)), kpt_cart, &
                                recip_lattice, dnn(shell), lmn, 1, multi(shell), num_kpts, timer, &
                                error, comm)
        if (allocated(error)) return
        do loop = 1, multi(shell)
          counter = counter + 1
          if (print_output%iprint > 0) write (stdout, '(a,I4,1x,a,2x,3f12.6,2x,a,2x,f12.6,a)') ' | b-vector  ', counter, ': (', &
            bvec_tmp(:, loop)/print_output%lenconfac, ')', dnn(shell)/print_output%lenconfac, '  |'
        end do
      end do
      if (print_output%iprint > 0) write (stdout, '(a)') ' '
      deallocate (bvec_tmp, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error deallocating bvec_tmp in kmesh_get', comm)
        return
      end if

      call set_error_fatal(error, 'kmesh_get: something wrong, found too many nearest neighbours', comm)
      return
    end if

    ! higher-order algorithm: include 2b, 3b, ..., Nb shells, and modify bweights
    if (kmesh_input%higher_order_nearest_shells) then
      if (print_output%iprint > 0) write (stdout, '(a)') &
        ' | WARNING: higher_order_nearest_shells is an experimental feature, and has   |', &
        ' | not been extensively tested.                                               |'
    else
      kmesh_info%nntot = kmesh_info%nntot*kmesh_input%higher_order_n
    end if

    allocate (kmesh_info%nnlist(num_kpts, kmesh_info%nntot), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating nnlist in kmesh_get', comm)
      return
    end if
    allocate (kmesh_info%neigh(num_kpts, kmesh_info%nntot/2), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating neigh in kmesh_get', comm)
      return
    end if
    allocate (kmesh_info%nncell(3, num_kpts, kmesh_info%nntot), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating nncell in kmesh_get', comm)
      return
    end if
    allocate (kmesh_info%wb(kmesh_info%nntot), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating wb in kmesh_get', comm)
      return
    end if
    allocate (kmesh_info%bka(3, kmesh_info%nntot/2), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating bka in kmesh_get', comm)
      return
    end if
    allocate (kmesh_info%bk(3, kmesh_info%nntot, num_kpts), stat=ierr)
    if (ierr /= 0) then
      call set_error_alloc(error, 'Error in allocating bk in kmesh_get', comm)
      return
    end if

    nnx = 0
    do loop_s = 1, kmesh_input%num_shells
      do loop_b = 1, multi(kmesh_input%shell_list(loop_s))
        nnx = nnx + 1
        wb_local(nnx) = bweight(loop_s)
      end do
    end do

    ! Now build up the list of nearest-neighbour shells for each k-point.
    ! nnlist(nkp,1...nnx) points to the nnx neighbours (ordered along increa
    ! shells) of the k-point nkp. nncell(i,nkp,nnth) tells us in which BZ is
    ! nnth nearest-neighbour of the k-point nkp. Construct the nnx b-vectors
    ! go from k-point nkp to each neighbour bk(1:3,nkp,1...nnx).
    ! Comment: Now we have bk(3,nntot,num_kps) 09/04/2006

    if (print_output%iprint > 0) then
      write (stdout, '(1x,a)') '+----------------------------------------------------------------------------+'
      write (stdout, '(1x,a)') '|                        Shell   # Nearest-Neighbours                        |'
      write (stdout, '(1x,a)') '|                        -----   --------------------                        |'
    end if
    !if (index(print_output%devel_flag, 'kmesh_degen') == 0) then
    !
    ! Standard routine
    !
    nnshell = 0
    do nkp = 1, num_kpts
      nnx = 0
      ok: do ndnnx = 1, kmesh_input%num_shells
        ndnn = kmesh_input%shell_list(ndnnx)
        do loop = 1, (2*kmesh_input%search_supcell_size + 1)**3
          l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop)
          vkpp2 = matmul(lmn(:, loop), recip_lattice)
          do nkp2 = 1, num_kpts
            vkpp = vkpp2 + kpt_cart(:, nkp2)
            dist = sqrt((kpt_cart(1, nkp) - vkpp(1))**2 &
                        + (kpt_cart(2, nkp) - vkpp(2))**2 + (kpt_cart(3, nkp) - vkpp(3))**2)
            if ((dist .ge. dnn(ndnn) - kmesh_input%tol) .and. (dist .le. dnn(ndnn) + kmesh_input%tol)) then
              nnx = nnx + 1
              nnshell(nkp, ndnn) = nnshell(nkp, ndnn) + 1
              kmesh_info%nnlist(nkp, nnx) = nkp2
              kmesh_info%nncell(1, nkp, nnx) = l
              kmesh_info%nncell(2, nkp, nnx) = m
              kmesh_info%nncell(3, nkp, nnx) = n
              bk_local(:, nnx, nkp) = vkpp(:) - kpt_cart(:, nkp)
            end if
            !if we have the right number of neighbours we can exit
            if (nnshell(nkp, ndnn) == multi(ndnn)) cycle ok
          end do
        end do
        ! check to see if too few neighbours here
      end do ok
    end do

    ! higher-order algorithm: include 2b, 3b, ..., Nb shells, and modify bweights
    if (.not. kmesh_input%higher_order_nearest_shells) then
      ! update num_shells, shell_list, dnn(distance to shells), multi, etc.
      call kmesh_shell_reconstruct(kmesh_input, num_kpts, multi, dnn, nnshell, bweight)

      ! update bk_local
      ! update nnlist(neighboring kpts), and nncell(G_lmn vectors of neighbors)
      do nkp = 1, num_kpts
        do nn = 2, kmesh_input%higher_order_n
          multi_cumulative = 0
          do ndnn = 1, kmesh_input%num_shells/kmesh_input%higher_order_n !first-order shells
            ! ndnn: index of shells (first order)
            ! ndnn2: index of shells (nn-th-order)
            ndnn2 = (nn - 1)*kmesh_input%num_shells/kmesh_input%higher_order_n + ndnn
            do nnx = 1 + multi_cumulative, multi(ndnn) + multi_cumulative
              counter = 0
              ! nnx: index of bvectors (first order)
              ! nnx2: index of bvectors (nn-th-order)
              nnx2 = (nn - 1)*kmesh_info%nntot/kmesh_input%higher_order_n + nnx
              bk_local(:, nnx2, nkp) = nn*bk_local(:, nnx, nkp)
              ! find nnlist and nncell
              !do loop = 1, (2*kmesh_input%search_supcell_size + 1)**3
              !  l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop)
              !  vkpp2 = matmul(lmn(:, loop), recip_lattice) !G_lmn vector
              !  do nkp2 = 1, num_kpts
              !    vkpp = vkpp2 + kpt_cart(:, nkp2) - kpt_cart(:, nkp)
              !    ! if kp2 - kp1 == bk_local(:, nnx2, :)
              !    call utility_compar(vkpp(1), bk_local(1, nnx2, nkp), ifpos, ifneg)
              !    if (ifpos .eq. 1) then
              !      counter = counter + 1
              !      kmesh_info%nnlist(nkp, nnx2) = nkp2
              !      kmesh_info%nncell(1, nkp, nnx2) = l
              !      kmesh_info%nncell(2, nkp, nnx2) = m
              !      kmesh_info%nncell(3, nkp, nnx2) = n
              !    endif
              !  enddo
              !enddo
              ! do not search supcell
              ! find nnlist(nkp2) and nncell(lmn)
              call utility_cart_to_frac(bk_local(:, nnx2, nkp), bk_latt, inv_lattice)
              lmn_temp(1) = floor(kpt_latt(1, nkp) + bk_latt(1) + 1.e-6_dp) ! e.g. 3.999999999 is 4
              lmn_temp(2) = floor(kpt_latt(2, nkp) + bk_latt(2) + 1.e-6_dp)
              lmn_temp(3) = floor(kpt_latt(3, nkp) + bk_latt(3) + 1.e-6_dp)
              vkpp2 = matmul(lmn_temp, recip_lattice) !G_lmn vector
              vkpp = kpt_cart(:, nkp) + bk_local(:, nnx2, nkp) - vkpp2 ! k_2 = k_1 + Nb - G_lmn
              do nkp2 = 1, num_kpts
                call utility_compar(kpt_cart(:, nkp2), vkpp, ifpos, ifneg)
                if (ifpos .eq. 1) then
                  counter = counter + 1
                  kmesh_info%nnlist(nkp, nnx2) = nkp2
                  kmesh_info%nncell(1, nkp, nnx2) = lmn_temp(1)
                  kmesh_info%nncell(2, nkp, nnx2) = lmn_temp(2)
                  kmesh_info%nncell(3, nkp, nnx2) = lmn_temp(3)
                end if
              end do
              if (counter == 0) then
                call set_error_fatal(error, 'Could not find Nb vectors in kmesh_get', comm)
              end if
              if (counter >= 2) then
                call set_error_fatal(error, 'Error in kmesh_get, try to modify tolerance in utility_compar', comm)
              end if
            end do
            multi_cumulative = multi_cumulative + multi(ndnn)
          end do
        end do
      end do
    end if
    nnx = 0
    do loop_s = 1, kmesh_input%num_shells
      do loop_b = 1, multi(kmesh_input%shell_list(loop_s))
        nnx = nnx + 1
        wb_local(nnx) = bweight(loop_s)
      end do
    end do

    if (kmesh_input%kmesh_shell_from_file) then
      ! should this be moved to the kmesh_shell_from_file as a function to simplify?
      nnshell = 0
      do nkp = 1, num_kpts
        nnx = 0
        ok2: do loop = 1, (2*kmesh_input%search_supcell_size + 1)**3
          l = lmn(1, loop); m = lmn(2, loop); n = lmn(3, loop)
          vkpp2 = matmul(lmn(:, loop), recip_lattice)
          do nkp2 = 1, num_kpts
            vkpp = vkpp2 + kpt_cart(:, nkp2)
            bnum = 0
            do ndnnx = 1, kmesh_input%num_shells
              do nbvec = 1, multi(ndnnx)
                bnum = bnum + 1
                kpbvec = kpt_cart(:, nkp) + bvec_inp(:, nbvec, ndnnx)
                dist = sqrt((kpbvec(1) - vkpp(1))**2 &
                            + (kpbvec(2) - vkpp(2))**2 + (kpbvec(3) - vkpp(3))**2)
                if (abs(dist) < kmesh_input%tol) then
                  nnx = nnx + 1
                  nnshell(nkp, ndnnx) = nnshell(nkp, ndnnx) + 1
                  kmesh_info%nnlist(nkp, bnum) = nkp2
                  kmesh_info%nncell(1, nkp, bnum) = l
                  kmesh_info%nncell(2, nkp, bnum) = m
                  kmesh_info%nncell(3, nkp, bnum) = n
                  bk_local(:, bnum, nkp) = bvec_inp(:, nbvec, ndnnx)
                end if
              end do
            end do
            if (nnx == sum(multi)) exit ok2
          end do
        end do ok2
      end do
    end if ! kmesh_shell_from_file

    if (kmesh_input%higher_order_n .eq. 1 .or. kmesh_input%higher_order_nearest_shells) then
      do ndnnx = 1, kmesh_input%num_shells
        ndnn = kmesh_input%shell_list(ndnnx)
        if (print_output%iprint > 0) then
          write (stdout, '(1x,a,24x,i3,13x,i3,33x,a)') '|', ndnn, nnshell(1, ndnn), '|'
        end if
      end do
    else
      do ndnnx = 1, num_first_shells
        ndnn = kmesh_input%shell_list(ndnnx)
        if (print_output%iprint > 0) then
          write (stdout, '(1x,a,24x,i3,13x,i3,33x,a)') '|', ndnn, nnshell(1, ndnn), '|'
        end if
      end do
      do i = 2, kmesh_input%higher_order_n
        do ndnnx = 1, num_first_shells
          ndnn = kmesh_input%shell_list(ndnnx)
          if (print_output%iprint > 0) then
            write (stdout, '(1x,a,20x,i3,a,i2,13x,i3,33x,a)') '|', i, ' x', ndnn, nnshell(1, ndnn), '|'
          end if
        end do
      end do
    end if
    if (print_output%iprint > 0) write (stdout, '(1x,"+",76("-"),"+")')

    do nkp = 1, num_kpts
      nnx = 0
      do ndnnx = 1, kmesh_input%num_shells
        ndnn = kmesh_input%shell_list(ndnnx)
        do nnsh = 1, nnshell(nkp, ndnn)
          bb1 = 0.0_dp
          bbn = 0.0_dp
          nnx = nnx + 1
          do i = 1, 3
            bb1 = bb1 + bk_local(i, nnx, 1)*bk_local(i, nnx, 1)
            bbn = bbn + bk_local(i, nnx, nkp)*bk_local(i, nnx, nkp)
          end do
          if (abs(sqrt(bb1) - sqrt(bbn)) .gt. kmesh_input%tol) then
            if (print_output%iprint > 0) write (stdout, '(1x,2f10.6)') bb1, bbn
            call set_error_fatal(error, 'Non-symmetric k-point neighbours in kmesh_get', comm)
            return
          end if
        end do
      end do
    end do

    ! now check that the completeness relation is satisfied for every kpoint
    ! We know it is true for kpt=1; but we check the rest to be safe.
    ! Eq. B1 in Appendix B PRB 56 12847 (1997)

    if ((.not. kmesh_input%skip_B1_tests) .and. kmesh_input%higher_order_nearest_shells) then
      do nkp = 1, num_kpts
        do i = 1, kmesh_input%higher_order_n
          if ((.not. kmesh_input%higher_order_nearest_shells) .and. i > 1) exit
          do j = 1, (1 + i)*(1 + 2*i) ! multiset coefficient ((3, 2i)) = (1 + i)*(1 + 2*i), 3: x,y,z, 2i: num. of b
            ! separate cartesian components
            ! e.g. for i=1
            ! j=1: num_x=2, num_y=0, num_z=0
            ! j=2: num_x=1, num_y=1, num_z=0
            ! j=3: num_x=0, num_y=2, num_z=0
            ! j=4: num_x=1, num_y=0, num_z=1
            ! j=5: num_x=0, num_y=1, num_z=1
            ! j=6: num_x=0, num_y=0, num_z=2
            do l = 0, 2*i
              if ((2*i + 1)*l - l*(l - 1)/2 <= j - 1 &
                  .and. (2*i + 1)*(l + 1) - (l + 1)*l/2 > j - 1) then
                num_z(j) = l
                exit
              end if
            end do
            num_y(j) = j - 1 - ((2*i + 1)*num_z(j) - num_z(j)*(num_z(j) - 1)/2)
            num_x(j) = 2*i - num_y(j) - num_z(j)
            ddelta = 0.0_dp
            nnx = 0
            do ndnnx = 1, kmesh_input%num_shells
              ndnn = kmesh_input%shell_list(ndnnx)
              do nnsh = 1, nnshell(1, ndnn)
                nnx = nnx + 1
                ddelta = ddelta + wb_local(nnx)*(bk_local(1, nnx, nkp)**num_x(j)) &
                         *(bk_local(2, nnx, nkp)**num_y(j))*(bk_local(3, nnx, nkp)**num_z(j))
              end do
            end do
            if (i .eq. 1 .and. (j .eq. 1 .or. j .eq. 3 .or. j .eq. 6)) then
              if (abs(ddelta - 1.0_dp) .gt. kmesh_input%tol) then
                if (print_output%iprint > 0) write (stdout, '(1x,3i3,f12.8)') num_x(j), num_y(j), num_z(j), ddelta
                call set_error_fatal(error, 'Eq. (B1) not satisfied in kmesh_get (1)', comm)
              end if
            else
              if (abs(ddelta) .gt. kmesh_input%tol) then
                if (print_output%iprint > 0) write (stdout, '(1x,3i3,f12.8)') num_x(j), num_y(j), num_z(j), ddelta
                call set_error_fatal(error, 'Eq. (B1) not satisfied in kmesh_get (2)', comm)
              end if
            end if
          end do
        end do
      end do
    end if

    if (print_output%iprint > 0) then
      write (stdout, '(1x,a)') '| Completeness relation is fully satisfied [Eq. (B1), PRB 56, 12847 (1997)]  |'
      if ((kmesh_input%higher_order_nearest_shells) .and. (kmesh_input%higher_order_n .gt. 1)) then
        write (stdout, '(1x,a)') '| Completeness relations for higher-order are fully satisfied                |'
      end if
      write (stdout, '(1x,"+",76("-"),"+")')
    end if

    kmesh_info%wbtot = 0.0_dp
    nnx = 0
    do ndnnx = 1, kmesh_input%num_shells
      ndnn = kmesh_input%shell_list(ndnnx)
      do nnsh = 1, nnshell(1, ndnn)
        nnx = nnx + 1
        kmesh_info%wbtot = kmesh_info%wbtot + wb_local(nnx)
      end do
    end do

    kmesh_info%nnh = kmesh_info%nntot/2
    ! make list of bka vectors from neighbours of first k-point
    ! delete any inverse vectors as you collect them
    na = 0
    do nn = 1, kmesh_info%nntot
      ifound = 0
      if (na .ne. 0) then
        do nap = 1, na
          call utility_compar(kmesh_info%bka(1, nap), bk_local(1, nn, 1), ifpos, ifneg)
          if (ifneg .eq. 1) ifound = 1
        end do
      end if
      if (ifound .eq. 0) then
        !         found new vector to add to set
        na = na + 1
        kmesh_info%bka(1, na) = bk_local(1, nn, 1)
        kmesh_info%bka(2, na) = bk_local(2, nn, 1)
        kmesh_info%bka(3, na) = bk_local(3, nn, 1)
      end if
    end do
    if (na .ne. kmesh_info%nnh) then
      call set_error_fatal(error, 'Did not find right number of bk directions kmesh_get', comm)
      return
    end if

    if (print_output%iprint > 0) then
      if (trim(print_output%length_unit) == 'Ang') then
        write (stdout, '(1x,a)') '|                  b_k Vectors (Ang^-1) and Weights (Ang^2)                  |'
        write (stdout, '(1x,a)') '|                  ----------------------------------------                  |'
      else
        write (stdout, '(1x,a)') '|                 b_k Vectors (Bohr^-1) and Weights (Bohr^2)                 |'
        write (stdout, '(1x,a)') '|                 ------------------------------------------                 |'
      end if
      write (stdout, '(1x,a)') '|            No.         b_k(x)      b_k(y)      b_k(z)        w_b           |'
      write (stdout, '(1x,a)') '|            ---        --------------------------------     --------        |'
      do i = 1, kmesh_info%nntot
        write (stdout, '(1x,"|",11x,i3,5x,3f12.6,3x,f10.6,8x,"|")') &
          i, (bk_local(j, i, 1)/print_output%lenconfac, j=1, 3), wb_local(i)*print_output%lenconfac**2
      end do
      write (stdout, '(1x,"+",76("-"),"+")')
      if (trim(print_output%length_unit) == 'Ang') then
        write (stdout, '(1x,a)') '|                           b_k Directions (Ang^-1)                          |'
        write (stdout, '(1x,a)') '|                           -----------------------                          |'
      else
        write (stdout, '(1x,a)') '|                           b_k Directions (Bohr^-1)                         |'
        write (stdout, '(1x,a)') '|                           ------------------------                         |'
      end if
      write (stdout, '(1x,a)') '|            No.           x           y           z                         |'
      write (stdout, '(1x,a)') '|            ---        --------------------------------                     |'
      do i = 1, kmesh_info%nnh
        write (stdout, '(1x,"|",11x,i3,5x,3f12.6,21x,"|")') i, (kmesh_info%bka(j, i)/print_output%lenconfac, j=1, 3)
      end do
      write (stdout, '(1x,"+",76("-"),"+")')
      write (stdout, *) ' '
    end if

    ! find index array
    do nkp = 1, num_kpts
      do na = 1, kmesh_info%nnh
        ! first, zero the index array so we can check it gets filled
        kmesh_info%neigh(nkp, na) = 0
        ! now search through list of neighbours of this k-point
        do nn = 1, kmesh_info%nntot
          call utility_compar(kmesh_info%bka(1, na), bk_local(1, nn, nkp), ifpos, ifneg)
          if (ifpos .eq. 1) kmesh_info%neigh(nkp, na) = nn
        end do
        ! check found
        if (kmesh_info%neigh(nkp, na) .eq. 0) then
          if (print_output%iprint > 0) write (stdout, *) ' nkp,na=', nkp, na
          call set_error_fatal(error, 'kmesh_get: failed to find neighbours for this kpoint', comm)
          return
        end if
      end do
    end do

    !fill in the global arrays from the local ones

    do loop = 1, kmesh_info%nntot
      kmesh_info%wb(loop) = wb_local(loop)
    end do

    do loop_s = 1, num_kpts
      do loop = 1, kmesh_info%nntot
        kmesh_info%bk(:, loop, loop_s) = bk_local(:, loop, loop_s)
      end do
    end do

![ysl-b]

    if (gamma_only) then
      ! use half of the b-vectors
      if (num_kpts .ne. 1) then
        call set_error_input(error, 'Error in kmesh_get: wrong choice of gamma_only option', comm)
        return
      end if

      ! reassign nnlist, nncell, wb, bk
      allocate (nnlist_tmp(num_kpts, kmesh_info%nntot), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating nnlist_tmp in kmesh_get', comm)
        return
      end if
      allocate (nncell_tmp(3, num_kpts, kmesh_info%nntot), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating nncell_tmp in kmesh_get', comm)
        return
      end if

      nnlist_tmp(:, :) = kmesh_info%nnlist(:, :)
      nncell_tmp(:, :, :) = kmesh_info%nncell(:, :, :)

      deallocate (kmesh_info%nnlist, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating nnlist in kmesh_get', comm)
        return
      end if
      deallocate (kmesh_info%nncell, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating nncell in kmesh_get', comm)
        return
      end if
      deallocate (kmesh_info%wb, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating wb in kmesh_get', comm)
        return
      end if
      deallocate (kmesh_info%bk, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating bk in kmesh_get', comm)
        return
      end if

      kmesh_info%nntot = kmesh_info%nntot/2

      allocate (kmesh_info%nnlist(num_kpts, kmesh_info%nntot), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating nnlist in kmesh_get', comm)
        return
      end if
      allocate (kmesh_info%nncell(3, num_kpts, kmesh_info%nntot), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating nncell in kmesh_get', comm)
        return
      end if
      allocate (kmesh_info%wb(kmesh_info%nntot), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating wb in kmesh_get', comm)
        return
      end if
      allocate (kmesh_info%bk(3, kmesh_info%nntot, num_kpts), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating bk in kmesh_get', comm)
        return
      end if

      na = 0
      do nn = 1, 2*kmesh_info%nntot
        ifound = 0
        if (na .ne. 0) then
          do nap = 1, na
            call utility_compar(kmesh_info%bk(1, nap, 1), bk_local(1, nn, 1), ifpos, ifneg)
            if (ifneg .eq. 1) ifound = 1
          end do
        end if
        if (ifound .eq. 0) then
          !         found new vector to add to set
          na = na + 1
          kmesh_info%bk(1, na, 1) = bk_local(1, nn, 1)
          kmesh_info%bk(2, na, 1) = bk_local(2, nn, 1)
          kmesh_info%bk(3, na, 1) = bk_local(3, nn, 1)
          kmesh_info%wb(na) = 2.0_dp*wb_local(nn)
          kmesh_info%nnlist(1, na) = nnlist_tmp(1, nn)
          kmesh_info%nncell(1, 1, na) = nncell_tmp(1, 1, nn)
          kmesh_info%nncell(2, 1, na) = nncell_tmp(2, 1, nn)
          kmesh_info%nncell(3, 1, na) = nncell_tmp(3, 1, nn)
          kmesh_info%neigh(1, na) = na
          ! check bk.eq.bka
          call utility_compar(kmesh_info%bk(1, na, 1), kmesh_info%bka(1, na), ifpos, ifneg)
          if (ifpos .ne. 1) then
            call set_error_input(error, 'Error in kmesh_get: bk is not identical to bka in gamma_only option', comm)
            return
          end if
        end if
      end do

      if (na .ne. kmesh_info%nnh) then
        call set_error_fatal(error, 'kmesh_get: Did not find right number of b-vectors in gamma_only option', comm)
        return
      end if

      if (print_output%iprint > 0) then
        write (stdout, '(1x,"+",76("-"),"+")')
        write (stdout, '(1x,a)') '|        Gamma-point: number of the b-vectors is reduced by half             |'
        write (stdout, '(1x,"+",76("-"),"+")')
        if (trim(print_output%length_unit) == 'Ang') then
          write (stdout, '(1x,a)') '|                  b_k Vectors (Ang^-1) and Weights (Ang^2)                  |'
          write (stdout, '(1x,a)') '|                  ----------------------------------------                  |'
        else
          write (stdout, '(1x,a)') '|                 b_k Vectors (Bohr^-1) and Weights (Bohr^2)                 |'
          write (stdout, '(1x,a)') '|                 ------------------------------------------                 |'
        end if
        write (stdout, '(1x,a)') '|            No.         b_k(x)      b_k(y)      b_k(z)        w_b           |'
        write (stdout, '(1x,a)') '|            ---        --------------------------------     --------        |'
        do i = 1, kmesh_info%nntot
          write (stdout, '(1x,"|",11x,i3,5x,3f12.6,3x,f10.6,8x,"|")') &
            i, (kmesh_info%bk(j, i, 1)/print_output%lenconfac, j=1, 3), kmesh_info%wb(i)*print_output%lenconfac**2
        end do
        write (stdout, '(1x,"+",76("-"),"+")')
        write (stdout, *) ' '
      end if

      deallocate (nnlist_tmp, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating nnlist_tmp in kmesh_get', comm)
        return
      end if
      deallocate (nncell_tmp, stat=ierr)
      if (ierr /= 0) then
        call set_error_dealloc(error, 'Error in deallocating nncell_tmp in kmesh_get', comm)
        return
      end if

    end if
![ysl-e]

    ! JJ, is use_ss_functional necessarily defined here, or must it be moved to "special"
    !if (wann_control%use_ss_functional) then
    if (.not. gamma_only) then
      allocate (kmesh_info%nnord(kmesh_info%nntot, num_kpts), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating kmesh_info%nnord in kmesh_get', comm)
        return
      end if
      allocate (kmesh_info%nninv(kmesh_info%nntot, num_kpts), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating kmesh_info%nninv in kmesh_get', comm)
        return
      end if
      allocate (kmesh_info%nnrev(kmesh_info%nntot, num_kpts), stat=ierr)
      if (ierr /= 0) then
        call set_error_alloc(error, 'Error in allocating kmesh_info%nnrev in kmesh_get', comm)
        return
      end if
      call kmesh_bvectors_perm(kmesh_info%bk(:, :, :), kmesh_info%bk(:, :, 1), num_kpts, &
                               kmesh_info%nntot, kmesh_info%nnord, kmesh_info%nninv, &
                               kmesh_info%nnrev, error, comm)
      if (allocated(error)) return
    end if

    deallocate (kpt_cart, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error deallocating kpt_cart in kmesh_get', comm)
      return
    end if

    deallocate (bk_local, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating bk_local in kmesh_get', comm)
      return
    end if

    deallocate (lmn, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating lmn in kmesh_get', comm)
      return
    end if

    deallocate (nnshell, stat=ierr)
    if (ierr /= 0) then
      call set_error_dealloc(error, 'Error in deallocating nnshell in kmesh_get', comm)
      return
    end if

    if (print_output%timing_level > 0) call io_stopwatch_stop('kmesh: get', timer)

    return

  end subroutine kmesh_get