mixedFullStop Subroutine

private subroutine mixedFullStop(ilut, csf_i, excitInfo, tempExcits, nExcits, excitations, t_no_singles_opt)

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
integer(kind=n_int), intent(inout), allocatable :: tempExcits(:,:)
integer, intent(inout) :: nExcits
integer(kind=n_int), intent(out), allocatable :: excitations(:,:)
logical, intent(in), optional :: t_no_singles_opt

Contents

Source Code


Source Code

    subroutine mixedFullStop(ilut, csf_i, excitInfo, tempExcits, nExcits, excitations, &
                             t_no_singles_opt)
        ! full stop routine for mixed generators
        integer(n_int), intent(in) :: ilut(0:nifguga)
        type(CSF_Info_t), intent(in) :: csf_i
        type(ExcitationInformation_t), intent(in) :: excitInfo
        integer(n_int), intent(inout), allocatable :: tempExcits(:, :)
        integer, intent(inout) :: nExcits
        integer(n_int), intent(out), allocatable :: excitations(:, :)
        logical, intent(in), optional :: t_no_singles_opt
        character(*), parameter :: this_routine = "mixedFullStop"

        integer :: iEx, ende, cnt, ierr, deltaB
        real(dp) :: bVal, tempWeight_0, tempWeight_1
        integer(n_int) :: t(0:nifguga)
        logical :: t_no_singles

        ASSERT(.not. isZero(ilut, excitInfo%fullEnd))
        ! not sure if i should check if ilut is not 3... since i could handle
        ! that differently

        if (present(t_no_singles_opt)) then
            t_no_singles = t_no_singles_opt
        else
            t_no_singles = .false.
        end if

        ende = excitInfo%fullEnd
        bVal = csf_i%B_real(ende)

        select case (csf_i%stepvector(ende))
        case (1)
            ! ending for -2 and 0 branches
            do iEx = 1, nExcits

                t = tempExcits(:, iEx)
                deltaB = getDeltaB(t)

                if (t_no_singles) then
                    ! check if the x0 is greater than 0, which indicates
                    ! purely delta-b = 0 route.., which should be disregarded
                    ! if we do not want to take singles into account.
                    if (.not. near_zero(extract_matrix_element(t, 1))) then
                        ! just to be sure this should also mean:
                        ASSERT(DeltaB == 0)

                        ! encode 0 matrix element, which then the loop at
                        ! the bottom should take care of and throw it out
                        call encode_matrix_element(t, 0.0_dp, 1)
                        call encode_matrix_element(t, 0.0_dp, 2)

                        tempExcits(:, iEx) = t

                        cycle

                    end if
                end if

                ! a +2 branch shouldnt arrive here
                ASSERT(deltaB /= 2)

                if (deltaB == 0) then
                    ! no change in stepvector so only calc matrix element
                    call getMixedFullStop(1, 1, 0, bVal, tempWeight_0, tempWeight_1)

                    ! finalize added up matrix elements in final step
                    tempWeight_0 = extract_matrix_element(t, 1) * tempWeight_0 + &
                                   extract_matrix_element(t, 2) * tempWeight_1

                else
                    ! -2 branch: change 1->2
                    set_orb(t, 2 * ende)
                    clr_orb(t, 2 * ende - 1)

                    ! matrix element in this case: only the x1 element is 1,
                    ! so just take the x1_element from the ilut
                    tempWeight_0 = extract_matrix_element(t, 2)

                end if
                ! and encode that into real part
                call encode_matrix_element(t, tempWeight_0, 1)
                call encode_matrix_element(t, 0.0_dp, 2)

                tempExcits(:, iEx) = t

            end do

        case (2)
            ! ending for 0 and +2 branches
            do iEx = 1, nExcits
                t = tempExcits(:, iEx)
                deltaB = getDeltaB(t)

                if (t_no_singles) then
                    ! check if the x0 is greater than 0, which indicates
                    ! purely delta-b = 0 route.., which should be disregarded
                    ! if we do not want to take singles into account.
                    if (.not. near_zero(extract_matrix_element(t, 1))) then
                        ! just to be sure this should also mean:
                        ASSERT(DeltaB == 0)

                        ! encode 0 matrix element, which then the loop at
                        ! the bottom should take care of and throw it out
                        call encode_matrix_element(t, 0.0_dp, 1)
                        call encode_matrix_element(t, 0.0_dp, 2)

                        tempExcits(:, iEx) = t

                        cycle

                    end if
                end if

                ASSERT(deltaB /= -2)

                if (deltaB == 0) then
                    call getMixedFullStop(2, 2, 0, bVal, tempWeight_0, tempWeight_1)

                    tempWeight_0 = extract_matrix_element(t, 1) * tempWeight_0 + &
                                   extract_matrix_element(t, 2) * tempWeight_1

                else
                    ! set 2-> 1 for +2 branch
                    set_orb(t, 2 * ende - 1)
                    clr_orb(t, 2 * ende)

                    tempWeight_0 = extract_matrix_element(t, 2)

                end if
                call encode_matrix_element(t, tempWeight_0, 1)
                call encode_matrix_element(t, 0.0_dp, 2)

                tempExcits(:, iEx) = t

            end do

        case (3)
            ! not sure if i ever access this function with 3 at end but
            ! nevertheless implement it for now...
            ! here only the 0 branches arrive, and whatever happend to switch
            ! at some point in the ovelap region has 0 matrix element
            do iEx = 1, nExcits
                t = tempExcits(:, iEx)
                deltaB = getDeltaB(t)

                ASSERT(deltaB == 0)

                call update_matrix_element(t, Root2, 1)
                call encode_matrix_element(t, 0.0_dp, 2)

                tempExcits(:, iEx) = t
            end do
        end select

        ! check again if there are maybe zero matrix elements
        cnt = 1
        do iEx = 1, nExcits
            if (near_zero(extract_matrix_element(tempExcits(:, iEx), 1))) cycle

            if (t_mixed_hubbard) then
                ! for mixed hubbard (and maybe other lattice systems,
                ! i do not want single excitations here (i guess.)
                if (.not. near_zero(extract_matrix_element(tempExcits(:, iEx), 1))) cycle
            end if

            tempExcits(:, cnt) = tempExcits(:, iEx)

            cnt = cnt + 1
        end do

        nExcits = cnt - 1

        ! do finishing up stuff..
        allocate(excitations(0:nifguga, nExcits), stat=ierr)

        excitations = tempExcits(:, 1:nExcits)

        deallocate(tempExcits)

    end subroutine mixedFullStop