run_excit_gen_tester_function Subroutine

private subroutine run_excit_gen_tester_function(excit_gen, excit_gen_name, opt_nI, opt_n_dets, gen_all_excits, calc_pgen, problem_filter, i_unit, successful)

@brief Test if an excitation generator generates all and only expected states with the correct pgen.

@author Werner Dobrautz, Oskar Weser

@details The pgen_diagnostic is given by \f[\sum_i \frac{1}{p_i N}\f] and should be roughly one. All problematic states with respect to that diagnostic are printed in the end with their pgen_diagnostic, excitation type (ic), matrix element, and pgen. @param[in] excit_gen, An excitation generator. @param[in] excit_gen_name, The name of the excitation generator. @param[in] opt_nI, An optional reference state. @param[in] opt_n_dets, An optional number of valid determinants to generate. Defaults to 100000. (Does not count determinants with nJ(1) == 0!) @param[in] gen_all_excits, An optional subroutine to generate all states that can be reached from the reference state. @param[in] calc_pgen, An optional function that calculates the pgen for a given reference an excitation. If a state is never generated, the pgen cannot be taken from the excitation generator. Adds an additional column to the output table. @param[in] problem_filter, An optional predicate function. Return true, if an excitation exc from determinant det_I and a given pgen_diagnostic (sum 1/pgen) is considered to be problematic. If it returns true, an entry in the final table is printed. If there is any problematic excitation, the out parameter successful will become .false.. By default all states with a pgen_diagnostic that deviate with more than 5\,\% from 100\,\% and have nonzereo matrix element are printed.


Type IntentOptional Attributes Name
procedure(generate_excitation_t) :: excit_gen
character(len=*), intent(in) :: excit_gen_name
integer, intent(in), optional :: opt_nI(nel)
integer, intent(in), optional :: opt_n_dets
procedure(generate_all_excits_t), optional :: gen_all_excits
procedure(calc_pgen_t), optional :: calc_pgen
procedure(problem_filter_t), optional :: problem_filter
integer, intent(in), optional :: i_unit
logical, intent(out), optional :: successful


Source Code

    subroutine run_excit_gen_tester_function(&
                        excit_gen, excit_gen_name, opt_nI, opt_n_dets, &
                        gen_all_excits, calc_pgen, problem_filter, i_unit, successful)
        procedure(generate_excitation_t) :: excit_gen
        character(*), intent(in) :: excit_gen_name
        integer, intent(in), optional :: opt_nI(nel), opt_n_dets
        procedure(generate_all_excits_t), optional :: gen_all_excits
        procedure(calc_pgen_t), optional :: calc_pgen
        procedure(problem_filter_t), optional :: problem_filter
        integer, intent(in), optional :: i_unit
        logical, intent(out), optional :: successful
        character(*), parameter :: this_routine = "run_excit_gen_tester"

        integer :: i, nI(nel), n_dets_target
        integer :: i_unit_
        integer, parameter :: default_n_dets_target = 100000

        integer(n_int) :: ilut(0:niftot), tgt_ilut(0:niftot)
        integer :: nJ(nel), n_excits, ex(2, maxExcit), ic, ex_flag, i_unused = 0
        type(excit_gen_store_type) :: store
        logical :: tPar, found_all
        real(dp) :: pgen, contrib
        real(dp), allocatable :: pgen_list(:)
        HElement_t(dp) :: hel
        integer(n_int), allocatable :: det_list(:, :)
        real(dp), allocatable :: contrib_list(:)
        logical, allocatable :: generated_list(:)
        integer :: n_generated, pos, n_iter

        procedure(problem_filter_t), pointer :: problem_filter_

        ! and also nbasis and stuff..
        ASSERT(nbasis > 0)
        ASSERT(nel <= nbasis)

        if (present(i_unit)) then
            i_unit_ = i_unit
            i_unit_ = stdout
        end if

        if (present(problem_filter)) then
            problem_filter_ => problem_filter
            problem_filter_ => default_predicate
        end if

        call init_excit_gen_store(store)
        ! use some default values if not provided:
        ! nel must be set!
        if (present(opt_nI)) then
            nI = opt_nI
            ! use HF-det as default
            nI = [(i, i=1, nel)]
        end if

        if (present(opt_n_dets)) then
            n_dets_target = opt_n_dets
            n_dets_target = default_n_dets_target
        end if

        ! Calculate all possible excitations.
        ! For special systems (hubbard, UEG, GAS, etc.) the calling site
        ! has to support a function.
        if (present(gen_all_excits)) then
            call gen_all_excits(nI, n_excits, det_list)
            call gen_all_excits_default(nI, n_excits, det_list)
        end if

        allocate(generated_list(n_excits), source=.false.)
        allocate(contrib_list(n_excits), source=0.0_dp)
        allocate(pgen_list(n_excits), source=0.0_dp)

        write(i_unit_, *) "---------------------------------"
        write(i_unit_, *) "testing: ", excit_gen_name
        write(i_unit_, *) "for ", size(det_list, 2), " configurations"
        write(i_unit_, *) " and ", n_dets_target, " determinants to generate"

        call EncodeBitDet(nI, ilut)
        write(i_unit_, *) "for starting determinant: ", nI

        write(i_unit_, *) ! linebreak
        write(i_unit_, '(A)') 'Progressbar'

        ! call this below now for the number of specified determinants
        ! (also use excitations of the first inputted, to be really
        !   consistent)

            integer(int64) :: L
            type(timer) :: loop_timer
            loop_timer%timer_name = 'loop '//excit_gen_name

            L = n_dets_target .div. 100
            n_generated = 0
            n_iter = 0
            contrib = 0.0_dp
            do while (n_generated < n_dets_target)
                n_iter = n_iter + 1

                call set_timer(loop_timer)
                call excit_gen(nI, ilut, nJ, tgt_ilut, ex_flag, ic, ex, tpar, pgen, &
                               hel, store)
                call halt_timer(loop_timer)

                if (nJ(1) == 0) cycle
                call EncodeBitDet(nJ, tgt_ilut)
                pos = binary_search_ilut(det_list, tgt_ilut, nifd + 1)
                if (pos < 0) then
                    write(i_unit_, *) "nJ: ", nJ
                    write(i_unit_, *) "ilutJ:", tgt_ilut
                    call stop_all(this_routine, 'Unexpected determinant generated')
                    if (mod(n_generated, L) == 0_int64) then
                        write(i_unit_, '(A)', advance='no') '#'
                        flush (i_unit_)
                    end if

                    generated_list(pos) = .true.
                    n_generated = n_generated + 1

                    contrib = contrib + 1.0_dp / pgen
                    contrib_list(pos) = contrib_list(pos) + 1.0_dp / pgen
                    pgen_list(pos) = pgen
                end if
            end do

            write(i_unit_, *) ! linebreak
            write(i_unit_, *) ! linebreak
            write(i_unit_, *) "=================================================================="
            write(i_unit_, *) n_generated, " dets generated in ", n_iter, " iterations "
            write(i_unit_, *) n_generated, (n_iter - n_generated), n_iter, 'valid  invalid and total excitations'
            write(i_unit_, *) 100.0_dp * real(n_iter - n_generated, dp) / real(n_iter, dp), "% abortion rate"
            write(i_unit_, *) "Averaged contribution: ", contrib / real(n_excits * n_dets_target, dp)
            write(i_unit_, *) ! linebreak
            write(i_unit_, *) 'Elapsed Time in seconds:', get_total_time(loop_timer)
            write(i_unit_, *) 'Elapsed Time in micro seconds per valid excitation:', get_total_time(loop_timer) / dble(n_dets_target)  * 10.**6
            write(i_unit_, *) 'Elapsed Time in micro seconds per excitation:', get_total_time(loop_timer) / dble(n_iter)  * 10.**6
            write(i_unit_, *) "=================================================================="
            write(i_unit_, *) ! linebreak
            write(i_unit_, *) ! linebreak
        end block

            real(dp) :: pgen_diagnostic
            integer :: nJ(nEl), exc(2, maxExcit)
            logical :: tParity

            write(i_unit_, *) "=================================="
            write(i_unit_, *) "Problematic contribution List: "
            write(i_unit_, *) "=================================="
            write(i_unit_, '("|       Determinant         |   Sum{1 / pgen} / n_iter |  ic |    <psi_I H psi_J >        |    pgen    |")', advance='no')
            if (present(calc_pgen)) write(i_unit_, '("   calc_pgen |")', advance='no')
            write(i_unit_, *) ! linebreak

            successful = .true.
            do i = 1, n_excits
                logical :: pgens_differ
                pgen_diagnostic = contrib_list(i) / real(n_iter, dp)
                ic = findbitexcitlevel(ilut, det_list(:, i))
                call decode_bit_det(nJ, det_list(:, i))
                call get_excitation(nI, nJ, ic, exc, tParity)

                if (present(calc_pgen)) then
                    pgens_differ = .not. (calc_pgen(nI, ilut, exc, ic, store%ClassCountOcc, store%ClassCountUnocc) .isclose. pgen_list(i))
                    pgens_differ = .false.
                end if

                if (problem_filter_(nI, exc, ic, pgen_diagnostic) .or. pgens_differ) then
                    successful = .false.
                    call write_det(i_unit_, nJ, .false.)
                    write(i_unit_, '("|"F10.5"|"I2"|"F10.5"|"F15.10"|")', advance='no') &
                        pgen_diagnostic, ic, get_helement(nI, nJ), pgen_list(i)
                    if (present(calc_pgen)) then
                        write(i_unit_, '(F15.10"|")', advance='no') &
                            calc_pgen(nI, det_list(:, i), exc, ic, store%ClassCountOcc, store%ClassCountUnocc)
                    end if
                    write(i_unit_, *)
                end if
            end block
            end do
        end block

        call clean_excit_gen_store(store)


        logical function default_predicate(nI, exc, ic, pgen_diagnostic)
            integer, intent(in) :: nI(nEl), exc(2, maxExcit), ic
            real(dp), intent(in) :: pgen_diagnostic
            default_predicate = &
                (abs(1._dp - pgen_diagnostic) >= 0.05_dp) &
                .and. .not. near_zero(dyn_sltcnd_excit_old(nI, ic, exc, .true.))
        end function

    end subroutine run_excit_gen_tester_function