subroutine print_frequency_histograms
! try to write one general one and not multiple as in my GUGA branch
use constants, only: int64
use util_mod, only: get_free_unit, get_unique_filename
use MPI_wrapper, only: root
character(255) :: filename, exname
integer(int64) :: all_frequency_bins_spec(n_frequency_bins)
integer(int64) :: all_frequency_bins(n_frequency_bins)
integer :: iunit, i
! sashas tip: why do i not just use a real as summation?
real(dp) :: sum_all
integer :: tmp_int, j, k
integer(int64) :: tmp_int_64
real(dp) :: cnt, threshold
real(dp) :: max_tmp, min_tmp
real(dp) :: temp_bins(n_frequency_bins)
all_frequency_bins = 0
! also do additional output here.. to get information about the
! number of ratios above the threshold and about the number of
! zero excitations and stuff
! maybe first check if we have only singles or only doubles like in
! the real-space or momentum space hubbard:
if (tHub .or. tUEG .or. &
(t_new_real_space_hubbard .and. .not. t_trans_corr_hop) .or. &
(t_k_space_hubbard .and. .not. t_3_body_excits)) then
! we only need to print one frequency_histogram:
! i need to change the rest of the code to use this
! frequency_bins histogram in this case! TODO
call MPIAllReduce(frequency_bins, MPI_SUM, all_frequency_bins)
max_tmp = 0.0_dp
min_tmp = huge(0.0_dp)
call mpiallreduce(gamma_doub, MPI_MAX, max_tmp)
call MPIAllReduce(min_doub, MPI_MIN, min_tmp)
if (iProcIndex == root) then
! also outout the obtained ratio integration here for now!
temp_bins = real(all_frequency_bins, dp)
sum_all = sum(temp_bins)
threshold = frq_ratio_cutoff * sum_all
iunit = get_free_unit()
call get_unique_filename('frequency_histogram', .true., &
.true., 1, filename)
open(iunit, file=filename, status='unknown')
cnt = 0.0_dp
write(stdout, *) "writing frequency histogram..."
do i = 1, n_frequency_bins
if (all_frequency_bins(i) == 0) cycle
write(iunit, "(f16.7)", advance='no') frq_step_size * i
write(iunit, "(i12)") all_frequency_bins(i)
if (cnt < threshold) then
cnt = cnt + temp_bins(i)
j = i
end if
! also print only as far as necessary until largest
! H_ij/pgen ratio
if (frq_step_size * i > max_tmp) exit
end do
close(iunit)
write(stdout, *) "Done!"
end if
tmp_int_64 = 0
call mpisum(zero_doubles, tmp_int_64)
if (iProcIndex == root) then
write(stdout, *) "Number of zero-valued excitations: ", tmp_int_64
write(stdout, *) "Number of valid excitations: ", sum_all
write(stdout, *) "ratio of zero-valued excitations: ", &
real(tmp_int_64, dp) / sum_all
! i guess i should also output the number of excitations
! above the threshold!
! this is not really working..
! because sasha had these cases
end if
tmp_int = 0
call mpisum(above_max_doubles, tmp_int)
if (iProcIndex == root) then
write(stdout, *) "Number of excitations above threshold: ", tmp_int
write(stdout, *) "ratio of excitations above threshold: ", &
real(tmp_int, dp) / sum_all
! also output the obtained integrated threshold and the
! maximum values H_ij/pgen ratio for this type of excitation:
write(stdout, *) "integrated H_ij/pgen ratio: ", j * frq_step_size
write(stdout, *) "for ", frq_ratio_cutoff, " percent coverage!"
! also output the maximum H_ij/pgen ratio.. and maybe the
! improvement of the integrated tau search?
write(stdout, *) "maximum H_ij/pgen ratio: ", max_tmp
write(stdout, *) "maximum/integrated ratio: ", max_tmp / (j * frq_step_size)
write(stdout, *) "minimum H_ij/pgen ratio: ", min_tmp
end if
else
if (iProcIndex == root) all_frequency_bins = 0
! in the other cases we definetly have singles and more
! first the singles:
if (iProcIndex == root) then
! change the name in case of the 2-body transcorrelated k-space hubbard
if (t_3_body_excits .and. .not. t_mol_3_body) then
call get_unique_filename('frequency_histogram_triples', .true., &
.true., 1, filename)
else
call get_unique_filename('frequency_histogram_singles', .true., &
.true., 1, filename)
end if
end if
exname = "singles"
call output_histogram(exname, frequency_bins_singles, gamma_sing, &
min_sing, zero_singles, above_max_singles, below_thresh_singles)
all_frequency_bins_spec = 0
! output the three-body histogram (if existing)
if (t_mol_3_body) then
call get_unique_filename('frequency_histogram_triples', .true., &
.true., 1, filename)
! output the triples in the same fashion as the singles
exname = "triples"
call output_histogram(exname, frequency_bins_triples, gamma_trip, &
min_trip, zero_triples, above_max_triples, below_thresh_triples)
all_frequency_bins_spec = 0
end if
! do the cases where there is antiparallel or parallel
if (t_consider_par_bias) then
! then the parallel:
call MPIAllReduce(frequency_bins_para, MPI_SUM, all_frequency_bins_spec)
max_tmp = 0.0_dp
min_tmp = huge(0.0_dp)
call mpiallreduce(gamma_par, MPI_MAX, max_tmp)
call MPIAllReduce(min_par, MPI_MIN, min_tmp)
if (iProcIndex == root) then
temp_bins = real(all_frequency_bins_spec, dp)
sum_all = sum(temp_bins)
threshold = frq_ratio_cutoff * sum_all
iunit = get_free_unit()
call get_unique_filename('frequency_histogram_para', .true., &
.true., 1, filename)
open(iunit, file=filename, status='unknown')
cnt = 0.0_dp
write(stdout, *) "writing parallel frequency histogram..."
do i = 1, n_frequency_bins
if (all_frequency_bins_spec(i) == 0) cycle
write(iunit, "(f16.7)", advance='no') frq_step_size * i
write(iunit, "(i12)") all_frequency_bins_spec(i)
if (cnt < threshold) then
cnt = cnt + temp_bins(i)
j = i
end if
if (frq_step_size * i > max_tmp) exit
end do
close(iunit)
write(stdout, *) "Done!"
end if
tmp_int_64 = 0
call MPISum(zero_para, tmp_int_64)
if (iProcIndex == root) then
write(stdout, *) "Number of zero-valued parallel excitations: ", tmp_int_64
write(stdout, *) "Number of valid parallel excitations: ", sum_all
write(stdout, *) "ratio of zero-valued parallel excitations: ", &
real(tmp_int_64, dp) / sum_all
end if
tmp_int = 0
call mpisum(above_max_para, tmp_int)
if (iProcIndex == root) then
write(stdout, *) "Number of parallel excitations above threshold: ", tmp_int
write(stdout, *) "ratio of parallel excitations above threshold: ", &
real(tmp_int, dp) / sum_all
end if
tmp_int = 0
call MPISum(below_thresh_para, tmp_int)
if (iProcIndex == root) then
write(stdout, *) "Number of parallel excitations below threshold: ", tmp_int
write(stdout, *) "ratio of parallel excitations below threshold: ", &
real(tmp_int, dp) / sum_all
write(stdout, *) "integrated parallel H_ij/pgen ratio: ", j * frq_step_size
write(stdout, *) "for ", frq_ratio_cutoff, " percent coverage!"
write(stdout, *) "maximum parallel H_ij/pgen ratio: ", max_tmp
write(stdout, *) "maximum/integrated parallel ratio: ", &
max_tmp / (j * frq_step_size)
write(stdout, *) "minimum parallel H_ij/pgen ratio: ", min_tmp
! and add them up for the final normed one
all_frequency_bins = all_frequency_bins + all_frequency_bins_spec
end if
all_frequency_bins_spec = 0
! then anti:
call MPIAllReduce(frequency_bins_anti, MPI_SUM, all_frequency_bins_spec)
max_tmp = 0.0_dp
min_tmp = huge(0.0_dp)
call mpiallreduce(gamma_opp, MPI_MAX, max_tmp)
call MPIAllReduce(min_opp, MPI_MIN, min_tmp)
if (iProcIndex == root) then
temp_bins = real(all_frequency_bins_spec, dp)
sum_all = sum(temp_bins)
threshold = frq_ratio_cutoff * sum_all
iunit = get_free_unit()
call get_unique_filename('frequency_histogram_anti', .true., &
.true., 1, filename)
open(iunit, file=filename, status='unknown')
cnt = 0.0_dp
write(stdout, *) "writing anti-parallel frequency histogram..."
do i = 1, n_frequency_bins
if (all_frequency_bins_spec(i) == 0) cycle
write(iunit, "(f16.7)", advance='no') frq_step_size * i
write(iunit, "(i12)") all_frequency_bins_spec(i)
if (cnt < threshold) then
cnt = cnt + temp_bins(i)
j = i
end if
if (frq_step_size * i > max_tmp) exit
end do
close(iunit)
write(stdout, *) "Done!"
end if
tmp_int_64 = 0
call MPISum(zero_anti, tmp_int_64)
if (iProcIndex == root) then
write(stdout, *) "Number of zero-valued anti-parallel excitations: ", tmp_int_64
write(stdout, *) "Number of valid anti-parallel excitations: ", sum_all
write(stdout, *) "ratio of zero-valued anti-parallel excitations: ", &
real(tmp_int_64, dp) / sum_all
end if
tmp_int = 0
call mpisum(above_max_anti, tmp_int)
if (iProcIndex == root) then
write(stdout, *) "Number of anti-parallel excitations above threshold: ", &
tmp_int
write(stdout, *) "ratio of anti-parallel excitations above threshold: ", &
real(tmp_int, dp) / sum_all
end if
tmp_int = 0
call mpisum(below_thresh_anti, tmp_int)
if (iProcIndex == root) then
write(stdout, *) "Number of anti-parallel excitations below threshold: ", &
tmp_int
write(stdout, *) "ratio of anti-parallel excitations below threshold: ", &
real(tmp_int, dp) / sum_all
write(stdout, *) "integrated anti-parallel H_ij/pgen ratio: ", j * frq_step_size
write(stdout, *) "for ", frq_ratio_cutoff, " percent coverage!"
write(stdout, *) "maximum anti-parallel H_ij/pgen ratio: ", max_tmp
write(stdout, *) "maximum/integrated anti-parallel ratio: ", &
max_tmp / (j * frq_step_size)
write(stdout, *) "minimum anti-parallel H_ij/pgen ratio: ", min_tmp
! and add them up for the final normed one
all_frequency_bins = all_frequency_bins + all_frequency_bins_spec
end if
all_frequency_bins_spec = 0
else
! we only have additional doubles:
call MPIAllReduce(frequency_bins_doubles, MPI_SUM, all_frequency_bins_spec)
max_tmp = 0.0_dp
min_tmp = huge(0.0_dp)
call mpiallreduce(gamma_doub, MPI_MAX, max_tmp)
call MPIAllReduce(min_doub, MPI_MIN, min_tmp)
if (iProcIndex == root) then
temp_bins = real(all_frequency_bins_spec, dp)
sum_all = sum(temp_bins)
threshold = frq_step_size * sum_all
iunit = get_free_unit()
call get_unique_filename('frequency_histogram_doubles', .true., &
.true., 1, filename)
open(iunit, file=filename, status='unknown')
cnt = 0.0_dp
write(stdout, *) "writing doubles frequency histogram..."
do i = 1, n_frequency_bins
if (all_frequency_bins_spec(i) == 0) cycle
write(iunit, "(f16.7)", advance='no') frq_step_size * i
write(iunit, "(i12)") all_frequency_bins_spec(i)
if (cnt < threshold) then
cnt = cnt + temp_bins(i)
j = i
end if
if (frq_step_size * i > max_tmp) exit
end do
close(iunit)
write(stdout, *) "Done!"
end if
tmp_int_64 = 0
call MPISUM(zero_doubles, tmp_int_64)
if (iprocindex == root) then
write(stdout, *) "Number of zero-valued double excitations: ", tmp_int_64
write(stdout, *) "Number of valid double excitations: ", sum_all
write(stdout, *) "ratio of zero-valued double excitations: ", &
real(tmp_int_64, dp) / sum_all
end if
tmp_int = 0
call mpisum(above_max_doubles, tmp_int)
if (iprocindex == root) then
write(stdout, *) "Number of excitations above threshold: ", tmp_int
write(stdout, *) "ratio of excitations above threshold: ", &
real(tmp_int, dp) / sum_all
end if
tmp_int = 0
call MPISUM(below_thresh_doubles, tmp_int)
if (iprocindex == root) then
write(stdout, *) "Number of excitations below threshold: ", tmp_int
write(stdout, *) "ratio of excitations below threshold: ", &
real(tmp_int, dp) / sum_all
write(stdout, *) "integrated doubles H_ij/pgen ratio: ", j * frq_step_size
write(stdout, *) "for ", frq_ratio_cutoff, " percent coverage!"
write(stdout, *) "maximum doubles H_ij/pgen ratio: ", max_tmp
write(stdout, *) "maximum/integrated doubles ratio: ", max_tmp / (j * frq_step_size)
write(stdout, *) "minimum doubles H_ij/pgen ratio: ", min_tmp
! and add them up for the final normed one
all_frequency_bins = all_frequency_bins + all_frequency_bins_spec
end if
end if
end if
! and the norm then is the same for all cases:
if (iprocindex == root) then
! also print out a normed version for comparison
temp_bins = real(all_frequency_bins, dp)
sum_all = sum(temp_bins)
! check if integer overflow:
if (.not. sum_all < 0.0_dp) then
iunit = get_free_unit()
call get_unique_filename('frequency_histogram_normed', .true., &
.true., 1, filename)
open(iunit, file=filename, status='unknown')
do i = 1, n_frequency_bins
! only print if above threshold
if (temp_bins(i) / sum_all < EPS) cycle
write(iunit, "(f16.7)", advance='no') frq_step_size * i
write(iunit, "(f16.7)") temp_bins(i) / sum_all
end do
close(iunit)
else
write(stdout, *) "Integer overflow in normed frequency histogram!"
write(stdout, *) "DO NOT PRINT IT!"
end if
end if
! why am i misusing the hist-tau-search also for these method?
! because i want to have even more info on the dead-end excitations!
if (t_log_ija) then
! i hope this mpiallreduce works..
allocate(all_ija_bins_sing(nBasis / 2))
all_ija_bins_sing = 0
call MPIAllReduce(ija_bins_sing, MPI_SUM, all_ija_bins_sing)
allocate(all_ija_orbs_sing(nBasis / 2))
all_ija_orbs_sing = 0
call MPIAllReduce(ija_orbs_sing, MPI_MAX, all_ija_orbs_sing)
if (iprocindex == root) then
iunit = get_free_unit()
call get_unique_filename('ija_bins_sing', .true., .true., 1, filename)
open(iunit, file=filename, status='unknown')
do i = 1, nBasis / 2
if (all_ija_bins_sing(i) == 0) cycle
write(iunit, '(i12, i6, i12)') all_ija_bins_sing(i), i, all_ija_orbs_sing(i)
end do
close(iunit)
end if
deallocate(all_ija_bins_sing)
deallocate(all_ija_orbs_sing)
deallocate(ija_bins_sing)
deallocate(ija_orbs_sing)
allocate(all_ija_bins(nBasis / 2, nBasis / 2, nBasis / 2))
allocate(all_ija_orbs(nBasis / 2, nBasis / 2, nBasis / 2))
all_ija_bins = 0
all_ija_orbs = 0
do i = 1, nBasis / 2
call MPIAllReduce(ija_bins_para(i, :, :), MPI_SUM, all_ija_bins(i, :, :))
! can i also communicate the number of symmetry allowed orbitals?
! and can i store it in spatial orbital information?
! i could make seperate ones for parallel and anti-parallel
! excitations.. and also for singles or?
call MPIAllReduce(ija_orbs_para(i, :, :), MPI_MAX, all_ija_orbs(i, :, :))
end do
if (iprocindex == root) then
iunit = get_free_unit()
call get_unique_filename('ija_bins_para', .true., .true., 1, filename)
open(iunit, file=filename, status='unknown')
! i know it is slower but loop over first row to get it
! nicely ordered in the output
! and change the number of occurences to the first line
! so it can be easily sorted with bash..
do i = 1, nBasis / 2
do j = 1, nbasis / 2
do k = 1, nBasis / 2
if (all_ija_bins(i, j, k) == 0) cycle
write(iunit, '(i12,i6, 2i3, i12)') all_ija_bins(i, j, k), i, j, k, all_ija_orbs(i, j, k)
end do
end do
end do
close(iunit)
end if
deallocate(ija_bins_para)
deallocate(ija_orbs_para)
all_ija_bins = 0
all_ija_orbs = 0
do i = 1, nBasis / 2
call MPIAllReduce(ija_bins_anti(i, :, :), MPI_SUM, all_ija_bins(i, :, :))
! can i also communicate the number of symmetry allowed orbitals?
! and can i store it in spatial orbital information?
! i could make seperate ones for parallel and anti-parallel
! excitations.. and also for singles or?
call MPIAllReduce(ija_orbs_anti(i, :, :), MPI_MAX, all_ija_orbs(i, :, :))
end do
if (iprocindex == root) then
iunit = get_free_unit()
call get_unique_filename('ija_bins_anti', .true., .true., 1, filename)
open(iunit, file=filename, status='unknown')
do i = 1, nBasis / 2
do j = 1, nbasis / 2
do k = 1, nBasis / 2
if (all_ija_bins(i, j, k) == 0) cycle
write(iunit, '(i12, i6, 2i3, i12)') all_ija_bins(i, j, k), i, j, k, all_ija_orbs(i, j, k)
end do
end do
end do
close(iunit)
end if
deallocate(ija_bins_anti)
deallocate(ija_orbs_anti)
end if
contains
subroutine output_histogram(histname, frequency_bins_loc, gamma_loc, min_loc, &
zero_loc, above_max_loc, below_thresh_loc)
implicit none
character(255), intent(in) :: histname
integer(int64), intent(in) :: frequency_bins_loc(:)
real(dp), intent(in) :: gamma_loc, min_loc
integer(int64), intent(in) :: zero_loc
integer, intent(in) :: above_max_loc, below_thresh_loc
all_frequency_bins_spec = 0
call MPIAllReduce(frequency_bins_loc, MPI_SUM, all_frequency_bins_spec)
max_tmp = 0.0_dp
min_tmp = huge(0.0_dp)
call mpiallreduce(gamma_loc, MPI_MAX, max_tmp)
call MPIAllReduce(min_loc, MPI_MIN, min_tmp)
if (iProcIndex == root) then
temp_bins = real(all_frequency_bins_spec, dp)
sum_all = sum(temp_bins)
threshold = frq_ratio_cutoff * sum_all
iunit = get_free_unit()
write(stdout, *) "writing "//trim(histname)//" frequency histogram..."
open(iunit, file=filename, status='unknown')
cnt = 0.0_dp
! also output here the integrated ratio!
do i = 1, n_frequency_bins
if (all_frequency_bins_spec(i) == 0) cycle
write(iunit, "(f16.7)", advance='no') frq_step_size * i
write(iunit, "(i12)") all_frequency_bins_spec(i)
if (cnt < threshold) then
! cnt = cnt + all_frequency_bins_spec(i)
cnt = cnt + temp_bins(i)
j = i
end if
if (frq_step_size * i > max_tmp) exit
end do
close(iunit)
write(stdout, *) "Done!"
end if
tmp_int_64 = 0
call mpisum(zero_loc, tmp_int_64)
if (iProcIndex == root) then
write(stdout, *) "Number of zero-valued "//trim(histname)//" excitations: ", tmp_int_64
! maybe also check the number of valid excitations
write(stdout, *) "Number of valid "//trim(histname)//" excitations: ", sum_all
if (abs(sum_all) > 0._dp) write(stdout, *) "ratio of zero-valued "//trim(histname)//" excitations: ", &
real(tmp_int_64, dp) / sum_all
end if
tmp_int = 0
call mpisum(above_max_loc, tmp_int)
if (iProcIndex == root) then
write(stdout, *) "Number of "//trim(histname)//" excitations above threshold: ", tmp_int
if (abs(sum_all) > 0._dp) write(stdout, *) "ratio of "//trim(histname)//" excitations above threshold: ", &
real(tmp_int, dp) / sum_all
end if
tmp_int = 0
call MPISum(below_thresh_loc, tmp_int)
if (iProcIndex == root) then
write(stdout, *) "Number of "//trim(histname)//" excitations below threshold: ", tmp_int
if (abs(sum_all) > 0._dp) write(stdout, *) "ratio of "//trim(histname)//" excitations below threshold: ", &
real(tmp_int, dp) / sum_all
write(stdout, *) "integrated "//trim(histname)//" H_ij/pgen ratio: ", j * frq_step_size
write(stdout, *) "for ", frq_ratio_cutoff, " percent coverage!"
write(stdout, *) "maximum "//trim(histname)//" H_ij/pgen ratio: ", max_tmp
write(stdout, *) ""//trim(histname)//" maximum/integrated ratio: ", max_tmp / (j * frq_step_size)
write(stdout, *) "minimum "//trim(histname)//" H_ij/pgen ratio: ", min_tmp
! and add them up for the final normed one
all_frequency_bins = all_frequency_bins + all_frequency_bins_spec
end if
end subroutine output_histogram
end subroutine print_frequency_histograms