subroutine calc_rdm_trace(rdm, rdm_trace)
! Calculate trace of the 2-RDM in the rdm object, and output it to
! rdm_trace.
! This trace is defined as
!
! rdm_trace = \sum_{ij} \Gamma_{ij,ij},
!
! where \Gamma_{ij,kl} is the 2-RDM stored in rdm, and i and j are
! spin orbital labels.
use rdm_data, only: rdm_spawn_t
type(rdm_list_t), intent(in) :: rdm
real(dp), intent(out) :: rdm_trace(rdm%sign_length)
integer(int_rdm) :: ijkl
integer :: ielem
integer :: ij, kl, i, j, k, l ! spin orbitals
integer(int_rdm) :: ij_, kl_
real(dp) :: rdm_sign(rdm%sign_length)
rdm_trace = 0.0_dp
! Loop over all RDM elements.
do ielem = 1, rdm%nelements
ijkl = rdm%elements(0, ielem)
if (tGUGA) then
call extract_2_rdm_ind(ijkl, i, j, k, l)
if ((i == j .and. k == l)) then! .or. (i == l .and. j == k)) then
call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)
rdm_trace = rdm_trace + rdm_sign
end if
else
call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l)
! If this is a diagonal element, add the element to the trace.
if (ij == kl) then
call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)
rdm_trace = rdm_trace + rdm_sign
end if
end if
end do
end subroutine calc_rdm_trace