function get_spin_square(dmat_a,dmat_b,ta,tb,abxc,Smat,nocb) result(s2)
! dmat_a / dmat_b -- alpha/beta density of the excited state
! ta / tb -- alpha/beta difference density matrix
use precision, only : dp
use mathlib, only: symmetrize_matrix, traceprod_sym_packed
use mathlib, only: pack_matrix, unpack_matrix
use messages, only: show_message, with_abort
implicit none
real(kind=dp), intent(in), dimension(:) :: &
dmat_a, dmat_b, ta, tb
real(kind=dp), intent(in), dimension(:,:) :: abxc
real(kind=dp), intent(in), dimension(:) :: smat
integer, intent(in) :: nocb
real(kind=dp) :: s2
real(kind=dp), allocatable :: scr1(:), dmat_t(:), &
dmat_t_sq(:,:), smat_sq(:,:), tmp1(:,:), tmp2(:,:)
integer :: nbf, nbf_tri, ok
real(kind=dp) :: dum1, dum2, dum3, dum4
nbf = ubound(abxc, 1)
nbf_tri = ubound(dmat_a, 1)
allocate(scr1(nbf_tri), &
dmat_t(nbf_tri), &
dmat_t_sq(nbf,nbf), &
smat_sq(nbf,nbf), &
tmp1(nbf,nbf), &
tmp2(nbf,nbf), &
source=0.0_dp, stat=ok)
if (ok/=0) call show_message('Cannot allocate memory in qet_spin_square',with_abort)
! Calculate spin expectation values
dum1 = nocb+1
! Symmetric matrix scr1 = Smat*Dmat_a*Smat
dmat_t = dmat_a + ta
call unpack_matrix(dmat_t,dmat_t_sq)
call unpack_matrix(smat,smat_sq)
call dgemm('n', 'n', nbf, nbf, nbf, &
1.0_dp, smat_sq, nbf, &
dmat_t_sq, nbf, &
0.0_dp, tmp1, nbf)
call dgemm('n', 'n', nbf, nbf, nbf, &
1.0_dp, tmp1, nbf, &
smat_sq, nbf, &
0.0_dp, tmp2, nbf)
call pack_matrix(tmp2,scr1)
! -tr[ Dmat_b*Smat*Dmat_a*Smat ]
dmat_t = dmat_b + tb
dum2 = -traceprod_sym_packed(dmat_t,scr1,nbf)
! Symmetric matrix scr1 = Smat*Ta*Smat
call unpack_matrix(Ta,tmp1)
call dgemm('n', 'n', nbf, nbf, nbf, &
1.0_dp, smat_sq, nbf, &
tmp1, nbf, &
0.0_dp, tmp2, nbf)
call dgemm('n', 'n', nbf, nbf, nbf, &
1.0_dp, tmp2, nbf, &
smat_sq, nbf, &
0.0_dp, tmp1, nbf)
call pack_matrix(tmp1,scr1)
! -tr[ Tb*Smat*Ta*Smat ])
dum3 =-traceprod_sym_packed(tb,scr1,nbf)
! +tr[ abxc*Smat ]
tmp1 = abxc
call symmetrize_matrix(tmp1, nbf)
call pack_matrix(tmp1, scr1)
dum4 = traceprod_sym_packed(scr1, smat, nbf)/2.0_dp
s2 = dum1 + dum2 - dum3 + dum4**2
end function get_spin_square