get_umat_el_heisenberg Function

    pure function get_umat_el_heisenberg(i, j, k, l) result(hel)
        integer, intent(in) :: i, j, k, l
        HElement_t(dp) :: hel
#ifdef DEBUG_
        character(*), parameter :: this_routine = "get_umat_el_heisenberg"
#endif
        ! how exactly do i do that? i access umat with spatial orbitals
        ! so i lose the info if the spins fit.. al i know is that
        ! <ij|kl> i and j are the electrons which must be neighors
        ! and k and l must be in the same orbital as one of the electrons
        ! and what about <ij|kl> = -<ij|lk> symmetry? hm..
        ! i have to think about that and the influence on the matrix element
        ASSERT(allocated(exchange_matrix))

        if (i == j) then
            hel = h_cast(0.0_dp)
        else
            if (i == k .and. j == l) then
                hel = h_cast(exchange_matrix(2 * i, 2 * j - 1))
            else if (i == l .and. j == k) then
                hel = h_cast(exchange_matrix(2 * i, 2 * j - 1))
            else
                hel = h_cast(0.0_dp)
            end if
        end if

    end function get_umat_el_heisenberg