subroutine extract_bit_rep_avsign_norm(rdm_defs, iLutnI, j, nI, SignI, FlagsI, IterRDMStartI, AvSignI, store)
! The following extract_bit_rep_avsign routine extracts the bit
! representation of the current determinant, and calculates the average
! sign since this determinant became occupied.
! In double run, we have to be particularly careful -- we need to start
! a new average when the determinant becomes newly occupied or
! unoccupied in either population (see CMO thesis). Additionally, we're
! also setting it up so that averages get restarted whenever we
! calculate the energy which saves a lot of faffing about, and storage
! of an extra set of RDMs, and is still unbiased. This is called for
! each determinant in the occupied list at the beginning of its FCIQMC
! cycle. It is used if we're calculating the RDMs with or without HPHF.
! Input: iLutnI (bit rep of current determinant).
! j - Which element in the CurrentDets array are we considering?
! Output: nI, SignI, FlagsI after extract.
! IterRDMStartI - new iteration the determinant became occupied (as a real).
! AvSignI - the new average walker population during this time (also real).
use FciMCData, only: PreviousCycles, Iter, IterRDMStart, excit_gen_store_type
use global_det_data, only: get_iter_occ_tot, get_av_sgn_tot
use global_det_data, only: len_av_sgn_tot, len_iter_occ_tot
use rdm_data, only: rdm_definitions_t
use LoggingData, only: RDMEnergyIter
type(rdm_definitions_t), intent(in) :: rdm_defs
integer(n_int), intent(in) :: iLutnI(0:nIfTot)
integer, intent(out) :: nI(nel), FlagsI
integer, intent(in) :: j
real(dp), intent(out) :: SignI(lenof_sign)
real(dp), intent(out) :: IterRDMStartI(len_iter_occ_tot), AvSignI(len_av_sgn_tot)
type(excit_gen_store_type), intent(inout), optional :: store
integer :: irdm
integer :: av_ind_1, av_ind_2
unused_var(store)
! This is the iteration from which this determinant has been occupied.
IterRDMStartI = get_iter_occ_tot(j)
! This extracts everything.
call extract_bit_rep(iLutnI, nI, SignI, FlagsI, j, store)
associate(ind => rdm_defs%sim_labels)
if (((Iter + PreviousCycles - IterRDMStart) > 0) .and. &
& (mod(((Iter - 1) + PreviousCycles - IterRDMStart + 1), &
RDMEnergyIter) == 0)) then
! The previous iteration was one where we added in diagonal elements
! To keep things unbiased, we need to set up a new averaging block now.
! NB: if doing single run cutoff, note that doing things this way is now
! NOT the same as the technique described in CMO (and DMC's) thesis.
! Would expect diagonal elements to be slightly worse quality, improving
! as one calculates the RDM energy less frequently. As this method is
! biased anyway, I'm not going to lose sleep over it.
do irdm = 1, rdm_defs%nrdms
av_ind_1 = irdm * 2 - 1
av_ind_2 = irdm * 2
AvSignI(av_ind_1) = SignI(ind(1, irdm))
IterRDMStartI(av_ind_1) = real(Iter + PreviousCycles, dp)
AvSignI(av_ind_2) = SignI(ind(2, irdm))
IterRDMStartI(av_ind_2) = real(Iter + PreviousCycles, dp)
end do
else
! Now let's consider other instances in which we need to start a new block:
do irdm = 1, rdm_defs%nrdms
! The indicies of the first and second replicas in this
! particular pair, in the *average* sign arrays.
av_ind_1 = irdm * 2 - 1
av_ind_2 = irdm * 2
if (near_zero(SignI(ind(1, irdm))) .and. (.not. near_zero(IterRDMStartI(av_ind_1)))) then
! The population has just gone to zero on population 1.
! Therefore, we need to start a new averaging block.
AvSignI(av_ind_1) = 0
IterRDMStartI(av_ind_1) = 0
AvSignI(av_ind_2) = SignI(ind(2, irdm))
IterRDMStartI(av_ind_2) = real(Iter + PreviousCycles, dp)
else if (near_zero(SignI(ind(2, irdm))) .and. (.not. near_zero(IterRDMStartI(av_ind_2)))) then
! The population has just gone to zero on population 2.
! Therefore, we need to start a new averaging block.
AvSignI(av_ind_2) = 0
IterRDMStartI(av_ind_2) = 0
AvSignI(av_ind_1) = SignI(ind(1, irdm))
IterRDMStartI(av_ind_1) = real(Iter + PreviousCycles, dp)
else if (.not. near_zero(SignI(ind(1, irdm))) .and. near_zero(IterRDMStartI(av_ind_1))) then
! Population 1 has just become occupied.
IterRDMStartI(av_ind_1) = real(Iter + PreviousCycles, dp)
IterRDMStartI(av_ind_2) = real(Iter + PreviousCycles, dp)
AvSignI(av_ind_1) = SignI(ind(1, irdm))
AvSignI(av_ind_2) = SignI(ind(2, irdm))
if (near_zero(SignI(ind(2, irdm)))) IterRDMStartI(av_ind_2) = 0
else if (.not. near_zero(SignI(ind(2, irdm))) .and. near_zero(IterRDMStartI(av_ind_2))) then
! Population 2 has just become occupied.
IterRDMStartI(av_ind_1) = real(Iter + PreviousCycles, dp)
IterRDMStartI(av_ind_2) = real(Iter + PreviousCycles, dp)
AvSignI(av_ind_1) = SignI(ind(1, irdm))
AvSignI(av_ind_2) = SignI(ind(2, irdm))
if (near_zero(SignI(ind(1, irdm)))) IterRDMStartI(av_ind_1) = 0
else
! Nothing unusual has happened so update both average
! populations as normal.
AvSignI(av_ind_1) = &
(((real(Iter + PreviousCycles, dp) - IterRDMStartI(av_ind_1)) * get_av_sgn_tot(j, av_ind_1)) &
+ SignI(ind(1, irdm))) / (real(Iter + PreviousCycles, dp) - IterRDMStartI(av_ind_1) + 1.0_dp)
AvSignI(av_ind_2) = &
(((real(Iter + PreviousCycles, dp) - IterRDMStartI(av_ind_2)) * get_av_sgn_tot(j, av_ind_2)) &
+ SignI(ind(2, irdm))) / (real(Iter + PreviousCycles, dp) - IterRDMStartI(av_ind_2) + 1.0_dp)
end if
end do
end if
end associate
end subroutine extract_bit_rep_avsign_norm