Not sure if it’s relevant anymore (perhaps you already figured it out), but in case it is (and also for record as this thread might come up in search results for many of us looking for similar questions) –
for nonadjacent sites such as for example a 3-site RDM with config. AB-C (where C is separated from AB by a distance of one site), just do the following :
orthogonalize!(psi1,A)
ket1 = psi1[A]
for i=A+1:C
ket1 *= psi1[i]
end
bra1 = dag(ket1)
rho = prime(ket1,"Site")*bra1
#@show inds(rho) <-- this shows that the indices of rho corresponding to the site between B and C are the 3rd and 7th element of the inds(rho) vector.
rho_ABC = rho*delta(inds(rho)[3],inds(rho)[7]) # delta(inds,inds) traces out the given inds (see the documentation)
This will yield the correct sum of eigenvalues = 1 and the correct number of eigenvalues 2^i, where i=3 in this example.
This is now readily generalizable to setups like AB–CD, AB–C–D and so on.