subroutine fill_RDM_offdiag_deterministic(rdm_defs, spawn, one_rdms)
use bit_rep_data, only: NIfD
use bit_reps, only: decode_bit_det
use DetBitOps, only: get_bit_excitmat, FindBitExcitLevel
use FciMCData, only: Iter, IterRDMStart, PreviousCycles, iLutHF_True, core_run
use core_space_util, only: cs_replicas
use LoggingData, only: RDMExcitLevel, RDMEnergyIter
use Parallel_neci, only: iProcIndex
use rdm_data, only: one_rdm_t, rdm_definitions_t
use rdm_data_utils, only: add_to_rdm_spawn_t
use SystemData, only: nel, tHPHF
type(rdm_definitions_t), intent(in) :: rdm_defs
type(rdm_spawn_t), intent(inout) :: spawn
type(one_rdm_t), intent(inout) :: one_rdms(:)
integer :: i, j, irdm
integer :: SingEx(2, 1), Ex(2, maxExcit)
real(dp) :: AvSignI(spawn%rdm_send%sign_length), AvSignJ(spawn%rdm_send%sign_length)
real(dp) :: full_sign(spawn%rdm_send%sign_length)
logical :: tParity
integer(n_int) :: iLutI(0:niftot), iLutJ(0:niftot)
type(CSF_Info_t) :: csf_i, csf_j
integer :: nI(nel), nJ(nel), IC, n
integer :: IterRDM, connect_elem, num_j
type(ExcitationInformation_t) :: excitInfo
character(*), parameter :: this_routine = "fill_rdm_offdiag_deterministic"
! IterRDM will be the number of iterations that the contributions are
! ech weighted by.
if (mod((iter + PreviousCycles - IterRDMStart + 1), RDMEnergyIter) == 0) then
! This numer of iterations is how regularly the energy is printed
! out.
IterRDM = RDMEnergyIter
else
! This must be the final iteration, as we've got tFill_RDM=.true.
! for an iteration where we wouldn't normally need the energy
IterRDM = mod((Iter + PreviousCycles - IterRDMStart + 1), RDMEnergyIter)
end if
Ex = 0
associate(rep => cs_replicas(core_run))
if (t_full_core_rdms) then
! num_j = rep%determ_sizes(iProcIndex)
num_j = rep%determ_space_size
end if
! Cycle over all core dets on this process.
do i = 1, rep%determ_sizes(iProcIndex)
if (.not. t_full_core_rdms) then
num_j = rep%sparse_core_ham(i)%num_elements - 1
end if
iLutI = rep%core_space(:, rep%determ_displs(iProcIndex) + i)
! Connections to the HF are added in elsewhere, so skip them here.
if (DetBitEq(iLutI, iLutHF_True, nifd)) cycle
do irdm = 1, rdm_defs%nrdms
AvSignI(irdm) = rep%full_determ_vecs_av(rdm_defs%sim_labels(2, irdm), &
rep%determ_displs(iProcIndex) + i)
end do
call decode_bit_det(nI, iLutI)
if (tGUGA) csf_i = CSF_Info_t(ilutI)
! if (tGUGA) call fill_csf_i(ilutI(0:nifd))
do j = 1, num_j
! Running over all non-zero off-diag matrix elements
! Connections to whole space (1 row), excluding diagonal elements
! Note: determ_displs holds sum(determ_sizes(0:proc-1))
! Core space holds all the core determinants on every processor,
! so we need to shuffle up to the range of indices corresponding
! to this proc (using determ_displs) and then select the
! correct one, i.
if (t_full_core_rdms) then
iLutJ = rep%core_space(:, j)
else
iLutJ = rep%core_space(:, rep%core_connections(i)%positions(j))
end if
! Connections to the HF are added in elsewhere, so skip them here.
if (DetBitEq(iLutJ, iLutHF_True, nifd)) cycle
if (DetBitEq(iLutJ, iLutI, nifd)) cycle
if (tGUGA) csf_j = CSF_Info_t(ilutJ)
if (t_full_core_rdms) then
if (.not. tGUGA) then
ic = FindBitExcitLevel(iLutI, ilutJ)
ex(1, 1) = ic
call GetBitExcitation(ilutI, ilutJ, ex, tParity)
end if
do irdm = 1, rdm_defs%nrdms
AvSignJ(irdm) = rep%full_determ_vecs_av(rdm_defs%sim_labels(1, irdm), j)
end do
else
do irdm = 1, rdm_defs%nrdms
AvSignJ(irdm) = rep%full_determ_vecs_av(rdm_defs%sim_labels(1, irdm), &
rep%core_connections(i)%positions(j))
end do
connect_elem = rep%core_connections(i)%elements(j)
IC = abs(connect_elem)
if (sign(1, connect_elem) > 0) then
tParity = .false.
else
tParity = .true.
end if
end if
! if (ic > 2) cycle
if (tParity) then
full_sign = -AvSignI * AvSignJ * IterRDM
else
full_sign = AvSignI * AvSignJ * IterRDM
end if
if (tHPHF) then
call decode_bit_det(nJ, iLutJ)
call Fill_Spin_Coupled_RDM(spawn, one_rdms, iLutI, iLutJ, &
nI, nJ, AvSignI * IterRDM, AvSignJ)
else if (tGUGA) then
call add_rdm_from_ij_pair_guga_exact(spawn, one_rdms, &
ilutI, csf_i, ilutJ, csf_j, AvSignI * IterRDM, AvSignJ)
else
if (IC == 1) then
! Single excitation - contributes to 1- and 2-RDM
! (if calculated).
! Note: get_bit_excitmat may be buggy (DetBitOps),
! but will do for now as we need the Ex...
call get_bit_excitmat(iLutI(0:NIfD), iLutJ(0:NIfD), SingEx, IC)
Ex(:, 1) = SingEx(:, 1)
! No need to explicitly fill symmetrically as we'll
! generate pairs of determinants both ways around using
! the connectivity matrix.
if (RDMExcitLevel /= 1) call fill_spawn_rdm_singles(spawn, nI, Ex, full_sign)
if (RDMExcitLevel == 1) then
call fill_sings_1rdm(one_rdms, Ex, tParity, AvSignI * IterRDM, AvSignJ, .false.)
end if
else if ((IC == 2) .and. (RDMExcitLevel /= 1)) then
! Note: get_bit_excitmat may be buggy (DetBitOps),
! but will do for now as we need the Ex...
call get_bit_excitmat(iLutI(0:NIfD), iLutJ(0:NIfD), Ex, IC)
call add_to_rdm_spawn_t(spawn, Ex(2, 1), Ex(2, 2), Ex(1, 1), Ex(1, 2), full_sign, .false.)
end if
end if
end do
end do
end associate
end subroutine fill_RDM_offdiag_deterministic