Hi Nick,
I think you may want to extract the itensor on the ith site form the MPS, multiply the itensor by the op function, and then plug it back. For example like this,
psi.position(i);
auto newpsi = noprime(psi(i) * op(sites,"Sz", i));
psi.set(i, newpsi);
Best,
Yixuan