I think that there might be two problems in your function rdm
:
for i=a+1:N
should befor i=a+1:b-1
prime(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? .