pickOrbitals_FastFast Subroutine

private subroutine pickOrbitals_FastFast(this, nI, ilut, csf_i, excitInfo, pgen)

Type Bound

PropVec_PCHB_FastFast_t

Arguments

Type IntentOptional Attributes Name
class(PropVec_PCHB_FastFast_t), intent(in) :: this
integer, intent(in) :: nI(nel)
integer(kind=n_int), intent(in) :: ilut(0:GugaBits%len_tot)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(out) :: excitInfo
real(kind=dp), intent(out) :: pgen

Source Code

    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