subroutine calc_rdm_energy(rdm, rdm_energy_1, rdm_energy_2)
! Calculate both the 1- and 2-electron contributions of the
! (unnormalised) energy from the 2-RDM object in rdm, and output them
! to rdm_energy_1 and rdm_energy_2.
use rdm_data, only: rdm_list_t
use rdm_integral_fns, only: one_elec_int, two_elec_int
use SystemData, only: nel
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, ij, kl, 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)
! Decode pqrs label into p, q, r and s labels.
call calc_separate_rdm_labels(ijkl, ij, kl, i, j, k, l)
! The 2-RDM contribution to the energy:
rdm_energy_2 = rdm_energy_2 + rdm_sign * two_elec_int(i, j, k, l)
! The 1-RDM contribution to the energy:
if (i == k) rdm_energy_1 = rdm_energy_1 + rdm_sign * one_elec_int(j, l) / (nel - 1)
if (j == l) rdm_energy_1 = rdm_energy_1 + rdm_sign * one_elec_int(i, k) / (nel - 1)
if (i == l) rdm_energy_1 = rdm_energy_1 - rdm_sign * one_elec_int(j, k) / (nel - 1)
if (j == k) rdm_energy_1 = rdm_energy_1 - rdm_sign * one_elec_int(i, l) / (nel - 1)
end do
end subroutine calc_rdm_energy