get_umat_heisenberg_spin_free Function

    pure function get_umat_heisenberg_spin_free(i, j, k, l) result(hel)
        ! for the spin-free form, I do not need information about
        ! the spin-orbitals
        integer, intent(in) :: i, j, k, l
        HElement_t(dp) :: hel
#ifdef DEBUG_
        character(*), parameter :: this_routine = "get_umat_heisenberg_spin_free"
#endif

        ASSERT(allocated(spin_free_exchange))

        if (i == j) then
            hel = h_cast(0.0_dp)
        else
            if (i == k .and. j == l) then
                hel = h_cast(spin_free_exchange(i, j))
            else if (i == l .and. j == k) then
                hel = h_cast(spin_free_exchange(i, j))
            else
                hel = h_cast(0.0_dp)
            end if
        end if

    end function get_umat_heisenberg_spin_free