subroutine fill_frequency_histogram_4ind(mat_ele, pgen, ic, t_parallel, ex)
! this is the specific routine to fill up the frequency histograms
! for the 4ind-weighted excitation generators, which use pParallel too
! (do they always by default?? check that! todo!
! i just realised, due to the linear and ordered bins, i actually dont
! need to binary search in them but can determine the index just
! through the ratio and frequency step size
! i want to stop filling the histograms after a certain number of
! excitations is stored. 1.: since i dont want integer overflows and
! i think it is unnecessary to store it as an 64bit integer array
! and 2.: the distribution shouldn't change too much after a
! certain point..
! but this would also mean, that the tau-update would not give
! any new values after i dont change the histograms anymore
! so i could stop the tau-adaptation after that point too..
! so maybe first experiment a bit with 64bit array to see if
! the time does change more after the 32 bit limit is reached
! and then change it back to 32 bits..
real(dp), intent(in) :: mat_ele, pgen
integer, intent(in) :: ic
logical, intent(in) :: t_parallel
integer, intent(in), optional :: ex(2, ic)
character(*), parameter :: this_routine = "fill_frequency_histogram_4ind"
! i think in my histogramming tau-search i have to increase this
! threshold by a LOT to avoid fluctuating behavior at the
! beginning of the calculation
integer, parameter :: cnt_threshold = 50
real(dp) :: ratio
integer :: ind
integer :: indi, indj, inda, indb
#ifdef DEBUG_
if (pgen < EPS) then
print *, "pgen: ", pgen
print *, "matrix element: ", mat_ele
end if
#endif
if (pgen < EPS) return
ASSERT(pgen > 0.0_dp)
ASSERT(ic == 1 .or. ic == 2 .or. (ic == 3 .and. t_3_body_excits))
! i can't make this exception here without changing alot in the
! other parts of the code.. and i have to talk to ali first about
! that
! but with my cutoff change i should change this here to keep
! track of the excitations above matele cutoff or?? hm..
! and in the rest of the code i have to abort these excitations really!
! talk to ali about that!
if (mat_ele < matele_cutoff) then
! but i still should keep track of these events!
select case (ic)
case (1)
! what if i get an int overflow here? should i also
! stop histogramming?
zero_singles = zero_singles + 1
case (2)
if (t_parallel) then
zero_para = zero_para + 1
else
zero_anti = zero_anti + 1
end if
case (3)
! adapt this to possible triples now
ASSERT(t_3_body_excits)
! but still re-use the singles counters in this case
zero_singles = zero_singles + 1
end select
return
end if
if (mat_ele < thresh) then
! maybe it would be better to measure if the ratio is
! below the thresh
select case (ic)
case (1)
below_thresh_singles = below_thresh_singles + 1
case (2)
if (t_parallel) then
below_thresh_para = below_thresh_para + 1
else
below_thresh_anti = below_thresh_anti + 1
end if
case (3)
ASSERT(t_3_body_excits)
below_thresh_singles = below_thresh_singles + 1
end select
end if
ratio = mat_ele / pgen
! then i have to decide which histogram to fill
if (ic == 1 .or. (ic == 3 .and. t_3_body_excits)) then
! if i ignore the ratios above the upper limit i also can
! only count these excitations if they do not get ignored..
! i have to change the way how we store those excitations in the
! histograms! i have to unbias it against the psingles, etc.
! quantities, so the histograms do not have feedback!
! i have to ensure pSingles is also set to 1.0_dp - pDoubles
! in the transcorr hubbard case..
ratio = ratio * pSingles
if (ratio < max_frequency_bound) then
! also keep track of the number of done excitations
if (.not. enough_sing_hist) then
cnt_sing_hist = cnt_sing_hist + 1
if (cnt_sing_hist > cnt_threshold) enough_sing_hist = .true.
end if
! find where to put the ratio
! since ordered and linear bin bounds i can just divide..
! i just hope there are no numerical errors creeping in ..
! if the ratio happens to be exactly the upper bound, it
! still has to be put in the last bin or? yes i guess..
! so take the minimum of the ind and the max bin index
ind = int(ratio / frq_step_size) + 1
frequency_bins_singles(ind) = frequency_bins_singles(ind) + 1
! for now also test if the actually ratio is below the
! threshold, although.. i already now the minumum ratio,
! which is not so small.. just too many of those happen!
! we have to avoid that!
else
! store the number of excitation which exceed the upper limit!
above_max_singles = above_max_singles + 1
print *, "Warning: single excitation H_ij/pgen above max_frequency_bound!"
print *, " H_ij: ", mat_ele, ", pgen: ", pgen, ", pSingles: ", pSingles
print *, "ex-maxtrix: ", get_src(ex), " -> ", get_tgt(ex)
print *, " H_ij/pgen: ", ratio, " ; bound: ", max_frequency_bound
print *, " Consider increasing the bound!"
#ifdef DEBUG_
indi = gtid(ex(1, 1))
indj = gtid(ex(1, 2))
inda = gtid(ex(2, 1))
indb = gtid(ex(2, 2))
print *, "umat (ij|ab) ", get_umat_el(indi, indj, inda, indb)
print *, "umat (ij|ba) ", get_umat_el(indi, indj, indb, inda)
print *, "diff: ", abs(get_umat_el(indi, indj, inda, indb) - &
get_umat_el(indi, indj, indb, inda))
print *, "(ia|ia): ", abs(get_umat_el(indi, inda, indi, inda))
print *, "(ja|ja): ", abs(get_umat_el(indj, inda, indj, inda))
print *, "(ib|ib): ", abs(get_umat_el(indi, indb, indi, indb))
print *, "(jb|jb): ", abs(get_umat_el(indj, indb, indj, indb))
print *, "******************"
#endif
end if
! also start to store the maximum values anyway..
if (ratio > gamma_sing) gamma_sing = ratio
! also start to store the smallest allowed ratio:
if (ratio < min_sing) min_sing = ratio
else
! check if parallel or anti-parallel
if (t_parallel) then
ratio = ratio * (pDoubles * pParallel)
#ifdef DEBUG_
! analyse the really low and really high ratios:
if (ratio < 0.001_dp) then
print *, "******************"
print *, "parallel excitation:"
print *, "ratio: ", ratio
print *, "mat_ele: ", mat_ele
print *, "pgen: ", pgen
print *, "ex-maxtrix: ", get_src(ex), " -> ", get_tgt(ex)
indi = gtid(ex(1, 1))
indj = gtid(ex(1, 2))
inda = gtid(ex(2, 1))
indb = gtid(ex(2, 2))
print *, "umat (ij|ab) ", get_umat_el(indi, indj, inda, indb)
print *, "umat (ij|ba) ", get_umat_el(indi, indj, indb, inda)
print *, "diff: ", abs(get_umat_el(indi, indj, inda, indb) - &
get_umat_el(indi, indj, indb, inda))
print *, "(ia|ia): ", abs(get_umat_el(indi, inda, indi, inda))
print *, "(ja|ja): ", abs(get_umat_el(indj, inda, indj, inda))
print *, "(ib|ib): ", abs(get_umat_el(indi, indb, indi, indb))
print *, "(jb|jb): ", abs(get_umat_el(indj, indb, indj, indb))
print *, "******************"
end if
#endif
if (ratio < max_frequency_bound) then
if (.not. enough_par_hist) then
cnt_par_hist = cnt_par_hist + 1
if (cnt_par_hist > 2 * cnt_threshold) enough_par_hist = .true.
end if
ind = int(ratio / frq_step_size) + 1
frequency_bins_para(ind) = frequency_bins_para(ind) + 1
else
above_max_para = above_max_para + 1
#ifdef DEBUG_
print *, "mat_ele: ", mat_ele
print *, "pgen: ", pgen
print *, "ex-maxtrix: ", get_src(ex), " -> ", get_tgt(ex)
indi = gtid(ex(1, 1))
indj = gtid(ex(1, 2))
inda = gtid(ex(2, 1))
indb = gtid(ex(2, 2))
print *, "umat (ij|ab) ", get_umat_el(indi, indj, inda, indb)
print *, "umat (ij|ba) ", get_umat_el(indi, indj, indb, inda)
print *, "diff: ", abs(get_umat_el(indi, indj, inda, indb) - &
get_umat_el(indi, indj, indb, inda))
print *, "(ia|ia): ", abs(get_umat_el(indi, inda, indi, inda))
print *, "(ja|ja): ", abs(get_umat_el(indj, inda, indj, inda))
print *, "(ib|ib): ", abs(get_umat_el(indi, indb, indi, indb))
print *, "(jb|jb): ", abs(get_umat_el(indj, indb, indj, indb))
print *, "******************"
#endif
end if
if (ratio > gamma_par) gamma_par = ratio
if (ratio < min_par) min_par = ratio
else
ratio = ratio * (pDoubles * (1.0_dp - pParallel))
! analyse the really low and really high ratios:
#ifdef DEBUG_
if (ratio < 0.001_dp) then
print *, "******************"
print *, "anti-parallel excitation:"
print *, "ratio: ", ratio
print *, "mat_ele: ", mat_ele
print *, "pgen: ", pgen
print *, "ex-maxtrix: ", get_src(ex), " -> ", get_tgt(ex)
indi = gtid(ex(1, 1))
indj = gtid(ex(1, 2))
inda = gtid(ex(2, 1))
indb = gtid(ex(2, 2))
print *, "umat (ij|ab) ", get_umat_el(indi, indj, inda, indb)
print *, "umat (ij|ba) ", get_umat_el(indi, indj, indb, inda)
print *, "diff: ", abs(get_umat_el(indi, indj, inda, indb) - &
get_umat_el(indi, indj, indb, inda))
print *, "(ia|ia): ", abs(get_umat_el(indi, inda, indi, inda))
print *, "(ja|ja): ", abs(get_umat_el(indj, inda, indj, inda))
print *, "(ib|ib): ", abs(get_umat_el(indi, indb, indi, indb))
print *, "(jb|jb): ", abs(get_umat_el(indj, indb, indj, indb))
print *, "******************"
end if
#endif
if (ratio < max_frequency_bound) then
if (.not. enough_opp_hist) then
cnt_opp_hist = cnt_opp_hist + 1
if (cnt_opp_hist > cnt_threshold) enough_opp_hist = .true.
end if
ind = int(ratio / frq_step_size) + 1
frequency_bins_anti(ind) = frequency_bins_anti(ind) + 1
else
above_max_anti = above_max_anti + 1
#ifdef DEBUG_
print *, "mat_ele: ", mat_ele
print *, "pgen: ", pgen
print *, "ex-maxtrix: ", get_src(ex), " -> ", get_tgt(ex)
indi = gtid(ex(1, 1))
indj = gtid(ex(1, 2))
inda = gtid(ex(2, 1))
indb = gtid(ex(2, 2))
print *, "umat (ij|ab) ", get_umat_el(indi, indj, inda, indb)
print *, "umat (ij|ba) ", get_umat_el(indi, indj, indb, inda)
print *, "diff: ", abs(get_umat_el(indi, indj, inda, indb) - &
get_umat_el(indi, indj, indb, inda))
print *, "(ia|ia): ", abs(get_umat_el(indi, inda, indi, inda))
print *, "(ja|ja): ", abs(get_umat_el(indj, inda, indj, inda))
print *, "(ib|ib): ", abs(get_umat_el(indi, indb, indi, indb))
print *, "(jb|jb): ", abs(get_umat_el(indj, indb, indj, indb))
print *, "******************"
#endif
end if
if (ratio > gamma_opp) gamma_opp = ratio
if (ratio < min_opp) min_opp = ratio
end if
! combine par and opp into enough_doub
if (enough_par_hist .and. enough_opp_hist) enough_doub_hist = .true.
end if
end subroutine fill_frequency_histogram_4ind