subroutine DiDj_Found_FillRDM(rdm_defs, spawn, one_rdms, Spawned_No, iLutJ, real_sign_j_all, &
tNonInits)
! This routine is called when we have found a Di (or multiple Di's)
! spawning onto a Dj with sign /= 0 (i.e. occupied). We then want to
! run through all the Di, Dj pairs and add their coefficients
! (with appropriate de-biasing factors) into the 1 and 2 electron RDM.
use bit_reps, only: decode_bit_det
use FciMCData, only: Spawned_Parents, Spawned_Parents_Index, iLutHF_True
use rdm_data, only: one_rdm_t, rdm_definitions_t
use SystemData, only: nel, tHPHF
use bit_reps, only: all_runs_are_initiator
type(rdm_definitions_t), intent(in) :: rdm_defs
type(rdm_spawn_t), intent(inout) :: spawn
type(one_rdm_t), intent(inout) :: one_rdms(:)
integer, intent(in) :: Spawned_No
integer(n_int), intent(in) :: iLutJ(0:NIfTot)
real(dp), intent(in) :: real_sign_j_all(lenof_sign)
logical, intent(in) :: tNonInits
integer :: i, irdm, rdm_ind, nI(nel), nJ(nel)
real(dp) :: realSignI
real(dp) :: input_sign_i(rdm_defs%nrdms), input_sign_j(rdm_defs%nrdms)
integer :: dest_part_type, source_part_type, run
integer(n_int) :: source_flags
logical :: spawning_from_ket_to_bra, tNonInitParent
integer(int_rdm) :: guga_rdm_ind
real(dp) :: x0, x1
! Spawning from multiple parents, to iLutJ, which has SignJ.
! We are at position Spawned_No in the SpawnedParts array.
! Spawned_Parents_Index(1,Spawned_No) is therefore the start position
! of the list of parents (Di's) which spawned on the Dj in
! SpawnedParts(Spawned_No). There are Spawned_Parents_Index(2,Spawned_No)
! of these parent Di's. Spawned_Parents(0:IlutBitsParent%len_orb,x)
! is the determinant Di, Spawned_Parents(IlutBitsParent%ind_pop,x)
! is the un-biased ci. Spawned_Parents(IlutBitsParent%ind_flag,x)
! is the flag (only initiator i believe. WD) and
! Spawned_Parents(IlutBitsParent%ind_source,x) is the simulation index,
! where the spawn was coming from
! Run through all Di's.
! in the GUGA RDMs it somehow can happen that a all 0 real_sign_j_all
! gets in here.. to abort if that happens
if (all(near_zero(real_sign_j_all))) return
do i = Spawned_Parents_Index(1, Spawned_No), &
Spawned_Parents_Index(1, Spawned_No) + Spawned_Parents_Index(2, Spawned_No) - 1
if (DetBitEQ(iLutHF_True, &
Spawned_Parents(0:IlutBitsParent%len_orb, i), IlutBits%len_orb)) then
! We've already added HF - S, and HF - D symmetrically.
! Any connection with the HF has therefore already been added.
cycle
end if
call decode_bit_det(nI, Spawned_Parents(0:IlutBits%len_orb, i))
call decode_bit_det(nJ, iLutJ)
realSignI = transfer(Spawned_Parents(IlutBits%ind_pop, i), realSignI)
! if the sign is 0, there is nothing to do, i.e. all contributions will be 0
if (abs(realSignI) < eps) cycle
! The original spawning event (and the RealSignI) came from this
! population.
source_part_type = int(Spawned_Parents(IlutBitsParent%ind_source, i))
! WD: due to some weird convention in Annihilation (TODO check that!)
! it can happen that a source_part_type = 0 can end up here..
! so cycle if that happens
if (source_part_type == 0) cycle
! if we only sum in initiator contriubtions, check the flags here
! This block is entered if
! a) tNonInits == .false. -> we are looking at the inits-rdms
! b) tNonInitsForRDMs == .false. -> we are calculating inits-only rdms
! In both cases, an initiator check is required
if (.not. (tNonInits .and. tNonInitsForRDMs)) then
tNonInitParent = .false.
! Only initiators are used
if (.not. (all_runs_are_initiator(ilutJ) .or. &
! or the non-variational inits-only version (only require an init in the ket)
(tNonVariationalRDMs .and. .not. tNonInits))) return
do run = 1, inum_runs
if (.not. btest(Spawned_Parents(IlutBitsParent%ind_flag, i), &
get_initiator_flag_by_run(run))) tNonInitParent = .true.
! if a non-initiator is participating, do not sum in
! that contribution
end do
if (tNonInitParent) cycle
end if
! Loop over all RDMs to which the simulation with label
! source_part_type contributes to.
do irdm = 1, rdm_defs%nrdms_per_sim(source_part_type)
! Get the label of the simulation that is paired with this,
! replica, for this particular RDM.
dest_part_type = rdm_defs%sim_pairs(irdm, source_part_type)
! The label of the RDM that this is contributing to.
rdm_ind = rdm_defs%rdm_labels(irdm, source_part_type)
input_sign_i = 0.0_dp
input_sign_j = 0.0_dp
input_sign_i(rdm_ind) = realSignI
input_sign_j(rdm_ind) = real_sign_j_all(dest_part_type)
! If the two FCIQMC simulations involved are different then we
! need a factor of a half, since spawning can occur from
! from either of the two simulations.
if (dest_part_type /= source_part_type) input_sign_i = 0.5_dp * input_sign_i
! sim_labels(1,irdm) holds the state in the 'bra' of the RDM.
! sim_labels(2,irdm) holds the state in the 'ket' of the RDM.
! If source_part_type is equal to the part type of the ket,
! then we can say that we are spawning from the ket to the bra.
spawning_from_ket_to_bra = source_part_type == rdm_defs%sim_labels(2, rdm_ind)
! If we are spawning from the ket of the RDM, then the RDM
! element to be added in is:
!
! \psi_j^b* \psi_i^a < nJ | ... | nI > (1)
!
! (where a and b are state labels), which is what we want.
!
! If we are spawning from the bra to the ket then we cannot
! just pass input_sign_i, input_sign_j, nI and nJ into the
! following routines in the same order, because then we would
! be adding in the same element as in Eq. (1), but we actually
! want the Hermitian conjugate of that term - <nJ| and |nI>
! will be in the wrong order, and so the RDM indices will
! be calculated the wrong way around, and the element will
! be added in on the 'wrong side' of the RDM's diagonal.
! As such, we instead swap the order of things so that we are
! adding in the term
!
! \psi_i^a* \psi_j^b < nI | ... | nJ > (2)
!
! which is clearly as it should be if the determinant we are
! spawning from (nI) is the bra of the RDM.
if (tGUGA) then
call extract_stochastic_rdm_info(IlutBitsParent, &
Spawned_Parents(:, i), guga_rdm_ind, x0, x1)
end if
if (spawning_from_ket_to_bra) then
if (tHPHF) then
call Fill_Spin_Coupled_RDM(spawn, one_rdms, &
Spawned_Parents(0:IlutBitsParent%len_orb, i), &
iLutJ, nI, nJ, input_sign_i, input_sign_j)
else if (tGUGA) then
call Add_RDM_From_IJ_Pair_GUGA(spawn, one_rdms, &
nI, nJ, input_sign_i, input_sign_j, spawning_from_ket_to_bra, &
rdm_ind_in=guga_rdm_ind, x0=x0, x1=x1)
else
call Add_RDM_From_IJ_Pair(spawn, one_rdms, nI, nJ, &
input_sign_i, input_sign_j)
end if
else
! Spawning from the bra to the ket - swap the order of the
! signs and states in the routine interface.
if (tHPHF) then
call Fill_Spin_Coupled_RDM(spawn, one_rdms, iLutJ, &
Spawned_Parents(0:IlutBitsParent%len_orb, i), &
nJ, nI, input_sign_j, input_sign_i)
else if (tGUGA) then
call Add_RDM_From_IJ_Pair_GUGA(spawn, one_rdms, &
nJ, nI, input_sign_j, input_sign_i, spawning_from_ket_to_bra, &
rdm_ind_in=guga_rdm_ind, x0=x0, x1=x1)
else
call Add_RDM_From_IJ_Pair(spawn, one_rdms, nJ, nI, &
input_sign_j, input_sign_i)
end if
end if
end do
end do
end subroutine DiDj_Found_FillRDM