Manual construction of 2-sites local MPO

Question

Hi ITensor team,

I am trying to manually construct a simple MPO, whose local MPO tensors acts on 2 spin-1/2 sites. Denoting the basis of 2-sites with sigma = 0,1,2,3, each local MPO tensor W^{sigma,sigma'}_{a,b} is diagonal in sigma basis:

W{ab}^{sigma, sigma'} = \delta{sigma, sigma'} A_{ab}^{sigma}

where A_{ab}^sigma are matrices of dimension D, therefore a,b = 1, ... ,D .

I was wondering what is the most efficient way to write down this object in ITensor (C++, v3) and if there are some related examples in the code.

Thanks a lot for your help!

All the best,
Jacopo.

Answer 1

Hi Jacopo,
Good question: yes this is a perfect thing to use ITensor for. There are a few different ways you can proceed. I'll outline some and we could chat more in comments below.

One way is to use the AutoMPO system if you MPO can be expressed as a sum of local operators. If so, then this would be the easiest way to construct your MPO.

A more flexible and powerful way is to just make your MPO directly, by forming the individual tensors and setting their elements. Essentially all you have to do is pre-define indices, construct ITensors, then set their elements in the usual way. Let us know if you need help with any of these steps. But a different pattern I have used is to "outer product" smaller tensors together and add them up to make MPO tensors. A specific example of that pattern is shown in the file sample/Heisenberg.h which is included with the C++ ITensor source code under the sample folder. It makes an MPO for the Heisenberg spin chain Hamiltonian. If that pattern is close to what you are looking for, please let me know if you have some questions about how that code works.

Best,
Miles

Comments

Hi Miles,
thanks a lot for the reply. Actually the MPO to write is really simple, it can be written locally (on 2 -sites) as
W = P_1 A_{ab} + (1-P_1) B_{ab}
where P_1 is the projector on the identity 1/sqrt[2] ( |++> + | - - > ).
I was wondering if I can construct it by accessing the elements of an MPO H, creating two itensor A and B and then doing SVD to assign H.position(b) and H.position(b+1), as for example you do to construct singlet states in the ancilla code?
Thanks and apologies for the possibly trivial question!
Jacopo

---

Hi Jacopo,
Thanks for explaining your question a bit more. So yes, here are some operations you can do with ITensor, and by combining them you can make an arbitrary MPO. You can:
1. create ITensors and set their elements to any values
2. SVD an ITensor to factorize it
3. assign ITensors to elements of an MPO

To do number 3, you actually do this:
MPO H(N);
H.ref(1) = H1;
H.ref(2) = H2;
...
H.ref(N) = HN;

The .ref(n) function of an MPS or MPO allows you to write to the n'th tensor of that MPS or MPO. (The .position function has a different purpose: it changes the 'gauge' of an MPS or MPO.)

To use an SVD in combination with making an MPO, you would SVD the ITensor you want into 3 pieces: U,S,V, then multiply S into either U or V, then finally assign the two MPO tensors one with (U*S) and one with V, say. So like:

auto [U,S,V] = svd(T,uinds);
H.ref(n) = (U*S);
H.ref(n+1) = V;

Please let me know if you have more questions.

Miles

---

Dear Miles,

thanks a lot for your reply and apologizes for my late response.
Here again just to know if I am doing things correctly:

suppose I want my MPO to be product of 2-sites projectors P_{Id} on |00>  + |11> with a matrix B_{a,b} = \delta_{a,b} auxiliary indices with dim = 2 (contracted with (1,1) on the left and (1,1) on the right), namely
W = W^1 ... W^N
with
W^i_{a,b} = P_{Id} B_{a,b}.

I copy my (draft) code here :



    auto W =  MPO(sites);
    
   
    for(int n = 1; n <= 2*N; n += 2)
    {
        auto s1 = sites(n);
        auto s2 = sites(n+1);
        auto s1p = prime(sites(n));
        auto s2p = prime(sites(n+1));
        
        auto i = Index(2,"Link");
        auto j = Index(2,"Link");
        
        auto AA = ITensor(s1,s2,s1p,s2p);
  // define the projector on | 00 >  + | 11 >
        
       AA.set(s1(1),s2(1),s1p(1),s2p(1), 1./2.);
        AA.set(s1(2),s2(2),s1p(2),s2p(2), 1./2. );
        AA.set(s1(1),s2(1),s1p(2),s2p(2), 1./2. );
        AA.set(s1(2),s2(2),s1p(1),s2p(1), 1./2. );
        
  // define the auxiliary indices for matrix BB, tensor product with AA and SVD to get one site MPOs

        if(n>1 && n< 2*N){
        
        auto BB =ITensor(i,j);
        
        BB.set(i(1),j(1), 1);
        BB.set(i(2),j(2), 1);
        
        AA = AA * BB;
            
        ITensor D;
        W.ref(n) = ITensor(s1,s1p,i);
        W.ref(n+1) = ITensor(s2,s2p,j);
        svd(AA,W.ref(n),D,W.ref(n+1));
        W.ref(n) *= D;
        }
        
        if(n==1){
        
        auto BB =ITensor(i);
            
        BB.set(i(1), 1);
        BB.set(i(2),1);
            
        AA = AA * BB;
            
            
        ITensor D;
        W.ref(n) = ITensor(s1,s1p);
        W.ref(n+1) = ITensor(s2,s2p,i);
        svd(AA,W.ref(n),D,W.ref(n+1));
        W.ref(n) *= D;
        }

        if(n==2*N){
        auto BB =ITensor(j);
            
        BB.set(i(1), 1);
        BB.set(i(2), 1);
            
        AA = AA * BB;
        
        ITensor D;
        W.ref(n) = ITensor(s1,s1p,i);
        W.ref(n+1) = ITensor(s2,s2p);
        svd(AA,W.ref(n),D,W.ref(n+1));
        W.ref(n) *= D;
            
            
        }
        
    }

I would be super grateful if you can let me know if this code does what I think and if now what is wrong :)
Thanks a lot!
Jacopo