Add classical_spin_temperature function

[Lazersmoke]

Re: https://www.sciencedirect.com/science/article/pii/S0010465516301412

We should steal the equation for "spin temperature" used in that paper (and derived in [26]). Here's a reference implementation:

function classical_spin_temperature(sys)
  dHdS = -Sunny.energy_grad_dipoles(sys)
  numer = sum([norm(cross(sys.dipoles[i],dHdS[i]))^2 for i = eachindex(sys.dipoles)])
  denom = 2 * sum([dot(sys.dipoles[i],dHdS[i]) for i = eachindex(sys.dipoles)])
  numer / denom
end

Using this, I annealed a system (with exponentially decreasing temperature schedule) and got the attached plot of "spin temperature" (y axis) vs Langevin temperature (x axis) [both log scale]. As you can see, the spin temperature formula accurately tracks the annealing temperature
[26] https://journals.aps.org/pre/pdf/10.1103/PhysRevE.82.031111

These two annealing runs were done insufficiently slowly, and you can see that the spin temperature doesn't track the thermostat temperature, reflecting this!

So this is a useful tool for people to check their annealing and also would be good for an example in the docs

[kbarros]

This is a nice idea. What should the behavior be for SU(N) mode?

[ddahlbom]

I think this generalizes using our standard tricks (cross product of S and dH/dT becomes commutator of the density matrix and \frakH, inner product becomes trace).

The following code reproduces Sam's results for an SU(2) system (in :SUN mode) up to numerical precision. It seems to correctly estimate the temperature for SU(N) systems.

function classical_spin_temperature(sys::System{N}) where N
    HZ = Sunny.energy_grad_coherents(sys)
    Z = sys.coherents
    numer = sum([real((tr((Z[i]*(HZ[i]') - HZ[i]*Z[i]')^2))) for i in eachsite(sys)])
    denom = 2N*sum([real(tr(Z[i]*(HZ[i])') - tr(sys.interactions_union[i.I[4]].onsite)/N) for i in eachsite(sys)])
    numer / denom
end

There's sensitivity to the fact that \frakH is not necessarily traceless (when there's an arbitrary onsite anisotropy). This needs to be worked out in detail, but the above works empirically.

[EDIT: Refined the original code snippet so that it actually generalizes.]