Compute entanglement entropy of an MPS

//Run DMRG to get a ground state `psi`
auto N = 20;
auto sites = SpinHalf(N);
auto ampo = AutoMPO(sites);
for( auto j : range1(N-1) )
    {
    ampo += 0.5,"S+",j,"S-",j+1;
    ampo += 0.5,"S-",j,"S+",j+1;
    ampo +=     "Sz",j,"Sz",j+1;
    }
auto H = toMPO(ampo);
auto state = InitState(sites,"Up");
for( auto n : range1(N) ) if( n%2==0 ) state.set(n,"Dn");
auto sweeps = Sweeps(5); //number of sweeps is 5
sweeps.maxdim() = 10,20,100,100,200;
sweeps.cutoff() = 1E-10;
auto [energy,psi] = dmrg(H,randomMPS(state),sweeps,"Silent");

//Given an MPS called "psi",
//and some particular bond "b" (1 <= b < length(psi))
//across which we want to compute the von Neumann entanglement
auto b = 10;

//"Gauge" the MPS to site b
psi.position(b); 

//SVD this wavefunction to get the spectrum
//of density-matrix eigenvalues
auto l = leftLinkIndex(psi,b);
auto s = siteIndex(psi,b);
auto [U,S,V] = svd(psi(b),{l,s});
auto u = commonIndex(U,S);

//Apply von Neumann formula
//to the squares of the singular values
Real SvN = 0.;
for(auto n : range1(dim(u)))
    {
    auto Sn = elt(S,n,n);
    auto p = sqr(Sn);
    if(p > 1E-12) SvN += -p*log(p);
    }
printfln("Across bond b=%d, SvN = %.10f",b,SvN);