Defining the spin glass susceptibility

[mgallegoq]

I wanted to implement the following operator $\chi_{SG} = \frac{1}{N} \sum_{ij}\langle \sigma_i^z \sigma_j^z \rangle^2$.

I have noticed that the implementation of a similar operator $\sum_{ij}\langle \sigma_i^z \sigma_j^z \rangle$ is much faster when you define the operator as the sum of local operators and then calculate the expectation value, rather than calculating each individual term and then performing the summation.

Of course if you naively try add the square as sum(sum((sigmaz(hi, i) * sigmaz(hi, j))**2 for j in range(i, L)) for i in range(L)) wouldn't work.

Is there a way to implement this operator in a similar, more efficient way?

[PhilipVinc]

What error do you get exactly?
If I run your code I see a

TypeError: unsupported operand type(s) for ** or pow(): 'LocalOperator' and 'int'

which is just because we have not defined th **operator. you can just do op*op and it will work fine

In [10]: sum(sum((sigmaz(hi, i) * sigmaz(hi, j))*(sigmaz(hi, i) * sigmaz(hi, j)) for j in range(i, L)) for i in range(L))
Out[10]: LocalOperator(dim=20, #acting_on=210 locations, constant=0.0, dtype=float64)

however this operator is trivial? it's just the identity..