subroutine calcFullStartRaisingStochastic(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 = "calcFullStartRaisingStochastic"
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(isZero(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
! but anyway need it for the final matrix element
! here in this case only this and the version with both index pairs
! exchange correspond to the reached excitation -> so the 1/2 in fron
! of the 2 electron part cancels in the hamiltonian.
umat = (get_umat_el(start, start, semi, ende) + &
get_umat_el(start, start, ende, semi)) / 2.0_dp
! but its not only this matrix element isnt it? have to also take the
! one with exchanged indices into account, but where only the type
! excitation determines the sign between the 2...
! also have to correctly index the routine with phycisist notation..
! check if i can deduce that by only the order entries of excitInfo..
! better do a excitation abortion here !TODO!! change that to the
! correct sum of umats
if (near_zero(umat)) then
branch_pgen = 0.0_dp
t = 0_n_int
return
end if
! set t
t = ilut
! set 0->3
set_orb(t, 2 * start)
set_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
select case (csf_i%stepvector(semi))
case (3)
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 both branches are zero, i have to abort the excitation
! altough that shouldnt happen...
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 3->1 first -1 branch first
clr_orb(t, 2 * semi)
! order does not matter since only x0 matrix element
call getDoubleMatrixElement(1, 3, deltaB, gen_type%R, gen_type%R, bVal, &
1.0_dp, tempWeight)
call setDeltaB(-1, t)
else
! do +1 branch
! then do 3->2: +1 branch(change from already set 1
clr_orb(t, 2 * semi - 1)
call getDoubleMatrixElement(2, 3, deltaB, gen_type%R, gen_type%R, bVal, &
1.0_dp, tempWeight)
call setDeltaB(1, t)
branch_pgen = 1.0_dp - branch_pgen
end if
else
! only 3->1, -1 branch possible
clr_orb(t, 2 * semi)
! order does not matter since only x0 matrix element
call getDoubleMatrixElement(1, 3, deltaB, gen_type%R, gen_type%R, 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 -> 0
clr_orb(t, 2 * semi - 1)
call getDoubleMatrixElement(0, 1, 0, gen_type%R, gen_type%R, bVal, &
1.0_dp, tempWeight)
call setDeltaB(1, t)
branch_pgen = 1.0_dp
case (2)
! here we have
! 2 -> 0
clr_orb(t, 2 * semi)
call getDoubleMatrixElement(0, 2, 0, gen_type%R, gen_type%R, bVal, &
1.0_dp, tempWeight)
call setDeltaB(-1, t)
branch_pgen = 1.0_dp
! encode fullstart contribution and pseudo overlap region here
! too in one go. one additional -1 due to semistop
end select
call encode_matrix_element(t, 0.0_dp, 2)
call encode_matrix_element(t, Root2 * tempWeight * (-1.0_dp)**nOpen, 1)
! x0 matrix elements cant be 0 at a RR semistop
! 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)
branch_pgen = branch_pgen * temp_pgen
! in the single overlap regions matrix elements cant actually
! become 0! -> so no check for matrix element is necessary here
! but i could be set to zero, due to both branches having
! zero probabilistic weight, which also should not happen, but just
! to be sure! this gets indicated by probWeight = 0 and t = 0
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 calcFullStartRaisingStochastic