subroutine createStochasticStart_single(ilut, csf_i, excitInfo, weights, posSwitches, &
negSwitches, t, probWeight)
! create a stochastic start for a single generator
! essentially exactly the same as in the full excitation calculation
! except that at a 0/3 start we have to do it stochastically ...
integer(n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(in) :: excitInfo
type(WeightObj_t), intent(in) :: weights
real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
integer(n_int), intent(out) :: t(0:nifguga) ! to use macros, have to use short names..
real(dp), intent(out) :: probWeight
character(*), parameter :: this_routine = "createStochasticStart_single"
real(dp) :: tempWeight, minusWeight, plusWeight, bVal
integer :: st, gen
! and also the possibility of the excitation should be ensured due to
! the correct choosing of excitation indices..
#ifdef DEBUG_
! also assert we are not calling it for a weight gen. accidently
ASSERT(excitInfo%currentGen /= 0)
ASSERT(isProperCSF_ilut(ilut))
! also check if calculated b vector really fits to ilut
ASSERT(all(csf_i%B_real .isclose. calcB_vector_ilut(ilut(0:nifd))))
if (excitInfo%currentGen == gen_type%R) then
ASSERT(.not. isThree(ilut, excitInfo%fullStart))
else if (excitInfo%currentGen == gen_type%L) then
ASSERT(.not. isZero(ilut, excitInfo%fullStart))
end if
! also have checked if atleast on branch way can lead to an excitaiton
#endif
! for more oversight
st = excitInfo%fullStart
gen = excitInfo%currentGen
bVal = csf_i%B_real(st)
t = ilut
select case (csf_i%stepvector(st))
case (1)
! set corresponding orbital to 0 or 3 depending on generator type
if (gen == gen_type%R) then ! raising gen case
! set the alpha orbital also to 1 to make d=1 -> d'=3
set_orb(t, 2 * st)
! does it work like that:
! would have all necessary values and if statements to directly
! calculated matrix element here... maybe do it here after all
! to be more efficient...
tempWeight = getSingleMatrixElement(3, 1, +1, gen, bVal)
else ! lowering gen case
! clear beta orbital to make d=1 -> d'=0
clr_orb(t, 2 * st - 1)
tempWeight = getSingleMatrixElement(0, 1, +1, gen, bVal)
end if
! set deltaB:
! for both cases deltaB = +1, set that as signed weight
! update: use previously unused flags to encode deltaB
call setDeltaB(1, t)
! have to set probabilistic weight to 1, since only 1 possibility
probWeight = 1.0_dp
case (2)
if (gen == gen_type%R) then
set_orb(t, 2 * st - 1)
! matrix elements
tempWeight = getSingleMatrixElement(3, 2, -1, gen, bVal)
else
clr_orb(t, 2 * st)
! matrix elements
tempWeight = getSingleMatrixElement(0, 2, -1, gen, bVal)
end if
call setDeltaB(-1, t)
probWeight = 1.0_dp
case (0, 3)
! if the deltaB = -1 weights is > 0 this one always works
! write weight calculation function! is a overkill here,
! but makes it easier to use afterwards with stochastic excitaions
! use new weight objects here.
minusWeight = weights%proc%minus(negSwitches(st), bVal, weights%dat)
plusWeight = weights%proc%plus(posSwitches(st), bVal, weights%dat)
! if both weights an b value allow both excitations, have to choose
! one stochastically
! also have to consider, that i might not actually can assure
! possible excitations, since im not yet sure if i check for
! possible switches... -> check that todo
! for now assume i checked for that, with checkCompatibility() or
! something similar..
! do not need the previous if structs since the weights
! already care for the facts if some excitation is incompatible..
! I only have to avoid a completely improper excitation...
! if b == 0 and minusWeight == 0 -> no excitation possible todo!
if (near_zero(plusWeight + minusWeight) .or. &
near_zero(bVal + minusWeight)) then
probWeight = 0.0_dp
t = 0
return
end if
! calc starting weight to go for the minusBranch
probWeight = calcStartProb(minusWeight, plusWeight)
if (genrand_real2_dSFMT() < probWeight) then
! do the -1 branch
if (gen == gen_type%R) then
ASSERT(isZero(ilut, st))
! set beta bit
set_orb(t, 2 * st - 1)
! matrix element
tempWeight = getSingleMatrixElement(1, 0, -1, gen, bVal)
else
ASSERT(isThree(ilut, st))
! clear alpha bit
clr_orb(t, 2 * st)
! matrix element
tempWeight = getSingleMatrixElement(1, 3, -1, gen, bVal)
end if
call setDeltaB(-1, t)
else
! do the +1 branch
if (gen == gen_type%R) then
ASSERT(isZero(ilut, st))
! set alpha bit
set_orb(t, 2 * st)
! matrix element
tempWeight = getSingleMatrixElement(2, 0, +1, gen, bVal)
else
ASSERT(isThree(ilut, st))
! cleat beta bit
clr_orb(t, 2 * st - 1)
! matrix element
tempWeight = getSingleMatrixElement(2, 3, +1, gen, bVal)
end if
call setDeltaB(+1, t)
! update the probabilistic weight with the actual use one
probWeight = 1.0_dp - probWeight
end if
end select
call encode_matrix_element(t, tempWeight, 1)
end subroutine createStochasticStart_single