From da1444ec0331fc2d093593250011f644176c57a1 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Thu, 4 Jul 2024 08:34:35 -0600 Subject: [PATCH] Change magnetic moment conventions (#284) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * `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 | 5 +++++ examples/05_MC_Ising.jl | 10 ++++----- examples/06_CP2_Skyrmions.jl | 8 +++---- examples/08_Momentum_Conventions.jl | 4 ++-- examples/extra/Advanced_MC/PT_WHAM_ising2d.jl | 2 +- examples/extra/Advanced_MC/REWL_ising2d.jl | 2 +- examples/extra/Advanced_MC/WL_ising2d.jl | 2 +- ext/PlottingExt.jl | 2 +- src/MCIF.jl | 2 +- src/Spiral/SpiralEnergy.jl | 2 +- src/System/Ewald.jl | 2 +- src/System/Interactions.jl | 3 ++- src/System/System.jl | 17 ++++---------- src/Units.jl | 22 +++++++++---------- test/test_chirality.jl | 2 +- test/test_ewald.jl | 7 +++--- test/test_lswt.jl | 17 +++++++------- test/test_spiral.jl | 2 +- test/test_units_scaling.jl | 11 ++++------ 19 files changed, 58 insertions(+), 64 deletions(-) diff --git a/docs/src/versions.md b/docs/src/versions.md index efd9939ab..0a47a0874 100644 --- a/docs/src/versions.md +++ b/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) diff --git a/examples/05_MC_Ising.jl b/examples/05_MC_Ising.jl index 0dac4ecce..900ad2682 100644 --- a/examples/05_MC_Ising.jl +++ b/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 diff --git a/examples/06_CP2_Skyrmions.jl b/examples/06_CP2_Skyrmions.jl index 7135d90c3..9ba49e00a 100644 --- a/examples/06_CP2_Skyrmions.jl +++ b/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 diff --git a/examples/08_Momentum_Conventions.jl b/examples/08_Momentum_Conventions.jl index c0e8342a1..e50746fad 100644 --- a/examples/08_Momentum_Conventions.jl +++ b/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 diff --git a/examples/extra/Advanced_MC/PT_WHAM_ising2d.jl b/examples/extra/Advanced_MC/PT_WHAM_ising2d.jl index 93dbc32a3..e4b3c04b6 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, units=Units.theory, 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 95cdbdf9f..1ecc31af9 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, units=Units.theory, 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 67bc5cd3f..aefdda4cb 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, units=Units.theory, seed=0) polarize_spins!(sys, (0, 0, 1)) set_exchange!(sys, -1.0, Bond(1,1,(1,0,0))) diff --git a/ext/PlottingExt.jl b/ext/PlottingExt.jl index 9dd968793..289f0a2ba 100644 --- a/ext/PlottingExt.jl +++ b/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 diff --git a/src/MCIF.jl b/src/MCIF.jl index f96fd0c52..6c6284398 100644 --- a/src/MCIF.jl +++ b/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") diff --git a/src/Spiral/SpiralEnergy.jl b/src/Spiral/SpiralEnergy.jl index c1c08c452..2823797d0 100644 --- a/src/Spiral/SpiralEnergy.jl +++ b/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) diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl index 2f089a852..1b738742e 100644 --- a/src/System/Ewald.jl +++ b/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 diff --git a/src/System/Interactions.jl b/src/System/Interactions.jl index a7410dbd6..ec6f937d4 100644 --- a/src/System/Interactions.jl +++ b/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 diff --git a/src/System/System.jl b/src/System/System.jl index 6da302942..ec68fad88 100644 --- a/src/System/System.jl +++ b/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 diff --git a/src/Units.jl b/src/Units.jl index 6627cae17..3372d74dd 100644 --- a/src/Units.jl +++ b/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,15 +14,14 @@ 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(; @@ -31,7 +29,7 @@ const Units = (; μ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 ), ) diff --git a/test/test_chirality.jl b/test/test_chirality.jl index 7ef84e57c..5b548c20a 100644 --- a/test/test_chirality.jl +++ b/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])) diff --git a/test/test_ewald.jl b/test/test_ewald.jl index c423523fb..ff1cb0ff7 100644 --- a/test/test_ewald.jl +++ b/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) diff --git a/test/test_lswt.jl b/test/test_lswt.jl index 58693c0d6..1ac3e8a73 100644 --- a/test/test_lswt.jl +++ b/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))) diff --git a/test/test_spiral.jl b/test/test_spiral.jl index 562935c17..34932fb12 100644 --- a/test/test_spiral.jl +++ b/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]) diff --git a/test/test_units_scaling.jl b/test/test_units_scaling.jl index 4c269a10b..d3f1e089f 100644 --- a/test/test_units_scaling.jl +++ b/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