subroutine setup_pchb_sampler_conditional(this)
class(GugaAliasSampler_t), intent(inout) :: this
integer :: i, j, ij, a, b, ab
integer(int64) :: ijMax, abMax
integer(int64), allocatable :: excit_info(:)
real(dp), allocatable :: w(:)
ijMax = fuseIndex(nSpatOrbs, nSpatOrbs)
abMax = fuseIndex(nSpatOrbs, nSpatOrbs)
! just to be save also code up a setup function with the
! necessary if statements, so i can check if my optimized
! sampler is correct (and if it is even faster..)
! weights per pair
allocate(w(abMax), source = 0.0_dp)
! excitation information per pair:
allocate(excit_info(abMax))
! allocate: all samplers have the same size
call this%alias_sampler%shared_alloc(int(ijMax), int(abMax), "GUGA-pchb")
! todo: do the same for the excit_info array!
call this%new_info_table(ijMax, abMax)
! the encode_excit_info function encodes as: E_{ai}E_{bj)
! but for some index combinations we have to change the input so
! actually E_{aj}E_{bi} is encoded! convention to use only
! certain type of excit-types
do i = 1, nSpatOrbs
do j = i, nSpatOrbs
w = 0.0_dp
ij = fuseIndex(i,j)
excit_info = 0_int64
do a = 1, nSpatOrbs
do b = a, nSpatOrbs
! shoud i point group symmetry restrctions here?
! this would avoid unnecessary if statements..
if (RandExcitSymLabelProd(&
SpinOrbSymLabel(2*i), SpinOrbSymLabel(2*j)) &
/= RandExcitSymLabelProd(&
SpinOrbSymLabel(2*a), SpinOrbSymLabel(2*b))) then
cycle
end if
ab = fuseIndex(a,b)
call get_weight_and_info(i, j, a, b, w(ab), excit_info(ab))
end do
end do
call this%alias_sampler%setup_entry(ij, root, w)
call this%set_info_entry(ij, excit_info)
end do
end do
end subroutine setup_pchb_sampler_conditional