print_frequency_histograms Subroutine

    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(histogramming_spec%n_frequency_bins)
        integer(int64) :: all_frequency_bins(histogramming_spec%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(histogramming_spec%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 = histogramming_spec%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, histogramming_spec%n_frequency_bins
                    if (all_frequency_bins(i) == 0) cycle
                    write(iunit, "(f16.7)", advance='no') histogramming_spec%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 (histogramming_spec%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 * histogramming_spec%frq_step_size
                write(stdout, *) "for ", histogramming_spec%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 * histogramming_spec%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 = histogramming_spec%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, histogramming_spec%n_frequency_bins
                        if (all_frequency_bins_spec(i) == 0) cycle
                        write(iunit, "(f16.7)", advance='no') histogramming_spec%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 (histogramming_spec%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 * histogramming_spec%frq_step_size
                    write(stdout, *) "for ", histogramming_spec%frq_ratio_cutoff, " percent coverage!"

                    write(stdout, *) "maximum parallel H_ij/pgen ratio: ", max_tmp
                    write(stdout, *) "maximum/integrated parallel ratio: ", &
                        max_tmp / (j * histogramming_spec%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 = histogramming_spec%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, histogramming_spec%n_frequency_bins
                        if (all_frequency_bins_spec(i) == 0) cycle
                        write(iunit, "(f16.7)", advance='no') histogramming_spec%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 (histogramming_spec%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 * histogramming_spec%frq_step_size
                    write(stdout, *) "for ", histogramming_spec%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 * histogramming_spec%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 = histogramming_spec%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, histogramming_spec%n_frequency_bins
                        if (all_frequency_bins_spec(i) == 0) cycle
                        write(iunit, "(f16.7)", advance='no') histogramming_spec%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 (histogramming_spec%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 * histogramming_spec%frq_step_size
                    write(stdout, *) "for ", histogramming_spec%frq_ratio_cutoff, " percent coverage!"

                    write(stdout, *) "maximum doubles H_ij/pgen ratio: ", max_tmp
                    write(stdout, *) "maximum/integrated doubles ratio: ", max_tmp / (j * histogramming_spec%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, histogramming_spec%n_frequency_bins
                    ! only print if above threshold
                    if (temp_bins(i) / sum_all < EPS) cycle
                    write(iunit, "(f16.7)", advance='no') histogramming_spec%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 = histogramming_spec%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, histogramming_spec%n_frequency_bins
                    if (all_frequency_bins_spec(i) == 0) cycle
                    write(iunit, "(f16.7)", advance='no') histogramming_spec%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 (histogramming_spec%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 * histogramming_spec%frq_step_size
                write(stdout, *) "for ", histogramming_spec%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 * histogramming_spec%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