calc_fullstart_alike_ex Subroutine

    pure subroutine calc_fullstart_alike_ex(csf_i, csf_j, excitInfo, mat_ele, &
                                       t_hamil, rdm_ind, rdm_mat)
        type(CSF_Info_t), intent(in) :: csf_i, csf_j
        type(ExcitationInformation_t), intent(in) :: excitInfo
        HElement_t(dp), intent(out) :: mat_ele
        logical, intent(in), optional :: t_hamil
        integer(int_rdm), intent(out), allocatable, optional :: rdm_ind(:)
        real(dp), intent(out), allocatable, optional :: rdm_mat(:)
        character(*), parameter :: this_routine = "calc_fullstart_alike_ex"

        integer :: i, step1, step2, db, delta_b(nSpatOrbs)
        real(dp) :: bVal, nOpen, temp_mat, guga_mat
        HElement_t(dp) :: umat
        logical :: t_hamil_

        def_default(t_hamil_, t_hamil, .true.)
        ASSERT(present(rdm_ind) .eqv. present(rdm_mat))

        ! set defaults for early exits
        mat_ele = h_cast(0.0_dp)
        delta_b = csf_i%B_int - csf_j%B_int

        associate (ii => excitInfo%i, jj => excitInfo%j, kk => excitInfo%k, &
                   ll => excitInfo%l, start => excitInfo%fullstart, &
                   ende => excitInfo%fullEnd, semi => excitInfo%firstEnd, &
                   gen => excitInfo%firstGen, typ => excitInfo%typ)

            if (present(rdm_ind) .and. present(rdm_mat)) then
                allocate(rdm_ind(1), source=0_int_rdm)
                allocate(rdm_mat(1), source=0.0_dp)
                rdm_ind(1) = contract_2_rdm_ind(ii, jj, kk, ll)
            end if

            ! i think i can exclude every deltaB > 1 sinve only dB = 0 branch
            ! allowed for the alike..
            if (any(abs(delta_b) > 1)) return

            if (t_hamil_) then
                if (typ == excit_type%fullstart_lowering) then
                    ! LL
                    umat = (get_umat_el(ende, semi, start, start) + &
                            get_umat_el(semi, ende, start, start)) / 2.0_dp

                else if (typ == excit_type%fullstart_raising) then
                    ! RR
                    umat = (get_umat_el(start, start, semi, ende) + &
                            get_umat_el(start, start, ende, semi)) / 2.0_dp
                end if
            else
                umat = h_cast(1.0_dp)
            end if

            if (t_hamil_ .and. near_zero(umat)) return

            nOpen = real(count_open_orbs_ij(csf_i, start, semi - 1), dp)

            ! do semi-stop
            step1 = csf_i%stepvector(semi)
            step2 = csf_j%stepvector(semi)
            db = delta_b(semi - 1)
            bVal = csf_i%b_real(semi)

            call getDoubleMatrixElement(step2, step1, db, gen, gen, bVal, 1.0_dp, temp_mat)

            if (near_zero(temp_mat)) return

            guga_mat = Root2 * temp_mat * (-1.0_dp)**nOpen

            ! do single range
            do i = semi + 1, ende

                step1 = csf_i%stepvector(i)
                step2 = csf_j%stepvector(i)
                db = delta_b(i - 1)
                bVal = csf_i%b_real(i)

                guga_mat = guga_mat * getSingleMatrixElement(step2, step1, db, gen, bVal)

                if (near_zero(guga_mat)) return

            end do

            mat_ele = guga_mat * umat

            ! also no influence on coupling coefficient sign from generator order
            if (present(rdm_mat)) rdm_mat = guga_mat

        end associate

    end subroutine calc_fullstart_alike_ex