calcRaisingSemiStopStochastic Subroutine

public subroutine calcRaisingSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, posSwitches, t, probWeight)

Arguments

Type IntentOptional Attributes Name
integer(kind=n_int), intent(in) :: ilut(0:nifguga)
type(CSF_Info_t), intent(in) :: csf_i
type(ExcitationInformation_t), intent(in) :: excitInfo
type(WeightObj_t), intent(in) :: weights
real(kind=dp), intent(in) :: negSwitches(nSpatOrbs)
real(kind=dp), intent(in) :: posSwitches(nSpatOrbs)
integer(kind=n_int), intent(inout) :: t(0:nifguga)
real(kind=dp), intent(inout) :: probWeight

Contents


Source Code

    subroutine calcRaisingSemiStopStochastic(ilut, csf_i, excitInfo, weights, negSwitches, &
                                             posSwitches, t, probWeight)
        integer(n_int), intent(in) :: ilut(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        type(WeightObj_t), intent(in) :: weights
        real(dp), intent(in) :: negSwitches(nSpatOrbs), posSwitches(nSpatOrbs)
        integer(n_int), intent(inout) :: t(0:nifguga)
        real(dp), intent(inout) :: probWeight
        character(*), parameter :: this_routine = "calcRaisingSemiStopStochastic"

        integer :: semi, deltaB
        real(dp) :: tempWeight_0, tempWeight_1, minusWeight, &
                    plusWeight, bVal

        ASSERT(.not. isZero(ilut, excitInfo%firstEnd))

        semi = excitInfo%firstEnd
        bVal = csf_i%B_real(semi)

        ! in the stochastic case i am sure that at there is no 3 or 0 at the
        ! full start... so definetly all deltaB branches can arrive here
        ! first deal with the forced ones
        deltaB = getDeltaB(t)

        select case (csf_i%stepvector(semi))
        case (1)
            if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
                if (getDeltaB(t) == 2) then
                    t = 0_n_int
                    probWeight = 0.0_dp
                    return
                end if
            end if

            ASSERT(getDeltaB(t) /= 2)

            ! do the 1 -> 0 switch
            clr_orb(t, 2 * semi - 1)

            call setDeltaB(deltaB + 1, t)

            call getDoubleMatrixElement(0, 1, deltaB, excitInfo%gen1, &
                                        excitInfo%gen2, bVal, excitInfo%order1, &
                                        tempWeight_0, tempWeight_1)

        case (2)
            if (t_approx_exchange .or. (t_approx_exchange_noninits .and. (.not. is_init_guga))) then
                if (getDeltaB(t) == -2) then
                    t = 0_n_int
                    probWeight = 0.0_dp
                    return
                end if
            end if

            ASSERT(getDeltaB(t) /= -2)

            ! do 2 -> 0
            clr_orb(t, 2 * semi)

            call setDeltaB(deltaB - 1, t)

            call getDoubleMatrixElement(0, 2, deltaB, excitInfo%gen1, &
                                        excitInfo%gen2, bVal, excitInfo%order1, &
                                        tempWeight_0, tempWeight_1)

        case (3)
            ! its a 3 -> have to check b if a 0 branch arrives
            select case (deltaB)
            case (-2)
                ! do 3 -> 2
                clr_orb(t, 2 * semi - 1)

                call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
                                            excitInfo%gen2, bVal, excitInfo%order1, &
                                            x1_element=tempWeight_1)

                tempWeight_0 = 0.0_dp

                call setDeltaB(-1, t)

            case (2)
                ! do 3 -> 1
                clr_orb(t, 2 * semi)

                call setDeltaB(+1, t)

                call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
                                            excitInfo%gen2, bVal, excitInfo%order1, &
                                            x1_element=tempWeight_1)

                tempWeight_0 = 0.0_dp

            case (0)
                ! deltaB = 0 branch arrives -> check b
                if (csf_i%B_int(semi) == 0) then
                    ! only 3 -> 1 possble
                    clr_orb(t, 2 * semi)

                    call setDeltaB(-1, t)

                    call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
                                                excitInfo%gen2, bVal, excitInfo%order1, &
                                                tempWeight_0, tempWeight_1)

                else
                    ! finally have to check the probs
                    minusWeight = weights%proc%minus(negSwitches(semi), &
                                                     bVal, weights%dat)
                    plusWeight = weights%proc%plus(posSwitches(semi), &
                                                   bVal, weights%dat)

                    if (near_zero(minusWeight + plusWeight)) then
                        probWeight = 0.0_dp
                        t = 0
                        return
                    end if

                    ! here i cant directly save it to probWeight ...
                    minusWeight = calcStartProb(minusWeight, plusWeight)

                    if (genrand_real2_dSFMT() < minusWeight) then
                        ! do -1 branch:
                        ! set 3 -> 1
                        clr_orb(t, 2 * semi)

                        call setDeltaB(-1, t)

                        call getDoubleMatrixElement(1, 3, deltaB, excitInfo%gen1, &
                                                    excitInfo%gen2, bVal, excitInfo%order1, &
                                                    tempWeight_0, tempWeight_1)

                        probWeight = probWeight * minusWeight

                    else
                        ! do +1 branch
                        ! set 3 -> 2
                        clr_orb(t, 2 * semi - 1)

                        call setDeltaB(1, t)

                        call getDoubleMatrixElement(2, 3, deltaB, excitInfo%gen1, &
                                                    excitInfo%gen2, bVal, excitInfo%order1, &
                                                    tempWeight_0, tempWeight_1)

                        probWeight = probWeight * (1 - minusWeight)
                    end if
                end if
            end select
        end select

        ! if its a mixed full-start raising into lowering -> check i a
        ! switch happened in the double overlap region
        if (excitInfo%typ == excit_type%fullstart_R_to_L) then
            ! this is indicated by a non-zero x0-matrix element
            if (.not. near_zero(extract_matrix_element(t, 1) * tempWeight_0)) then
                probWeight = 0.0_dp
                t = 0
                return
            end if
        end if

        ! do not combine the x0 and x1 elements at this point, because its
        ! easier to deal with the contributing integrals if we keep them
        ! seperate until the end!

        if (near_zero(extract_matrix_element(t, 1) * tempWeight_0) .and. &
            near_zero(extract_matrix_element(t, 2) * tempWeight_1)) then
            probWeight = 0.0_dp
            t = 0_n_int
            return
        end if

        call update_matrix_element(t, tempWeight_0, 1)
        call update_matrix_element(t, tempWeight_1, 2)

    end subroutine calcRaisingSemiStopStochastic