subroutine singleStochasticUpdate(ilut, csf_i, s, excitInfo, weights, posSwitches, &
negSwitches, 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) :: 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 = "singleStochasticUpdate"
real(dp) :: tempWeight, bVal, minusWeight, plusWeight
integer :: gen, deltaB, step
! is also very similar to full single update
ASSERT(isProperCSF_ilut(ilut))
ASSERT(s > 0 .and. s <= nSpatOrbs)
! change! have to update the x1 matrix elements here too,
! to keep them seperated to use the individually for the contributing
! integrals
! set standard probWeight
probWeight = 1.0_dp
select case (csf_i%stepvector(s))
case (0)
! do nothin actually.. not even change matrix elements
return
case (3)
! only change matrix element to negative one
call update_matrix_element(t, -1.0_dp, 1)
call update_matrix_element(t, -1.0_dp, 2)
return
end select
! to generally use this function i need to define a current generator...
gen = excitInfo%currentGen
bVal = csf_i%B_real(s)
! only need weights in certain cases..
deltaB = getDeltaB(t)
! is it possible to only check for possible switches and else handle
! it in one go? try!
select case (csf_i%stepvector(s) + deltaB)
case (0)
! d=1 + b=-1 : 0
! no bValue restrictions here, so that should be handlebar
plusWeight = weights%proc%plus(posSwitches(s), bVal, weights%dat)
minusWeight = weights%proc%minus(negSwitches(s), bVal, weights%dat)
! here do a check if not both weights are 0...
if (near_zero(plusWeight + minusWeight)) then
probWeight = 0.0_dp
t = 0
return
end if
! calc. staying probabiliy
probWeight = calcStayingProb(minusWeight, plusWeight, bVal)
if (genrand_real2_dSFMT() < probWeight) then
! stay on -1 branch
tempWeight = getSingleMatrixElement(1, 1, deltaB, gen, bVal)
else
! switch to +1 branch 1 -> 2
set_orb(t, 2 * s)
clr_orb(t, 2 * s - 1)
call setDeltaB(1, t)
tempWeight = getSingleMatrixElement(2, 1, deltaB, gen, bVal)
probWeight = 1.0_dp - probWeight
end if
case (3)
! d=2 + b=1 : 3
! do i need bValue check here?
! probably... to distinguish forced switches
if (csf_i%B_int(s) > 0) then
! both excitations possible
plusWeight = weights%proc%plus(posSwitches(s), bVal, weights%dat)
minusWeight = weights%proc%minus(negSwitches(s), bVal, weights%dat)
! here do a check if not both weights are 0...
if (near_zero(plusWeight + minusWeight)) then
probWeight = 0.0_dp
t = 0
return
end if
probWeight = calcStayingProb(plusWeight, minusWeight, bVal)
if (genrand_real2_dSFMT() < probWeight) then
ASSERT(plusWeight > 0.0_dp)
! stay on +1 branch
tempWeight = getSingleMatrixElement(2, 2, deltaB, gen, bVal)
else
ASSERT(minusWeight > 0.0_dp)
! do the 2 -> 1 switch
clr_orb(t, 2 * s)
set_orb(t, 2 * s - 1)
! change deltaB
call setDeltaB(-1, t)
! and calc matrix elements
tempWeight = getSingleMatrixElement(1, 2, deltaB, gen, bVal)
probWeight = 1.0_dp - probWeight
end if
else
! forced switch
#ifdef DEBUG_
minusWeight = weights%proc%minus(negSwitches(s), bVal, weights%dat)
if (.not. minusWeight > 0.0_dp) then
print *, "+", plusWeight
print *, "-", minusWeight
print *, negSwitches
print *, "overlap:", excitInfo%overlap
call write_det_guga(stdout, ilut)
call write_det_guga(stdout, t)
call print_excitInfo(excitInfo)
end if
ASSERT(minusWeight > 0.0_dp)
#endif
! do 2 -> 1forced switch
clr_orb(t, 2 * s)
set_orb(t, 2 * s - 1)
! change deltaB
call setDeltaB(-1, t)
! and calc matrix elements
tempWeight = getSingleMatrixElement(1, 2, deltaB, gen, bVal)
end if
case (1, 2)
! d=1 + b=1: 2
! d=2 + b=-1: 1
! just staying possibility...
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
plusWeight = weights%proc%plus(posSwitches(s), bVal, weights%dat)
minusWeight = weights%proc%minus(negSwitches(s), bVal, weights%dat)
if ((deltaB == -1 .and. near_zero(minusWeight)) &
.or. (deltaB == 1 .and. near_zero(plusWeight))) then
t = 0_n_int
probWeight = 0.0_dp
return
end if
end if
#ifdef DEBUG_
! for sanity check for weights here too just to be sure to not get
! into non compatible excitations
plusWeight = weights%proc%plus(posSwitches(s), bVal, weights%dat)
minusWeight = weights%proc%minus(negSwitches(s), bVal, weights%dat)
if (deltaB == -1) then
ASSERT(minusWeight > 0.0_dp)
end if
if (deltaB == +1) then
ASSERT(plusWeight > 0.0_dp)
end if
#endif
! can i efficiently code that up?
! deltaB stays the same.., stepvector stays the same.. essentiall
! only matrix element changes.. -> need deltaB, generators and stepvalue
step = csf_i%stepvector(s)
! to get correct bValue use it for next spatial orbital..
tempWeight = getSingleMatrixElement(step, step, deltaB, gen, bVal)
end select
call update_matrix_element(t, tempWeight, 1)
call update_matrix_element(t, tempWeight, 2)
if (t_trunc_guga_matel) then
if (abs(extract_matrix_element(t, 1)) < trunc_guga_matel) then
probWeight = 0.0_dp
t = 0_n_int
end if
end if
end subroutine singleStochasticUpdate