subroutine update_shift(iter_data, replica_pairs)
use CalcData, only: tInstGrowthRate, tL2GrowRate
type(fcimc_iter_data), intent(in) :: iter_data
logical, intent(in) :: replica_pairs
integer(int64) :: tot_walkers
logical, dimension(inum_runs) :: tReZeroShift
real(dp), dimension(inum_runs) :: AllGrowRateRe, AllGrowRateIm
real(dp), dimension(inum_runs) :: AllHFGrowRate, AllWalkers
real(dp), dimension(lenof_sign) :: denominator, all_denominator
real(dp), dimension(inum_runs) :: rel_tot_trial_numerator
integer :: error, i, proc, pos, run, lb, ub
logical, dimension(inum_runs) :: defer_update
logical :: start_varying_shift
character(*), parameter :: this_routine = 'update_shift'
! Normally we allow the shift to vary depending on the conditions
! tested. Sometimes we want to defer this to the next cycle...
defer_update(:) = .false.
! collate_iter_data --> The values used are only valid on Root
i_am_root : if (iProcIndex == Root) then
if (tL2GrowRate) then
! use the L2 norm to determine the growrate
AllGrowRate(:) = norm_psi(:) / old_norm_psi(:)
AllWalkers(:) = norm_psi(:)
else if (tInstGrowthRate) then
! Calculate the growth rate simply using the two points at
! the beginning and the end of the update cycle.
do run = 1, inum_runs
lb = min_part_type(run)
ub = max_part_type(run)
AllGrowRate(run) = (sum(iter_data%update_growth_tot(lb:ub) &
+ iter_data%tot_parts_old(lb:ub))) &
/ real(sum(iter_data%tot_parts_old(lb:ub)), dp)
AllWalkers(run) = (sum(iter_data%update_growth_tot(lb:ub) &
+ iter_data%tot_parts_old(lb:ub)))
end do
else
! Instead attempt to calculate the average growth over every
! iteration over the update cycle
if (all(.not. near_zero(OldAllAvWalkersCyc))) then
AllGrowRate(:) = AllSumWalkersCyc(:) / real(StepsSft, dp) &
/ OldAllAvWalkersCyc(:)
end if
AllWalkers(:) = AllSumWalkersCyc(:) / real(StepsSft, dp)
end if
! For complex case, obtain both Re and Im parts
#ifdef CMPLX_
do run = 1, inum_runs
lb = min_part_type(run)
ub = max_part_type(run)
if (.not. near_zero(iter_data%tot_parts_old(lb))) then
AllGrowRateRe(run) = &
(iter_data%update_growth_tot(lb) + iter_data%tot_parts_old(lb)) &
/ iter_data%tot_parts_old(lb)
end if
if (.not. near_zero(iter_data%tot_parts_old(ub))) then
AllGrowRateIm(run) = &
(iter_data%update_growth_tot(ub) + iter_data%tot_parts_old(ub)) &
/ iter_data%tot_parts_old(ub)
end if
end do
#endif
! If any run uses the fixtrial option, we need to add the offset to the
! trial numerator
if (tTrialWavefunction .and. tTrialShift) &
rel_tot_trial_numerator = real(relative_trial_numerator( &
tot_trial_numerator, tot_trial_denom, replica_pairs), dp)
! Exit the single particle phase if the number of walkers exceeds
! the value in the input file. If particle no has fallen, re-enter
! it.
tReZeroShift = .false.
do run = 1, inum_runs
lb = min_part_type(run)
ub = max_part_type(run)
if (tTrialShift .and. .not. tFixTrial(run) .and. tTrialWavefunction .and. abs(tot_trial_denom(run)) >= TrialTarget) then
!When reaching target overlap with trial wavefunction, set flag to keep it fixed.
tFixTrial(run) = .True.
write(stdout, '(a,i13,a,i1)') 'Exiting the varaible shift phase on iteration: ' &
, iter + PreviousCycles, ' - overlap with trial wavefunction of the following run is now fixed: ', run
end if
if (tFixedN0) then
if (tModShiftTarget) tModShiftPhase = .true.
if (.not. tSkipRef(run) .and. abs(AllHFCyc(run)) >= N0_Target) then
!When reaching target N0, set flag to keep the population of reference det fixed.
tSkipRef(run) = .True.
write(stdout, '(a,i13,a,i1)') 'Exiting the fixed shift phase on iteration: ' &
, iter + PreviousCycles, ' - reference population of the following run is now fixed: ', run
!Set these parameters because other parts of the code depends on them
VaryShiftIter(run) = Iter
iBlockingIter(run) = Iter + PreviousCycles
tSinglePartPhase(run) = .false.
end if
if (tSkipRef(run)) then
!Use the projected energy as the shift to fix the
!population of the reference det and thus reduce the
!fluctuations of the projected energy.
!ToDo: Make DiagSft complex
DiagSft(run) = real((AllENumCyc(run)) &
/ (AllHFCyc(run)) + proje_ref_energy_offsets(run), dp)
if (tModShiftTarget) then
tot_walkers = int(InitWalkers, int64) * int(nNodes, int64)
TargetWalkers2 = tot_walkers - AllNoModShiftWalk(run)
if (TargetWalkers2 < 1.0_dp) TargetWalkers2 = 1.0_dp
if (iter < (ModShiftTargetIter+ModShiftConstIters)) then
DiagSft2(run) = DiagSft2Init
else
DiagSft2(run) = &
DiagSft2(run) &
- (log(AllNoNonModShiftWalk(run)/OldAllNoNonModShiftWalk(run)) * SftDamp &
+ log(AllNoNonModShiftWalk(run)/TargetWalkers2) * SftDamp2) &
/ (Tau * StepsSft)
end if
end if
! Update the shift averages
if ((iter - VaryShiftIter(run)) >= nShiftEquilSteps) then
if ((iter - VaryShiftIter(run) - nShiftEquilSteps) < StepsSft) &
write(stdout, '(a,i14)') 'Beginning to average shift value on iteration: ', iter + PreviousCycles
VaryShiftCycles(run) = VaryShiftCycles(run) + 1
SumDiagSft(run) = SumDiagSft(run) + DiagSft(run)
AvDiagSft(run) = SumDiagSft(run) / real(VaryShiftCycles(run), dp)
end if
else
!Keep shift equal to input till target reference population is reached.
DiagSft(run) = InputDiagSft(run)
if (tModShiftTarget) DiagSft2(run) = DiagSft2Init
end if
else if (tFixTrial(run)) then
!Use the trial energy as the shift to fix the
!overlap with trial wavefunction and thus reduce the
!fluctuations of the trial energy.
!ToDo: Make DiafSft complex
DiagSft(run) = real((rel_tot_trial_numerator(run) / tot_trial_denom(run)) - Hii, dp)
! Update the shift averages
if ((iter - VaryShiftIter(run)) >= nShiftEquilSteps) then
if ((iter - VaryShiftIter(run) - nShiftEquilSteps) < StepsSft) &
write(stdout, '(a,i14)') 'Beginning to average shift value on iteration: ', iter + PreviousCycles
VaryShiftCycles(run) = VaryShiftCycles(run) + 1
SumDiagSft(run) = SumDiagSft(run) + DiagSft(run)
AvDiagSft(run) = SumDiagSft(run) / real(VaryShiftCycles(run), dp)
end if
else if(.not. tRenorm) then!not Fixed-N0 and not Trial-Shift
tot_walkers = int(InitWalkers, int64) * int(nNodes, int64)
single_part_phase : if (TSinglePartPhase(run)) then
#ifdef CMPLX_
if ((sum(AllTotParts(lb:ub)) > tot_walkers) .or. &
(abs_sign(AllNoatHF(lb:ub)) > MaxNoatHF)) then
write(stdout, '(a,i13,a)') 'Exiting the single particle growth phase on iteration: ', iter + PreviousCycles, &
' - Shift can now change'
VaryShiftIter(run) = Iter
iBlockingIter(run) = Iter + PreviousCycles
tSinglePartPhase(run) = .false.
if (abs(TargetGrowRate(run)) > EPS) then
write(stdout, "(A)") "Setting target growth rate to 1."
TargetGrowRate = 0.0_dp
end if
! If enabled, jump the shift to the value preducted by the
! projected energy!
if (tJumpShift) then
if (tJumpShift .and. &
(.not. (isnan(real(proje_iter(run), dp))) .or. &
.not. (is_inf(real(proje_iter(run), dp))))) then
DiagSft(run) = real(proje_iter(run), dp)
defer_update(run) = .true.
end if
end if
end if
#else
start_varying_shift = .false.
if (tLetInitialPopDie) then
if (AllTotParts(run) < tot_walkers) start_varying_shift = .true.
else if (tTargetShiftdamp .or. tBosonNoSpawn) then
start_varying_shift = .true.
else
if ((AllTotParts(run) > tot_walkers) .or. &
(abs(AllNoatHF(run)) > MaxNoatHF)) start_varying_shift = .true.
end if
if (start_varying_shift) then
write(stdout, '(a,i13,a,i1)') 'Exiting the single particle growth phase on iteration: ' &
, iter + PreviousCycles, ' - Shift can now change for population', run
VaryShiftIter(run) = Iter
iBlockingIter(run) = Iter + PreviousCycles
tSinglePartPhase(run) = .false.
! [W.D. 15.5.2017]
! change equal 0 comps
if (abs(TargetGrowRate(run)) > EPS) then
write(stdout, "(A)") "Setting target growth rate to 1."
TargetGrowRate(run) = 0.0_dp
end if
! If enabled, jump the shift to the value preducted by the
! projected energy!
if (tJumpShift) then
DiagSft(run) = real(proje_iter(run), dp)
defer_update(run) = .true.
end if
if (tNonInitModShift .and. .not. tModShiftTarget) then
tModShiftPhase = .true.
end if
end if
#endif
else ! .not.tSinglePartPhase(run)
#ifdef CMPLX_
if (abs_sign(AllNoatHF(lb:ub)) < MaxNoatHF - HFPopThresh) then
#else
if (abs(AllNoatHF(run)) < MaxNoatHF - HFPopThresh) then
#endif
write(stdout, '(a,i13,a)') &
'No at HF has fallen too low - reentering the single particle growth phase on iteration', &
iter + PreviousCycles, ' - particle number may grow again.'
tSinglePartPhase(run) = .true.
tReZeroShift(run) = .true.
end if
end if single_part_phase
! How should the shift change for the entire ensemble of walkers
! over all processors.
if (tModShiftTarget .and. iter > ModShiftTargetIter) then
tModShiftPhase = .true.
end if
if (.not. tRenorm .and. .not. (tSinglePartPhase(run) &
.and. near_zero(TargetGrowRate(run)) &
.or. defer_update(run))) then
!In case we want to continue growing, TargetGrowRate > 0.0_dp
! New shift value
! if(TargetGrowRate(run).ne.0.0_dp) then
! [W.D. 15.5.2017]
if (abs(TargetGrowRate(run)) > EPS) then
#ifdef CMPLX_
if (sum(AllTotParts(lb:ub)) > TargetGrowRateWalk(run)) then
#else
if (AllTotParts(run) > TargetGrowRateWalk(run)) then
#endif
if (tTargetShiftdamp) then
call stop_all(this_routine, &
"Target-shiftdamp not compatible with targetgrowrate!")
end if
!Only allow targetgrowrate to kick in once we have > TargetGrowRateWalk walkers.
DiagSft(run) = DiagSft(run) - (log(AllGrowRate(run) - TargetGrowRate(run)) * SftDamp) / &
(Tau * StepsSft)
! Same for the info shifts for complex walkers
#ifdef CMPLX_
DiagSftRe(run) = DiagSftRe(run) - (log(AllGrowRateRe(run) - TargetGrowRate(run)) * SftDamp) / &
(Tau * StepsSft)
DiagSftIm(run) = DiagSftIm(run) - (log(AllGrowRateIm(run) - TargetGrowRate(run)) * SftDamp) / &
(Tau * StepsSft)
#endif
end if
else
if (tShiftonHFPop) then
!Calculate the shift required to keep the HF population constant
AllHFGrowRate(run) = abs(AllHFCyc(run) / real(StepsSft, dp)) / abs(OldAllHFCyc(run))
DiagSft(run) = DiagSft(run) - (log(AllHFGrowRate(run)) * SftDamp) / &
(Tau * StepsSft)
else if (tInitShift) then
DiagSft(run) = DiagSft(run) - (log(AllNoInitWalk(run)/OldAllNoInitWalk(run)) * SftDamp) / &
(Tau * StepsSft)
else if (tModShiftTarget) then
DiagSft(run) = DiagSft(run) - (log(AllNoModShiftWalk(run)/OldAllNoModShiftWalk(run)) * SftDamp + &
log(AllNoModShiftWalk(run)/TargetWalkers1) * SftDamp2) / &
(Tau * StepsSft)
if (iter < (modShiftTargetIter + modShiftConstIters)) then
DiagSft2(run) = DiagSft2Init
else
DiagSft2(run) = DiagSft2(run) - (log(AllNoNonModShiftWalk(run)/OldAllNoNonModShiftWalk(run)) * SftDamp + &
log(AllNoNonModShiftWalk(run)/TargetWalkers2) * SftDamp2) / &
(Tau * StepsSft)
end if
else if (tBosonNoSpawn) then
DiagSft(run) = DiagSft(run) - (log(AllNoBosonActiveWalk(run)/OldAllNoBosonActiveWalk(run)) * SftDamp + &
log(AllNoBosonActiveWalk(run)/tot_walkers) * SftDamp2) / &
(Tau * StepsSft)
else if (tTargetShiftdamp) then
if (.not. near_zero(AllGrowRate(run)) .and. .not. near_zero(AllWalkers(run))) then
DiagSft(run) = DiagSft(run) - (log(AllGrowRate(run)) * SftDamp + &
log(AllWalkers(run)/tot_walkers) * SftDamp2) / &
(Tau * StepsSft)
else
call stop_all(this_routine, "Shift undefined because walker growth rate is zero. Aborting.")
end if
else
if (.not. near_zero(AllGrowRate(run))) then
DiagSft(run) = DiagSft(run) - (log(AllGrowRate(run)) * SftDamp) / &
(Tau * StepsSft)
else
call stop_all(this_routine, "Shift undefined because walker growth rate is zero. Aborting.")
end if
end if
end if
! Update the shift averages
if ((iter - VaryShiftIter(run)) >= nShiftEquilSteps) then
if ((iter - VaryShiftIter(run) - nShiftEquilSteps) < StepsSft) &
write(stdout, '(a,i14)') 'Beginning to average shift value on iteration: ', iter + PreviousCycles
VaryShiftCycles(run) = VaryShiftCycles(run) + 1
SumDiagSft(run) = SumDiagSft(run) + DiagSft(run)
AvDiagSft(run) = SumDiagSft(run) / real(VaryShiftCycles(run), dp)
end if
end if
end if !tFixedN0 or not
! only update the shift this way if possible
if (abs_sign(AllNoatHF(lb:ub)) > EPS) then
#ifdef CMPLX_
! Calculate the instantaneous 'shift' from the HF population
HFShift(run) = -1.0_dp / abs_sign(AllNoatHF(lb:ub)) * &
(abs_sign(AllNoatHF(lb:ub)) - abs_sign(OldAllNoatHF(lb:ub)) / &
(Tau * real(StepsSft, dp)))
InstShift(run) = -1.0_dp / sum(AllTotParts(lb:ub)) * &
((sum(AllTotParts(lb:ub)) - sum(AllTotPartsOld(lb:ub))) / &
(Tau * real(StepsSft, dp)))
#else
! Calculate the instantaneous 'shift' from the HF population
HFShift(run) = -1.0_dp / abs(AllNoatHF(run)) * &
(abs(AllNoatHF(run)) - abs(OldAllNoatHF(run)) / &
(Tau * real(StepsSft, dp)))
InstShift(run) = -1.0_dp / AllTotParts(run) * &
((AllTotParts(run) - AllTotPartsOld(run)) / &
(Tau * real(StepsSft, dp)))
#endif
end if
! When using a linear combination, the denominator is summed
! directly.
all_sum_proje_denominator(run) = ARR_RE_OR_CPLX(AllSumNoatHF, run)
all_cyc_proje_denominator(run) = AllHFCyc(run)
! Calculate the projected energy.
if (.not. near_zero(AllSumNoatHF(run))) then
ProjectionE(run) = (AllSumENum(run)) / (all_sum_proje_denominator(run)) &
+ proje_ref_energy_offsets(run)
end if
if (abs(AllHFCyc(run)) > EPS) then
proje_iter(run) = (AllENumCyc(run)) / (all_cyc_proje_denominator(run)) &
+ proje_ref_energy_offsets(run)
AbsProjE(run) = (AllENumCycAbs(run)) / (all_cyc_proje_denominator(run)) &
+ proje_ref_energy_offsets(run)
inits_proje_iter(run) = (AllInitsENumCyc(run)) / (all_cyc_proje_denominator(run)) &
+ proje_ref_energy_offsets(run)
end if
! If we are re-zeroing the shift
if (tReZeroShift(run)) then
DiagSft(run) = 0.0_dp
DiagSft2(run) = 0.0_dp
VaryShiftCycles(run) = 0
SumDiagSft(run) = 0.0_dp
AvDiagSft(run) = 0.0_dp
end if
end do
! Get some totalled values
if (abs(sum(all_sum_proje_denominator(1:inum_runs))) > EPS) then
projectionE_tot = sum(AllSumENum(1:inum_runs)) &
/ sum(all_sum_proje_denominator(1:inum_runs))
end if
if (abs(sum(all_cyc_proje_denominator(1:inum_runs))) > EPS) then
proje_iter_tot = sum(AllENumCyc(1:inum_runs)) &
/ sum(all_cyc_proje_denominator(1:inum_runs))
inits_proje_iter_tot = sum(AllInitsENumCyc(1:inum_runs)) &
/ sum(all_cyc_proje_denominator(1:inum_runs))
end if
end if i_am_root
! Broadcast the shift from root to all the other processors
call MPIBcast(tSinglePartPhase)
call MPIBCast(tModShiftPhase)
call MPIBcast(VaryShiftIter)
call MPIBcast(DiagSft)
call MPIBcast(DiagSft2)
call MPIBcast(tSkipRef)
call MPIBcast(tFixTrial)
call MPIBcast(VaryShiftCycles)
call MPIBcast(SumDiagSft)
call MPIBcast(AvDiagSft)
do run = 1, inum_runs
if (.not. tSinglePartPhase(run)) then
TargetGrowRate(run) = 0.0_dp
end if
end do
if (tau_search_method /= possible_tau_search_methods%off) then
if (end_of_search_reached(tau_search_method, tau_stop_method)) then
call stop_tau_search(tau_stop_method)
end if
end if
end subroutine update_shift