ws_apply_ndegen Subroutine

  subroutine ws_apply_ndegen(ws_distance, use_ws_distance, num_wann, nrpts, ndegen, &
                             nrpts_full, ir_map, op_R, op_R_full)
    !================================================!
    !! Divide a real-space operator by its degeneracy weights and spread it over
    !! the expanded lattice-vector list built by ws_expand_rvec, so that it can be
    !! Fourier transformed with a plain sum over exp(i k.R), irrespective of
    !! use_ws_distance.
    !================================================!

    use w90_constants, only: cmplx_0
    use w90_types, only: ws_distance_type

    implicit none

    type(ws_distance_type), intent(in) :: ws_distance

    logical, intent(in) :: use_ws_distance
    integer, intent(in) :: num_wann
    integer, intent(in) :: nrpts
    integer, intent(in) :: ndegen(nrpts)
    integer, intent(in) :: nrpts_full
    integer, intent(in) :: ir_map(:, :, :, :)

    complex(kind=dp), intent(in) :: op_R(num_wann, num_wann, nrpts)
    !! operator on the folded grid, before applying the degeneracy weights
    complex(kind=dp), intent(out) :: op_R_full(num_wann, num_wann, nrpts_full)
    !! operator on the expanded grid, after applying the degeneracy weights

    integer :: ir, jr, i, j, ideg

    if (use_ws_distance) then
      op_R_full = cmplx_0
      do ir = 1, nrpts
        do j = 1, num_wann
          do i = 1, num_wann
            do ideg = 1, ws_distance%ndeg(i, j, ir)
              jr = ir_map(ideg, i, j, ir)
              op_R_full(i, j, jr) = op_R_full(i, j, jr) &
                                    + op_R(i, j, ir)/real(ndegen(ir)*ws_distance%ndeg(i, j, ir), dp)
            end do
          end do
        end do
      end do
    else
      ! nrpts_full == nrpts in this case
      do ir = 1, nrpts
        op_R_full(:, :, ir) = op_R(:, :, ir)/real(ndegen(ir), dp)
      end do
    end if
    !================================================!
  end subroutine ws_apply_ndegen