subroutine calc_en_pert_energy(en_pert, rdm_energy_num, rdm_norm, &
energy_pert)
! Calculate the Epstein-Nesbet perturbation energy.
! This is defined as
!
! E_{PN} = \sum_{a} \frac{ ( \sum_i H_{ai} \psi_i )^2 }{ E_{RDM} - E_{aa} }.
!
! The values ( \sum_i H_{ai} \psi_i )^2 are what is stored as the signs
! in the en_pert%dets objects.
use bit_rep_data, only: nifd, NIfTot
use bit_reps, only: decode_bit_det
use determinants, only: get_helement
use FciMCData, only: Hii
use hphf_integrals, only: hphf_diag_helement
use rdm_data, only: en_pert_t
use rdm_data_utils, only: extract_sign_EN
use SystemData, only: nel, tHPHF
type(en_pert_t), intent(in) :: en_pert
real(dp), intent(in) :: rdm_energy_num(en_pert%sign_length)
real(dp), intent(in) :: rdm_norm(en_pert%sign_length)
real(dp), intent(out) :: energy_pert(en_pert%sign_length)
integer(n_int) :: ilut(0:NIfTot)
integer :: nI(nel)
integer :: idet, istate
real(dp) :: h_aa
real(dp) :: contrib(en_pert%sign_length)
real(dp) :: contrib_rdm(en_pert%sign_length)
energy_pert = 0.0_dp
ilut = 0_n_int
! Loop over all determinants.
do idet = 1, en_pert%ndets
ilut(0:nifd) = en_pert%dets(0:nifd, idet)
call decode_bit_det(nI, ilut)
if (tHPHF) then
h_aa = hphf_diag_helement(nI, ilut)
else
h_aa = get_helement(nI, nI, 0)
end if
call extract_sign_EN(en_pert%sign_length, en_pert%dets(:, idet), contrib)
do istate = 1, en_pert%sign_length
contrib_rdm(istate) = contrib(istate) / ((rdm_energy_num(istate) / rdm_norm(istate)) - h_aa)
end do
energy_pert = energy_pert + contrib_rdm
end do
end subroutine calc_en_pert_energy