Hi Zach,
Thanks for the follow-up question here. So for the case of a pure state which is an MPS, there is a shortcut to getting the kind of rdm you want (where sites 1:k are traced over). The shortcut is to “gauge” the MPS so that the tensors on sites 1,2,…,k are all “left-orthogonal”. This can be done in ITensor by calling orthogonalize!(psi,k+1) on an MPS psi. Afterward, the “orthogonality center” will be at site k+1, meaning that the tensors 1,2,…,k will all be left-orthogonal (and tensors k+2,k+3,...,N will also be right-orthogonal).
After that step, you can show that tracing sites 1:k just amounts to leaving out those tensors (they cancel) and just contracting MPS tensor number k+1 with its primed conjugate over the virtual or link index on bond k. If we call the virtual index for bond k+1 b_{k+1} and the site index for site k+1 s_{k+1} then the resulting rdm will have indices s_{k+1},b_{k+1}, s'_{k+1},b'_{k+1}. Then you can input this rdm to our eigen function in ITensor and factorize it to get its eigenvalues.
I’d have to play around with the code to remember how to extract the QNs. Once you have the QN objects themselves, you can call the val function on them to get their numerical values. Say the name of the QN you want is "Sz" so like inside of a QN object like q-QN("Sz",-1,"Nf",1). Then you do val(q,"Sz") and it will return the number -1.