Change convention from forces to energy_grad

src/Integrators.jl

@@ -63,22 +63,22 @@ such as [`LocalSampler`](@ref).
 function step! end
 
 function step!(sys::System{0}, integrator::Langevin)
-    (B, s₁, f₁, r₁, ξ) = get_dipole_buffers(sys, 5)
+    (∇E, s₁, f₁, r₁, ξ) = get_dipole_buffers(sys, 5)
     (; kT, λ, Δt) = integrator
     s = sys.dipoles
 
     randn!(sys.rng, ξ)
     ξ .*= √(2λ*kT)
 
     # Euler step
-    set_forces_dipoles!(B, s, sys)
-    @. f₁ = rhs_dipole(s, B, λ)
+    set_energy_grad_dipoles!(∇E, s, sys)
+    @. f₁ = rhs_dipole(s, -∇E, λ)
     @. r₁ = rhs_dipole(s, ξ)   # note absence of λ argument -- noise only appears once in rhs.
     @. s₁ = s + Δt * f₁ + √Δt * r₁
 
     # Corrector step
-    set_forces_dipoles!(B, s₁, sys)
-    @. s = s + 0.5 * Δt * (f₁ + rhs_dipole(s₁, B, λ)) + 0.5 * √Δt * (r₁ + rhs_dipole(s₁, ξ))
+    set_energy_grad_dipoles!(∇E, s₁, sys)
+    @. s = s + 0.5 * Δt * (f₁ + rhs_dipole(s₁, -∇E, λ)) + 0.5 * √Δt * (r₁ + rhs_dipole(s₁, ξ))
     @. s = normalize_dipole(s, sys.κs)
     nothing
 end
@@ -93,7 +93,7 @@ function step!(sys::System{0}, integrator::ImplicitMidpoint)
     s = sys.dipoles
     (; Δt, atol) = integrator
 
-    (B, s̄, ŝ, s̄′) = get_dipole_buffers(sys, 4)
+    (∇E, s̄, ŝ, s̄′) = get_dipole_buffers(sys, 4)
 
     # Initial guess for midpoint
     @. s̄ = s
@@ -103,8 +103,8 @@ function step!(sys::System{0}, integrator::ImplicitMidpoint)
         # Integration step for current best guess of midpoint s̄. Produces
         # improved midpoint estimator s̄′.
         @. ŝ = normalize_dipole(s̄, sys.κs)
-        set_forces_dipoles!(B, ŝ, sys)
-        @. s̄′ = s + 0.5 * Δt * rhs_dipole(ŝ, B)
+        set_energy_grad_dipoles!(∇E, ŝ, sys)
+        @. s̄′ = s + 0.5 * Δt * rhs_dipole(ŝ, -∇E)
 
         # If converged, then we can return
         if fast_isapprox(s̄, s̄′,atol=atol* √length(s̄))
@@ -143,22 +143,22 @@ end
 
 function step!(sys::System{N}, integrator::Langevin) where N
     (Z′, ΔZ₁, ΔZ₂, ξ, HZ) = get_coherent_buffers(sys, 5)
-    B = get_dipole_buffers(sys, 1) |> only
+    ∇E = get_dipole_buffers(sys, 1) |> only
     Z = sys.coherents
 
     randn!(sys.rng, ξ)
 
     # Prediction
     @. sys.dipoles = expected_spin(Z) # temporarily desyncs dipoles and coherents
-    set_forces_dipoles!(B, sys.dipoles, sys)
-    set_forces_coherents!(HZ, B, Z, sys)
+    set_energy_grad_dipoles!(∇E, sys.dipoles, sys)
+    set_energy_grad_coherents!(HZ, ∇E, Z, sys)
     rhs_langevin!(ΔZ₁, Z, ξ, HZ, integrator, sys)
     @. Z′ = normalize_ket(Z + ΔZ₁, sys.κs)
 
     # Correction
     @. sys.dipoles = expected_spin(Z′) # temporarily desyncs dipoles and coherents
-    set_forces_dipoles!(B, sys.dipoles, sys)
-    set_forces_coherents!(HZ, B, Z′, sys)
+    set_energy_grad_dipoles!(∇E, sys.dipoles, sys)
+    set_energy_grad_coherents!(HZ, ∇E, Z′, sys)
     rhs_langevin!(ΔZ₂, Z′, ξ, HZ, integrator, sys)
     @. Z = normalize_ket(Z + (ΔZ₁+ΔZ₂)/2, sys.κs)
 
@@ -186,7 +186,7 @@ end
 function step!(sys::System{N}, integrator::ImplicitMidpoint; max_iters=100) where N
     (; atol) = integrator
     (ΔZ, Z̄, Z′, Z″, HZ) = get_coherent_buffers(sys, 5)
-    B = get_dipole_buffers(sys, 1) |> only
+    ∇E = get_dipole_buffers(sys, 1) |> only
     Z = sys.coherents
 
     @. Z′ = Z 
@@ -196,8 +196,8 @@ function step!(sys::System{N}, integrator::ImplicitMidpoint; max_iters=100) wher
         @. Z̄ = (Z + Z′)/2
 
         @. sys.dipoles = expected_spin(Z̄) # temporarily desyncs dipoles and coherents
-        set_forces_dipoles!(B, sys.dipoles, sys)
-        set_forces_coherents!(HZ, B, Z̄, sys)
+        set_energy_grad_dipoles!(∇E, sys.dipoles, sys)
+        set_energy_grad_coherents!(HZ, ∇E, Z̄, sys)
         rhs_ll!(ΔZ, HZ, integrator, sys)
 
         @. Z″ = Z + ΔZ

src/Optimization.jl

@@ -2,13 +2,13 @@
 # projective parameterization is formally the same for both dipoles and coherent
 # states (the element types are just real in the first case, compelx in the
 # second).
-function set_forces_optim!(Hgrad, B, sys::System{N}) where {N}
-    Sunny.set_forces_dipoles!(B, sys.dipoles, sys)
-    Sunny.set_forces_coherents!(Hgrad, B, sys.coherents, sys)
+function set_forces_optim!(∇H, ∇H_dip, sys::System{N}) where {N}
+    Sunny.set_energy_grad_dipoles!(∇H_dip, sys.dipoles, sys)
+    Sunny.set_energy_grad_coherents!(∇H, ∇H_dip, sys.coherents, sys)
 end
 
-function set_forces_optim!(Hgrad, _, sys::System{0}) 
-    Sunny.set_forces_dipoles!(Hgrad, sys.dipoles, sys)
+function set_forces_optim!(∇H, _, sys::System{0}) 
+    Sunny.set_energy_grad_dipoles!(∇H, sys.dipoles, sys)
 end
 
 @inline function set_spin_optim!(sys::System{N}, α, z, site) where N
@@ -64,7 +64,6 @@ end
 # Calculate dH(u(α))/dα 
 function optim_gradient!(buf, αs, zs, B, sys::System{N}, halted, quickmode=false) where N
     T = N == 0 ? Vec3 : CVec{N} 
-    jacsign = N == 0 ? -1 : 1  # Because set_forces_dipoles gives -dH/dS whereas set_forces_coherents gives dH/dZ*
     αs = reinterpret(reshape, T, αs)
     Hgrad = reinterpret(reshape, T, buf)
 
@@ -86,7 +85,7 @@ function optim_gradient!(buf, αs, zs, B, sys::System{N}, halted, quickmode=fals
     set_forces_optim!(Hgrad, B, sys)
 
     for site in all_sites(sys)
-        Hgrad[site] = jacsign * apply_projective_jacobian(Hgrad[site], αs[site], zs[site])
+        Hgrad[site] = apply_projective_jacobian(Hgrad[site], αs[site], zs[site])
     end
 end
 

src/System/Ewald.jl

@@ -154,7 +154,7 @@ end
 
 # Use FFT to accumulate the entire field -dE/ds for long-range dipole-dipole
 # interactions
-function accum_ewald_force!(B::Array{Vec3, 4}, dipoles::Array{Vec3, 4}, sys::System{N}) where N
+function accum_ewald_grad!(∇E::Array{Vec3, 4}, dipoles::Array{Vec3, 4}, sys::System{N}) where N
     (; ewald, units, gs) = sys
     (; μ, FA, Fμ, Fϕ, ϕ, plan, ift_plan) = ewald
 
@@ -178,7 +178,7 @@ function accum_ewald_force!(B::Array{Vec3, 4}, dipoles::Array{Vec3, 4}, sys::Sys
     mul!(ϕr, ift_plan, Fϕ)
 
     for site in all_sites(sys)
-        B[site] -= 2 * units.μB * (gs[site]' * ϕ[site])
+        ∇E[site] += 2 * units.μB * (gs[site]' * ϕ[site])
     end
 end
 

src/System/Interactions.jl

@@ -259,47 +259,47 @@ function energy_aux(sys::System{N}, ints::Interactions, i::Int, cells, foreachbo
     return E
 end
 
-# Updates B in-place to hold negative energy gradient, -dE/ds, for each spin.
-function set_forces_dipoles!(B, dipoles::Array{Vec3, 4}, sys::System{N}) where N
+# Updates ∇E in-place to hold energy gradient, dE/ds, for each spin.
+function set_energy_grad_dipoles!(∇E, dipoles::Array{Vec3, 4}, sys::System{N}) where N
     (; crystal, latsize, extfield, ewald) = sys
 
-    fill!(B, zero(Vec3))
+    fill!(∇E, zero(Vec3))
 
     # Zeeman coupling
     for site in all_sites(sys)
-        B[site] += sys.units.μB * (sys.gs[site]' * extfield[site])
+        ∇E[site] -= sys.units.μB * (sys.gs[site]' * extfield[site])
     end
 
     # Anisotropies and exchange interactions
     for i in 1:natoms(crystal)
         if is_homogeneous(sys)
             # Interaction is the same at every cell
             interactions = sys.interactions_union[i]
-            set_forces_aux!(B, dipoles, interactions, sys, i, all_cells(sys), homog_bond_iterator(latsize))
+            set_energy_grad_aux!(∇E, dipoles, interactions, sys, i, all_cells(sys), homog_bond_iterator(latsize))
         else
             for cell in all_cells(sys)
                 # There is a different interaction at every cell
                 interactions = sys.interactions_union[cell,i]
-                set_forces_aux!(B, dipoles, interactions, sys, i, (cell,), inhomog_bond_iterator(latsize, cell))
+                set_energy_grad_aux!(∇E, dipoles, interactions, sys, i, (cell,), inhomog_bond_iterator(latsize, cell))
             end
         end
     end
 
     if !isnothing(ewald)
-        accum_ewald_force!(B, dipoles, sys)
+        accum_ewald_grad!(∇E, dipoles, sys)
     end
 end
 
-# Calculate the force `B' for the sublattice `i' at all elements of `cells`. The
-# function `foreachbond` enables efficient iteration over neighboring cell
-# pairs.
-function set_forces_aux!(B, dipoles::Array{Vec3, 4}, ints::Interactions, sys::System{N}, i::Int, cells, foreachbond) where N
+# Calculate the energy gradient `∇H' for the sublattice `i' at all elements of
+# `cells`. The function `foreachbond` enables efficient iteration over
+# neighboring cell pairs.
+function set_energy_grad_aux!(∇E, dipoles::Array{Vec3, 4}, ints::Interactions, sys::System{N}, i::Int, cells, foreachbond) where N
     # Single-ion anisotropy only contributes in dipole mode. In SU(N) mode, the
     # anisotropy matrix will be incorporated directly into ℌ.
     if N == 0
         for cell in cells
             s = dipoles[cell, i]
-            B[cell, i] -= energy_and_gradient_for_classical_anisotropy(s, ints.onsite.stvexp)[2]
+            ∇E[cell, i] += energy_and_gradient_for_classical_anisotropy(s, ints.onsite.stvexp)[2]
         end
     end
 
@@ -309,8 +309,8 @@ function set_forces_aux!(B, dipoles::Array{Vec3, 4}, ints::Interactions, sys::Sy
 
         # Bilinear
         J = coupling.bilin
-        B[site1] -= J  * sⱼ
-        B[site2] -= J' * sᵢ
+        ∇E[site1] += J  * sⱼ
+        ∇E[site2] += J' * sᵢ
 
         # Biquadratic
         if !iszero(coupling.biquad)
@@ -321,46 +321,46 @@ function set_forces_aux!(B, dipoles::Array{Vec3, 4}, ints::Interactions, sys::Sy
                 Sⱼ = (sys.Ns[site2]-1)/2
                 S = √(Sᵢ*Sⱼ)
                 r = (1 - 1/S + 1/4S^2)
-                B[site1] -= J * (2r*sⱼ*(sᵢ⋅sⱼ) - sⱼ/2)
-                B[site2] -= J * (2r*sᵢ*(sᵢ⋅sⱼ) - sᵢ/2)
+                ∇E[site1] += J * (2r*sⱼ*(sᵢ⋅sⱼ) - sⱼ/2)
+                ∇E[site2] += J * (2r*sᵢ*(sᵢ⋅sⱼ) - sᵢ/2)
             elseif sys.mode == :large_S
-                B[site1] -= J * 2sⱼ*(sᵢ⋅sⱼ)
-                B[site2] -= J * 2sᵢ*(sᵢ⋅sⱼ)
+                ∇E[site1] += J * 2sⱼ*(sᵢ⋅sⱼ)
+                ∇E[site2] += J * 2sᵢ*(sᵢ⋅sⱼ)
             elseif sys.mode == :SUN
                 error("Biquadratic currently unsupported in SU(N) mode.")
             end
         end
     end
 end
 
-# Updates `HZ` in-place to hold `dH/dZ^*`, the analog in the Schroedinger formulation
-# of the classical quantity `dE/ds` (note sign).
-function set_forces_coherents!(HZ, B, Z, sys::System{N}) where N
+# Updates `HZ` in-place to hold `dH/dZ^*`, the analog in the Schroedinger
+# formulation of the quantity `dE/ds`.
+function set_energy_grad_coherents!(HZ, ∇E, Z, sys::System{N}) where N
     if is_homogeneous(sys)
         ints = interactions_homog(sys)
         for site in all_sites(sys)
             Λ = ints[to_atom(site)].onsite.matrep
-            HZ[site] = mul_spin_matrices(Λ, -B[site], Z[site])  
+            HZ[site] = mul_spin_matrices(Λ, ∇E[site], Z[site])  
         end 
     else
         ints = interactions_inhomog(sys)
         for site in all_sites(sys)
             Λ = ints[site].onsite.matrep
-            HZ[site] = mul_spin_matrices(Λ, -B[site], Z[site])  
+            HZ[site] = mul_spin_matrices(Λ, ∇E[site], Z[site])  
         end 
     end
 end
 
-# Returns (Λ + B⋅S) Z
-@generated function mul_spin_matrices(Λ, B::Sunny.Vec3, Z::Sunny.CVec{N}) where N
+# Returns (Λ - ∇E⋅S) Z
+@generated function mul_spin_matrices(Λ, ∇E::Sunny.Vec3, Z::Sunny.CVec{N}) where N
     S = spin_matrices(; N)
     out = map(1:N) do i
         out_i = map(1:N) do j
             terms = Any[:(Λ[$i,$j])]
             for α = 1:3
                 S_αij = S[α][i,j]
                 if !iszero(S_αij)
-                    push!(terms, :(B[$α] * $S_αij))
+                    push!(terms, :(∇E[$α] * $S_αij))
                 end
             end
             :(+($(terms...)) * Z[$j])
@@ -407,7 +407,7 @@ end
 Returns the effective local field (force) at each site, ``𝐁 = -∂E/∂𝐬``.
 """
 function forces(sys::System{N}) where N
-    B = zero(sys.dipoles)
-    set_forces_dipoles!(B, sys.dipoles, sys)
-    return B
+    ∇E = zero(sys.dipoles)
+    set_energy_grad_dipoles!(∇E, sys.dipoles, sys)
+    return -∇E
 end