Default seed for a System is generated from the task-local RNG

docs/make.jl

@@ -37,6 +37,7 @@ function build_examples(example_sources, destdir)
             """
             # Download this example as [Julia file]($assetsdir/scripts/$name.jl) or [Jupyter notebook]($assetsdir/notebooks/$name.ipynb).
 
+            import Random; Random.seed!(0); #hide
             """ * str
         end
         # Write to `src/$destpath/$name.md`

docs/src/parallelism.md

@@ -15,19 +15,20 @@ copied and pasted into your preferred Julia development environment.
 
 The serial approach to calculating a structure factor, covered in the [FeI₂
 tutorial](@ref "4. Generalized spin dynamics of FeI₂ at finite *T*"), involves
-thermalizing a spin `System` and then calling [`add_sample!`](@ref).
+thermalizing a spin [`System`](@ref) and then calling [`add_sample!`](@ref).
 `add_sample!` uses the state of the `System` as an initial condition for the
 calculation of a dynamical trajectory. The correlations of the trajectory are
 calculated and accumulated into a running average of the ``\mathcal{S}(𝐪,ω)``.
 This sequence is repeated to generate additional samples.
 
 To illustrate, we'll set up a a simple model: a spin-1 antiferromagnet on a BCC
-crystal. 
+crystal. Constructing the `System` with a specific random number `seed` ensures
+full reproducibility of the simulation.
 
 ```julia
 using Sunny, GLMakie
 
-function make_system(; seed=nothing)
+function make_system(seed)
     latvecs = lattice_vectors(1, 1, 1, 90, 90, 90)
     positions = [[0, 0, 0]/2, [1, 1, 1]/2]
     cryst = Crystal(latvecs, positions)
@@ -36,7 +37,7 @@ function make_system(; seed=nothing)
     return sys
 end
 
-sys = make_system()
+sys = make_system(0)
 ```
 
 A serial calculation of [`SampledCorrelations`](@ref) involving the
@@ -80,11 +81,12 @@ Before going further, make sure that `Threads.nthreads()` returns a number great
 
 We will use multithreading in a very simple way, essentially employing a
 distributed memory approach to avoid conflicts around memory access. First
-preallocate a number of systems and correlations.
+preallocate a number of systems and correlations. The integer `id` of each
+system is used as its random number seed.
 
 ```julia
 npar = Threads.nthreads()
-systems = [make_system(; seed=id) for id in 1:npar]
+systems = [make_system(id) for id in 1:npar]
 scs = [SampledCorrelations(sys; dt=0.1, energies, measure) for _ in 1:npar]
 ```
 
@@ -151,7 +153,7 @@ environments. This is easily achieved with the `@everywhere` macro.
 ```julia
 @everywhere using Sunny
 
-@everywhere function make_system(; seed=nothing)
+@everywhere function make_system(seed)
     latvecs = lattice_vectors(1, 1, 1, 90, 90, 90)
     positions = [[0, 0, 0]/2, [1, 1, 1]/2]
     cryst = Crystal(latvecs, positions)
@@ -181,7 +183,7 @@ called `scs`.
 
 ```julia
 scs = pmap(1:ncores) do id
-    sys = make_system(; seed=id)
+    sys = make_system(id)
     sc = SampledCorrelations(sys; dt=0.1, energies, measure)
     integrator = Langevin(0.05; damping=0.2, kT=0.5)
 

docs/src/versions.md

@@ -6,6 +6,9 @@
 * Better error message when a $g$-tensor is symmetry disallowed.
 * Higher-precision convergence in [`minimize_energy!`](@ref).
 * Fix [`minimize_energy!`](@ref) when used with [`set_vacancy_at!`](@ref).
+* The `System` constructor now, by default, seeds its internal random number
+  generator with `seed=rand(UInt)`. Note that Julia's global random number
+  generator can itself be seeded with `Random.seed!`.
 
 ## v0.7.3
 (Nov 12, 2024)

examples/02_LLD_CoRh2O4.jl

@@ -9,9 +9,8 @@
 # classical-to-quantum correction factor, the resulting intensities can be
 # compared to inelastic neutron scattering data.
 
-# Construct the system as in the [previous tutorial](@ref "1. Spin wave
-# simulations of CoRh₂O₄"). For this antiferromagnetic model on the diamond
-# cubic lattice, the ground state is unfrustrated Néel order.
+# Construct the CoRh₂O₄ antiferromagnet as before. Energy minimization yields
+# the expected Néel order.
 
 using Sunny, GLMakie
 
@@ -22,7 +21,7 @@ cryst = Crystal(latvecs, [[1/8, 1/8, 1/8]], 227)
 
 sys = System(cryst, [1 => Moment(s=3/2, g=2)], :dipole)
 J = 0.63 # (meV)
-set_exchange!(sys, J, Bond(2, 3, [0,0,0]))
+set_exchange!(sys, J, Bond(2, 3, [0, 0, 0]))
 randomize_spins!(sys)
 minimize_energy!(sys)
 plot_spins(sys; color=[S[3] for S in sys.dipoles])

examples/03_LSWT_SU3_FeI2.jl

@@ -84,11 +84,11 @@ print_symmetry_table(cryst, 8.0)
 # magnitude. This physics is, however, well captured with a theory of SU(_N_)
 # coherent states, where ``N = 2S+1 = 3`` is the number of levels. Activate this
 # generalized theory by selecting `:SUN` mode instead of `:dipole` mode.
-# 
-# An optional random number `seed` will make the calculations exactly
-# reproducible.
+#
+# An optional `seed` for random number generation can be used to to make the
+# calculation exactly reproducible.
 
-sys = System(cryst, [1 => Moment(s=1, g=2)], :SUN, seed=2)
+sys = System(cryst, [1 => Moment(s=1, g=2)], :SUN; seed=2)
 
 # Set the exchange interactions for FeI₂ following the fits of [Bai et
 # al.](https://doi.org/10.1038/s41567-020-01110-1)
@@ -156,42 +156,37 @@ sys = resize_supercell(sys, (4, 4, 4))
 randomize_spins!(sys)
 minimize_energy!(sys)
 
-# A positive number above indicates that the procedure has converged to a local
-# energy minimum. The configuration, however, may still have defects. This can
-# be checked by visualizing the expected spin dipoles, colored according to
-# their ``z``-components.
+# The positive step-count above indicates successful convergence to a local
+# energy minimum. Defects, however, are visually apparent.
 
 plot_spins(sys; color=[S[3] for S in sys.dipoles])
 
-# To better understand the spin configuration, we could inspect the static
-# structure factor ``\mathcal{S}(𝐪)`` in the 3D space of momenta ``𝐪``. The
-# general tool for this analysis is [`SampledCorrelationsStatic`](@ref). For the
-# present purposes, however, it is more convenient to use
-# [`print_wrapped_intensities`](@ref), which reports ``\mathcal{S}(𝐪)`` with
-# periodic wrapping of all commensurate ``𝐪`` wavevectors into the first
-# Brillouin zone.
+# 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.
 
 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]``. Small systems can spontaneously break this symmetry, but
-# for larger systems, defects and competing domains are to be expected.
-# Nonetheless, `print_wrapped_intensities` shows large intensities consistent
-# with a subset of the known ordering wavevectors.
+# ``[-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.
 #
 # Let's break the three-fold symmetry by hand. The function
-# [`suggest_magnetic_supercell`](@ref) takes one or more ``𝐤`` modes, and
-# suggests a magnetic cell shape that is commensurate.
+# [`suggest_magnetic_supercell`](@ref) takes any number of ``𝐤`` modes and
+# suggests a magnetic cell shape that is commensurate with all of them.
 
 suggest_magnetic_supercell([[0, -1/4, 1/4]])
 
-# Calling [`reshape_supercell`](@ref) yields a much smaller system, making it
-# much easier to find the global energy minimum. Plot the system again, now
-# including "ghost" spins out to 12Å, to verify that the magnetic order is
-# consistent with FeI₂.
+# 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₂.
 
 sys_min = reshape_supercell(sys, [1 0 0; 0 2 1; 0 -2 1])
 randomize_spins!(sys_min)

examples/05_MC_Ising.jl

@@ -16,7 +16,7 @@ crystal = Crystal(latvecs, [[0, 0, 0]])
 # dipole ``𝐒`` is ``-𝐁⋅𝐒``. The system size is 128×128.
 
 L = 128
-sys = System(crystal, [1 => Moment(s=1, g=-1)], :dipole; dims=(L, L, 1), seed=0)
+sys = System(crystal, [1 => Moment(s=1, g=-1)], :dipole; dims=(L, L, 1))
 polarize_spins!(sys, (0, 0, 1))
 
 # Use [`set_exchange!`](@ref) to include a ferromagnetic Heisenberg interaction

examples/07_Dipole_Dipole.jl

@@ -18,8 +18,10 @@ positions = [[0, 0, 0]]
 cryst = Crystal(latvecs, positions, 227)
 view_crystal(cryst)
 
-# Create a system and reshape to the primitive cell, which contains four atoms.
-# Add antiferromagnetic nearest neighbor exchange interactions.
+# Create a [`System`](@ref) with a random number `seed` that was empirically
+# selected to produce the desired type of spontaneous symmetry breaking. Reshape
+# to the primitive cell, which contains four atoms. Add antiferromagnetic
+# nearest neighbor exchange interactions.
 
 sys = System(cryst, [1 => Moment(s=7/2, g=2)], :dipole; seed=0)
 sys = reshape_supercell(sys, primitive_cell(cryst))

examples/spinw_tutorials/SW10_Energy_cut.jl

@@ -33,7 +33,7 @@ grid = q_space_grid(cryst, [1, 0, 0], range(0, 2, 201), [0, 1, 0], range(0, 2, 2
 # Apply a line broadening with a full-width half-max of 0.2 meV to approximately
 # capture intensities between 3.5 and 4.0 meV.
 
-swt = SpinWaveTheory(sys; measure=ssf_perp(sys))
+swt = SpinWaveTheory(sys; measure=ssf_trace(sys))
 res = intensities(swt, grid; energies=[3.75], kernel=gaussian(fwhm=0.2))
 plot_intensities(res; units)
 

examples/spinw_tutorials/SW15_Ba3NbFe3Si2O14.jl

@@ -30,7 +30,7 @@ view_crystal(cryst)
 # [Loire et al., Phys. Rev. Lett. **106**, 207201
 # (2011)](http://dx.doi.org/10.1103/PhysRevLett.106.207201).
 
-sys = System(cryst, [1 => Moment(s=5/2, g=2)], :dipole; seed=0)
+sys = System(cryst, [1 => Moment(s=5/2, g=2)], :dipole)
 J₁ = 0.85
 J₂ = 0.24
 J₃ = 0.053

examples/spinw_tutorials/SW18_Distorted_kagome.jl

@@ -25,7 +25,7 @@ view_crystal(cryst)
 # Define the interactions.
 
 moments = [1 => Moment(s=1/2, g=2), 3 => Moment(s=1/2, g=2)]
-sys = System(cryst, moments, :dipole, seed=0)
+sys = System(cryst, moments, :dipole)
 J   = -2
 Jp  = -1
 Jab = 0.75

examples/spinw_tutorials/SW19_Different_Ions.jl

@@ -25,7 +25,7 @@ J_Cu_Cu = 1.0
 J_Fe_Fe = 1.0
 J_Cu_Fe = -0.1
 moments = [1 => Moment(s=1/2, g=2), 2 => Moment(s=2, g=2)]
-sys = System(cryst, moments, :dipole; dims=(2, 1, 1), seed=0)
+sys = System(cryst, moments, :dipole; dims=(2, 1, 1))
 set_exchange!(sys, J_Cu_Cu, Bond(1, 1, [-1, 0, 0]))
 set_exchange!(sys, J_Fe_Fe, Bond(2, 2, [-1, 0, 0]))
 set_exchange!(sys, J_Cu_Fe, Bond(2, 1, [0, 1, 0]))
@@ -50,17 +50,17 @@ path = q_space_path(cryst, qs, 400)
 fig = Figure(size=(768,600))
 
 formfactors = [1 => FormFactor("Cu2"), 2 => FormFactor("Fe2")]
-swt = SpinWaveTheory(sys; measure=ssf_perp(sys; formfactors))
+swt = SpinWaveTheory(sys; measure=ssf_trace(sys; formfactors))
 res = intensities_bands(swt, path)
 plot_intensities!(fig[1, 1], res; units, title="All correlations")
 
 formfactors = [1 => FormFactor("Cu2"), 2 => zero(FormFactor)]
-swt = SpinWaveTheory(sys; measure=ssf_perp(sys; formfactors))
+swt = SpinWaveTheory(sys; measure=ssf_trace(sys; formfactors))
 res = intensities_bands(swt, path)
 plot_intensities!(fig[1, 2], res; units, title="Cu-Cu correlations")
 
 formfactors = [1 => zero(FormFactor), 2 => FormFactor("Fe2")]
-swt = SpinWaveTheory(sys; measure=ssf_perp(sys; formfactors))
+swt = SpinWaveTheory(sys; measure=ssf_trace(sys; formfactors))
 res = intensities_bands(swt, path)
 plot_intensities!(fig[2, 2], res; units, title="Fe-Fe correlations")
 

src/System/System.jl

@@ -23,13 +23,15 @@ be renormalized to improve accuracy. To disable this renormalization, use the
 mode `:dipole_uncorrected`, which corresponds to the formal ``s → ∞`` limit. For
 details, see the documentation page: [Interaction Renormalization](@ref).
 
-An integer `seed` for the random number generator can optionally be specified to
-enable reproducible calculations.
+Stochastic operations on this system can be made reproducible by selecting a
+`seed` for this system's pseudo-random number generator. The default system seed
+will be generated from Julia's task-local RNG, which can itself be seeded using
+`Random.seed!`.
 
 All spins are initially polarized in the global ``z``-direction.
 """
 function System(crystal::Crystal, moments::Vector{Pair{Int, Moment}}, mode::Symbol;
-                dims::NTuple{3,Int}=(1, 1, 1), seed::Union{Int,Nothing}=nothing, units=nothing)
+                dims::NTuple{3,Int}=(1, 1, 1), seed=nothing, units=nothing)
     if !isnothing(units)
         @warn "units argument to System is deprecated and will be ignored!"
     end
@@ -47,7 +49,7 @@ function System(crystal::Crystal, moments::Vector{Pair{Int, Moment}}, mode::Symb
     if !isnothing(crystal.root)
         @assert crystal.latvecs == crystal.root.latvecs
     end
-    
+
     na = natoms(crystal)
 
     moments = propagate_moments(crystal, moments)
@@ -78,7 +80,8 @@ function System(crystal::Crystal, moments::Vector{Pair{Int, Moment}}, mode::Symb
     coherents = fill(zero(CVec{N}), 1, 1, 1, na)
     dipole_buffers = Array{Vec3, 4}[]
     coherent_buffers = Array{CVec{N}, 4}[]
-    rng = isnothing(seed) ? Random.Xoshiro() : Random.Xoshiro(seed)
+
+    rng = isnothing(seed) ? Random.Xoshiro(rand(UInt64, 4)...) : Random.Xoshiro(seed)
 
     ret = System(nothing, mode, crystal, (1, 1, 1), Ns, κs, gs, interactions, ewald,
                  extfield, dipoles, coherents, dipole_buffers, coherent_buffers, rng)