subroutine fill_frequency_histogram_sd(mat_ele, pgen, ic)
! this is the routine, which for now i can use in the symmetry
! adapted GUGA code, where i only distinguish between single and
! double excitation for now..
! if i adapt that to the already used method in nosym_guga to
! distinguish between more excitaiton i have to write an additional
! subroutine which deals with that
! and this can also be used in the case where we dont differentiate
! between parallel or antiparallel double excitations in the old
! excitation generators
real(dp), intent(in) :: mat_ele, pgen
integer, intent(in) :: ic
character(*), parameter :: this_routine = "fill_frequency_histogram_sd"
integer, parameter :: cnt_threshold = 50
real(dp) :: ratio
integer :: ind
ASSERT(pgen > EPS)
ASSERT(ic == 1 .or. ic == 2 .or. ic == 3)
if (mat_ele < matele_cutoff) then
select case (ic)
case (1)
zero_singles = zero_singles + 1
case (2)
zero_doubles = zero_doubles + 1
case (3)
if (mat_ele < lMatEps) zero_triples = zero_triples + 1
end select
return
end if
if (mat_ele < thresh) then
select case (ic)
case (1)
below_thresh_singles = below_thresh_singles + 1
case (2)
below_thresh_doubles = below_thresh_doubles + 1
case (3)
below_thresh_triples = below_thresh_triples + 1
end select
end if
if (pgen < EPS) return
ratio = mat_ele / pgen
if (t_mol_3_body .and. ic < 3) ratio = ratio * (1.0 - pTriples)
if (ic == 1) then
! change to unbias the stored ratios:
ratio = ratio * pSingles
if (ratio < max_frequency_bound) then
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
ind = int(ratio / frq_step_size) + 1
frequency_bins_singles(ind) = frequency_bins_singles(ind) + 1
else
above_max_singles = above_max_singles + 1
print *, "Warning: single excitation H_ij/pgen above max_frequency_bound!"
print *, " H_ij/pgen: ", ratio, " ; bound: ", max_frequency_bound
print *, " Consider increasing the bound!"
end if
if (ratio > gamma_sing) gamma_sing = ratio
if (ratio < min_sing) min_sing = ratio
else if (ic == 2) then
ratio = ratio * pDoubles
if (ratio < max_frequency_bound) then
if (.not. enough_doub_hist) then
cnt_doub_hist = cnt_doub_hist + 1
if (cnt_doub_hist > cnt_threshold) enough_doub_hist = .true.
end if
ind = int(ratio / frq_step_size) + 1
frequency_bins_doubles(ind) = frequency_bins_doubles(ind) + 1
else
above_max_doubles = above_max_doubles + 1
print *, "Warning: double excitation H_ij/pgen above max_frequency_bound!"
print *, " H_ij/pgen: ", ratio, " ; bound: ", max_frequency_bound
print *, " Consider increasing the bound!"
end if
if (ratio > gamma_doub) gamma_doub = ratio
if (ratio < min_doub) min_doub = ratio
else
ratio = ratio * pTriples
! check if ratio is within range
if (ratio < max_frequency_bound) then
! check if enough triple spawns have been logged
if (.not. enough_trip_hist) then
cnt_trip_hist = cnt_trip_hist + 1
if (cnt_trip_hist > cnt_threshold) enough_trip_hist = .true.
end if
! get the corresponding bin in the histogram
ind = int(ratio / frq_step_size) + 1
frequency_bins_triples(ind) = frequency_bins_triples(ind) + 1
else
! the ratio is not within range, log this event
above_max_triples = above_max_triples + 1
end if
! update largest/smallest ratio for triples
if (ratio > gamma_trip) gamma_trip = ratio
if (ratio < min_trip) min_trip = ratio
end if
end subroutine fill_frequency_histogram_sd