subroutine calcRaisingSemiStartStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
posSwitches, t, probWeight)
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) :: negSwitches(nSpatOrbs), posSwitches(nSpatOrbs)
integer(n_int), intent(inout) :: t(0:nifguga)
real(dp), intent(inout) :: probWeight
character(*), parameter :: this_routine = "calcRaisingSemiStartStochastic"
integer :: se, deltaB
real(dp) :: tempWeight_0, tempWeight_1, minusWeight, &
plusWeight, zeroWeight, bVal
ASSERT(.not. isThree(ilut, excitInfo%secondStart))
! i can be sure that there is no 3 or 0 at the fullEnd, or otherwise
! this would be single-excitation like.
se = excitInfo%secondStart
bVal = csf_i%B_real(se)
deltaB = getDeltaB(t)
! do non-choosing possibs first
! why does this cause a segfault on compilation with gfortran??
! do some debugging:
! fix for gfortran compilation for some reasono
! i can probably fix it when i finally get to this point in
! test running
select case (csf_i%stepvector(se))
case (1)
! 1 -> 3
set_orb(t, 2 * se)
call setDeltaB(deltaB + 1, t)
call getDoubleMatrixElement(3, 1, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, &
tempWeight_0, tempWeight_1)
case (2)
! 2 -> 3
set_orb(t, 2 * se - 1)
call setDeltaB(deltaB - 1, t)
call getDoubleMatrixElement(3, 2, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, &
tempWeight_0, tempWeight_1)
case (0)
! 0:
! for arriving -1 branch branching is always possible
if (deltaB == -1) then
! here the choice is between 0 and -2 branch
minusWeight = weights%proc%minus(negSwitches(se), bVal, weights%dat)
zeroWeight = weights%proc%zero(negSwitches(se), posSwitches(se), &
bVal, weights%dat)
if (near_zero(minusWeight + zeroWeight)) then
probWeight = 0.0_dp
t = 0
return
end if
! cant directly cant assign it to probWeight
zeroWeight = calcStartProb(zeroWeight, minusWeight)
if (genrand_real2_dSFMT() < zeroWeight) then
! go to 0 branch
! 0 -> 2
set_orb(t, 2 * se)
call setDeltaB(0, t)
call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, &
tempWeight_0, tempWeight_1)
probWeight = probWeight * zeroWeight
else
! to to -2 branch
! 0 -> 1
set_orb(t, 2 * se - 1)
call setDeltaB(-2, t)
call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, &
x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
probWeight = probWeight * (1.0_dp - zeroWeight)
end if
else
! +1 branch arriving -> have to check b values
! UPDATE: include b value check into probWeight calculation
! todo
if (csf_i%B_int(se) < 2) then
! only 0 branch possible
! todo: in this forced cases due to the b value, have to
! think about, how that influences the probWeight...
! 0 -> 1
set_orb(t, 2 * se - 1)
call setDeltaB(0, t)
call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, &
tempWeight_0, tempWeight_1)
else
! need probs
plusWeight = weights%proc%plus(posSwitches(se), &
bVal, weights%dat)
zeroWeight = weights%proc%zero(negSwitches(se), posSwitches(se), &
bVal, weights%dat)
if (near_zero(plusWeight + zeroWeight)) then
probWeight = 0.0_dp
t = 0
return
end if
! cant directly cant assign it to probWeight
zeroWeight = calcStartProb(zeroWeight, plusWeight)
if (genrand_real2_dSFMT() < zeroWeight) then
! do 0 branch
! 0 -> 1
set_orb(t, 2 * se - 1)
call setDeltaB(0, t)
call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, &
tempWeight_0, tempWeight_1)
probWeight = probWeight * zeroWeight
else
! do +2 branch
! 0 -> 2
set_orb(t, 2 * se)
call setDeltaB(2, t)
call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order, &
x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
probWeight = probWeight * (1.0_dp - zeroWeight)
end if
end if
end if
end select
call encode_matrix_element(t, extract_matrix_element(t, 1) * tempWeight_1, 2)
call update_matrix_element(t, tempWeight_0, 1)
if (near_zero(tempWeight_0) .and. near_zero(tempWeight_1)) then
probWeight = 0.0_dp
t = 0
end if
end subroutine calcRaisingSemiStartStochastic