Hi Javad,
Thanks for the discussion above. I thought of a solution that could work well for what you want to do. I'll just give an outline, though happy to discuss some technical aspects of implementing it in ITensor (C++ v3).
First of all, prepare an MPS which is a product of singlets on sites (1,2), (3,4), (5,6), etc. To do so, you can use code similar to the following:
https://github.com/ITensor/ITensor/blob/a1254d56a180d5e0b992c9b7d9a463bb06a5333d/tutorial/finiteT/ancilla.cc#L58
which is from the "tutorial/finiteT/ancilla.cc" code included with ITensor.
Then use "swap moves" to exchange MPS sites to obtain the pattern of singlets that you want. For example, if you merge the MPS tensors on sites 2 and 3:
auto wf = psi(2)*psi(3);
then SVD these tensors back apart, but now grouping the 3rd physical index with the MPS bond index connecting MPS tensors 1 and 2:
auto b1 = commonIndex(psi(1),psi(2));
auto [U,S,V] = svd(wf,{b1,sites(2)});
Finally, reassign U and (S*V) into the MPS to update it:
psi.ref(2) = U;
psi.ref(3) = S*V;
Now the pattern of singlets will be (1,3),(2,4),(5,6),... and the bond dimensions of the MPS will be increased accordingly, and should now be: 2,4,2,1,2,1,... However note that the site indices will be out of order from the original ones you may want: they will be 1,3,2,4,5,...
Keep repeating the above steps to swap pairs of neighboring physical sites. It's ok if you ignore the MPS gauge (orthogonality) conditions as long as all of the SVD's you do are exact and don't truncate.
Finally, the site indices of the MPS will be very scrambled up compared to the original order. You can either fix this by doing code lik psi.ref(j).replaceInds({sites(n)},{sites(j)}); on each MPS tensor to put the physical index sites(j) back on MPS tensor j. Or you can just work with the new order of site indices if that's convenient to do.
Best regards,
Miles