calcRaisingSemiStartStochastic Subroutine

public subroutine calcRaisingSemiStartStochastic(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 calcRaisingSemiStartStochastic(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 = "calcRaisingSemiStartStochastic"

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

        ASSERT(.not. isThree(ilut, excitInfo%secondStart))

        ! i can be sure that there is no 3 or 0 at the fullEnd, or otherwise
        ! this would be single-excitation like.
        se = excitInfo%secondStart
        bVal = csf_i%B_real(se)

        deltaB = getDeltaB(t)

        ! do non-choosing possibs first
        ! why does this cause a segfault on compilation with gfortran??
        ! do some debugging:

        ! fix for gfortran compilation for some reasono
        ! i can probably fix it when i finally get to this point in
        ! test running

        select case (csf_i%stepvector(se))
        case (1)
            ! 1 -> 3
            set_orb(t, 2 * se)

            call setDeltaB(deltaB + 1, t)

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

        case (2)
            ! 2 -> 3
            set_orb(t, 2 * se - 1)

            call setDeltaB(deltaB - 1, t)

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

        case (0)
            ! 0:
            ! for arriving -1 branch branching is always possible
            if (deltaB == -1) then
                ! here the choice is between 0 and -2 branch
                minusWeight = weights%proc%minus(negSwitches(se), bVal, weights%dat)
                zeroWeight = weights%proc%zero(negSwitches(se), posSwitches(se), &
                                               bVal, weights%dat)

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

                ! cant directly cant assign it to probWeight
                zeroWeight = calcStartProb(zeroWeight, minusWeight)

                if (genrand_real2_dSFMT() < zeroWeight) then
                    ! go to 0 branch
                    ! 0 -> 2
                    set_orb(t, 2 * se)

                    call setDeltaB(0, t)

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

                    probWeight = probWeight * zeroWeight

                else
                    ! to to -2 branch
                    ! 0 -> 1
                    set_orb(t, 2 * se - 1)

                    call setDeltaB(-2, t)

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

                    tempWeight_0 = 0.0_dp

                    probWeight = probWeight * (1.0_dp - zeroWeight)
                end if
            else
                ! +1 branch arriving -> have to check b values
                ! UPDATE: include b value check into probWeight calculation
                ! todo
                if (csf_i%B_int(se) < 2) then
                    ! only 0 branch possible
                    ! todo: in this forced cases due to the b value, have to
                    ! think about, how that influences the probWeight...
                    ! 0 -> 1
                    set_orb(t, 2 * se - 1)

                    call setDeltaB(0, t)

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

                else
                    ! need probs
                    plusWeight = weights%proc%plus(posSwitches(se), &
                                                   bVal, weights%dat)
                    zeroWeight = weights%proc%zero(negSwitches(se), posSwitches(se), &
                                                   bVal, weights%dat)

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

                    ! cant directly cant assign it to probWeight
                    zeroWeight = calcStartProb(zeroWeight, plusWeight)

                    if (genrand_real2_dSFMT() < zeroWeight) then
                        ! do 0 branch
                        ! 0 -> 1
                        set_orb(t, 2 * se - 1)

                        call setDeltaB(0, t)

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

                        probWeight = probWeight * zeroWeight

                    else
                        ! do +2 branch
                        ! 0 -> 2
                        set_orb(t, 2 * se)

                        call setDeltaB(2, t)

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

                        tempWeight_0 = 0.0_dp

                        probWeight = probWeight * (1.0_dp - zeroWeight)
                    end if
                end if
            end if
        end select

        call encode_matrix_element(t, extract_matrix_element(t, 1) * tempWeight_1, 2)
        call update_matrix_element(t, tempWeight_0, 1)

        if (near_zero(tempWeight_0) .and. near_zero(tempWeight_1)) then
            probWeight = 0.0_dp
            t = 0
        end if

    end subroutine calcRaisingSemiStartStochastic