subroutine test_excitation_generator(sampleSize, pTot, pNull, numEx, nFound, t_calc_pgen, start_nI)
! Test an excitation generator by creating a large number of excitations and
! compare the generated excits with a precomputed list of all excitations
! We thus make sure that
! a) all possible excitations are generated with some weight
! b) no invalid excitations are obtained
integer, intent(in) :: sampleSize
real(dp), intent(out) :: pTot, pNull
integer, intent(out) :: numEx, nFound
logical, intent(in) :: t_calc_pgen
integer, intent(in), optional :: start_nI(nEl)
integer :: nI(nel), nJ(nel)
integer :: i, ex(2, maxExcit), exflag
integer(n_int) :: ilut(0:NIfTot), ilutJ(0:NIfTot)
real(dp) :: pgen
logical :: tPar, tAllExFound, tFound
integer :: j, nSingles, nDoubles
integer(n_int), allocatable :: allEx(:, :)
real(dp) :: pgenArr(lenof_sign)
real(dp) :: matel, matelN, pgenCalc
logical :: exDoneDouble(0:nBasis, 0:nBasis, 0:nBasis, 0:nBasis)
logical :: exDoneSingle(0:nBasis, 0:nBasis)
integer :: ic, part, nullExcits
integer :: ClassCountOcc(scratchSize), ClassCountUnocc(scratchSize)
integer(int64) :: start, finish, rate
character(*), parameter :: t_r = "test_excitation_generator"
HElement_t(dp) :: HEl
exDoneDouble = .false.
exDoneSingle = .false.
call system_clock(count_rate=rate)
! some starting det - do NOT use the reference for the pcpp test, that would
! defeat the purpose
if (present(start_nI)) then
nI = start_nI
else
do i = 1, nel
if (2 * i < nBasis) then
nI(i) = 2 * i - mod(i, 2)
else
nI(i) = i
end if
end do
call sort(nI)
end if
call EncodeBitDet(nI, ilut)
exflag = 3
ex = 0
! create a list of all singles and doubles for reference
call CountExcitations3(nI, exflag, nSingles, nDoubles)
allocate(allEx(0:(NIfTot + 1), nSingles + nDoubles), source=0_n_int)
numEx = 0
tAllExFound = .false.
do
call GenExcitations3(nI, ilut, nJ, exflag, ex(:,1:2), tPar, tAllExFound, .false.)
if (tAllExFound) exit
call encodeBitDet(nJ, ilutJ)
numEx = numEx + 1
allEx(0:nifd, numEx) = ilutJ(0:nifd)
end do
call sort(nI)
write(stdout, *) "In total", numEx, "excits, (", nSingles, nDoubles, ")"
write(stdout, *) "Exciting from", nI
call EncodeBitDet(nI, ilut)
! set the biases for excitation generation
pParallel = 0.5_dp
pSingles = 0.1_dp
pDoubles = 0.9_dp
pNull = 0.0_dp
nullExcits = 0
call system_clock(start)
do i = 1, sampleSize
fcimc_excit_gen_store%tFilled = .false.
call generate_excitation(nI, ilut, nJ, ilutJ, exFlag, ic, ex, tPar, pgen, HEl, fcimc_excit_gen_store, part)
! lookup the excitation
tFound = .false.
do j = 1, numEx
if (DetBitEQ(ilutJ, allEx(:, j))) then
pgenArr = pgen
call encode_sign(allEx(:, j), pgenArr)
! increase its counter by 1
allEx(NIfTot + 1, j) = allEx(NIfTot + 1, j) + 1
tFound = .true.
exit
end if
end do
! if it was not found, and is not marked as invalid, we have a problem: this is not
! an excitaion
if (.not. tFound .and. .not. nJ(1) == 0) then
call decode_bit_det(nJ, ilutJ)
write(stdout, *) "Created excitation", nJ
call stop_all(t_r, "Error: Invalid excitation")
end if
! check if the generated excitation is invalid, if it is, mark this specific constellation
! so we do not double-count when calculating pNull
if (nJ(1) == 0) then
nullExcits = nullExcits + 1
if (ic == 2) then
if (.not. exDoneDouble(ex(1, 1), ex(1, 2), ex(2, 1), ex(2, 2))) then
exDoneDouble(ex(1, 1), ex(1, 2), ex(2, 1), ex(2, 2)) = .true.
pNull = pNull + pgen
end if
else if (ic == 1) then
if (.not. exDoneSingle(ex(1, 1), ex(1, 2))) then
exDoneSingle(ex(1, 1), ex(1, 2)) = .true.
pNull = pNull + pgen
end if
end if
end if
end do
call system_clock(finish)
! check that all excits have been generated and all pGens are right
! probability normalization
pTot = pNull
matelN = 0.0
do i = 1, numEx
call decode_bit_det(nJ, allEx(:, i))
matelN = matelN + abs(get_helement(nI, nJ))
end do
nFound = 0
! class counts might be required for comparing the pgen
call construct_class_counts(nI, classCountOcc, classCountUnocc)
do i = 1, numEx
call extract_sign(allEx(:, i), pgenArr)
call decode_bit_det(nJ, allEx(:, i))
matel = get_helement(nI, nJ)
if (pgenArr(1) > eps) then
nFound = nFound + 1
write(stdout, *) i, pgenArr(1), real(allEx(NIfTot + 1, i)) / real(sampleSize), &
abs(matel) / (pgenArr(1) * matelN)
! compare the stored pgen to the directly computed one
if (t_calc_pgen) then
if (i > nSingles) then
ic = 2
else
ic = 1
end if
ex(1, 1) = 2
call getBitExcitation(ilut, allEx(:, i), ex, tPar)
pgenCalc = calc_pgen(nI, ilut, ex, ic, ClassCountOcc, ClassCountUnocc)
if (abs(pgenArr(1) - pgenCalc) > eps) then
write(stdout, *) "Stored: ", pgenArr(1), "calculated:", pgenCalc
write(stdout, *) "For excit", nJ
call stop_all(t_r, "Incorrect pgen")
end if
end if
else
! excitations with zero matrix element are not required to be found
if (abs(matel) < eps) then
nFound = nFound + 1
else if (i < nSingles) then
write(stdout, *) "Unfound single excitation", nJ
else
write(stdout, *) "Unfound double excitation", nJ, matel
end if
end if
pTot = pTot + pgenArr(1)
end do
write(stdout, *) "Total prob. ", pTot
write(stdout, *) "pNull ", pNull
write(stdout, *) "Null ratio", nullExcits / real(sampleSize)
write(stdout, *) "In total", numEx, "excitations"
write(stdout, *) "With", nSingles, "single excitation"
write(stdout, *) "Found", nFound, "excitations"
write(stdout, *) 'Elapsed Time in seconds:', dble(finish - start) / dble(rate)
write(stdout, *) 'Elapsed Time in micro seconds per excitation:', dble(finish - start) * 1e6_dp / dble(sampleSize* rate)
end subroutine test_excitation_generator