Improve Hamiltonian evaluation for vectorized systems

So far, the new system.VectorizedSystem offers a hamiltonian_submatrix method that offers a subset of the functionality of the eponymous method of system.System.

This is a good first step, but we are not bound by backwards compatibility here (the new hamiltonian_submatrix already breaks it), since there is no legacy code that enables vectorization.

This gives as a unique opportunity to modernize the API for evaluating Hamiltonians.

Here are some problems with the old API:

• It's incompatible with ND-systems.
• It's not possible to selectively evaluate terms (for example only terms that depend on a particular parameter could be of interest).
• By default, a dense matrix is returned (#244).
• It's impossible to act with the Hamiltonian on a vector without creating a sparse matrix in memory.
• inter_cell_hopping may return non-square matrices.