subroutine pickOrbitals_FastFast(this, nI, ilut, csf_i, excitInfo, pgen)
class(PropVec_PCHB_FastFast_t), intent(in) :: this
integer, intent(in) :: nI(nel)
integer(n_int), intent(in) :: ilut(0:GugaBits%len_tot)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(out) :: excitInfo
real(dp), intent(out) :: pgen
integer :: src(2), a, b, orbs(2), ij, ab, i_sg
real(dp) :: pgen_elec, pgen_orbs
routine_name("pickOrbitals_FastFast")
if (this%use_lookup) then
i_sg = this%indexer%lookup_prop_vec_idx(csf_I)
else
i_sg = this%indexer%idx_nI(nI)
end if
call this%particle_selector%draw(nI, ilut, csf_i, i_sg, src, pgen_elec)
ASSERT(src(1) <= src(2))
if (src(1) == 0) then
excitInfo%valid = .false.
pgen = 0.0_dp
return
end if
! use the sampler for this electron pair -> order of src electrons
! does not matter
ij = fuse_symm_idx(src(1), src(2))
! get a pair of orbitals using the precomputed weights
call this%alias_sampler%sample(ij, i_sg, ab, pgen_orbs)
! can happen for PropVec and is ok
if (ab == 0) then
excitInfo%valid = .false.
pgen = 0.0_dp
return
! unfortunately, there is a super-rare case when, due to floating point error,
! an excitation with pGen=0 is created. Those are invalid, too
else if (near_zero(pgen_orbs)) then
excitInfo%valid = .false.
pgen = 0.0_dp
! Yes, print. Those events are signficant enough to be always noted in the output
print *, "WARNING: Generated excitation with probability of 0"
return
end if
! split the index ab (using a table containing mapping ab -> (a,b))
orbs = this%tgtOrbs(:, ab)
a = orbs(1); b = orbs(2)
! check if the picked orbs are a valid choice - if they are the same, match one
! occupied orbital or are zero (maybe because there are no allowed picks for
! the given source) abort
! these are the 'easy' checks for GUGA.. more checks need to be done
! to see if it is actually a valid combination..
if (any(orbs == 0) &
.or. (csf_i%stepvector(a) == 3) &
.or. (csf_i%stepvector(b) == 3) &
.or. (a == b .and. csf_i%stepvector(b) /= 0)) then
excitInfo%valid = .false.
return
end if
call extract_excit_info(this%get_info_entry(ij, ab), excitInfo)
excitInfo%i_sg_start = i_sg
pGen = pgen_elec * pgen_orbs
ASSERT(pgen .isclose. this%calc_pgen(nI, ilut, csf_i, excitInfo))
end subroutine pickOrbitals_FastFast