Hi Wangwei -
While one approach would be to create a new site set (similar to the Hubbard site set but with extra "labels" on the fermions) a simpler and likely better approach is to view your multi-band model as just an extended Hubbard model with a unit cell of more than one site (physically this is often the microscopic origin of multiple bands anyway!).
For example, you might be able to rewrite your model as a two-leg ladder with different hoppings on the rungs versus along the legs.
An easy way to set this up is to use AutoMPO. You can use the sample/exthubbard.cc code as an example to get started from. It's a useful example because it shows how to make a second-neighbor hopping. In general, you just need to choose a 1d ordering of your sites (even if they form a ladder or a 2d system) and then just add the terms you want to the AutoMPO as appropriate.
Let me know if you have additional questions about the specific model you want to simulate -