Hi Miles,
Thanks you so much. It worked. I used it following way:
psi=psi0;
orthogonalize!(psi,j)
psidag_j = dag(prime(psi[j], "Site"))
SPIN+ = scalar(psidag_j * op(sites, "S+", j) * psi[j])
SPIN- = scalar(psidag_j * op(sites, "S-", j) * psi[j])
SPINz= scalar(psidag_j * op(sites, "Sz", j) * psi[j])
Also, is there any command to get complex output rather than scalar ( e.g. to get expectation value of S_y or S_x) like "InnerC or eltC" defined in ITensor C++?