Hi all,

I don't know if it is allowed to answer in the same post but I reconsidered the problem after a while an came with the following idea.

Considering a chain of L sites, with two species of bosons, I wanted to construct the reduced density matrix of the first 4 sites with a given flux QN(("Na",2),("Nb",2)), Na and Nb are the number of bosons of type a and type b. The values of Na and Nb are taken arbitrarly in this example.

I tried this:

ψ=orthogonalize(starting_MPS,1)

ψdag = dag(ψ)

r1=linkind(ψ,1)

site1=1 #first site

site2=4 #last site

ρA = prime(ψ[site1],r1)*prime(ψdag[site1],"Site")

for j in site1+1:site2-1

lj = linkind(ψ,j-1)

rj = linkind(ψ,j)

ρA *= prime(prime(ψ[j],lj),rj)

ρA *= prime(ψdag[j],"Site")

end

l2 = linkind(ψ,site2-1)

l2p = prime(linkind(ψdag,site2-1))

r2 = prime(linkind(ψ,site2))

r2p = prime(linkind(ψdag,site2))

ρA *= prime(ψ[site2],l2)

ρA *= delta(QN(("Na",2),("Nb",2)),l2,l2p)

ρA *= delta(QN(("Na",2),("Nb",2)),r2p,r2)

ρA *= prime(ψdag[site2])

and realized that delta(QN(("Na",2),("Nb",2)),l2,l2p) and delta(QN(("Na",2),("Nb",2)),r2,r2p) were empty.

At first I thought I was selecting an incorrect values of the flux but it is compatible with l2, which is the following

(dim=64|id=983|"Link,fact")

1: QN(("Na",4),("Nb",4)) => 1

2: QN(("Na",4),("Nb",3)) => 3

3: QN(("Na",4),("Nb",2)) => 3

4: QN(("Na",4),("Nb",1)) => 1

5: QN(("Na",3),("Nb",4)) => 3

6: QN(("Na",3),("Nb",3)) => 9

7: QN(("Na",3),("Nb",2)) => 9

8: QN(("Na",3),("Nb",1)) => 3

9: QN(("Na",2),("Nb",4)) => 3

10: QN(("Na",2),("Nb",3)) => 9

11: QN(("Na",2),("Nb",2)) => 9

12: QN(("Na",2),("Nb",1)) => 3

13: QN(("Na",1),("Nb",4)) => 1

14: QN(("Na",1),("Nb",3)) => 3

15: QN(("Na",1),("Nb",2)) => 3

16: QN(("Na",1),("Nb",1)) => 1

for this reason I assumed the problem was in l2p in the sense that dag(index) inverts the sign of the flux as I understood correctly looking at the documentation.

However delta(QN(("Na",2),("Nb",-2)),l2,l2p) does not work as well and I'm a bit confused about the construction of this delta tensor.

For example delta(QN(("Na",2),("Nb",-2)),l2,l2) works and selects the right sector but the indexes are both outwards and the contraction with the other tensors becomes cumbersome.

Thank you for your time and I apologize if this answer goes against the policy of this forum. It is not written in the guidelines.

Best,

Vittorio

EDIT: The idea of multiplying by the delta tensor is to have a tensor which is different from zero only in the subspace I want to construct the reduced density matrix.

EDIT2: I managed to do what I wanted using setelt in some fashion like this one:

temp=setelt(r2p=>96,r2=>96)

for n in 97:124

temp+=setelt(r2p=>n,r2=>n)

end

It is not pretty nor automatic but it works. I put it here in case anyone is interested.