1-RDM from Julia version

This is about the question of computing a RDM from an MPS here. Given an MPS over n spin sites , I want to compute the 1-RDM of the kth site only (all n-1 sites with be traced except site index k). I want to retrieve the 1-RDM and its spectrum thereafter. I followed the post below

However was not able to understand the last response completely. I understand the orthogonalize part but not specifically the lines beyond as mentioned here

What did I try ?

The following is the code I tried

#Function for density matrix
function density_matrix(sites, psi::MPS)
    prime!(psi_diag, "Link")
    return outer(psi,psi_diag)

#Computing 1st order reduced density matrix at each site
function one_RDM(sites, state, index)
    rho=density_matrix(sites, state)
    #rho_tilde = copy(rho)
    #orthogonalize!(rho, index)
    for (i,s) in zip(sort!(setdiff(collect(1:length(sites)), index)), sites)
        rho[i] *= delta(s, prime(s))

    L = ITensor(1.0)
    L = prod(rho)   # conversion to Itensor object to apply eigen methods
    #for i in 1:sort!(setdiff(collect(1:length(sites)),index))
    #    L *= rho[i]      #same as prod above
    return L

It works , but is remarkably slow for large bond dimensions. Following the comment in the above linked post, I am not sure how to proceed beyond the following :

using ITensors
n = 6
s = siteinds("S=1",n;conserve_qns=true)
state = ["Dn" for i = 1:n]
psi = randomMPS(s,state)
k = 4  # k is the physical index where i want the 1RDM not traced (the ^1\rho^k_k' are the indices)

orthogonalize!(psi, k)
for i in k+1:length(s)
     psi[i] = psi[i]*prime(dag(psi[i]),siteinds(psi))
     # What needs to be done next following the discussion to get correct 1-RDM and its spectrum ?

Any help would be appreciated on how to efficiently compute the quantity above (even for larger bond dimensions) Also what would change if I wanted to retrieve the 2-RDM at say site l and m. I can find a C++ post about this ITensor but not a Julia post. Certain commands like rightLinkIndex or leftLinkIndex has no analogue it seems

Hi Manas, here are some notes about how to compute the 1-RDM and the reasoning for the code. It helps to draw it as diagrams:

The main points are:

  • you only need to prime the site index, as you want to contract over the Link indices
  • you should just use the * operator (ITensor contraction) and not the outer function

Regarding your other question, yes in Julia there’s no analogue of leftLinkIndex and rightLinkIndex currently. For simplicity there is just linkind(psi,j) which gives you the Link index between tensors j and j+1.