get_helement_heisenberg_general Function

    function get_helement_heisenberg_general(nI, nJ, ic_ret) result(hel)
        integer, intent(in) :: nI(nel), nJ(nel)
        integer, intent(inout), optional :: ic_ret
        HElement_t(dp) :: hel

        integer :: ic, ex(2, maxExcit)
        logical :: tpar
        integer(n_int) :: ilutI(0:NIfTot), ilutJ(0:NIfTot)

        if (present(ic_ret)) then
            if (ic_ret == 0) then
                hel = get_diag_helement_heisenberg(nI)

            else if (ic_ret == 2) then
                ex(1, 1) = 2
                call GetExcitation(nI, nJ, nel, ex, tpar)
                hel = get_offdiag_helement_heisenberg(nI, ex, tpar)

            else if (ic_ret == -1) then
                call EncodeBitDet(nI, ilutI)
                call EncodeBitDet(nJ, ilutJ)

                ic_ret = FindBitExcitLevel(ilutI, ilutJ)

                if (ic_ret == 0) then
                    hel = get_diag_helement_heisenberg(nI)

                else if (ic_ret == 2) then
                    ex(1, 1) = 2
                    call GetBitExcitation(ilutI, ilutJ, ex, tpar)

                    hel = get_offdiag_helement_heisenberg(nI, ex, tpar)

                else
                    hel = h_cast(0.0_dp)
                end if
            else
                hel = h_cast(0.0_dp)
            end if
        else
            call EncodeBitDet(nI, ilutI)
            call EncodeBitDet(nJ, ilutJ)

            ic = FindBitExcitLevel(ilutI, ilutJ)

            if (ic == 0) then
                hel = get_diag_helement_heisenberg(nI)

            else if (ic == 2) then
                ex(1, 1) = 2
                call GetBitExcitation(ilutI, ilutJ, ex, tpar)
                hel = get_offdiag_helement_heisenberg(nI, ex, tpar)

            else
                hel = h_cast(0.0_dp)
            end if
        end if

    end function get_helement_heisenberg_general