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