subroutine doubleUpdateStochastic(ilut, csf_i, s, excitInfo, weights, negSwitches, &
posSwitches, t, probWeight)
integer(n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
integer, intent(in) :: s
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 = "doubleUpdateStochastic"
real(dp) :: minusWeight, plusWeight, zeroWeight
integer :: gen1, gen2, deltaB
real(dp) :: bVal, tempWeight_0, tempWeight_1, order
ASSERT(isProperCSF_ilut(ilut))
ASSERT(s > 0 .and. s <= nSpatOrbs)
if (csf_i%Occ_int(s) /= 1) then
! no change in stepvector or matrix element in this case
return
end if
! and for more readibility extract certain values:
gen1 = excitInfo%gen1
gen2 = excitInfo%gen2
bVal = csf_i%B_real(s)
! stupid! only need at order at semistarts and semistops and not for
! the overlap region
order = 1.0_dp
deltaB = getDeltaB(t)
! new idea: make the combination stepvalue + deltaB !
! this give me 6 distinct integer quantities which i can choose
! from in a select case statement!
select case (csf_i%stepvector(s) + deltaB)
! depending on the deltaB value different possibs
case (3)
! d=1 + b=2 = 3
! only staying
call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, bVal, &
order, x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
case (-1)
! d=1 + b=-2 = -1
! both 0 and -2 branches are possible
minusWeight = weights%proc%minus(negSwitches(s), &
bVal, weights%dat)
zeroWeight = weights%proc%zero(negSwitches(s), &
posSwitches(s), bVal, weights%dat)
if (near_zero(minusWeight + zeroWeight)) then
probWeight = 0.0_dp
t = 0
return
end if
minusWeight = calcStayingProb(minusWeight, zeroWeight, bVal)
if (genrand_real2_dSFMT() < minusWeight) then
! stay on -2 branch
call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
probWeight = probWeight * minusWeight
else
! switch to 0 branch
! 1 -> 2
clr_orb(t, 2 * s - 1)
set_orb(t, 2 * s)
call setDeltaB(0, t)
call getDoubleMatrixElement(2, 1, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
probWeight = probWeight * (1.0_dp - minusWeight)
end if
case (1)
! d=1 + b=0 = 1
! 0 branch arrives -> have to check b value
! hm... is this a bug?? if d = 1, b cant be 0, its atleast
! 1.. and then a switch to +2 cant happen..
! but that should have been dealt with the weights below
! probably.. so thats why it didnt matter probably..
if (csf_i%B_int(s) == 1) then
! only staying branch
call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, bVal, &
order, tempWeight_0, tempWeight_1)
else
! both 0 and +2 branch are possible
plusWeight = weights%proc%plus(posSwitches(s), &
bVal, weights%dat)
zeroWeight = weights%proc%zero(negSwitches(s), &
posSwitches(s), bVal, weights%dat)
if (near_zero(plusWeight + zeroWeight)) then
probWeight = 0.0_dp
t = 0
end if
zeroWeight = calcStayingProb(zeroWeight, plusWeight, bVal)
if (genrand_real2_dSFMT() < zeroWeight) then
! stay on 0 branch
call getDoubleMatrixElement(1, 1, deltaB, gen1, gen2, &
bVal, order, tempWeight_0, tempWeight_1)
probWeight = probWeight * zeroWeight
else
! switch to +2 branch
! 1 -> 2
clr_orb(t, 2 * s - 1)
set_orb(t, 2 * s)
call setDeltaB(2, t)
call getDoubleMatrixElement(2, 1, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
probWeight = probWeight * (1.0_dp - zeroWeight)
end if
end if
case (0)
! d=2 + b=-2 : 0
! only staying
call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, bVal, &
order, x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
case (2)
! d=2 + b=0 : 2
! always -2 and 0 branching possible
minusWeight = weights%proc%minus(negSwitches(s), &
bVal, weights%dat)
zeroWeight = weights%proc%zero(negSwitches(s), &
posSwitches(s), bVal, weights%dat)
! here should be the only case where b*w_0 + w_- = 0 ->
! have to avoid divison by 0 then
! but in this case only the stayin on 0 branch is valid
! since it should be always > 0 and the -branch has 0 weight
if (near_zero(bVal * zeroWeight + minusWeight)) then
zeroWeight = 1.0_dp
else
zeroWeight = calcStayingProb(zeroWeight, minusWeight, bVal)
end if
if (near_zero(zeroWeight + minusWeight)) then
probWeight = 0.0_dp
t = 0
return
end if
if (genrand_real2_dSFMT() < zeroWeight) then
! stay on 0 branch
call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, bVal, &
order, tempWeight_0, tempWeight_1)
probWeight = probWeight * zeroWeight
else
! switch to -2 branch
! 2 -> 1
set_orb(t, 2 * s - 1)
clr_orb(t, 2 * s)
call setDeltaB(-2, t)
call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
probWeight = probWeight * (1.0_dp - zeroWeight)
end if
case (4)
! d=2 + b=2 : 4
! have to check b value if branching is possible
if (csf_i%B_int(s) < 2) then
! only switch possible
! 2 -> 1
set_orb(t, 2 * s - 1)
clr_orb(t, 2 * s)
call setDeltaB(0, t)
call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
else
! both 0 and +2 branches possible
plusWeight = weights%proc%plus(posSwitches(s), &
bVal, weights%dat)
zeroWeight = weights%proc%zero(negSwitches(s), &
posSwitches(s), bVal, weights%dat)
if (near_zero(plusWeight + zeroWeight)) then
probWeight = 0.0_dp
t = 0
return
end if
plusWeight = calcStayingProb(plusWeight, zeroWeight, bVal)
if (genrand_real2_dSFMT() < plusWeight) then
! stay on +2 branch
call getDoubleMatrixElement(2, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
probWeight = probWeight * plusWeight
else
! switch to 0 branch
! 2 -> 1
set_orb(t, 2 * s - 1)
clr_orb(t, 2 * s)
call setDeltaB(0, t)
call getDoubleMatrixElement(1, 2, deltaB, gen1, gen2, &
bVal, order, x1_element=tempWeight_1)
tempWeight_0 = 0.0_dp
probWeight = probWeight * (1.0_dp - plusWeight)
end if
end if
end select
call update_matrix_element(t, tempWeight_0, 1)
call update_matrix_element(t, tempWeight_1, 2)
if (near_zero(tempWeight_0) .and. near_zero(tempWeight_1)) then
probWeight = 0.0_dp
t = 0
return
end if
if (t_trunc_guga_pgen .or. &
(t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
if (probWeight < trunc_guga_pgen) then
probWeight = 0.0_dp
t = 0_n_int
end if
end if
if (t_trunc_guga_matel) then
if (abs(extract_matrix_element(t, 1)) < trunc_guga_matel .and. &
abs(extract_matrix_element(t, 2)) < trunc_guga_matel) then
probWeight = 0.0_dp
t = 0_n_int
end if
end if
end subroutine doubleUpdateStochastic