The relax() stopping criterion seems wrong to me. The current relax() procedure is:

1. Run until the total energy stops decreasing.
2. Run until the average torque stops decreasing.
3. Repeat step 2, decreasing MaxErr until MaxErr<=1e-9.

There are some issues that seem wrong to me:

• The torque criterion should consider the maximum torque and not the average torque. Often, only a very small subset of the spins are moving, and relax() should keep running until those spins relax. The maximum torque will measure those fastest spins, independently of the number of already relaxed spins. Moreover, the average torque is noisier, as it's a sum of lots of near-zero values.
• The torque criterion should not be that the torque stops decreasing, it should be that the torque nearly reaches zero (within the numerical noise). A slowly relaxing system may show a constant (non-zero) torque for a long time, such as a domain wall moving in a small constant field. Also, in many relaxation situations the torque may increase, before the system is relaxed. For example, a depinning domain wall at its critical field will start by deforming, slower and slower, up to a point when it depins and the torques increase -- relax() should not stop in this case.
• if the initial MaxErr is <=1e-9, the step 2 will not run.

I suggest that the maximum torque be taken instead of the average torque, and that relax() stop when maxTorque reaches a user-definable threshold. I'm not very familiar with the code, but something like changing lines 62-76 of relax.go by:

    var thresholdTorque := 1e-3; //user-settable
    maxTorque := func() float32 {
        return cuda.MaxVecNorm(solver.k1)
    }
     for maxTorque() < thresholdTorque && !pause {
            relaxSteps(N)
        }
    for MaxErr > 1e-9 && !pause {
        MaxErr /= math.Sqrt2
        for maxTorque() > thresholdTorque && !pause {
            relaxSteps(N)
        }
    }

ed: I've written the code implementing these changes in pull request #148 .