Suggest integration timestep

[kbarros]

Collect statistics from numerical integration to provide suggestions about the choice of dt.

[kbarros]

Intended to address first part of https://github.com/SunnySuite/Sunny.jl/issues/149.

As currently implemented, this PR collects statistics of grad E over the course of a dynamical trajectory and stores them in the Langevin struct. An alternative implementation would have check_timestep take a System, and suggest a timestep based on a single grad E calculation. Switching to that alternative implementation would be simpler, and allows the possibility of getting a dt estimate prior to running dynamics. A disadvantage to switching is that the estimates might become worse for small systems where there are less statistics (?).

Another question is how to define an interpretable error tolerance parameter, tol. Since Heun is weak first-order accurate, the error in statistical observables likely decays as 1 / dt. It might be the case that if we switch to implicit-midpoint as our Langevin integrator, the error decays faster, as 1 / dt^2. We could determine this with some numerical experiments.

Prior to merging, we should also add support for ImplicitMidpoint.

[ddahlbom]

I think that collecting statistics for a chosen time steps is reasonable, even if in an ideal world a user would like to have a function that takes one step and gets a suggestion. I assume the statistics will vary a lot depending on whether the system is equilibrated or not, and the user should probably be encouraged to do a small study for the situation that they are interested in studying. I think the report that this provides is very useful.

[kbarros]

I assume the statistics will vary a lot depending on whether the system is equilibrated or not, and the user should probably be encouraged to do a small study for the situation that they are interested in studying. I think the report that this provides is very useful.

A worst case scenario is maybe a single-site problem, where norm(HZ) could be any eigenvalue of the single-ion Hamiltonian, depending on state Z.

Probably for a system with multiple sites, a random spin configuration will give a good sense of the typical dt required.

Conversely, for a local Monte-Carlo sampler, I don't see anyway around the need to collect statistics over a trajectory. Given that this type of interface must exist, it might make sense to use it consistently for Langevin sampling.

[kbarros]

I tested the alternate implementation, which suggests Δt based on the energy gradient of a single snapshot.

Results for the dynamics in our two main tutorials. The suggested Δt below varies according to the type of spin snapshots used in the estimation:

Tutorial03, CoRh2O4

Δt ~ 0.02528 (Polarized in z) Δt ~ 0.04959 (Randomized spins) Δt ~ 0.02528 (Energy minimum) Δt ~ 0.03933 (T=16K equilibrium snapshot) Δt ~ 0.04201 (T=16K, equilibrium averaged)

Tutorial04, FeI2

Δt ~ 0.1482 (Polarized in z) Δt ~ 0.05107 (Randomized spins) Δt ~ 0.02734 (Energy minimum) Δt ~ 0.03794 (T=3.5meV, equilibrium snapshot) Δt ~ 0.03885 (T=3.5meV, equilibrium averaged)

Key points:

1. Energy gradient is larger for configurations at lower temperatures, leading to smaller suggested dt.
2. At finite temperatures, one should in principle average over many equilibrium snapshots. For these particular systems, with 64 sites, this time averaging isn't overkill, because a single snapshot already includes an average over sites.
3. The suggested dt at T=0 is approximately 30% lower than at temperatures around than transition.

I think a reasonable advice to users is "first run minimize_energy! and then call suggest_timestep. This is a simple interface, and can be done prior to any time-integration.

[kbarros]

@ddahlbom pointed out that there might be a different trend with frustrated systems. This can be verified in the SpinW tutorial:

SW08, Kagome antiferro:

Δt ~ 0.025 (Polarized in z) Δt ~ 0.045 ± 0.01 (Randomized spins) Δt ~ 0.05 (Energy minimum)

In this case, both energy minimizing and random spin configurations lead to about the same suggested dt, but curiously the "spin polarized" configuration requires dt smaller by a factor of 1/2.