Just to say that the approach posted in Is this a correct way to calculate one-site and two-site entanglement spectrum and entropy? for calculating the n-site reduced density matrix works correctly. I used it to calculate the mutual information for an example study involving the transverse field Ising model and it successfully benchmarked with a paper that I was following for the results.
I’d found this thread from a search related to my query about correctly extracting the reduced density matrices. Hence thought of posting about what worked for me.