Change magnetic moment conventions (#284)

* `magnetic_moment` return value in units of Bohr magneton μB.
* Must now use `g=-1` to favor spin dipole aligned with field, even in "theory" units.

docs/src/versions.md

@@ -3,6 +3,11 @@
 ## v0.6.1
 (In progress)
 
+* **Breaking changes**: [`magnetic_moment`](@ref) is now reported in units of
+  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)).
+
 ## v0.6.0
 (June 18, 2024)
 

examples/05_MC_Ising.jl

@@ -18,12 +18,12 @@ crystal = Crystal(latvecs, [[0,0,0]])
 # [`polarize_spins!`](@ref). Following the Ising convention, we will restrict
 # these spins to the $z$-axis and give them magnitude $S=1$.
 # 
-# By default, Sunny uses physical units, e.g. magnetic field in tesla. Here we
-# specify an alternative [`Units`](@ref) system, so that the Zeeman coupling
-# between the spin dipole $𝐬$ and an external field $𝐁$ has the dimensionless
-# form $-𝐁⋅𝐬$.
+# By default, Sunny expects the magnetic field in tesla. Selecting
+# [`Units.theory`](@ref Units) instead allows for dimensionless units. Following
+# 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, units=Units.theory, seed=0)
 polarize_spins!(sys, (0,0,1))
 
 # Use [`set_exchange!`](@ref) to include a ferromagnetic Heisenberg interaction

examples/06_CP2_Skyrmions.jl

@@ -29,13 +29,13 @@ lat_vecs = lattice_vectors(1, 1, 10, 90, 90, 120)
 basis_vecs = [[0,0,0]]
 cryst = Crystal(lat_vecs, basis_vecs)
 
-# The crystal is then used to create a spin [`System`](@ref). All parameters in
-# this model system are dimensionless, so we select "theory" units and set the
-# g-factor to one. 
+# Create a spin [`System`](@ref) containing $L×L$ cells. Selecting
+# [`Units.theory`](@ref Units) with $g=-1$ provides a dimensionless Zeeman
+# coupling of the form $-𝐁⋅𝐬$.
 
 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, units=Units.theory)
 
 # 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

examples/08_Momentum_Conventions.jl

@@ -22,8 +22,8 @@ cryst = Crystal(latvecs, [[0,0,0]], "P1")
 
 # Construct a 1D chain system that extends along ``𝐚₃``. The Hamiltonian
 # includes DM and Zeeman coupling terms, ``ℋ = ∑_j D ẑ ⋅ (𝐒_j × 𝐒_{j+1}) -
-# ∑_j 𝐁 ⋅ 𝐌_j``, where ``𝐌_j = - μ_B g 𝐒_j`` is the
-# [`magnetic_moment`](@ref) and ``𝐁 ∝ ẑ``.
+# ∑_j 𝐁 ⋅ μ_j``, where ``μ_j = - μ_B g 𝐒_j`` is the [`magnetic_moment`](@ref)
+# and ``𝐁 ∝ ẑ``.
 
 sys = System(cryst, (1, 1, 25), [SpinInfo(1; S=1, g=2)], :dipole; seed=0)
 D = 0.1 # meV

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, units=Units.theory, 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, units=Units.theory, 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, units=Units.theory, seed=0)
 polarize_spins!(sys, (0, 0, 1))
 set_exchange!(sys, -1.0, Bond(1,1,(1,0,0)))
 

ext/PlottingExt.jl

@@ -615,7 +615,7 @@ function draw_atoms_or_dipoles(; ax, full_crystal_toggle, dipole_menu, cryst, sy
                 N0 = norm(sys.Ns) / sqrt(L)
                 S0 = (N0 - 1) / 2
                 spin_dipoles = sys.dipoles[sites] / S0
-                magn_dipoles = magnetic_moment.(Ref(sys), sites) / (S0*g0*sys.units.μB)
+                magn_dipoles = magnetic_moment.(Ref(sys), sites) / (S0*g0)
                 for (dipoles, obs) in [(spin_dipoles, show_spin_dipoles), (magn_dipoles, show_magn_dipoles)]
                     a0 = 5ionradius
                     arrowsize = 0.4a0

src/MCIF.jl

@@ -114,7 +114,7 @@ function set_dipoles_from_mcif_aux!(sys; positions, moments, magn_operations, ma
             end
 
             # Get spin dipole by inverting the `magnetic_moment` transformation
-            dipole = inv(- sys.units.μB * sys.gs[site]) * μ_new
+            dipole = - sys.gs[site] \ μ_new
             s_prev = sys.dipoles[site]
             set_dipole!(sys, dipole, site)
             s_prev == zero(Vec3) || s_prev ≈ sys.dipoles[site] || error("Conflicting dipoles at site $site")

src/Spiral/SpiralEnergy.jl

@@ -80,7 +80,7 @@ function spiral_energy_and_gradient_aux!(dEds, sys::System{0}; k, axis)
 
     # See "spiral_energy.lyx" for derivation
     if !isnothing(sys.ewald)
-        μ = [magnetic_moment(sys, site) for site in eachsite(sys)]
+        μ = [sys.units.μB*magnetic_moment(sys, site) for site in eachsite(sys)]
 
         A0 = sys.ewald.A
         A0 = reshape(A0, Na, Na)

src/System/Ewald.jl

@@ -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] = magnetic_moment(sys, site)
+        μ[site] = sys.units.μB*magnetic_moment(sys, site)
     end
 
     E = 0.0

src/System/Interactions.jl

@@ -76,9 +76,10 @@ 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 system of [`Units`](@ref).
+determined by the system of [`Units`](@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)
     return
 end

src/System/System.jl

@@ -241,22 +241,13 @@ end
 """
     magnetic_moment(sys::System, site::Site)
 
-Get the magnetic moment for a [`Site`](@ref). This is the spin dipole ``𝐒``
-multiplied by the local ``g``-tensor and the Bohr magneton ``μ_B``.
-
-```math
-μ = - μ_B g 𝐒.
-```
-
-The constant ``μ_B`` depends on the choice of [`Units`](@ref). By default,
-energy is in meV and external field is in tesla, such that ``μ_B = 0.05788…``
-(meV/T). If `Units.theory` is selected instead, the constant ``μ_B = -1``
-effectively aligns magnetic and spin dipoles, which is a common convention when
-specifying theoretical model systems.
+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.
 """
 function magnetic_moment(sys::System, site)
     site = to_cartesian(site)
-    return - sys.units.μB * sys.gs[site] * sys.dipoles[site]
+    return - sys.gs[site] * sys.dipoles[site]
 end
 
 # Total volume of system

src/Units.jl

@@ -1,8 +1,7 @@
 """
     meV_per_K = 0.086173332621451774
 
-A physical constant. Useful for converting kelvin into the default energy units,
-meV.
+Boltzmann's constant ``k_B`` in units of meV per kelvin.
 """
 const meV_per_K = 0.086173332621451774
 
@@ -15,23 +14,22 @@ end
     Units.meV
     Units.theory
 
-The unit system is implicitly determined by the definition of two physical
-constants: the vacuum permeability ``μ₀`` and the Bohr magneton ``μ_B``.
-Temperatures are effectively measured in units of energy (``k_B = 1``) and time
-is effectively measured in units of inverse energy (``ħ = 1``). The default unit
-system, `Units.meV`, employs (meV, Å, tesla). Select alternatively
-`Units.theory` for a units system defined so that ``μ₀ = 1`` and ``μ_B = -1``,
-which produces a Zeeman coupling of ``-g 𝐁⋅𝐒``.
+The default units system is `Units.meV`, which employs (meV, Å, tesla). Time is
+measured as an inverse energy, where factors of ``ħ`` are implicit.
 
-See also [`meV_per_K`](@ref).
+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.
+
+See also [`meV_per_K`](@ref) to convert between temperature and energy.
 """
 const Units = (;
     meV = PhysicalConsts(;
         μ0 = 201.33545383470705041,   # T^2 Å^3 / meV
         μB = 0.057883818060738013331, # meV / T
     ),
     theory = PhysicalConsts(;
-        μ0 = 1.0,
-        μB = -1.0,
+        μ0 = NaN, # dipole-dipole interactions are invalid
+        μB = 1.0, # arbitrary energy units
     ),
 )

test/test_chirality.jl

@@ -31,7 +31,7 @@ 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; units=Units.theory)
     D = 1
     B = 10.0
     set_exchange!(sys, dmvec([0, 0, D]), Bond(1, 1, [0, 0, 1]))

test/test_ewald.jl

@@ -18,13 +18,14 @@
     positions = [[0,0,0]]
     cryst = Crystal(latvecs, positions)
     infos = [SpinInfo(1, S=1, g=1)]
-    sys = System(cryst, (1,1,1), infos, :dipole; units=Units.theory)
+    units = Sunny.PhysicalConsts(; μ0=1, μB=1)
+    sys = System(cryst, (1,1,1), infos, :dipole; units)
     enable_dipole_dipole!(sys)
     @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=Units.theory)
+    sys = System(cryst, (2,3,4), infos, :dipole; units)
     enable_dipole_dipole!(sys)
     @test isapprox(energy_per_site(sys), -1/6; atol=1e-13)
 
@@ -38,7 +39,7 @@
         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=Units.theory)
+    sys = System(cryst, (1,1,1), infos, :dipole; units)
     enable_dipole_dipole!(sys)
     randomize_spins!(sys)
 

test/test_lswt.jl

@@ -249,7 +249,7 @@ 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; units=Units.theory)
         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])
@@ -289,8 +289,8 @@ 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; units=Units.theory)
+    sys_SUN = System(cryst, dims, [SpinInfo(1; S, g=-1)], :SUN; units=Units.theory)
 
     # The strengths of single-ion anisotropy (must be negative to favor the dipolar ordering under consideration)
     Ds = -rand(3)
@@ -358,18 +358,19 @@ end
     cryst = Crystal(latvecs, [[0,0,0]])
 
     for mode in (:dipole, :SUN)
-        sys = System(cryst, (1,1,1), [SpinInfo(1; S=1, g=1)], mode; units=Units.theory)
+        sys = System(cryst, (1,1,1), [SpinInfo(1; S=1, g=1)], mode)
         enable_dipole_dipole!(sys)
 
         polarize_spins!(sys, (0,0,1))
-        @test energy_per_site(sys) ≈ -0.1913132980155851
+        @test energy_per_site(sys) ≈ -0.1913132980155851 * (sys.units.μ0 * sys.units.μB^2)
 
         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)
-        @test disps[:,end] ≈ [0.5689399140467553, 0.23914164251944922, 0.23914164251948083, 0.5689399140467553]
+        disps_ref = [0.5689399140467553, 0.23914164251944922, 0.23914164251948083, 0.5689399140467553] * (sys.units.μ0 * sys.units.μB^2)
+        @test isapprox(disps[:,end], disps_ref; atol=1e-7)
         @test is[:,end] ≈ [1, 1, 201/202, 1]
     end
 
@@ -432,7 +433,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, units=Units.theory, seed=0)
 
     J₁ = -1
     J₃ = 1.6234898018587323
@@ -725,7 +726,7 @@ end
         latvecs = lattice_vectors(a, a, a, 90, 90, 90)
         positions = [[0, 0, 0]]
         fcc = Crystal(latvecs, positions, 225)
-        sys_afm1 = System(fcc, (1, 1, 1), [SpinInfo(1; S, g=1)], mode; units=Units.theory)
+        sys_afm1 = System(fcc, (1, 1, 1), [SpinInfo(1; S, g=1)], mode)
         set_exchange!(sys_afm1, J, Bond(1, 2, [0, 0, 0]))
         set_dipole!(sys_afm1, (0, 0,  1), position_to_site(sys_afm1, (0, 0, 0)))
         set_dipole!(sys_afm1, (0, 0, -1), position_to_site(sys_afm1, (1/2, 1/2, 0)))

test/test_spiral.jl

@@ -65,7 +65,7 @@ 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; units=Units.theory)
         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])

test/test_units_scaling.jl

@@ -10,7 +10,10 @@
     latsize = (4, 4, 4)
     infos = [SpinInfo(1, S=1, g=2)]
 
-    units = [Sunny.Units.meV, Sunny.Units.theory]
+    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)
@@ -33,10 +36,8 @@
         @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])
@@ -45,10 +46,8 @@
         @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()
@@ -62,7 +61,5 @@
         @test E1 / c1 ≈ E2 / c2
         @test ∇E1 / c1 ≈ ∇E2 / c2
     end
-
     validate_dipole_scaling()
-
 end