Hi Miles,
I thought the singular values are sorted from the largest to the smallest after "svd". But it seems they are not if we conserve the parity. It is easy to sort by ourselves, but I'm not sure if you missed it for ConserveParity=true. Thanks.
#include "itensor/all.h"
using namespace itensor;
int
main()
{
int N = 8;
auto sites = SpinHalf(N,{"ConserveSz=",false, "ConserveParity",true});
auto ampo = AutoMPO(sites);
for(int j = 1; j < N; ++j)
{
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 sweeps = Sweeps(20); //number of sweeps is 5
sweeps.maxdim() = 10,20,100,100,200;
sweeps.cutoff() = 1E-10;
auto state0 = InitState(sites);
for(int is = 1; is <= N; ++is)
{
state0.set(is,is%2==1 ? "Up" : "Dn");
//state0.set(is,"Dn");
}
auto psi0 = randomMPS(state0);
auto [energy,psi] = dmrg(H,psi0,sweeps,{"Quiet",true,"EnergyErrgoal=",1e-6,"EntropyErrgoal=",1e-5});
println("Ground State Energy = ",energy);
auto b = 4;
//"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(n<=4) println(n, " ", p);
if(p > 1E-14) SvN += -p*log(p);
}
printfln("Across bond b=%d, SvN = %.10f",b,SvN);
return 0;
}
Jin