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)
psi_diag=dag(prime(psi,sites))
prime!(psi_diag, "Link")
return outer(psi,psi_diag)
end
#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))
end
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
#end
@show(tr(L))
return L
end
```

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 **rightLinkInde**x or **leftLinkIndex** has no analogue it seems