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