diff --git a/src/Integrators.jl b/src/Integrators.jl
index 34fb275c7..28bc850e2 100644
--- a/src/Integrators.jl
+++ b/src/Integrators.jl
@@ -259,10 +259,12 @@ end
     (; dt, damping) = integrator
     λ = damping
 
-    if iszero(λ)
-        @. Δs = -s × (- dt*∇E)
-    else
-        @. Δs = -s × (- dt*∇E + ξ - dt*λ*(s × ∇E))
+    for i = eachindex(Δs)
+        if iszero(λ)
+            Δs[i] = -s[i] × (- dt*∇E[i])
+        else
+            Δs[i] = -s[i] × (- dt*∇E[i] + ξ[i] - dt*λ*(s[i] × ∇E[i]))
+        end
     end
 end
 
@@ -383,18 +385,24 @@ function step!(sys::System{0}, integrator::ImplicitMidpoint; max_iters=100)
 
     for _ in 1:max_iters
         # Current guess for midpoint ŝ
-        @. ŝ = normalize_dipole((s + s′)/2, sys.κs)
+        for i = eachindex(ŝ )
+            ŝ[i] = normalize_dipole((s[i] + s′[i])/2, sys.κs[i])
+        end
 
         set_energy_grad_dipoles!(∇E, ŝ, sys)
         rhs_dipole!(Δs, ŝ, ξ, ∇E, integrator)
 
-        @. s″ = s + Δs
+        for i = eachindex(s″)
+            s″[i] = s[i] + Δs[i]
+        end
 
         # If converged, then we can return
         if fast_isapprox(s′, s″; atol)
             # Normalization here should not be necessary in principle, but it
             # could be useful in practice for finite `atol`.
-            @. s = normalize_dipole(s″, sys.κs)
+            for i = eachindex(s)
+                s[i] = normalize_dipole(s″[i], sys.κs[i])
+            end
             return
         end