subroutine calcLoweringSemiStopStochastic(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 = "calcLoweringSemiStopStochastic"
integer :: semi, deltaB
real(dp) :: tempWeight_0, tempWeight_1, minusWeight, &
plusWeight, bVal
ASSERT(.not. isThree(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
return
end if
end if
ASSERT(getDeltaB(t) /= 2)
! do the 1 -> 3 switch
set_orb(t, 2 * semi)
call setDeltaB(deltaB + 1, t)
call getDoubleMatrixElement(3, 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
return
end if
end if
ASSERT(getDeltaB(t) /= -2)
! do 2 -> 3
set_orb(t, 2 * semi - 1)
call setDeltaB(deltaB - 1, t)
call getDoubleMatrixElement(3, 2, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
tempWeight_0, tempWeight_1)
case (0)
! its a 0 -> have to check b if a 0 branch arrives
select case (deltaB)
case (-2)
! do 0 -> 2
set_orb(t, 2 * semi)
call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
call setDeltaB(-1, t)
case (2)
! do 0 -> 1
set_orb(t, 2 * semi - 1)
call setDeltaB(+1, t)
call getDoubleMatrixElement(1, 0, 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 0 -> 1 possble
set_orb(t, 2 * semi - 1)
call setDeltaB(-1, t)
call getDoubleMatrixElement(1, 0, 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 0 -> 1
set_orb(t, 2 * semi - 1)
call setDeltaB(-1, t)
call getDoubleMatrixElement(1, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
tempWeight_0, tempWeight_1)
probWeight = probWeight * minusWeight
else
! do +1 branch
! set 0 -> 2
set_orb(t, 2 * semi)
call setDeltaB(1, t)
call getDoubleMatrixElement(2, 0, deltaB, excitInfo%gen1, &
excitInfo%gen2, bVal, excitInfo%order1, &
tempWeight_0, tempWeight_1)
probWeight = probWeight * (1.0_dp - minusWeight)
end if
end if
end select
end select
! for mixed fullstart check if no switch happened in the double
! overlap region, indicated by a non-zero-x0 matrix element
if (excitInfo%typ == excit_type%fullStart_L_to_R) then
if (.not. near_zero(extract_matrix_element(t, 1) * tempWeight_0)) then
probWeight = 0.0_dp
t = 0
return
end if
end if
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
return
end if
call update_matrix_element(t, tempWeight_0, 1)
call update_matrix_element(t, tempWeight_1, 2)
end subroutine calcLoweringSemiStopStochastic