calcFullStartRaising Subroutine

public subroutine calcFullStartRaising(ilut, csf_i, excitInfo, excitations, nExcits, posSwitches, negSwitches)

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(inout) :: excitInfo
integer(kind=n_int), intent(out), allocatable :: excitations(:,:)
integer, intent(out) :: nExcits
real(kind=dp), intent(in) :: posSwitches(nSpatOrbs)
real(kind=dp), intent(in) :: negSwitches(nSpatOrbs)

Contents

Source Code


Source Code

    subroutine calcFullStartRaising(ilut, csf_i, excitInfo, excitations, nExcits, &
                                    posSwitches, negSwitches)
        integer(n_int), intent(in) :: ilut(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(inout) :: excitInfo
        integer(n_int), intent(out), allocatable :: excitations(:, :)
        integer, intent(out) :: nExcits
        real(dp), intent(in) :: posSwitches(nSpatOrbs), negSwitches(nSpatOrbs)
        character(*), parameter :: this_routine = "calcFullStartRaising"

        integer :: nMax, ierr, iOrb, start, ende, semi, gen
        integer(n_int) :: t(0:nifguga)
        real(dp) :: tempWeight, minusWeight, plusWeight, bVal, nOpen
        type(WeightObj_t) :: weights
        integer(n_int), allocatable :: tempExcits(:, :)

        ASSERT(isZero(ilut, excitInfo%fullStart))
        ASSERT(isProperCSF_ilut(ilut))

        ! create the fullStart
        start = excitInfo%fullStart
        ende = excitInfo%fullEnd
        semi = excitInfo%firstEnd
        gen = excitInfo%lastGen
        bVal = csf_i%B_real(semi)
        ! first have to allocate both arrays for the determinant and weight list
        ! worse then worst case for single excitations are 2^|i-j| excitations
        ! for a given CSF -> for now use that to allocate the arrays first
        ! update only need the number of open orbitals between i and j, and
        ! some additional room if 0/3 at start
        ! use already provided open orbital counting function.
        ! nMax = 2**(ende - start)
        nMax = 2 + 2**count_open_orbs_ij(csf_i, start, ende)
        allocate(tempExcits(0:nifguga, nMax), stat=ierr)

        t = ilut

        ! have to make the switch and store the first matrix element in it too
        ! additionally also already calculate the sign coming from the
        ! pseudo double excitation which only depends on the number of open
        ! orbitals in the overlap region
        nOpen = real(count_open_orbs_ij(csf_i, start, semi - 1), dp)

        ! set 0->3
        set_orb(t, 2 * start)
        set_orb(t, 2 * start - 1)

        ! could encode matrix element here, but do it later for more efficiency
        ! call encode_part_sign(t,Root2 * (-1.0)**nMax, 1)

        ! also dont need to store it here, also later
        ! tempExcits(:,1) = t

        ! since only the 0 branch is taken, where there are now changes in
        ! the stepvector and the matrix element is just a sign i can directly
        ! go to the lowering semi stop
        ! i could write it specifically here, and i should for efficiency

        ! have to calc. weights here, which are only the normel single
        ! excitation weights
        weights = init_singleWeight(csf_i, ende)
        minusWeight = weights%proc%minus(negSwitches(semi), csf_i%B_real(semi), weights%dat)
        plusWeight = weights%proc%plus(posSwitches(semi), csf_i%B_real(semi), weights%dat)

        ASSERT(.not. isZero(ilut, semi))
        ASSERT(minusWeight + plusWeight > 0.0_dp)

        call encode_matrix_element(t, 0.0_dp, 2)

        select case (csf_i%stepvector(semi))
        case (3)
            ! have to check a few things with 3 semi stop

            ! did something wrong with matrix elements before here.. the 3
            ! semi-stop has an aditional sign...but all matrix elements
            ! the same..
            call encode_matrix_element(t, (-1.0_dp)**(nOpen + 1.0_dp), 1)

            ! bullshit...
            if (minusWeight > 0.0_dp .and. plusWeight > 0.0_dp) then
                ! both +1 and -1 excitations possible
                ! do 3->1 first -1 branch first
                clr_orb(t, 2 * semi)

                call setDeltaB(-1, t)

                tempExcits(:, 1) = t

                ! then do 3->2: +1 branch(change from already set 1
                clr_orb(t, 2 * semi - 1)
                set_orb(t, 2 * semi)

                call setDeltaB(1, t)

                tempExcits(:, 2) = t

                nExcits = 2

            else if (near_zero(plusWeight)) then
                ! only -1 branch possible
                ! do 3->1 first -1 branch first
                clr_orb(t, 2 * semi)

                call setDeltaB(-1, t)

                tempExcits(:, 1) = t

                nExcits = 1
            else if (near_zero(minusWeight)) then
                ! only +1 branch possible
                ! then do 3->2: +1 branch
                clr_orb(t, 2 * semi - 1)

                call getDoubleMatrixElement(2, 3, 0, gen_type%R, gen_type%R, bVal, 1.0_dp, tempWeight)

                call encode_matrix_element(t, Root2 * tempWeight * ((-1.0_dp)**nOpen), 1)

                call setDeltaB(1, t)

                tempExcits(:, 1) = t

                nExcits = 1

            else
                ! something went wrong, or i have to cancel exciation.. todo how
                call stop_all(this_routine, "something went wrong. shouldnt be here!")

            end if

        case (1)
            ! only one excitation possible, which also has to have
            ! non-zero weight or otherwise i wouldnt even be here
            ! and reuse already provided t
            ! 1 -> 0
            clr_orb(t, 2 * semi - 1)

            call getDoubleMatrixElement(0, 1, 0, gen_type%R, gen_type%R, bVal, 1.0_dp, tempWeight)

            call setDeltaB(1, t)

            ! encode fullstart contribution and pseudo overlap region here
            ! too in one go. one additional -1 due to semistop

            ! no Root2 here, since it cancels with t from matEle,
            call encode_matrix_element(t, Root2 * tempWeight * ((-1.0_dp)**(nOpen)), 1)

            tempExcits(:, 1) = t

            nExcits = 1

        case (2)
            ! here we have
            ! 2 -> 0
            clr_orb(t, 2 * semi)

            call getDoubleMatrixElement(0, 2, 0, gen_type%R, gen_type%R, bVal, 1.0_dp, tempWeight)

            call setDeltaB(-1, t)

            ! encode fullstart contribution and pseudo overlap region here
            ! too in one go. one additional -1 due to semistop

            ! no Root2 here, since it cancels with t from matEle,
            call encode_matrix_element(t, Root2 * tempWeight * ((-1.0_dp)**(nOpen)), 1)

            tempExcits(:, 1) = t

            nExcits = 1

        end select

        ! and then we have to do just a regular single excitation
        excitInfo%currentGen = excitInfo%lastGen

        do iOrb = semi + 1, ende - 1
            call singleUpdate(ilut, csf_i, iOrb, excitInfo, posSwitches, negSwitches, &
                              weights, tempExcits, nExcits)
        end do

        ! and do a regular end step
        call singleEnd(ilut, csf_i, excitInfo, tempExcits, nExcits, excitations)

        ! that should be it...

    end subroutine calcFullStartRaising