subroutine PerformFCIMCycPar(iter_data, err)
use mpi
use global_det_data, only: get_iter_occ_tot, get_av_sgn_tot
use global_det_data, only: set_av_sgn_tot, set_iter_occ_tot
use global_det_data, only: len_av_sgn_tot, len_iter_occ_tot
use rdm_data, only: two_rdm_spawn, two_rdm_recv, two_rdm_main, one_rdms
use rdm_data, only: rdm_definitions
use rdm_data_utils, only: communicate_rdm_spawn_t, add_rdm_1_to_rdm_2, clear_rdm_list_t
use symrandexcit_Ex_Mag, only: test_sym_excit_ExMag
! Iteration specific data
type(fcimc_iter_data), intent(inout) :: iter_data
integer, intent(out) :: err
! Now the local, iteration specific, variables
integer :: j, p, error, proc_temp, i, HFPartInd, isym
integer :: DetCurr(nel), nJ(nel), FlagsCurr, parent_flags
real(dp), dimension(lenof_sign) :: SignCurr, child, child_for_stats, SpawnSign
integer(kind=n_int) :: iLutnJ(0:niftot)
integer :: IC, walkExcitLevel, walkExcitLevel_toHF, ex(2, 3), TotWalkersNew, part_type, run
integer(int64) :: tot_parts_tmp(lenof_sign)
logical :: tParity, tSuccess, tCoreDet(inum_runs), tGlobalCoreDet
real(dp) :: prob, HDiagCurr, EnergyCurr, hdiag_bare, TempTotParts, Di_Sign_Temp
real(dp) :: RDMBiasFacCurr
real(dp) :: lstart
real(dp) :: AvSignCurr(len_av_sgn_tot), IterRDMStartCurr(len_iter_occ_tot)
real(dp) :: av_sign(len_av_sgn_tot), iter_occ(len_iter_occ_tot)
HElement_t(dp) :: HDiagTemp, HElGen, HOffDiagCurr
character(*), parameter :: this_routine = 'PerformFCIMCycPar'
HElement_t(dp), dimension(inum_runs) :: delta
integer :: proc, pos, determ_index, irdm
real(dp) :: r, sgn(lenof_sign), prob_extra_walker
integer :: DetHash, FinalVal, clash, PartInd, k, y, MaxIndex
type(ll_node), pointer :: TempNode
integer :: ms, allErr
real(dp) :: precond_fac
! average number of excitations per walker for a given determinant
real(dp) :: AvMCExcitsLoc, scale, max_spawn
HElement_t(dp) :: hdiag_spawn, h_diag_correct
logical :: signChanged, newlyOccupied, flag_mixed
real(dp) :: currArg, spawnArg
integer :: scaleFactor
! how many entries were added to (the end of) CurrentDets in the last iteration
integer, save :: detGrowth = 0
real(dp) :: inst_rdm_occ
call set_timer(Walker_Time, 30)
err = 0
allErr = 0
MaxInitPopPos = 0.0_dp
MaxInitPopNeg = 0.0_dp
HighPopNeg = 1
HighPopPos = 1
FlagsCurr = 0
! Reset iteration variables
! Next free position in newly spawned list.
ValidSpawnedList = InitialSpawnedSlots
FreeSlot(1:iEndFreeSlot) = 0 !Does this cover enough?
iStartFreeSlot = 1
iEndFreeSlot = 0
! Clear the hash table for the spawning array.
if (use_spawn_hash_table) call clear_hash_table(spawn_ht)
! Index for counting deterministic states.
determ_index = 0
call rezero_iter_stats_each_iter(iter_data, rdm_definitions)
! The processor with the HF determinant on it will have to check
! through each determinant until it's found. Once found, tHFFound is
! true, and it no longer needs to be checked.
precond_fac = 1.0_dp
IFDEBUGTHEN(FCIMCDebug, stdout)
write(stdout, "(A)") "Hash Table: "
write(stdout, "(A)") "Position in hash table, Position in CurrentDets"
do j = 1, nWalkerHashes
TempNode => HashIndex(j)
if (TempNode%Ind /= 0) then
write(stdout, '(i9)', advance='no') j
do while (associated(TempNode))
write(stdout, '(i9)', advance='no') TempNode%Ind
TempNode => TempNode%Next
end do
write(stdout, '()', advance='yes')
end if
end do
ENDIFDEBUG
IFDEBUG(FCIMCDebug, 3) write(stdout, "(A,I12)") "Walker list length: ", TotWalkers
IFDEBUG(FCIMCDebug, 3) write(stdout, "(A)") "TW: Walker Det"
! This block decides whether or not to calculate the contribution to the RDMs from
! the diagonal elements (and explicit connections to the HF) for each occupied determinant.
! For efficiency, this is only done on the final iteration, or one where the RDM energy is
! being printed.
tFill_RDM = .false.
if (tFillingStochRDMonFly) then
if (mod((Iter + PreviousCycles - IterRDMStart + 1), RDMEnergyIter) == 0) then
! RDM energy is being printed, calculate the diagonal elements for
! the last RDMEnergyIter iterations.
tFill_RDM = .true.
IterLastRDMFill = RDMEnergyIter
else if(Iter == NMCyc .or. Iter - maxval(VaryShiftIter) == eq_cyc) then
! Last iteration, calculate the diagonal element for the iterations
! since the last time they were included.
tFill_RDM = .true.
IterLastRDMFill = mod((Iter + PreviousCycles - IterRDMStart + 1), RDMEnergyIter)
end if
end if
lstart = mpi_wtime()
loop_over_determinants: do j = 1, int(TotWalkers)
! N.B. j indicates the number of determinants, not the number
! of walkers.
! Indicate that the scratch storage used for excitation generation
! from the same walker has not been filled (it is filled when we
! excite from the first particle on a determinant).
fcimc_excit_gen_store%tFilled = .false.
! let the interested excitation generator know about the index in the CurrentDets array.
fcimc_excit_gen_store%idx_curr_dets = j
! Make sure that the parent flags from the last walker don't through.
parent_flags = 0
! If we're not calculating the RDM (or we're calculating some HFSD combination of the
! RDM) this just extracts info from the bit representation like normal.
! IterRDMStartCurr and AvSignCurr just come out as 1.0_dp.
! Otherwise, it extracts the Curr info, and calculates the iteration this determinant
! became occupied (IterRDMStartCurr) and the average population during that time
! (AvSignCurr).
! Is this state is in the deterministic space?
do run = 1, inum_runs
tCoreDet(run) = check_determ_flag(CurrentDets(:, j), run)
end do
tGlobalCoreDet = all(tCoreDet)
call extract_bit_rep_avsign(rdm_definitions, CurrentDets(:, j), j, DetCurr, SignCurr, FlagsCurr, &
IterRDMStartCurr, AvSignCurr, fcimc_excit_gen_store)
! if we CurrentDets is almost full, make some space
if (t_activate_decay) then
if (test_flag(CurrentDets(:, j), flag_prone)) then
! kill a prone determinant with probability given by the
! ratio of how many space we probably need to how many prone dets exist
! (prone dets always have only a single scaled walker)
r = genrand_real2_dSFMT()
if (n_prone_dets > 0) then
if (r < detGrowth / n_prone_dets) then
! log the removal
iter_data%nremoved = iter_data%nremoved + abs(SignCurr)
! remove all walkers here (this will be counted as unocc. later
! and thus become an empty slot)
SignCurr = 0.0_dp
! kill all walkers on the determinant
call nullify_ilut(CurrentDets(:, j))
! and remove it from the hashtable
if (.not. tAccumEmptyDet(CurrentDets(:, j))) &
call RemoveHashDet(HashIndex, DetCurr, j)
cycle
end if
end if
end if
end if
! We only need to find out if determinant is connected to the
! reference (so no ex. level above 2 required,
! truncated etc.)
walkExcitLevel = FindBitExcitLevel(iLutRef(:, 1), CurrentDets(:, j), &
max_calc_ex_level, .true.)
if (tRef_Not_HF) then
walkExcitLevel_toHF = FindBitExcitLevel(iLutHF_true, CurrentDets(:, j), &
max_calc_ex_level, .true.)
else
walkExcitLevel_toHF = walkExcitLevel
end if
if (tFillingStochRDMonFly) then
! Set the average sign and occupation iteration which were
! found in extract_bit_rep_avsign.
call set_av_sgn_tot(j, AvSignCurr)
call set_iter_occ_tot(j, IterRDMStartCurr)
! If this is an iteration where we print out the RDM energy,
! add in the diagonal contribution to the RDM for this
! determinant, for each rdm.
if (tFill_RDM .and. (.not. tNoNewRDMContrib)) then
av_sign = get_av_sgn_tot(j)
iter_occ = get_iter_occ_tot(j)
if (tInitsRDM .and. all_runs_are_initiator(CurrentDets(:, j))) &
call fill_rdm_diag_currdet(two_rdm_inits_spawn, inits_one_rdms, &
CurrentDets(:, j), DetCurr, walkExcitLevel_toHF, av_sign, iter_occ, &
tGlobalCoreDet, .false.)
call fill_rdm_diag_currdet(two_rdm_spawn, one_rdms, CurrentDets(:, j), &
DetCurr, walkExcitLevel_toHF, av_sign, iter_occ, tGlobalCoreDet, tApplyLC)
end if
end if
! This if-statement is only entered when using semi-stochastic and
! only if this determinant is in the core space.
! Potential optimization: This list does not change between iterations, only set it up once
if (allocated(cs_replicas)) then
do run = 1, size(cs_replicas)
associate(rep => cs_replicas(run))
if (tCoreDet(run)) then
! Store the index of this state, for use in annihilation later.
! A global core-space is ordered in the same fashion
! in currentdets as internally in the core_space
if (t_global_core_space) then
determ_index = determ_index + 1
else
! In general, that is not true however, the
! core-space index has to be looked up
determ_index = core_space_pos(CurrentDets(:, j), DetCurr, &
run) - rep%determ_displs(iProcIndex)
end if
rep%indices_of_determ_states(determ_index) = j
! Add this amplitude to the deterministic vector.
rep%partial_determ_vecs(:, determ_index) = SignCurr(rep%min_part():rep%max_part())
end if
end associate
end do
end if
! As the main list (which is storing a hash table) no longer needs
! to be contiguous, we need to skip sites that are empty.
if (IsUnoccDet(SignCurr)) then
! The deterministic states are always kept in CurrentDets, even when
! the amplitude is zero. Hence we must check if the amplitude is zero,
! and if so, skip the state.
if (any(tCoreDet) .or. tAccumEmptyDet(CurrentDets(:, j))) cycle
!It has been removed from the hash table already
!Just add to the "freeslot" list
iEndFreeSlot = iEndFreeSlot + 1
FreeSlot(iEndFreeSlot) = j
cycle
end if
! sum in (fmu-1)*cmu^2 for the purpose of RDMs
if (tAdaptiveShift .and. all(.not. tSinglePartPhase) .and. tFillingStochRDMOnFly) then
! Only add the contribution from the corespace if it is explicitly demanded
if ((.not. tGlobalCoreDet) .or. tCoreAdaptiveShift) &
call SumCorrectionContrib(SignCurr, j)
end if
! The current diagonal matrix element is stored persistently.
HDiagCurr = det_diagH(j)
HOffDiagCurr = det_offdiagH(j)
EnergyCurr = det_diagH(j) + Hii
if (tSeniorInitiators) then
SpawnSign = get_all_spawn_pops(j)
do run = 1, inum_runs
if (.not. is_run_unnocc(SignCurr, run) .and. .not. is_run_unnocc(SpawnSign, run)) then
#ifdef CMPLX_
!For complex walkers, we consider the sign changed when the argument of the complex
!number changes more than pi/2.
CurrArg = DATAN2(SignCurr(max_part_type(run)), SignCurr(min_part_type(run)))
SpawnArg = DATAN2(SpawnSign(max_part_type(run)), SpawnSign(min_part_type(run)))
signChanged = mod(abs(CurrArg - SpawnArg), PI) > PI / 2.0_dp
#else
signChanged = SpawnSign(min_part_type(run)) * SignCurr(min_part_type(run)) < 0.0_dp
#endif
else
signChanged = .false.
end if
newlyOccupied = is_run_unnocc(SignCurr, run) .and. .not. is_run_unnocc(SpawnSign, run)
if (signChanged .or. newlyOccupied) then
call reset_tau_int(j, run)
call reset_shift_int(j, run)
call set_spawn_pop(j, min_part_type(run), SignCurr(min_part_type(run)))
#ifdef CMPLX_
call set_spawn_pop(j, max_part_type(run), SignCurr(max_part_type(run)))
#endif
else
call update_tau_int(j, run, tau)
call update_shift_int(j, run, DiagSft(run) * tau)
end if
end do
end if
if (tTruncInitiator) then
call CalcParentFlag(j, DetCurr, WalkExcitLevel, parent_flags)
end if
!Debug output.
IFDEBUGTHEN(FCIMCDebug, 3)
write(stdout, "(A,I10,a)", advance='no') 'TW:', j, '['
do part_type = 1, lenof_sign
write(stdout, "(f10.5)", advance='no') SignCurr(part_type)
end do
write(stdout, '(a,i7)', advance='no') '] ', FlagsCurr
call WriteBitDet(stdout, CurrentDets(:, j), .true.)
call neci_flush(stdout)
ENDIFDEBUG
if (walkExcitLevel_toHF == 0) HFInd = j
IFDEBUGTHEN(FCIMCDebug, 1)
if (j > 1) then
if (DetBitEQ(CurrentDets(:, j - 1), CurrentDets(:, j), nifd)) then
call stop_all(this_routine, "Shouldn't have the same determinant twice")
end if
end if
ENDIFDEBUG
! in the guga simulation it is probably better to initialize
! the csf_irmation out here, so it can be used in the
! new way to calculate the reference energy and then the flag
! does not have to be checked each time we loop over the walkers..
if (tGUGA) call fill_csf_i(CurrentDets(0:nifd, j), current_csf_i)
! Sum in any energy contribution from the determinant, including
! other parameters, such as excitlevel info.
! This is where the projected energy is calculated.
call SumEContrib(DetCurr, WalkExcitLevel, SignCurr, CurrentDets(:, j), &
HDiagCurr, HOffDiagCurr, 1.0_dp, tPairedReplicas, j)
if (t_calc_double_occ) then
#ifdef CMPLX_
call stop_all(this_routine, "not implemented for complex")
#else
inst_double_occ = inst_double_occ + &
get_double_occupancy(CurrentDets(:, j), SignCurr)
#endif
end if
if (t_measure_local_spin) then
if (tGUGA) then
call measure_local_spin(SignCurr, current_csf_i)
else
call stop_all(this_routine, "measure_local_spin works only for GUGA")
end if
end if
if (t_spin_measurements) then
call measure_double_occ_and_spin_diff(CurrentDets(:, j), &
DetCurr, SignCurr)
end if
! If we're on the Hartree-Fock, and all singles and doubles are in
! the core space, then there will be no stochastic spawning from
! this determinant, so we can the rest of this loop.
if (tSemiStochastic .and. ss_space_in%tDoubles .and. walkExcitLevel_toHF == 0 .and. tDetermHFSpawning) cycle
! For HPHFs, get the diagonal Hamiltonian element without the
! cross-term correction. This will make calculating the diagonal
! elements for new spawnins more efficient
if (tPreCond .or. tReplicaEstimates) then
if (tHPHF) then
hdiag_bare = get_hdiag_bare_hphf(DetCurr, CurrentDets(:, j), EnergyCurr)
else
hdiag_bare = EnergyCurr
end if
end if
! Loop over the 'type' of particle.
! lenof_sign == 1 --> Only real particles
! lenof_sign == 2 --> complex walkers
! --> part_type == 1, 2; real and complex walkers
! --> OR double run
! --> part_type == 1, 2; population sets 1 and 2, both real
! alis additional idea to skip the number of attempted excitations
! for noninititators in the back-spawning approach
! remove that for now
! try a mixed excitation scheme for guga, where we only do a full
! excitation for initiators and the crude one for non-inits
if (tGen_guga_mixed) then
if (t_guga_mixed_init) then
flag_mixed = any_run_is_initiator(CurrentDets(:, j))
else if (t_guga_mixed_semi) then
if (tSemiStochastic) then
flag_mixed = tGlobalCoreDet
else
flag_mixed = any_run_is_initiator(CurrentDets(:, j))
end if
end if
if (flag_mixed) then
generate_excitation => generate_excitation_guga
else
if (tGen_sym_guga_mol) then
generate_excitation => gen_excit_4ind_weighted2
else if (tGen_sym_guga_ueg) then
if (t_new_real_space_hubbard) then
generate_excitation => gen_excit_rs_hubbard
else if (t_k_space_hubbard) then
generate_excitation => gen_excit_k_space_hub
end if
end if
end if
end if
if (t_crude_exchange_noninits .or. t_approx_exchange_noninits) then
is_init_guga = any_run_is_initiator(CurrentDets(:, j))
end if
if (t_trunc_guga_pgen_noninits) then
is_init_guga = any_run_is_initiator(CurrentDets(:, j))
end if
if (t_guga_back_spawn) then
is_init_guga = any_run_is_initiator(CurrentDets(:, j))
end if
loop_over_type : do part_type = 1, lenof_sign
run = part_type_to_run(part_type)
TempSpawnedPartsInd = 0
! Loop over all the particles of a given type on the
! determinant. CurrentSign gives number of walkers. Multiply
! up by AvMCExcits if attempting multiple excitations from
! each walker (default 1.0_dp).
AvMCExcitsLoc = AvMCExcits
! optional: Adjust the number of spawns to the expected maximum
! Hij/pgen ratio of this determinant -> prevent blooms
! Only done while not updating tau (to prevent interdependencies)
if (tScaleBlooms .and. tau_search_method == possible_tau_search_methods%off) then
max_spawn = tau * get_max_ratio(j)
if (max_spawn > max_allowed_spawn) then
scale = max_spawn / max_allowed_spawn
AvMCExcitsLoc = AvMCExcitsLoc * scale
end if
end if
call decide_num_to_spawn(SignCurr(part_type), AvMCExcitsLoc, WalkersToSpawn)
loop_over_walkers : do p = 1, WalkersToSpawn
! Zero the bit representation, to ensure no extraneous
! data gets through.
ilutnJ = 0_n_int
child = 0.0_dp
! for the 3-body excitations i really do not want to change
! all the interfaces to the other excitation generators,
! which all just assume ex(2,2) as size.. so use a
! if here..
! Generate a (random) excitation
call generate_excitation(DetCurr, CurrentDets(:,j), nJ, &
ilutnJ, exFlag, IC, ex, tParity, prob, &
HElGen, fcimc_excit_gen_store, part_type)
!If we are fixing the population of reference det, skip spawing into it.
if (tSkipRef(run) .and. all(nJ == projEdet(:, run))) then
!Set nJ to null
nJ(1) = 0
end if
! If a valid excitation, see if we should spawn children.
if (IsNullDet(nJ)) then
nInvalidExcits = nInvalidExcits + 1
else
nValidExcits = nValidExcits + 1
if (tSemiStochastic) then
call encode_child(CurrentDets(:, j), iLutnJ, ic, ex)
! Temporary fix: FindExcitBitDet copies the flags of the parent onto the
! child, which causes semi-stochastic simulations to crash. Should it copy
! these flags? There are comments questioning this in create_particle, too.
iLutnJ(IlutBits%ind_flag) = 0_n_int
! If the parent state in the core space.
if (check_determ_flag(CurrentDets(:, j), run)) then
! Is the spawned state in the core space?
tInDetermSpace = is_core_state(iLutnJ, nJ, run)
! If spawning is from and to the core space, cancel it.
if (tInDetermSpace) cycle
! Set the flag to specify that the spawning is occuring
! from the core space.
call set_flag(iLutnJ, flag_determ_parent)
end if
end if
if (tPreCond .or. tReplicaEstimates) then
hdiag_spawn = get_hdiag_from_excit(DetCurr, nJ, &
iLutnJ, ic, ex, hdiag_bare)
if (tPreCond) then
precond_fac = hdiag_spawn - proj_e_for_precond(part_type) - &
proje_ref_energy_offsets(part_type) - Hii
end if
end if
child = attempt_create(DetCurr, &
CurrentDets(:, j), SignCurr, &
nJ, iLutnJ, Prob, HElGen, IC, ex, &
tParity, walkExcitLevel, part_type, &
AvSignCurr, AvMCExcitsLoc, RDMBiasFacCurr, precond_fac)
! Note these last two, AvSignCurr and
! RDMBiasFacCurr are not used unless we're
! doing an RDM calculation.
end if
IFDEBUG(FCIMCDebug, 3) then
write(stdout, '(a)', advance='no') 'SP: ['
do y = 1, lenof_sign
write(stdout, '(f12.5)', advance='no') &
child(y)
end do
write(stdout, '("] ")', advance='no')
call write_det(stdout, nJ, .true.)
call neci_flush(stdout)
end if
! Children have been chosen to be spawned.
is_child_created : if (.not. all(near_zero(child))) then
! Encode child if not done already.
if (.not. tSemiStochastic) then
call encode_child(CurrentDets(:, j), iLutnJ, ic, ex)
end if
! FindExcitBitDet copies the parent flags so that unwanted flags must be unset.
! Should it really do this?
if (tTrialWavefunction) then
call clr_flag(iLutnJ, flag_trial)
call clr_flag(iLutnJ, flag_connected)
end if
! If using a preconditioner, update the child weight for statistics
! (mainly for blooms and hence updating the time step).
! Note, the final preconiditoner is applied in annihilation, once
! the exact projected energy is know.
child_for_stats = child
if (tPreCond) child_for_stats = child_for_stats / precond_fac
if (tScaleBlooms) call update_max_ratio(abs(HElGen) / prob, j)
call new_child_stats_normal(iter_data, CurrentDets(:, j), &
nJ, iLutnJ, ic, walkExcitLevel, &
child_for_stats, parent_flags, part_type)
if (use_spawn_hash_table) then
call create_particle_with_hash_table(nJ, ilutnJ, child, part_type, &
CurrentDets(:, j), iter_data, err)
if (err /= 0) call stop_all(this_routine, 'failure from create_particle_with_hash_table')
else
call create_particle(nJ, iLutnJ, child, part_type, hdiag_spawn, err, &
CurrentDets(:, j), SignCurr, p, &
RDMBiasFacCurr, WalkersToSpawn, abs(HElGen), j)
if (err /= 0) call stop_all(this_routine, 'failure from create_particle')
end if
end if is_child_created ! (child /= 0), Child created.
end do loop_over_walkers
end do loop_over_type ! Cycling over 'type' of particle on a given determinant.
! If we are performing a semi-stochastic simulation and this state
! is in the deterministic space, then the death step is performed
! deterministically later. Otherwise, perform the death step now.
! If using a preconditioner, then death is done in the annihilation
! routine, after the energy has been calculated.
if (tDeathBeforeComms) then
call walker_death(iter_data, DetCurr, CurrentDets(:, j), &
HDiagCurr, SignCurr, j, WalkExcitLevel, &
t_core_die_=.false.)
end if
end do loop_over_determinants
!loop timing for this iteration on this MPI rank
lt_arr(mod(Iter - 1, 100) + 1) = mpi_wtime() - lstart
! if any proc ran out of memory, terminate
call MPISumAll(err, allErr)
if (allErr /= 0) then
err = allErr
return
end if
IFDEBUGTHEN(FCIMCDebug, 2)
write(stdout, *) 'Finished loop over determinants'
write(stdout, *) "Holes in list: ", iEndFreeSlot
ENDIFDEBUG
if (tSemiStochastic) then
! For semi-stochastic calculations only: Gather together the parts
! of the deterministic vector stored on each processor, and then
! perform the multiplication of the exact projector on this vector.
if (tDeathBeforeComms) then
call determ_projection()
else
call determ_projection_no_death()
end if
if (tFillingStochRDMonFly) then
! For RDM calculations, add the current core amplitudes into the
! running average.
call average_determ_vector()
! If this is an iteration where the RDM energy is printed then
! add the off-diagonal contributions from the core determinants
! (the diagonal contributions are done in the same place for
! all determinants, regardless of whether they are core or not,
! so are not added in here).
if (tFill_RDM) then
call fill_RDM_offdiag_deterministic(rdm_definitions, two_rdm_spawn, one_rdms)
! deterministic space is always only initiators, so it fully counts towards
! the initiator-only RDMs
if (tInitsRDM) call fill_RDM_offdiag_deterministic(rdm_inits_defs, &
two_rdm_inits_spawn, inits_one_rdms)
end if
end if
end if
! With this algorithm, the determinants do not move, and therefore
! TotWalkersNew is simply equal to TotWalkers
TotWalkersNew = int(TotWalkers)
! Update the statistics for the end of an iteration.
! Why is this done here - before annihilation!
call end_iter_stats(TotWalkersNew)
! Print bloom/memory warnings
call end_iteration_print_warn(totWalkersNew)
call halt_timer(walker_time)
! For the direct annihilation algorithm. The newly spawned
! walkers should be in a seperate array (SpawnedParts) and the other
! list should be ordered.
call set_timer(annihil_time, 30)
!HolesInList is returned from direct annihilation with the number of unoccupied determinants in the list
!They have already been removed from the hash table though.
call communicate_and_merge_spawns(MaxIndex, iter_data, .false.)
if (tSimpleInit) then
call get_ests_from_spawns_simple(MaxIndex, proj_e_for_precond)
else
call get_ests_from_spawns(MaxIndex, proj_e_for_precond)
end if
if (tSimpleInit) call rm_non_inits_from_spawnedparts(MaxIndex, iter_data)
! If performing FCIQMC with preconditioning, then apply the
! the preconditioner to the spawnings, and perform the death step.
if (tPreCond) then
call set_timer(rescale_time, 30)
call rescale_spawns(MaxIndex, proj_e_for_precond, iter_data)
call halt_timer(rescale_time)
end if
! If we haven't performed the death step yet, then do it now.
if (.not. tDeathBeforeComms) then
call set_timer(death_time, 30)
call perform_death_all_walkers(iter_data)
call halt_timer(death_time)
end if
call DirectAnnihilation(TotWalkersNew, MaxIndex, iter_data, err)
! if any proc ran out of memory, terminate
call MPISumAll(err, allErr)
if (allErr /= 0) then
err = allErr
return
end if
! The growth in the size of the occupied part of CurrentDets
! this is important for the purpose of prone_walkers
detGrowth = int(TotWalkersNew - TotWalkers)
! This indicates the number of determinants in the list + the number
! of holes that have been introduced due to annihilation.
TotWalkers = TotWalkersNew
! if we still have plenty of empty slots in the list, deactivate the decay
if (t_activate_decay .and. TotWalkers < 0.95_dp * real(MaxWalkersPart, dp)) then
t_activate_decay = .false.
end if
! The superinitiators are now the same as they will be at the beginning of
! the next cycle (this flag is reset if they change)
tReferenceChanged = .false.
CALL halt_timer(Annihil_Time)
IFDEBUG(FCIMCDebug, 2) write(stdout, *) "Finished Annihilation step"
! If we are orthogonalising the replica wavefunctions, to generate
! excited states, then do that here.
if (tOrthogonaliseReplicas .and. iter > orthogonalise_iter) then
call orthogonalise_replicas(iter_data)
else if (tPrintReplicaOverlaps .and. inum_runs > 1) then
call calc_replica_overlaps()
end if
if (tFillingStochRDMonFly) then
! if we use the initiator-only rdms as gamma_0, get them in their own entity
if (tInitsRDM) call fill_rdm_diag_wrapper(rdm_inits_defs, two_rdm_inits_spawn, &
inits_one_rdms, CurrentDets, int(TotWalkers), .false., .false.)
call fill_rdm_diag_wrapper(rdm_definitions, two_rdm_spawn, one_rdms, &
CurrentDets, int(TotWalkers), .true., tApplyLC)
end if
if (tTrialWavefunction .and. tTrialShift) then
call fix_trial_overlap(iter_data)
end if
call update_iter_data(iter_data)
! This routine will take the CurrentDets and search the array to find all single and double
! connections - adding them into the RDM's.
! This explicit way of doing this is very expensive, but o.k for very small systems.
if (tFillingExplicRDMonFly) then
if (tHistSpawn) THEN
call Fill_Hist_ExplicitRDM_this_Iter()
else
call Fill_ExplicitRDM_this_Iter(TotWalkers)
end if
end if
if (tFillingStochRDMonFly .or. tFillingExplicRDMonFly) then
! Fill the receiving RDM list from the beginning.
two_rdm_recv%nelements = 0
call communicate_rdm_spawn_t(two_rdm_spawn, two_rdm_recv)
call add_rdm_1_to_rdm_2(two_rdm_recv, two_rdm_main)
two_rdm_recv%nelements = 0
if (tInitsRDM) then
call communicate_rdm_spawn_t(two_rdm_inits_spawn, two_rdm_recv)
call add_rdm_1_to_rdm_2(two_rdm_recv, two_rdm_inits)
end if
end if
end subroutine PerformFCIMCycPar