subroutine createStochasticExcitation_double(ilut, nI, csf_i, excitation, pgen, excit_typ)
! calculate one possible double excitation and the corresponding
! probabilistic weight. and hamilton matrix element for a given CSF
integer(n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
integer, intent(in) :: nI(nel)
integer(n_int), intent(out) :: excitation(0:nifguga)
real(dp), intent(out) :: pgen
integer, intent(out) :: excit_typ(2)
character(*), parameter :: this_routine = "createStochasticExcitation_double"
type(ExcitationInformation_t) :: excitInfo
integer(n_int), allocatable :: excitations(:, :)
integer(n_int) :: ilutI(0:niftot), ilutJ(0:niftot)
real(dp) :: orb_pgen, branch_pgen
HElement_t(dp) :: mat_ele
type(WeightObj_t) :: weights
logical :: compFlag
real(dp) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
ASSERT(isProperCSF_ilut(ilut))
! do that in the calling interface funcition already..
! or maybe even in the FciMCPar to use the same b and occvector
! for and occupied determinant/CSF
call pickOrbitals_double(ilut, nI, csf_i, excitInfo, orb_pgen)
! check if orbitals were correctly picked
if ( .not. excitInfo%valid ) then
excitation = 0_n_int
pgen = 0.0_dp
return
end if
! do i need the checkcomp flag here?? how often does it happen that
! i create a wrong excitation information? can i avoid to create
! a wrong excitation information and thus not use checkComp here?
! profiler tells me this takes a lot of time..
! and if i have to do it, i could get rid of all the other uncessary
! call of calcRemainingSwitches ..
!todo: since the weights also get initialized within the
!checkcompatibility function, i should maybe output it also, similar to
!the posSwitches and negSwitches quantitites..
! or just dont do a compatibility check, since it will get aborted
! anyway if it does not work in the excitation generation..
call checkCompatibility(&
csf_i, excitInfo, compFlag, posSwitches, negSwitches, weights)
if (.not. compFlag) then
excitation = 0
pgen = 0.0_dp
return
end if
if (t_guga_back_spawn) then
if (increase_ex_levl(csf_i, excitInfo) .and. .not. is_init_guga) then
if (t_guga_back_spawn_trunc) then
pgen = 0.0_dp
excitation = 0_n_int
return
end if
call create_crude_guga_double(ilut, nI, csf_i, excitation, branch_pgen, excitInfo)
pgen = orb_pgen * branch_pgen
return
end if
end if
if (tgen_guga_crude .and. .not. tgen_guga_mixed) then
! do the crude approximation here, where i do not switch
! in the excitation range, but only at the picked electrons
! and orbitals
! this includes the change, that I always switch at mixed
! start and ends too!
call create_crude_double(ilut, excitation, branch_pgen, excitInfo)
if (near_zero(branch_pgen)) then
excitation = 0_n_int
pgen = 0.0_dp
return
end if
call convert_ilut_toNECI(ilut, ilutI)
call convert_ilut_toNECI(excitation, ilutJ)
call calc_guga_matrix_element(ilutI, csf_i, ilutJ, CSF_Info_t(ilutJ), excitInfo, mat_ele, .true.)
if (near_zero(mat_ele)) then
excitation = 0
pgen = 0.0_dp
return
end if
call encode_matrix_element(excitation, 0.0_dp, 2)
call encode_matrix_element(excitation, mat_ele, 1)
pgen = orb_pgen * branch_pgen
return
end if
! depending on the excitation chosen -> call specific stochastic
! excitation calculators. similar to the exact calculation
! only certain cases, where the excitation really can not be
! described by a single excitation, should arrive here.
select case (excitInfo%typ)
case (excit_type%single_overlap_L_to_R)
! single overlap lowering into raising
! similar to a single excitation except the (predetermined)
! single overlap site.
call calcSingleOverlapMixedStochastic(ilut, csf_i, excitInfo, excitation, &
branch_pgen, posSwitches, negSwitches, weights)
pgen = orb_pgen * branch_pgen
case (excit_type%single_overlap_R_to_L) ! single overlap raising into lowering
! similar to a single excitation except the (predetermined)
! single overlap site.
!todo: mention on the weight input here: in some routines below,
! the weights get reinitialized in the different sectors! so be
! careful to just input the weights everywhere, and also check the
! checkCompatibility function, if the weights get reinitialized
! there correctly!
call calcSingleOverlapMixedStochastic(ilut, csf_i, excitInfo, excitation, &
branch_pgen, posSwitches, negSwitches, weights)
pgen = orb_pgen * branch_pgen
case (excit_type%double_lowering) ! normal double two lowering
call calcDoubleLoweringStochastic(ilut, csf_i, excitInfo, excitation, branch_pgen, &
posSwitches, negSwitches, weights)
pgen = orb_pgen * branch_pgen
case (excit_type%double_raising) ! normal double two raising
call calcDoubleRaisingStochastic(ilut, csf_i, excitInfo, excitation, branch_pgen, &
posSwitches, negSwitches, weights)
pgen = orb_pgen * branch_pgen
case (excit_type%double_L_to_R_to_L) ! lowergin into raising into lowering
! should be able to use the same general function as above to
! calculate the excitation, but the matrix element calculation
! should be a little bit different... maybe additional input needed
call calcDoubleL2R2L_stochastic(ilut, csf_i, excitInfo, excitation, branch_pgen, &
posSwitches, negSwitches, weights)
pgen = orb_pgen * branch_pgen
case (excit_type%double_R_to_L_to_R) ! raising into lowering into raising
! dito about matrix elements as above...
call calcDoubleR2L2R_stochastic(ilut, csf_i, excitInfo, excitation, branch_pgen, &
posSwitches, negSwitches, weights)
pgen = orb_pgen * branch_pgen
case (excit_type%double_L_to_R) ! lowering into raising double
call calcDoubleL2R_stochastic(ilut, csf_i, excitInfo, excitation, branch_pgen, &
posSwitches, negSwitches, weights)
pgen = orb_pgen * branch_pgen
case (excit_type%double_R_to_L) ! raising into lowering double
call calcDoubleR2L_stochastic(ilut, csf_i, excitInfo, excitation, branch_pgen, &
posSwitches, negSwitches, weights)
pgen = orb_pgen * branch_pgen
case (excit_type%fullstop_lowering) ! full stop 2 lowering
! again the double overlap part is easy to deal with, since its
! only the deltaB=0 branch
call calcFullstopLoweringStochastic(ilut, csf_i, excitInfo, excitation, &
branch_pgen, posSwitches, negSwitches, weights)
pgen = orb_pgen * branch_pgen
case (excit_type%fullstop_raising) ! full-stop 2 raising
! again only deltaB = 0 branch in DE overlap region
call calcFullstopRaisingStochastic(ilut, csf_i, excitInfo, excitation, &
branch_pgen, posSwitches, negSwitches, weights)
pgen = orb_pgen * branch_pgen
case (excit_type%fullstop_L_to_R) ! full-stop lowering into raising
if (t_crude_exchange .or. (t_crude_exchange_noninits .and. (.not. is_init_guga))) then
! in case of "crude" exchange excitation perform a
! determinant-like excitation without spin-recouplings
! but only for non-inits.. so this information has to
! be passed in here!
call perform_crude_excitation(ilut, csf_i, excitInfo, excitation, compFlag)
! in this case the pgen is just the orbital pgen, as only
! on CSF can be created from it..
! but be sure if the excitation is then actually possible
if (.not. compFlag) then
excitation = 0_n_int
pgen = 0.0_dp
return
else
pgen = orb_pgen
end if
else
call calcFullStopL2R_stochastic(ilut, csf_i, excitInfo, excitation, &
branch_pgen, posSwitches, negSwitches, weights)
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
pgen = branch_pgen * orb_pgen
else
pgen = branch_pgen
end if
end if
case (excit_type%fullstop_R_to_L) ! full-stop raising into lowering
if (t_crude_exchange .or. (t_crude_exchange_noninits .and. (.not. is_init_guga))) then
call perform_crude_excitation(ilut, csf_i, excitInfo, excitation, compFlag)
! in this case the pgen is just the orbital pgen, as only
! on CSF can be created from it..
! but be sure if the excitation is then actually possible
if (.not. compFlag) then
excitation = 0_n_int
pgen = 0.0_dp
return
else
pgen = orb_pgen
end if
else
call calcFullStopR2L_stochastic(ilut, csf_i, excitInfo, excitation, &
branch_pgen, posSwitches, negSwitches, weights)
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
pgen = branch_pgen * orb_pgen
else
pgen = branch_pgen
end if
end if
case (excit_type%fullstart_lowering) ! full-start 2 lowering
! again only deltaB = 0 branch in DE overlap region
call calcFullStartLoweringStochastic(ilut, csf_i, excitInfo, excitation, &
branch_pgen, posSwitches, negSwitches, weights)
pgen = orb_pgen * branch_pgen
case (excit_type%fullstart_raising) ! full-start 2 raising
! the double overlap part is again really easy here, since only
! the deltaB=0 branch is non-zero -> and the second part can be
! treated as a single excitation
call calcFullStartRaisingStochastic(ilut, csf_i, excitInfo, excitation, &
branch_pgen, posSwitches, negSwitches, weights)
pgen = orb_pgen * branch_pgen
case (excit_type%fullStart_L_to_R) ! full-start lowering into raising
if (t_crude_exchange .or. (t_crude_exchange_noninits .and. (.not. is_init_guga))) then
call perform_crude_excitation(ilut, csf_i, excitInfo, excitation, compFlag)
! in this case the pgen is just the orbital pgen, as only
! on CSF can be created from it..
! but be sure if the excitation is then actually possible
if (.not. compFlag) then
excitation = 0_n_int
pgen = 0.0_dp
return
else
pgen = orb_pgen
end if
else
call calcFullStartL2R_stochastic(ilut, csf_i, excitInfo, excitation, &
branch_pgen, posSwitches, negSwitches, weights)
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
pgen = branch_pgen * orb_pgen
else
pgen = branch_pgen
end if
end if
case (excit_type%fullstart_R_to_L) ! full-start raising into lowering
if (t_crude_exchange .or. (t_crude_exchange_noninits .and. (.not. is_init_guga))) then
call perform_crude_excitation(ilut, csf_i, excitInfo, excitation, compFlag)
! in this case the pgen is just the orbital pgen, as only
! on CSF can be created from it..
! but be sure if the excitation is then actually possible
if (.not. compFlag) then
excitation = 0_n_int
pgen = 0.0_dp
return
else
pgen = orb_pgen
end if
else
call calcFullStartR2L_stochastic(ilut, csf_i, excitInfo, excitation, &
branch_pgen, posSwitches, negSwitches, weights)
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
pgen = branch_pgen * orb_pgen
else
pgen = branch_pgen
end if
end if
! start, case by case how they appear in the documentary
case (excit_type%fullstart_stop_alike) ! full-start into full-stop alike
! since there is only the deltaB = 0 branch with non-zero weight
! only 1 excitation is possible, which is calculated by the
! exact version
! check matrix element before calculating anything
if (near_zero(get_umat_el(excitInfo%i, excitInfo%i, excitInfo%j, excitInfo%j))) then
excitation = 0_n_int
pgen = 0.0_dp
return
end if
call calcFullStartFullStopAlike(ilut, csf_i, excitInfo, excitations)
excitation = excitations(:, 1)
! have to encode the umat/2 matrix element here to keep this
! above function compatible with the exact routines
call update_matrix_element(excitation, get_umat_el(excitInfo%i, &
excitInfo%i, excitInfo%j, excitInfo%j) / 2.0_dp, 1)
! probWeight is just the index choosing probability
! multiply further down
pgen = orb_pgen
! can a full-start-full stop alike excitation have zero
! matrix element?
! yes since D(b=1,0) and d(b=0,1) can be zero
! to use t_trunc_guga_pgen
branch_pgen = 1.0_dp
deallocate(excitations)
case (excit_type%fullstart_stop_mixed) ! full-start into full-stop mixed
! here it is garuanteed that its a open orbital at the start and
! end to allow for an non-diagonal excitation.
! but somehow i have to ensure, that the deltaB=0 path is left
! at somepoint... -> chances are slim, but there.. maybe can bias
! for that.. since i know how many switche possibilities are
! left and could include that somehow... todo
! on another note... since i know, that i have to switch at some
! point, i could totally ignore the x0-matrix element since it
! will be zero in the event of switching...
! i also need this behavior in the case of mixed full starts and
! full stops... -> so maybe write a specific double update function
! for that ...
! would have to include some alreadySwitched flag in the excitation
! to see and save if a switch already happened to enforce if
! necessary..
! or use the x0 matrix element as a kind of flag...
! since if its 0 it means a switch happended at some point, but
! thats seems a bit inefficient.
if (t_crude_exchange .or. (t_crude_exchange_noninits .and. (.not. is_init_guga))) then
call perform_crude_excitation(ilut, csf_i, excitInfo, excitation, compFlag)
! in this case the pgen is just the orbital pgen, as only
! on CSF can be created from it..
! but be sure if the excitation is then actually possible
if (.not. compFlag) then
excitation = 0_n_int
pgen = 0.0_dp
return
else
pgen = orb_pgen
end if
else
call calcFullStartFullStopMixedStochastic(ilut, csf_i, excitInfo, &
excitation, branch_pgen, posSwitches, negSwitches, weights)
if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
pgen = branch_pgen * orb_pgen
else
pgen = branch_pgen
end if
end if
! random notes:
! those are the "easy" ones..., which dont actually have something
! to do in the DE overlap region.. well except the full-start into
! full-stop mixed excitation, where we actually HAVE to do a switch or
! else we get a already dealt with single excitation
! UPDATE: see the picking of the full-start or stop orbital as picking
! the first or last, necessary, stepvector switch. since it has to be
! switched somewhere anyway
! these are all the possibilities for (ii,jj), (ii,jk) index picking.
! a thinking still is to maybe combine all of them into one orbital picker
! and allow the second picked orbital to be the same as the first , so
! also this would account for the correct probability to pick to similar
! ones.. (a little bit of renormalization is needed, since the order
! how to pick the orbitals shouldnt matter and thus has to be accounted
! for. after here 4 differeing (ij,kl) indices were picked:
! that should be about all...
end select
! what if probWeight is 0 for some reason? shouldnt be..
! yes it could be since i indicate zero-values excitations in this way
if (t_trunc_guga_pgen .or. (t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
if (branch_pgen < trunc_guga_pgen) then
pgen = 0.0_dp
excitation = 0_n_int
return
end if
end if
! check if for some reason the matrix element of the excitation is 0
if (near_zero(extract_h_element(excitation))) then
pgen = 0.0_dp
excitation = 0
else
! also store information of type of excitation for automated tau-search
! for the non-weighted guga-excitation-generator
select case (excitInfo%excitLvl)
case (0)
! (ii,jj) RR/LL excitation
excit_typ(1) = 2
excit_typ(2) = 1
case (1)
! (ii,jj) RL excitation
excit_typ(1) = 2
excit_typ(2) = 0
case (2)
! (ii,jk) RR/LL excitation
excit_typ(1) = 3
excit_typ(2) = 1
case (3)
! (ii,jk) RL excitation
excit_typ(1) = 3
excit_typ(2) = 0
case (4)
! (ij,kl) excitation
excit_typ(1) = 4
excit_typ(2) = 0
case default
call print_excitInfo(excitInfo)
call stop_all(this_routine, "wrong excit level info!")
end select
end if
select case(excitInfo%typ)
case (excit_type%single_overlap_L_to_R, &
excit_type%single_overlap_R_to_L, &
excit_type%double_lowering, &
excit_type%double_raising, &
excit_type%double_L_to_R_to_L, &
excit_type%double_R_to_L_to_R, &
excit_type%double_L_to_R, &
excit_type%double_R_to_L, &
excit_type%fullstop_raising, &
excit_type%fullstop_lowering, &
excit_type%fullstart_lowering, &
excit_type%fullstart_raising, &
excit_type%fullstart_stop_alike)
! in the other cases global_excitInfo gets assigned before it
! gets changed
global_excitInfo = excitInfo
end select
end subroutine createStochasticExcitation_double