Since exponentiating a matrix is often needed for tensor networks, the scipy.linalg.expm function is also
implemented for np_conserved Arrays.
When you talk about string order parameters, do you really want to calculate that?
What you wrote suggests that you try to calculate something like
\( \langle A_i \rangle + \sum_{i < k < j} \langle B_k \rangle + \langle C_j \rangle \)
This is different from the usual definition of String order parameters, which would look something like
\( \langle A_i \prod_{i < k < j} B_k C_j \rangle \)
Note that the latter is just a
single expectation value!
To calculate the latter, use the
MPS.correlation_function and provide \(B_k\) as argument "opstr".
Moreover, let me note that if you look at Spin-1/2, choosing \(B_k = \exp(i \pi S_z) \) is probably not what you want, because then B_k would just be minus the identity matrix.