@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 | Intent | Optional | 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 |
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
else
i_unit_ = stdout
end if
if (present(problem_filter)) then
problem_filter_ => problem_filter
else
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
else
! use HF-det as default
nI = [(i, i=1, nel)]
end if
if (present(opt_n_dets)) then
n_dets_target = opt_n_dets
else
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)
else
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)
block
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')
else
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
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
block
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))
else
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)
contains
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