subroutine calc_rdm_spin(rdm, rdm_norm, rdm_spin)
! Return the (unnormalised) estimate of <S^2> from the instantaneous
! 2-RDM estimates.
use rdm_data, only: rdm_spawn_t
use SystemData, only: nel
use UMatCache, only: spatial
type(rdm_list_t), intent(in) :: rdm
real(dp), intent(in) :: rdm_norm(rdm%sign_length)
real(dp), intent(out) :: rdm_spin(rdm%sign_length)
integer(int_rdm) :: ijkl
integer :: ielem
integer :: ij, kl, i, j, k, l ! spin orbitals
integer :: p, q, r, s ! spatial orbitals
real(dp) :: rdm_sign(rdm%sign_length)
rdm_spin = 0.0_dp
! Loop over all RDM elements.
do ielem = 1, rdm%nelements
ijkl = rdm%elements(0, ielem)
! Obtain spin orbital labels.
call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l)
! Obtain spatial orbital labels.
p = spatial(i); q = spatial(j);
r = spatial(k); s = spatial(l);
! Note to the reader for the following code: if mod(i,2) == 1 then
! i is a beta (b) orbital, if mod(i,2) == 0 then it is an alpha (a)
! obrital.
! The following if-statement allows IJIJ spatial combinations.
if (p == r .and. q == s) then
! If we get to this point then we definitely have a contribution
! to add in, so extract the sign.
call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)
! If all labels have the same spatial part (IIII):
if (p == q) then
if (is_beta(i) .and. is_alpha(j) .and. is_beta(k) .and. is_alpha(l)) then
rdm_spin = rdm_spin - 6.0_dp * rdm_sign
end if
else
! We only get here if the spatial parts obey IJIJ, for I /= J:
! The following if-statement allows the following spin combinations:
! aaaa, bbbb, abab and baba.
if (mod(i, 2) == mod(k, 2) .and. mod(j, 2) == mod(l, 2)) then
if (mod(i, 2) == mod(j, 2)) then
! aaaa and bbbb.
rdm_spin = rdm_spin + 2.0_dp * rdm_sign
else
! abab and baba.
rdm_spin = rdm_spin - 2.0_dp * rdm_sign
end if
else
! We only get here if the spin parts are abba or baab.
rdm_spin = rdm_spin + 4.0_dp * rdm_sign
end if
end if
end if
end do
rdm_spin = rdm_spin + 3.0_dp * real(nel, dp) * rdm_norm
rdm_spin = rdm_spin / 4.0_dp
end subroutine calc_rdm_spin