calcFullstopLowering Subroutine

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

        ! not sure if single or double weight is necessary here...
        type(WeightObj_t) :: weights
        integer(n_int), allocatable :: tempExcits(:, :)
        integer(n_int) :: t(0:nifguga)
        integer :: iOrb, iEx, deltaB, cnt, ierr
        real(dp) :: tempWeight, sig, bVal

        ! create weight object here
        ! i think i only need single excitations weights here, since
        ! the semi stop in this case is like an end step...
        weights = init_singleWeight(csf_i, excitInfo%secondStart)

        ! create start
        call createSingleStart(ilut, csf_i, excitInfo, posSwitches, negSwitches, &
                               weights, tempExcits, nExcits)

        excitInfo%currentGen = excitInfo%firstGen
        ! loop until semi-start
        do iOrb = excitInfo%fullStart + 1, excitInfo%secondStart - 1
            call singleUpdate(ilut, csf_i, iOrb, excitInfo, posSwitches, negSwitches, &
                              weights, tempExcits, nExcits)
        end do

        ! at semi start with alike generator full stops only the deltaB=0
        ! are compatible, i hope this is correctly accounted by the weights..
        iOrb = excitInfo%secondStart
        bVal = csf_i%B_real(iOrb)

        select case (csf_i%stepvector(iOrb))
        case (1)
            ! only delta -1 branches can lead to deltaB=0 in overlap here
            do iEx = 1, nExcits
                ! correct use of weight gens should exclude this, but check:

                ASSERT(getDeltaB(tempExcits(:, iEx)) == -1)

                t = tempExcits(:, iEx)

                ! change 1 -> 0
                clr_orb(t, 2 * iOrb - 1)

                call getDoubleMatrixElement(0, 1, -1, gen_type%L, gen_type%L, &
                                            bVal, excitInfo%order, tempWeight)

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

                call setDeltaB(0, t)

                tempExcits(:, iEx) = t
            end do

        case (2)
            ! only deltaB=+1 branches lead to a 0 branch
            do iEx = 1, nExcits
                ASSERT(getDeltaB(tempExcits(:, iEx)) == +1)

                t = tempExcits(:, iEx)

                ! change 2 -> 0
                clr_orb(t, 2 * iOrb)

                call getDoubleMatrixElement(0, 2, +1, gen_type%L, gen_type%L, &
                                            bVal, 1.0_dp, tempWeight)

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

                call setDeltaB(0, t)

                tempExcits(:, iEx) = t
            end do

        case (3)
            ! has to be a 3 in the lowering case.
            ASSERT(isThree(ilut, iOrb))
            ! -1 branches always go to 0 branch
            ! +1 branches might go to 0 branch if bValue high enough
            ! do not have to check weights here, since 0 branch has to have
            ! non-zero weight or it wouldnt even be compatible
            ! no! bullshit! switchin to 0 branch always works!
            ! both branches survive
            do iEx = 1, nExcits
                t = tempExcits(:, iEx)

                deltaB = getDeltaB(t)

                if (deltaB == -1) then

                    ! change 3->2
                    clr_orb(t, 2 * iOrb - 1)

                    call getDoubleMatrixElement(2, 3, deltaB, gen_type%L, gen_type%L, &
                                                bVal, 1.0_dp, tempWeight)

                else
                    ! change 3->1
                    clr_orb(t, 2 * iOrb)

                    call getDoubleMatrixElement(1, 3, deltaB, gen_type%L, gen_type%L, &
                                                bVal, 1.0_dp, tempWeight)

                end if

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

                call setDeltaB(0, t)

                tempExcits(:, iEx) = t

            end do
        end select

        ! continue on with double excitation region, only the 0 branch
        ! valid here, where there is no change in stepvector and matrix
        ! element only a sign dependent on the number of open orbitals
        sig = (-1.0_dp)**real(count_open_orbs_ij(csf_i, excitInfo%secondStart + 1, &
                                                 excitInfo%fullEnd - 1), dp)

        ! do the ending
        ASSERT(isZero(ilut, excitInfo%fullEnd))

        iOrb = excitInfo%fullEnd

        cnt = 1

        do iEx = 1, nExcits
            t = tempExcits(:, iEx)

            if (near_zero(extract_matrix_element(t, 1))) cycle
            ! change 0 -> 3
            set_orb(t, 2 * iOrb)
            set_orb(t, 2 * iOrb - 1)

            ! have to cancel 0 weighted excitations

            ! also include the sign of the open orbitals here
            call update_matrix_element(t, sig * Root2, 1)

            tempExcits(:, cnt) = t
            cnt = cnt + 1

        end do

        nExcits = cnt - 1

        ! have to put it in excitations stupid!
        allocate(excitations(0:nifguga, nExcits), stat=ierr)
        excitations = tempExcits(:, 1:nExcits)
        deallocate(tempExcits)

    end subroutine calcFullstopLowering