subroutine calcRaisingSemiStopStochastic(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 = "calcRaisingSemiStopStochastic"
integer :: semi, deltaB
real(dp) :: tempWeight_0, tempWeight_1, minusWeight, &
plusWeight, bVal
ASSERT(.not. isZero(ilut, excitInfo%firstEnd))
semi = excitInfo%firstEnd
bVal = csf_i%B_real(semi)
! in the stochastic case i am sure that at there is no 3 or 0 at the
! full start... so definetly all deltaB branches can arrive here
! first deal with the forced ones
deltaB = getDeltaB(t)
select case (csf_i%stepvector(semi))
case (1)
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
if (getDeltaB(t) == 2) then
t = 0_n_int
probWeight = 0.0_dp
return
end if
end if
ASSERT(getDeltaB(t) /= 2)
! do the 1 -> 0 switch
clr_orb(t, 2 * semi - 1)
call setDeltaB(deltaB + 1, t)
call getDoubleMatrixElement(0, 1, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
tempWeight_0, tempWeight_1)
case (2)
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
if (getDeltaB(t) == -2) then
t = 0_n_int
probWeight = 0.0_dp
return
end if
end if
ASSERT(getDeltaB(t) /= -2)
! do 2 -> 0
clr_orb(t, 2 * semi)
call setDeltaB(deltaB - 1, t)
call getDoubleMatrixElement(0, 2, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
tempWeight_0, tempWeight_1)
case (3)
! its a 3 -> have to check b if a 0 branch arrives
select case (deltaB)
case (-2)
! do 3 -> 2
clr_orb(t, 2 * semi - 1)
call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
call setDeltaB(-1, t)
case (2)
! do 3 -> 1
clr_orb(t, 2 * semi)
call setDeltaB(+1, t)
call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
case (0)
! deltaB = 0 branch arrives -> check b
if (csf_i%B_int(semi) == 0) then
! only 3 -> 1 possble
clr_orb(t, 2 * semi)
call setDeltaB(-1, t)
call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
tempWeight_0, tempWeight_1)
else
! finally have to check the probs
minusWeight = weights%proc%minus(negSwitches(semi), &
bVal, weights%dat)
plusWeight = weights%proc%plus(posSwitches(semi), &
bVal, weights%dat)
if (near_zero(minusWeight + plusWeight)) then
probWeight = 0.0_dp
t = 0
return
end if
! here i cant directly save it to probWeight ...
minusWeight = calcStartProb(minusWeight, plusWeight)
if (genrand_real2_dSFMT() < minusWeight) then
! do -1 branch:
! set 3 -> 1
clr_orb(t, 2 * semi)
call setDeltaB(-1, t)
call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
tempWeight_0, tempWeight_1)
probWeight = probWeight * minusWeight
else
! do +1 branch
! set 3 -> 2
clr_orb(t, 2 * semi - 1)
call setDeltaB(1, t)
call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
tempWeight_0, tempWeight_1)
probWeight = probWeight * (1 - minusWeight)
end if
end if
end select
end select
! if its a mixed full-start raising into lowering -> check i a
! switch happened in the double overlap region
if (excitInfo%typ == excit_type%fullstart_R_to_L) then
! this is indicated by a non-zero x0-matrix element
if (.not. near_zero(extract_matrix_element(t, 1) * tempWeight_0)) then
probWeight = 0.0_dp
t = 0
return
end if
end if
! do not combine the x0 and x1 elements at this point, because its
! easier to deal with the contributing integrals if we keep them
! seperate until the end!
if (near_zero(extract_matrix_element(t, 1) * tempWeight_0) .and. &
near_zero(extract_matrix_element(t, 2) * tempWeight_1)) then
probWeight = 0.0_dp
t = 0_n_int
return
end if
call update_matrix_element(t, tempWeight_0, 1)
call update_matrix_element(t, tempWeight_1, 2)
end subroutine calcRaisingSemiStopStochastic