Hi Chengshu, ok now I see what you are trying to measure, thanks. The reason measuring the density did not work is that the operator expinentiaion does not commute with taking an expectation value. I.e. the parity operator is sensitive to higher order correlations in the wavefunction not captured by the sum of the local densities.
But to measure an operator like this is relatively straightforward with an MPS. First write the operator as a product of local site operators, using the fact that the n_i operators commute. (I see you've already done this above.) Then sketch out the tensor diagram for the expectation value of the operator. It will look like the diagram for the overlap of the wavefunction with itself, except at every site one of the (1-2n_i) operators is inserted. Finally, compute this diagram by contracting tensor together from left to right, say.
auto P = psi.A(1) * op1 * dag(prime(psi.A(1)));
P *= psi.A(2);
P *= op2;
P *= dag(prime(psi.A(2)));
Etc.
Here by op1, op2, etc I mean the (1-2*n_i) operators which you can make either on the fly or ahead of time. You could make them from operators provided by the site set or by explicitly setting tensor elements.
At the end, P will be a scalar tensor if all goes well and you can call P.real() to extract its value.
Let me know if you have some more questions about it.