calc_mixed_start_contr_pgen Subroutine

public subroutine calc_mixed_start_contr_pgen(ilut, csf_i, t, excitInfo, branch_pgen, pgen)

Arguments

Type IntentOptional Attributes Name
integer(kind=n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
integer(kind=n_int), intent(in) :: t(0:nifguga)
type(ExcitationInformation_t), value :: excitInfo
real(kind=dp), value :: branch_pgen
real(kind=dp), intent(out) :: pgen

Contents


Source Code

    subroutine calc_mixed_start_contr_pgen(ilut, csf_i, t, excitInfo, branch_pgen, &
            pgen)
        integer(n_int), intent(in) :: ilut(0:nifguga), t(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), value :: excitInfo
        real(dp), value :: branch_pgen
        real(dp), intent(out) :: pgen
        character(*), parameter :: this_routine = "calc_mixed_start_contr_pgen"
        integer :: st, se, en, elecInd, holeInd, sw, step, i
        real(dp) :: orb_pgen, zero_weight, switch_weight, stay_weight
        real(dp) :: bot_cont, posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
        real(dp) :: new_pgen, start_mat, stay_mat, start_weight
        type(WeightObj_t) :: weights
        logical :: below_flag

        st = excitInfo%fullStart
        se = excitInfo%firstEnd
        en = excitInfo%fullEnd
        ! depending on the type of excitaiton, calculation of orbital pgens
        ! change
        if (excitInfo%typ == excit_type%fullStart_L_to_R) then
            elecInd = en
            holeInd = se
        else if (excitInfo%typ == excit_type%fullstart_R_to_L) then
            elecInd = se
            holeInd = en
        else
            call stop_all(this_routine, "should not be here!")
        end if

        sw = findFirstSwitch(ilut, t, st, se)

        ! what can i precalculate beforehand?
        step = csf_i%stepvector(st)

        ! do i actually deal with the actual start orbital influence??
        ! fuck i don't think so.. wtf..
        call calc_orbital_pgen_contrib_start(csf_i, [2 * st, 2 * elecInd], &
            holeInd, orb_pgen)

        pgen = orb_pgen * branch_pgen

        ! since weights only depend on the number of switches at the
        ! semistop and semistop and full-end index i can calculate
        ! it beforehand for all?
        excitInfo%fullStart = 1
        excitInfo%secondStart = 1
        call calcRemainingSwitches_excitInfo_double(csf_i, excitInfo, posSwitches, negSwitches)

        weights = init_fullStartWeight(csf_i, se, en, negSwitches(se), &
                                       posSwitches(se), csf_i%B_real(se))

        ! determine the original starting weight
        zero_weight = weights%proc%zero(negSwitches(st), posSwitches(st), &
                                        csf_i%B_real(st), weights%dat)

        if (step == 1) then

            switch_weight = weights%proc%plus(posSwitches(st), csf_i%B_real(st), &
                                              weights%dat)

            if (isOne(t, st)) then

                bot_cont = Root2 * sqrt((csf_i%B_real(st) - 1.0_dp) / &
                                        (csf_i%B_real(st) + 1.0_dp))

                stay_weight = calcStayingProb(zero_weight, switch_weight, &
                                              csf_i%B_real(st))

                start_weight = zero_weight / (zero_weight + switch_weight)

            else

                bot_cont = -sqrt(2.0_dp / ((csf_i%B_real(st) - 1.0_dp) * &
                                           (csf_i%B_real(st) + 1.0_dp)))

                stay_weight = 1.0_dp - calcStayingProb(zero_weight, switch_weight, &
                                                       csf_i%B_real(st))

                start_weight = switch_weight / (zero_weight + switch_weight)

            end if
        else

            switch_weight = weights%proc%minus(negSwitches(st), &
                                               csf_i%B_real(st), weights%dat)

            if (isOne(t, st)) then
                bot_cont = -sqrt(2.0_dp / ((csf_i%B_real(st) + 1.0_dp) * &
                                           (csf_i%B_real(st) + 3.0_dp)))

                stay_weight = 1.0_dp - calcStayingProb(zero_weight, switch_weight, &
                                                       csf_i%B_real(st))

                start_weight = switch_weight / (zero_weight + switch_weight)

            else
                bot_cont = -Root2 * sqrt((csf_i%B_real(st) + 3.0_dp) / &
                                         (csf_i%B_real(st) + 1.0_dp))

                stay_weight = calcStayingProb(zero_weight, switch_weight, &
                                              csf_i%B_real(st))

                start_weight = zero_weight / (zero_weight + switch_weight)
            end if
        end if

        ASSERT(.not. near_zero(start_weight))
        ! update the pgen stumbs here to reuse start_weight variable
        new_pgen = stay_weight * branch_pgen / start_weight

        ! divide out the original starting weight:
        branch_pgen = branch_pgen / start_weight

        ! loop from start backwards so i can abort at a d=1 & b=1 stepvalue
        ! also consider if bot_cont < EPS to avoid unnecarry calculations
        if (.not. near_zero(bot_cont)) then

            below_flag = .false.

            do i = st - 1, 1, -1
                if (csf_i%Occ_int(i) /= 1) cycle

                ! then check if thats the last stepvalue to consider
                if (csf_i%stepvector(i) == 1 .and. csf_i%B_int(i) == 1) then
                    below_flag = .true.
                end if

                ! then i need to calculate the orbital probability
                ! from the fact that this is a lowering into raising fullstart
                ! i know, more about the restrictions...
                ! and that fullend is the electron eg.
                ! depening on the type of excitation (r2l or l2r) the electron
                ! orbitals change here
                call calc_orbital_pgen_contrib_start(csf_i, [2*i, 2*elecInd], &
                    holeInd, orb_pgen)

                ! then deal with the matrix element and branching probabilities
                step = csf_i%stepvector(i)

                ! get both start and staying matrix elements -> and update
                ! matrix element contributions on the fly to avoid second loop!
                call getDoubleMatrixElement(step, step, -1, gen_type%R, gen_type%L, &
                    csf_i%B_real(i), 1.0_dp, x1_element=start_mat)

                call getDoubleMatrixElement(step, step, 0, gen_type%R, gen_type%L, &
                    csf_i%B_real(i), 1.0_dp, x1_element=stay_mat)

                ! another check.. although this should not happen
                ! except the other d = 1 & b = 1 condition is already met
                ! above, to not continue:
                if (near_zero(stay_mat)) below_flag = .true.

                zero_weight = weights%proc%zero(negSwitches(i), posSwitches(i), &
                    csf_i%B_real(i), weights%dat)

                if (step == 1) then
                    switch_weight = weights%proc%plus(posSwitches(i), &
                                                      csf_i%B_real(i), weights%dat)
                else
                    switch_weight = weights%proc%minus(negSwitches(i), &
                                                       csf_i%B_real(i), weights%dat)
                end if

                start_weight = zero_weight / (zero_weight + switch_weight)
                stay_weight = calcStayingProb(zero_weight, switch_weight, &
                                              csf_i%B_real(i))

                ! i think i could avoid the second loop over j
                ! if i express everything in terms of already calculated
                ! quantities!
                ! "normally" matrix element shouldnt be 0 anymore... still check
                if (.not. near_zero(start_mat)) then
                    pgen = pgen + orb_pgen * start_weight * new_pgen
                end if

                if (below_flag) exit

                if (t_trunc_guga_pgen .or. &
                    (t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
                    if (new_pgen < trunc_guga_pgen) then
                        new_pgen = 0.0_dp
                    end if
                end if

                ! update new_pgen for next cycle
                new_pgen = stay_weight * new_pgen

            end do

        end if

        ! start to switch loop: here matrix elements are not 0!
        ! and its only db = 0 branch and no stepvalue change!
        ! if the start is the switch nothing happens

        step = csf_i%stepvector(st)

        ! calculate the necarry values needed to formulate everything in terms
        ! of the already calculated quantities:
        call getDoubleMatrixElement(step, step, -1, gen_type%L, gen_type%R, &
            csf_i%B_real(st), 1.0_dp, x1_element=start_mat)

        if (.not. near_zero(abs(start_mat))) then

            do i = st + 1, sw - 1
                ! the good thing here is, i do not need to loop a second time,
                ! since i can recalc. the matrix elements and pgens on-the fly
                ! here the matrix elements should not be 0 or otherwise the
                ! excitation wouldnt have happended anyways
                if (csf_i%Occ_int(i) /= 1) cycle

                ! calculate orbitals pgen first and cycle if 0
                call calc_orbital_pgen_contrib_start(csf_i, [2 * i, 2 * elecInd], &
                    holeInd, orb_pgen)

                step = csf_i%stepvector(i)

                ! update inverse product
#ifdef DEBUG_
                call getDoubleMatrixElement(step, step, 0, gen_type%L, gen_type%R, &
                    csf_i%B_real(i), 1.0_dp, x1_element=stay_mat)
                ASSERT(.not. near_zero(stay_mat))
#endif

                ! and update pgens also
                zero_weight = weights%proc%zero(negSwitches(i), posSwitches(i), &
                    csf_i%B_real(i), weights%dat)

                if (step == 1) then
                    switch_weight = weights%proc%plus(posSwitches(i), &
                                                csf_i%B_real(i), weights%dat)
                else
                    switch_weight = weights%proc%minus(negSwitches(i), &
                                                csf_i%B_real(i), weights%dat)
                end if

                stay_weight = calcStayingProb(zero_weight, switch_weight, &
                                              csf_i%B_real(i))

                start_weight = zero_weight / (zero_weight + switch_weight)

                branch_pgen = branch_pgen / stay_weight

                if (t_trunc_guga_pgen .or. &
                    (t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then

                    if (branch_pgen * start_weight > trunc_guga_pgen) then
                        pgen = pgen + orb_pgen * branch_pgen * start_weight
                    end if
                else
                    ! and add up correctly
                    pgen = pgen + orb_pgen * branch_pgen * start_weight
                end if

            end do

            ! handle switch seperately (but only if switch > start)
            if (sw > st) then

                ! check orb_pgen otherwise no influencce
                call calc_orbital_pgen_contrib_start(csf_i, [2 * sw, 2 * elecInd], &
                    holeInd, orb_pgen)

                if (.not. near_zero(orb_pgen)) then

                    step = csf_i%stepvector(sw)

                    zero_weight = weights%proc%zero(negSwitches(sw), posSwitches(sw), &
                                                    csf_i%B_real(sw), weights%dat)

                    ! on the switch the original probability is:
                    if (step == 1) then
                        switch_weight = weights%proc%plus(posSwitches(sw), &
                                                csf_i%B_real(sw), weights%dat)

#ifdef DEBUG_
                        call getDoubleMatrixElement(2, 1, 0, gen_type%L, gen_type%R, &
                            csf_i%B_real(sw), 1.0_dp, x1_element=stay_mat)
#endif

                    else
                        switch_weight = weights%proc%minus(negSwitches(sw), &
                                                csf_i%B_real(sw), weights%dat)

#ifdef DEBUG_
                        call getDoubleMatrixElement(1, 2, 0, gen_type%L, gen_type%R, &
                            csf_i%B_real(sw), 1.0_dp, x1_element=stay_mat)
#endif

                    end if

                    ! update inverse product
                    ! and also get starting contribution
                    ASSERT(.not. near_zero(stay_mat))

                    stay_weight = 1.0_dp - calcStayingProb(zero_weight, switch_weight, &
                                                           csf_i%B_real(sw))

                    ! and the new startProb is also the non-b=0 branch
                    start_weight = switch_weight / (zero_weight + switch_weight)

                    if (t_trunc_guga_pgen .or. &
                        (t_trunc_guga_pgen_noninits .and. .not. is_init_guga)) then
                        if (branch_pgen * start_weight / start_weight > trunc_guga_pgen) then
                            pgen = pgen + orb_pgen * branch_pgen * start_weight / stay_weight
                        end if

                    else
                        pgen = pgen + orb_pgen * branch_pgen * start_weight / stay_weight
                    end if

                end if
            end if
        end if

        ! i also need to consider the electron pair picking probability..
        pgen = pgen / real(ElecPairs, dp)
        ! and if the second electron is in a double occupied orbital I have
        ! to modify it with 2
        if (csf_i%stepvector(elecInd) == 3) pgen = pgen * 2.0_dp

    end subroutine calc_mixed_start_contr_pgen