diff --git a/docs/src/versions.md b/docs/src/versions.md
index 0a47a0874..70af906cd 100644
--- a/docs/src/versions.md
+++ b/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
diff --git a/examples/02_LSWT_CoRh2O4.jl b/examples/02_LSWT_CoRh2O4.jl
index e73ebf87d..1f8078661 100644
--- a/examples/02_LSWT_CoRh2O4.jl
+++ b/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]))
 
diff --git a/examples/03_LLD_CoRh2O4.jl b/examples/03_LLD_CoRh2O4.jl
index 2d5fb8f52..cad40b327 100644
--- a/examples/03_LLD_CoRh2O4.jl
+++ b/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
diff --git a/examples/04_GSD_FeI2.jl b/examples/04_GSD_FeI2.jl
index 712cbe36b..8cecde105 100644
--- a/examples/04_GSD_FeI2.jl
+++ b/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)
diff --git a/examples/05_MC_Ising.jl b/examples/05_MC_Ising.jl
index 900ad2682..0ff8cf101 100644
--- a/examples/05_MC_Ising.jl
+++ b/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)
diff --git a/examples/06_CP2_Skyrmions.jl b/examples/06_CP2_Skyrmions.jl
index 9ba49e00a..c74e574dd 100644
--- a/examples/06_CP2_Skyrmions.jl
+++ b/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,
 
diff --git a/examples/07_Dipole_Dipole.jl b/examples/07_Dipole_Dipole.jl
index ea3002395..71af85656 100644
--- a/examples/07_Dipole_Dipole.jl
+++ b/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
diff --git a/examples/08_Momentum_Conventions.jl b/examples/08_Momentum_Conventions.jl
index e50746fad..3b807910d 100644
--- a/examples/08_Momentum_Conventions.jl
+++ b/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.
diff --git a/examples/extra/Advanced_MC/PT_WHAM_ising2d.jl b/examples/extra/Advanced_MC/PT_WHAM_ising2d.jl
index e4b3c04b6..593836b03 100644
--- a/examples/extra/Advanced_MC/PT_WHAM_ising2d.jl
+++ b/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)))
 
diff --git a/examples/extra/Advanced_MC/REWL_ising2d.jl b/examples/extra/Advanced_MC/REWL_ising2d.jl
index 1ecc31af9..b424060db 100644
--- a/examples/extra/Advanced_MC/REWL_ising2d.jl
+++ b/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)))
 
diff --git a/examples/extra/Advanced_MC/WL_ising2d.jl b/examples/extra/Advanced_MC/WL_ising2d.jl
index aefdda4cb..139ff67a7 100644
--- a/examples/extra/Advanced_MC/WL_ising2d.jl
+++ b/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)))
 
diff --git a/src/Reshaping.jl b/src/Reshaping.jl
index 082c89b95..c6701c767 100644
--- a/src/Reshaping.jl
+++ b/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
diff --git a/src/SpinWaveTheory/HamiltonianDipole.jl b/src/SpinWaveTheory/HamiltonianDipole.jl
index 009f72386..98a813885 100644
--- a/src/SpinWaveTheory/HamiltonianDipole.jl
+++ b/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
diff --git a/src/SpinWaveTheory/HamiltonianSUN.jl b/src/SpinWaveTheory/HamiltonianSUN.jl
index 77d3fcb4d..cce78442c 100644
--- a/src/SpinWaveTheory/HamiltonianSUN.jl
+++ b/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)
diff --git a/src/SpinWaveTheory/SpinWaveTheory.jl b/src/SpinWaveTheory/SpinWaveTheory.jl
index 70df3571a..60f4f3b5d 100644
--- a/src/SpinWaveTheory/SpinWaveTheory.jl
+++ b/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
diff --git a/src/Spiral/LuttingerTisza.jl b/src/Spiral/LuttingerTisza.jl
index 8c4530540..2f4c81054 100644
--- a/src/Spiral/LuttingerTisza.jl
+++ b/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)
diff --git a/src/Spiral/SpiralEnergy.jl b/src/Spiral/SpiralEnergy.jl
index 2823797d0..50e4d389d 100644
--- a/src/Spiral/SpiralEnergy.jl
+++ b/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
diff --git a/src/Spiral/SpiralSWT.jl b/src/Spiral/SpiralSWT.jl
index d88a04462..7466a3263 100644
--- a/src/Spiral/SpiralSWT.jl
+++ b/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′)
diff --git a/src/Sunny.jl b/src/Sunny.jl
index 774ff15d7..9905299f2 100644
--- a/src/Sunny.jl
+++ b/src/Sunny.jl
@@ -44,7 +44,7 @@ export Crystal, subcrystal, standardize, lattice_vectors, lattice_params, primit
     reference_bonds, print_site, print_bond, print_symmetry_table, print_suggested_frame
 
 include("Units.jl")
-export meV_per_K, Units
+export Units
 
 include("System/SpinInfo.jl")
 include("System/Types.jl")
@@ -56,7 +56,7 @@ include("System/Interactions.jl")
 export SpinInfo, System, Site, clone_system, eachsite, position_to_site, global_position, magnetic_moment,
     set_coherent!, set_dipole!, polarize_spins!, randomize_spins!, set_spin_rescaling!, energy, energy_per_site,
     spin_label, set_onsite_coupling!, set_pair_coupling!, set_exchange!, dmvec, enable_dipole_dipole!,
-    set_external_field!, to_inhomogeneous, set_external_field_at!, set_vacancy_at!, set_onsite_coupling_at!,
+    set_field!, to_inhomogeneous, set_field_at!, set_vacancy_at!, set_onsite_coupling_at!,
     set_exchange_at!, set_pair_coupling_at!, symmetry_equivalent_bonds, remove_periodicity!,
     modify_exchange_with_truncated_dipole_dipole!
 
@@ -125,6 +125,7 @@ include("MonteCarlo/ParallelWangLandau.jl")
 export propose_uniform, propose_flip, propose_delta, @mix_proposals, LocalSampler
 
 include("deprecated.jl")
+export set_external_field!, set_external_field_at!, meV_per_K
 
 isloaded(pkg::String) = any(k -> k.name == pkg, keys(Base.loaded_modules))
 
diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl
index 1b738742e..f3cc4f634 100644
--- a/src/System/Ewald.jl
+++ b/src/System/Ewald.jl
@@ -1,8 +1,8 @@
 
-function Ewald(sys::System{N}) where N
-    (; crystal, latsize, units) = sys
+function Ewald(sys::System{N}, μ0_μB²) where N
+    (; crystal, latsize) = sys
 
-    A = precompute_dipole_ewald(crystal, latsize, units.μ0)
+    A = precompute_dipole_ewald(crystal, latsize) * μ0_μB²
 
     na = natoms(crystal)
     μ = zeros(Vec3, latsize..., na)
@@ -18,7 +18,7 @@ function Ewald(sys::System{N}) where N
     plan     = FFTW.plan_rfft(mock_spins, 2:4; flags=FFTW.MEASURE)
     ift_plan = FFTW.plan_irfft(Fμ, latsize[1], 2:4; flags=FFTW.MEASURE)
 
-    return Ewald(A, μ, ϕ, FA, Fμ, Fϕ, plan, ift_plan)
+    return Ewald(μ0_μB², A, μ, ϕ, FA, Fμ, Fϕ, plan, ift_plan)
 end
 
 # Ideally, this would clone all mutable state within Ewald. Note that `A`, `FA`
@@ -29,23 +29,23 @@ end
 # https://github.com/JuliaMath/FFTW.jl/issues/261.
 function clone_ewald(ewald::Ewald)
     error("Not supported")
-    (; A, μ, ϕ, FA, Fμ, Fϕ, plan, ift_plan) = ewald
-    return Ewald(A, copy(μ), copy(ϕ), FA, copy(Fμ), copy(Fϕ), copy(plan), copy(ift_plan))
+    (; μ0_μB², A, μ, ϕ, FA, Fμ, Fϕ, plan, ift_plan) = ewald
+    return Ewald(μ0_μB², A, copy(μ), copy(ϕ), FA, copy(Fμ), copy(Fϕ), copy(plan), copy(ift_plan))
 end
 
 # Tensor product of 3-vectors
 (⊗)(a::Vec3,b::Vec3) = reshape(kron(a,b), 3, 3)
 
 
-function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int}, μ0)
-    precompute_dipole_ewald_aux(cryst, latsize, μ0, Vec3(0,0,0), cos, Val{Float64}())
+function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int})
+    precompute_dipole_ewald_aux(cryst, latsize, Vec3(0,0,0), cos, Val{Float64}())
 end
 
-function precompute_dipole_ewald_at_wavevector(cryst::Crystal, latsize::NTuple{3,Int}, μ0, q_reshaped)
-    precompute_dipole_ewald_aux(cryst, latsize, μ0, q_reshaped, cis, Val{ComplexF64}())
+function precompute_dipole_ewald_at_wavevector(cryst::Crystal, latsize::NTuple{3,Int}, q_reshaped)
+    precompute_dipole_ewald_aux(cryst, latsize, q_reshaped, cis, Val{ComplexF64}())
 end
 
-# Precompute the pairwise interaction matrix A beween magnetic moments μ. For
+# Precompute the pairwise interaction matrix A between magnetic moments μ. For
 # q_reshaped = 0, this yields the usual Ewald energy, E = μᵢ Aᵢⱼ μⱼ / 2. Nonzero
 # q_reshaped is useful in spin wave theory. Physically, this amounts to a
 # modification of the periodic boundary conditions, such that μ(q) can be
@@ -58,7 +58,7 @@ end
 # part cancels in the symmetric sum over ±k. Specifically, replace `cis(x) ≡
 # exp(i x) = cos(x) + i sin(x)` with just `cos(x)` for efficiency. The parameter
 # `T ∈ {Float64, ComplexF64}` controls the return type in a type-stable way. 
-function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, μ0, q_reshaped, cis, ::Val{T}) where T
+function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, q_reshaped, cis, ::Val{T}) where T
     na = natoms(cryst)
     A = zeros(SMatrix{3, 3, T, 9}, latsize..., na, na)
 
@@ -107,7 +107,7 @@ function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, μ0
                 erfc0 = erfc(r/(√2*σ))
                 gauss0 = √(2/π) * (r/σ) * exp(-r²/2σ²)    
                 phase = cis(2π * dot(q_reshaped, n))
-                acc += phase * (μ0/4π) * ((I₃/r³) * (erfc0 + gauss0) - (3(rhat⊗rhat)/r³) * (erfc0 + (1+r²/3σ²) * gauss0))
+                acc += phase * (1/4π) * ((I₃/r³) * (erfc0 + gauss0) - (3(rhat⊗rhat)/r³) * (erfc0 + (1+r²/3σ²) * gauss0))
             end
         end
 
@@ -126,14 +126,14 @@ function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, μ0
                 # acc += (μ0/2V) * I₃
             elseif ϵ² < k² <= kmax*kmax
                 phase = cis(-k⋅Δr)
-                acc += phase * (μ0/V) * (exp(-σ²*k²/2) / k²) * (k⊗k)
+                acc += phase * (1/V) * (exp(-σ²*k²/2) / k²) * (k⊗k)
             end
         end
 
         #####################################################
         ## Remove self energies
         if iszero(Δr)
-            acc += - μ0*I₃/(3(2π)^(3/2)*σ³)
+            acc += - I₃/(3(2π)^(3/2)*σ³)
         end
 
         # For sites site1=(cell1, i) and site2=(cell2, j) offset by an amount
@@ -154,7 +154,7 @@ function ewald_energy(sys::System{N}) where N
     even_rft_size = latsize[1] % 2 == 0
 
     for site in eachsite(sys)
-        μ[site] = sys.units.μB*magnetic_moment(sys, site)
+        μ[site] = magnetic_moment(sys, site)
     end
 
     E = 0.0
@@ -182,12 +182,12 @@ end
 # Use FFT to accumulate the entire field dE/ds for long-range dipole-dipole
 # interactions
 function accum_ewald_grad!(∇E, dipoles, sys::System{N}) where N
-    (; ewald, units, gs) = sys
+    (; ewald, gs) = sys
     (; μ, FA, Fμ, Fϕ, ϕ, plan, ift_plan) = ewald
 
     # Fourier transformed magnetic moments
     for site in eachsite(sys)
-        μ[site] = units.μB * gs[site] * dipoles[site]
+        μ[site] = gs[site] * dipoles[site]
     end
     mul!(Fμ, plan, reinterpret(reshape, Float64, μ))
 
@@ -205,20 +205,20 @@ function accum_ewald_grad!(∇E, dipoles, sys::System{N}) where N
     mul!(ϕr, ift_plan, Fϕ)
     
     for site in eachsite(sys)
-        ∇E[site] += units.μB * (gs[site]' * ϕ[site])
+        ∇E[site] += gs[site]' * ϕ[site]
     end
 end
 
 # Calculate the field dE/ds at site1 generated by a dipole at site2.
 function ewald_pairwise_grad_at(sys::System{N}, site1, site2) where N
-    (; gs, ewald, units) = sys
+    (; gs, ewald) = sys
     latsize = size(ewald.ϕ)[1:3]
     cell_offset = mod.(Tuple(to_cell(site2)-to_cell(site1)), latsize)
     cell = CartesianIndex(cell_offset .+ (1,1,1))
 
     # The factor of 1/2 in the energy formula `E = μ (A ⋆ μ) / 2` disappears due
-    # to quadratic appearance of μ = - μB g S.
-    return units.μB^2 * gs[site1]' * ewald.A[cell, to_atom(site1), to_atom(site2)] * gs[site2] * sys.dipoles[site2]
+    # to quadratic appearance of μ = - g S.
+    return gs[site1]' * ewald.A[cell, to_atom(site1), to_atom(site2)] * gs[site2] * sys.dipoles[site2]
 end
 
 # Calculate the field dE/ds at `site` generated by all `dipoles`.
@@ -233,39 +233,40 @@ end
 # Calculate the change in dipole-dipole energy when the spin at `site` is
 # updated to `s`
 function ewald_energy_delta(sys::System{N}, site, s::Vec3) where N
-    (; dipoles, ewald, units) = sys
+    (; dipoles, ewald) = sys
     Δs = s - dipoles[site]
-    Δμ = units.μB * (sys.gs[site] * Δs)
+    Δμ = sys.gs[site] * Δs
     i = to_atom(site)
     ∇E = ewald_grad_at(sys, site)
     return Δs⋅∇E + dot(Δμ, ewald.A[1, 1, 1, i, i], Δμ) / 2
 end
 
 """
-    modify_exchange_with_truncated_dipole_dipole!(sys::System, cutoff)
-
-This *experimental* function is subject to change.
+    modify_exchange_with_truncated_dipole_dipole!(sys::System, cutoff, μ0_μB²)
 
 Like [`enable_dipole_dipole!`](@ref), the purpose of this function is to
 introduce long-range dipole-dipole interactions between magnetic moments.
 Whereas `enable_dipole_dipole!` employs Ewald summation, this function instead
-employs real-space pair couplings, with truncation at some user-specified
-`cutoff` distance. If the cutoff is relatively small, then this function can be
-faster than `enable_dipole_dipole!`, especially for spin wave theory
-calculations.
-
-**Caution**. This function will modify existing bilinear couplings between spins
-by adding dipole-dipole interactions. It must therefore be called _after_ all
-other pair couplings have been specified. Conversely, any calls to
-`set_exchange!`, `set_pair_coupling!`, etc. will irreversibly delete the
-dipole-dipole interactions that have been introduced by this function.
+employs real-space pair couplings with truncation at the specified `cutoff`
+distance. If the cutoff is relatively small, then this function may be faster
+than `enable_dipole_dipole!`.
+
+!!! warning "Mutation of existing couplings"
+
+    This function will modify existing bilinear couplings between spins by
+    adding dipole-dipole interactions. It must therefore be called _after_
+    all other pair couplings have been specified. Conversely, any calls to
+    `set_exchange!`, `set_pair_coupling!`, etc. will irreversibly delete the
+    dipole-dipole interactions that have been introduced by this function.
 """
-function modify_exchange_with_truncated_dipole_dipole!(sys::System{N}, cutoff) where N
-    (; gs) = sys
-    (; μB, μ0) = sys.units
+function modify_exchange_with_truncated_dipole_dipole!(sys::System{N}, cutoff, μ0_μB²=nothing) where N
+    if isnothing(μ0_μB²)
+        @warn "Deprecated syntax! Consider `modify_exchange_with_truncated_dipole_dipole!(sys, cutoff, units.vacuum_permeability)` where `units = Units(:meV)`."
+        μ0_μB² = Units(:meV).vacuum_permeability
+    end
 
     if !isnothing(sys.origin)
-        modify_exchange_with_truncated_dipole_dipole!(sys.origin, cutoff)
+        modify_exchange_with_truncated_dipole_dipole!(sys.origin, cutoff, μ0_μB²)
         transfer_interactions!(sys, sys.origin)
         return
     end
@@ -280,7 +281,7 @@ function modify_exchange_with_truncated_dipole_dipole!(sys::System{N}, cutoff) w
                 r = global_displacement(sys.crystal, bond′)
                 iszero(r) && continue
                 r̂ = normalize(r)
-                bilin = (μ0/4π) * μB^2 * gs[i]' * ((I - 3r̂⊗r̂) / norm(r)^3) * gs[j]
+                bilin = (μ0_μB²/4π) * sys.gs[i]' * ((I - 3r̂⊗r̂) / norm(r)^3) * sys.gs[j]
                 replace_coupling!(ints[i].pair, PairCoupling(bond′, 0.0, Mat3(bilin), 0.0, zero(TensorDecomposition)); accum=true)
             end
         end
diff --git a/src/System/Interactions.jl b/src/System/Interactions.jl
index ec6f937d4..d3f54b005 100644
--- a/src/System/Interactions.jl
+++ b/src/System/Interactions.jl
@@ -64,7 +64,7 @@ end
 
 
 """
-    enable_dipole_dipole!(sys::System)
+    enable_dipole_dipole!(sys::System, μ0_μB²)
 
 Enables long-range interactions between magnetic dipole moments,
 
@@ -73,39 +73,71 @@ Enables long-range interactions between magnetic dipole moments,
 ```
 
 where the sum is over all pairs of spins (singly counted), including periodic
-images, regularized using the Ewald summation convention. See
-[`magnetic_moment`](@ref) for the relationship between ``μ_i`` and the spin
-angular momentum. The vacuum permeability ``μ_0`` is a physical constant
-determined by the system of [`Units`](@ref).
+images, regularized using the Ewald summation convention. The
+[`magnetic_moment`](@ref) is defined as ``μ = -g μ_B 𝐒``, where ``𝐒`` is the
+spin angular momentum dipole. The parameter `μ0_μB²` specifies the physical
+constant ``μ_0 μ_B^2``, which has dimensions of length³-energy. Obtain this
+constant for a given system of [`Units`](@ref) via its `vacuum_permeability`
+property.
+
+# Example
+
+```julia
+# Valid for a system with lengths in Å and energies in meV
+units = Units(:meV)
+enable_dipole_dipole!(sys, units.vacuum_permeability)
+```
+
+See also [`modify_exchange_with_truncated_dipole_dipole!`](@ref).
 """
-function enable_dipole_dipole!(sys::System{N}) where N
-    isnan(sys.units.μ0) && error("Dipole-dipole interactions incompatible with Units.theory")
-    sys.ewald = Ewald(sys)
+function enable_dipole_dipole!(sys::System{N}, μ0_μB²=nothing) where N
+    if isnothing(μ0_μB²)
+        @warn "Deprecated syntax! Consider `enable_dipole_dipole!(sys, units.vacuum_permeability)` where `units = Units(:meV)`."
+        μ0_μB² = Units(:meV).vacuum_permeability
+    end
+    sys.ewald = Ewald(sys, μ0_μB²)
     return
 end
 
 """
-    set_external_field!(sys::System, B::Vec3)
-
-Sets the external field ``𝐁`` that couples to all magnetic moments, ``- ∑_i
-𝐁⋅μ_i``. See [`magnetic_moment`](@ref) for the relationship between ``μ_i`` and
-the spin angular momentum.
+    set_field!(sys::System, B_μB)
+
+Sets the external magnetic field ``𝐁`` scaled by the Bohr magneton ``μ_B``.
+This scaled field has units of energy and couples directly to the dimensionless
+[`magnetic_moment`](@ref). At every site, the Zeeman coupling contributes an
+energy ``+ (𝐁 μ_B) ⋅ (g 𝐒)``, involving the local ``g``-tensor and spin
+angular momentum ``𝐒``. Commonly, ``g ≈ +2`` such that ``𝐒`` is favored to
+anti-align with the applied field ``𝐁``. Note that a given system of
+[`Units`](@ref) will implicitly use the Bohr magneton to convert between field
+and energy dimensions.
+
+# Example
+
+```julia
+# In units of meV, apply a 2 tesla field in the z-direction
+units = Units(:meV)
+set_field!(sys, [0, 0, 2] * units.T)
+```
 """
-function set_external_field!(sys::System, B)
+function set_field!(sys::System, B_μB)
     for site in eachsite(sys)
-        set_external_field_at!(sys, B, site)
+        set_field_at!(sys, B_μB, site)
     end
 end
 
 """
-    set_external_field_at!(sys::System, B::Vec3, site::Site)
+    set_field_at!(sys::System, B_μB, site::Site)
+
+Sets the external magnetic field ``𝐁`` scaled by the Bohr magneton ``μ_B`` for
+a single [`Site`](@ref). This scaled field has units of energy and couples
+directly to the dimensionless [`magnetic_moment`](@ref). Note that a given
+system of [`Units`](@ref) will implicitly use the Bohr magneton to convert
+between field and energy dimensions.
 
-Sets a local field ``𝐁`` that couples to a single magnetic moment, ``-𝐁⋅μ_i``.
-See [`magnetic_moment`](@ref) for the relationship between ``μ_i`` and the spin
-angular momentum. [`Site`](@ref) includes a unit cell and a sublattice index.
+See the documentation of [`set_field!`](@ref) for more information.
 """
-function set_external_field_at!(sys::System, B, site)
-    sys.extfield[to_cartesian(site)] = Vec3(B)
+function set_field_at!(sys::System, B_μB, site)
+    sys.extfield[to_cartesian(site)] = Vec3(B_μB)
 end
 
 """
@@ -140,7 +172,7 @@ function local_energy_change(sys::System{N}, site, state::SpinState) where N
     ΔE = 0.0
 
     # Zeeman coupling to external field
-    ΔE += sys.units.μB * dot(extfield[site], sys.gs[site], Δs)
+    ΔE += dot(extfield[site], sys.gs[site], Δs)
 
     # Single-ion anisotropy, dipole or SUN mode
     if N == 0
@@ -217,7 +249,7 @@ function energy(sys::System{N}) where N
 
     # Zeeman coupling to external field
     for site in eachsite(sys)
-        E += sys.units.μB * sys.extfield[site] ⋅ (sys.gs[site] * sys.dipoles[site])
+        E += sys.extfield[site] ⋅ (sys.gs[site] * sys.dipoles[site])
     end
 
     # Anisotropies and exchange interactions
@@ -328,7 +360,7 @@ function set_energy_grad_dipoles!(∇E, dipoles::Array{Vec3, 4}, sys::System{N})
 
     # Zeeman coupling
     for site in eachsite(sys)
-        ∇E[site] += sys.units.μB * (sys.gs[site]' * sys.extfield[site])
+        ∇E[site] += sys.gs[site]' * sys.extfield[site]
     end
 
     # Anisotropies and exchange interactions
diff --git a/src/System/SpinInfo.jl b/src/System/SpinInfo.jl
index d6663fa39..92f30cae9 100644
--- a/src/System/SpinInfo.jl
+++ b/src/System/SpinInfo.jl
@@ -4,7 +4,7 @@
 Characterizes the spin at a given `atom` index within the crystal unit cell. `S`
 is an integer multiple of 1/2 and gives the spin angular momentum in units of ħ.
 `g` is the g-factor or tensor, such that an angular momentum dipole ``s``
-produces a magnetic moment ``g s`` in units of the Bohr magneton.
+produces a magnetic moment ``+ g s`` in units of the Bohr magneton.
 """
 struct SpinInfo
     atom   :: Int     # Index of atom in unit cell
diff --git a/src/System/System.jl b/src/System/System.jl
index ec68fad88..8b0e8d7e8 100644
--- a/src/System/System.jl
+++ b/src/System/System.jl
@@ -1,10 +1,10 @@
 """
-    System(crystal::Crystal, latsize, infos, mode; units=Units.meV, seed::Int)
+    System(crystal::Crystal, latsize, infos, mode; seed::Int)
 
 Construct a `System` of spins for a given [`Crystal`](@ref) symmetry. The
 `latsize` parameter determines the number of unit cells in each lattice vector
-direction. The `infos` parameter is a list of [`SpinInfo`](@ref) objects, which
-determine the magnitude ``S`` and ``g``-tensor of each spin.
+direction. The `infos` parameter is a list of [`SpinInfo`](@ref) objects, one
+for each symmetry-inequivalent sublattice.
 
 The two primary options for `mode` are `:SUN` and `:dipole`. In the former, each
 spin-``S`` degree of freedom is described as an SU(_N_) coherent state, i.e. a
@@ -18,16 +18,17 @@ be renormalized to improve accuracy. To disable this renormalization, use the
 mode `:dipole_large_S` which applies the ``S → ∞`` classical limit. For details,
 see the documentation page: [Interaction Strength Renormalization](@ref).
 
-The default units system of (meV, Å, tesla) can be overridden by with the
-`units` parameter; see [`Units`](@ref). 
-
 An optional `seed` may be provided to achieve reproducible random number
 generation.
 
 All spins are initially polarized in the ``z``-direction.
 """
 function System(crystal::Crystal, latsize::NTuple{3,Int}, infos::Vector{SpinInfo}, mode::Symbol;
-                units=Units.meV, seed=nothing)
+                seed=nothing, units=nothing)
+    if !isnothing(units)
+        @warn "units argument to System is deprecated and will be ignored!"
+    end
+    
     if !in(mode, (:SUN, :dipole, :dipole_large_S))
         error("Mode must be `:SUN`, `:dipole`, or `:dipole_large_S`.")
     end
@@ -74,7 +75,7 @@ function System(crystal::Crystal, latsize::NTuple{3,Int}, infos::Vector{SpinInfo
     rng = isnothing(seed) ? Random.Xoshiro() : Random.Xoshiro(seed)
 
     ret = System(nothing, mode, crystal, latsize, Ns, κs, gs, interactions, ewald,
-                 extfield, dipoles, coherents, dipole_buffers, coherent_buffers, units, rng)
+                 extfield, dipoles, coherents, dipole_buffers, coherent_buffers, rng)
     polarize_spins!(ret, (0,0,1))
     return ret
 end
@@ -132,7 +133,7 @@ Creates a full clone of the system, such that mutable updates to one copy will
 not affect the other, and thread safety is guaranteed.
 """
 function clone_system(sys::System{N}) where N
-    (; origin, mode, crystal, latsize, Ns, gs, κs, extfield, interactions_union, ewald, dipoles, coherents, units, rng) = sys
+    (; origin, mode, crystal, latsize, Ns, gs, κs, extfield, interactions_union, ewald, dipoles, coherents, rng) = sys
 
     origin_clone = isnothing(origin) ? nothing : clone_system(origin)
     ewald_clone = nothing # TODO: use clone_ewald(ewald)
@@ -147,13 +148,13 @@ function clone_system(sys::System{N}) where N
 
     ret = System(origin_clone, mode, crystal, latsize, Ns, copy(κs), copy(gs),
                  interactions_clone, ewald_clone, copy(extfield), copy(dipoles), copy(coherents),
-                 empty_dipole_buffers, empty_coherent_buffers, units, copy(rng))
+                 empty_dipole_buffers, empty_coherent_buffers, copy(rng))
 
     if !isnothing(ewald)
         # At the moment, clone_ewald is unavailable, so instead rebuild the
         # Ewald data structures from scratch. This might be fixed eventually.
         # See https://github.com/JuliaMath/FFTW.jl/issues/261.
-        enable_dipole_dipole!(ret)
+        enable_dipole_dipole!(ret, ewald.μ0_μB²)
     end
 
     return ret
@@ -169,7 +170,7 @@ atom within the unit cell).
 
 This object can be used to index `dipoles` and `coherents` fields of a `System`.
 A `Site` is also required to specify inhomogeneous interactions via functions
-such as [`set_external_field_at!`](@ref) or [`set_exchange_at!`](@ref).
+such as [`set_field_at!`](@ref) or [`set_exchange_at!`](@ref).
 
 Note that the definition of a cell may change when a system is reshaped. In this
 case, it is convenient to construct the `Site` using [`position_to_site`](@ref),
@@ -241,9 +242,8 @@ end
 """
     magnetic_moment(sys::System, site::Site)
 
-Returns ``μ/μ_B = - g 𝐒``, the local magnetic moment ``μ`` in units of the Bohr
-magneton ``μ_B``. The spin dipole ``𝐒`` and ``g``-tensor may both be
-[`Site`](@ref) dependent.
+Returns ``- g 𝐒``, the local magnetic moment in units of the Bohr magneton. The
+spin dipole ``𝐒`` and ``g``-tensor may both be [`Site`](@ref) dependent.
 """
 function magnetic_moment(sys::System, site)
     site = to_cartesian(site)
diff --git a/src/System/Types.jl b/src/System/Types.jl
index 2cbba9912..81fca5296 100644
--- a/src/System/Types.jl
+++ b/src/System/Types.jl
@@ -71,6 +71,7 @@ const rBFTPlan = FFTW.rFFTWPlan{ComplexF64, 1, false, 5, UnitRange{Int64}}
 const rIFTPlan = FFTW.AbstractFFTs.ScaledPlan{ComplexF64, rBFTPlan, Float64}
 
 struct Ewald
+    μ0_μB²   :: Float64               # Strength of dipole-dipole interactions
     A        :: Array{Mat3, 5}        # Interaction matrices in real-space         [offset+1,i,j]
     μ        :: Array{Vec3, 4}        # Magnetic moments μ = g s                   [cell,i]
     ϕ        :: Array{Vec3, 4}        # Cross correlation, ϕ = A⋆μ                 [cell,i]
@@ -110,6 +111,5 @@ mutable struct System{N}
     const coherent_buffers :: Vector{Array{CVec{N}, 4}} # Buffers for dynamics routines
 
     # Global data
-    const units            :: PhysicalConsts
     const rng              :: Random.Xoshiro
 end
diff --git a/src/Units.jl b/src/Units.jl
index 3372d74dd..e97364d77 100644
--- a/src/Units.jl
+++ b/src/Units.jl
@@ -1,35 +1,68 @@
-"""
-    meV_per_K = 0.086173332621451774
+# Angstrom in selected length unit
+const angstrom_in = Base.ImmutableDict(
+    :angstrom => 1.0, # == Å / Å
+    :nm => 0.1,       # == Å / nm
+)
 
-Boltzmann's constant ``k_B`` in units of meV per kelvin.
-"""
-const meV_per_K = 0.086173332621451774
+# meV in selected energy unit
+const meV_in = Base.ImmutableDict(
+    :meV => 1.0,                   # meV / meV
+    :K   => 801088317/69032450,    # meV / kB K
+    :THz => 267029439/1104345025,  # meV / h THz
+    :T   => 17.275985474052367519, # meV / μB T
+)
 
-Base.@kwdef struct PhysicalConsts
-    μ0::Float64    # Vacuum permeability
-    μB::Float64    # Bohr magneton 
-end
 
 """
-    Units.meV
-    Units.theory
+    Units(energy, length=:angstrom)
+
+Physical constants in units of reference `energy` and `length` scales. Possible
+lengths are `[:angstrom, :nm]` and possible energies are `[:meV, :K, :THz]`.
+Kelvin is converted to energy via the Boltzmann constant ``k_B``. Similarly,
+hertz is converted via the Planck constant ``h``, and tesla (field strength) via
+the Bohr magneton ``μ_B``. For a given `Units` system, one can access length
+scales (`angstrom`, `nm`) and energy scales (`meV`, `K`, `THz`, `T`).
+
+# Examples
+
+```julia
+# Create a unit system where energies are measured in meV
+units = Units(:meV)
 
-The default units system is `Units.meV`, which employs (meV, Å, tesla). Time is
-measured as an inverse energy, where factors of ``ħ`` are implicit.
+# Use the Boltzmann constant ``k_B`` to convert 1 kelvin into meV
+@assert units.K ≈ 0.0861733326
 
-An arbitrary unit system can be selected via two physical constants: the Bohr
-magneton ``μ_B`` and the vacuum permeability ``μ₀``. The choice `Units.theory`
-selects ``μ_B = 1``, such that the external magnetic field has energy units.
+# Use the Planck constant ``h`` to convert 1 THz into meV
+@assert units.THz ≈ 4.135667696
 
-See also [`meV_per_K`](@ref) to convert between temperature and energy.
+# Use the Bohr magneton ``μ_B`` to convert 1 tesla into meV
+@assert units.T ≈ 0.05788381806
+
+# The physical constant ``μ_0 μ_B²`` in units of ų meV.
+@assert u.vacuum_permeability ≈ 0.6745817653
+```
 """
-const Units = (;
-    meV = PhysicalConsts(;
-        μ0 = 201.33545383470705041,   # T^2 Å^3 / meV
-        μB = 0.057883818060738013331, # meV / T
-    ),
-    theory = PhysicalConsts(;
-        μ0 = NaN, # dipole-dipole interactions are invalid
-        μB = 1.0, # arbitrary energy units
-    ),
-)
+struct Units{E, L}
+    function Units(energy, length=:angstrom)
+        length in keys(angstrom_in) || error("`length` must be one of $(keys(angstrom_in))")
+        energy in keys(meV_in) || error("`energy` must be one of $(keys(meV_in))")
+        return new{energy, length}()
+    end
+end
+
+function Base.getproperty(u::Units{E, L}, name::Symbol) where {E, L}
+    if name in (:meV, :K, :THz, :T)
+        return meV_in[E] / meV_in[name]
+    end
+    
+    if name in (:angstrom, :nm)
+        return angstrom_in[L] / angstrom_in[name]
+    end
+
+    if name == :vacuum_permeability
+        # 0.6745... = μ0 μB² / ų meV
+        return 0.6745817653324668 * u.angstrom^3 * u.meV
+    end
+
+    error("type Units has no constant $name")
+end
diff --git a/src/deprecated.jl b/src/deprecated.jl
index a596134aa..495f818e2 100644
--- a/src/deprecated.jl
+++ b/src/deprecated.jl
@@ -34,26 +34,53 @@ function Base.setproperty!(value::Langevin, name::Symbol, x)
 end
 
 
-function lorentzian(x, η)
+Base.@deprecate lorentzian(x, η) let
     @warn "`lorentzian(x, η)` is deprecated! Use `lorentzian(; fwhm=2η)(x)` instead."
     return lorentzian(; fwhm=2η)(x)
 end
 
-function lorentzian(η)
+Base.@deprecate lorentzian(η) let
     @warn "`lorentzian(η)` is deprecated! Use `lorentzian(; fwhm=2η)` instead."
     return lorentzian(; fwhm=2η)
 end
 
-function integrated_lorentzian(η::Float64)
+Base.@deprecate integrated_lorentzian(η::Float64) let
     @warn "`integrated_lorentzian(η)` is deprecated! Use `integrated_lorentzian(; fwhm=2η)` instead."
     return integrated_lorentzian(; fwhm=2η)
 end
 
+Base.@deprecate set_external_field!(sys::System, B) let
+    @warn "`set_external_field!(sys, B)` is deprecated! Consider `set_field!(sys, B*units.T)` where `units = Units(:meV)`."
+    set_field!(sys, Vec3(B) * Units(:meV).T)
+end
+
+Base.@deprecate set_external_field_at!(sys::System, B, site) let
+    @warn "`set_external_field_at!(sys, B, site)` is deprecated! Consider `set_field_at!(sys, B*units.T, site)` where `units = Units(:meV)`."
+    set_field_at!(sys, Vec3(B) * Units(:meV).T, site)
+end
+
+
+# Consider `units.K` where `units = Units(:meV)`.
+Base.@deprecate_binding meV_per_K Units(:meV).K
 
-# REMEMBER TO ALSO DELETE:
 
-# view_crystal(cryst, max_dist)
-# λ argument in Langevin constructor
-# Δt argument in dynamical_correlations
-# large_S argument in set_exchange! and set_exchange_at!
-# Argument `q` in set_spiral_order*
+function Base.getproperty(x::Type{Units}, name::Symbol)
+    if name in (:theory, :meV)
+        @warn "Units.$name is deprecated! See `Units` docs for new interface."
+        return nothing
+    end
+    return getfield(x, name)
+end
+
+
+# REMEMBER TO ALSO DELETE:
+#
+# * view_crystal(cryst, max_dist)
+# * λ argument in Langevin constructor
+# * Δt argument in dynamical_correlations
+# * large_S argument in set_exchange! and set_exchange_at!
+# * Argument q in set_spiral_order*
+# * Argument units to System
+# * Missing μ0_μB² in enable_dipole_dipole! and
+#   modify_exchange_with_truncated_dipole_dipole!
+# * hermitianpart[!] for VERSION < v"1.10"
diff --git a/test/Project.toml b/test/Project.toml
index 2f3482295..e1f951077 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -12,3 +12,6 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b"
+
+[compat]
+TestItemRunner = "0.2.3"
diff --git a/test/shared.jl b/test/shared.jl
index be9c46db8..e00001e3c 100644
--- a/test/shared.jl
+++ b/test/shared.jl
@@ -9,7 +9,7 @@ using Random, LinearAlgebra, IOCapture
 # Various possible interactions appropriate to diamond crystal
 
 function add_linear_interactions!(sys, mode)
-    set_external_field!(sys, (0.0, 1.0, 1.0))
+    set_field!(sys, [0.0, 0.2, 0.3])
     if mode == :SUN
         # Kets scale as z → √κ z, so ⟨Λ⟩ → κ ⟨Λ⟩ is linear in κ
         set_onsite_coupling!(sys, S -> 0.2*(S[1]^4+S[2]^4+S[3]^4), 1)
diff --git a/test/test_chirality.jl b/test/test_chirality.jl
index 5b548c20a..232b65c8d 100644
--- a/test/test_chirality.jl
+++ b/test/test_chirality.jl
@@ -1,41 +1,39 @@
 @testitem "Spin precession handedness" begin
     using LinearAlgebra
 
-    for units in (Units.meV, Units.theory)
-        crystal = Crystal(lattice_vectors(1, 1, 1, 90, 90, 90), [[0, 0, 0]])
-        sys_dip = System(crystal, (1, 1, 1), [SpinInfo(1; S=1, g=2)], :dipole; units)
-        sys_sun = System(crystal, (1, 1, 1), [SpinInfo(1; S=1, g=2)], :SUN; units)
-    
-        B = [0, 0, 1] / abs(units.μB)
-        set_external_field!(sys_dip, B)
-        set_external_field!(sys_sun, B)
-    
-        ic = [1/√2, 0, 1/√2]
-        set_dipole!(sys_dip, ic, (1, 1, 1, 1))
-        set_dipole!(sys_sun, ic, (1, 1, 1, 1))
-    
-        integrator = ImplicitMidpoint(0.05)
-        for _ in 1:5
-            step!(sys_dip, integrator)
-            step!(sys_sun, integrator)
-        end
-    
-        dip_is_lefthanded = B ⋅ (ic × magnetic_moment(sys_dip, (1,1,1,1))) < 0
-        sun_is_lefthanded = B ⋅ (ic × magnetic_moment(sys_sun, (1,1,1,1))) < 0
-    
-        @test dip_is_lefthanded == sun_is_lefthanded == true
+    crystal = Crystal(lattice_vectors(1, 1, 1, 90, 90, 90), [[0, 0, 0]])
+    sys_dip = System(crystal, (1, 1, 1), [SpinInfo(1; S=1, g=2)], :dipole)
+    sys_sun = System(crystal, (1, 1, 1), [SpinInfo(1; S=1, g=2)], :SUN)
+
+    B = [0, 0, 1]
+    set_field!(sys_dip, B)
+    set_field!(sys_sun, B)
+
+    ic = [1/√2, 0, 1/√2]
+    set_dipole!(sys_dip, ic, (1, 1, 1, 1))
+    set_dipole!(sys_sun, ic, (1, 1, 1, 1))
+
+    integrator = ImplicitMidpoint(0.05)
+    for _ in 1:5
+        step!(sys_dip, integrator)
+        step!(sys_sun, integrator)
     end
+
+    dip_is_lefthanded = B ⋅ (ic × magnetic_moment(sys_dip, (1,1,1,1))) < 0
+    sun_is_lefthanded = B ⋅ (ic × magnetic_moment(sys_sun, (1,1,1,1))) < 0
+
+    @test dip_is_lefthanded == sun_is_lefthanded == true
 end
 
 
 @testitem "DM chain" begin
     latvecs = lattice_vectors(2, 2, 1, 90, 90, 90)
     cryst = Crystal(latvecs, [[0,0,0]], "P1")
-    sys = System(cryst, (1,1,1), [SpinInfo(1,S=1,g=-1)], :dipole; units=Units.theory)
+    sys = System(cryst, (1,1,1), [SpinInfo(1,S=1,g=-1)], :dipole)
     D = 1
     B = 10.0
     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])
 
     # Above the saturation field, the ground state is fully polarized, with no
     # energy contribution from the DM term.
@@ -57,7 +55,7 @@ end
 
     sys2 = resize_supercell(sys, (1, 1, 4))
     B = 1
-    set_external_field!(sys2, [0, 0, B])
+    set_field!(sys2, [0, 0, B])
     randomize_spins!(sys2)
     minimize_energy!(sys2)
     @test energy_per_site(sys2) ≈ -5/4
@@ -93,7 +91,7 @@ end
     # Finally, test fully polarized state
 
     B = 10
-    set_external_field!(sys3, [0, 0, B])
+    set_field!(sys3, [0, 0, B])
     polarize_spins!(sys3, [0, 0, 1])
     @test energy_per_site(sys3) ≈ -B
     swt = SpinWaveTheory(sys3)
diff --git a/test/test_energy_consistency.jl b/test/test_energy_consistency.jl
index 8c0ab0a49..850b2be88 100644
--- a/test/test_energy_consistency.jl
+++ b/test/test_energy_consistency.jl
@@ -8,11 +8,11 @@
         add_linear_interactions!(sys, mode)
         add_quadratic_interactions!(sys, mode)
         add_quartic_interactions!(sys, mode)
-        enable_dipole_dipole!(sys)
+        enable_dipole_dipole!(sys, 0.5)
 
         # Random field
         for site in eachsite(sys)
-            set_external_field_at!(sys, randn(sys.rng, 3), site)
+            set_field_at!(sys, 0.1*randn(sys.rng, 3), site)
         end
         # Random spin rescaling
         rand!(sys.rng, sys.κs)
diff --git a/test/test_ewald.jl b/test/test_ewald.jl
index ff1cb0ff7..d048cc3e6 100644
--- a/test/test_ewald.jl
+++ b/test/test_ewald.jl
@@ -10,7 +10,7 @@
         # magnetic moments
         dipoles = [magnetic_moment(sys, site) for site in eachsite(sys)][:]
         # energy from traditional Ewald summation
-        Ewalder.energy(Ewalder.System(; latvecs, pos); dipoles) / (4π*sys.units.μ0)
+        Ewalder.energy(Ewalder.System(; latvecs, pos); dipoles) / 4π
     end
 
     # Long-range energy of single dipole in cubic box with PBC
@@ -18,15 +18,14 @@
     positions = [[0,0,0]]
     cryst = Crystal(latvecs, positions)
     infos = [SpinInfo(1, S=1, g=1)]
-    units = Sunny.PhysicalConsts(; μ0=1, μB=1)
-    sys = System(cryst, (1,1,1), infos, :dipole; units)
-    enable_dipole_dipole!(sys)
+    sys = System(cryst, (1,1,1), infos, :dipole)
+    enable_dipole_dipole!(sys, 1.0)
     @test ewalder_energy(sys) ≈ -1/6
     @test isapprox(energy(sys), -1/6; atol=1e-13)
 
     # Same thing, with multiple unit cells
-    sys = System(cryst, (2,3,4), infos, :dipole; units)
-    enable_dipole_dipole!(sys)
+    sys = System(cryst, (2,3,4), infos, :dipole)
+    enable_dipole_dipole!(sys, 1.0)
     @test isapprox(energy_per_site(sys), -1/6; atol=1e-13)
 
     # Create a random box
@@ -39,8 +38,8 @@
         SpinInfo(2, S=3/2, g=rand(3,3)),
         SpinInfo(3, S=2, g=rand(3,3)),
     ]
-    sys = System(cryst, (1,1,1), infos, :dipole; units)
-    enable_dipole_dipole!(sys)
+    sys = System(cryst, (1,1,1), infos, :dipole)
+    enable_dipole_dipole!(sys, 1.0)
     randomize_spins!(sys)
 
     # Energy per site is independent of resizing
diff --git a/test/test_lswt.jl b/test/test_lswt.jl
index 1ac3e8a73..030167ce5 100644
--- a/test/test_lswt.jl
+++ b/test/test_lswt.jl
@@ -141,10 +141,11 @@ end
     positions = [[0, 0, 0]]
     fcc = Crystal(latvecs, positions, 225)
 
+    units = Units(:meV)
     S = 5/2
     g = 2
-    J = 22.06 * meV_per_K
-    K = 0.15  * meV_per_K
+    J = 22.06 * units.K
+    K = 0.15  * units.K
     C = J + K
     J₁ = diagm([J, J, C])
     D = 25/24
@@ -249,10 +250,10 @@ end
         cryst = Crystal(latvecs, positions)
         q = [0.12, 0.23, 0.34]
         
-        sys = System(cryst, (1, 1, 1), [SpinInfo(1; S, g=-1)], :dipole; units=Units.theory)
+        sys = System(cryst, (1, 1, 1), [SpinInfo(1; S, g=-1)], :dipole)
         set_exchange!(sys, J, Bond(1, 1, [1, 0, 0]))
         set_onsite_coupling!(sys, S -> D*S[3]^2, 1)
-        set_external_field!(sys, [0, 0, h])
+        set_field!(sys, [0, 0, h])
 
         # Numerical
         sys_swt_dip = reshape_supercell(sys, [1 -1 0; 1 1 0; 0 0 1])
@@ -289,12 +290,12 @@ end
 
     dims = (1, 1, 1)
     S = 2
-    sys_dip = System(cryst, dims, [SpinInfo(1; S, g=-1)], :dipole; units=Units.theory)
-    sys_SUN = System(cryst, dims, [SpinInfo(1; S, g=-1)], :SUN; units=Units.theory)
+    sys_dip = System(cryst, dims, [SpinInfo(1; S, g=-1)], :dipole)
+    sys_SUN = System(cryst, dims, [SpinInfo(1; S, g=-1)], :SUN)
 
     # The strengths of single-ion anisotropy (must be negative to favor the dipolar ordering under consideration)
     Ds = -rand(3)
-    h  = rand()
+    h  = 0.1*rand()
     M = normalize(rand(3))
     SM = M' * spin_matrices(S)
     aniso = Ds[1]*SM^2 + Ds[2]*SM^4 + Ds[3]*SM^6
@@ -302,8 +303,8 @@ end
     set_onsite_coupling!(sys_dip, aniso, 1)
     set_onsite_coupling!(sys_SUN, aniso, 1)
 
-    set_external_field!(sys_dip, h*M)
-    set_external_field!(sys_SUN, h*M)
+    set_field!(sys_dip, h*M)
+    set_field!(sys_SUN, h*M)
 
     set_dipole!(sys_dip, M, (1,1,1,1))
     set_dipole!(sys_SUN, M, (1,1,1,1))
@@ -359,17 +360,17 @@ end
 
     for mode in (:dipole, :SUN)
         sys = System(cryst, (1,1,1), [SpinInfo(1; S=1, g=1)], mode)
-        enable_dipole_dipole!(sys)
+        enable_dipole_dipole!(sys, 1.0)
 
         polarize_spins!(sys, (0,0,1))
-        @test energy_per_site(sys) ≈ -0.1913132980155851 * (sys.units.μ0 * sys.units.μB^2)
+        @test energy_per_site(sys) ≈ -0.1913132980155851
         
         swt = SpinWaveTheory(sys)
         formula = intensity_formula(swt, :perp; kernel=delta_function_kernel)
         
         qpoints = [[0, 0, 0], [0, 0, 1/2], [0, 1/2, 1/2], [0, 0, 0]]
         disps, is = intensities_bands(swt, qpoints, formula)
-        disps_ref = [0.5689399140467553, 0.23914164251944922, 0.23914164251948083, 0.5689399140467553] * (sys.units.μ0 * sys.units.μB^2)
+        disps_ref = [0.5689399140467553, 0.23914164251944922, 0.23914164251948083, 0.5689399140467553]
         @test isapprox(disps[:,end], disps_ref; atol=1e-7)
         @test is[:,end] ≈ [1, 1, 201/202, 1]
     end
@@ -377,7 +378,7 @@ end
     begin
         cryst = Sunny.bcc_crystal()
         sys = System(cryst, (1, 1, 1), [SpinInfo(1, S=1, g=2)], :dipole, seed=2)
-        enable_dipole_dipole!(sys)
+        enable_dipole_dipole!(sys, Units(:meV).vacuum_permeability)
         polarize_spins!(sys, (1,2,3)) # arbitrary direction
         
         R = hcat([1,1,-1], [-1,1,1], [1,-1,1]) / 2
@@ -433,7 +434,7 @@ end
     tol = 1e-7
     latvecs = lattice_vectors(1, 1, 10, 90, 90, 120)
     cryst = Crystal(latvecs, [[0, 0, 0]])
-    sys = System(cryst, (7,7,1), [SpinInfo(1; S=1, g=-1)], :dipole, units=Units.theory, seed=0)
+    sys = System(cryst, (7,7,1), [SpinInfo(1; S=1, g=-1)], :dipole; seed=0)
 
     J₁ = -1
     J₃ = 1.6234898018587323
@@ -444,7 +445,7 @@ end
     set_onsite_coupling!(sys, S -> D*S[3]^2, 1)
 
     h = 2.0*J₃
-    set_external_field!(sys, [0, 0, h])
+    set_field!(sys, [0, 0, h])
     
     dipole_array = zeros(Float64,3,7,7)
     dipole_array[:,:,1]= [0.856769   0.811926  -0.485645   -0.950288   -0.429787  -0.0800935  0.217557;
@@ -538,7 +539,7 @@ end
         K2 = Sunny.tracelesspart(K2 + K2')
 
         set_pair_coupling!(sys, (S1, S2) -> S1'*R*J*R'*S2 + (S1'*R*K1*R'*S1)*(S2'*R*K2*R'*S2), Bond(1, 2, [0,0,0]))
-        set_external_field!(sys, R*h)
+        set_field!(sys, R*h)
 
         return sys
     end
@@ -629,7 +630,7 @@ end
         K2 = diagm([-1, -1, 2])
         set_pair_coupling!(sys, (Si, Sj) -> -Si'*Sj + (Si'*K1*Si)*(Sj'*K2*Sj), Bond(1,1,[1,0,0]); extract_parts=true)
         set_onsite_coupling!(sys, S -> S[3]^2, 1)
-        set_external_field!(sys, [0,0,0.1])
+        set_field!(sys, [0,0,0.1])
 
         randomize_spins!(sys)
         minimize_energy!(sys; maxiters=1_000)
diff --git a/test/test_optimization.jl b/test/test_optimization.jl
index 9c9b0b6c7..82bd04c38 100644
--- a/test/test_optimization.jl
+++ b/test/test_optimization.jl
@@ -73,7 +73,7 @@ end
         sys = System(cryst, dims, [SpinInfo(1; S, g=2)], mode; seed) 
         set_exchange!(sys, -1, Bond(1,1,[1,0,0]))
         set_onsite_coupling!(sys, S -> -S[3]^2, 1)
-        enable_dipole_dipole!(sys)
+        enable_dipole_dipole!(sys, 1.0)
         sys
     end
 
diff --git a/test/test_spiral.jl b/test/test_spiral.jl
index 34932fb12..b1f0c283c 100644
--- a/test/test_spiral.jl
+++ b/test/test_spiral.jl
@@ -17,7 +17,7 @@
     # compute ewald energy using J(k)
     sys = System(cryst, (1, 1, 1), [SpinInfo(1, S=1, g=1)], :dipole, seed=0)
     randomize_spins!(sys)
-    enable_dipole_dipole!(sys)
+    enable_dipole_dipole!(sys, 1.0)
     E1 = Sunny.spiral_energy_per_site(sys; k, axis)
 
     # compute ewald energy using supercell
@@ -25,7 +25,7 @@
     for i in 1:Sunny.natoms(sys.crystal)
         set_spiral_order_on_sublattice!(sys_large, i; k, axis, S0=sys.dipoles[i])
     end
-    enable_dipole_dipole!(sys_large)
+    enable_dipole_dipole!(sys_large, 1.0)
     E2 = energy_per_site(sys_large)
 
     @test E1 ≈ E2
@@ -40,7 +40,7 @@
     # compute ewald energy using J(Q)
     sys = System(cryst, (1, 1, 1), [SpinInfo(1, S=1, g=1)], :dipole, seed=0)
     randomize_spins!(sys)
-    enable_dipole_dipole!(sys)
+    enable_dipole_dipole!(sys, 1.0)
     E1 = Sunny.spiral_energy_per_site(sys; k, axis)
 
     # compute ewald energy using supercell
@@ -48,7 +48,7 @@
     for i in 1:Sunny.natoms(sys.crystal)
         set_spiral_order_on_sublattice!(sys_large, i; k, axis, S0=sys.dipoles[i])
     end
-    enable_dipole_dipole!(sys_large)
+    enable_dipole_dipole!(sys_large, 1.0)
     E2 = energy_per_site(sys_large)
 
     @test E1 ≈ E2
@@ -65,10 +65,10 @@ end
         cryst = Crystal(latvecs, positions)
 
         dims = (1, 1, 1)
-        sys = System(cryst, dims, [SpinInfo(1; S, g=-1)], :dipole; units=Units.theory)
+        sys = System(cryst, dims, [SpinInfo(1; S, g=-1)], :dipole)
         set_exchange!(sys, J, Bond(1, 1, [1, 0, 0]))
         set_onsite_coupling!(sys, S -> (D/rcs)*S[3]^2, 1)
-        set_external_field!(sys, [0, 0, h])
+        set_field!(sys, [0, 0, h])
 
         k = Sunny.minimize_luttinger_tisza_exchange(sys; k_guess=randn(3))
         @test k[1:2] ≈ [0.5, 0.5]
diff --git a/test/test_symmetry.jl b/test/test_symmetry.jl
index c9a0d17cb..68ed5d890 100644
--- a/test/test_symmetry.jl
+++ b/test/test_symmetry.jl
@@ -140,18 +140,20 @@ end
 
 @testitem "Standardize" begin
     using LinearAlgebra
-
-    function test_standardize(latvecs, r)
-        cryst = Crystal(latvecs, [r])
+    
+    function test_standardize(cryst)
         cryst2 = standardize(cryst; idealize=false)
-        @test cryst2.latvecs * cryst2.positions[1] ≈ cryst.latvecs * r
+        @test cryst2.latvecs * cryst2.positions[1] ≈ cryst.latvecs * cryst.positions[1]
         cryst3 = standardize(cryst)
         @test norm(cryst3.positions[1]) < 1e-12    
     end
     
-    r = [0.1, 0.2, 0.3]
-    test_standardize([1 0 1; 1 1 0; 0 1 1], r)
-    test_standardize(lattice_vectors(1, 1, 1, 90, 90, 60), r)
+    cryst = Crystal([1 0 1; 1 1 0; 0 1 1], [[0.1, 0.2, 0.3]])
+    test_standardize(cryst)
+    
+    msg = "Found a nonconventional hexagonal unit cell. Consider using `lattice_vectors(a, a, c, 90, 90, 120)`."
+    @test_warn msg cryst = Crystal(lattice_vectors(1, 1, 1, 90, 90, 60), [[0.1, 0.2, 0.3]])
+    test_standardize(cryst)
 end
 
 
diff --git a/test/test_units.jl b/test/test_units.jl
new file mode 100644
index 000000000..68ea9b683
--- /dev/null
+++ b/test/test_units.jl
@@ -0,0 +1,19 @@
+@testitem "Units" begin
+    function f(units)
+        crystal = Sunny.diamond_crystal(; a=units.angstrom)
+        sys = System(crystal, (4, 4, 4), [SpinInfo(1, S=1, g=2)], :dipole; seed=0)
+        randomize_spins!(sys)
+        set_exchange!(sys, 2 * units.K, Bond(1, 2, [0,0,0]))
+        set_field!(sys, [0, 0, 1] * units.T)
+        enable_dipole_dipole!(sys, units.vacuum_permeability)
+        E = energy(sys) / units.meV
+        ∇E = Sunny.energy_grad_dipoles(sys) / units.THz
+        return (E, ∇E)
+    end
+
+    (E1, ∇E1) = f(Units(:meV, :angstrom))
+    (E2, ∇E2) = f(Units(:THz, :nm))
+
+    @test E1 ≈ E2
+    @test ∇E1 ≈ ∇E2
+end
diff --git a/test/test_units_scaling.jl b/test/test_units_scaling.jl
deleted file mode 100644
index d3f1e089f..000000000
--- a/test/test_units_scaling.jl
+++ /dev/null
@@ -1,65 +0,0 @@
-@testitem "Dimensional Scaling" begin
-    include("shared.jl")
-
-    # These tests checks that varying μ0, μB results in the correct scale
-    # transformations in various energies / fields. Note:
-    # does not check that the actual absolute values are
-    # correct.
-
-    crystal = Sunny.diamond_crystal()
-    latsize = (4, 4, 4)
-    infos = [SpinInfo(1, S=1, g=2)]
-
-    units = [
-        Sunny.PhysicalConsts(; μ0=2, μB=3)
-        Sunny.PhysicalConsts(; μ0=1, μB=1)
-    ]
-
-    function collect_energy_and_grad(test, units)
-        sys = System(crystal, latsize, infos, :dipole; units, seed=0)
-        if test == :zeeman
-            set_external_field!(sys, randn(sys.rng, Sunny.Vec3))
-        elseif test == :exchange
-            add_exchange_interactions!(sys, false)
-        else test == :dipdip
-            enable_dipole_dipole!(sys)
-        end
-        randomize_spins!(sys)
-        return energy(sys), Sunny.energy_grad_dipoles(sys)
-    end
-
-    # All exchange interactions should be invariant under changes to physical
-    # constants
-    function validate_exchanges_scaling()
-        E1, ∇E1 = collect_energy_and_grad(:exchange, units[1])
-        E2, ∇E2 = collect_energy_and_grad(:exchange, units[2])
-        @test E1 ≈ E2
-        @test ∇E1 ≈ ∇E2
-    end
-    validate_exchanges_scaling()
-
-    # Zeeman energies/forces should scale linearly with μB, invariant to μ0
-    function validate_zeeman_scaling()
-        E1, ∇E1 = collect_energy_and_grad(:zeeman, units[1])
-        E2, ∇E2 = collect_energy_and_grad(:zeeman, units[2])
-
-        @test E1 / units[1].μB ≈ E2 / units[2].μB
-        @test ∇E1 / units[1].μB ≈ ∇E2 / units[2].μB
-    end
-    validate_zeeman_scaling()
-
-    # Dipole-dipole interactions should scale linearly with μ0, and
-    # quadratically with μB
-    function validate_dipole_scaling()
-        E1, ∇E1 = collect_energy_and_grad(:dipdip, units[1])
-        E2, ∇E2 = collect_energy_and_grad(:dipdip, units[2])
-
-        (; μB, μ0) = units[1]
-        c1 = μ0*μB^2
-        (; μB, μ0) = units[2]
-        c2 = μ0*μB^2
-        @test E1 / c1 ≈ E2 / c2
-        @test ∇E1 / c1 ≈ ∇E2 / c2
-    end
-    validate_dipole_scaling()
-end