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