Two-site reduced density matrix

I think that there might be two problems in your function rdm:

  1. for i=a+1:N should be for i=a+1:b-1
  2. prime(psi1[b],linkinds(psi1,b+1)) should be prime(psi1[b],linkinds(psi1,b))

And before obtain psi1dag, I feel like you should do orthogonalize!(psi1,a) according to the original post Is this a correct way to calculate one-site and two-site entanglement spectrum and entropy? .