Clarify internal nature of energy_grad_*

src/System/Interactions.jl

@@ -379,6 +379,19 @@ function set_energy_grad_coherents!(HZ, Z, sys::System{N}) where N
     @. dE_ds = dipoles = Vec3(0,0,0)
 end
 
+# Internal testing functions
+function energy_grad_dipoles(sys::System{N}) where N
+    ∇E = zero(sys.dipoles)
+    set_energy_grad_dipoles!(∇E, sys.dipoles, sys)
+    return ∇E
+end
+function energy_grad_coherents(sys::System{N}) where N
+    ∇E = zero(sys.coherents)
+    set_energy_grad_coherents!(∇E, sys.coherents, sys)
+    return ∇E
+end
+
+
 # Returns (Λ + (dE/ds)⋅S) Z
 @generated function mul_spin_matrices(Λ, dE_ds::Sunny.Vec3, Z::Sunny.CVec{N}) where N
     S = spin_matrices(; N)
@@ -428,9 +441,3 @@ function homog_bond_iterator(latsize)
         end
     end
 end
-
-function energy_grad(sys::System{N}) where N
-    ∇E = zero(sys.dipoles)
-    set_energy_grad_dipoles!(∇E, sys.dipoles, sys)
-    return ∇E
-end

test/test_ewald.jl

@@ -48,7 +48,7 @@
     # Calculate energy gradient using a sum over pairs, or using an FFT-based
     # convolution
     ∇E = [Sunny.ewald_grad_at(sys, site) for site in eachsite(sys)]
-    @test isapprox(Sunny.energy_grad(sys), ∇E; atol=1e-12)
+    @test isapprox(Sunny.energy_grad_dipoles(sys), ∇E; atol=1e-12)
 
     # Calculation of energy as a sum over pairs
     E = sum((1/2)d⋅b for (d, b) in zip(sys.dipoles, ∇E))

test/test_units_scaling.jl

@@ -22,7 +22,7 @@
             enable_dipole_dipole!(sys)
         end
         randomize_spins!(sys)
-        return energy(sys), Sunny.energy_grad(sys)
+        return energy(sys), Sunny.energy_grad_dipoles(sys)
     end
 
     # All exchange interactions should be invariant under changes to physical