Constructing 3-site reduced density matrix.

ITensor Julia Questions
1

Can someone point out what mistake am I doing in the below code snippet to construct a 3-site reduced density matrix (RDM) ? This method, taken from another post in this forum long time ago, works correctly when it’s a 2-site RDM that I’ve been using in several problems but not when it’s a 3-site (or 4-site) RDM… I am using this within a TEBD evolution, and noticed that the sum of RDM eigvals is > 1 and number of eigvals is \neq 2^3 (after a few time steps, if not immediately so)… I’ve printed the indices, and played around with the prime function, but haven’t succeeded in finding the issue. Perhaps it’s something silly, but let me know what’s gone wrong here.

function RDM3(ψ,A,B,C)
    
    psi1 = deepcopy(ψ)
    
    orthogonalize!(psi1,A)
    
    # 3-site RDM construction for subsystem = {A} U {B} U {C}.
    # A, B, C denote the sites.

    psi1dag=prime(dag(psi1),linkinds(psi1)) 

    rho_ABC = prime(psi1[A],linkinds(psi1,A-1))*prime(psi1dag[A],s[A])
    
    for i=A+1:B-1                  
     rho_ABC *= psi1[i]             
     rho_ABC *= psi1dag[i]
    end

    rho_ABC *= prime(psi1[B],linkinds(psi1,B))*prime(psi1dag[B],s[B])
    
    for i=B+1:C-1
        rho_ABC *= psi1[i]
        rho_ABC *= psi1dag[i]
    end
    
    rho_ABC *= prime(psi1[C],linkinds(psi1,C))*prime(psi1dag[C],s[C])
    
    DABC,UABC=eigen(rho_ABC)
    
    @show length(diag(DABC))
    cc = 0
    for i=1:length(diag(DABC))
        cc += real(diag(DABC))[i]
    end
    @show cc       
    
    return DABC
end
2

Two comments I can see - you still have a link inds in rho_ABC that shouldnt be there (they should be contracted), and for eigen you can specify rho is hermitian eigen(rho_ABC,ishermitian=true)

1 Like
3

Hi Ryan,

thanks for your reply. I modified my code on your suggestion (shown below), and a part of the issue is solved - the sum of RDM eigvals. is now = 1 and the obtained results for quantities depending on this (such as mutual info involving 3 sites, and so on) agrees with a few benchmarking cases available in the literature.

However, the number of RDM eigvals is still \neq 2^3 in most cases. When printing out diag(DABC) I am seeing that in those cases a lot of zeros (or rather, numbers typically of 1e-10 or smaller) are getting included in the DABC, thus not affecting the sum = 1 and other physical results, but I’d still like to figure out the source of this ambiguity but unable to see why this is happening and finally asking back here. This is especially the case when the three sites are not adjacent but there is some distance between two or all of them… When the three sites are adjacent, then aside from the first few time-steps eventually it shows length(diag(DABC)) = 8 but even this appears to depend on the initial state fed into the TEBD as per the experiments I tested out so far. This is weird.

I must mention that the previously posted manner of constructing the RDMs were perfectly fine in all aspects (no. of eigvals, sum of eigvals, agreement with literature results) when the RDMs were 2-site It was apparently only a problem with 3 (or more)-site RDMs). However the same issues come up with 2-RDMs also with this same modification as shown below.

Also, having ishermitian=true makes no difference to these issues.

Another thing worth mentioning is that the rho_ABC here gets constructed as an ITensor object, not an MPO.

function RDM3(ψ,A,B,C)
    
    psi1 = deepcopy(ψ)
    
    orthogonalize!(psi1,A)
    
    # 3-site RDM construction for subsystem = {A} U {B} U {C}.
    # A, B, C denote the sites.

    psi1dag=prime(dag(psi1))

    rho_ABC = prime(psi1[A],s[A])*prime(psi1dag[A],s[A])
    
    for i=A+1:B-1                  
     rho_ABC *= psi1[i]             
     rho_ABC *= psi1dag[i]
    end

    rho_ABC *= prime(psi1[B],s[B])*prime(psi1dag[B],s[B])
    
    for i=B+1:C-1
        rho_ABC *= psi1[i]
        rho_ABC *= psi1dag[i]
    end
    
    rho_ABC *= prime(psi1[C],s[C])*prime(psi1dag[C],s[C])
    
    DABC,UABC=eigen(rho_ABC)
    
    @show length(diag(DABC))
    cc = 0
    for i=1:length(diag(DABC))
        cc += real(diag(DABC))[i]
    end
    @show cc       
    
    return DABC
end

Please let me know if you can see the source of the issue.

4

It seems the RDM you obtain is still not a correct RDM.

s = siteinds("S=1/2",10)
p = randomMPS(s,linkdims=4);
RDM3(p,2,4,6)

I added a @show inds(rho_ABC) inside your funciton and it gives me

inds(rho_ABC) = ((dim=2|id=322|"Link,l=1"), (dim=2|id=322|"Link,l=1")', (dim=2|id=320|"S=1/2,Site,n=3"), (dim=2|id=320|"S=1/2,Site,n=3")', (dim=2|id=693|"S=1/2,Site,n=5"), (dim=2|id=693|"S=1/2,Site,n=5")', (dim=4|id=950|"Link,l=6"), (dim=4|id=950|"Link,l=6")')

A RDM should have entirely link indices or entirely site indices (here you probably want site indices). That is why you are getting the wrong number of eigvals. To help, you can add @assert hastags(rho_ABC,"Link") == false to your code before the eigen call while debugging. This will also speed things up and make the tensor contractions smaller.

5

Putting that assert gives an assertion error, so yeah that’s probably where it’s going wrong though not sure how to correct it yet (funny how this doesn’t seem to matter to 2-RDMs which are carrying both site and link indices). Thanks for this tip.