calc_pgen_back_spawn Function

public function calc_pgen_back_spawn(nI, ilutI, ex, ic, part_type) result(pgen)


Type IntentOptional Attributes Name
integer, intent(in) :: nI(nel)
integer(kind=n_int), intent(in) :: ilutI(0:niftot)
integer, intent(in) :: ex(2,2)
integer, intent(in) :: ic
integer, intent(in) :: part_type

Return Value real(kind=dp)


Source Code

Source Code

    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)

            ! 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)

                    ! 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), &

                    ! reuse the routine..
                    call pick_occupied_orbital_single(ilutI, ssrc, cc_index, &
                                                      part_type, orb_pgen, dummy, .true.)


                    ! 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)), &

                ! 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)


                    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)

                        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
                            cpt_pair = 0.0_dp
                            sum_pair = 1.0_dp
                        end if

                        ! 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), &

                        ! 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))

                            cpt_pair = 0.0_dp
                            sum_pair = 1.0_dp
                        end if

                    end if


                    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.
                        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))
                        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))


                pgen = 0.0_dp

            end if
        end if

    end function calc_pgen_back_spawn