Thanks for the question. There isn't exactly a built-in feature for what you are asking, but it's straightforward to do. I made a sample code below that hopefully show you the steps you need.

To transform psi into a "vector" (tensor with N indices in the full Hilbert space), you can contract all of the MPS tensors into a single tensor using either a for loop or the `prod`

function (prod does something similar to the capital-Pi symbol in mathematics).

From what I understand of your question, there should be no need to really form the entire density matrix from psi. This is because it's a pure state so just an outer product of psi with itself and therefore much more efficient to deal with each copy of psi separately.

In the code below, after transforming psi into a single tensor psiT, I make another (arbitrary) tensor AT which has two copies of the "site" indices: one set with prime level 0 and the other with prime level 1. Performing the tensor contraction using the ITensor * operation contracts over the unprimed indices, and then the noprime function removes the primes from the remaining indices.

In the end, the resulting density matrix is formally defined as this new state Apsi outer producted with your original state psi (being sloppy here about conjugation, but that's easy to include if dealing with complex numbers or QN-conserving tensors).

Please reply below if you have any follow-up questions.

```
using ITensors
let
N = 10
chi = 2
s = siteinds("S=1/2",N)
# Make a random MPS with sites s
# and bond dimension chi:
psi = randomMPS(s,chi)
# psi could come from a DMRG
# calculation instead
# Contract MPS tensors
# together to make full
# wavefunction tensor:
psiT = ITensor(1.0)
for j=1:N
psiT = psiT*psi[j]
end
# s1 s2 s3 ...
# | | | |
# psiT = ============== ...
#
# Note:
# A shorter way to do the above
# is just `psiT = prod(psi)`
#
# Make a random tensor with
# s and s' indices:
A = randn(2^N,2^N)
AT = ITensor(A,prime.(s)...,s...)
#
# s1' s2' s3' ...
# | | | |
# AT = ============== ...
# | | | |
# s1 s2 s3 ...
#
ApsiT = noprime(AT*psiT)
# Note: afterward can turn ApsiT
# back into an MPS by calling
# `Apsi = MPS(ApsiT,s...;cutoff=1E-10)`
return
end
```