Solving fermionic Hamiltonians

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:

#include <mps/lattice.h>

...

int L = J.rows();
Lattice problem(L, N); // [1]
CSparse H = problem.Hamiltonian(J, U, 0.0, Lattice::FERMIONS); // [2]
CTensor psi;
int neig = 1; // Number of eigenstates
CTensor energies = linalg::eigs(H, linalg::SmallestReal, neig, &psi); // [3]

1. Line [1] creates an object that scans the list of possible configurations that N fermions moving in L lattice sites can have. This set is typically much smaller than the full set of configurations, 2^N, and leads to a significant memory saving.
2. Using this idea we can now build a Hamiltonian H that is restricted to the subspace of N particles and which is rather sparse (provided that J and U contain a lot of zeros).
3. Finally, we use the Arnoldi method to compute the lowest energy state (linalg::SmallestReal) of this Hamiltonian. Note that we could ask for two, three or more eigenstates replacing neig with the right value.

Note that this whole code is independent from the MPS algorithms but it is distributed jointly with them. If you are interested in the implementation, you may have a look at the github repository. Othewise just follow the instructions in this blog, install the libraries and get going.

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.