subroutine calcFullStartLoweringStochastic(ilut, csf_i, excitInfo, t, branch_pgen, &
posSwitches, negSwitches, opt_weight)
! in this case there is no ambiguity in the matrix elements, as they
! are uniquely determined and thus can be efficiently calculated on
! the fly
integer(n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(in) :: excitInfo
integer(n_int), intent(out) :: t(0:nifguga)
real(dp), intent(out) :: branch_pgen
real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
type(WeightObj_t), intent(in), optional :: opt_weight
character(*), parameter :: this_routine = "calcFullStartLoweringStochastic"
real(dp) :: tempWeight, minusWeight, plusWeight, nOpen, bVal, temp_pgen
HElement_t(dp) :: umat
integer :: start, ende, semi, gen, iOrb, deltaB
type(WeightObj_t) :: weights
ASSERT(isThree(ilut, excitInfo%fullStart))
ASSERT(isProperCSF_ilut(ilut))
! create the fullStart
start = excitInfo%fullStart
ende = excitInfo%fullEnd
semi = excitInfo%firstEnd
gen = excitInfo%firstGen
bVal = csf_i%B_real(semi)
! could check U matrix element here.. although in the final version
! of the orbital picker, it probably should be accessed already
! for the cauchy-schwarz criteria
! in fullstart lowering the x1-elements are always zero -> so no
! minus sign in the integral contribution -> so just add the two
! relevant elements!
umat = (get_umat_el(ende, semi, start, start) + &
get_umat_el(semi, ende, start, start)) / 2.0_dp
! todo! correct combination of umat sum terms.. and correct indexing
if (near_zero(umat)) then
branch_pgen = 0.0_dp
t = 0_n_int
return
end if
! set t
t = ilut
! set 3->0
clr_orb(t, 2 * start)
clr_orb(t, 2 * start - 1)
! double overlap region influence only determined by number of open
! orbitals
nOpen = real(count_open_orbs_ij(csf_i, start, semi - 1), dp)
deltaB = getDeltaB(t)
if (present(opt_weight)) then
weights = opt_weight
else
weights = init_singleWeight(csf_i, ende)
end if
! could encode matrix element here, but do it later for more efficiency
tempWeight = 0.0_dp
select case (csf_i%stepvector(semi))
case (0)
if (csf_i%B_int(semi) > 0) then
minusWeight = weights%proc%minus(negSwitches(semi), bVal, weights%dat)
plusWeight = weights%proc%plus(posSwitches(semi), bVal, weights%dat)
if (near_zero(minusWeight + plusWeight)) then
branch_pgen = 0.0_dp
t = 0_n_int
return
end if
branch_pgen = calcStartProb(minusWeight, plusWeight)
if (genrand_real2_dSFMT() < branch_pgen) then
! do -1 branch
! do 0->1 first -1 branch first
set_orb(t, 2 * semi - 1)
! order does not matter since only x0 matrix element
call getDoubleMatrixElement(1, 0, deltaB, gen_type%L, gen_type%L, bVal, &
1.0_dp, tempWeight)
call setDeltaB(-1, t)
else
! do +1 branch
! then do 0->2: +1 branch(change from already set 1
set_orb(t, 2 * semi)
call getDoubleMatrixElement(2, 0, deltaB, gen_type%L, gen_type%L, bVal, &
1.0_dp, tempWeight)
call setDeltaB(1, t)
branch_pgen = 1.0_dp - branch_pgen
end if
else
! only 0->1, -1 branch possible
set_orb(t, 2 * semi - 1)
! order does not matter since only x0 matrix element
call getDoubleMatrixElement(1, 0, deltaB, gen_type%L, gen_type%L, bVal, &
1.0_dp, tempWeight)
call setDeltaB(-1, t)
branch_pgen = 1.0_dp
end if
case (1)
! only one excitation possible, which also has to have
! non-zero weight or otherwise i wouldnt even be here
! and reuse already provided t
! 1 -> 3
set_orb(t, 2 * semi)
call setDeltaB(1, t)
! still get the matrix elements here to deal with this
! generally
call getDoubleMatrixElement(3, 1, 0, gen_type%L, gen_type%L, bVal, &
1.0_dp, tempWeight)
branch_pgen = 1.0_dp
case (2)
! here we have
! 2 -> 3
set_orb(t, 2 * semi - 1)
call setDeltaB(-1, t)
call getDoubleMatrixElement(3, 2, 0, gen_type%L, gen_type%L, bVal, &
1.0_dp, tempWeight)
branch_pgen = 1.0_dp
! encode fullstart contribution and pseudo overlap region here
end select
call encode_matrix_element(t, 0.0_dp, 2)
call encode_matrix_element(t, Root2 * tempWeight * (-1.0_dp)**nOpen, 1)
! and then we have to do just a regular single excitation
do iOrb = semi + 1, ende - 1
call singleStochasticUpdate(ilut, csf_i, iOrb, excitInfo, weights, &
posSwitches, negSwitches, t, temp_pgen)
! matrix element cant be 0 but maybe both branch weights were..
branch_pgen = branch_pgen * temp_pgen
check_abort_excit(branch_pgen, t)
end do
call singleStochasticEnd(csf_i, excitInfo, t)
if (tFillingStochRDMOnFly) then
call encode_stochastic_rdm_info(GugaBits, t, rdm_ind= &
contract_2_rdm_ind(excitInfo%i, excitInfo%j, excitInfo%k, excitInfo%l, &
excit_lvl=2, excit_typ=excitInfo%typ), x1=0.0_dp, &
x0=extract_matrix_element(t, 1))
end if
call update_matrix_element(t, umat, 1)
end subroutine calcFullStartLoweringStochastic