R_wz_sc Subroutine

  subroutine R_wz_sc(R_in, R0, ndeg, R_out, shifts, mp_grid, real_lattice, inv_lattice, &
                     ws_search_size, ws_distance_tol, error, comm)
    !================================================!
    !! Put R_in in the Wigner-Seitz cell centered around R0,
    !! and find all equivalent vectors to this (i.e., with same distance).
    !! Return their coordinates and the degeneracy, as well as the integer
    !! shifts needed to get the vector (these are always multiples of
    !! the mp_grid, i.e. they are supercell displacements in the large supercell)
    !================================================!

    use w90_utility, only: utility_cart_to_frac, utility_frac_to_cart

    implicit none

    ! arguments
    integer, intent(in) :: mp_grid(3)
    integer, intent(in) :: ws_search_size(3)
    real(kind=dp), intent(in) :: real_lattice(3, 3)
    real(kind=dp), intent(in) :: inv_lattice(3, 3)
    real(kind=dp), intent(in) :: ws_distance_tol
    real(DP), intent(in) :: R_in(3)
    real(DP), intent(in) :: R0(3)
    integer, intent(out) :: ndeg
    real(DP), intent(out) :: R_out(3, ndegenx)
    integer, intent(out) :: shifts(3, ndegenx)
    type(w90_error_type), allocatable, intent(out) :: error
    type(w90_comm_type), intent(in) :: comm

    ! local variables
    real(DP) :: R(3), R_f(3), R_in_f(3), R_bz(3), mod2_R_bz
    integer :: i, j, k

    ! init
    ndeg = 0
    R_out = 0._dp
    shifts = 0
    R_bz = R_in
    mod2_R_bz = SUM((R_bz - R0)**2)
    !
    ! take R_bz to cryst(frac) coord for translating
    call utility_cart_to_frac(R_bz, R_in_f, inv_lattice)

    ! In this first loop, I just look for the shortest vector that I obtain
    ! by trying to displace the second Wannier function by all
    ! 'large-supercell' vectors
    ! The size of the supercell, controlled by ws_search_size,
    ! is incremented by one unit in order to account for WFs whose centre
    ! wanders away from the original reference unit cell
    do i = -ws_search_size(1) - 1, ws_search_size(1) + 1
      do j = -ws_search_size(2) - 1, ws_search_size(2) + 1
        do k = -ws_search_size(3) - 1, ws_search_size(3) + 1

          R_f = R_in_f + REAL((/i*mp_grid(1), j*mp_grid(2), k*mp_grid(3)/), &
                              kind=DP)
          call utility_frac_to_cart(R_f, R, real_lattice)

          if (SUM((R - R0)**2) < mod2_R_bz) then
            R_bz = R
            mod2_R_bz = SUM((R_bz - R0)**2)
            ! I start to set a first shift that is applied to get R_bz.
            ! Note: I reset these every time I find a smaller vector.
            !
            ! At this stage, this is the same for all potentially degenerate
            ! points (hence the use of : in shifts(1,:), for instance)
            ! In the second loop below, this shift will be added to the
            ! additional shift that differs for each degenerate but
            ! equivalent point
            shifts(1, :) = i*mp_grid(1)
            shifts(2, :) = j*mp_grid(2)
            shifts(3, :) = k*mp_grid(3)
          end if
        end do
      end do
    end do

    ! Now, second loop to find the list of R_out that differ from R_in
    ! by a large-supercell lattice vector and are equally distant from R0
    ! (i.e. that are on the edges of the WS cell centered on R0)
    ! As above, the size of the supercell, controlled by ws_search_size,
    ! is incremented by one unit in order to account for WFs whose centre
    ! wanders away from the original reference unit cell

    ! I start from the last R_bz found
    mod2_R_bz = SUM((R_bz - R0)**2)
    ! check if R0 and R_in are the same vector
    if (mod2_R_bz < ws_distance_tol**2) then
      ndeg = 1
      R_out(:, 1) = R0
      ! I can safely return as 'shifts' is already set
      return
    end if
    !
    ! take R_bz to cryst(frac) coord for translating
    call utility_cart_to_frac(R_bz, R_in_f, inv_lattice)

    do i = -ws_search_size(1) - 1, ws_search_size(1) + 1
      do j = -ws_search_size(2) - 1, ws_search_size(2) + 1
        do k = -ws_search_size(3) - 1, ws_search_size(3) + 1

          r_f = r_in_f + real((/i*mp_grid(1), j*mp_grid(2), k*mp_grid(3)/), &
                              kind=DP)
          call utility_frac_to_cart(R_f, R, real_lattice)

          if (abs(sqrt(sum((r - r0)**2)) - sqrt(mod2_r_bz)) < ws_distance_tol) then
            ndeg = ndeg + 1
            if (ndeg > ndegenx) then
              call set_error_fatal(error, "surprising ndeg, I wouldn't expect a degeneracy larger than 8...", comm)
              return
            end if
            R_out(:, ndeg) = R
            ! I return/update also the shifts. Note that I have to sum these
            ! to the previous value since in this second loop I am using
            ! R_bz (from the first loop) as the 'central' reference point,
            ! that is already shifted by shift(:,ndeg)
            shifts(1, ndeg) = shifts(1, ndeg) + i*mp_grid(1)
            shifts(2, ndeg) = shifts(2, ndeg) + j*mp_grid(2)
            shifts(3, ndeg) = shifts(3, ndeg) + k*mp_grid(3)
          end if

        end do
      end do
    end do
    !================================================!
  end subroutine R_wz_sc