Given the initial determinant (both as nI and ilut), create a random double excitation using the hamiltonian matrix elements as weights
GAS_PCHB_DoublesSpinorbFullyWeightedExcGenerator_t
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(GAS_PCHB_DoublesSpinorbFullyWeightedExcGenerator_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)
!! Given the initial determinant (both as nI and ilut), create a random
!! double excitation using the hamiltonian matrix elements as weights
class(GAS_PCHB_DoublesSpinorbFullyWeightedExcGenerator_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), tgt(2), IJ
real(dp) :: p_IJ, p_first(2), p_second(2), renorm_first, renorm_second(2)
integer :: unoccupied(nBasis - nEl)
integer(n_int) :: ilut_unoccupied(0 : nIfD)
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, p_IJ)
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. (p_IJ .isclose. this%particle_selector%get_pgen(nI, i_sg, src(1), src(2)))) then
write(stderr, *) ""
write(stderr, *) "Assertion p_IJ .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_spinorb_fullyweighted.fpp:179"
call stop_all (this_routine, "Assert fail: p_IJ .isclose. this%particle_selector%get_pgen(nI, i_sg, src(1), src(2))")
end if
end block
#endif
call this%get_unoccupied(ilutI(0 : nIfD), ilut_unoccupied, unoccupied)
IJ = fuseIndex(src(1), src(2))
renorm_first = 1._dp - sum(this%A_sampler%get_prob(IJ, i_sg, nI))
if (do_direct_calculation(renorm_first)) then
renorm_first = sum(this%A_sampler%get_prob(IJ, i_sg, unoccupied))
end if
block
integer :: dummy
call this%A_sampler%constrained_sample(&
IJ, i_sg, unoccupied, ilut_unoccupied, renorm_first, dummy, tgt(1), p_first(1))
end block
if (tgt(1) == 0) then
call invalidate()
return
end if
renorm_second(1) = 1._dp - sum(this%B_sampler%get_prob(tgt(1), IJ, i_sg, nI))
if (do_direct_calculation(renorm_second(1))) then
renorm_second(1) = sum(this%B_sampler%get_prob(tgt(1), IJ, i_sg, unoccupied))
end if
if (near_zero(renorm_second(1))) then
call invalidate()
return
end if
block
integer :: dummy
call this%B_sampler%constrained_sample(&
tgt(1), IJ, i_sg, unoccupied, ilut_unoccupied, &
renorm_second(1), dummy, tgt(2), p_second(1))
end block
if (tgt(2) == 0) then
call invalidate()
return
end if
! We could have picked them the other way round.
! Account for that.
! The renormalization for the first hole is the same
p_first(2) = this%A_sampler%constrained_getProb(&
IJ, i_sg, unoccupied, renorm_first, tgt(2))
renorm_second(2) = 1._dp - sum(this%B_sampler%get_prob(tgt(2), IJ, i_sg, nI))
if (do_direct_calculation(renorm_second(2))) then
renorm_second(2) = sum(this%B_sampler%get_prob(tgt(2), IJ, i_sg, unoccupied))
end if
p_second(2) = this%B_sampler%constrained_getProb(&
tgt(2), IJ, i_sg, unoccupied, renorm_second(2), tgt(1))
pGen = p_IJ * sum(p_first * p_second)
call make_double(nI, nJ, elecs(1), elecs(2), tgt(1), tgt(2), ex, tParity)
ilutJ = exciteIlut(ilutI, src, tgt)
block
integer :: ClassCount2(ScratchSize), ClassCountUnocc2(ScratchSize)
#ifdef DEBUG_
block
use util_mod, only: stop_all
use constants, only: stderr
if (.not. (isclose(pgen, this%get_pgen(nI, ilutI, ex, ic, ClassCount2, ClassCountUnocc2), atol=1e-10_dp))) then
write(stderr, *) ""
write(stderr, *) "Assertion isclose(pgen, this%get_pgen(nI, ilutI, ex, ic, ClassCount2, ClassCountUnocc2),&
& atol=1e-10_dp)"
write(stderr, *) "failed in /scratch/jenkins/jobs/existing_branch_doc/workspace/build_config/gfortran-doc/src/gasci_pch&
&b_doubles_spinorb_fullyweighted.fpp:236"
write(stderr, *) "Values leading to this error:"
write(stderr, *) "pgen = ", pgen
write(stderr, *) "this%get_pgen(nI, ilutI, ex, ic, ClassCount2, ClassCountUnocc2) = ", this%get_pgen(nI, ilutI, ex,&
& ic, ClassCount2, ClassCountUnocc2)
write(stderr, *) "abs(pgen - this%get_pgen(nI, ilutI, ex, ic, ClassCount2, ClassCountUnocc2)) = ", abs(pgen -&
& this%get_pgen(nI, ilutI, ex, ic, ClassCount2, ClassCountUnocc2))
call stop_all (this_routine, "Assert fail: isclose(pgen, this%get_pgen(nI, ilutI, ex, ic, ClassCount2,&
& ClassCountUnocc2), atol=1e-10_dp)")
end if
end block
#endif
end block
contains
subroutine invalidate()
nJ = 0
ilutJ = 0_n_int
ex(1, 1 : 2) = src
ex(2, 1 : 2) = tgt
end subroutine invalidate
end subroutine GAS_doubles_PCHB_gen_exc