pickOrbs_sym_uniform_mol_double Subroutine

public subroutine pickOrbs_sym_uniform_mol_double(ilut, nI, csf_i, excitInfo, pgen)

Arguments

Type IntentOptional Attributes Name
integer(kind=n_int), intent(in) :: ilut(0:GugaBits%len_tot)
integer, intent(in) :: nI(nel)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(out) :: excitInfo
real(kind=dp), intent(out) :: pgen

Contents


Source Code

    subroutine pickOrbs_sym_uniform_mol_double(ilut, nI, csf_i, excitInfo, pgen)
        ! new orbital picking routine, which is closer to simons already
        ! implemented one for the determinant version
        integer(n_int), intent(in) :: ilut(0:GugaBits%len_tot)
        integer, intent(in) :: nI(nel)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(out) :: excitInfo
        real(dp), intent(out) :: pgen
        character(*), parameter :: this_routine = "pickOrbs_sym_uniform_mol_double"

        integer :: occ_orbs(2), sym_prod, cc_a, cc_b, sum_ml, ind_res, &
                   st, en, a, b, i, j
        real(dp) :: int_contrib(2), cum_sum(2), cum_arr(nSpatOrbs), &
                    int_switch(2), cum_switch(2)
        logical :: range_flag

        ! has to be in the interface for function pointers
        unused_var(ilut)

        ! pick 2 ocupied orbitals randomly:
        call pick_elec_pair_uniform_guga(nI, occ_orbs, sym_prod, sum_ml, &
                                         pgen)
        ! the 2 picked electrons are ALWAYS ordered too! so i already know
        ! that I < J

        ! also store spatial orbitals:
        i = gtID(occ_orbs(1))
        j = gtID(occ_orbs(2))

        ! the occ_orbs are given in "spin-orbital" form.. not sure yet which
        ! implementation to choose!...

        ! then pick orbital a, weighted with FCIDUMP integrals
        call pick_a_orb_guga_mol(csf_i, occ_orbs, int_contrib(1), cum_sum(1), &
                                 cum_arr, a)

        ! changed that orbital a is now a spatial orbital already!!
        ! pick_a_orb now gives me a spatial orbital again!
        ! kill excitation if a == 0
        if (a == 0) then
            excitInfo%valid = .false.
            return
        end if

        ! now the additional GUGA and symmetry restrictions kick in..

        ! TODO: check if symmetries still work with GUGA so easily
        ! here always index it with the beta_spin part by default! and figure
        ! out how to later access it to not restrict any excitations on spin!
        ! change it to 2*a only.. since from now on i am ignoring spin
        ! symmetry anyways
        cc_a = ClassCountInd(2 * a)
        ! cc_a = ClassCountInd(1, SpinOrbSymLabel(2*a), G1(2*a)%Ml)
        ! hm have to think on how to use the symmetries of the spin-orbitals
        ! in the GUGA case..
        ! since original quanities sum_ml and iSpn are not really meaning
        ! anything in the CSF approach

        ! picking i and j in determinental case gives me information if my
        ! electrons have parallel or antiparallel spin ->
        ! depending on that my classcount for b depends on iSpn and the "spin"
        ! of orbitals b -> so i have to somehow exclude this info in the
        ! orbital choosing of b, and only consider spatial symmetries
        ! the spin restriction is not valid in the guga case since i
        ! essentially pick spatial orbitals -> to take into account both
        ! spin-orbital types -> choose cc_a to always take the beta
        ! orbital -> and choose both ispn = 2 and 3 to get all the
        ! possible spin_orbitals -> and then consruct a list of
        ! fitting spatial orbitals in the orbital picker!

        ! for (symorbs) quantity should not matter which spin the picked
        ! electrons have! (still have to check that thoroughly though!)
        ! so it should be enough to pick only one cc_index for b and
        ! convert to spatial orbitals then, and pick randomly from those!
        cc_b = get_paired_cc_ind(cc_a, sym_prod, sum_ml, 2)

        ! determine the GUGA restrictions on the orbitals!
        ! + to determine those restrictions also allows partially determining
        ! the typ of excitation already!
        ! restrctions on b can be:
        ! only exclude orbital (a) -> this is always the case! -> no additional
        ! exclude spin-orbital belonging to I or J
        ! b being strictly lower/higher than a certain orbital

        ! those two are essentially the only restrictions
        ! how to control that? do optional orbital input to the b orbital
        ! picker:
        ! if no additional input -> no additional restriction
        ! if positive integer input: exclude this specific orbital
        ! if integer + logical: if integer positive : excude everything above
        !                       if negative exclude everything below!

        ! do this + the symmetry and integral contributions, then i should
        ! be finished..

        ! and write the corresponding logic out here to determine the cases
        ! and to pre-determine the type of excitation

        ! now first check if the two spin_orb indices belong to the same
        ! spatial orbitals

        ! have to do the logic to determine the type of guga restrictions
        ! and determine the type of excitation here!
        ! and this also influences the pgen calculation!

        ! set default
        ind_res = 0
        range_flag = .false.

        if (is_in_pair(occ_orbs(1), occ_orbs(2))) then
            ! this also implies that there is no x1 matrix element contribution
            ! for the generator matrix elements -> always bias towards
            ! U + U' for orbital b
            ! this implies d(i) == 3 and there is NO additional restriction
            ! on hole a, excpet symmetries and stuff..
            ! in this case:
            ! full-start lowering and
            ! full-stop raising
            ! mixed raising into lowering single overlap
            ! excitations are possible
            ! depending where a and b are compared to I/J

            ! additional there are no restrictions on picking b
            ! todo: i could move the determination of the excitation
            ! information into the pick_b orb function. since there i
            ! already have to check the relation to the other already picked
            ! orbitals.. but not now i guess..
            call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, int_contrib(2), &
                                     cum_sum(2), b)
            ! TODO: still have to decide if i output SPATIAL or SPIN orbital
            ! for picked b... so i might have to convert at some point here
            ! and below!

            if (b == 0) then
                excitInfo%valid = .false.
                return
            end if
            ! changed, that b ouput is already a spatial orb..

            ! short names for stupid reasons
            en = max(a, b)
            st = min(a, b)

            ! also have to consider is a and b are from same spatial orb or?

            if (a == b) then
                if (a > i) then
                    !_LL_(ij) > ^LL^(ab)
                    excitInfo = assign_excitInfo_values_double( &
                                excit_type%fullstart_stop_alike, &
                                gen_type%L, gen_type%L, gen_type%L, gen_type%L, gen_type%L, &
                                b, i, b, i, i, i, b, b, 0, 2, 1.0_dp, 1.0_dp)

                else
                    ! _RR_(ab) > ^RR^(ij)
                    excitInfo = assign_excitInfo_values_double( &
                                excit_type%fullstart_stop_alike, &
                                gen_type%R, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                a, i, a, i, a, a, i, i, 0, 2, 1.0_dp, 1.0_dp)
                end if
            else
                if (a > i .and. b > i) then
                    ! _LL_(ij) > ^LL(min(a,b)) -> ^L(max(a,b))
                    excitInfo = assign_excitInfo_values_double( &
                                excit_type%fullstart_lowering, &
                                gen_type%L, gen_type%L, gen_type%L, gen_type%L, gen_type%L, &
                                en, i, st, i, i, i, st, en, 0, 2, 1.0_dp, 1.0_dp, 2)

                else if (a < i .and. b < i) then
                    ! _R(min(a,b)) > _RR(max(a,b) > ^RR^(ij)
                    excitInfo = assign_excitInfo_values_double( &
                                excit_type%fullstop_raising, &
                                gen_type%R, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                st, i, en, i, st, en, i, i, 0, 2, 1.0_dp, 1.0_dp, 2)

                else
                    ! _R(min(a,b)) > ^RL_(i) > ^L(max(a,b))
                    excitInfo = assign_excitInfo_values_double( &
                                excit_type%single_overlap_R_to_L, &
                                gen_type%R, gen_type%L, gen_type%R, gen_type%R, gen_type%L, &
                                st, i, en, i, st, i, i, en, 0, 2, 1.0_dp, 1.0_dp, 1)

                end if

            end if
            ! could have picked a and b in either way around..
            ! have p(a)*p(b|a) have to determine p(b)*p(a|b)
            ! can resuse cum_arrays form both picking
            ! no.. only first cum_arr -> have to reconstruct second one..
            ! but actually dont need an array
            ! do that below at the end since its always the same
            call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, a, &
                                          int_switch(2), cum_switch(2))

            ! should not happen but assert here that the cummulative
            ! probabilty is not 0
            ASSERT(.not. near_zero(cum_switch(2)))

        else
            ! if they are not in a pair there are more possibilites
            ! i know through the triangular mapping, that the 2 picked
            ! electrons are always ordered!
            if (csf_i%stepvector(i) == 3) then
                ! if the picked stepvector is doubly occupied, since i pick
                ! based on spinorbitals but then only consider spatial orbitals
                ! there are more chances to pick those orbitals then...
                ! but thats probably a not valid bias towards those type
                ! of excitations... -> have to think of a new way to pick
                ! occupied spatial orbitals
                !TODO: yes definetly not pick based on spin-orbitals then.
                ! i unnaturally pick towards doubly occupied sites.. hm..
                pgen = 2.0_dp * pgen
                if (csf_i%stepvector(j) == 3) then

                    ! see above for description
                    pgen = 2.0_dp * pgen

                    ! no additional restrictions in picking b
                    call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, &
                                             int_contrib(2), cum_sum(2), b)

                    if (b == 0) then
                        excitInfo%valid = .false.
                        return
                    end if

                    ! the good thing -> independent of the ordering, the pgens
                    ! have the same restrictions independent of the order how
                    ! a and b are picked!
                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, a, &
                                                  int_switch(2), cum_switch(2))

                    ! now determine the type of excitation:
                    if (a == b) then

                        if (a > j) then
                            ! _L(i) -> _LL(j) > ^LL^(ab)
                            excitInfo = assign_excitInfo_values_double( &
                                        excit_type%fullstop_lowering, &
                                        gen_type%L, gen_type%L, gen_type%L, gen_type%L, gen_type%L, &
                                        a, i, a, j, i, j, a, a, 0, 2, 1.0_dp, 1.0_dp, 2)

                        else if (a < i) then
                            ! _RR_(ab) > ^RR(i) > ^R(j)
                            excitInfo = assign_excitInfo_values_double( &
                                        excit_type%fullstart_raising, &
                                        gen_type%R, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                        a, i, a, j, a, a, i, j, 0, 2, 1.0_dp, 1.0_dp, 2)

                        else
                            ! _L(i) > ^LR_(ab) > ^R(j)
                            excitInfo = assign_excitInfo_values_double( &
                                        excit_type%single_overlap_L_to_R, &
                                        gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%R, &
                                        a, i, a, j, i, a, a, j, 0, 2, 1.0_dp, 1.0_dp, 1)

                        end if
                    else
                        ! only the maximums count..
                        en = max(a, b)
                        st = min(a, b)
                        if (en > j) then
                            if (st > j) then
                                ! _L(i) > _LL(j) > ^LL(min) > ^L(max)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%double_lowering, &
                                            gen_type%L, gen_type%L, gen_type%L, gen_type%L, gen_type%L, &
                                            st, i, en, j, i, j, st, en, 0, 4, 1.0_dp, 1.0_dp)
                            else if (st < i) then
                                ! _R(min) > _LR(i) > ^RL(j) > ^L(max)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%double_R_to_L, &
                                            gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%L, &
                                            st, j, en, i, st, i, j, en, 0, 4, 1.0_dp, 1.0_dp)
                            else
                                ! _L(i) > _RL(min) > ^RL(j) > ^L(max)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%double_L_to_R_to_L, &
                                            gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
                                            en, i, st, j, i, st, j, en, 0, 4, 1.0_dp, 1.0_dp)

                            end if
                        else if (en < i) then
                            ! this implies that st < en < I !
                            ! _R(min) > _RR(max) > ^RR(i) > ^R(j)
                            excitInfo = assign_excitInfo_values_double( &
                                        excit_type%double_raising, &
                                        gen_type%R, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                        st, j, en, i, st, en, i, j, 0, 4, 1.0_dp, 1.0_dp)

                        else
                            if (st < i) then
                                ! _R(min) > _LR(i) > ^LR(max) > ^R(j)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%double_R_to_L_to_R, &
                                            gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                            st, j, en, i, st, i, en, j, 0, 4, 1.0_dp, 1.0_dp)
                            else
                                ! _L(i) > _RL(min) > ^LR(max) > ^R(j)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%double_L_to_R, &
                                            gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%R, &
                                            en, i, st, j, i, st, en, j, 0, 4, 1.0_dp, 1.0_dp)
                            end if
                        end if
                    end if
                else
                    ! so its a [3 1] configuration to start with
                    ! depending on the position of a there are restrictions
                    ! on b
                    ! have to check if A == J
                    if (a == j) then
                        ! then b has to be strictly lower then j!
                        call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, &
                                                 int_contrib(2), cum_sum(2), b, -j, .true.)

                        if (b == 0) then
                            excitInfo%valid = .false.
                            return
                        end if

!
                        ! but for picking it the other way around there is no
                        ! restriction on the orbitals
                        ! although for the nasty mixed full-stops i have to
                        ! recalculate the pgens anyway..
                        call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, &
                                                      a, int_switch(2), cum_switch(2))

                        ! determine excit
                        ! ATTENTION for the 2 below orbital J must not be
                        ! doubly occupied or else its just a single-like
                        ! excitation!! -> check if thats the case!
                        ! since a only can be a empty orbital -> it only can be
                        ! a singly occupied orbital!
                        if (b < i) then
                            ! _R(b) > _LR(i) > ^RL^(ja)
                            excitInfo = assign_excitInfo_values_double( &
                                        excit_type%fullstop_R_to_L, &
                                        gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                        b, j, j, i, b, i, j, j, 0, 4, 1.0_dp, 1.0_dp)

                        else
                            ! _L(i) > _RL(b) > ^RL^(ja)
                            excitInfo = assign_excitInfo_values_double( &
                                        excit_type%fullstop_L_to_R, &
                                        gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
                                        j, i, b, j, i, b, j, j, 0, 4, 1.0_dp, 1.0_dp)

                        end if

                    else
                        ! if its not j
                        if (a > j) then
                            ! b can not be J!
                            call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, &
                                                     int_contrib(2), cum_sum(2), b, j)

                            if (b == 0) then
                                excitInfo%valid = .false.
                                return
                            end if

                            ! depending on where b is the excitation is defined
                            ! and the pgen influence from picking a <-> b
                            if (b > j) then
                                ! both are on top -> same pgen
                                ! p(b) is always determinable from cum_arr... do
                                ! it outside!
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, &
                                                              a, int_switch(2), cum_switch(2), j)

                                if (a == b) then
                                    ! _L(i) > _LL(j) > ^LL^(ab)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%fullstop_lowering, &
                                                gen_type%L, gen_type%L, gen_type%L, gen_type%L, gen_type%L, &
                                                a, i, a, j, i, j, a, a, 0, 4, 1.0_dp, 1.0_dp)

                                else
                                    st = min(a, b)
                                    en = max(a, b)
                                    ! _L(i) > _LL(j) > ^LL(min) > ^L(max)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%double_lowering, &
                                                gen_type%L, gen_type%L, gen_type%L, gen_type%L, gen_type%L, &
                                                st, i, en, j, i, j, st, en, 0, 4, 1.0_dp, 1.0_dp)

                                end if

                            else
                                ! pgen is not restricted
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, &
                                                              a, int_switch(2), cum_switch(2))
                                if (b < i) then
                                    ! _R(b) > _LR(i) > ^RL(j) > ^L(a)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%double_R_to_L, &
                                                gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%L, &
                                                b, j, a, i, b, i, j, a, 0, 4, 1.0_dp, 1.0_dp)

                                else
                                    ! _L(i) > _RL(b) > ^RL(j) > ^L(a)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%double_L_to_R_to_L, &
                                                gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
                                                a, i, b, j, i, b, j, a, 0, 4, 1.0_dp, 1.0_dp)

                                end if
                            end if
                        else
                            ! there is no restriction in picking b
                            call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, &
                                                     int_contrib(2), cum_sum(2), b)

                            if (b == 0) then
                                excitInfo%valid = .false.
                                return
                            end if
                            ! check where b is
                            if (b == j) then
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                              b, a, int_switch(2), cum_switch(2), &
                                                              -j, .true.)

                                if (a < i) then
                                    ! _R(a) > _LR(i) > ^RL^(jb)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%fullstop_R_to_L, &
                                                gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                a, j, j, i, a, i, j, j, 0, 2, 1.0_dp, 1.0_dp)

                                else
                                    ! _L(i) > _RL(a) > ^RL^(jb)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%fullstop_L_to_R, &
                                                gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
                                                j, i, a, j, i, a, j, j, 0, 2, 1.0_dp, 1.0_dp)

                                end if

                            else
                                if (b > j) then
                                    ! a could not have been j
                                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                  b, a, int_switch(2), cum_switch(2), j)
                                    if (a < i) then
                                        ! _R(a) > _LR(i) > ^RL(j) > ^L(b)
                                        excitInfo = assign_excitInfo_values_double( &
                                                    excit_type%double_R_to_L, &
                                                    gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%L, &
                                                    a, j, b, i, a, i, j, b, 0, 4, 1.0_dp, 1.0_dp)

                                    else
                                        ! _L(i) > _RL(a) > ^RL(j) > ^L(b)
                                        excitInfo = assign_excitInfo_values_double( &
                                                    excit_type%double_L_to_R_to_L, &
                                                    gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
                                                    b, i, a, j, i, a, j, b, 0, 4, 1.0_dp, 1.0_dp)

                                    end if
                                else
                                    ! no restric, on other order
                                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                  b, a, int_switch(2), cum_switch(2))

                                    if (a == b) then
                                        if (a < i) then
                                            ! _RR_(ab) > ^RR(i) > ^R(j)
                                            excitInfo = assign_excitInfo_values_double( &
                                                        excit_type%fullstart_raising, &
                                                        gen_type%R, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                        a, i, a, j, a, a, i, j, 0, 2, 1.0_dp, 1.0_dp)
                                        else
                                            ! _L(i) > ^LR_(ab) > ^R(j)
                                            excitInfo = assign_excitInfo_values_double( &
                                                        excit_type%single_overlap_L_to_R, &
                                                        gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%R, &
                                                        a, i, a, j, i, a, a, j, 0, 2, 1.0_dp, 1.0_dp, 1)
                                        end if
                                    else
                                        ! only min and max importtant
                                        st = min(a, b)
                                        en = max(a, b)

                                        if (en < i) then
                                            ! _R(min) > _RR(max) > ^RR(i) > ^R(j)
                                            excitInfo = assign_excitInfo_values_double( &
                                                        excit_type%double_raising, &
                                                        gen_type%R, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                        st, i, en, j, st, en, i, j, 0, 4, 1.0_dp, 1.0_dp)

                                        else
                                            if (st < i) then
                                                ! _R(min) > _LR(i) > ^LR(max) > ^R(j)
                                                excitInfo = assign_excitInfo_values_double( &
                                                            excit_type%double_R_to_L_to_R, &
                                                            gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                            st, j, en, i, st, i, en, j, 0, 4, 1.0_dp, 1.0_dp)

                                            else
                                                ! _L(i) > _RL(min) > ^LR(max) > ^R(j)
                                                excitInfo = assign_excitInfo_values_double( &
                                                            excit_type%double_L_to_R, &
                                                            gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%R, &
                                                            en, i, st, j, i, st, en, j, 0, 4, 1.0_dp, 1.0_dp)
                                            end if
                                        end if
                                    end if
                                end if
                            end if
                        end if
                    end if
                end if
            else
                if (csf_i%stepvector(j) == 3) then
                    pgen = 2.0_dp * pgen
                    ! its a [1 3] configuration -> very similar to [3 1]
                    if (a == i) then
                        ! b has to be strictly higher then I
                        call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, &
                                                 int_contrib(2), cum_sum(2), b, i, .true.)

                        if (b == 0) then
                            excitInfo%valid = .false.
                            return
                        end if
                        ! no restriction the other way around
                        call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, a, &
                                                      int_switch(2), cum_switch(2))

                        ! ATTENTION: here there have to be the additional
                        ! constraint that spat orb I must not be doubly
                        ! occupied! or otherwise its a single-like excitation!!
                        ! did i consider that?`
                        ! yes! it cant be a doubly occupied orbital since
                        ! orbital a and b can only be non -occupied orbitals!
                        if (b > j) then
                            ! _RL_(ia) > ^RL(j) > ^L(b)
                            excitInfo = assign_excitInfo_values_double( &
                                        excit_type%fullstart_R_to_L, &
                                        gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
                                        b, i, i, j, i, i, j, b, 0, 2, 1.0_dp, 1.0_dp)

                        else
                            ! _RL_(ia) > ^LR(b) > ^R(j)
                            excitInfo = assign_excitInfo_values_double( &
                                        excit_type%fullstart_L_to_R, &
                                        gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                        b, i, i, j, i, i, b, j, 0, 2, 1.0_dp, 1.0_dp)
                        end if

                    else
                        ! if its not i
                        if (a < i) then
                            ! b cannot be I
                            call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, &
                                                     int_contrib(2), cum_sum(2), b, i)

                            if (b == 0) then
                                excitInfo%valid = .false.
                                return
                            end if

                            if (b < i) then
                                ! same pgen restrictions
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                              b, a, int_switch(2), cum_switch(2), i)

                                if (a == b) then
                                    ! _RR_(ab) > ^RR(i) > ^R(j)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%fullstart_raising, &
                                                1, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                a, i, a, j, a, a, i, j, 0, 2, 1.0_dp, 1.0_dp)
                                else
                                    st = min(a, b)
                                    en = max(a, b)

                                    ! _R(min) > _RR(max) > ^RR(i) > ^R(j)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%double_raising, &
                                                gen_type%R, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                en, i, st, j, st, en, i, j, 0, 4, 1.0_dp, 1.0_dp)
                                end if
                            else
                                ! pgen is not restricted
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                              b, a, int_switch(2), cum_switch(2))

                                if (b > j) then
                                    ! _R(a) > _LR(i) > ^RL(j) > ^L(b)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%double_R_to_L, &
                                                gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%L, &
                                                a, j, b, i, a, i, j, b, 0, 4, 1.0_dp, 1.0_dp)

                                else
                                    ! _R(a) > _LR(i)> ^LR(b) > ^R(j)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%double_R_to_L_to_R, &
                                                gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                a, j, b, i, a, i, b, j, 0, 4, 1.0_dp, 1.0_dp)

                                end if
                            end if
                        else
                            ! there is no restriction on b
                            call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, &
                                                     int_contrib(2), cum_sum(2), b)

                            if (b == 0) then
                                excitInfo%valid = .false.
                                return
                            end if

                            ! check where b is:
                            if (b == i) then
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                              b, a, int_switch(2), cum_switch(2), i, .true.)

                                if (a > j) then
                                    ! _RL_(ib) > ^RL(j) > ^L(a)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%fullstart_R_to_L, &
                                                gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
                                                i, j, a, i, i, i, j, a, 0, 2, 1.0_dp, 1.0_dp)

                                else
                                    ! _RL_(ib) > ^LR(a) > ^R(j)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%fullstart_L_to_R, &
                                                gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                i, j, a, i, i, i, a, j, 0, 2, 1.0_dp, 1.0_dp)

                                end if

                            else
                                if (b < i) then
                                    ! a coud not have been i
                                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                  b, a, int_switch(2), cum_switch(2), i)

                                    if (a > j) then
                                        ! _R(b) > _LR(i) > ^RL(j) > ^L(a)
                                        excitInfo = assign_excitInfo_values_double( &
                                                    excit_type%double_R_to_L, &
                                                    gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%L, &
                                                    b, j, a, i, b, i, j, a, 0, 4, 1.0_dp, 1.0_dp)

                                    else
                                        ! _R(b) > _LR(i) > ^LR(a) > ^R(j)
                                        excitInfo = assign_excitInfo_values_double( &
                                                    excit_type%double_R_to_L_to_R, &
                                                    gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                    b, j, a, i, b, i, a, j, 0, 4, 1.0_dp, 1.0_dp)

                                    end if

                                else
                                    ! no restrictions on pgen
                                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                  b, a, int_switch(2), cum_switch(2))

                                    if (a == b) then
                                        if (a > j) then
                                            ! _L(i) > _LL(j) > ^LL^(ab)
                                            excitInfo = assign_excitInfo_values_double( &
                                                        excit_type%fullstop_lowering, &
                                                        gen_type%L, gen_type%L, gen_type%L, gen_type%L, gen_type%L, &
                                                        a, i, a, j, i, j, a, a, 0, 2, 1.0_dp, 1.0_dp)

                                        else
                                            ! _L(i) > ^LR_(ab) > ^R(j)
                                            excitInfo = assign_excitInfo_values_double( &
                                                        excit_type%single_overlap_L_to_R, &
                                                        gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%R, &
                                                        a, i, a, j, i, a, a, j, 0, 2, 1.0_dp, 1.0_dp, 1)

                                        end if

                                    else
                                        ! only extremes coun
                                        st = min(a, b)
                                        en = max(a, b)

                                        if (st > j) then
                                            ! _L(i) > _LL(j) > ^LL(min) > ^L(max)
                                            excitInfo = assign_excitInfo_values_double( &
                                                        excit_type%double_lowering, &
                                                        gen_type%L, gen_type%L, gen_type%L, gen_type%L, gen_type%L, &
                                                        st, i, en, j, i, j, st, en, 0, 4, 1.0_dp, 1.0_dp)

                                        else
                                            if (en > j) then
                                                ! _L(i) > _RL(min) > ^RL(j) > ^L(max)
                                                excitInfo = assign_excitInfo_values_double( &
                                                            excit_type%double_L_to_R_to_L, &
                                                            gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
                                                            en, i, st, j, i, st, j, en, 0, 4, 1.0_dp, 1.0_dp)

                                            else
                                                ! _L(i) > _RL(min) ^LR(max) > ^R(j)
                                                excitInfo = assign_excitInfo_values_double( &
                                                            excit_type%double_L_to_R, &
                                                            gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%R, &
                                                            en, i, st, j, i, st, en, j, 0, 4, 1.0_dp, 1.0_dp)

                                            end if
                                        end if
                                    end if
                                end if
                            end if
                        end if
                    end if
                else
                    ! [ 1 1 ] config
                    ! in this case, if i dont exclude these cases beforehand
                    ! i have to check if there is a possible switch if
                    ! (a) is (i) or (j)
                    if (a == j) then
                        if (csf_i%stepvector(i) == 1 .and. &
                            csf_i%stepvector(j) == 1) then
                            if (count_alpha_orbs_ij(csf_i, i, j) == 0) then
                                ! no valid excitation
                                excitInfo%valid = .false.
                                return
                            end if
                        else if (csf_i%stepvector(i) == 2 .and. &
                                 csf_i%stepvector(j) == 2) then
                            if (count_beta_orbs_ij(csf_i, i, j) == 0) then
                                excitInfo%valid = .false.
                                return
                            end if
                        end if
                        ! b has to be lower than J
                        call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, &
                                                 int_contrib(2), cum_sum(2), b, -j, .true.)

                        if (b == 0) then
                            excitInfo%valid = .false.
                            return
                        end if

                        ! have to check where b is to determine switchen
                        ! pgen contribution!
                        if (b == i) then
                            ! a would have to have been > I
                            call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, &
                                                          a, int_switch(2), cum_switch(2), i, .true.)

                            ! _RL_(ib) > ^RL^(ja)
                            excitInfo = assign_excitInfo_values_double( &
                                        excit_type%fullstart_stop_mixed, &
                                        gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                        i, j, j, i, i, i, j, j, 0, 2, 1.0_dp, 1.0_dp)

                        else
                            if (b < i) then
                                ! I is restricted for a
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, &
                                                              a, int_switch(2), cum_switch(2), i)

                                ! why am i never here??

                                ! _R(b) > _LR(i) > ^RL^(ja)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%fullstop_R_to_L, &
                                            gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                            b, j, j, i, b, i, j, j, 0, 2, 1.0_dp, 1.0_dp)

                            else
                                ! no restrictions
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, &
                                                              a, int_switch(2), cum_switch(2))

                                ! _L(i) > _RL(b) > ^RL^(ja)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%fullstop_L_to_R, &
                                            gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%R, &
                                            j, i, b, j, i, b, j, j, 0, 2, 1.0_dp, 1.0_dp)

                            end if
                        end if
                    else if (a == i) then
                        ! check if there is a possible switch if both i and j
                        ! have the same stepvalue
                        if (csf_i%stepvector(i) == 1 .and. &
                            csf_i%stepvector(j) == 1) then
                            if (count_alpha_orbs_ij(csf_i, i, j) == 0) then
                                ! no valid excitation
                                excitInfo%valid = .false.
                                return
                            end if
                        else if (csf_i%stepvector(i) == 2 .and. &
                                 csf_i%stepvector(j) == 2) then
                            if (count_beta_orbs_ij(csf_i, i, j) == 0) then
                                excitInfo%valid = .false.
                                return
                            end if
                        end if
                        ! b has to be higher than I
                        call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, &
                                                 int_contrib(2), cum_sum(2), b, i, .true.)

                        if (b == 0) then
                            excitInfo%valid = .false.
                            return
                        end if

                        ! check where b is

                        if (b == j) then
                            ! a would have to have been < J
                            call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, &
                                                          a, int_switch(2), cum_switch(2), -j, .true.)

                            ! _RL_(ia) > ^RL^(jb)
                            excitInfo = assign_excitInfo_values_double( &
                                        excit_type%fullstart_stop_mixed, &
                                        gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                        j, i, i, j, i, i, j, j, 0, 2, 1.0_dp, 1.0_dp)
                        else
                            if (b > j) then
                                ! J would be restricted
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, &
                                                              a, int_switch(2), cum_switch(2), j)

                                ! _RL_(ia) > ^RL(j) > ^L(b)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%fullstart_R_to_L, &
                                            gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
                                            i, j, b, i, i, i, j, b, 0, 2, 1.0_dp, 1.0_dp)

                            else
                                ! no restrictions
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, &
                                                              a, int_switch(2), cum_switch(2))

                                ! _RL_(ia) > ^LR(b) > ^R(j)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%fullstart_L_to_R, &
                                            gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                            i, j, b, i, i, i, b, j, 0, 2, 1.0_dp, 1.0_dp)

                            end if
                        end if
                    else
                        ! check were a is
                        if (a > j) then
                            ! b cant be J
                            call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, &
                                                     int_contrib(2), cum_sum(2), b, j)

                            if (b == 0) then
                                excitInfo%valid = .false.
                                return
                            end if

                            ! check where b is
                            if (b == i) then
                                ! check if there is a possible switch if both i and j
                                ! have the same stepvalue
                                if (csf_i%stepvector(i) == 1 .and. &
                                    csf_i%stepvector(j) == 1) then
                                    if (count_alpha_orbs_ij(csf_i, i, j) == 0) then
                                        ! no valid excitation
                                        excitInfo%valid = .false.
                                        return
                                    end if
                                else if (csf_i%stepvector(i) == 2 .and. &
                                         csf_i%stepvector(j) == 2) then
                                    if (count_beta_orbs_ij(csf_i, i, j) == 0) then
                                        excitInfo%valid = .false.
                                        return
                                    end if
                                end if

                                ! everything below I is restricted
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, &
                                                              a, int_switch(2), cum_switch(2), i, .true.)

                                ! _RL_(ib) > ^RL(j) > ^L(a)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%fullstart_R_to_L, &
                                            gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
                                            i, j, a, i, i, i, j, a, 0, 2, 1.0_dp, 1.0_dp)

                            else
                                if (b < i) then
                                    ! I would have been off limits
                                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                  b, a, int_switch(2), cum_switch(2), i)

                                    ! _R(b) > _LR(i) > ^RL(j) > ^L(a)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%double_R_to_L, &
                                                gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%L, &
                                                b, j, a, i, b, i, j, a, 0, 4, 1.0_dp, 1.0_dp)

                                else if (b > j) then

                                    if (a == b) then
                                        ! J would have been off limits, although in
                                        ! this case i could just copy the probs
                                        ! since a == b .. duh
                                        int_switch(2) = int_contrib(2)
                                        cum_switch(2) = cum_sum(2)

                                        ! and its a:
                                        ! _L(i) -> _LL(j) > ^LL^(a,b)
                                        excitInfo = assign_excitInfo_values_double( &
                                                    excit_type%fullstop_lowering, &
                                                    gen_type%L, gen_type%L, gen_type%L, gen_type%L, gen_type%L, &
                                                    a, i, a, j, i, j, a, a, 0, 2, 1.0_dp, 1.0_dp, 2)

                                    else

                                        ! only extremes count.. and J would have
                                        ! been off-limits
                                        call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                      b, a, int_switch(2), cum_switch(2), j)

                                        st = min(a, b)
                                        en = max(a, b)

                                        ! _L(i) > _LL(j) > ^LL(min) > ^L(max)
                                        excitInfo = assign_excitInfo_values_double( &
                                                    excit_type%double_lowering, &
                                                    gen_type%L, gen_type%L, gen_type%L, gen_type%L, gen_type%L, &
                                                    st, i, en, j, i, j, st, en, 0, 4, 1.0_dp, 1.0_dp)
                                    end if

                                else
                                    ! no restrictions
                                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                  b, a, int_switch(2), cum_switch(2))

                                    ! _L(i) > _RL(b) > ^RL(j) > ^L(a)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%double_L_to_R_to_L, &
                                                gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
                                                a, i, b, j, i, b, j, a, 0, 4, 1.0_dp, 1.0_dp)

                                end if
                            end if
                        else if (a < i) then
                            ! b cant be I
                            call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, &
                                                     int_contrib(2), cum_sum(2), b, i)

                            if (b == 0) then
                                excitInfo%valid = .false.
                                return
                            end if

                            ! check where b is
                            if (b == j) then
                                ! check if there is a possible switch if both i and j
                                ! have the same stepvalue
                                if (csf_i%stepvector(i) == 1 .and. &
                                    csf_i%stepvector(j) == 1) then
                                    if (count_alpha_orbs_ij(csf_i, i, j) == 0) then
                                        ! no valid excitation
                                        excitInfo%valid = .false.
                                        return
                                    end if
                                else if (csf_i%stepvector(i) == 2 .and. &
                                         csf_i%stepvector(j) == 2) then
                                    if (count_beta_orbs_ij(csf_i, i, j) == 0) then
                                        excitInfo%valid = .false.
                                        return
                                    end if
                                end if

                                ! wrong: I would have been off limits
                                ! the above comment is not right i would have
                                ! had to picked something below j, but why
                                ! is it 0? but thats atleast consistent with
                                ! above.. is the umat read in wrong?
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, b, &
                                                              a, int_switch(2), cum_switch(2), -j, .true.)

                                ! _R(a) > _LR(i) > ^RL^(jb)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%fullstop_R_to_L, &
                                            gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                            a, j, j, i, a, i, j, j, 0, 2, 1.0_dp, 1.0_dp)

                            else
                                if (b > j) then
                                    ! J would have been off-limits
                                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                  b, a, int_switch(2), cum_switch(2), j)
                                    ! _R(a) > _LR(i) > ^RL(j) > ^L(b)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%double_R_to_L, &
                                                gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%L, &
                                                a, j, b, i, a, i, j, b, 0, 4, 1.0_dp, 1.0_dp)

                                else if (b < i) then

                                    ! a can be b! why did i forget that...
                                    if (a == b) then
                                        int_switch(2) = int_contrib(2)
                                        cum_switch(2) = cum_sum(2)

                                        ! _RR_(ab) > ^RR(i) > ^R(j)
                                        excitInfo = assign_excitInfo_values_double( &
                                                    excit_type%fullstart_raising, &
                                                    gen_type%R, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                    a, i, a, j, a, a, i, j, 0, 2, 1.0_dp, 1.0_dp, 2)
                                    else

                                        ! only extremes and I would have been off
                                        ! lmits
                                        call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                      b, a, int_switch(2), cum_switch(2), i)

                                        st = min(a, b)
                                        en = max(a, b)

                                        ! _R(min) > _RR(max) > ^RR(i) > ^R(j)
                                        excitInfo = assign_excitInfo_values_double( &
                                                    excit_type%double_raising, &
                                                    gen_type%R, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                    en, j, st, i, st, en, i, j, 0, 4, 1.0_dp, 1.0_dp)
                                    end if

                                else
                                    ! no restrictions
                                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                  b, a, int_switch(2), cum_switch(2))

                                    ! _R(a) > _LR(i) > ^LR(b) > ^R(j)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%double_R_to_L_to_R, &
                                                gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                a, j, b, i, a, i, b, j, 0, 4, 1.0_dp, 1.0_dp)

                                end if
                            end if
                        else
                            ! a is between i and j -> no b restrictions
                            call pick_b_orb_guga_mol(csf_i, occ_orbs, a, cc_b, &
                                                     int_contrib(2), cum_sum(2), b)

                            if (b == 0) then
                                excitInfo%valid = .false.
                                return
                            end if
                            if (b < i) then
                                ! I off limits
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                              b, a, int_switch(2), cum_switch(2), i)

                                ! _R(b) > _LR(i) > ^LR(a) > ^R(j)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%double_R_to_L_to_R, &
                                            gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                            b, j, a, i, b, i, a, j, 0, 4, 1.0_dp, 1.0_dp)

                            else if (b > j) then
                                ! J off limits
                                call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                              b, a, int_switch(2), cum_switch(2), j)

                                ! _L(i) > _RL(a) > ^RL(j) > ^L(b)
                                excitInfo = assign_excitInfo_values_double( &
                                            excit_type%double_L_to_R_to_L, &
                                            gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%L, &
                                            b, i, a, j, i, a, j, b, 0, 4, 1.0_dp, 1.0_dp)

                            else

                                ! check where b is
                                if (a == b) then
                                    ! no restrictions on pgen
                                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                  b, a, int_switch(2), cum_switch(2))

                                    ! _L(i) > ^LR_(ab) > ^R(j)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%single_overlap_L_to_R, &
                                                gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%R, &
                                                a, i, a, j, i, a, a, j, 0, 2, 1.0_dp, 1.0_dp, 1)

                                else if (b == i) then
                                    ! a > I
                                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                  b, a, int_switch(2), cum_switch(2), i, .true.)

                                    ! _RL_(ib) > ^LR(a) > ^R(j)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%fullstart_L_to_R, &
                                                gen_type%L, gen_type%R, gen_type%R, gen_type%R, gen_type%R, &
                                                i, j, a, i, i, i, a, j, 0, 2, 1.0_dp, 1.0_dp)

                                else if (b == j) then
                                    ! a < J
                                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                  b, a, int_switch(2), cum_switch(2), -j, .true.)

                                    ! _L(i) > _RL(a) > ^RL^(jb)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%fullstop_L_to_R, &
                                                gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%R, &
                                                j, i, a, j, i, a, j, j, 0, 2, 1.0_dp, 1.0_dp)

                                else
                                    ! no restrictions
                                    call pgen_select_orb_guga_mol(csf_i, occ_orbs, &
                                                                  b, a, int_switch(2), cum_switch(2))

                                    ! only extremes count
                                    st = min(a, b)
                                    en = max(a, b)

                                    ! _L(i) > _RL(min) > ^LR(max) > ^R(j)
                                    excitInfo = assign_excitInfo_values_double( &
                                                excit_type%double_L_to_R, &
                                                gen_type%L, gen_type%R, gen_type%L, gen_type%L, gen_type%R, &
                                                en, i, st, j, i, st, en, j, 0, 4, 1.0_dp, 1.0_dp)
                                end if
                            end if
                        end if
                    end if
                end if
            end if
        end if

        ! if the holes are the same spatial orbital it makes no sense of
        ! considering the other order of picking..
        if (a == b) then
            pgen = pgen / 2.0_dp
        end if

        if (b == 1) then
            int_switch(1) = cum_arr(1)
        else
            int_switch(1) = cum_arr(b) - cum_arr(b - 1)
        end if
        cum_switch(1) = cum_sum(1)

        ! have to correctly adress if doubly occupied orbital was chosen for i,j
        ! and think about other stuff too!
        ! but in general it should look like:
        if (any(near_zero(cum_sum)) .or. any(near_zero(cum_switch))) then
            pgen = 0.0_dp
        else
            pgen = pgen * (product(int_contrib) / product(cum_sum) + &
                           product(int_switch) / product(cum_switch))
        end if

    end subroutine pickOrbs_sym_uniform_mol_double