does anybody know an _easy_ way to compute the square of a MPO/model? I want to compute something like \((\Delta H)^2 = <H>^2 - <H^2>\)
to check how good an Eigenstate my result is.

Defining a new model that is the square of the Hamiltonian by hand seems rather cumbersome (and doesnt lead to a nice nn-coupling model) and I was hoping there is a quick and dirty way to calculate it directly.

I finally took the time to implement at least the "naive" approach of contracting the network with the squared MPO for <psi| H H |psi> in the new variance.
Right now it works only for finite systems - for infinite MPS, this is a whole new story, since there is no direct notion of overlaps any more.
This is not the most efficient way possible, I think it scales as \(\mathcal{O}(\chi_{MPS}^3 \chi_{MPO}^2 d + \chi_{MPS}^2 \chi_{MPO}^2 d^2)\) and is hence more expensive than DMRG itself in the limit of large MPS and MPO bond dimensions!
The more efficient way would be to try to compress the squared MPO, but that's not implemented and also non-trivial...

I actually need it for infinite systems, do you (or anybody else) know of another way to check "how good an eigenstate" the result is other than the variance?