Flexible Units (#285)

* More general Units constants
* Replace `set_external_field!` for `set_field``, which takes an energy
* Vacuum permittivity required for `enable_dipole_dipole!`

docs/src/versions.md

@@ -7,6 +7,10 @@
   the Bohr magneton, ``μ_B``. For model systems where the Zeeman coupling aligns
   spin dipole with field (e.g., the Ising model convention), create a `SpinInfo`
   with `g=-1`. ([PR 284](https://github.com/SunnySuite/Sunny.jl/pull/284)).
+* More flexible [`Units`](@ref) system. `set_external_field!` is deprecated in
+  favor of [`set_field!`](@ref), which now expects a field in energy units.
+  [`enable_dipole_dipole!`](@ref) now expects a scale parameter ``μ_0 μ_B^2``
+  that can be obtained from `units.vacuum_permeability`.
 
 ## v0.6.0
 (June 18, 2024)
@@ -335,7 +339,7 @@ The parameter `mode` is one of `:SUN` or `:dipole`.
 ### Setting interactions
 
 Interactions are now added mutably to an existing `System` using the following
-functions: [`set_external_field!`](@ref), [`set_exchange!`](@ref),
+functions: `set_external_field!`, [`set_exchange!`](@ref),
 [`set_onsite_coupling!`](@ref), [`enable_dipole_dipole!`](@ref).
 
 As a convenience, one can use [`dmvec(D)`](@ref) to convert a DM vector to a
@@ -355,8 +359,7 @@ combination of Stevens operators. To see this expansion use
 
 ### Inhomogeneous field
 
-An external field can be applied to a single site with
-[`set_external_field_at!`](@ref). 
+An external field can be applied to a single site with `set_external_field_at!`. 
 
 
 ### Structure factor rewrite

examples/02_LSWT_CoRh2O4.jl

@@ -32,7 +32,7 @@ view_crystal(cryst)
 
 latsize = (1, 1, 1)
 S = 3/2
-J = 7.5413*meV_per_K # (~ 0.65 meV)
+J = 0.63 # (meV)
 sys = System(cryst, latsize, [SpinInfo(1; S, g=2)], :dipole; seed=0)
 set_exchange!(sys, J, Bond(1, 3, [0,0,0]))
 

examples/03_LLD_CoRh2O4.jl

@@ -11,6 +11,8 @@ using Sunny, GLMakie, Statistics
 a = 8.5031 # (Å)
 latvecs = lattice_vectors(a, a, a, 90, 90, 90)
 cryst = Crystal(latvecs, [[0,0,0]], 227, setting="1")
+
+units = Units(:meV)
 latsize = (2, 2, 2)
 S = 3/2
 J = 0.63 # (meV)
@@ -29,7 +31,7 @@ sys = resize_supercell(sys, (10, 10, 10))
 # This dynamics involves a dimensionless `damping` magnitude and target
 # temperature `kT` for thermal fluctuations.
 
-kT = 16 * meV_per_K  # 16K, a temperature slightly below ordering
+kT = 16 * units.K  # 16 K ≈ 1.38 meV, slightly below ordering temperature
 langevin = Langevin(; damping=0.2, kT)
 
 # Use [`suggest_timestep`](@ref) to select an integration timestep for the given

examples/04_GSD_FeI2.jl

@@ -44,6 +44,7 @@ positions = [[0,0,0], [1/3, 2/3, 1/4], [2/3, 1/3, 3/4]]#hide
 types = ["Fe", "I", "I"]#hide
 FeI2 = Crystal(latvecs, positions; types)#hide
 cryst = subcrystal(FeI2, "Fe")#hide
+units = Units(:meV)#hide
 sys = System(cryst, (4,4,4), [SpinInfo(1,S=1,g=2)], :SUN, seed=2)#hide
 J1pm   = -0.236#hide
 J1pmpm = -0.161#hide
@@ -133,10 +134,9 @@ sys_large = resize_supercell(sys, (16,16,4)) # 16x16x4 copies of the original un
 plot_spins(sys_large; color=[s[3] for s in sys_large.dipoles])
 
 # Now we will re-thermalize the system to a configuration just above the
-# ordering temperature. Sunny expects energies in meV by default, so
-# we use `meV_per_K` to convert from kelvin.
+# ordering temperature.
 
-kT = 3.5 * meV_per_K     # 3.5K ≈ 0.30 meV
+kT = 3.5 * units.K # 3.5 K ≈ 0.30 meV
 langevin.kT = kT
 for _ in 1:10_000
     step!(sys_large, langevin)

examples/05_MC_Ising.jl

@@ -23,7 +23,7 @@ crystal = Crystal(latvecs, [[0,0,0]])
 # Ising conventions, select $g=-1$ so that the Zeeman coupling between external
 # field $𝐁$ and spin dipole $𝐬$ is $-𝐁⋅𝐬$.
 L = 128
-sys = System(crystal, (L,L,1), [SpinInfo(1, S=1, g=-1)], :dipole, units=Units.theory, seed=0)
+sys = System(crystal, (L,L,1), [SpinInfo(1, S=1, g=-1)], :dipole, seed=0)
 polarize_spins!(sys, (0,0,1))
 
 # Use [`set_exchange!`](@ref) to include a ferromagnetic Heisenberg interaction
@@ -33,10 +33,9 @@ polarize_spins!(sys, (0,0,1))
 # the symmetries of the square lattice.
 set_exchange!(sys, -1.0, Bond(1,1,(1,0,0)))
 
-# If an external field is desired, it can be set using
-# [`set_external_field!`](@ref).
+# If an external field is desired, it can be set using [`set_field!`](@ref).
 B = 0
-set_external_field!(sys, (0,0,B))
+set_field!(sys, (0, 0, B))
 
 # The critical temperature for the Ising model is known analytically.
 Tc = 2/log(1+√2)

examples/06_CP2_Skyrmions.jl

@@ -35,7 +35,7 @@ cryst = Crystal(lat_vecs, basis_vecs)
 
 L = 40
 dims = (L, L, 1)
-sys = System(cryst, dims, [SpinInfo(1, S=1, g=-1)], :SUN; seed=101, units=Units.theory)
+sys = System(cryst, dims, [SpinInfo(1, S=1, g=-1)], :SUN; seed=101)
 
 # We proceed to implement each term of the Hamiltonian, selecting our parameters
 # so that the system occupies a region of the phase diagram that supports
@@ -57,7 +57,7 @@ set_exchange!(sys, ex2, Bond(1, 1, [1, 2, 0]))
 # Next we add the external field,
 
 h = 15.5
-field = set_external_field!(sys, [0, 0, h])
+field = set_field!(sys, [0, 0, h])
 
 # and finally an easy-plane single-ion anisotropy,
 

examples/07_Dipole_Dipole.jl

@@ -13,9 +13,9 @@ latvecs = lattice_vectors(10.19, 10.19, 10.19, 90, 90, 90)
 positions = [[0, 0, 0]]
 cryst = Crystal(latvecs, positions, 227, setting="2")
 
+units = Units(:meV)
 sys = System(cryst, (1, 1, 1), [SpinInfo(1, S=7/2, g=2)], :dipole, seed=2)
-
-J1 = 0.304 * meV_per_K
+J1 = 0.304 * units.K
 set_exchange!(sys, J1, Bond(1, 2, [0,0,0]))
 
 # Reshape to the primitive cell with four atoms. To facilitate indexing, the
@@ -30,7 +30,7 @@ set_dipole!(sys_prim, [-1, +1, 0], position_to_site(sys_prim, [1/4, 1/4, 0]))
 set_dipole!(sys_prim, [+1, +1, 0], position_to_site(sys_prim, [1/4, 0, 1/4]))
 set_dipole!(sys_prim, [-1, -1, 0], position_to_site(sys_prim, [0, 1/4, 1/4]))
 
-plot_spins(sys_prim; ghost_radius=8, color=[:red,:blue,:yellow,:purple])
+plot_spins(sys_prim; ghost_radius=8, color=[:red, :blue, :yellow, :purple])
 
 # Calculate dispersions with and without long-range dipole interactions. The
 # high-symmetry k-points are specified with respect to the conventional cubic
@@ -67,18 +67,18 @@ ax = Axis(fig[1,1]; xlabel="", ylabel="Energy (K)", xticks, xticklabelrotation=0
 ylims!(ax, 0, 2)
 xlims!(ax, 1, size(disp1, 1))
 for i in axes(disp1, 2)
-    lines!(ax, 1:length(disp1[:,i]), fudge_factor*disp1[:,i]/meV_per_K)
+    lines!(ax, 1:length(disp1[:,i]), fudge_factor*disp1[:,i]/units.K)
 end
 
 ax = Axis(fig[1,2]; xlabel="", ylabel="Energy (K)", xticks, xticklabelrotation=0)
 ylims!(ax, 0.0, 3)
 xlims!(ax, 1, size(disp2, 1))
 for i in axes(disp2, 2)
-    lines!(ax, 1:length(disp2[:,i]), fudge_factor*disp2[:,i]/meV_per_K)
+    lines!(ax, 1:length(disp2[:,i]), fudge_factor*disp2[:,i]/units.K)
 end
 
 for i in axes(disp3, 2)
-    lines!(ax, 1:length(disp3[:,i]), fudge_factor*disp3[:,i]/meV_per_K; color=:gray, linestyle=:dash)
+    lines!(ax, 1:length(disp3[:,i]), fudge_factor*disp3[:,i]/units.K; color=:gray, linestyle=:dash)
 end
 
 fig

examples/08_Momentum_Conventions.jl

@@ -27,9 +27,9 @@ cryst = Crystal(latvecs, [[0,0,0]], "P1")
 
 sys = System(cryst, (1, 1, 25), [SpinInfo(1; S=1, g=2)], :dipole; seed=0)
 D = 0.1 # meV
-B = 10D / (sys.units.μB*2) # ~ 8.64T
+B = 5D # ~ 8.64T
 set_exchange!(sys, dmvec([0, 0, D]), Bond(1, 1, [0, 0, 1]))
-set_external_field!(sys, [0, 0, B])
+set_field!(sys, [0, 0, B])
 
 # The large external field fully polarizes the system. Here, the DM coupling
 # contributes nothing, leaving only Zeeman coupling.

examples/extra/Advanced_MC/PT_WHAM_ising2d.jl

@@ -6,7 +6,7 @@ crystal = Crystal(latvecs, [[0,0,0]])
 
 # 20×20 sites, ferromagnetic exchange
 L = 20
-sys = System(crystal, (L,L,1), [SpinInfo(1, S=1, g=-1)], :dipole, units=Units.theory, seed=0)
+sys = System(crystal, (L,L,1), [SpinInfo(1, S=1, g=-1)], :dipole, seed=0)
 polarize_spins!(sys, (0,0,1))
 set_exchange!(sys, -1.0, Bond(1,1,(1,0,0)))
 

examples/extra/Advanced_MC/REWL_ising2d.jl

@@ -6,7 +6,7 @@ crystal = Crystal(latvecs, [[0,0,0]])
 
 # 20×20 sites, ferromagnetic exchange
 L = 20
-sys = System(crystal, (L,L,1), [SpinInfo(1, S=1, g=-1)], :dipole, units=Units.theory, seed=0)
+sys = System(crystal, (L,L,1), [SpinInfo(1, S=1, g=-1)], :dipole, seed=0)
 polarize_spins!(sys, (0,0,1))
 set_exchange!(sys, -1.0, Bond(1,1,(1,0,0)))
 

examples/extra/Advanced_MC/WL_ising2d.jl

@@ -6,7 +6,7 @@ crystal = Crystal(latvecs, [[0,0,0]])
 
 # 20×20 sites, ferromagnetic exchange
 L = 20
-sys = System(crystal, (L,L,1), [SpinInfo(1, S=1, g=-1)], :dipole, units=Units.theory, seed=0)
+sys = System(crystal, (L,L,1), [SpinInfo(1, S=1, g=-1)], :dipole, seed=0)
 polarize_spins!(sys, (0, 0, 1))
 set_exchange!(sys, -1.0, Bond(1,1,(1,0,0)))
 

src/Reshaping.jl

@@ -81,7 +81,7 @@ function reshape_supercell_aux(sys::System{N}, new_latsize::NTuple{3, Int}, new_
     orig_sys = @something sys.origin sys
 
     new_sys = System(orig_sys, sys.mode, new_cryst, new_latsize, new_Ns, new_κs, new_gs, new_ints, new_ewald,
-        new_extfield, new_dipoles, new_coherents, new_dipole_buffers, new_coherent_buffers, sys.units, copy(sys.rng))
+        new_extfield, new_dipoles, new_coherents, new_dipole_buffers, new_coherent_buffers, copy(sys.rng))
 
     # Transfer interactions. In the case of an inhomogeneous system, the
     # interactions in `sys` have detached from `orig`, so we use the latest
@@ -102,7 +102,7 @@ function reshape_supercell_aux(sys::System{N}, new_latsize::NTuple{3, Int}, new_
     # Restore dipole-dipole interactions if present. This involves pre-computing
     # an interaction matrix that depends on `new_latsize`.
     if !isnothing(sys.ewald)
-        enable_dipole_dipole!(new_sys)
+        enable_dipole_dipole!(new_sys, sys.ewald.μ0_μB²)
     end
 
     return new_sys

src/SpinWaveTheory/HamiltonianDipole.jl

@@ -2,7 +2,7 @@
 function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_reshaped::Vec3)
     (; sys, data) = swt
     (; local_rotations, stevens_coefs, sqrtS) = data
-    (; extfield, gs, units) = sys
+    (; extfield, gs) = sys
 
     L = nbands(swt) 
     @assert size(H) == (2L, 2L)
@@ -17,7 +17,7 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r
 
     for (i, int) in enumerate(sys.interactions_union)
         # Zeeman term
-        B = units.μB * (gs[1, 1, 1, i]' * extfield[1, 1, 1, i]) 
+        B = gs[1, 1, 1, i]' * extfield[1, 1, 1, i]
         B′ = - dot(B, local_rotations[i][:, 3]) / 2
         H11[i, i] += B′
         H22[i, i] += B′
@@ -97,15 +97,14 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r
 
     # Add long-range dipole-dipole
     if !isnothing(sys.ewald)
-        μB² = units.μB^2
         Rs = local_rotations
 
         # Interaction matrix for wavevector q
-        A = precompute_dipole_ewald_at_wavevector(sys.crystal, (1,1,1), units.μ0, q_reshaped)
+        A = precompute_dipole_ewald_at_wavevector(sys.crystal, (1,1,1), q_reshaped) * sys.ewald.μ0_μB²
         A = reshape(A, L, L)
 
         # Interaction matrix for wavevector (0,0,0). It could be recalculated as:
-        # precompute_dipole_ewald(sys.crystal, (1,1,1), units.μ0)
+        # precompute_dipole_ewald(sys.crystal, (1,1,1)) * sys.ewald.μ0_μB²
         A0 = sys.ewald.A
         A0 = reshape(A0, L, L)
 
@@ -114,8 +113,8 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r
             # An ordered pair of magnetic moments contribute (μᵢ A μⱼ)/2 to the
             # energy. A symmetric contribution will appear for the bond reversal
             # (i, j) → (j, i).  Note that μ = -μB g S.
-            J = μB² * gs[i]' * A[i, j] * gs[j] / 2
-            J0 = μB² * gs[i]' * A0[i, j] * gs[j] / 2
+            J = gs[i]' * A[i, j] * gs[j] / 2
+            J0 = gs[i]' * A0[i, j] * gs[j] / 2
 
             # Perform same transformation as appears in usual bilinear exchange.
             # Rⱼ denotes a rotation from ẑ to the ground state dipole Sⱼ.
@@ -179,14 +178,14 @@ function multiply_by_hamiltonian_dipole!(y::Array{ComplexF64, 2}, x::Array{Compl
 
     # Add Zeeman and single-ion anisotropy. These entries are q-independent and
     # could be precomputed.
-    (; units, extfield, gs) = sys
+    (; extfield, gs) = sys
     for i in 1:L
         (; c2, c4, c6) = stevens_coefs[i]
         S = sqrtS[i]^2
         A1 = -3S*c2[3] - 40*S^3*c4[5] - 168*S^5*c6[7]
         A2 = im*(S*c2[5] + 6S^3*c4[7] + 16S^5*c6[9]) + (S*c2[1] + 6S^3*c4[3] + 16S^5*c6[5])
 
-        B = units.μB * (gs[1, 1, 1, i]' * extfield[1, 1, 1, i]) 
+        B = gs[1, 1, 1, i]' * extfield[1, 1, 1, i]
         B′ = - dot(B, view(local_rotations[i], :, 3)) / 2 
 
         # Seems to be no benefit to breaking this into two loops acting on

src/SpinWaveTheory/HamiltonianSUN.jl

@@ -3,7 +3,7 @@
 function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_reshaped::Vec3)
     (; sys, data) = swt
     (; spins_localized) = data
-    (; gs, units) = sys
+    (; gs) = sys
 
     N = sys.Ns[1]
     Na = natoms(sys.crystal)
@@ -70,10 +70,9 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh
 
     if !isnothing(sys.ewald)
         N = sys.Ns[1]
-        μB² = units.μB^2
 
         # Interaction matrix for wavevector q
-        A = precompute_dipole_ewald_at_wavevector(sys.crystal, (1,1,1), units.μ0, q_reshaped)
+        A = precompute_dipole_ewald_at_wavevector(sys.crystal, (1,1,1), q_reshaped) * sys.ewald.μ0_μB²
         A = reshape(A, Na, Na)
 
         # Interaction matrix for wavevector (0,0,0)
@@ -82,10 +81,10 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh
 
         for i in 1:Na, j in 1:Na
             # An ordered pair of magnetic moments contribute (μᵢ A μⱼ)/2 to the
-            # energy. A symmetric contribution will appear for the bond reversal
-            # (i, j) → (j, i).  Note that μ = -μB g S.
-            J = μB² * gs[i]' * A[i, j] * gs[j] / 2
-            J0 = μB² * gs[i]' * A0[i, j] * gs[j] / 2
+            # energy, where μ = - g S. A symmetric contribution will appear for
+            # the bond reversal (i, j) → (j, i).
+            J = gs[i]' * A[i, j] * gs[j] / 2
+            J0 = gs[i]' * A0[i, j] * gs[j] / 2
 
             for α in 1:3, β in 1:3
                 Ai = view(spins_localized, :, :, α, i)

src/SpinWaveTheory/SpinWaveTheory.jl

@@ -150,7 +150,7 @@ function swt_data_sun(sys::System{N}, obs) where N
 
         # Accumulate Zeeman terms into OnsiteCoupling
         S = spin_matrices_of_dim(; N)
-        B = sys.units.μB * (sys.gs[atom]' * sys.extfield[atom])
+        B = sys.gs[atom]' * sys.extfield[atom]
         int.onsite += B' * S
 
         # Rotate onsite anisotropy

src/Spiral/LuttingerTisza.jl

@@ -22,11 +22,11 @@ function luttinger_tisza_exchange(sys::System; k, ϵ=0)
     end
 
     if !isnothing(sys.ewald)
-        A = Sunny.precompute_dipole_ewald_at_wavevector(sys.crystal, (1,1,1), sys.units.μ0, k)
+        A = Sunny.precompute_dipole_ewald_at_wavevector(sys.crystal, (1,1,1), k) * sys.ewald.μ0_μB²
         A = reshape(A, Na, Na)
         for i in 1:Na, j in 1:Na
-            J_k[:, i, :, j] += μB² * gs[i]' * A[i, j] * gs[j] / 2
-        end    
+            J_k[:, i, :, j] += gs[i]' * A[i, j] * gs[j] / 2
+        end
     end
 
     J_k = reshape(J_k, 3*Na, 3*Na)

src/Spiral/SpiralEnergy.jl

@@ -70,22 +70,22 @@ function spiral_energy_and_gradient_aux!(dEds, sys::System{0}; k, axis)
         E += E_aniso
 
         # Zeeman coupling
-        E += sys.extfield[i]' * (sys.units.μB * sys.gs[i] * Si)
+        E += sys.extfield[i]' * (sys.gs[i] * Si)
 
         if accum_grad
             dEds[i] += dEds_aniso
-            dEds[i] += sys.units.μB * sys.gs[i]' * sys.extfield[i]
+            dEds[i] += sys.gs[i]' * sys.extfield[i]
         end
     end
 
     # See "spiral_energy.lyx" for derivation
     if !isnothing(sys.ewald)
-        μ = [sys.units.μB*magnetic_moment(sys, site) for site in eachsite(sys)]
+        μ = [magnetic_moment(sys, site) for site in eachsite(sys)]
 
         A0 = sys.ewald.A
         A0 = reshape(A0, Na, Na)
 
-        Ak = Sunny.precompute_dipole_ewald_at_wavevector(sys.crystal, (1,1,1), sys.units.μ0, k)
+        Ak = Sunny.precompute_dipole_ewald_at_wavevector(sys.crystal, (1,1,1), k) * sys.ewald.μ0_μB²
         Ak = reshape(Ak, Na, Na)
 
         ϵ = 1e-8

src/Spiral/SpiralSWT.jl

@@ -58,9 +58,8 @@ function swt_hamiltonian_dipole_spiral!(H::Matrix{ComplexF64}, swt::SpinWaveTheo
     H[:,:] = H / 2
 
     # Add Zeeman term
-    (; extfield, gs, units) = sys
     for i in 1:L
-        B = units.μB * (Transpose(extfield[1, 1, 1, i]) * gs[1, 1, 1, i]) 
+        B = sys.extfield[1, 1, 1, i]' * sys.gs[1, 1, 1, i]
         B′ = - (B * local_rotations[i][:, 3]) / 2 
         H[i, i]     += B′
         H[i+L, i+L] += conj(B′)