From ac70cd3a73fdf179771e680a9a8ebdbdd8f0769e Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Tue, 16 Jan 2024 10:12:57 -0700 Subject: [PATCH 1/6] Add `check_timestep` This should be used after a time-integration to verify expected accuracy. --- docs/src/versions.md | 2 + examples/03_LLD_CoRh2O4.jl | 5 ++ examples/04_GSD_FeI2.jl | 7 +++ src/Integrators.jl | 120 +++++++++++++++++++++++++++++++++++-- src/OnlineStatistics.jl | 27 +++++++++ src/Sunny.jl | 4 +- 6 files changed, 158 insertions(+), 7 deletions(-) create mode 100644 src/OnlineStatistics.jl diff --git a/docs/src/versions.md b/docs/src/versions.md index e570f6576..db62a3061 100644 --- a/docs/src/versions.md +++ b/docs/src/versions.md @@ -3,6 +3,8 @@ ## v0.5.9 (In development) +* New function [`check_timestep`](@ref) to verify accuracy of spin dynamics + integration. ## v0.5.8 (Jan 4, 2024) diff --git a/examples/03_LLD_CoRh2O4.jl b/examples/03_LLD_CoRh2O4.jl index 373f03181..b7e3054eb 100644 --- a/examples/03_LLD_CoRh2O4.jl +++ b/examples/03_LLD_CoRh2O4.jl @@ -50,6 +50,11 @@ for _ in 1:1000 end plot(energies, color=:blue, figure=(size=(600,300),), axis=(xlabel="Time steps", ylabel="Energy (meV)")) +# After running a Langevin trajectory, it is a good practice to call +# [`check_timestep`](@ref). + +check_timestep(langevin; tol=1e-2) + # Thermal fluctuations are apparent in the spin configuration. S_ref = sys.dipoles[1,1,1,1] diff --git a/examples/04_GSD_FeI2.jl b/examples/04_GSD_FeI2.jl index 09c50cfd6..117324520 100644 --- a/examples/04_GSD_FeI2.jl +++ b/examples/04_GSD_FeI2.jl @@ -102,6 +102,11 @@ for _ in 1:20_000 step!(sys, langevin) end +# After running a Langevin trajectory, it is a good practice to call +# [`check_timestep`](@ref). + +check_timestep(langevin; tol=1e-2) + # Although thermal fluctuations are present, the correct antiferromagnetic order # (2 up, 2 down) is apparent. @@ -131,6 +136,8 @@ for _ in 1:10_000 step!(sys_large, langevin) end +check_timestep(langevin; tol=1e-2) + # 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 diff --git a/src/Integrators.jl b/src/Integrators.jl index 88283bc5f..9b0128378 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -53,10 +53,14 @@ mutable struct Langevin λ :: Float64 kT :: Float64 + # Averages the square of the dynamical drift term for use in + # `check_timestep`. + drift_squared :: OnlineStatistics{Float64} + function Langevin(Δt; λ, kT) Δt <= 0 && error("Select positive Δt") - return new(Δt, λ, kT) - end + return new(Δt, λ, kT, OnlineStatistics{Float64}()) + end end """ @@ -80,13 +84,95 @@ mutable struct ImplicitMidpoint Δt :: Float64 atol :: Float64 + # Averages the square of the dynamical drift term for use in + # `check_timestep`. + drift_squared :: OnlineStatistics{Float64} + function ImplicitMidpoint(Δt; atol=1e-12) Δt <= 0 && error("Select positive Δt") - return new(Δt, atol) + return new(Δt, atol, OnlineStatistics{Float64}()) end end +function reset_statistics(integrator::Union{Langevin, ImplicitMidpoint}) + integrator.drift_squared = OnlineStatistics{Float64}() +end + +""" + check_timestep(integrator; tol) + +Prints a suggested timestep scale for the relative error tolerance `tol`, based +on statistics collected from previous time integration. A reasonable tolerance +might be of order `tol = 1e-2`, but should be selected on a case-by-case basis. +Note that `tol` is interpreted as the error for the deterministic part of a +single integration timestep, which scales like `Δt²`. Because the stochastic +Heun integration scheme is weakly convergent of order-1, errors in the estimates +of averaged observables may instead scale like ` √tol ~ Δt`. It is the +responsibility of the user to choose `tol` accordingly. + +The suggested timestep will be inversely proportional to the largest energy +scale that plays a role in the dynamics. If there is a strong Langvin coupling +``λ`` to the thermal bath, it will effectively rescale the energy scale. + +Timestep statistics stored in `integrator` will be reset after calling this +function. +""" +function check_timestep(langevin::Langevin; tol, reset=true) + (; Δt, λ, kT) = langevin + + drift_magnitude = sqrt(mean(langevin.drift_squared)) + if isnan(drift_magnitude) + error("Run trajectory first to collect statistics") + end + + # Heun integration is second-order accurate, so the local error from each + # deterministic timestep scales as dθ². Angular displacement per timestep dθ + # scales like dt |drift|, yielding err1 ~ (dt |drift|)^2 + # + # The "relevant" error introduced by thermal noise should be quantified in + # terms of the effect on statistical observables. To avoid subtleties in + # defining this error, we instead naïvely assume it 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|² + c₂² λ² kT²)) + # + # for some empirical constants c₁ and c₂. + + c1 = 1.0 + c2 = 1.0 + Δt_bound = sqrt(tol / ((c1*drift_magnitude)^2 + (c2*λ*kT)^2)) + + if Δt_bound/2 < Δt < 2Δt_bound + opinion = "reasonable." + elseif Δt_bound/5 < Δt < 5Δt_bound + opinion = "_marginal_." + elseif Δt <= Δt_bound/5 + opinion = "SMALL!" + elseif Δt >= 5Δt_bound + opinion = "LARGE!" + end + + λkTstr, tolstr, Δtstr, boundstr = number_to_simple_string.((λ*kT, tol, Δt, Δt_bound); digits=4) + + println("Suggest Δt ~ $boundstr for tol = $tolstr.") + println("Current Δt = $Δtstr is $opinion") + if c2*λ*kT > c1*drift_magnitude + println("Thermal noise contributes significantly, with λ*kT = $λkTstr.") + end + + reset && reset_statistics(langevin) + + return +end + + ################################################################################ # Dipole integration ################################################################################ @@ -122,7 +208,19 @@ 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 + + # Collect statistics about the magnitude of the drift term squared. For the + # energy-conserving part, it is important to use |∇E|² rather than |s×∇E|², + # because the former gives direct information about the frequency of + # oscillation, while the latter may artificially vanish (consider, e.g., a + # spin nearly aligned with an external field). + for i in eachsite(sys) + s0 = sys.κs[i] # == norm(s[i]) + ∇E² = ∇E[i]' * ∇E[i] + accum!(integrator.drift_squared, (1 + (s0*λ)^2) * ∇E²) + end + + return end # The spherical midpoint method, Phys. Rev. E 89, 061301(R) (2014) @@ -178,8 +276,9 @@ 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 @@ -200,6 +299,15 @@ function step!(sys::System{N}, integrator::Langevin) where N # Coordinate dipole data @. sys.dipoles = expected_spin(Z) + + # Collect statistics about the magnitude of the drift term squared. Here, + # |HZ|² plays the role of |∇E|² in the dipole case. + for i in eachsite(sys) + HZ² = HZ[i]' * HZ[i] + accum!(integrator.drift_squared, (1 + integrator.λ^2) * HZ²) + end + + return end function rhs_langevin!(ΔZ::Array{CVec{N}, 4}, Z::Array{CVec{N}, 4}, ξ::Array{CVec{N}, 4}, 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..624de97ed 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -39,6 +39,8 @@ const HermitianC64 = Hermitian{ComplexF64, Matrix{ComplexF64}} end end +include("OnlineStatistics.jl") + include("Operators/Spin.jl") include("Operators/Rotation.jl") include("Operators/Stevens.jl") @@ -82,7 +84,7 @@ include("Reshaping.jl") export reshape_supercell, resize_supercell, repeat_periodically include("Integrators.jl") -export Langevin, ImplicitMidpoint, step! +export Langevin, ImplicitMidpoint, step!, check_timestep include("Optimization.jl") export minimize_energy! From ef96abb855c868243e56f512edee1b2b42c096d4 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Tue, 16 Jan 2024 10:12:57 -0700 Subject: [PATCH 2/6] Stateless version --- docs/src/versions.md | 4 +- examples/03_LLD_CoRh2O4.jl | 74 +++++++++------- examples/04_GSD_FeI2.jl | 94 +++++++++++--------- src/Integrators.jl | 175 +++++++++++++++++++------------------ src/Sunny.jl | 2 +- 5 files changed, 190 insertions(+), 159 deletions(-) diff --git a/docs/src/versions.md b/docs/src/versions.md index db62a3061..5b55dcdb7 100644 --- a/docs/src/versions.md +++ b/docs/src/versions.md @@ -3,8 +3,8 @@ ## v0.5.9 (In development) -* New function [`check_timestep`](@ref) to verify accuracy of spin dynamics - integration. +* New functions [`suggest_timestep`](@ref) and [`check_timestep`](@ref) to + assist in the accurate 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 b7e3054eb..1413fbcd8 100644 --- a/examples/03_LLD_CoRh2O4.jl +++ b/examples/03_LLD_CoRh2O4.jl @@ -25,35 +25,44 @@ 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 + +# Use [`suggest_timestep`](@ref) to obtain a reasonable integration timestep. +# The spin configuration in `sys` should ideally be equilibrated for the target +# temperature `kT`, but the current, energy-minimized configuration will also +# work reasonably well. Usually `tol=1e-2` is good tolerance to numerical error. + +suggest_timestep(sys; tol=1e-2, λ, kT) + +# We now have all data needed to construct a Langevin integrator. + +Δt = 0.025 langevin = Langevin(Δt; λ, kT); -# Because the magnetic order has been initialized correctly, relatively few -# additional Langevin time-steps are required to reach thermal equilibrium. +# 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)")) -# After running a Langevin trajectory, it is a good practice to call -# [`check_timestep`](@ref). +# Calling [`check_timestep`](@ref) indicates that our choice of `Δt` was a +# little smaller than necessary; increasing it will improve efficiency. + +check_timestep(sys, langevin; tol=1e-2) +langevin.Δt = 0.042 -check_timestep(langevin; tol=1e-2) +# 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. @@ -63,8 +72,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` @@ -101,7 +110,7 @@ formfactors = [FormFactor("Co2")] instant_formula = intensity_formula(sc, :perp; formfactors) iq = instant_intensities_interpolated(sc, qs, instant_formula); -# Plot the resulting intensity data ``I(𝐪)``. The color scale is clipped to 50% +# Plot the resulting intensity data ``𝒮(𝐪)``. The color scale is clipped to 50% # of the maximum intensity. heatmap(q1s, q2s, iq; @@ -117,18 +126,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 @@ -137,7 +149,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 117324520..51a5c5778 100644 --- a/examples/04_GSD_FeI2.jl +++ b/examples/04_GSD_FeI2.jl @@ -74,46 +74,54 @@ 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 +# As previously observed, direct energy minimization is susceptible to trapping +# in a local energy minimum. + +randomize_spins!(sys) +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 + +# Use [`suggest_timestep`](@ref) to obtain a reasonable integration timestep. It +# is important that the system has already been initialized to a low-energy +# configuration. Usually `tol=1e-2` is good tolerance to numerical error. + +suggest_timestep(sys; tol=1e-2, λ, kT) + +# This information is sufficient to define the Langevin integrator. + +Δt = 0.027 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. +# Langevin dynamics can be used to search for a magnetically ordered state. This +# works well because the temperature `kT = 0.2` has been carefully selected. It +# is below the ordering temperature, but large enough that the dynamical +# trajectory can overcome local energy barriers and annihilate defects. -randomize_spins!(sys) -for _ in 1:20_000 +for _ in 1:10_000 step!(sys, langevin) end -# After running a Langevin trajectory, it is a good practice to call -# [`check_timestep`](@ref). +# Calling [`check_timestep`](@ref) shows that thermalization has not +# substantially altered the suggested `Δt`. -check_timestep(langevin; tol=1e-2) +check_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$ # @@ -136,7 +144,10 @@ for _ in 1:10_000 step!(sys_large, langevin) end -check_timestep(langevin; tol=1e-2) +# With this increase in temperature, the suggested timestep has increased slightly. + +check_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 @@ -148,18 +159,21 @@ check_timestep(langevin; tol=1e-2) # 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 @@ -169,7 +183,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: @@ -202,9 +216,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 9b0128378..c791a3987 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -53,13 +53,9 @@ mutable struct Langevin λ :: Float64 kT :: Float64 - # Averages the square of the dynamical drift term for use in - # `check_timestep`. - drift_squared :: OnlineStatistics{Float64} - function Langevin(Δt; λ, kT) Δt <= 0 && error("Select positive Δt") - return new(Δt, λ, kT, OnlineStatistics{Float64}()) + return new(Δt, λ, kT) end end @@ -84,92 +80,119 @@ mutable struct ImplicitMidpoint Δt :: Float64 atol :: Float64 - # Averages the square of the dynamical drift term for use in - # `check_timestep`. - drift_squared :: OnlineStatistics{Float64} - function ImplicitMidpoint(Δt; atol=1e-12) Δt <= 0 && error("Select positive Δt") - return new(Δt, atol, OnlineStatistics{Float64}()) + return new(Δt, atol) end end - -function reset_statistics(integrator::Union{Langevin, ImplicitMidpoint}) - integrator.drift_squared = OnlineStatistics{Float64}() -end - """ - check_timestep(integrator; tol) - -Prints a suggested timestep scale for the relative error tolerance `tol`, based -on statistics collected from previous time integration. A reasonable tolerance -might be of order `tol = 1e-2`, but should be selected on a case-by-case basis. -Note that `tol` is interpreted as the error for the deterministic part of a -single integration timestep, which scales like `Δt²`. Because the stochastic -Heun integration scheme is weakly convergent of order-1, errors in the estimates -of averaged observables may instead scale like ` √tol ~ Δt`. It is the -responsibility of the user to choose `tol` accordingly. - -The suggested timestep will be inversely proportional to the largest energy -scale that plays a role in the dynamics. If there is a strong Langvin coupling -``λ`` to the thermal bath, it will effectively rescale the energy scale. - -Timestep statistics stored in `integrator` will be reset after calling this -function. + suggest_timestep(sys; tol, λ=0, kT=0) + +Suggests a timestep for the numerical integration of spin dynamics according to +a given error tolerance `tol`. 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 is appropriate for numerical integration schemes that are +second-order accurate. + +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). + +The optional parameters ``λ`` and ``kT`` will tighten the timestep bound by +accounting for [`Langevin`](@ref) coupling to a thermal bath. If the damping +magnitude ``λ`` is order 1 or larger, it will effectively rescale the suggested +timestep by ``1/λ``. The target temperature ``kT`` controls the magnitude of the +noise term, and is treated as an additional energy scale for purposes of +estimating ``Δt``. + +When the dynamics includes a noise term, quantifying numerical error becomes +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 adjusted accordingly. """ -function check_timestep(langevin::Langevin; tol, reset=true) - (; Δt, λ, kT) = langevin +function suggest_timestep(sys::System{N}; tol, λ=0, kT=0) where N + Δt_bound = suggest_timestep_aux(sys; tol, λ, kT) + Δt_str, tol_str = number_to_simple_string.((Δt_bound, tol); digits=4) + println("Suggested timestep Δt ≲ $Δt_str for tol = $tol_str at the given spin configuration.") +end - drift_magnitude = sqrt(mean(langevin.drift_squared)) - if isnan(drift_magnitude) - error("Run trajectory first to collect statistics") +function suggest_timestep_aux(sys::System{N}; tol, λ=0, kT=0) 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 - # Heun integration is second-order accurate, so the local error from each - # deterministic timestep scales as dθ². Angular displacement per timestep dθ - # scales like dt |drift|, yielding err1 ~ (dt |drift|)^2 + # `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 # - # The "relevant" error introduced by thermal noise should be quantified in - # terms of the effect on statistical observables. To avoid subtleties in - # defining this error, we instead naïvely assume it 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)². + # 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|² + c₂² λ² kT²)) + # 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_magnitude)^2 + (c2*λ*kT)^2)) - - if Δt_bound/2 < Δt < 2Δt_bound - opinion = "reasonable." - elseif Δt_bound/5 < Δt < 5Δt_bound - opinion = "_marginal_." - elseif Δt <= Δt_bound/5 - opinion = "SMALL!" - elseif Δt >= 5Δt_bound - opinion = "LARGE!" - end + return sqrt(tol / ((c1*drift_rms)^2 + (c2*λ*kT)^2)) +end - λkTstr, tolstr, Δtstr, boundstr = number_to_simple_string.((λ*kT, tol, Δt, Δt_bound); digits=4) +""" + check_timestep(sys, integrator; tol) - println("Suggest Δt ~ $boundstr for tol = $tolstr.") - println("Current Δt = $Δtstr is $opinion") - if c2*λ*kT > c1*drift_magnitude - println("Thermal noise contributes significantly, with λ*kT = $λkTstr.") - end +Compares the the timestep in `integrator` to that of [`suggest_timestep`](@ref), +and prints an informative message. +""" +function check_timestep(sys::System{N}, langevin::Langevin; tol) where N + (; Δt, λ, kT) = langevin - reset && reset_statistics(langevin) + Δt_bound = suggest_timestep_aux(sys; tol, λ, kT) + Δtstr, boundstr = number_to_simple_string.((Δt, Δt_bound); digits=4) - return + print("Current Δt = $Δtstr vs suggested Δt = $(boundstr).") + if Δt <= Δt_bound/2 + println("\nTimestep looks small! Increasing it will make the simulation faster.") + elseif Δt >= 2Δt_bound + println("\nWARNING: Timestep looks large! Decreasing it will improve accuracy.") + else + println(" Agreement seems reasonable.") + end end @@ -209,17 +232,6 @@ function step!(sys::System{0}, integrator::Langevin) @. s = s + 0.5 * Δt * (f₁ + rhs_dipole(s₁, -∇E, λ)) + 0.5 * √Δt * (r₁ + rhs_dipole(s₁, ξ)) @. s = normalize_dipole(s, sys.κs) - # Collect statistics about the magnitude of the drift term squared. For the - # energy-conserving part, it is important to use |∇E|² rather than |s×∇E|², - # because the former gives direct information about the frequency of - # oscillation, while the latter may artificially vanish (consider, e.g., a - # spin nearly aligned with an external field). - for i in eachsite(sys) - s0 = sys.κs[i] # == norm(s[i]) - ∇E² = ∇E[i]' * ∇E[i] - accum!(integrator.drift_squared, (1 + (s0*λ)^2) * ∇E²) - end - return end @@ -300,13 +312,6 @@ function step!(sys::System{N}, integrator::Langevin) where N # Coordinate dipole data @. sys.dipoles = expected_spin(Z) - # Collect statistics about the magnitude of the drift term squared. Here, - # |HZ|² plays the role of |∇E|² in the dipole case. - for i in eachsite(sys) - HZ² = HZ[i]' * HZ[i] - accum!(integrator.drift_squared, (1 + integrator.λ^2) * HZ²) - end - return end diff --git a/src/Sunny.jl b/src/Sunny.jl index 624de97ed..e83862b21 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -84,7 +84,7 @@ include("Reshaping.jl") export reshape_supercell, resize_supercell, repeat_periodically include("Integrators.jl") -export Langevin, ImplicitMidpoint, step!, check_timestep +export Langevin, ImplicitMidpoint, step!, suggest_timestep, check_timestep include("Optimization.jl") export minimize_energy! From 5a06a7bb1cbb6b3a6b3e9b1f17a10be6b21b8998 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Tue, 16 Jan 2024 10:12:57 -0700 Subject: [PATCH 3/6] Simplify interface --- docs/src/versions.md | 4 +- examples/03_LLD_CoRh2O4.jl | 24 ++++----- examples/04_GSD_FeI2.jl | 35 ++++++------ src/Integrators.jl | 107 +++++++++++++++++++++---------------- src/Sunny.jl | 2 +- 5 files changed, 91 insertions(+), 81 deletions(-) diff --git a/docs/src/versions.md b/docs/src/versions.md index 5b55dcdb7..6057f12e1 100644 --- a/docs/src/versions.md +++ b/docs/src/versions.md @@ -3,8 +3,8 @@ ## v0.5.9 (In development) -* New functions [`suggest_timestep`](@ref) and [`check_timestep`](@ref) to - assist in the accurate simulation of spin dynamics. +* 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 1413fbcd8..6ec241cbe 100644 --- a/examples/03_LLD_CoRh2O4.jl +++ b/examples/03_LLD_CoRh2O4.jl @@ -31,18 +31,15 @@ sys = resize_supercell(sys, (10, 10, 10)) λ = 0.2 # Magnitude of damping (dimensionless) kT = 16 * meV_per_K # 16K, a temperature slightly below ordering +langevin = Langevin(; λ, kT) -# Use [`suggest_timestep`](@ref) to obtain a reasonable integration timestep. -# The spin configuration in `sys` should ideally be equilibrated for the target -# temperature `kT`, but the current, energy-minimized configuration will also -# work reasonably well. Usually `tol=1e-2` is good tolerance to numerical error. +# 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; tol=1e-2, λ, kT) - -# We now have all data needed to construct a Langevin integrator. - -Δt = 0.025 -langevin = Langevin(Δt; λ, kT); +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. @@ -53,10 +50,11 @@ for _ in 1:1000 push!(energies, energy_per_site(sys)) end -# Calling [`check_timestep`](@ref) indicates that our choice of `Δt` was a -# little smaller than necessary; increasing it will improve efficiency. +# 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. -check_timestep(sys, langevin; tol=1e-2) +suggest_timestep(sys, langevin; tol=1e-2) langevin.Δt = 0.042 # The energy per site has converged, which suggests that the system has reached diff --git a/examples/04_GSD_FeI2.jl b/examples/04_GSD_FeI2.jl index 51a5c5778..b69c5fdb3 100644 --- a/examples/04_GSD_FeI2.jl +++ b/examples/04_GSD_FeI2.jl @@ -74,8 +74,9 @@ sys # ## Finding a ground state -# As previously observed, direct energy minimization is susceptible to trapping -# in a local energy minimum. +# 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) minimize_energy!(sys) @@ -89,31 +90,29 @@ plot_spins(sys; color=[s[3] for s in sys.dipoles]) λ = 0.2 # Dimensionless damping time-scale kT = 0.2 # Temperature in meV +langevin = Langevin(; λ, kT) -# Use [`suggest_timestep`](@ref) to obtain a reasonable integration timestep. It -# is important that the system has already been initialized to a low-energy -# configuration. Usually `tol=1e-2` is good tolerance to numerical error. +# 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; tol=1e-2, λ, kT) +suggest_timestep(sys, langevin; tol=1e-2) +langevin.Δt = 0.027 -# This information is sufficient to define the Langevin integrator. - -Δt = 0.027 -langevin = Langevin(Δt; kT, λ); - -# Langevin dynamics can be used to search for a magnetically ordered state. This -# works well because the temperature `kT = 0.2` has been carefully selected. It -# is below the ordering temperature, but large enough that the dynamical -# trajectory can overcome local energy barriers and annihilate defects. +# 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 [`check_timestep`](@ref) shows that thermalization has not +# Calling [`suggest_timestep`](@ref) shows that thermalization has not # substantially altered the suggested `Δt`. -check_timestep(sys, langevin; tol=1e-2) +suggest_timestep(sys, langevin; tol=1e-2) # Although thermal fluctuations are present, the correct antiferromagnetic order # (2 up, 2 down) has been found. @@ -146,7 +145,7 @@ end # With this increase in temperature, the suggested timestep has increased slightly. -check_timestep(sys_large, langevin; tol=1e-2) +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 diff --git a/src/Integrators.jl b/src/Integrators.jl index c791a3987..90e0233c5 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 @@ -58,6 +60,8 @@ mutable struct Langevin return new(Δt, λ, kT) 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,42 +91,47 @@ 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; tol, λ=0, kT=0) + suggest_timestep(sys, integrator; tol) Suggests a timestep for the numerical integration of spin dynamics according to -a given error tolerance `tol`. 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 is appropriate for numerical integration schemes that are -second-order accurate. +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). -The optional parameters ``λ`` and ``kT`` will tighten the timestep bound by -accounting for [`Langevin`](@ref) coupling to a thermal bath. If the damping -magnitude ``λ`` is order 1 or larger, it will effectively rescale the suggested -timestep by ``1/λ``. The target temperature ``kT`` controls the magnitude of the -noise term, and is treated as an additional energy scale for purposes of -estimating ``Δt``. - -When the dynamics includes a noise term, quantifying numerical error becomes +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 adjusted accordingly. +true numerical error, and should be selected with this in mind. """ -function suggest_timestep(sys::System{N}; tol, λ=0, kT=0) where N - Δt_bound = suggest_timestep_aux(sys; tol, λ, kT) - Δt_str, tol_str = number_to_simple_string.((Δt_bound, tol); digits=4) - println("Suggested timestep Δt ≲ $Δt_str for tol = $tol_str at the given spin configuration.") +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, λ=0, kT=0) where N +function suggest_timestep_aux(sys::System{N}; tol, Δt, λ, kT) where N acc = 0.0 if N == 0 ∇Es, = get_dipole_buffers(sys, 1) @@ -170,28 +181,22 @@ function suggest_timestep_aux(sys::System{N}; tol, λ=0, kT=0) where N c1 = 1.0 c2 = 1.0 - return sqrt(tol / ((c1*drift_rms)^2 + (c2*λ*kT)^2)) -end - -""" - check_timestep(sys, integrator; tol) - -Compares the the timestep in `integrator` to that of [`suggest_timestep`](@ref), -and prints an informative message. -""" -function check_timestep(sys::System{N}, langevin::Langevin; tol) where N - (; Δt, λ, kT) = langevin - - Δt_bound = suggest_timestep_aux(sys; tol, λ, kT) - Δtstr, boundstr = number_to_simple_string.((Δt, Δt_bound); digits=4) - - print("Current Δt = $Δtstr vs suggested Δt = $(boundstr).") - if Δt <= Δt_bound/2 - println("\nTimestep looks small! Increasing it will make the simulation faster.") - elseif Δt >= 2Δt_bound - println("\nWARNING: Timestep looks large! Decreasing it will improve accuracy.") - else - println(" Agreement seems reasonable.") + Δ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 @@ -214,6 +219,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 @@ -242,6 +249,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 @@ -294,6 +303,8 @@ end 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 @@ -333,6 +344,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/Sunny.jl b/src/Sunny.jl index e83862b21..fa8158b2d 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -84,7 +84,7 @@ include("Reshaping.jl") export reshape_supercell, resize_supercell, repeat_periodically include("Integrators.jl") -export Langevin, ImplicitMidpoint, step!, suggest_timestep, check_timestep +export Langevin, ImplicitMidpoint, step!, suggest_timestep include("Optimization.jl") export minimize_energy! From 6f6577fe9862c7d664f166c3c23b73b0af1b25a6 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Tue, 16 Jan 2024 10:12:57 -0700 Subject: [PATCH 4/6] Don't suppress print statements --- examples/03_LLD_CoRh2O4.jl | 4 ++-- examples/04_GSD_FeI2.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/03_LLD_CoRh2O4.jl b/examples/03_LLD_CoRh2O4.jl index 6ec241cbe..28d42db1a 100644 --- a/examples/03_LLD_CoRh2O4.jl +++ b/examples/03_LLD_CoRh2O4.jl @@ -39,7 +39,7 @@ langevin = Langevin(; λ, kT) # configuration will also work reasonably well. suggest_timestep(sys, langevin; tol=1e-2) -langevin.Δt = 0.025 +langevin.Δt = 0.025; # Now run a dynamical trajectory to sample spin configurations. Keep track of # the energy per site at each time step. @@ -55,7 +55,7 @@ end # simulations faster. suggest_timestep(sys, langevin; tol=1e-2) -langevin.Δt = 0.042 +langevin.Δt = 0.042; # The energy per site has converged, which suggests that the system has reached # thermal equilibrium. diff --git a/examples/04_GSD_FeI2.jl b/examples/04_GSD_FeI2.jl index b69c5fdb3..ae377a44b 100644 --- a/examples/04_GSD_FeI2.jl +++ b/examples/04_GSD_FeI2.jl @@ -98,7 +98,7 @@ langevin = Langevin(; λ, kT) # configuration will also work reasonably well. suggest_timestep(sys, langevin; tol=1e-2) -langevin.Δt = 0.027 +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 @@ -146,7 +146,7 @@ end # With this increase in temperature, the suggested timestep has increased slightly. suggest_timestep(sys_large, langevin; tol=1e-2) -langevin.Δt = 0.040 +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 From a22591fad82577c2a9de1e39c182ddb5f8b88411 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Tue, 16 Jan 2024 10:12:57 -0700 Subject: [PATCH 5/6] Better show() on integrators. --- src/Integrators.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/Integrators.jl b/src/Integrators.jl index 90e0233c5..636ac5e33 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -201,6 +201,19 @@ function suggest_timestep_aux(sys::System{N}; tol, Δt, λ, kT) where N end +function Base.show(io::IO, integrator::Langevin) + (; Δt, λ, kT) = integrator + Δt = isnan(integrator.Δt) ? "" : repr(Δt) + println(io, "Langevin($Δt; λ=$λ, kT=$kT)") +end + +function Base.show(io::IO, integrator::ImplicitMidpoint) + (; Δt, atol) = integrator + Δt = isnan(integrator.Δt) ? "" : repr(Δt) + println(io, "ImplicitMidpoint($Δt; atol=$atol)") +end + + ################################################################################ # Dipole integration ################################################################################ From a2f474e81929acb8bdfa65c81935d22dabc5d600 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Tue, 16 Jan 2024 10:44:49 -0700 Subject: [PATCH 6/6] Tweaks --- examples/03_LLD_CoRh2O4.jl | 2 +- src/Sunny.jl | 2 -- 2 files changed, 1 insertion(+), 3 deletions(-) diff --git a/examples/03_LLD_CoRh2O4.jl b/examples/03_LLD_CoRh2O4.jl index 28d42db1a..e70d479b3 100644 --- a/examples/03_LLD_CoRh2O4.jl +++ b/examples/03_LLD_CoRh2O4.jl @@ -108,7 +108,7 @@ formfactors = [FormFactor("Co2")] instant_formula = intensity_formula(sc, :perp; formfactors) iq = instant_intensities_interpolated(sc, qs, instant_formula); -# Plot the resulting intensity data ``𝒮(𝐪)``. The color scale is clipped to 50% +# Plot the resulting intensity data ``I(𝐪)``. The color scale is clipped to 50% # of the maximum intensity. heatmap(q1s, q2s, iq; diff --git a/src/Sunny.jl b/src/Sunny.jl index fa8158b2d..1da32746a 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -39,8 +39,6 @@ const HermitianC64 = Hermitian{ComplexF64, Matrix{ComplexF64}} end end -include("OnlineStatistics.jl") - include("Operators/Spin.jl") include("Operators/Rotation.jl") include("Operators/Stevens.jl")