subroutine calc_rdm_energy_guga(rdm, rdm_energy_1, rdm_energy_2)
! the GUGA RDM version of the energy calculator
type(rdm_list_t), intent(in) :: rdm
real(dp), intent(out) :: rdm_energy_1(rdm%sign_length)
real(dp), intent(out) :: rdm_energy_2(rdm%sign_length)
integer(int_rdm) :: ijkl
integer :: ielem, i, j, k, l
real(dp) :: rdm_sign(rdm%sign_length)
rdm_energy_1 = 0.0_dp
rdm_energy_2 = 0.0_dp
! Loop over all elements in the 2-RDM.
do ielem = 1, rdm%nelements
ijkl = rdm%elements(0, ielem)
call extract_sign_rdm(rdm%elements(:, ielem), rdm_sign)
call extract_2_rdm_ind(ijkl, i, j, k, l)
! TODO maybe the indices are not correctly accessed here!
rdm_energy_2 = rdm_energy_2 + rdm_sign * get_umat_el(i, k, j, l) / 2.0_dp
if (i == j) then
rdm_energy_1 = rdm_energy_1 + &
rdm_sign * GetTMatEl(2 * k, 2 * l) / real(nel - 1, dp) / 2.0_dp
end if
if (k == l) then
rdm_energy_1 = rdm_energy_1 + &
rdm_sign * GetTMatEl(2 * i, 2 * j) / real(nel - 1, dp) / 2.0_dp
end if
end do
end subroutine calc_rdm_energy_guga