gen_double_4ind_ex2 Subroutine

public subroutine gen_double_4ind_ex2(nI, ilutI, nJ, ilutJ, ex, par, pgen)


Type IntentOptional Attributes Name
integer, intent(in) :: nI(nel)
integer(kind=n_int), intent(in) :: ilutI(0:NIfTot)
integer, intent(out) :: nJ(nel)
integer(kind=n_int), intent(out) :: ilutJ(0:NIfTot)
integer, intent(out) :: ex(2,maxExcit)
logical, intent(out) :: par
real(kind=dp), intent(out) :: pgen


Source Code

Source Code

    subroutine gen_double_4ind_ex2(nI, ilutI, nJ, ilutJ, ex, par, pgen)

        use GenRandSymExcitNUMod, only: RandExcitSymLabelProd
        use SymExcitDataMod, only: SpinOrbSymLabel
        use lattice_models_utils, only: pick_spin_opp_elecs

        integer, intent(in) :: nI(nel)
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
        integer, intent(out) :: nJ(nel), ex(2, maxExcit)
        integer(n_int), intent(out) :: ilutJ(0:NIfTot)
        logical, intent(out) :: par
        real(dp), intent(out) :: pgen
        character(*), parameter :: this_routine = 'gen_double_4ind_ex2'

        integer :: elecs(2), src(2), orbs(2), ispn, sum_ml, cc_a, cc_b
        integer :: sym_product
        real(dp) :: int_cpt(2), cum_sum(2), cpt_pair(2), sum_pair(2)

        ! This array stores the cumulative list used in selecting the first
        ! orbital so that it can be used for calculating the reverse
        ! probability if required.
        real(dp) :: cum_arr(nbasis)

        real(dp) :: scratch_cpt, scratch_sm
        integer :: scratch_orb
        integer :: loc

        ! Pick the electrons in a weighted fashion
        if (t_mixed_hubbard .or. t_olle_hubbard) then
            call pick_biased_elecs(nI, elecs, src, sym_product, ispn, sum_ml, pgen)
            call pick_weighted_elecs(nI, elecs, src, sym_product, ispn, sum_ml, &
        end if

        ! then first pick (a) orbital:
        ! for opposite spin excitations (a) is restricted to be a beta orbital!
        ! and the probability is split p(a|ij) = p(j)*p(a|i)
        ! except in the symmetric excitation generator, which isn't used
        ! ever anyway..
        orbs(1) = pick_a_orb(ilutI, src, iSpn, int_cpt(1), cum_sum(1), cum_arr)

        ! Select the B orbital, in the same way as before!!
        ! The symmetry of this second orbital depends on that of the first.
        if (orbs(1) /= 0) then
            cc_a = ClassCountInd(orbs(1))
            cc_b = get_paired_cc_ind(cc_a, sym_product, sum_ml, iSpn)

            ! pick the last orbitals weighted with the exact matrix
            ! element
            orbs(2) = select_orb(ilutI, src, cc_b, orbs(1), int_cpt(2), &
        end if

        ! what does this assert do?  do i have to pick the electrons in a
        ! certain order??

        ASSERT((.not. (is_beta(orbs(2)) .and. .not. is_beta(orbs(1)))) .or. tGen_4ind_2_symmetric)
        if (any(orbs == 0)) then
            nJ(1) = 0
            pgen = 0.0_dp
        end if

        ! can i exit right away if this happens??
        ! i am pretty sure this means
        if (any(cum_sum < EPS)) then
            cum_sum = 1.0_dp
            int_cpt = 0.0_dp
        end if

        ! Calculate the pgens. Note that all of these excitations can be
        ! selected as both A--B or B--A. So these need to be calculated
        ! explicitly.
        if (.not. tGen_guga_crude) then
        end if
        if ((is_beta(orbs(1)) .eqv. is_beta(orbs(2))) .or. tGen_4ind_2_symmetric) then

            ! in the case of parallel spin excitations or symmetrice excitation
            ! generation(but does actually someone use that?) we have to
            ! calculate the probability of picking the holes the other
            ! way around
            call pgen_select_a_orb(ilutI, src, orbs(2), iSpn, cpt_pair(1), &
                                   sum_pair(1), cum_arr, .false.)
            call pgen_select_orb(ilutI, src, orbs(2), orbs(1), &
                                 cpt_pair(2), sum_pair(2))

            cpt_pair = 0.0_dp
            sum_Pair = 1.0_dp
        end if

        if (any(sum_pair < EPS)) then
            cpt_pair = 0.0_dp
            sum_pair = 1.0_dp
        end if

        pgen = pgen * (product(int_cpt) / product(cum_sum) + &
                       product(cpt_pair) / product(sum_pair))

        ! And generate the actual excitation
        call make_double(nI, nJ, elecs(1), elecs(2), orbs(1), orbs(2), &
                         ex, par)
        ilutJ = ilutI
        clr_orb(ilutJ, src(1))
        clr_orb(ilutJ, src(2))
        set_orb(ilutJ, orbs(1))
        set_orb(ilutJ, orbs(2))

    end subroutine gen_double_4ind_ex2