implement "models" module, optionally incorporate units into it
See here for a rough implementation.
The overall idea is to use some symbolic definition of tight binding models, e.g. using a string "p**2 / 2 m + V", interpret it and process it to produce efficient value functions that can be used in Kwant systems. Adding unit support would be an extra benefit. Most likely this will rely on sympy
for symbolic manipulation.