subroutine init_get_helement_heisenberg_guga
if (associated(tmat2d)) deallocate(tmat2d)
allocate(tmat2d(nbasis, nbasis), source=h_cast(0.0_dp))
call setup_exchange_matrix(lat)
call setup_spin_free_exchange(lat)
get_umat_el => get_umat_heisenberg_spin_free
end subroutine init_get_helement_heisenberg_guga