Ok, so in the file Heisenberg.h from line 94, the multiplication of two IQIndices (one In, one Out) with an IQTensor/ITensor yields a 4-tensor with only the block corresponding to that particular operator and all other blocks empty? And these are then added to form the full MPO? I am familiar with this form of the MPO but not familiar with McCullough's articles (yet).

What I cannot get is that is the Heisenberg example the MPO has a dimension of (5,5,(2*1 + 1), (2*1 + 1)) since S=1. But you already know that the MPO will have dimension 5 from the number and size of the sectors. For example in the Ising model (spin 1/2), the MPO should have dimension 3 (I suppose?), but I dont know where to put the elements. I may be missing something crucial theory-wise, since I am very new to MPS!

Also, I'm not sure from what you write how I should find the flux. Each "block" has two indices (I assume that you are talking about the operators in each "cell" of the matrix of operators, when the MPO is viewed as a matrix given by doing a kronecker product over indices 1,3 and 2,4 respectively).