# 1-RDM from Julia version

**URL:** <https://itensor.discourse.group/t/1-rdm-from-julia-version/354>\
**Category:** ITensor Julia Questions\
**Created:** [August 18, 2022, 1:22am UTC](https://itensor.discourse.group/t/1-rdm-from-julia-version/354 "2022-08-18T01:22:40Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Manas](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/manas/32/376_2.png) [@Manas](https://itensor.discourse.group/u/Manas)\
**Post date:** [August 18, 2022, 1:22am UTC](https://itensor.discourse.group/t/1-rdm-from-julia-version/354/1 "2022-08-18T01:22:40Z")

</div>

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

> [@How to calculate flux of a reduced density matrix?](https://itensor.discourse.group/t/how-to-calculate-flux-of-a-reduced-density-matrix/243/6):
>
> Hi Miles, Sorry for the late update. It seems i have some troubles for obtaining rdm and its factorization. For obtaining rdm from a MPS |\psi\rangle, a naive way is to call outer(psi',psi) (for example, [How to get density matrix operator from MPS wave function? - ITensor Support Q&A](http://itensor.org/support/3217/how-to-get-density-matrix-operator-from-mps-wave-function?show=3217#q3217)) and partially traced out half of the system. (for example, the discussion here [trace and partial trace of MPO - #3 by jisutich](https://itensor.discourse.group/t/trace-and-partial-trace-of-mpo/79/3)) However, this method seems to be very slow when the bond dimension is large. The…

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

> [@How to calculate flux of a reduced density matrix?](https://itensor.discourse.group/t/how-to-calculate-flux-of-a-reduced-density-matrix/243/7):
>
> Hi Zach, Thanks for the follow-up question here. So for the case of a pure state which is an MPS, there is a shortcut to getting the kind of rdm you want (where sites 1:k are traced over). The shortcut is to “gauge” the MPS so that the tensors on sites 1,2,…,k are all “left-orthogonal”. This can be done in ITensor by calling orthogonalize!(psi,k+1) on an MPS psi. Afterward, the “orthogonality center” will be at site k+1, meaning that the tensors 1,2,…,k will all be left-orthogonal (and tensors k…

**What did I try ?**

The following is the code I tried

```julia
#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 :

```julia
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](http://itensor.org/docs.cgi?vers=cppv3&page=formulas/mps_two_rdm) but not a Julia post. Certain commands like **rightLinkInde** x or **leftLinkIndex** has no analogue it seems

---

<div class="post-metadata">

**Author:** ![miles](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/miles/32/6_2.png) [@miles](https://itensor.discourse.group/u/miles)\
**Post date:** [August 18, 2022, 3:07pm UTC](https://itensor.discourse.group/t/1-rdm-from-julia-version/354/2 "2022-08-18T15:07:01Z")

</div>

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:

 ![IMG_FA4A50D42792-1](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/048df29f1230bf9482806ef8507864713456cfb6.jpeg)

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`.