diff --git a/examples/spinw_tutorials/SW05_Simple_kagome_FM.jl b/examples/spinw_tutorials/SW05_Simple_kagome_FM.jl index 4355beae2..93a9a016a 100644 --- a/examples/spinw_tutorials/SW05_Simple_kagome_FM.jl +++ b/examples/spinw_tutorials/SW05_Simple_kagome_FM.jl @@ -69,6 +69,6 @@ res2 = powder_average(cryst, radii, 1000) do qs end fig = Figure(size=(768, 800)) -plot_intensities!(fig[1, 1], res1; units, colorrange=(0,10)) -plot_intensities!(fig[2, 1], res2; units, colorrange=(0,10)) +plot_intensities!(fig[1, 1], res1; units, colorrange=(0,10), title="FWHM 0.02 meV") +plot_intensities!(fig[2, 1], res2; units, colorrange=(0,10), title="FWHM 0.25 meV") fig diff --git a/examples/spinw_tutorials/SW18_Distorted_kagome.jl b/examples/spinw_tutorials/SW18_Distorted_kagome.jl index 2a2f05690..c849c8b57 100644 --- a/examples/spinw_tutorials/SW18_Distorted_kagome.jl +++ b/examples/spinw_tutorials/SW18_Distorted_kagome.jl @@ -24,8 +24,8 @@ view_crystal(cryst) # Define the interactions. -spininfos = [1 => Moment(s=1/2, g=2), 3 => Moment(s=1/2, g=2)] -sys = System(cryst, spininfos, :dipole, seed=0) +moments = [1 => Moment(s=1/2, g=2), 3 => Moment(s=1/2, g=2)] +sys = System(cryst, moments, :dipole, seed=0) J = -2 Jp = -1 Jab = 0.75 diff --git a/ext/PlottingExt.jl b/ext/PlottingExt.jl index 4006b026a..0d4f30715 100644 --- a/ext/PlottingExt.jl +++ b/ext/PlottingExt.jl @@ -609,9 +609,8 @@ function draw_atoms_or_dipoles(; ax, full_crystal_toggle, dipole_menu, cryst, sy # Show dipoles. Mostly consistent with code in plot_spins. if !isnothing(sys) && ismagnetic sites = Sunny.position_to_site.(Ref(sys), rs) - L = length(eachsite(sys)) - g0 = norm(sys.gs) / sqrt(L * 3) - N0 = norm(sys.Ns) / sqrt(L) + g0 = norm(sys.gs) / sqrt(length(sys.gs) * 3) + N0 = norm(sys.Ns) / sqrt(length(sys.Ns)) s0 = (N0 - 1) / 2 spin_dipoles = sys.dipoles[sites] / s0 magn_dipoles = magnetic_moment.(Ref(sys), sites) / (s0*g0) diff --git a/src/Integrators.jl b/src/Integrators.jl index f23f9cafe..fe4710ddb 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -210,7 +210,7 @@ function suggest_timestep_aux(sys::System{N}, integrator; tol) where N # information about the frequency of oscillation. Consider, e.g., a spin # approximately aligned with an external field: the precession frequency is # given by |∇E| = |B|. - drift_rms = sqrt(acc/length(eachsite(sys))) + drift_rms = sqrt(acc/nsites(sys)) if iszero(drift_rms) error("Cannot suggest a timestep without an energy scale!") end diff --git a/src/KPM/SpinWaveTheoryKPM.jl b/src/KPM/SpinWaveTheoryKPM.jl index e7a61b462..421683d22 100644 --- a/src/KPM/SpinWaveTheoryKPM.jl +++ b/src/KPM/SpinWaveTheoryKPM.jl @@ -70,7 +70,7 @@ function intensities!(data, swt_kpm::SpinWaveTheoryKPM, qpts; energies, kernel:: @assert eltype(data) == eltype(measure) @assert size(data) == (length(energies), length(qpts.qs)) - Na = length(eachsite(sys)) + Na = nsites(sys) Ncells = Na / natoms(cryst) Nf = nflavors(swt) L = Nf*Na diff --git a/src/Measurements/IntensitiesTypes.jl b/src/Measurements/IntensitiesTypes.jl index 5c33076b3..75340b1a3 100644 --- a/src/Measurements/IntensitiesTypes.jl +++ b/src/Measurements/IntensitiesTypes.jl @@ -43,7 +43,12 @@ struct PowderIntensities{T} <: AbstractIntensities end function Base.show(io::IO, res::AbstractIntensities) - sz = join([size(res.data, 1), sizestr(res.qpts)], "×") + sz = string(size(res.data, 1)) * "×" * sizestr(res.qpts) + print(io, string(typeof(res)) * " ($sz elements)") +end + +function Base.show(io::IO, res::StaticIntensities) + sz = sizestr(res.qpts) print(io, string(typeof(res)) * " ($sz elements)") end diff --git a/src/MonteCarlo/WangLandau.jl b/src/MonteCarlo/WangLandau.jl index f3b8f5cfe..99857e8a9 100644 --- a/src/MonteCarlo/WangLandau.jl +++ b/src/MonteCarlo/WangLandau.jl @@ -31,8 +31,7 @@ function drive_system_to_energy_bounds!(sys::System{N}, bounds, propose::PR, max sampler = LocalSampler(; kT, nsweeps=0.01, propose) for _ in 1:max_sweeps step!(sys, sampler) - nsites = length(eachsite(sys)) - E = E0 + sampler.ΔE/nsites + E = E0 + sampler.ΔE/nsites(sys) inbounds(E, bounds) && return end @@ -64,17 +63,17 @@ end # sample for specified number of sweeps function step_ensemble!(WL::WangLandau, nsweeps::Int64) (; sys) = WL - nsites = length(eachsite(sys)) - WL.E = energy(sys) / nsites + N = nsites(sys) + WL.E = energy_per_site(sys) @assert inbounds(WL.E, WL.bounds) # try to fulfill the histogram criterion in set number of MC sweeps accepts = 0 - for _ in 1:nsweeps*nsites + for _ in 1:nsweeps*N # perform single-spin update and calculate proposal energy site = rand(sys.rng, eachsite(sys)) state = WL.propose(sys, site) - E_next = WL.E + local_energy_change(sys, site, state) / nsites + E_next = WL.E + local_energy_change(sys, site, state) / N if inbounds(E_next, WL.bounds) iszero(WL.ln_g[E_next]) && add_new!(WL, E_next) diff --git a/src/Operators/Rotation.jl b/src/Operators/Rotation.jl index 8e05badbe..cf89f976d 100644 --- a/src/Operators/Rotation.jl +++ b/src/Operators/Rotation.jl @@ -118,7 +118,8 @@ function unitary_for_rotation(R::Mat3, gen) return exp(-im*θ*(n'*gen)) end -# Unitary for a rotation matrix in the N-dimensional irrep of SU(2). +# Unitary for a rotation matrix in the N-dimensional irrep of SU(2). TODO: Use +# polynomial expansions for speed: http://www.emis.de/journals/SIGMA/2014/084/. function unitary_irrep_for_rotation(R::Mat3; N::Int) gen = spin_matrices_of_dim(; N) unitary_for_rotation(R, gen) diff --git a/src/SampledCorrelations/CorrelationUtils.jl b/src/SampledCorrelations/CorrelationUtils.jl index 5c4f791e0..845d19482 100644 --- a/src/SampledCorrelations/CorrelationUtils.jl +++ b/src/SampledCorrelations/CorrelationUtils.jl @@ -1,13 +1,15 @@ """ - available_wave_vectors(sc::SampledCorrelations; bzsize=(1,1,1)) + available_wave_vectors(sc::SampledCorrelations; counts=(1,1,1)) -Returns all wave vectors for which `sc` contains exact values. `bsize` specifies -the number of Brillouin zones to be included. +Returns the grid of wave vectors for which `sc` contains exact values. +Optionally extend by a given number of `counts` along each grid axis. If the +system was not reshaped, then the number of Brillouin zones included is +`prod(counts)`. """ -function available_wave_vectors(sc::SampledCorrelations; bzsize=(1,1,1)) +function available_wave_vectors(sc::SampledCorrelations; counts=(1,1,1)) Ls = sc.sys_dims offsets = map(L -> isodd(L) ? 1 : 0, Ls) - up = Ls .* bzsize + up = Ls .* counts hi = map(L -> L - div(L, 2), up) .- offsets lo = map(L -> 1 - div(L, 2), up) .- offsets diff --git a/src/SampledCorrelations/DataRetrieval.jl b/src/SampledCorrelations/DataRetrieval.jl index 2c08e46c6..1b89fd46d 100644 --- a/src/SampledCorrelations/DataRetrieval.jl +++ b/src/SampledCorrelations/DataRetrieval.jl @@ -70,7 +70,8 @@ end contains_dynamic_correlations(sc) = !isnan(sc.Δω) -# Documented under intensities function for LSWT. +# Documented under intensities function for LSWT. TODO: As a hack, this function +# is also being used as the back-end to intensities_static. function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT) if !isnothing(kT) && kT <= 0 error("Positive `kT` required for classical-to-quantum corrections, or set `kT=nothing` to disable.") @@ -92,9 +93,11 @@ function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT # Prepare memory and configuration variables for actual calculation qpts = Base.convert(AbstractQPoints, qpts) + qs_reshaped = [to_reshaped_rlu(sc, q) for q in qpts.qs] + ffs = sc.measure.formfactors[1, :] # FIXME intensities = zeros(eltype(sc.measure), isnan(sc.Δω) ? 1 : length(ωs), length(qpts.qs)) # N.B.: Inefficient indexing order to mimic LSWT - q_idx_info = pruned_wave_vector_info(sc, qpts.qs) + q_idx_info = pruned_wave_vector_info(sc, qs_reshaped) crystal = @something sc.origin_crystal sc.crystal NCorr = Val{size(sc.data, 1)}() NAtoms = Val{size(sc.data, 2)}() @@ -102,6 +105,9 @@ function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT # Intensities calculation intensities_rounded!(intensities, sc.data, sc.crystal, sc.measure, ffs, q_idx_info, ωidcs, NCorr, NAtoms) + # Convert to a q-space density in original (not reshaped) RLU. + intensities .*= det(sc.crystal.recipvecs) / det(crystal.recipvecs) + # Post-processing steps for dynamical correlations if contains_dynamic_correlations(sc) # Convert time axis to a density. @@ -120,7 +126,6 @@ function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT intensities = reshape(intensities, length(ωs), size(qpts.qs)...) - # TODO: Refactor this logic so that return value is uniform. return if contains_dynamic_correlations(sc) Intensities(crystal, qpts, collect(ωs), intensities) else @@ -202,6 +207,3 @@ end =# - - - diff --git a/src/SampledCorrelations/SampledCorrelations.jl b/src/SampledCorrelations/SampledCorrelations.jl index 9c8804ab9..284548b0d 100644 --- a/src/SampledCorrelations/SampledCorrelations.jl +++ b/src/SampledCorrelations/SampledCorrelations.jl @@ -3,7 +3,7 @@ mutable struct SampledCorrelations const data :: Array{ComplexF64, 7} # Raw SF with sublattice indices (ncorrs × natoms × natoms × sys_dims × nω) const M :: Union{Nothing, Array{Float64, 7}} # Running estimate of (nsamples - 1)*σ² (where σ² is the variance of intensities) const crystal :: Crystal # Crystal for interpretation of q indices in `data` - const origin_crystal :: Union{Nothing,Crystal} # Original user-specified crystal (if different from above) -- needed for FormFactor accounting + const origin_crystal :: Union{Nothing,Crystal} # Original user-specified crystal (if different from above) const Δω :: Float64 # Energy step size (could make this a virtual property) measure :: MeasureSpec # Observable, correlation pairs, and combiner diff --git a/src/SpinWaveTheory/DispersionAndIntensities.jl b/src/SpinWaveTheory/DispersionAndIntensities.jl index eb648dfc0..6c2684b15 100644 --- a/src/SpinWaveTheory/DispersionAndIntensities.jl +++ b/src/SpinWaveTheory/DispersionAndIntensities.jl @@ -154,7 +154,7 @@ function intensities_bands(swt::SpinWaveTheory, qpts; kT=0, with_negative=false) # Number of (magnetic) atoms in magnetic cell @assert sys.dims == (1,1,1) - Na = length(eachsite(sys)) + Na = nsites(sys) # Number of chemical cells in magnetic cell Ncells = Na / natoms(cryst) # Number of quasiparticle modes diff --git a/src/SpinWaveTheory/SpinWaveTheory.jl b/src/SpinWaveTheory/SpinWaveTheory.jl index 8dda1fd84..3e37804f9 100644 --- a/src/SpinWaveTheory/SpinWaveTheory.jl +++ b/src/SpinWaveTheory/SpinWaveTheory.jl @@ -43,7 +43,7 @@ function SpinWaveTheory(sys::System; measure::Union{Nothing, MeasureSpec}, regul end measure = @something measure empty_measurespec(sys) - if length(eachsite(sys)) != prod(size(measure.observables)[2:5]) + if nsites(sys) != prod(size(measure.observables)[2:5]) error("Size mismatch. Check that measure is built using consistent system.") end @@ -162,7 +162,7 @@ end # rotating these into the local reference frame determined by the ground state. function swt_data(sys::System{N}, measure) where N # Calculate transformation matrices into local reference frames - Na = length(eachsite(sys)) + Na = nsites(sys) Nobs = size(measure.observables, 1) observables = reshape(measure.observables, Nobs, Na) @@ -228,7 +228,7 @@ end # Compute Stevens coefficients in the local reference frame function swt_data(sys::System{0}, measure) - Na = length(eachsite(sys)) + Na = nsites(sys) Nobs = size(measure.observables, 1) # Operators for rotating vectors into local frame diff --git a/src/Spiral/SpinWaveTheorySpiral.jl b/src/Spiral/SpinWaveTheorySpiral.jl index 89f87fa0a..e4d1a5b11 100644 --- a/src/Spiral/SpinWaveTheorySpiral.jl +++ b/src/Spiral/SpinWaveTheorySpiral.jl @@ -189,9 +189,9 @@ function intensities_bands(sswt::SpinWaveTheorySpiral, qpts; kT=0) # TODO: branc # Number of atoms in magnetic cell @assert sys.dims == (1,1,1) - Na = length(eachsite(sys)) + Na = nsites(sys) # Number of chemical cells in magnetic cell - Ncells = Na / natoms(cryst) # TODO check invariance + Ncells = Na / natoms(cryst) # Number of quasiparticle modes L = nbands(swt) # Number of wavevectors diff --git a/src/Spiral/SpiralEnergy.jl b/src/Spiral/SpiralEnergy.jl index e508cbd2c..4e8a89dc1 100644 --- a/src/Spiral/SpiralEnergy.jl +++ b/src/Spiral/SpiralEnergy.jl @@ -33,7 +33,7 @@ The [`spiral_energy`](@ref) divided by the number of sites in `sys`. The special case ``𝐤 = 0`` yields a result identical to [`energy_per_site`](@ref). """ function spiral_energy_per_site(sys::System{0}; k, axis) - return spiral_energy(sys; k, axis) / length(eachsite(sys)) + return spiral_energy(sys; k, axis) / nsites(sys) end function spiral_energy_and_gradient_aux!(dEds, sys::System{0}; k, axis) diff --git a/src/Sunny.jl b/src/Sunny.jl index 181012f9a..d0c07ea43 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -126,7 +126,7 @@ export BinningParameters, load_nxs, generate_mantid_script_from_binning_paramete include("deprecated.jl") export set_external_field!, set_external_field_at!, meV_per_K, dynamic_correlations, instant_correlations, intensity_formula, reciprocal_space_path, - set_spiral_order_on_sublattice!, set_spiral_order! + delta_function_kernel, set_spiral_order_on_sublattice!, set_spiral_order! isloaded(pkg::String) = any(k -> k.name == pkg, keys(Base.loaded_modules)) diff --git a/src/System/Interactions.jl b/src/System/Interactions.jl index 3ff015153..100d8b266 100644 --- a/src/System/Interactions.jl +++ b/src/System/Interactions.jl @@ -234,7 +234,7 @@ end The total system [`energy`](@ref) divided by the number of sites. """ function energy_per_site(sys::System{N}) where N - return energy(sys) / length(eachsite(sys)) + return energy(sys) / nsites(sys) end """ diff --git a/src/System/System.jl b/src/System/System.jl index 2f99b5470..5c42c23af 100644 --- a/src/System/System.jl +++ b/src/System/System.jl @@ -194,7 +194,7 @@ const Site = Union{NTuple{4, Int}, CartesianIndex{4}} end # Offset a `cell` by `ncells` -@inline offsetc(cell::CartesianIndex{3}, ncells, dims) = CartesianIndex(altmod1.(Tuple(cell) .+ Tuple(ncells), dims)) +@inline offsetc(cell::CartesianIndex{3}, delta, dims) = CartesianIndex(altmod1.(Tuple(cell) .+ Tuple(delta), dims)) # Split a site `site` into its cell and sublattice parts @inline to_cell(site) = CartesianIndex((site[1],site[2],site[3])) @@ -227,6 +227,15 @@ An iterator over all [`Site`](@ref)s in the system. """ @inline eachsite(sys::System) = CartesianIndices(sys.dipoles) +""" +nsites(sys::System) = length(eachsite(sys)) +""" +nsites(sys::System) = length(eachsite(sys)) + +# Number of (original) crystal cells in the system +ncells(sys::System) = nsites(sys) / natoms(orig_crystal(sys)) + + """ global_position(sys::System, site::Site) diff --git a/src/deprecated.jl b/src/deprecated.jl index 2427878df..aa0617368 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -85,7 +85,7 @@ Base.@deprecate instant_correlations(sys; opts...) let end Base.@deprecate intensity_formula(opts1...; opts2...) let - error("Formulas have been removed. Call `intensities` directly.") + error("Formulas have been removed. See revised Sunny docs for new interface.") end Base.@deprecate reciprocal_space_path(cryst::Crystal, qs, density) let @@ -110,6 +110,9 @@ Base.@deprecate SpinInfo(i; S, g) let i => Moment(; s=S, g) end +Base.@deprecate_binding delta_function_kernel nothing + + # REMEMBER TO ALSO DELETE: # # * view_crystal(cryst, max_dist) diff --git a/test/test_correlation_sampling.jl b/test/test_correlation_sampling.jl index 61c2543f3..a2e24480d 100644 --- a/test/test_correlation_sampling.jl +++ b/test/test_correlation_sampling.jl @@ -26,7 +26,6 @@ step!(sys, langevin) end end - # Test classical sum rule in SU(N) mode sys = simple_model_fcc(; mode=:SUN) @@ -35,13 +34,13 @@ sc = SampledCorrelations(sys; dt=0.08, energies, measure=ssf_trace(sys; apply_g=false)) Δω = sc.Δω add_sample!(sc, sys) - qgrid = Sunny.QPoints(Sunny.available_wave_vectors(sc)[:]) + qgrid = Sunny.available_wave_vectors(sc)[:] Δq³ = 1/prod(sys.dims) # Fraction of a BZ Sqw = intensities(sc, qgrid; energies=:available_with_negative, kT=nothing) expected_sum_rule = Sunny.norm2(sys.dipoles[1]) # S^2 classical sum rule @test isapprox(sum(Sqw.data) * Δq³ * Δω, expected_sum_rule; atol=1e-12) - # Test classical sum rule in dipole mode + # Test classical sum rule in dipole mode sys = simple_model_fcc(; mode=:dipole) thermalize_simple_model!(sys; kT=0.1) sc = SampledCorrelations(sys; dt=0.08, energies, measure=ssf_trace(sys; apply_g=false)) @@ -57,7 +56,6 @@ total_intensity_unpolarized = sum(Sqw_perp.data) @test total_intensity_unpolarized < total_intensity_trace - # Test diagonal elements are approximately real (at one wave vector) function ssf_diag_imag(sys::System{N}; apply_g=false) where N return ssf_custom(sys; apply_g) do _, ssf @@ -65,36 +63,34 @@ end end sc.measure = ssf_diag_imag(sys; apply_g=false) - is = intensities(sc, Sunny.QPoints([[0.25, 0.5, 0]]); energies=:available, kT=nothing) - is.data - @test sum(is.data) < 1e-14 - + res = intensities(sc, [[0.25, 0.5, 0]]; energies=:available, kT=nothing) + res.data + @test sum(res.data) < 1e-14 # Test form factor correction works and is doing something. formfactors = [1 => FormFactor("Fe2")] sc.measure = ssf_trace(sys; apply_g=false, formfactors) - is = intensities(sc, qgrid; energies=:available_with_negative, kT=nothing) - total_intensity_ff = sum(is.data) + res = intensities(sc, qgrid; energies=:available_with_negative, kT=nothing) + total_intensity_ff = sum(res.data) @test total_intensity_ff != total_intensity_trace - # Test static from dynamic intensities working sc.measure = ssf_trace(sys; apply_g=false) - is_static = intensities_static(sc, qgrid; kT=nothing) - total_intensity_static = sum(is_static.data) + res_static = intensities_static(sc, qgrid; kT=nothing) + total_intensity_static = sum(res_static.data) @test isapprox(total_intensity_static, total_intensity_trace * sc.Δω; atol=1e-9) # Order of summation can lead to very small discrepancies # Test quantum-to-classical increases intensity - is_static_c2q = intensities_static(sc, qgrid; kT=0.1) - total_intensity_static_c2q = sum(is_static_c2q.data) - @test total_intensity_static_c2q > total_intensity_static + res_static_c2q = intensities_static(sc, qgrid; kT=0.1) + total_intensity_static_c2q = sum(res_static_c2q.data) + @test total_intensity_static_c2q > total_intensity_static # Test static intensities working sys = simple_model_fcc(; mode=:dipole) thermalize_simple_model!(sys; kT=0.1) - ic = SampledCorrelationsStatic(sys; measure=ssf_trace(sys; apply_g=false)) - add_sample!(ic, sys) - true_static_vals = intensities_static(ic, qgrid) + res = SampledCorrelationsStatic(sys; measure=ssf_trace(sys; apply_g=false)) + add_sample!(res, sys) + true_static_vals = intensities_static(res, qgrid) true_static_total = sum(true_static_vals.data) @test isapprox(true_static_total / prod(sys.dims), 1.0; atol=1e-12) end @@ -106,7 +102,7 @@ end randomize_spins!(sys) # Set up Langevin sampler. - dt_langevin = 0.07 + dt_langevin = 0.07 langevin = Langevin(dt_langevin; damping=0.1, kT=0.1723) measure = ssf_trace(sys) @@ -141,10 +137,10 @@ end sc = SampledCorrelations(sys; energies=range(0.0, 5.5, 10), dt=0.12, measure=ssf_perp(sys)) add_sample!(sc, sys) - qs = Sunny.QPoints([[0.0, 0.0, 0.0], [-0.2, 0.4, -0.1]]) - is = intensities(sc, qs; energies=:available, kT=nothing) + qs = [[0.0, 0.0, 0.0], [-0.2, 0.4, -0.1]] + res = intensities(sc, qs; energies=:available, kT=nothing) # Compare with reference - data_golden = [33.52245944537883 31.523781055757002; 16.76122972268928 16.188214427928443; 1.6337100022968968e-14 5.3112747876921045; 1.590108516238891e-13 1.8852219123773621; -1.5194916032483394e-14 0.06400012935688847; -6.217248937900877e-15 -0.01803103943766904; 7.863406895322359e-14 -0.04445301974061088; 7.014086281484114e-14 0.0025512102338097653; 3.195591939212742e-14 -0.02515685630480813; 3.6201681383269686e-14 0.023924996100518413] - @test is.data ≈ data_golden + data_golden = [33.52245944537883 31.523781055757002; 16.76122972268928 16.188214427928443; 1.6337100022968968e-14 5.3112747876921045; 1.590108516238891e-13 1.8852219123773621; -1.5194916032483394e-14 0.06400012935688847; -6.217248937900877e-15 -0.01803103943766904; 7.863406895322359e-14 -0.04445301974061088; 7.014086281484114e-14 0.0025512102338097653; 3.195591939212742e-14 -0.02515685630480813; 3.6201681383269686e-14 0.023924996100518413] + @test res.data ≈ data_golden end diff --git a/test/test_ewald.jl b/test/test_ewald.jl index da504e2c3..9be398c48 100644 --- a/test/test_ewald.jl +++ b/test/test_ewald.jl @@ -17,14 +17,14 @@ latvecs = lattice_vectors(1,1,1,90,90,90) positions = [[0,0,0]] cryst = Crystal(latvecs, positions) - infos = [1 => Moment(s=1, g=1)] - sys = System(cryst, infos, :dipole) + moments = [1 => Moment(s=1, g=1)] + sys = System(cryst, moments, :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, infos, :dipole; dims=(2, 3, 4)) + sys = System(cryst, moments, :dipole; dims=(2, 3, 4)) enable_dipole_dipole!(sys, 1.0) @test isapprox(energy_per_site(sys), -1/6; atol=1e-13) @@ -33,12 +33,12 @@ positions = [[0,0,0], [0.1,0,0], [0.6,0.4,0.5]] cryst = Crystal(latvecs, positions) Random.seed!(0) # Don't have sys.rng yet - infos = [ + moments = [ 1 => Moment(s=1, g=rand(3,3)), 2 => Moment(s=3/2, g=rand(3,3)), 3 => Moment(s=2, g=rand(3,3)), ] - sys = System(cryst, infos, :dipole) + sys = System(cryst, moments, :dipole) enable_dipole_dipole!(sys, 1.0) randomize_spins!(sys) diff --git a/test/test_intensities_setup.jl b/test/test_intensities_setup.jl index 2befbf6f3..d6eea6d0b 100644 --- a/test/test_intensities_setup.jl +++ b/test/test_intensities_setup.jl @@ -1,14 +1,16 @@ +# TODO: Optimize performance, perhaps via `unitary_irrep_for_rotation`. Runtime +# is 1.1s on Sunny 0.7. @testitem "Magnetization Observables" begin using LinearAlgebra positions = [[0.,0,0], [0.1, 0.2, 0.3]] cryst = Crystal(I(3), positions, 1) g1 = [0.8 3.2 5.1; -0.3 0.6 -0.1; 1.0 0. 0.] - g2 = [1.8 2.2 -3.1; -0.3 2.6 0.1; 1.0 0. 0.3] + g2 = [1.8 2.2 -3.1; -0.3 2.6 0.1; 1.0 0. 0.3] for mode = [:SUN, :dipole] - infos = [1 => Moment(s=3/2, g=g1), 2 => Moment(s=(mode == :SUN ? 3/2 : 1/2), g=g2)] - sys = System(cryst, infos, mode) + moments = [1 => Moment(s=3/2, g=g1), 2 => Moment(s=(mode == :SUN ? 3/2 : 1/2), g=g2)] + sys = System(cryst, moments, mode) set_dipole!(sys, [1,3,2], (1,1,1,1)) set_dipole!(sys, [3,4,5], (1,1,1,2)) @@ -16,8 +18,8 @@ # Dipole magnetization observables (classical) sc = SampledCorrelationsStatic(sys; measure=ssf_custom((q, ssf) -> ssf, sys; apply_g=true)) add_sample!(sc, sys) - is = intensities_static(sc, Sunny.QPoints([[0,0,0]])) - corr_mat = is.data[1] + res = intensities_static(sc, [[0,0,0]]) + corr_mat = res.data[1] # Compute magnetization correlation "by hand", averaging over sites mag_corr = sum([sys.gs[i] * sys.dipoles[i] * (sys.gs[j] * sys.dipoles[j])' for i = 1:2, j = 1:2]) / Sunny.natoms(cryst) @@ -26,10 +28,10 @@ # For spin wave theory, check that `apply_g=true` is equivalent to # setting `apply_g` and manually contracting spin indices with the - # g-tensor. This only works when g is homogeneous. TODO: Test with - # inhomogeneous g-tensors. - infos_homog = [1 => Moment(s=3/2, g=g1), 2 => Moment(s=(mode == :SUN ? 3/2 : 1/2), g=g1)] - sys_homog = System(cryst, infos_homog, mode) + # g-tensor. This only works when g is equal among sites. TODO: Test with + # anisotropic g-tensors. + moments_homog = [1 => Moment(s=3/2, g=g1), 2 => Moment(s=(mode == :SUN ? 3/2 : 1/2), g=g1)] + sys_homog = System(cryst, moments_homog, mode) measure = ssf_custom((q, ssf) -> ssf, sys_homog; apply_g=false) swt = SpinWaveTheory(sys_homog; measure) @@ -60,46 +62,82 @@ end @test all(abs.(vals[2:end]) .< 1e-12) end -@testitem "Polyatomic sum rule" begin - sys = System(Sunny.diamond_crystal(), [1 => Moment(s=1/2, g=2)], :SUN; dims=(4, 1, 1), seed=1) +@testitem "Sum rule with reshaping" begin + s = 1/2 + g = 2.3 + cryst = Sunny.diamond_crystal() + dims = (1, 2, 1) # TODO: why not (1,1,1)? + sys = System(cryst, [1 => Moment(; s, g)], :SUN; dims) randomize_spins!(sys) - sc = SampledCorrelations(sys; dt=0.8, energies=range(0.0, 1.0, 3), measure=ssf_trace(sys; apply_g=true)) + sc = SampledCorrelationsStatic(sys; measure=ssf_trace(sys; apply_g=true)) add_sample!(sc, sys) - sum_rule_ixs = [1, 4, 6] # indices for zz, yy, xx - sub_lat_sum_rules = sum(sc.data[sum_rule_ixs,:,:,:,:,:,:], dims=[1,4,5,6,7])[1,:,:,1,1,1,1] - - Δq³ = 1/prod(sys.dims) # Fraction of a BZ - n_all_ω = size(sc.data, 7) - # Intensities in sc.data are a density in q, but already integrated over dω - # bins, and then scaled by n_all_ω. Therefore, we need the factor below to - # convert the previous sum to an integral. (See same logic in - # intensities_interpolated function.) - sub_lat_sum_rules .*= Δq³ / n_all_ω - - # SU(N) sum rule for S = 1/2: - # ⟨∑ᵢSᵢ²⟩ = 3/4 on every site, but because we're classical, we - # instead compute ∑ᵢ⟨Sᵢ⟩² = (1/2)^2 = 1/4 since the ⟨Sᵢ⟩ form a vector with - # length (1/2). Since the actual observables are the magnetization M = gS, we - # need to include the g factor. This is the equal-space-and-time correlation value: - gS_squared = (2 * 1/2)^2 - - expected_sum = gS_squared - # This sum rule should hold for each sublattice, independently, and only - # need to be taken over a single BZ (which is what sc.data contains) to hold: - [sub_lat_sum_rules[i,i] for i in 1:Sunny.natoms(sc.crystal)] ≈ expected_sum * ones(ComplexF64,Sunny.natoms(sc.crystal)) - - # The polyatomic sum rule demands going out 4 BZ's for the diamond crystal - # since there is an atom at relative position [1/4, 1/4, 1/4]. It also - # requires integrating over the full sampling frequency range, in this - # case by going over both positive and negative energies. - nbzs = (4, 4, 4) - qs = Sunny.available_wave_vectors(sc; bzsize=nbzs) - is = intensities(sc, Sunny.QPoints(qs[:]); energies=:available_with_negative, kT=nothing) - calculated_sum = sum(is.data) * Δq³ * sc.Δω - - # This tests that `negative_energies = true` spans exactly one sampling frequency - expected_multi_BZ_sum = gS_squared * prod(nbzs) # ⟨S⋅S⟩ - expected_multi_BZ_sum_times_natoms = expected_multi_BZ_sum * Sunny.natoms(sc.crystal) # Nₐ×⟨S⋅S⟩ - @test calculated_sum ≈ expected_multi_BZ_sum_times_natoms + # For the diamond cubic crystal, reciprocal space is periodic over a + # distance of 4 BZs. Verify that the average intensity matches the expected + # sum rule. + counts = (4, 4, 4) + qs = Sunny.available_wave_vectors(sc.parent; counts) + res = intensities_static(sc, qs[:]) + @test sum(res.data) / length(qs) ≈ Sunny.natoms(cryst) * s^2 * g^2 + + # Repeat the same calculation for a primitive cell. + shape = [0 1 1; 1 0 1; 1 1 0] / 2 + sys_prim = reshape_supercell(sys, shape) + sys_prim = repeat_periodically(sys_prim, (4, 4, 4)) + sc_prim = SampledCorrelationsStatic(sys_prim; measure=ssf_trace(sys_prim; apply_g=true)) + add_sample!(sc_prim, sys_prim) + + qs = Sunny.available_wave_vectors(sc_prim.parent; counts) + res_prim = intensities_static(sc_prim, qs[:]) + @test sum(res_prim.data) / length(qs) ≈ Sunny.natoms(cryst) * s^2 * g^2 +end + +@testitem "Sublattice sum rule" begin + latvecs = lattice_vectors(2, 1, 1, 90, 90, 90) + cryst = Crystal(latvecs, [[0, 0, 0], [0.5, 0, 0]]; types=["A", "B"]) + s = 1/2 + g1 = 2 + g2 = 2.3 + dims = (2, 1, 1) # TODO: Why won't (1,1,1) work? + sys = System(cryst, [1 => Moment(; s, g=g1), 2 => Moment(; s, g=g2)], :dipole; dims) + randomize_spins!(sys) + + sc = SampledCorrelationsStatic(sys; measure=ssf_trace(sys)) + qs = Sunny.available_wave_vectors(sc.parent; counts=(2, 1, 1)) + add_sample!(sc, sys) + res = intensities_static(sc, qs[:]) + @test sum(res.data) / length(qs) ≈ s^2 * (g1^2 + g2^2) + + # Just atoms A + formfactors = [1 => one(FormFactor), 2 => zero(FormFactor)] + sc = SampledCorrelationsStatic(sys; measure=ssf_trace(sys; formfactors)) + add_sample!(sc, sys) + res = intensities_static(sc, qs[:]) + @test sum(res.data) / length(qs) ≈ s^2 * g1^2 + + # Just atoms B + formfactors = [1 => zero(FormFactor), 2 => one(FormFactor)] + sc = SampledCorrelationsStatic(sys; measure=ssf_trace(sys; formfactors)) + add_sample!(sc, sys) + res = intensities_static(sc, qs[:]) + @test sum(res.data) / length(qs) ≈ s^2 * g2^2 +end + +@testitem "Dynamic sum rule" begin + s = 3/2 + g = 2.3 + cryst = Sunny.cubic_crystal() + sys = System(cryst, [1 => Moment(; s, g)], :dipole) + set_exchange!(sys, 1.0, Bond(1, 1, [1, 0, 0])) + randomize_spins!(sys) + + sc = SampledCorrelations(sys; dt=0.1, energies=range(0.0, 1.0, 3), measure=ssf_trace(sys)) + add_sample!(sc, sys) + + qs = Sunny.available_wave_vectors(sc) + res = intensities(sc, qs[:]; energies=:available_with_negative, kT=nothing) + + # Integrate over energies to get static intensity, then average over sampled + # q values to get sum rule. + @test sum(res.data) * sc.Δω / length(qs) ≈ (s*g)^2 end diff --git a/test/test_lswt.jl b/test/test_lswt.jl index 4a2f1f06d..358b910b8 100644 --- a/test/test_lswt.jl +++ b/test/test_lswt.jl @@ -163,8 +163,7 @@ end function test_biquad(mode, q, s) # System - infos = [1 => Moment(; s, g=2)] - sys = System(cryst, infos, mode; dims=(2, 2, 2)) + sys = System(cryst, [1 => Moment(; s, g=2)], mode; dims=(2, 2, 2)) α = -0.4π J = 1.0 JL, JQ = J * cos(α), J * sin(α) / s^2 @@ -464,8 +463,8 @@ end function build_system(R, D1, D2, J, K1, K2, h, g) latvecs = lattice_vectors(1, 1, 1, 92, 93, 94) cryst = Crystal(latvecs, [[0,0,0], [0.4,0,0]]; types=["A", "B"]) - infos = [1 => Moment(s=1, g=2), 2 => Moment(s=2, g=R*g*R')] - sys = System(cryst, infos, :dipole; seed=101) + moments = [1 => Moment(s=1, g=2), 2 => Moment(s=2, g=R*g*R')] + sys = System(cryst, moments, :dipole; seed=101) set_onsite_coupling!(sys, S -> S'*R*(D1+D1')*R'*S, 1) set_onsite_coupling!(sys, S -> S'*R*(D2+D2')*R'*S, 2)