Expose clone_system in public API (#259)

src/Sunny.jl

@@ -52,7 +52,7 @@ include("System/PairExchange.jl")
 include("System/OnsiteCoupling.jl")
 include("System/Ewald.jl")
 include("System/Interactions.jl")
-export SpinInfo, System, Site, eachsite, position_to_site, global_position, magnetic_moment, 
+export SpinInfo, System, Site, clone_system, eachsite, position_to_site, global_position, magnetic_moment,
     set_coherent!, set_dipole!, polarize_spins!, randomize_spins!, set_spin_rescaling!, energy, energy_per_site,
     spin_label, set_onsite_coupling!, set_pair_coupling!, set_exchange!, dmvec, enable_dipole_dipole!,
     set_external_field!, to_inhomogeneous, set_external_field_at!, set_vacancy_at!, set_onsite_coupling_at!,

src/System/Ewald.jl

@@ -21,14 +21,16 @@ function Ewald(sys::System{N}) where N
     return Ewald(A, μ, ϕ, FA, Fμ, Fϕ, plan, ift_plan)
 end
 
-# Clone all mutable state within Ewald. Note that `A`, `FA` are immutable data.
-# We're also treating FFTW plans as immutable. This is not 100% correct, because
-# plans mutably cache inverse plans, which may possibly lead to data races in a
-# multithreaded context. Unfortunately, FFTW does not yet provide a copy
-# function. For context, see https://github.com/JuliaMath/FFTW.jl/issues/261.
+# Ideally, this would clone all mutable state within Ewald. Note that `A`, `FA`
+# are immutable data. A blocker is that FFTW plans cannot currently be copied,
+# and it is not 100% clear whether they can be treated as immutable. For
+# example, they cache inverse plans, which may possibly lead to data races in a
+# multithreaded context. See discussion at
+# https://github.com/JuliaMath/FFTW.jl/issues/261.
 function clone_ewald(ewald::Ewald)
+    error("Not supported")
     (; A, μ, ϕ, FA, Fμ, Fϕ, plan, ift_plan) = ewald
-    return Ewald(A, copy(μ), copy(ϕ), FA, copy(Fμ), copy(Fϕ), plan, ift_plan)
+    return Ewald(A, copy(μ), copy(ϕ), FA, copy(Fμ), copy(Fϕ), copy(plan), copy(ift_plan))
 end
 
 # Tensor product of 3-vectors

src/System/System.jl

@@ -125,14 +125,17 @@ end
 # stack of `clone` functions instead.
 Base.deepcopy(_::System) = error("Use `clone_system` instead of `deepcopy`.")
 
-# Creates a clone of the system where all the mutable internal data is copied.
-# It is intended to be thread-safe to use the original and the copied systems,
-# without any restrictions, but see caveats in `clone_ewald()`.
+"""
+    clone_system(sys::System)
+
+Creates a full clone of the system, such that mutable updates to one copy will
+not affect the other, and thread safety is guaranteed.
+"""
 function clone_system(sys::System{N}) where N
     (; origin, mode, crystal, latsize, Ns, gs, κs, extfield, interactions_union, ewald, dipoles, coherents, units, rng) = sys
 
     origin_clone = isnothing(origin) ? nothing : clone_system(origin)
-    ewald_clone  = isnothing(ewald)  ? nothing : clone_ewald(ewald)
+    ewald_clone = nothing # TODO: use clone_ewald(ewald)
 
     # Dynamically dispatch to the correct `map` function for either homogeneous
     # (Vector) or inhomogeneous interactions (4D Array)
@@ -142,9 +145,18 @@ function clone_system(sys::System{N}) where N
     empty_dipole_buffers = Array{Vec3, 4}[]
     empty_coherent_buffers = Array{CVec{N}, 4}[]
 
-    System(origin_clone, mode, crystal, latsize, Ns, copy(κs), copy(gs),
-           interactions_clone, ewald_clone, copy(extfield), copy(dipoles), copy(coherents),
-           empty_dipole_buffers, empty_coherent_buffers, units, copy(rng))
+    ret = System(origin_clone, mode, crystal, latsize, Ns, copy(κs), copy(gs),
+                 interactions_clone, ewald_clone, copy(extfield), copy(dipoles), copy(coherents),
+                 empty_dipole_buffers, empty_coherent_buffers, units, copy(rng))
+
+    if !isnothing(ewald)
+        # At the moment, clone_ewald is unavailable, so instead rebuild the
+        # Ewald data structures from scratch. This might be fixed eventually.
+        # See https://github.com/JuliaMath/FFTW.jl/issues/261.
+        enable_dipole_dipole!(ret)
+    end
+
+    return ret
 end