subroutine annihilate_rdm_list(rdm)
! Perform annihilation on RDM elements in rdm%elements.
! * WARNING * rdm%hash_table is not used by this routine, and will not
! be updated by it after reordering and compressing rdm%elements.
use sort_mod, only: sort
type(rdm_list_t), intent(inout) :: rdm
integer :: ielem, counter
real(dp) :: rdm_sign(rdm%sign_length), summed_rdm_sign(rdm%sign_length)
if (rdm%nelements > 0) then
call sort(rdm%elements(:, 1:rdm%nelements))
summed_rdm_sign = 0.0_dp
counter = 0
do ielem = 1, rdm%nelements - 1
call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)
summed_rdm_sign = summed_rdm_sign + rdm_sign
! Is the next RDM element the same as this one? If so then
! don't keep this element yet, but wait until all signs on
! this RDM element have been summed (annihilated).
if (.not. (rdm%elements(0, ielem) == rdm%elements(0, ielem + 1))) then
! If this element is zero for all RDMs, then don't keep it.
if (any(abs(summed_rdm_sign) > 1.0e-12_dp)) then
counter = counter + 1
rdm%elements(0, counter) = rdm%elements(0, ielem)
call encode_sign_rdm(rdm%elements(:, counter), summed_rdm_sign)
end if
summed_rdm_sign = 0.0_dp
end if
end do
! Account for the final element separately.
call extract_sign_rdm(rdm%elements(:, rdm%nelements), rdm_sign)
summed_rdm_sign = summed_rdm_sign + rdm_sign
if (any(abs(summed_rdm_sign) > 1.0 - 12_dp)) then
counter = counter + 1
rdm%elements(0, counter) = rdm%elements(0, rdm%nelements)
call encode_sign_rdm(rdm%elements(:, counter), summed_rdm_sign)
end if
rdm%nelements = counter
end if
end subroutine annihilate_rdm_list