How do we take into account QN conservation in self-writing an MPO?

Question

I want to use the iDMRG code in a 2D system, so I am trying to write a 2D version of the Heisenberg.h file for my 2D spinless fermionic system. I managed to write the MPO with the W-matrices, similarly as in the Heisenberg.h code, and also taken into account the Jordan-Wigner transformation.

Still, I do not know how to handle QN conservation (I want the number of fermions to be conserved), and when I run the code I get an error related to the divergence of the IQTensors.

I understand that I should arrange the tensor in non-zero blocks, but in practice how can I do it when adding terms to the W-matrices? I do not see any special line in the Heisenberg.h file related to this (I see that from lines 84-94 the different terms are added and that's all), but I guess there is something more involved in my system, which is 2D and requires the Jordan-Wigner strings mediating interactions.

Thanks a lot for your time,

Sergi

Answer 1

Hi Sergi,
Thanks for the question. This is definitely a rather more complicated, technical topic that we have documentation for at the moment.

To begin with, please look at the question and answer here, which is similar to yours:
http://itensor.org/support/2083/properly-handle-quantum-number-for-system-spin-half-moments

If you understand everything in the discussion there, then you will be a long way toward understanding all of the aspects of the code in Heisenber.h. If not, then please ask some more questions in the comments below about the parts you do not yet understand.

To answer some of your questions above more specifically:

• the convention for the blocking of the tensors in the Heisenberg.h code (and in QN conserving ITensors more generally) is set by how the Index objects are constructed. The ordering and arrangement of QN blocks (QN subspaces) determines, ahead of time, the blocking pattern of the tensors

• when various terms are added into the W tensors in that code, the addition operator of ITensor automatically figures out where the terms should go into the existing W tensor. So if the term being added in has certain non-zero blocks, then because the Index objects are the same (and carry all of the QN details), the ITensor system can figure out on its own how to correctly place the tensor being added into the block structure of the W tensor. So you don’t really have to manage this yourself; you just have to add ITensors together that carry the same total QN flux and it always works correctly (or you get an error if the flux is different).

• for your Jordan-Wigner strings, it is a complication but you may find it’s not too bad of one. In a normal (non fermionic) MPO, there will be identity operators that go before a set of non-identity operators, or after. These will remain identity operators in your case. But there will be some identity operators which go in between a pair of fermionic operators (operators which change the fermion parity, such as C and Cdag). These should be replaced by F operators.

• the 2D aspect is more of a serious complication, but I’m sure you know how to handle it, just mapping terms in 2D into various 1D distances. E.g. on a square lattice of “width” Ny, horizontal, nearest-neighbor terms (which have a 2D distance of 1) have a 1D distance of Ny (assuming you are using the “zig zag” MPS path which always goes up each column, rather than a snaking path which goes up odd columns and down even columns. I recommend the zig zag path.)

Hope that helps you make some progress. One thing I like to do is to devise a bunch of tests of a new MPO I’m making by using the ITensor function inner to compute matrix elements between various product states. Like you can make a product state |phij> which has a single fermion on some site j and |phik> which has a single fermion on site k, then compute <phik|H|phij> and work out by hand what it should be for your H. Then you can check whether you get the right thing from inner. Often you can do this for every single term in the Hamiltonian and thus be nearly sure that it’s bug free. (I say nearly because checking the JW string requires at least two particles to be present, but you can do checks for that too.) Also of course you can do some DMRG calculations in limits where you know what the answer should be.

Best,
Miles

P.S. as mentioned please ask more question in the comments below if you have them

Comments

Dear Miles,

thank you very much for your fast and extended answer, it really helped me a lot. I think that the key is in your second point. If I understood correctly, regarding QN conservation, the main modification that I have to do to Heisenberg.h is:

links.at(l) = Index(QN({"Nf", 0}),1,
                            QN({"Nf",-1}),1,
                            QN({"Nf",+1}),1,
                            Out,
                            ts);

and the rest will be carried out automatically by the ITensor system, am I right? I wasn't understanding properly how the link Index works until I read the related post, and I had this part of the code wrong. Then, I think that I have the other points (JW strings, 2D system) under control, but for sure I will need to check things with this trick that you said in order to debug the Hamiltonian, so thanks again.

---

Glad the answer was helpful! Yes what you said sounds right, depending on what Hamiltonian you’re making of course, in terms of the QN sub spaces needed and their sizes. You can work that out by determining the fluxes of each operator then enforcing the convention of making each MPO a zero-flux ITensor. Then the role of the virtual/link Index QNs are to cancel off those of the operators.

If you make each term that you are adding in out of the same link indices, then yes the ITensor system (the ITensor “+” operator specifically) will handle all the details of summing terms into your MPO tensors. You can print out each step using PrintData to see the contents of the tensors changing as you go if you’d like.

---

Actually I think I would still need to figure out some details for the part:

links.at(l) = Index(QN({"Nf", 0}),?,
                            QN({"Nf",-1}),?,
                            QN({"Nf",+1}),?,
                            Out,
                            ts);

I guess that the 0, +1, -1 are okay, as my Hamiltonian only contains density-density interactions ,which do not change the QN, and hoppings which carry +-1. Still, even after reading carefully your explanation in the related post, I do not how to figure out the value of the dimensions ?'s above. Naively, I would say that for the QNs +-1 I would have dimension 1, as I only have one
type of term for each in the Hamiltonian (c^tc +h.c.).

For the QN(0) I would say that I have dimension 2 for identities, plus 1 for density-density interactions (analogous to S_zS_z), plus 1 for F operators appearing due to the Jordan-Wigner transformation, which would make it 4 overall. But this is not working. I guess my reasoning is wrong and the long-range nature of interactions, with F operators appearing in between make these dimensions more difficult to count, right?

---

Further comment: I now understand that the link indexes are very much related to the matrix representation of my 2D MPO, which is much more involved than the Heisenberg.h one. Therefore, even though it is maybe too technical to ask/answer, should I write the link indexes according to the operators on the first row of my matrix representation? If with Heisenberg.h the first operators on the left are I, I, S_z, Sm, Sp and this gives {0,0,0,-2, 2}, then if the left part of my matrix is something like I, I, N, C, Cdag, F, should it be {0, 0, 0, -1, 1, 0} ?

---

Hi, so to hopefully answer your questions, there is a great deal of freedom of how you arrange the operators inside your MPO construction, and thus the pattern of the link indices can vary quite a lot depending on these choices too. So there's no one standard choice.

For example, in a lot MPO's that you see written in papers, by authors such as McCulloch, Zaletel, etc. you will see an identity in the upper left and lower right corners. This is pretty common. But for QN conserving MPOs I often group those two identities together into the same block. Which leads me to my main point:

A key thing when constructing QN conserving MPOs is to arrange the operators so that when you compute the QN fluxes of the link index subspaces, the subspaces (settings, or values of the Index) which have the *same flux get grouped together*. This allows for maximum blocking of non-zero elements together. So that's the pattern I always follow: first the flux-zero sector which includes things like the starting and ending identity operators, as well as on-site operators or density-density operators. Then the rest of the flux sectors after that.

So I would do the following:
1. write out the matrix form of your MPO the best way you know how.
2. Now label the rows and columns by the fluxes, which are uniquely determined by requiring every MPO tensor to have flux zero (so once you start on the left hand side, and use the fluxes of the individual operators inside each tensor, the fluxes of the links will get determined).
3. Then ask: are all of the fluxes of the same value grouped together?
4. If not, swap rows and corresponding columns so that they do get grouped more togther. Such a swap will be a gauge transformation of the MPO and will not change the overall operator it makes.
5. Continue doing these row/column swaps until all of the same-flux Index values are grouped together, then you're done.

If you follow the above steps it should give you the MPO that you want. It's also technically ok to have more than one subspace of a QN Index in ITensor with the same QN value. It's just that it will lead to ITensors with more blocks than necessary.

Miles