subroutine gyrotropic_get_NOA_Bnl_spin(S_h, num_occ, occ_list, num_unocc, unocc_list, Bnl)
!================================================!
!
! Calculating the matrix
! B_{nl,ac}^spin(num_occ,num_unocc,3,3)=
! -i eps_{abc} < u_n | sigma_b | u_l >
! ( dimensionless )
!================================================!
use w90_constants, only: dp, cmplx_i, cmplx_0
implicit none
! arguments
integer, intent(in) :: num_occ, num_unocc
integer, intent(in) :: occ_list(:), unocc_list(:)
complex(kind=dp), intent(in) :: S_h(:, :, :) ! n,l,a
complex(kind=dp), intent(out) :: Bnl(:, :, :, :) ! n,l,a,c
! local variables
integer n, l, a, b, c, n1, l1
Bnl(:, :, :, :) = cmplx_0
do b = 1, 3
c = alpha_A(b)
a = beta_A(b)
do n1 = 1, num_occ
n = occ_list(n1)
do l1 = 1, num_unocc
l = unocc_list(l1)
Bnl(n1, l1, a, c) = S_h(n, l, b)
end do !l1
end do !n1
end do !b
Bnl = Bnl*(-cmplx_i)
end subroutine gyrotropic_get_NOA_Bnl_spin