In our daily live as a quantum physicist, a very routine task is to solve fermionic problems: from computing ground states to time evolution. It happens to me often enough but not enough that I keep the same code around for tests. Moreover, sometimes we need to push the limits of the problems we want to solve and then you do not need any code but the right one.
Until this comes for you, you might consider using the tensor and MPS libraries, combining the exact diagonalization routine eigs() with the object in the MPS library that creates fermionic lattice Hamiltonians with up to about 12 particles in 24 sites (26 if you have some extra memory and processors).
The code is quite straightforward. Let's assume you have constructed two matrices, J and U, describing the hopping and interaction between fermions (or hard-core bosons)
Then you would write something like:
int neig = 1; // Number of eigenstates
Finally, while I am right now testing the code, it is provided as is, with no implicit or explicit warranty that it will suit your needs or fulfill your goals or make you coffee. In particular it is quite usual that for complex problems (i.e. topological order in 'J') the eigs() solver may get trapped or the lattice construction of J may be tricky. Be careful. And please report any issues you find.