diff --git a/docs/src/versions.md b/docs/src/versions.md
index e570f6576..6057f12e1 100644
--- a/docs/src/versions.md
+++ b/docs/src/versions.md
@@ -3,6 +3,8 @@
 ## v0.5.9
 (In development)
 
+* New function [`suggest_timestep`](@ref) to assist in performing accurate and
+  efficient simulation of spin dynamics.
 
 ## v0.5.8
 (Jan 4, 2024)
diff --git a/examples/03_LLD_CoRh2O4.jl b/examples/03_LLD_CoRh2O4.jl
index 373f03181..e70d479b3 100644
--- a/examples/03_LLD_CoRh2O4.jl
+++ b/examples/03_LLD_CoRh2O4.jl
@@ -25,30 +25,42 @@ minimize_energy!(sys)
 sys = resize_supercell(sys, (10, 10, 10))
 @assert energy_per_site(sys) ≈ -2J*S^2
 
-# Use the stochastic Landau-Lifshitz dynamics to thermalize system into
-# equilibrium at finite temperature via a [`Langevin`](@ref) integrator. The
-# timestep ``Δt`` controls integration accuracy. In `:dipole` mode, it should be
-# inversely proportional to the largest effective field in the system. For
-# CoRh₂O₄, this is the antiferromagnetic exchange ``J`` times the spin magnitude
-# ``S=3/2``. The dimensionless parameter ``λ`` determines the magnitude of
-# Langevin noise and damping terms. A reasonable choice is `λ = 0.2`. The
-# temperature `kT` is linked to the magnitude of the Langevin noise term via a
-# fluctuation-dissipation theorem.
-
-Δt = 0.05/abs(J*S)   # Integration timestep
-λ  = 0.2             # Dimensionless damping time-scale
+# We will be using a [`Langevin`](@ref) spin dynamics to thermalize the system.
+# This involves a damping term of strength `λ` and a noise term determined by
+# the target temperature `kT`.
+
+λ  = 0.2             # Magnitude of damping (dimensionless)
 kT = 16 * meV_per_K  # 16K, a temperature slightly below ordering
-langevin = Langevin(Δt; λ, kT);
+langevin = Langevin(; λ, kT)
+
+# Use [`suggest_timestep`](@ref) to select an integration timestep for the given
+# error tolerance, e.g. `tol=1e-2`. The spin configuration in `sys` should
+# ideally be relaxed into thermal equilibrium, but the current, energy-minimized
+# configuration will also work reasonably well.
 
-# Because the magnetic order has been initialized correctly, relatively few
-# additional Langevin time-steps are required to reach thermal equilibrium.
+suggest_timestep(sys, langevin; tol=1e-2)
+langevin.Δt = 0.025;
+
+# Now run a dynamical trajectory to sample spin configurations. Keep track of
+# the energy per site at each time step.
 
 energies = [energy_per_site(sys)]
 for _ in 1:1000
     step!(sys, langevin)
     push!(energies, energy_per_site(sys))
 end
-plot(energies, color=:blue, figure=(size=(600,300),), axis=(xlabel="Time steps", ylabel="Energy (meV)"))
+
+# Now that the spin configuration has relaxed, we can learn that `Δt` was a
+# little smaller than necessary; increasing it will make the remaining
+# simulations faster.
+
+suggest_timestep(sys, langevin; tol=1e-2)
+langevin.Δt = 0.042;
+
+# The energy per site has converged, which suggests that the system has reached
+# thermal equilibrium.
+
+plot(energies, color=:blue, figure=(size=(600,300),), axis=(xlabel="Timesteps", ylabel="Energy (meV)"))
 
 # Thermal fluctuations are apparent in the spin configuration.
 
@@ -58,8 +70,8 @@ plot_spins(sys; color=[s'*S_ref for s in sys.dipoles])
 # ### Instantaneous structure factor
 
 # To visualize the instantaneous (equal-time) structure factor, create an object
-# [`instant_correlations`](@ref) and use [`add_sample!`](@ref) to accumulated
-# data for the equilibrated spin configuration.
+# [`instant_correlations`](@ref). Use [`add_sample!`](@ref) to accumulate data
+# for each equilibrated spin configuration.
 
 sc = instant_correlations(sys)
 add_sample!(sc, sys)    # Accumulate the newly sampled structure factor into `sf`
@@ -112,18 +124,21 @@ heatmap(q1s, q2s, iq;
 # ### Dynamical structure factor
 
 # To collect statistics for the dynamical structure factor intensities
-# ``I(𝐪,ω)`` at finite temperature, use [`dynamical_correlations`](@ref). Now,
-# each call to `add_sample!` will run a classical spin dynamics trajectory.
-# Longer-time trajectories will be required to achieve greater energy
-# resolution, as controlled by `nω`. Here, we pick a moderate number of
-# energies, `nω = 50`, which will make the simulation run quickly.
-
-ωmax = 6.0  # Maximum  energy to resolve (meV)
+# ``I(𝐪,ω)`` at finite temperature, use [`dynamical_correlations`](@ref). The
+# integration timestep `Δt` used for measuring dynamical correlations can be
+# somewhat larger than that used by the Langevin dynamics. We must also specify
+# `nω` and `ωmax`, which determine the frequencies over which intensity data
+# will be collected.
+
+Δt = 2*langevin.Δt
+ωmax = 6.0  # Maximum energy to resolve (meV)
 nω = 50     # Number of energies to resolve
 sc = dynamical_correlations(sys; Δt, nω, ωmax, process_trajectory=:symmetrize)
 
-# Each sample requires running a full dynamical trajectory to measure
-# correlations, so we here restrict to 5 samples.
+# Use Langevin dynamics to sample spin configurations from thermal equilibrium.
+# For each sample, use [`add_sample!`](@ref) to run a classical spin dynamics
+# trajectory and measure dynamical correlations. Here we average over just 5
+# samples, but this number could be increased for better statistics.
 
 for _ in 1:5
     for _ in 1:100
@@ -132,7 +147,7 @@ for _ in 1:5
     add_sample!(sc, sys)
 end
 
-# Define points that define a piecewise-linear path through reciprocal space,
+# Select points that define a piecewise-linear path through reciprocal space,
 # and a sampling density.
 
 points = [[3/4, 3/4,   0],
diff --git a/examples/04_GSD_FeI2.jl b/examples/04_GSD_FeI2.jl
index 09c50cfd6..ae377a44b 100644
--- a/examples/04_GSD_FeI2.jl
+++ b/examples/04_GSD_FeI2.jl
@@ -74,41 +74,53 @@ sys
 
 # ## Finding a ground state
 
-# Sunny introduces a [Langevin dynamics of SU(_N_) coherent
-# states](https://arxiv.org/abs/2209.01265), which can be used to sample spin
-# configurations from the thermal equlibrium.
-#
-# The [`Langevin`](@ref) integrator requires several parameters. The timestep
-# ``Δt`` controls integration accuracy. In `:SUN` mode, it should be inversely
-# proportional to the largest energy scale in the system. For FeI₂, this is the
-# easy-axis anisotropy energy scale, ``D S^2``. The dimensionless parameter
-# ``λ`` determines the magnitude of Langevin noise and damping terms. A
-# reasonable choice is `λ = 0.2`. The temperature `kT` is linked to the
-# magnitude of the noise via a fluctuation-dissipation theorem.
-
-S = 1
-Δt = 0.05/abs(D*S^2)  # Integration timestep
-λ  = 0.2              # Dimensionless damping time-scale
-kT = 0.2              # Temperature in meV
-langevin = Langevin(Δt; kT, λ);
-
-# Langevin dynamics can be used to search for a magnetically ordered state. For
-# this, the temperature `kT` must be below the ordering temperature, but large
-# enough that the dynamical sampling procedure can overcome local energy
-# barriers and eliminate defects.
+# As [previously observed](@ref "1. Multi-flavor spin wave simulations of FeI₂
+# (Showcase)"), direct energy minimization is susceptible to trapping in a local
+# energy minimum.
 
 randomize_spins!(sys)
-for _ in 1:20_000
+minimize_energy!(sys)
+plot_spins(sys; color=[s[3] for s in sys.dipoles])
+
+# Alternatively, one can search for the ordered state by sampling spin
+# configurations from thermal equilibrium. Sunny supports this via a
+# [`Langevin`](@ref) dynamics of SU(_N_) coherent states. This dynamics involves
+# a damping term of strength `λ` and a noise term determined by the target
+# temperature `kT`.
+
+λ  = 0.2  # Dimensionless damping time-scale
+kT = 0.2  # Temperature in meV
+langevin = Langevin(; λ, kT)
+
+# Use [`suggest_timestep`](@ref) to select an integration timestep for the given
+# error tolerance, e.g. `tol=1e-2`. The spin configuration in `sys` should
+# ideally be relaxed into thermal equilibrium, but the current, energy-minimized
+# configuration will also work reasonably well.
+
+suggest_timestep(sys, langevin; tol=1e-2)
+langevin.Δt = 0.027;
+
+# Sample spin configurations using Langevin dynamics. We have carefully selected
+# a temperature of 0.2 eV that is below the ordering temperature, but large
+# enough to that the dynamics can overcome local energy barriers and annihilate
+# defects.
+
+for _ in 1:10_000
     step!(sys, langevin)
 end
 
+# Calling [`suggest_timestep`](@ref) shows that thermalization has not
+# substantially altered the suggested `Δt`.
+
+suggest_timestep(sys, langevin; tol=1e-2)
+
 # Although thermal fluctuations are present, the correct antiferromagnetic order
-# (2 up, 2 down) is apparent.
+# (2 up, 2 down) has been found.
 
 plot_spins(sys; color=[s[3] for s in sys.dipoles])
 
-# For other systems, it can be much harder to find the magnetic ordering in an
-# unbiased way, and more complicated sampling procedures may be necessary.
+# For other phases, it can be much harder to find thermal equilibrium, and more
+# complicated sampling procedures may be necessary.
 
 # ## Calculating Thermal-Averaged Correlations $\langle S^{\alpha\beta}(𝐪,ω)\rangle$
 #
@@ -131,6 +143,11 @@ for _ in 1:10_000
     step!(sys_large, langevin)
 end
 
+# With this increase in temperature, the suggested timestep has increased slightly.
+
+suggest_timestep(sys_large, langevin; tol=1e-2)
+langevin.Δt = 0.040;
+
 # The next step is to collect correlation data ``S^{\alpha\beta}``. This will
 # involve sampling spin configurations from thermal equilibrium, and then
 # integrating [an energy-conserving generalized classical spin
@@ -141,18 +158,21 @@ end
 # this a real-space calculation, data is only available for discrete ``q`` modes
 # (the resolution scales like inverse system size).
 #
-# To store the correlation data, we initialize a `SampledCorrelations` object by
-# calling [`dynamical_correlations`](@ref). It requires three keyword arguments:
-# an integration step size, a target number of ωs to retain, and a maximum
-# energy ω to resolve. For the time step, twice the value used for the Langevin
-# integrator is usually a good choice.
+# The function [`dynamical_correlations`](@ref) creates an object to store
+# sampled correlations. The integration timestep `Δt` used for measuring
+# dynamical correlations can be somewhat larger than that used by the Langevin
+# dynamics. We must also specify `nω` and `ωmax`, which determine the
+# frequencies over which intensity data will be collected.
 
-sc = dynamical_correlations(sys_large; Δt=2Δt, nω=120, ωmax=7.5)
+Δt = 2*langevin.Δt
+ωmax = 7.5  # Maximum energy to resolve (meV)
+nω = 120    # Number of energies to resolve
+sc = dynamical_correlations(sys_large; Δt, nω, ωmax)
 
 # The function [`add_sample!`](@ref) will collect data by running a dynamical
 # trajectory starting from the current system configuration. 
 
-add_sample!(sc, sys_large)        # Accumulate the sample into `sc`
+add_sample!(sc, sys_large)
 
 # To collect additional data, it is required to re-sample the spin configuration
 # from the thermal distribution. For efficiency, the dynamics should be run long
@@ -162,7 +182,7 @@ for _ in 1:2
     for _ in 1:1000               # Enough steps to decorrelate spins
         step!(sys_large, langevin)
     end
-    add_sample!(sc, sys_large)    # Accumulate the sample into `sc`
+    add_sample!(sc, sys_large)
 end
 
 # Now, `sc` has more samples included:
@@ -195,9 +215,9 @@ lines!(ωs, is[2,:]; label="(π,π,π)")
 axislegend()
 fig
 
-# The resolution in energy can be improved by increasing `nω` (and decreasing `Δt`),
-# and the general accuracy can be improved by collecting additional samples from the thermal
-# equilibrium.
+# The resolution in energy can be improved by increasing `nω`, and the
+# statistical accuracy can be improved by collecting additional samples from the
+# thermal equilibrium.
 #
 # For real calculations, one often wants to apply further corrections and more
 # accurate formulas. Here, we apply [`FormFactor`](@ref) corrections appropriate
diff --git a/src/Integrators.jl b/src/Integrators.jl
index 88283bc5f..636ac5e33 100644
--- a/src/Integrators.jl
+++ b/src/Integrators.jl
@@ -6,11 +6,12 @@ thermal bath, of strength `λ`. One call to the [`step!`](@ref) function will
 advance a [`System`](@ref) by `Δt` units of time. Can be used to sample from the
 Boltzmann distribution at temperature `kT`. An alternative approach to sampling
 states from thermal equilibrium is [`LocalSampler`](@ref), which proposes local
-Monte Carlo moves. For example, use `LocalSampler` to sample Ising-like spins.
+Monte Carlo moves. For example, use `LocalSampler` instead of `Langevin` to
+sample Ising-like spins.
 
 Setting `λ = 0` disables coupling to the thermal bath, yielding an
 energy-conserving spin dynamics. The `Langevin` integrator uses an explicit
-numerical integrator that allows energy drift. Alternatively, the
+numerical integrator which cannot prevent energy drift. Alternatively, the
 [`ImplicitMidpoint`](@ref) method can be used, which is more expensive but
 prevents energy drift through exact conservation of the symplectic 2-form.
 
@@ -46,7 +47,8 @@ constructor interprets its `λ` argument as either ``λ`` or ``λ̃``, for modes
 
 References:
 
-1. [D. Dahlbom et al., Phys. Rev. B 106, 235154 (2022)](https://arxiv.org/abs/2209.01265).
+1. [D. Dahlbom et al., Phys. Rev. B 106, 235154
+   (2022)](https://arxiv.org/abs/2209.01265).
 """
 mutable struct Langevin
     Δt  :: Float64
@@ -56,8 +58,10 @@ mutable struct Langevin
     function Langevin(Δt; λ, kT)
         Δt <= 0 && error("Select positive Δt")
         return new(Δt, λ, kT)
-    end    
+    end
 end
+Langevin(; λ, kT) = Langevin(NaN; λ, kT)
+
 
 """
     ImplicitMidpoint(Δt::Float64; atol=1e-12) where N
@@ -73,8 +77,10 @@ approach eliminates energy drift over long simulation trajectories.
 
 References:
 
-1. [H. Zhang and C. D. Batista, Phys. Rev. B 104, 104409 (2021)](https://arxiv.org/abs/2106.14125).
-2. [D. Dahlbom et al, Phys. Rev. B 106, 054423 (2022)](https://arxiv.org/abs/2204.07563).
+1. [H. Zhang and C. D. Batista, Phys. Rev. B 104, 104409
+   (2021)](https://arxiv.org/abs/2106.14125).
+2. [D. Dahlbom et al, Phys. Rev. B 106, 054423
+   (2022)](https://arxiv.org/abs/2204.07563).
 """
 mutable struct ImplicitMidpoint
     Δt   :: Float64
@@ -85,6 +91,127 @@ mutable struct ImplicitMidpoint
         return new(Δt, atol)
     end    
 end
+ImplicitMidpoint(; atol) = ImplicitMidpoint(NaN; atol)
+
+function check_timestep_available(integrator)
+    isnan(integrator.Δt) && error("Set integration timestep `Δt`.")
+end
+
+"""
+    suggest_timestep(sys, integrator; tol)
+
+Suggests a timestep for the numerical integration of spin dynamics according to
+a given error tolerance `tol`. The `integrator` should be [`Langevin`](@ref) or
+[`ImplicitMidpoint`](@ref). The suggested ``Δt`` will be inversely proportional
+to the magnitude of the effective field ``|dE/d𝐬|`` arising from the current
+spin configuration in `sys`. The recommended timestep ``Δt`` scales like `√tol`,
+which assumes second-order accuracy of the integrator.
+
+The system `sys` should be initialized to an equilibrium spin configuration for
+the target temperature. Alternatively, a reasonably timestep estimate can be
+obtained from any low-energy spin configuration. For this, one can use
+[`randomize_spins!`](@ref) and then [`minimize_energy!`](@ref).
+
+If `integrator` is of type [`Langevin`](@ref), then the damping magnitude ``λ``
+and target temperature ``kT`` will tighten the timestep bound. If ``λ`` exceeds
+1, it will rescale the suggested timestep by an approximate the factor ``1/λ``.
+If ``kT`` is the largest energy scale, then the suggested timestep will scale
+like ``1/λkT``. Quantification of numerical error for stochastic dynamics is
+subtle. The stochastic Heun integration scheme is weakly convergent of order-1,
+such that errors in the estimates of averaged observables may scale like ``Δt``.
+This implies that the `tol` argument may actually scale like the _square_ of the
+true numerical error, and should be selected with this in mind.
+"""
+function suggest_timestep(sys::System{N}, integrator::Langevin; tol) where N
+    (; Δt, λ, kT) = integrator
+    suggest_timestep_aux(sys; tol, Δt, λ, kT)
+end
+function suggest_timestep(sys::System{N}, integrator::ImplicitMidpoint; tol) where N
+    (; Δt) = integrator
+    suggest_timestep_aux(sys; tol, Δt, λ=0, kT=0)
+end
+
+function suggest_timestep_aux(sys::System{N}; tol, Δt, λ, kT) where N
+    acc = 0.0
+    if N == 0
+        ∇Es, = get_dipole_buffers(sys, 1)
+        set_energy_grad_dipoles!(∇Es, sys.dipoles, sys)
+        for (κ, ∇E) in zip(sys.κs, ∇Es)
+            # In dipole mode, the spin magnitude `κ = |s|` scales the effective
+            # damping rate.
+            acc += (1 + (κ*λ)^2) * norm(∇E)^2
+        end
+    else
+        ∇Es, = get_coherent_buffers(sys, 1)
+        set_energy_grad_coherents!(∇Es, sys.coherents, sys)
+        for ∇E in ∇Es
+            acc += (1 + λ^2) * norm(∇E)^2
+        end
+    end
+
+    # `drift_rms` gives the root-mean-squared of the drift term for one
+    # integration timestep of the Langevin dynamics. It is associated with the
+    # angular velocity dθ/dt where dθ ~ dS/|S| or dZ/|Z| for :dipole or :SUN
+    # mode, respectively. In calculating `drift_rms`, it is important to use the
+    # energy gradient |∇E| directly, rather than projecting out the component of
+    # ∇E aligned with the spin. Without projection, one obtains direct
+    # information about the frequency of oscillation. Consider, e.g., a spin
+    # approximately aligned with an external field: the precession frequency is
+    # given by |∇E| = |B|.
+    drift_rms = sqrt(acc/length(eachsite(sys)))
+
+    # In a second-order integrator, the local error from each deterministic
+    # timestep scales as dθ². Angular displacement per timestep dθ scales like
+    # dt drift_rms, yielding err1 ~ (dt drift_rms)^2
+    #
+    # Quantifying the "error" introduced by thermal noise is subtle. E.g., for
+    # weak convergence, we should consider the effect on statistical
+    # observables. We avoid all subtleties by naïvely assuming this error
+    # continues to be second order in `dt`. To determine the proportionality
+    # constant, consider the high-T limit, where each spin undergoes Brownian
+    # motion. Here, the diffusion constant D ~ λ kT sets an inverse time-scale.
+    # This implies err2 ~ (dt λ kT)².
+    #
+    # The total error (err1 + err2) should be less than the target tolerance.
+    # After some algebra, this implies,
+    #
+    # dt ≲ sqrt(tol / (c₁² drift_rms² + c₂² λ² kT²))
+    #
+    # for some empirical constants c₁ and c₂.
+
+    c1 = 1.0
+    c2 = 1.0
+    Δt_bound = sqrt(tol / ((c1*drift_rms)^2 + (c2*λ*kT)^2))
+
+    # Print suggestion
+    bound_str, tol_str = number_to_simple_string.((Δt_bound, tol); digits=4)
+    print("Consider Δt ≈ $bound_str for this spin configuration at tol = $tol_str.")
+
+    # Compare with existing Δt if present
+    if !isnan(Δt)
+        Δt_str = number_to_simple_string(Δt; digits=4)
+        if Δt <= Δt_bound/2
+            println("\nCurrent value Δt = $Δt_str seems small! Increasing it will make the simulation faster.")
+        elseif Δt >= 2Δt_bound
+            println("\nCurrent value Δt = $Δt_str seems LARGE! Decreasing it will improve accuracy.")
+        else
+            println(" Current value is Δt = $Δt_str.")
+        end
+    end
+end
+
+
+function Base.show(io::IO, integrator::Langevin)
+    (; Δt, λ, kT) = integrator
+    Δt = isnan(integrator.Δt) ? "<missing>" : repr(Δt)
+    println(io, "Langevin($Δt; λ=$λ, kT=$kT)")
+end
+
+function Base.show(io::IO, integrator::ImplicitMidpoint)
+    (; Δt, atol) = integrator
+    Δt = isnan(integrator.Δt) ? "<missing>" : repr(Δt)
+    println(io, "ImplicitMidpoint($Δt; atol=$atol)")
+end
 
 
 ################################################################################
@@ -105,6 +232,8 @@ such as [`LocalSampler`](@ref).
 function step! end
 
 function step!(sys::System{0}, integrator::Langevin)
+    check_timestep_available(integrator)
+
     (∇E, s₁, f₁, r₁, ξ) = get_dipole_buffers(sys, 5)
     (; kT, λ, Δt) = integrator
     s = sys.dipoles
@@ -122,7 +251,8 @@ function step!(sys::System{0}, integrator::Langevin)
     set_energy_grad_dipoles!(∇E, s₁, sys)
     @. s = s + 0.5 * Δt * (f₁ + rhs_dipole(s₁, -∇E, λ)) + 0.5 * √Δt * (r₁ + rhs_dipole(s₁, ξ))
     @. s = normalize_dipole(s, sys.κs)
-    nothing
+
+    return
 end
 
 # The spherical midpoint method, Phys. Rev. E 89, 061301(R) (2014)
@@ -132,6 +262,8 @@ end
 #   (s′ - s)/Δt = 2(s̄ - s)/Δt = - ŝ × B,
 # where B = -∂E/∂ŝ.
 function step!(sys::System{0}, integrator::ImplicitMidpoint)
+    check_timestep_available(integrator)
+
     s = sys.dipoles
     (; Δt, atol) = integrator
 
@@ -178,11 +310,14 @@ end
 ################################################################################
 # SU(N) integration
 ################################################################################
-@inline function proj(a::T, Z::T)  where T <: CVec
-    (a - ((Z' * a) * Z))  
+
+@inline function proj(a::T, Z::T) where T <: CVec
+    a - Z * ((Z' * a) / (Z' * Z))
 end
 
 function step!(sys::System{N}, integrator::Langevin) where N
+    check_timestep_available(integrator)
+
     (Z′, ΔZ₁, ΔZ₂, ξ, HZ) = get_coherent_buffers(sys, 5)
     Z = sys.coherents
 
@@ -200,6 +335,8 @@ function step!(sys::System{N}, integrator::Langevin) where N
 
     # Coordinate dipole data
     @. sys.dipoles = expected_spin(Z)
+
+    return
 end
 
 function rhs_langevin!(ΔZ::Array{CVec{N}, 4}, Z::Array{CVec{N}, 4}, ξ::Array{CVec{N}, 4},
@@ -220,6 +357,8 @@ end
 #   (Z′-Z)/Δt = - i H(Z̄) Z, where Z̄ = (Z+Z′)/2
 #
 function step!(sys::System{N}, integrator::ImplicitMidpoint; max_iters=100) where N
+    check_timestep_available(integrator)
+
     (; atol) = integrator
     (ΔZ, Z̄, Z′, Z″, HZ) = get_coherent_buffers(sys, 5)
     Z = sys.coherents
diff --git a/src/OnlineStatistics.jl b/src/OnlineStatistics.jl
new file mode 100644
index 000000000..b56cfae65
--- /dev/null
+++ b/src/OnlineStatistics.jl
@@ -0,0 +1,27 @@
+
+mutable struct OnlineStatistics{T}
+    n::Int
+    μ::T
+    M::Float64
+ 
+    function OnlineStatistics{T}() where T
+        new(0, T(NaN), T(NaN))
+    end
+end
+
+function accum!(o::OnlineStatistics{T}, x) where T
+    o.n += 1
+    if o.n == 1
+        o.μ = x
+        o.M = 0.0
+    else
+        μ_prev = o.μ
+        o.μ += (x - μ_prev) / o.n
+        o.M += real((x - o.μ)' * (x - μ_prev))
+    end
+    o
+end
+
+# If Statistics becomes in scope, define methods Statistics.{mean, var} instead
+mean(o::OnlineStatistics) = o.μ
+var(o::OnlineStatistics; corrected=true) = corrected ? o.M/(o.n-1) : o.M/o.n
diff --git a/src/Sunny.jl b/src/Sunny.jl
index a56e21ef1..1da32746a 100644
--- a/src/Sunny.jl
+++ b/src/Sunny.jl
@@ -82,7 +82,7 @@ include("Reshaping.jl")
 export reshape_supercell, resize_supercell, repeat_periodically
 
 include("Integrators.jl")
-export Langevin, ImplicitMidpoint, step!
+export Langevin, ImplicitMidpoint, step!, suggest_timestep
 
 include("Optimization.jl")
 export minimize_energy!