Hi Victor,
Interesting question. So yes, you can efficiently measure S^2 for a large system by representing the total S^2 operator as an MPO. (You can use the overlap(psi,S2,psi) function to efficiently compute the expectation value of an MPO S2 with respect to an MPS psi.)
Since this is something I've needed for projects in the past, I am pasting below a function which will generate an MPO for the total S^2 operator using a Hubbard site set to define the Hilbert space of the lattice. If you're doing a different kind of system, like a spin system, you should be able to adapt this code pretty straightforwardly for your case, but let me know if you encounter any problems. The main change you'd need to make is to define the QN objects going into the link/virtual indices to not have the "Nf" quantum number anymore.
The code below makes S2 as an IQMPO, but if you are not using quantum numbers you should be able just to convert the return value to an MPO like this:
MPO S2 = makeS2(sites);
Oh, and if this solves your problem, could I recruit you for a task? This would be a nice function to turn into a "code formula" to go on the ITensor website (http://itensor.org/docs.cgi?page=formulas). Could you do a brief writeup of it (see the source for the other code formulas here:
https://github.com/ITensor/website/tree/master/docs/formulas
Of course you don't have to if you're too busy. Just an idea & I'll put your name on the contributors page with a badge.
Ok here is the function:
IQMPO
makeS2(Hubbard const& sites)
{
auto N = sites.N();
auto S2 = IQMPO(sites);
auto links = std::vector<IQIndex>(N+1);
for(auto n : range(N+1))
{
links.at(n) = IQIndex(nameint("L",n),
Index("0",3),QN("Nf=",0,"Sz=",0),
Index("+",1),QN("Nf=",0,"Sz=",-2),
Index("-",1),QN("Nf=",0,"Sz=",+2));
}
for(auto n : range1(N))
{
auto row = dag(links.at(n-1));
auto col = links.at(n);
auto& W = S2.Aref(n);
W = IQTensor(row,col,dag(sites(n)),prime(sites(n)));
W += sites.op("Id",n) * row(1) * col(1);
W += sites.op("Id",n) * row(2) * col(2);
W += sites.op("S2",n) * row(2) * col(1);
W += 2*sites.op("Sz",n) * row(2) * col(3);
W += sites.op("Id",n) * row(3) * col(3);
W += sites.op("Sz",n) * row(3) * col(1);
W += sites.op("S+",n) * row(2) * col(4);
W += sites.op("Id",n) * row(4) * col(4);
W += sites.op("S-",n) * row(4) * col(1);
W += sites.op("S-",n) * row(2) * col(5);
W += sites.op("Id",n) * row(5) * col(5);
W += sites.op("S+",n) * row(5) * col(1);
W.scaleTo(1.);
}
S2.Aref(1) *= setElt(links.at(0)(2));
S2.Aref(N) *= setElt(dag(links.at(N))(1));
return S2;
}