subroutine InitFCIMC_pops(iPopAllTotWalkers, PopNIfSgn, PopNel, &
pops_nnodes, pops_walkers, perturbs, &
PopBalanceBlocks, PopDiagSft, source_name)
use CalcData, only: iReadWalkersRoot
use hash, only: clear_hash_table, fill_in_hash_table
use perturbations, only: apply_perturbation_array
use semi_stoch_procs, only: fill_in_diag_helements
integer(int64), intent(in) :: iPopAllTotWalkers
integer, intent(in) :: PopNIfSgn
integer, intent(inout) :: PopBalanceBlocks
! Note that PopNel might be recalculated in ReadFromPopsfile (and might
! not have even been set on input).
integer, intent(inout) :: PopNel
integer, intent(in) :: pops_nnodes
integer(int64), intent(in) :: pops_walkers(0:nProcessors - 1)
real(dp), intent(in) :: PopDiagSft(inum_runs)
! Perturbation operators to apply to the determinants after they have
! been read in.
type(perturbation), intent(in), allocatable, optional :: perturbs(:)
character(*), intent(in), optional :: source_name
integer(int64) :: tot_walkers
integer :: run, ReadBatch
logical :: apply_pert
integer :: TotWalkersIn
character(255) :: identifier
integer(int64) :: l
real(dp) :: TempSign(lenof_sign)
character(len=*), parameter :: this_routine = "InitFCIMC_pops"
type(gdata_io_t) :: gdata_read_handler
HElement_t(dp):: InstE(inum_runs)
if (iReadWalkersRoot == 0) then
! ReadBatch is the number of walkers to read in from the
! popsfile at one time. The larger it is, the fewer
! communictions will be needed to scatter the particles.
!
! By default, the new array (which is only created on the
! root processors) is the same length as the spawning
! arrays.
ReadBatch = MaxSpawned
else
ReadBatch = iReadWalkersRoot
end if
apply_pert = .false.
if (present(perturbs)) then
if (allocated(perturbs)) apply_pert = .true.
end if
! Only active the accumlation of populations if the popsfile has
! apvals and accum-pops is specified with a proper iteration
tAccumPopsActive = .false.
if (tAccumPops .and. tPopAccumPops) then
if (iAccumPopsIter <= PreviousCycles) then
tAccumPopsActive = .true.
iAccumPopsCounter = PopAccumPopsCounter
write (stdout, *) "Accumulated populations are found. Accumulation will continue."
else
write (stdout, *) "Accumulated populations are found, but will be discarded."
write (stdout, *) "Accumulation will restart at iteration: ", iAccumPopsIter
endif
endif
! decide which global det data is read
call gdata_read_handler%init_gdata_io(tPopAutoAdaptiveShift, tPopScaleBlooms, &
tPopAccumPops, fvals_size, max_ratio_size, apvals_size)
if (present(source_name)) then
identifier = source_name
else
if (tPopsAlias) then
identifier = aliasStem
else
identifier = 'POPSFILE'
endif
endif
! If applying perturbations, read the popsfile into the array
! popsfile_dets and then apply the perturbations to the determinants
! in this array.
! If not, then read the popsfile straight into CurrentDets.
if (apply_pert) then
call ReadFromPopsfile(iPopAllTotWalkers, ReadBatch, TotWalkers, TotParts, NoatHF, &
popsfile_dets, MaxWalkersPart, pops_nnodes, pops_walkers, PopNIfSgn, &
PopNel, PopBalanceBlocks, gdata_read_handler, tCalcExtraInfo=.false., &
filename_stem=identifier)
write (stdout, *) "Applying perturbation to read-in walker confguration!"
write (stdout, *) "Total number of walkers before perturbation: ", TotWalkers
! also store the original walker number in the real-time FCIQMC
TotWalkersIn = int(TotWalkers)
! also store this original value
if (t_real_time_fciqmc) then
TotWalkers_orig = TotWalkersIn
if (t_kspace_operators .and. tReal) then
call apply_perturbation_array(perturbs, TotWalkersIn, popsfile_dets, &
CurrentDets, phase_factors)
else
call apply_perturbation_array(perturbs, TotWalkersIn, popsfile_dets, &
CurrentDets)
endif
else
call apply_perturbation_array(perturbs, TotWalkersIn, popsfile_dets, CurrentDets)
endif
TotWalkers = int(TotWalkersIn, int64)
write (stdout, *) "Total number of walkers after perturbation: ", TotWalkers
else
call ReadFromPopsfile(iPopAllTotWalkers, ReadBatch, TotWalkers, TotParts, NoatHF, &
CurrentDets, MaxWalkersPart, pops_nnodes, pops_walkers, PopNIfSgn, &
PopNel, PopBalanceBlocks, gdata_read_handler, tCalcExtraInfo=.false., &
filename_stem=identifier)
if (t_real_time_fciqmc) TotWalkers_orig = TotWalkers
end if
if (abs(ScaleWalkers - 1) > 1.0e-12_dp) then
write (stdout, *) "Rescaling walkers by a factor of: ", ScaleWalkers
do l = 1, TotWalkers
call extract_sign(CurrentDets(:, l), TempSign)
!do j = 1, lenof_sign
!IntegerPart=int(ScaleWalkers*TempSign(j))
!FracPart=TempSign(j)-real(IntegerPart,dp)
!r = genrand_real2_dSFMT()
!if(r.lt.abs(FracPart)) then
!if(FracPart.lt.0) then
!TempSign(j)=TempSign(j)-1
!else
!TempSign(j)=TempSign(j)+1
!endif
!endif
!end do
call encode_sign(CurrentDets(:, l), TempSign * ScaleWalkers)
enddo
SumNoatHF = SumNoatHF * ScaleWalkers
AllSumNoatHF = AllSumNoatHF * ScaleWalkers
SumENum = SumENum * ScaleWalkers
AllSumENum = AllSumENum * ScaleWalkers
InstNoatHF = InstNoatHF * ScaleWalkers
end if
call fill_in_diag_helements()
if (associated(lookup_supergroup_indexer)) then
block
integer :: i, nI(nel)
do i = 1, int(TotWalkers)
call decode_bit_det(nI, CurrentDets(:, i))
call set_supergroup_idx( &
i, lookup_supergroup_indexer%idx_nI(nI))
end do
end block
end if
call clear_hash_table(HashIndex)
call fill_in_hash_table(HashIndex, nWalkerHashes, CurrentDets, int(TotWalkers), .true.)
! If we are doing load balancing, do an initial balance now that we
! have read the particles in
! n.b. This must be done after the hash tables have been initialised
! properly, or everything will break!
if (tLoadBalanceBlocks) then
! n.b. The free-slot list will be empty after reading from POPS
! (which are densely packed). Not otherwised initialised by here.
iStartFreeSlot = 1
iEndFreeSlot = 0
call pops_init_balance_blocks(PopBalanceBlocks)
end if
!Set initial global data after potentially load balancing
call set_initial_global_data(TotWalkers, CurrentDets)
! If tWalkContGrow has been set, and therefore the shift reset, in
! when the POPSFILE header was read, check that we aren't already
! beyond the walker threshold. Otherwise unexpected (but technically
! "correct") behaviour occurs.
!
! n.b. 0.95 makes this a soft limit. Otherwise on restarting a calc
! it might have (randomly) dropped just below the target.
if (tWalkContGrow) then
tot_walkers = int(0.95 * InitWalkers * nNodes, int64)
do run = 1, inum_runs
#ifdef CMPLX_
if ((tLetInitialPopDie .and. sum(AllTotParts(min_part_type(run):max_part_type(run))) < tot_walkers) .or. &
((.not. tLetInitialPopDie) .and. sum(AllTotParts(min_part_type(run):max_part_type(run))) > tot_walkers)) then
write (stdout, '("WALKCONTGROW set in input, but simulation already exceeds target number of particles")')
write (stdout, '("Continuing with DIAGSHIFT from POPSFILE")')
if (tHDF5PopsRead) then
root_print "diagshift:", hdf5_diagsft(run)
DiagSft(run) = hdf5_diagsft(run)
else
root_print "diagshift:", PopDiagSft
DiagSft(run) = PopDiagSft(run)
end if
tSinglePartPhase(run) = .false.
end if
#else
if ((tLetInitialPopDie .and. AllTotParts(run) < tot_walkers) .or. &
((.not. tLetInitialPopDie) .and. AllTotParts(run) > tot_walkers)) then
write (stdout, '("WALKCONTGROW set in input, but simulation already exceeds target number of particles")')
write (stdout, '("Continuing with DIAGSHIFT from POPSFILE for run ",i4)') run
if (tHDF5PopsRead) then
root_print "diagshift:", hdf5_diagsft(run)
DiagSft(run) = hdf5_diagsft(run)
else
root_print "diagshift:", PopDiagSft
DiagSft(run) = PopDiagSft(run)
end if
tSinglePartPhase(run) = .false.
end if
#endif
end do
end if
! If necessary, recalculate the instantaneous projected energy, and
! then update the shift to that value.
if (tPopsJumpShift .and. .not. tWalkContGrow) then
call calc_proje(InstE)
DiagSft = real(InstE, dp)
write (stdout, *) 'Calculated instantaneous projected energy', DiagSft
end if
end subroutine InitFCIMC_pops