1subroutine print_spin(qspin)
2  !
3  ! Prints spin in output and CML files
4  !
5  use m_spin,          only: spin
6  use atomlist,        only: no_l
7  use sparse_matrices, only: listhptr, numh
8  use sparse_matrices, only: S, Dscf   ! Dscf could have 1, 2, 4, or 8 components
9  use siesta_cml
10  use parallel,        only: IOnode
11  use precision,       only: dp
12#ifdef MPI
13  use m_mpi_utils,     only: globalize_sum
14#endif
15
16  implicit none
17
18  real(dp), intent(out)  :: qspin(spin%Grid)
19
20  integer  :: ispin, io, j, ind
21  real(dp) :: qaux
22  real(dp) :: Stot        ! Total spin magnitude
23  real(dp) :: Svec(3)     ! Total spin vector
24#ifdef MPI
25  real(dp) :: qtmp(spin%Grid)
26#endif
27
28  qspin(:) = 0.0_dp
29
30  if ( spin%Grid < 2 ) return
31
32  if ( spin%SO ) then
33
34    do io = 1,no_l
35      do j = 1,numh(io)
36        ind = listhptr(io)+j
37        ! In the SOC case, hermitify Dscf
38        qspin(1:2) = qspin(1:2) + dscf(ind,1:2) * S(ind)
39        qspin(3) = qspin(3) + 0.5_dp*(dscf(ind,3)+dscf(ind,7)) * S(ind)
40        qspin(4) = qspin(4) + 0.5_dp*(dscf(ind,4)+dscf(ind,8)) * S(ind)
41      end do
42    end do
43
44  else
45
46    do io = 1,no_l
47      do j = 1,numh(io)
48        ind = listhptr(io)+j
49        qspin(:) = qspin(:) + Dscf(ind,:) * S(ind)
50      end do
51    end do
52
53  endif
54#ifdef MPI
55  ! Global reduction of spin components
56  call globalize_sum(qspin(1:spin%Grid),qtmp(1:spin%Grid))
57  qspin(1:spin%Grid) = qtmp(1:spin%Grid)
58#endif
59
60  ! We are left with printing out to stdout
61  if ( .not. IONode ) return
62
63  if ( spin%Grid == 2) then
64
65    Svec(1) = 0.0_dp
66    Svec(2) = 0.0_dp
67    Svec(3) = qspin(1) - qspin(2)
68    Stot = Svec(3)
69    write(6,'(5x,a,f10.5,2f10.1,f10.5)') 'spin moment: S , {S} = ', Stot, Svec
70    if (cml_p) call cmlAddProperty(xf=mainXML,            &
71        value=qspin(1)-qspin(2), dictref='siesta:stot', &
72        units='siestaUnits:spin')
73
74  else if ( spin%Grid == 4 ) then
75
76    call spnvec( spin%Grid, qspin, qaux, Stot, Svec )
77    write(6,'(5x,a,4f10.5)') 'spin moment: S , {S} = ', Stot, Svec
78    if (cml_p) then
79      call cmlAddProperty(xf=mainXML, value=Stot,  &
80          dictref='siesta:stot', units='siestaUnits:spin')
81      call cmlAddProperty(xf=mainXML, value=Svec,  &
82          dictref='siesta:svec', units='siestaUnits:spin')
83    end if !cml_p
84
85  end if
86
87end subroutine print_spin
88