create_cum_list_heisenberg Subroutine

    subroutine create_cum_list_heisenberg(ilutI, src, neighbors, cum_arr, cum_sum, &
                                          tgt, cpt)
        integer(n_int), intent(in) :: ilutI(0:NIfTot)
        integer, intent(in) :: src, neighbors(:)
        real(dp), intent(out), allocatable :: cum_arr(:)
        real(dp), intent(out) :: cum_sum
        integer, intent(in), optional :: tgt
        real(dp), intent(out), optional :: cpt
#ifdef DEBUG_
        character(*), parameter :: this_routine = "create_cum_list_heisenberg"
#endif
        integer :: flip, i, temp_ex(2, maxExcit), nI(nel)
        real(dp) :: elem

        ASSERT(IsOcc(ilutI, src))

        allocate(cum_arr(size(neighbors)))
        cum_arr = 0.0_dp
        cum_sum = 0.0_dp

        ! for the heisenberg model, where we know that every site is singly
        ! occupied, we search for empty spin-orbital neighbors, to see if
        ! a spin-flip is possible! so we do not need to check for the type
        ! of spin of orbital src here!
        ! although for the matrix element i need the opposite spin!
        ! add the according flip to get the other spin!

        ! also use an excitation matrix array to effectively calculate the
        ! transcorrelation factor everywhere needed!
        temp_ex(1, 1) = src

        call decode_bit_det(nI, ilutI)

        if (is_beta(src)) then
            flip = +1
        else
            flip = -1
        end if

        ! add the spinflipped
        temp_ex(2, 1) = src + flip

        if (present(tgt)) then
            ASSERT(present(cpt))
            cpt = 0.0_dp

            do i = 1, ubound(neighbors, 1)
                elem = 0.0_dp
                if (IsNotOcc(ilutI, neighbors(i))) then
                    temp_ex(1, 2) = neighbors(i) + flip
                    temp_ex(2, 2) = neighbors(i)
                    elem = abs(get_offdiag_helement_heisenberg(nI, temp_ex, .false.))
!                     elem = abs(get_heisenberg_exchange(src, neighbors(i)+flip))
                end if
                if (neighbors(i) == tgt) then
                    cpt = elem
                end if
                cum_sum = cum_sum + elem
            end do
            if (cum_sum < EPS) then
                cpt = 0.0_dp
            else
                cpt = cpt / cum_sum
            end if

        else
            do i = 1, ubound(neighbors, 1)
                elem = 0.0_dp
                if (IsNotOcc(ilutI, neighbors(i))) then
                    ! this is a valid orbital to choose from
                    ! but for the matrix element calculation, we need to
                    ! have the opposite spin of the neighboring orbital!
                    temp_ex(1, 2) = neighbors(i) + flip
                    temp_ex(2, 2) = neighbors(i)
                    elem = abs(get_offdiag_helement_heisenberg(nI, temp_ex, .false.))
!                     elem = abs(get_heisenberg_exchange(src, neighbors(i)+flip))
                end if
                cum_sum = cum_sum + elem
                cum_arr(i) = cum_sum
            end do
        end if

    end subroutine create_cum_list_heisenberg