@brief Given the initial determinant (both as nI and ilut), create a random double excitation using the hamiltonian matrix elements as weights
@param[in] nI determinant to excite from @param[in] elec_map map to translate electron picks to orbitals @param[in] ilut determinant to excite from in ilut format @param[out] nJ on return, excited determinant @param[out] excitMat on return, excitation matrix nI -> nJ @param[out] tParity on return, the parity of the excitation nI -> nJ @param[out] pGen on return, the probability of generating the excitation nI -> nJ @param[in] idet Optional index of determinant in the CurrentDets array.
GAS_PCHB_DoublesSpatOrbFastWeightedExcGenerator_t
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(GAS_PCHB_DoublesSpatOrbFastWeightedExcGenerator_t), | intent(inout) | :: | this | |||
| integer, | intent(in) | :: | nI(nel) | |||
| integer(kind=n_int), | intent(in) | :: | ilutI(0:NIfTot) | |||
| integer, | intent(out) | :: | nJ(nel) | |||
| integer(kind=n_int), | intent(out) | :: | ilutJ(0:NifTot) | |||
| integer, | intent(in) | :: | exFlag | |||
| integer, | intent(out) | :: | ic | |||
| integer, | intent(out) | :: | ex(2,maxExcit) | |||
| logical, | intent(out) | :: | tParity | |||
| real(kind=dp), | intent(out) | :: | pGen | |||
| real(kind=dp), | intent(out) | :: | hel | |||
| type(excit_gen_store_type), | intent(inout), | target | :: | store | ||
| integer, | intent(in), | optional | :: | part_type |
subroutine GAS_doubles_PCHB_gen_exc(&
this, nI, ilutI, nJ, ilutJ, exFlag, ic, &
ex, tParity, pGen, hel, store, part_type)
class(GAS_PCHB_DoublesSpatOrbFastWeightedExcGenerator_t), intent(inout) :: this
integer, intent(in) :: nI(nel), exFlag
integer(n_int), intent(in) :: ilutI(0:NIfTot)
integer, intent(out) :: nJ(nel), ic, ex(2, maxExcit)
integer(n_int), intent(out) :: ilutJ(0:NifTot)
real(dp), intent(out) :: pGen
logical, intent(out) :: tParity
HElement_t(dp), intent(out) :: hel
type(excit_gen_store_type), intent(inout), target :: store
integer, intent(in), optional :: part_type
character(*), parameter :: this_routine = 'GAS_doubles_PCHB_gen_exc'
integer :: elecs(2), src(2), ij
integer :: orbs(2), ab
real(dp) :: pGenHoles
logical :: invalid
integer :: spin(2), samplerIndex
integer :: i_sg
#ifdef WARNING_WORKAROUND_
associate(exFlag => exFlag); end associate
associate(part_type => part_type); end associate
#endif
ic = 2
#ifdef WARNING_WORKAROUND_
hel = h_cast(0.0_dp)
#endif
if (this%use_lookup) then
i_sg = this%indexer%lookup_supergroup_idx(store%idx_curr_dets, nI)
else
i_sg = this%indexer%idx_nI(nI)
end if
! first, pick two random elecs
call this%particle_selector%draw(nI, ilutI, i_sg, elecs, src, pGen)
if (src(1) == 0) then
call invalidate()
return
end if
#ifdef DEBUG_
block
use util_mod, only: stop_all
use constants, only: stderr
if (.not. (pGen .isclose. this%particle_selector%get_pgen(nI, i_sg, src(1), src(2)))) then
write(stderr, *) ""
write(stderr, *) "Assertion pGen .isclose. this%particle_selector%get_pgen(nI, i_sg, src(1), src(2))"
write(stderr, *) "failed in /scratch/jenkins/jobs/existing_branch_doc/workspace/build_config/gfortran-doc/src/gasci_pch&
&b_doubles_spatorb_fastweighted.fpp:198"
call stop_all (this_routine, "Assert fail: pGen .isclose. this%particle_selector%get_pgen(nI, i_sg, src(1), src(2))")
end if
end block
#endif
invalid = .false.
! use the sampler for this electron pair -> order of src electrons does not matter
ij = fuseIndex(gtID(src(1)), gtID(src(2)))
! the spin of the electrons: 0 - alpha, 1 - beta
spin = getSpinIndex(src)
! determine type of spin-excitation: same-spin, opp spin w exchange, opp spin w/o exchange
if (spin(1) == spin(2)) then
! same spin
samplerIndex = SAME_SPIN
else
! else, pick exchange with...some ij-spin bias
if (genrand_real2_dSFMT() < this%pExch(ij, i_sg)) then
samplerIndex = OPP_SPIN_EXCH
! adjust pgen
pGen = pGen * this%pExch(ij, i_sg)
! the spins of the target are the opposite of the source spins
call swap(spin(1), spin(2))
else
samplerIndex = OPP_SPIN_NO_EXCH
! adjust pgen
pGen = pGen * (1.0_dp - this%pExch(ij, i_sg))
end if
end if
! get a pair of orbitals using the precomputed weights
call this%pchb_samplers%sample(ij, samplerIndex, i_sg, ab, pGenHoles)
#ifdef DEBUG_
block
use util_mod, only: stop_all
use constants, only: stderr
if (.not. (ab /= 0 .implies. (pGenHoles .isclose. this%pchb_samplers%get_prob(ij, samplerIndex, i_sg, ab)))) then
write(stderr, *) ""
write(stderr, *) "Assertion ab /= 0 .implies. (pGenHoles .isclose. this%pchb_samplers%get_prob(ij, samplerIndex,&
& i_sg, ab))"
write(stderr, *) "failed in /scratch/jenkins/jobs/existing_branch_doc/workspace/build_config/gfortran-doc/src/gasci_pch&
&b_doubles_spatorb_fastweighted.fpp:225"
call stop_all (this_routine, "Assert fail: ab /= 0 .implies. (pGenHoles .isclose. this%pchb_samplers%get_prob(ij,&
& samplerIndex, i_sg, ab))")
end if
end block
#endif
if (ab == 0) then
invalid = .true.
orbs = 0
else
! split the index ab (using a table containing mapping ab -> (a,b))
orbs = this%tgtOrbs(:, ab)
! convert orbs to spin-orbs with the same spin
orbs = 2 * orbs - spin
#ifdef DEBUG_
block
use util_mod, only: stop_all
use constants, only: stderr
if (.not. (all(orbs /= 0) .implies. orbs(1) /= orbs(2))) then
write(stderr, *) ""
write(stderr, *) "Assertion all(orbs /= 0) .implies. orbs(1) /= orbs(2)"
write(stderr, *) "failed in /scratch/jenkins/jobs/existing_branch_doc/workspace/build_config/gfortran-doc/src/gasci_pch&
&b_doubles_spatorb_fastweighted.fpp:236"
call stop_all (this_routine, "Assert fail: all(orbs /= 0) .implies. orbs(1) /= orbs(2)")
end if
end block
#endif
! note: nI is an array, not a scalar, so we need two `any` checks below
invalid = any(orbs == 0) .or. any(orbs(1) == nI) .or. any(orbs(2) == nI)
end if
! unfortunately, there is a super-rare case when, due to floating point error,
! an excitation with pGen=0 is created. Those are invalid, too
if (.not. invalid .and. near_zero(pGenHoles)) then
invalid = .true.
! Yes, print. Those events are signficant enough to be always noted in the output
write(stdout, *) "WARNING: Generated excitation with probability of 0"
end if
if (invalid) then
! if 0 is returned, there are no excitations for the chosen elecs
! -> return nulldet
call invalidate()
else
! else, construct the det from the chosen orbs/elecs
call make_double(nI, nJ, elecs(1), elecs(2), orbs(1), orbs(2), ex, tParity)
ilutJ = exciteIlut(ilutI, src, orbs)
pGen = pGen * pGenHoles
block
integer :: ClassCount2(ScratchSize), ClassCountUnocc2(ScratchSize)
#ifdef DEBUG_
block
use util_mod, only: stop_all
use constants, only: stderr
if (.not. (pgen .isclose. this%get_pgen(nI, ilutI, ex, ic, ClassCount2, ClassCountUnocc2))) then
write(stderr, *) ""
write(stderr, *) "Assertion pgen .isclose. this%get_pgen(nI, ilutI, ex, ic, ClassCount2, ClassCountUnocc2)"
write(stderr, *) "failed in /scratch/jenkins/jobs/existing_branch_doc/workspace/build_config/gfortran-doc/src/gasci_pch&
&b_doubles_spatorb_fastweighted.fpp:264"
call stop_all (this_routine, "Assert fail: pgen .isclose. this%get_pgen(nI, ilutI, ex, ic, ClassCount2,&
& ClassCountUnocc2)")
end if
end block
#endif
end block
end if
contains
subroutine invalidate()
nJ = 0
ilutJ = 0_n_int
ex(1, 1 : 2) = src
ex(2, 1 : 2) = orbs
end subroutine invalidate
end subroutine GAS_doubles_PCHB_gen_exc