subroutine init_pcpp_singles_excitgen()
! create the probability distributions for drawing single excitations
! The normalization 1/n_{jb} used by Neufeld and Thom seems to be just a constant, we
! omit it for now
implicit none
call setup_elecs_sampler()
call setup_holes_sampler()
contains
!------------------------------------------------------------------------------------------!
subroutine setup_elecs_sampler()
! the probability distribution for selection of electrons
! creates a sampler for electron indices within the reference determinant
! these later have to be transferred to the current determinant
! even though strictly speaking, these are sums of the hole probabilities,
! expressing them in terms of the latter would be an unwanted dependency,
! since this relation is merely a matter of choice of the algorithm and should not
! be reflected in the design of the implementation
implicit none
real(dp) :: w(nel)
integer :: i, a
integer :: iEl
do iEl = 1, nel
w(iEl) = 0
i = refDet(iEl)
do a = 1, nBasis
w(iEl) = w(iEl) + acc_doub_matel(i, a)
end do
end do
! load the probabilites for electron selection into the alias table
call single_elec_sampler%setup(root, w)
end subroutine setup_elecs_sampler
!------------------------------------------------------------------------------------------!
subroutine setup_holes_sampler()
! the probability distributions for selection of holes, given the electron orbital
implicit none
real(dp) :: w(nBasis)
integer :: i, a
integer :: aerr
! there is one table for given source orbital
allocate(single_hole_sampler(nBasis), stat=aerr)
! each table has probabilities for all given virtual orbitals a (some of them might
! be 0, this has no additional cost)
do i = 1, nBasis
w = 0.0_dp
do a = 1, nBasis
! we never want to sample the source orbital
if (i /= a) &
! store the accumulated matrix elements (= un-normalized probability) with
! the corresponding symmetry (if spins of a/i are different, w is 0)
w(a) = acc_doub_matel(i, a)
end do
call single_hole_sampler(i)%setup(root, w)
end do
end subroutine setup_holes_sampler
!------------------------------------------------------------------------------------------!
function acc_doub_matel(src, tgt) result(prob)
! Accumulate all single excitation matrix elements connecting
! D_(j,src)^(b,tgt) for all (j,b), where D is the reference
! Input: src - orbital we want to excite from (any orbital is allowed)
! tgt - orbital we want to excite to (any orbital is allowed)
! Output: prob - Accumulated matrix elements of singles
! If src/tgt have different spin, the result is 0, if they have
! different symmetry, it is not necessarily
!
! IMPORTANT: The matrix elements are calculated as
! <D_j^b|H|D_(j,src)^(b,tgt)> := h_src^tgt + sum_(k occ in D_j^b) <tgt k|src k> - <tgt k|k src>
! That is, the left side is to be understood symbolic, there is no actual excitation
! from src to tgt
! symmetry has to be preserved
implicit none
integer, intent(in) :: src, tgt
real(dp) :: prob
integer :: b, j
integer :: nI(nel)
type(Excite_1_t) :: ex
logical :: tPar
prob = 0
if (symAllowed(src, tgt)) then
do b = 1, nBasis
! loop over all non-occupied orbitals
if (.not. any(b == refDet(:))) then
do j = 1, nel
! get the excited determinant D_j^b used for the matrix element
if (symAllowed(refDet(j), b)) then
call make_single(refDet(:), nI, j, b, ex%val, tPar)
! this is a symbolic excitation, we do NOT require src to be occupied
! we just use the formula for excitation matrix elements
prob = prob + abs(sltcnd_excit(nI, Excite_1_t(src, tgt), tPar, assert_occupation=.false.))
end if
end do
end if
end do
end if
end function acc_doub_matel
!------------------------------------------------------------------------------------------!
end subroutine init_pcpp_singles_excitgen