Add copy_spins!. Expose jitter argument in minimize_energy! (#465)

Also includes various updates to the documentation. In particular, `minimize_energy!` suggests a code pattern to search for the ground state.

Author: kbarros

docs/src/library.md

@@ -31,6 +31,7 @@ Units
 add_sample!
 clone_correlations
 clone_system
+copy_spins!
 dispersion
 dmvec
 domain_average

docs/src/versions.md

@@ -3,7 +3,9 @@
 ## v0.8.1
 (In progress)
 
-* Fix plotting error with SCGA-calculated intensities ([#464](@ref)).
+* Fix plotting error with [`SCGA`](@ref)-calculated intensities ([#464](@ref)).
+* Add [`copy_spins!`](@ref). Expose `jitter` parameter in
+  [`minimize_energy!`](@ref) ([#465](@ref)).
 
 ## v0.8.0
 (Nov 30, 2025)

examples/03_LSWT_SU3_FeI2.jl

@@ -142,15 +142,15 @@ set_onsite_coupling!(sys, S -> -D*S[3]^2, 1)
 
 # ### Finding the ground state
 
-# This model has been fitted so that energy minimization yields the physically
-# correct ground state. Knowing this, we could set the magnetic configuration
-# manually by calling [`set_dipole!`](@ref) on each site in the system. Another
-# approach, as we will demonstrate, is to search for the ground-state via
-# [`minimize_energy!`](@ref).
+# The model parameters have already been fitted so that energy minimization
+# yields the physically correct ground state. Knowing this, one could manually
+# set the magnetic configuration by calling [`set_dipole!`](@ref) at each site.
+# Alternatively, Sunny provides tools to interactively search for the ground
+# state.
 #
-# To reduce bias in the search, use [`resize_supercell`](@ref) to create a
-# relatively large system of 4×4×4 chemical cells. Randomize all spins
-# (represented as SU(3) coherent states) and minimize the energy.
+# Use [`resize_supercell`](@ref) to create a relatively large system of 4×4×4
+# chemical cells. Call [`randomize_spins!`](@ref) and [`minimize_energy!`](@ref)
+# in sequence.
 
 sys = resize_supercell(sys, (4, 4, 4))
 randomize_spins!(sys)
@@ -161,36 +161,36 @@ minimize_energy!(sys)
 
 plot_spins(sys; color=[S[3] for S in sys.dipoles])
 
-# One could precisely quantify the Fourier-space static structure factor
-# ``\mathcal{S}(𝐪)`` of this spin configuration using
-# [`SampledCorrelationsStatic`](@ref). For the present purposes, however, it is
-# most convenient to use [`print_wrapped_intensities`](@ref), which effectively
-# averages ``\mathcal{S}(𝐪)`` over all Brillouin zones.
+# Finding the true ground state can be a challenging task. If the system is not
+# too large, a good strategy is repeated randomization and then energy
+# minimization. This particular system might require about 30 minimization runs
+# to find the defect-free ground state.
+#
+# Another tool is [`print_wrapped_intensities`](@ref). It reports weights
+# analogous to the static structure factor ``\mathcal{S}(𝐪)``, but without
+# accounting for phase interference between magnetic sublattices. To calculate
+# the true ``\mathcal{S}(𝐪)``, use instead [`SampledCorrelationsStatic`](@ref).
 
 print_wrapped_intensities(sys)
 
-# The known zero-field energy-minimizing magnetic structure of FeI₂ is a two-up,
-# two-down order. It can be described as a generalized spiral with a single
-# propagation wavevector ``𝐤``. Rotational symmetry allows for three equivalent
-# orientations: ``±𝐤 = [0, -1/4, 1/4]``, ``[1/4, 0, 1/4]``, or
-# ``[-1/4,1/4,1/4]``. Energy minimization of large systems can easily get
-# trapped in a metastable state with competing domains and defects. Nonetheless,
-# `print_wrapped_intensities` shows that a non-trivial fraction of the total
-# intensity is consistent with known ordering wavevectors.
+# The correct ground state for FeI₂ is known to be a generalized spiral with one
+# of three propagation wavevectors, ``[0, -1/4, 1/4]`` or ``[1/4, 0, 1/4]`` or
+# ``[-1/4, 1/4, 1/4]``, which are equivalent under 120° rotations. The result of
+# `print_wrapped_intensities` hints at this spiral phase.
 #
-# Let's break the three-fold symmetry by hand. The function
-# [`suggest_magnetic_supercell`](@ref) takes any number of ``𝐤`` modes and
-# suggests a magnetic cell shape that is commensurate with all of them.
+# Let's break the 3-fold symmetry by hand. The function
+# [`suggest_magnetic_supercell`](@ref) takes any number of propagation
+# wavevectors and suggests a commensurate magnetic cell.
 
 suggest_magnetic_supercell([[0, -1/4, 1/4]])
 
-# Using the minimal system returned by [`reshape_supercell`](@ref), it is now
-# easy to find the ground state. Plot the system again, now including "ghost"
-# spins out to 12Å, to verify that the magnetic order is consistent with FeI₂.
+# After calling [`reshape_supercell`](@ref), it becomes easy to find the ground
+# state. Plot the system again, now including "ghost" spins out to 12Å. The
+# correct two-up, two-down magnetic order is visually apparent.
 
-sys_min = reshape_supercell(sys, [1 0 0; 0 2 1; 0 -2 1])
+sys_min = reshape_supercell(sys, [1 0 0; 0 1 -2; 0 1 2])
 randomize_spins!(sys_min)
-minimize_energy!(sys_min);
+minimize_energy!(sys_min)
 plot_spins(sys_min; color=[S[3] for S in sys_min.dipoles], ghost_radius=12)
 
 # ### Spin wave theory

src/Integrators.jl

@@ -142,27 +142,24 @@ 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 ``dt`` will be inversely proportional
-to the magnitude of the effective field ``|dE/d𝐒|`` arising from the current
-spin configuration in `sys`. The recommended timestep ``dt`` 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).
-
-Large `damping` magnitude or target temperature `kT` will tighten the timestep
-bound. If `damping` exceeds 1, it will rescale the suggested timestep by an
-approximate the factor ``1/damping``. If `kT` is the largest energy scale, then
-the suggested timestep will scale like `1/(damping*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 `dt`. 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.
+Suggests a timestep `dt` for spin dynamics simulation at a given error tolerance
+`tol`. The `integrator` should be [`Langevin`](@ref) or
+[`ImplicitMidpoint`](@ref). Ideally, the spin configuration in `sys` would be
+equilibrated to the target thermodynamic conditions. In practice, a
+configuration obtained from [`minimize_energy!`](@ref) should give a reasonable,
+if conservative, `dt` suggestion.
+
+The suggested `dt` scales like `√tol`, consistent with a second order
+integration scheme. In most cases, `dt` will also be inversely proportional to
+the characteristic magnitude of the energy gradient, ``∂E/∂𝐒_i``. If the
+Langevin noise magnitude ``λ k_B T`` dominates, however, then its inverse will
+limit the `dt` scale.
+
+Analysis of error in Langevin dynamics can be subtle. Sunny uses the stochastic
+Heun scheme, which has a weak convergence rate of order 1. This means that
+errors in certain statistical observables may scale like `dt` rather than
+`dt^2`. In such cases, the `tol` parameter controls the _square_ of the
+numerical error, and can be tightened appropriately.
 """
 function suggest_timestep(sys::System, integrator::Union{Langevin, ImplicitMidpoint}; tol)
     (; dt) = integrator

src/MagneticOrdering.jl

@@ -1,19 +1,21 @@
 """
     print_wrapped_intensities(sys::System; nmax=10)
 
-For Bravais lattices: Prints up to `nmax` wavevectors according to their
-instantaneous (static) structure factor intensities, listed in descending order.
-For non-Bravais lattices: Performs the same analysis for each spin sublattice
-independently; the output weights are naïvely averaged over sublattices, without
-incorporating phase shift information. This procedure therefore wraps all
-wavevectors into the first Brillouin zone. Each wavevector coordinate is given
-between ``-1/2`` and ``1/2`` in reciprocal lattice units (RLU).  The output from
-this function will typically be used as input to
+Prints up to `nmax` wavevectors ``𝐪`` and their "wrapped" static structure
+factor weights. Each ``𝐪`` is exactly commensurate with the system volume and
+has components between ``-1/2`` and ``1/2`` in reciprocal lattice units (RLU).
+The output from this function will typically be used as input to
 [`suggest_magnetic_supercell`](@ref).
 
-Because this function does not incorporate phase information in its averaging
-over sublattices, the printed weights are not directly comparable with
-experiment. For that purpose, use [`SampledCorrelationsStatic`](@ref) instead.
+For simplicity, phase interference between sublattices is neglected. The
+reported weights are the sum of static structure factors
+``\\mathcal{S}_{jj}(𝐪)`` calculated independently for each sublattice ``j`` of
+the chemical cell. This is mathematically equivalent to averaging
+``\\mathcal{S}(𝐪)`` over all cells of the infinite reciprocal lattice. It is in
+this sense that the intensities are "wrapped" into the first reciprocal cell.
+
+To calculate the true ``\\mathcal{S}(𝐪)`` as an experimental observable, use
+[`SampledCorrelationsStatic`](@ref) instead.
 """
 function print_wrapped_intensities(sys::System{N}; nmax=10) where N
     sys.crystal == orig_crystal(sys) || error("Cannot perform this analysis on reshaped system.")

src/Optimization.jl

@@ -1,17 +1,22 @@
 struct MinimizationResult
     converged :: Bool
-    iterations :: Int
+    niters :: Int
+    energy :: Float64
     data :: Optim.OptimizationResults
+
+    function MinimizationResult(niters, res)
+        return new(Optim.converged(res), niters, Optim.minimum(res), res)
+    end
 end
 
 function Base.show(io::IO, ::MIME"text/plain", res::MinimizationResult)
-    (; converged, iterations, data) = res
+    (; converged, niters, data) = res
     Δx = number_to_simple_string(Optim.x_abschange(data); digits=3)
     g_res = number_to_simple_string(Optim.g_residual(data); digits=3)
     if converged
-        print(io, "Converged in $iterations iterations")
+        print(io, "Converged in $niters iterations")
     else
-        print(io, "Non-converged after $iterations iterations: |Δx|=$Δx, |∂E/∂x|=$g_res")
+        print(io, "Non-converged after $niters iterations: |Δx|=$Δx, |∂E/∂x|=$g_res")
     end
 end
 
@@ -52,34 +57,52 @@ function Optim.project_tangent!(sm::SpinManifold, g, x)
 end
 
 function optimize_with_restarts(; calc_f, calc_g!, x, method, maxiters, options_args)
-    iters = 0
+    niters = 0
     while true
-        options = Optim.Options(; iterations=maxiters-iters, options_args...)
+        options = Optim.Options(; iterations=maxiters-niters, options_args...)
         res = Optim.optimize(calc_f, calc_g!, x, method, options)
         x = Optim.minimizer(res)
-        iters += Optim.iterations(res)
-        if Optim.converged(res) || iters >= maxiters
-            return (res, iters)
+        niters += Optim.iterations(res)
+        if Optim.converged(res) || niters >= maxiters
+            return (; res, niters)
         end
     end
 end
 
 
 """
-    minimize_energy!(sys::System; maxiters=1000, kwargs...)
-
-Optimizes the spin configuration in `sys` to minimize energy. A total of
-`maxiters` iterations will be attempted. Any remaining `kwargs` will be included
-in the `Options` constructor of the [Optim.jl
-package](https://github.com/JuliaNLSolvers/Optim.jl)
-
-Convergence status is stored in the field `ret.converged` of the return value
-`ret`. Additional optimization statistics are stored in the field `ret.data`.
+    minimize_energy!(sys::System; maxiters=1000, jitter=1e-8, kwargs...)
+
+Optimizes the spin configuration in `sys` to find a local minimum of the energy.
+Large magnetic cells will be slower to converge; increase `maxiters` as needed.
+Prior to optimization, each spin will be randomly perturbed with the
+dimensionless `jitter` magnitude, which can be useful to break accidental
+symmetries. Any remaining `kwargs` will be included in the `Options` constructor
+of the [Optim.jl package](https://github.com/JuliaNLSolvers/Optim.jl).
+
+Returns an object that can be inspected for optimization statistics.
+
+!!! tip "Escaping local minima"  
+    To search for the global energy minimum, a simple strategy is to repeatedly
+    call [`randomize_spins!`](@ref) and then `minimize_energy!`. The example
+    below keeps the minimum energy spin state found within 100 runs.
+
+    ```julia
+    tmp_sys = clone_system(sys)
+    for i in 1:100
+        randomize_spins!(sys)
+        minimize_energy!(sys)
+        if energy(sys) < energy(tmp_sys)
+            copy_spins!(tmp_sys, sys)
+        end
+    end
+    copy_spins!(sys, tmp_sys)
+    ```
 """
-function minimize_energy!(sys::System{N}; maxiters=1000, δ=1e-8, kwargs...) where N
+function minimize_energy!(sys::System{N}; maxiters=1000, jitter=1e-8, kwargs...) where N
     # Small perturbation to destabilize an accidental stationary point (local
     # maximum or saddle).
-    perturb_spins!(sys, δ)
+    perturb_spins!(sys, jitter)
 
     # Optimization variables are normalized spins or coherent states. In case of
     # a vacancy, use an arbitrary representative on the sphere: [1, 1, …] / √N.
@@ -132,10 +155,10 @@ function minimize_energy!(sys::System{N}; maxiters=1000, δ=1e-8, kwargs...) whe
     manifold = SpinManifold(iszero(N) ? 3 : N, length(eachsite(sys)))
     method = Optim.ConjugateGradient(; alphaguess=LineSearches.InitialHagerZhang(; αmax=10.0), manifold)
     options_args = (; g_abstol, x_abstol, x_reltol=NaN, f_reltol=NaN, f_abstol=NaN, kwargs...)
-    (res, iters) = optimize_with_restarts(; calc_f, calc_g!, x, method, maxiters, options_args)
+    (; res, niters) = optimize_with_restarts(; calc_f, calc_g!, x, method, maxiters, options_args)
 
     load_spins!(Optim.minimizer(res))
-    mr = MinimizationResult(Optim.converged(res), iters, res)
+    mr = MinimizationResult(niters, res)
     if !mr.converged
         @warn repr("text/plain", mr)
     end

src/Sunny.jl

@@ -63,7 +63,7 @@ include("System/OnsiteCoupling.jl")
 include("System/Ewald.jl")
 include("System/Interactions.jl")
 export Moment, System, Site, clone_system, eachsite, position_to_site, global_position,
-    magnetic_moment, set_coherent!, set_dipole!, polarize_spins!, randomize_spins!,
+    magnetic_moment, set_coherent!, set_dipole!, polarize_spins!, copy_spins!, randomize_spins!,
     set_spin_rescaling!, set_spin_s_at!, set_spin_rescaling_for_static_sum_rule!,
     energy, energy_per_site, spin_label, set_onsite_coupling!, set_pair_coupling!,
     set_exchange!, dmvec, enable_dipole_dipole!, set_field!, to_inhomogeneous, set_field_at!,

src/Symmetry/Printing.jl

@@ -146,9 +146,9 @@ end
 """
     print_bond(cryst::Crystal, b::Bond; b_ref=b)
 
-Prints symmetry information for bond `b`. An optional symmetry-equivalent
-reference bond `b_ref` can be provided to keep a consistent meaning of the free
-parameters `A`, `B`, etc.
+Prints symmetry-allowed interactions for the [`Bond`](@ref) `b`. An optional
+symmetry-equivalent reference bond `b_ref` can be provided to keep a consistent
+meaning of the free parameters `A`, `B`, etc.
 """
 function print_bond(cryst::Crystal, b::Bond; b_ref=b, io=stdout)
     # How many digits to use in printing coefficients
@@ -193,10 +193,12 @@ end
 """
     print_symmetry_table(cryst::Crystal, max_dist)
 
-Print symmetry information for all equivalence classes of sites and bonds, up to
-a maximum bond distance of `max_dist`. Equivalent to calling `print_bond(cryst,
-b)` for every bond `b` in `reference_bonds(cryst, max_dist)`, where
-`Bond(i, i, [0,0,0])` refers to a single site `i`.
+Prints the allowed interactions, as constrained by the spacegroup symmetries of
+the provided [`Crystal`](@ref). The bond distance cutoff `max_dist` is in global
+length units, e.g. angstrom.
+
+The same information can be obtained from [`print_site`](@ref),
+[`print_bond`](@ref), and [`reference_bonds`](@ref).
 """
 function print_symmetry_table(cryst::Crystal, max_dist; io=stdout)
     for b in reference_bonds(cryst, max_dist)
@@ -206,11 +208,11 @@ end
 
 
 """
-    print_suggested_frame(cryst, i)
+    print_suggested_frame(cryst::Crystal, i)
 
-Print a suggested reference frame, as a rotation matrix `R`, that can be used as
-input to `print_site()`. The purpose is to simplify the description of allowed
-anisotropies.
+Prints a suggested reference frame for atom `i`. This is given as a rotation `R`
+of the Cartesian basis. When passed to [`print_site`](@ref), it may simplify the
+description of symmetry-allowed anisotropies.
 """
 function print_suggested_frame(cryst::Crystal, i::Int)
     R = suggest_frame_for_atom(cryst, i)
@@ -220,13 +222,13 @@ end
 
 
 """
-    print_site(cryst, i; i_ref=i, R=I)
+    print_site(cryst::Crystal, i; i_ref=i, R=I)
 
-Print symmetry information for the site `i`, including allowed g-tensor and
-allowed anisotropy operator.  An optional symmetry-equivalent reference atom
-`i_ref` can be provided to keep a consistent meaning of the free parameters. An
-optional rotation matrix `R` can map to a new Cartesian reference frame for
-expression of the allowed anisotropy.
+Print symmetry information for atom `i`, including allowed g-tensor and allowed
+anisotropy operator.  An optional symmetry-equivalent reference atom `i_ref` can
+be provided to keep a consistent meaning of the free parameters. An optional
+rotation matrix `R` will transform the Cartesian basis for expression of the
+allowed anisotropy.
 """
 function print_site(cryst::Crystal, i; i_ref=i, R=Mat3(I), ks=[2,4,6], io=stdout)
     R_global = convert(Mat3, R)

src/Symmetry/SymmetryAnalysis.jl

@@ -195,10 +195,10 @@ end
 
 """    reference_bonds(cryst::Crystal, max_dist)
 
-Returns a full list of bonds, one for each symmetry equivalence class, up to
-distance `max_dist`. The reference bond `b` for each equivalence class is
-selected according to a scoring system that prioritizes simplification of the
-elements in `basis_for_symmetry_allowed_couplings(cryst, b)`."""
+Returns a list of [`Bond`](@ref)s, one for each symmetry equivalence class, up
+to the `max_dist` cutoff in length units. These reference bonds are
+heuristically selected to simplify the expression of symmetry-allowed
+interactions."""
 function reference_bonds(cryst::Crystal, max_dist::Float64; min_dist=0.0)
     # Bonds, one for each equivalence class
     ref_bonds = Bond[]