From 6fbaf150153edd1ac5826bfa5a5df8bbabb4828f Mon Sep 17 00:00:00 2001 From: ddahlbom Date: Fri, 21 Jul 2023 16:09:41 -0400 Subject: [PATCH] Change convention from forces to energy_grad --- src/Integrators.jl | 32 ++++++++++----------- src/Optimization.jl | 13 ++++----- src/System/Ewald.jl | 4 +-- src/System/Interactions.jl | 58 +++++++++++++++++++------------------- 4 files changed, 53 insertions(+), 54 deletions(-) diff --git a/src/Integrators.jl b/src/Integrators.jl index 4d51412a6..4694e40c0 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -63,7 +63,7 @@ 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 @@ -71,14 +71,14 @@ function step!(sys::System{0}, integrator::Langevin) ξ .*= √(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 diff --git a/src/Optimization.jl b/src/Optimization.jl index 14ee8b9da..e4a04fd28 100644 --- a/src/Optimization.jl +++ b/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 diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl index 9cb66a2db..ca696bd5c 100644 --- a/src/System/Ewald.jl +++ b/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 diff --git a/src/System/Interactions.jl b/src/System/Interactions.jl index ec82ca5e9..25a83d380 100644 --- a/src/System/Interactions.jl +++ b/src/System/Interactions.jl @@ -259,15 +259,15 @@ 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 @@ -275,31 +275,31 @@ function set_forces_dipoles!(B, dipoles::Array{Vec3, 4}, sys::System{N}) where N 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,11 +321,11 @@ 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 @@ -333,26 +333,26 @@ function set_forces_aux!(B, dipoles::Array{Vec3, 4}, ints::Interactions, sys::Sy 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 @@ -360,7 +360,7 @@ end 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