Let’s say you want to measure the reduced density matrix of the sites c and c+1 of an MPS. Then you can do that with the following code:
orthogonalize!(psi,c)
ket = psi[c]*psi[c+1]
bra = dag(ket)
rho = prime(ket,"Site")*bra
I would strongly suggest you write a tensor diagram and think about the “orthogonality conditions” of an MPS to understand what the code above is doing. Also I give a longer answer in this post about how to get a reduced density matrix for the entire left half of a finite system that could be helpful background for you.
To do more sites, you can generalize the above code to make ket include more MPS tensors, such as ket = psi[c]*psi[c+1]*psi[c+2] for the three-site case. Note that this calculation scales exponentially in the number of sites, but should work well for a small number of sites.