subroutine calc_rho_ii_and_sum_n(one_rdms, norm_1rdm, SumN_Rho_ii)
! Calculate the rho_ii arrays, which hold the diagonals of the 1-RDMs,
! appropriately normalised and sorted in terms of the energy ordering
! of the orbitals. Then, calculate SumN_Rho_ii, which holds the sum of
! this diagonal for orbitals occupied in the HF determinant.
use FciMCData, only: HFDet_True
use LoggingData, only: tDiagRDM
use rdm_data, only: tOpenShell
use RotateOrbsData, only: SymLabelListInv_rot, NoOrbs
use SystemData, only: BRR, nel
use UMatCache, only: gtID
type(one_rdm_t), intent(inout) :: one_rdms(:)
real(dp), intent(in) :: norm_1rdm(:)
real(dp), intent(out) :: SumN_Rho_ii(:)
integer :: irdm, i, HFDet_ID, BRR_ID
! Want to sum the diagonal elements of the 1-RDM for the HF orbitals.
! Given the HF orbitals, SymLabelListInv_rot tells us their position
! in the 1-RDM.
SumN_Rho_ii = 0.0_dp
do irdm = 1, size(one_rdms)
do i = 1, NoOrbs
! rho_ii is the diagonal elements of the 1-RDM. We want this
! ordered according to the energy of the orbitals. Brr has the
! orbital numbers in order of energy... i.e Brr(2) = the orbital
! index with the second lowest energy. BRR is always in spin
! orbitals. i gives the energy level, BRR gives the orbital,
! SymLabelListInv_rot gives the position of this orbital in
! one_rdm.
associate(ind => SymLabelListInv_rot)
if (tDiagRDM) then
if (tOpenShell) then
one_rdms(irdm)%rho_ii(i) = one_rdms(irdm)%matrix(ind(BRR(i)), ind(BRR(i))) * norm_1rdm(irdm)
else
BRR_ID = gtID(BRR(2 * i))
one_rdms(irdm)%rho_ii(i) = one_rdms(irdm)%matrix(ind(BRR_ID), ind(BRR_ID)) * norm_1rdm(irdm)
end if
end if
if (i <= nel) then
if (tOpenShell) then
SumN_Rho_ii(irdm) = SumN_Rho_ii(irdm) + &
(one_rdms(irdm)%matrix(ind(HFDet_True(i)), ind(HFDet_True(i))) * norm_1rdm(irdm))
else
HFDet_ID = gtID(HFDet_True(i))
SumN_Rho_ii(irdm) = SumN_Rho_ii(irdm) + &
(one_rdms(irdm)%matrix(ind(HFDet_ID), ind(HFDet_ID)) * norm_1rdm(irdm)) / 2.0_dp
end if
end if
end associate
end do
end do
end subroutine calc_rho_ii_and_sum_n