function calc_pgen_back_spawn(nI, ilutI, ex, ic, part_type) result(pgen)
! to use HPHF keyword and also the test if the pgens are correct
! or just to be able and to be sure and more save i need a way to
! recalculate the generation probability also for a back-spawn
! excitation from a non-initiator determinant
! so although thats a hassle implement it, otherwise i cannot be
! quite sure about the method
integer, intent(in) :: nI(nel), ex(2, 2), ic, part_type
integer(n_int), intent(in) :: ilutI(0:niftot)
real(dp) :: pgen
integer :: dummy, ssrc, stgt, cc_index, src(2), tgt(2), dummy_elecs(2), &
dummy_orbs(2), ispn, loc, sum_ml, sym_prod, cc_a, cc_b, &
dummy_ispn, dummy_sum_ml
real(dp) :: elec_pgen, orb_pgen, int_cpt(2), cum_sum(2), cpt_pair(2), &
sum_pair(2), cum_arr(nbasis)
logical :: t_gen_list, t_in_ref, t_par
! i should only call this function when i am on a non-inititator or?
! in the HPHF framework i could end up here also on a non-initiator
! so here i should then go to 4ind-2 pgen calculator if it is a
! inititator
if (test_flag(ilutI, get_initiator_flag(part_type))) then
pgen = calc_pgen_4ind_weighted2(nI, ilutI, ex, ic)
else
! in the hphf framework we definetly only end up here if the
! back-spawn method is already turned on.. but is this
! ensured in the other places the function might get called?
! decision: if we call this function we want to get the pgen in
! the back-spawn method, so ensure outside if it is turned on
! or not.
if (ic == 1) then
! depending on the implementation
! here since i want to calculate the real pgen if back_spawn
! is or would be on, i should work with the
! _option keywords..
ssrc = ex(1, 1)
stgt = ex(2, 1)
if (t_back_spawn_flex_option) then
elec_pgen = 1.0_dp / real(nel, dp)
loc = check_electron_location([ssrc, 0], 1, part_type)
else
! it is slow anyway.. so just do the dummy implementation
! as Ali mentioned in the HPHF implementation it is not
! that common to have a doubly connected HPHF det
call pick_virtual_electron_single(nI, part_type, dummy, elec_pgen, .true.)
loc = -1
end if
if ((t_back_spawn_occ_virt) .or. (t_back_spawn_flex .and. &
((loc == 2 .and. occ_virt_level /= -1) .or. occ_virt_level == 2))) then
! the symmetry of the orbital is known after all
cc_index = ClassCountInd(get_spin(stgt), SpinOrbSymLabel(stgt), &
G1(stgt)%Ml)
! reuse the routine..
call pick_occupied_orbital_single(ilutI, ssrc, cc_index, &
part_type, orb_pgen, dummy, .true.)
else
! reuse the other pgen single function
! but i have to multiply out the electron pgen
orb_pgen = pgen_single_4ind(nI, ilutI, ssrc, stgt) * &
real(nel, dp)
end if
pgen = pSingles * elec_pgen * orb_pgen
else if (ic == 2) then
src = get_src(ex)
ispn = get_ispn(src)
sum_ml = sum(G1(src)%ml)
sym_prod = RandExcitSymLabelProd(SpinOrbSymLabel(src(1)), &
SpinOrbSymLabel(src(2)))
! i have to be careful with the orbitals..
! because i guess those get ordered.. and i have to check if
! it is possible to have picked the orbitals in either
! order..
! ok.. and there is the restriction to pick a beta orbital
! in src(1) first if it is a anti-parallel excitation
! ok this means atleast src(1) is definetly the first picked
! orbital.. is this also the case in HPHF?? argh i hate that
! stuff
! do this testing here once
! todo: i have to fix this since here i do not know anymore
! what was orb a and orb b since ex is sorted..
tgt = get_tgt(ex)
t_in_ref = (is_in_ref(tgt(1), part_type) .and. is_in_ref(tgt(2), part_type))
t_par = (is_beta(tgt(1)) .eqv. is_beta(tgt(2)))
! now i know.. for opposite spin excitations the beta orbital
! of tgt was the first picked one! so it would be best to
! switch it back, so i can correctly recalculate the
! pgen.. for parallel spin excitations the order does not
! matter
if (.not. t_par) then
if (.not. is_beta(tgt(1))) then
tgt = [tgt(2), tgt(1)]
end if
end if
if (t_back_spawn_flex) then
elec_pgen = pgen_weighted_elecs(nI, src)
loc = check_electron_location(src, 2, part_type)
else
call pick_virtual_electrons_double(nI, part_type, dummy_elecs, &
dummy_orbs, dummy_ispn, dummy_sum_ml, elec_pgen, .true.)
loc = -1
end if
! for some back_spawn_flex it can happen that we have no
! restrictions on the orbitals.. but thats rare i guess..
if (t_back_spawn_occ_virt .or. (t_back_spawn_flex .and. &
((loc == 1 .and. occ_virt_level /= -1) .or. loc == 2 .or. &
(loc == 0 .and. occ_virt_level >= 1)))) then
! in this case it matters which orbital was picked
! first..
! in this case we can be sure that atleast on of the
! orbitals is in the reference..
if (t_par) then
if (.not. is_in_ref(tgt(1), part_type)) then
! then it is incorrectly ordered.. reorder!
tgt = [tgt(2), tgt(1)]
end if
end if
call pick_occupied_orbital(ilutI, src, ispn, &
part_type, int_cpt(1), cum_sum(1), dummy_orbs(1), .true.)
! i can atleast do some stuff for picking it the other
! way or?
! because if it is a parallel spin-excitations and orbital
! (b) is in the reference the probs p(b|ij) are the same
if (t_par .and. t_in_ref) then
cpt_pair(1) = int_cpt(1)
sum_pair(1) = cum_sum(1)
else
cpt_pair = 0.0_dp
sum_pair = 1.0_dp
! maybe i could set a flag that it does not have to
! be recalced otherwise..
end if
! here i should go on to calculate p(b|aij) since in the
! other case we have everything..
if (t_back_spawn_flex .and. ( &
(loc == 2 .and. occ_virt_level /= -1) .or. occ_virt_level == 2)) then
! this is the case where orbital (b) is also restricted
cc_a = ClassCountInd(tgt(1))
cc_b = get_paired_cc_ind(cc_a, sym_prod, sum_ml, ispn)
call pick_second_occupied_orbital(ilutI, cc_b, tgt(1), &
ispn, part_type, int_cpt(2), cum_sum(2), dummy_orbs(2), .true.)
! and ofc both probs are the same if the spin-fits
if (t_par) then
cpt_pair = int_cpt
sum_pair = cum_sum
else
cpt_pair = 0.0_dp
sum_pair = 1.0_dp
end if
else
! the order should not matter or??
! or does it? and i always force a beta orbital
! to be picked first.. hm..
! the order does matter! and the excitation matrix
! is always ordered at this point but during picking
! the orbitals it is not!
! so we have to be more careful here!
call pgen_select_orb(ilutI, src, tgt(1), tgt(2), int_cpt(2), &
cum_sum(2))
! in this case (a) was restricted but (b) was not
! so check if the other way would have been
! possible
if (t_par .and. t_in_ref) then
! p(b|ij) already set above
call pgen_select_orb(ilutI, src, tgt(2), tgt(1), &
cpt_pair(2), sum_pair(2))
else
cpt_pair = 0.0_dp
sum_pair = 1.0_dp
end if
end if
else
call pgen_select_a_orb(ilutI, src, tgt(1), ispn, int_cpt(1), &
cum_sum(1), cum_arr, .true.)
! but i guess in this case i can already cover all the
! pgen..
if (int_cpt(1) > EPS) then
call pgen_select_orb(ilutI, src, tgt(1), tgt(2), &
int_cpt(2), cum_sum(2))
t_gen_list = .false.
else
t_gen_list = .false.
int_cpt = 0.0_dp
cum_sum = 1.0_dp
end if
! otherwise i can pick orbital (a) freely.. which also
! means that i definetly can pick (b) freely and do it
! the other way around..
call pgen_select_a_orb(ilutI, src, tgt(2), ispn, cpt_pair(1), &
sum_pair(1), cum_arr, t_gen_list)
if (cpt_pair(1) > EPS) then
call pgen_select_orb(ilutI, src, tgt(2), tgt(1), &
cpt_pair(2), sum_pair(2))
else
cpt_pair = 0.0_dp
sum_pair = 1.0_dp
end if
end if
! how do we handle this now??
! i only want to handle these cases, which i have not
! recalculated above already. especially for the p(b|ij)
! probability..
! duh.. i have to actually do something with the above
! calculated pgens..
! now i have to figure p(ab), p(ba) stuff.. and implement this
! above.. annoying..
if (any(cum_sum < EPS)) then
int_cpt = 0.0_dp
cum_sum = 1.0_dp
end if
if (any(sum_pair < EPS)) then
cpt_pair = 0.0_dp
sum_pair = 1.0_dp
end if
pgen = pDoubles * elec_pgen * (product(int_cpt) / product(cum_sum) + &
product(cpt_pair) / product(sum_pair))
else
pgen = 0.0_dp
end if
end if
end function calc_pgen_back_spawn