I think that there might be two problems in your function rdm:
for i=a+1:Nshould befor i=a+1:b-1prime(psi1[b],linkinds(psi1,b+1))should beprime(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? .