subroutine spin_get_S(kpt, S, num_wann, ws_region, wannier_data, real_lattice, mp_grid, &
ws_distance, HH_R, SS_R, wigner_seitz, error, comm)
!================================================!
!
! Computes <psi_{nk}^(H)|S|psi_{nk}^(H)> (n=1,...,num_wann)
! where S = (S_x,S_y,S_z) is the vector of Pauli matrices
!
!================================================ !
use w90_constants, only: dp
use w90_utility, only: utility_diagonalize, utility_rotate_diag
use w90_types, only: print_output_type, wannier_data_type, ws_region_type, &
ws_distance_type
use w90_postw90_common, only: pw90common_fourier_R_to_k
use w90_postw90_types, only: wigner_seitz_type
use w90_comms, only: w90_comm_type
! arguments
type(ws_region_type), intent(in) :: ws_region
type(wannier_data_type), intent(in) :: wannier_data
type(wigner_seitz_type), intent(in) :: wigner_seitz
type(ws_distance_type), intent(inout) :: ws_distance
type(w90_error_type), allocatable, intent(out) :: error
type(w90_comm_type), intent(in) :: comm
integer, intent(in) :: mp_grid(3)
integer, intent(in) :: num_wann
real(kind=dp), intent(in) :: kpt(3)
real(kind=dp), intent(in) :: real_lattice(3, 3)
real(kind=dp), intent(out) :: S(num_wann, 3)
complex(kind=dp), allocatable, intent(inout) :: HH_R(:, :, :) ! <0n|r|Rm>
complex(kind=dp), allocatable, intent(inout) :: SS_R(:, :, :, :) ! <0n|sigma_x,y,z|Rm>
! local variables
! Physics
complex(kind=dp), allocatable :: HH(:, :)
complex(kind=dp), allocatable :: UU(:, :)
complex(kind=dp), allocatable :: SS(:, :, :)
real(kind=dp) :: eig(num_wann)
! Misc/Dummy
integer :: i
allocate (HH(num_wann, num_wann))
allocate (UU(num_wann, num_wann))
allocate (SS(num_wann, num_wann, 3))
call pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, HH, HH_R, &
kpt, real_lattice, mp_grid, 0, num_wann, error, comm)
if (allocated(error)) return
call utility_diagonalize(HH, num_wann, eig, UU, error, comm)
if (allocated(error)) return
do i = 1, 3
call pw90common_fourier_R_to_k(ws_region, wannier_data, ws_distance, wigner_seitz, &
SS(:, :, i), SS_R(:, :, :, i), kpt, real_lattice, mp_grid, &
0, num_wann, error, comm)
if (allocated(error)) return
S(:, i) = real(utility_rotate_diag(SS(:, :, i), UU, num_wann), dp)
end do
end subroutine spin_get_S