subroutine update_tau_hist()
! split up the new tau-update routine from the old one!
character(*), parameter :: this_routine = "update_tau_hist"
real(dp) :: psingles_new, tau_new, mpi_tmp, tau_death, pParallel_new, pTriples_new
real(dp) :: ratio_singles, ratio_anti, ratio_para, ratio_doubles, ratio_triples, ratio
logical :: mpi_ltmp
ASSERT(tau_search_method == possible_tau_search_methods%HISTOGRAMMING)
! What needs doing depends on the number of parametrs that are being
! updated.
call MPIAllLORLogical(enough_sing_hist, mpi_ltmp)
enough_sing_hist = mpi_ltmp
call MPIAllLORLogical(enough_doub_hist, mpi_ltmp)
enough_doub_hist = mpi_ltmp
if (t_mol_3_body) then
call MPIAllLORLogical(enough_trip_hist, mpi_ltmp)
enough_trip_hist = mpi_ltmp
end if
! singles is always used..
! thats not quite right.. for the hubbard/UEG case it is not..
if (tUEG .or. tHub .or. &
(t_new_real_space_hubbard .and. .not. t_trans_corr_hop) .or. &
(t_k_space_hubbard .and. .not. t_trans_corr_2body)) then
! also use my new flags and exclude the 2-body transcorrelation
! in the k-space hubbard due to triple excitations and
! parallel doubles
! here i could implement the summing in the case of Hubbard
! and UEG models.. although I could "just" implement the
! optimal time-step in the case of Hubbard models!
! for UEG not i guess..
call integrate_frequency_histogram_spec(frequency_bins, ratio)
if (ratio < 0.0_dp) then
! TODO: and also print out in this case!
! this means i had an int overflow and should stop the tau-searching
root_print "The single excitation histogram is full!"
root_print "stop the hist_tau_search with last time-step: ", tau
tau_search_method = possible_tau_search_methods%OFF
return
end if
if (enough_doub_hist) then
! in this case the enough_doub_hist flag is (mis)used to
! indicate enough spawning events!
! the doubles flag is also used for single excitations in the
! real-space hubbard! be careful
tau_new = MaxWalkerBloom / ratio
! and use the mixing now:
if (t_mix_ratios) then
tau_new = (1.0_dp - mix_ratio) * tau + mix_ratio * tau_new
end if
else
! what do i do in the case if not enough spawns..
! nothing i guess..
tau_new = tau
end if
else
! NOTE: in the case of the 2-body transcorrelated k-space hubbard
! the singles histogram is used to store the triples!
call integrate_frequency_histogram_spec(frequency_bins_singles, &
ratio_singles)
! for analyis print this also out:
! if i have a integer overflow i should probably deal here with it..
if (ratio_singles < 0.0_dp) then
! TODO: and also print out in this case!
! this means i had an int overflow and should stop the tau-searching
root_print "The single excitation histogram is full!"
root_print "stop the hist_tau_search with last time-step: ", tau
root_print "and pSingles and pDoubles:", pSingles, pDoubles
tau_search_method = possible_tau_search_methods%OFF
return
end if
! for tc, also consider the triples
if (t_mol_3_body) then
call integrate_frequency_histogram_spec(frequency_bins_triples, ratio_triples)
! overflow error -> disable hist_tau_search
if (ratio_triples < 0.0_dp) then
tau_search_method = possible_tau_search_methods%OFF
return
end if
else
ratio_triples = 0.0
end if
if (consider_par_bias) then
! sum up all 3 ratios: single, parallel and anti-parallel
call integrate_frequency_histogram_spec(frequency_bins_para, &
ratio_para)
if (ratio_para < 0.0_dp) then
root_print "The parallel excitation histogram is full!"
root_print "stop the hist_tau_search with last time-step: ", tau
root_print "and pSingles and pParallel: ", pSingles, pParallel
tau_search_method = possible_tau_search_methods%OFF
return
end if
call integrate_frequency_histogram_spec(frequency_bins_anti, &
ratio_anti)
if (ratio_anti < 0.0_dp) then
root_print "The anti-parallel excitation histogram is full!"
root_print "stop the hist_tau_search with last time-step: ", tau
root_print "and pSingles and pParallel: ", pSingles, pParallel
tau_search_method = possible_tau_search_methods%OFF
return
end if
! to compare the influences on the time-step:
! change that the ratios are already stored in an unbiased
! way, so no feedback happens!
! also calculate new time-step through this method and
! check the difference to the old method
if (enough_sing_hist .and. enough_doub_hist) then
! puh.. for the mixing i am not sure how to exactly do
! that.. since all of the ratios depend on each other..
! test that!
pparallel_new = ratio_para / (ratio_anti + ratio_para)
if (t_mix_ratios) then
pparallel_new = (1.0_dp - mix_ratio) * pParallel + &
mix_ratio * pparallel_new
end if
psingles_new = ratio_singles * pparallel_new / &
(ratio_para + ratio_singles * pparallel_new)
if (t_mix_ratios) then
psingles_new = (1.0_dp - mix_ratio) * psingles + &
mix_ratio * psingles_new
end if
tau_new = psingles_new * MaxWalkerBloom / ratio_singles
if (t_mix_ratios) then
tau_new = (1.0_dp - mix_ratio) * tau + mix_ratio * tau_new
end if
if (psingles_new > 1e-5_dp .and. &
psingles_new < (1.0_dp - 1e-5_dp)) then
if (abs(psingles - psingles_new) / psingles > 0.0001_dp) then
root_print "Updating singles/doubles bias. pSingles = ", psingles_new
root_print " pDoubles = ", 1.0_dp - psingles_new
end if
pSingles = psingles_new
pDoubles = 1.0_dp - pSingles
end if
! although checking for enough doubles is probably more efficient
! than always checking the reals below..
if (pParallel_new > 1e-4_dp .and. pParallel_new < (1.0_dp - 1e-4_dp)) then
! enough_doub implies that both enough_opp and enough_par are
! true.. so this if statement makes no sense
! and otherwise pParallel_new is the same as before
if (abs(pparallel_new - pParallel) / pParallel > 0.0001_dp) then
root_print "Updating parallel-spin bias; new pParallel = ", pParallel_new
end if
! in this new implementation the weighting make pParallel
! smaller and smaller.. so also limit it to some lower bound
pParallel = pParallel_new
end if
else
pparallel_new = pParallel
psingles_new = pSingles
! im not sure if this possible division by 0 is cool...
! and if this is somehow cool in the histogramming
! approach.. this again is some kind of worst case
! adaptation.. hm.. todo
if (abs(ratio_singles) > EPS .or. abs(ratio_para) > EPS &
.or. abs(ratio_anti) > EPS) then
tau_new = MaxWalkerBloom * min( &
pSingles / max(EPS, ratio_singles), &
pDoubles * pParallel / max(EPS, ratio_para), &
pDoubles * (1.0_dp - pParallel) / max(EPS, ratio_anti))
else
tau_new = tau
end if
end if
else
! here i only use doubles for now
call integrate_frequency_histogram_spec(frequency_bins_doubles, &
ratio_doubles)
if (ratio_doubles < 0.0_dp) then
root_print "The double excitation histogram is full!"
root_print "stop the hist_tau_search with last time-step: ", tau
root_print "and pSingles and pDoubles: ", pSingles, pDoubles
tau_search_method = possible_tau_search_methods%OFF
return
end if
! to compare the influences on the time-step:
! change that so it does store the ratio unbiased
if (enough_sing_hist .and. enough_doub_hist) then
psingles_new = ratio_singles / (ratio_doubles + ratio_singles)
if (t_mix_ratios) then
psingles_new = (1.0_dp - mix_ratio) * pSingles + &
mix_ratio * psingles_new
end if
tau_new = MaxWalkerBloom / (ratio_doubles + ratio_singles + ratio_triples)
if (t_mix_ratios) then
tau_new = (1.0_dp - mix_ratio) * tau + mix_ratio * tau_new
end if
if (psingles_new > 1e-5_dp .and. &
psingles_new < (1.0_dp - 1e-5_dp)) then
if (abs(psingles - psingles_new) / psingles > 0.0001_dp) then
root_print "Updating singles/doubles bias. pSingles = ", psingles_new
root_print " pDoubles = ", 1.0_dp - psingles_new
end if
pSingles = psingles_new
pDoubles = 1.0_dp - pSingles
end if
! update pTriples analogously - note that as pTriples already modifies
! the effective pSingles/pDoubles, they do not have to be updated
! with ratio_triples
if (t_mol_3_body .and. enough_doub_hist) then
pTriples_new = ratio_triples / (ratio_doubles + ratio_singles + &
ratio_triples)
! update pTriples if it has a reasonable value
if (.not. t_exclude_3_body_excits .and. pTriples_new > 1e-5_dp &
.and. pTriples_new < (1.0_dp - 1e-5_dp)) then
if (abs(pTriples_new - pTriples) > 0.001_dp) then
root_print "Updating triples bias. pTriples = ", pTriples_new
end if
pTriples = pTriples_new
end if
end if
else
psingles_new = pSingles
if (abs(ratio_singles) > EPS .or. abs(ratio_doubles) > EPS) then
tau_new = MaxWalkerBloom * &
min(pSingles / max(EPS, ratio_singles), pDoubles / max(EPS, ratio_doubles))
else
tau_new = tau
end if
end if
end if
end if
! to deatch check again and finally update time-step
! The range of tau is restricted by particle death. It MUST be <=
! the value obtained to restrict the maximum death-factor to 1.0.
call MPIAllReduce(max_death_cpt, MPI_MAX, mpi_tmp)
max_death_cpt = mpi_tmp
! again, only count deaths if any occured
if(abs(max_death_cpt) > EPS) then
tau_death = 1.0_dp / max_death_cpt
if (tau_death < tau_new) then
if (tau_death < min_tau) then
root_print "time-step reduced, due to death events! reset min_tau to:", tau_death
min_tau = tau_death
end if
tau_new = tau_death
end if
end if
! If the calculated tau is less than the current tau, we should ALWAYS
! update it. Once we have a reasonable sample of excitations, then we
! can permit tau to increase if we have started too low.
! make the right if-statements here.. and remember enough_doub_hist is
! used for singles in the case of the real-space transcorrelated hubbard!
tau_new = clamp(tau_new, min_tau, max_tau)
if (tau_new < tau .or. &
(tUEG .or. tHub .or. enough_sing_hist .or. &
(t_k_space_hubbard .and. .not. t_trans_corr_2body) .and. enough_doub_hist) .or. &
(t_new_real_space_hubbard .and. (enough_doub_hist .and. &
(.not. t_trans_corr_hop .or. enough_sing_hist)))) then
! remove this constriction to just work in the transcorrelated
! case and let the user decide!
! Make the final tau smaller than tau_new by a small amount
! so that we don't get spawns exactly equal to the
! initiator threshold, but slightly below it instead.
! does this make sense in the new implmentation? this way
! i will always decrease the time-step even if its not necessary..
call assign_value_to_tau(tau_new, 'Histogramming tau search')
end if
end subroutine update_tau_hist