From a31ed309723a4846d3419d1a200588661ff18d59 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 13 Dec 2023 10:01:17 -0500 Subject: [PATCH 01/11] Implement KPM Capitalize KPM. Implement KPM intensity formula Add KPM element contraction. Implement KPM intensity formula (with Sam's help). Implement quasi-Hermitian Lanczos. Comment test. Finish integrating QH Lanczos. Fix kpm test Add plotting files. fix plotting2d Fix test Fix sqrt(M) factor. KPM fix scaling of Bose factor and regularization. Fix repetition post refactor. Minor edits to bring in line with new intensities. Squashed commit bringing up to date with Sam's branch commit f001c5ec324473126551c1f7457f1b546e559a9c Author: Sam Quinn Date: Sat Nov 4 11:45:44 2023 -0400 update commit 69ca7190971cac19be991f9f57095545ea632638 Author: Sam Quinn Date: Thu Nov 2 10:53:47 2023 -0400 update intensity_formula_kpm commit cc344468ed862eee33061344e0f0bde0936db903 Author: Harry Lane Date: Wed Nov 1 11:18:47 2023 +0000 KPM fix scaling of Bose factor and regularization. commit e3ea347bdf1cac26dd737ecdd506f86d86dfebc7 Author: Harry Lane Date: Mon Oct 30 15:26:55 2023 +0000 Remove symmetrization to get offdiag corrs. commit 1fc1790aff1fa675dc97041a9a4fddfbd7596cf3 Author: Sam Quinn Date: Wed Nov 1 11:13:06 2023 -0400 positive and negative frequency LSWT commit 6558c6ee67df942e265557053ae8dacd8ef16b76 Author: Sam Quinn Date: Tue Oct 31 12:15:43 2023 -0400 Add regularization options commit 72e056dd2b8f4e7f7b4ec71d9054e960279dde5f Author: Harry Lane Date: Mon Oct 30 15:11:27 2023 +0000 Introduce cubic regularization. commit cd6347e01d2bf4391a0fe04e3e904bcc98997b08 Author: Harry Lane Date: Mon Oct 30 10:31:19 2023 +0000 Fix sqrt(M) factor. commit 89123f2710a6f74693118878a8103ef9ffdfc16b Author: [Harry Lane] <[hlane34@gatech.edu]> Date: Mon Sep 11 23:16:50 2023 -0400 Fix test commit 9f522ca7a5ebd2136aefb8b7ff8b75957fd593b8 Author: [Harry Lane] <[hlane34@gatech.edu]> Date: Thu Sep 7 12:18:50 2023 -0400 fix plotting2d commit 39ed498bac677e08edd67d1941b287294462e02b Author: [Harry Lane] <[hlane34@gatech.edu]> Date: Thu Sep 7 10:03:47 2023 -0400 Add plotting files. commit 97549f46ace5f8aad57b68a884a7adc9882a2bbf Author: [Harry Lane] <[hlane34@gatech.edu]> Date: Tue Sep 5 22:29:50 2023 -0400 Fix kpm test commit 8694b091a610d8c16420d87bd85590f98c53b7de Author: [Harry Lane] <[hlane34@gatech.edu]> Date: Tue Sep 5 18:40:26 2023 -0400 Finish integrating QH Lanczos. commit 09ac14fd53e98087ebc75d5e4f31717fcf02093c Author: [Harry Lane] <[hlane34@gatech.edu]> Date: Wed Aug 30 20:39:39 2023 -0400 Comment test. commit 191f2f27f3decdb0a72e80910d20040dafdc156a Author: [Harry Lane] <[hlane34@gatech.edu]> Date: Tue Aug 29 22:24:35 2023 -0400 Implement quasi-Hermitian Lanczos. commit e1a13756aa1a388bd37fd37ffcc24005e8ac4ad2 Author: [Harry Lane] <[hlane34@gatech.edu]> Date: Fri Aug 11 15:29:34 2023 -0400 Implement KPM intensity formula (with Sam's help). commit be51cda0a41691e3e7025c8a0b4988650496b5f7 Author: [Harry Lane] <[hlane34@gatech.edu]> Date: Thu Aug 10 13:31:50 2023 -0400 Add KPM element contraction. commit 824ecf4414e0a14026d3dc34ff5445e891dbf67e Author: Sam Quinn Date: Fri Aug 4 15:15:40 2023 -0400 Implement KPM intensity formula commit 6f29102a8825b1692ab720903093df3b906684ec Author: [Harry Lane] <[hlane34@gatech.edu]> Date: Wed Aug 2 10:55:19 2023 -0400 Add kpm dipole mode. commit 1ce57e42b97bd33d34ca1b94d8bdf533b8f45124 Author: [Harry Lane] <[hlane34@gatech.edu]> Date: Tue Aug 1 12:45:23 2023 -0400 Capitalize KPM. commit 951722ce8dd115a4c405baa53d50c16c727f2e95 Author: [Harry Lane] <[hlane34@gatech.edu]> Date: Tue Aug 1 12:14:56 2023 -0400 KPM Implementation. Remove redundant `intensity_formula` patch `intensity_formula_kpm` to not depend on omega Add susceptibility regularization style Rebased onto inhomog_lswt and fix ref to rots Take max(abs()) for Lanczos. --- Project.toml | 1 + examples/extra/chebyshev_usage.jl | 124 ++++++ examples/extra/plotting2d.jl | 124 ++++++ examples/extra/plotting_2d.jl | 145 +++++++ src/Intensities/ElementContraction.jl | 12 + src/Intensities/LinearSpinWaveIntensities.jl | 4 +- src/SpinWaveTheory/Chebyshev.jl | 83 ++++ .../DispersionAndIntensities.jl | 17 +- src/SpinWaveTheory/HamiltonianSUN.jl | 2 +- src/SpinWaveTheory/KPM.jl | 376 ++++++++++++++++++ src/SpinWaveTheory/Lanczos.jl | 76 +++- src/SpinWaveTheory/kpm.jl | 44 -- src/Sunny.jl | 13 +- test/test_lswt.jl | 72 +++- 14 files changed, 1010 insertions(+), 83 deletions(-) create mode 100644 examples/extra/chebyshev_usage.jl create mode 100644 examples/extra/plotting2d.jl create mode 100644 examples/extra/plotting_2d.jl create mode 100644 src/SpinWaveTheory/Chebyshev.jl create mode 100644 src/SpinWaveTheory/KPM.jl delete mode 100644 src/SpinWaveTheory/kpm.jl diff --git a/Project.toml b/Project.toml index 1c8f91d23..c14717117 100644 --- a/Project.toml +++ b/Project.toml @@ -17,6 +17,7 @@ Optim = "429524aa-4258-5aef-a3af-852621145aeb" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RowEchelon = "af85af4c-bcd5-5d23-b03a-a909639aa875" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Spglib = "f761d5c5-86db-4880-b97f-9680a7cccfb5" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/examples/extra/chebyshev_usage.jl b/examples/extra/chebyshev_usage.jl new file mode 100644 index 000000000..ff4e42bf1 --- /dev/null +++ b/examples/extra/chebyshev_usage.jl @@ -0,0 +1,124 @@ +using Sunny, GLMakie, LinearAlgebra + +# Make the internal functions public for the purposes of illustration +import Sunny: cheb_coefs, cheb_eval, jackson_kernel + +# Specify the function we wish to approximate and related parameters +bose_reg(x, a) = 1/((2cosh(x/a))^a - 1) # "Regularized" Bose function +bounds = (-1.0, 1.0) # Function domain +nsamps = 4000 # Number of abscissa points used to sample the function (within bounds) +npolys = 1000 # Number of polynomials to use in approximation + +# Calculate the coefficients +a = 0.05 +coefs = cheb_coefs(npolys, nsamps, x -> bose_reg(x, a), bounds) + +# The coefficients may now be used to evaluate the approximant at arbitrary x +# within the bounds +x = 0.2 +y_exact = bose_reg(x, a) +y_approx = cheb_eval(x, bounds, coefs) + +# If we wish to apply the Jackson kernel to the coefficients, this may be passed +# as an additional argument. +kernel = jackson_kernel(npolys) +y_approx = cheb_eval(x, bounds, coefs; kernel) + +# Consider now an example of a function not defined on the interval [-1, 1]. +# We'll set the domain of our function to the interval [0, 10]. +bounds = (0.0, 10.0) + +# Define the function itself (a couple of Gaussian bumps) +somefunc(x) = 0.5exp(-0.5((x - 2.0)/0.05)^2) + exp(-0.5((x - 7.0)/0.1)^2) + +# Calculate the coefficients +coefs = cheb_coefs(npolys, nsamps, somefunc, bounds) + +# Construct the reference and approximation +xs_ref = range(bounds[1], bounds[2], 2000) +xs_rec = range(bounds[1], bounds[2], 400) +ref = map(somefunc, xs_ref) +rec = map(x -> cheb_eval(x, bounds, coefs), xs_rec) + +# Plot the results +begin + fig = Figure() + ax = Axis(fig[1,1]; yscale=log10, ylabel=L"c_n", xlabel=L"n", title="Coefficients") + scatter!(ax, abs.(coefs); markersize=5.0) + ax = Axis(fig[1,2]; ylabel=L"f(x)", xlabel=L"x", title="Example Function") + lines!(ax, xs_ref, ref; color=(:black, 0.6), label="Reference") + scatter!(ax, xs_rec, rec; markersize=8.0, color=:red, label="Approximant") + axislegend(ax) + fig +end + + + +# Examine approximations using different numbers of polynomials. +begin + maxN = 20 # Number of polynomials to use in reconstruction -- choose between 2 and 1000 + + # Calculate the coefficients + a = 0.05 + bounds = (-1, 1) + coefs = cheb_coefs(npolys, nsamps, x -> bose_reg(x, a), bounds) + + # Create reference (evaluate original function on 1000 points) + xs_ref = range(bounds[1], bounds[2], 1000) + ref = map(x -> bose_reg(x, 0.05), xs_ref) # "Ground truth" + + # Reconstruct from coefficients, without and with Jackson kernel + xs_rec = range(bounds[1], bounds[2], 200) # Points to evaluate reconstruction + kernel = jackson_kernel(maxN) # Coefficient weightings from Jackson kernel + + rec = map(x -> cheb_eval(x, bounds, coefs; maxN), xs_rec) + jac = map(x -> cheb_eval(x, bounds, coefs; kernel, maxN), xs_rec) + + # Plot results + p = lines(xs_ref, ref; color=(:black, 0.6), label="Reference") + scatter!(xs_rec, rec; marker=:circle, markersize=10.0, color=:orange, label="Reconstruction") + scatter!(xs_rec, jac; marker=:cross, markersize=10.0, color=:magenta, label="Reconstruction (Jackson)") + axislegend(p.axis) + p +end + + +# Example of how coefficients die off with varying width parameter +begin + fig = Figure(resolution=(1200,500)) + as = [0.1, 0.05, 0.025] # Width parameter of regularized Bose function + + # Plot coefficient magnitudes for each case + numtokeep = Int64[] + stored_coefs = [] + for (n, a) in enumerate(as) + # Calculate and save coefficients + coefs = cheb_coefs(npolys, nsamps, x -> bose_reg(x, a), bounds) + push!(stored_coefs, coefs) + + # Plot + axis = Axis(fig[1,n]; yscale=log10, ylabel = n == 1 ? L"c_n" : "", xlabel=L"n") + scatter!(axis, abs.(coefs); markersize=2.0) + + # Figure out how many coefficients to use for reconstruction (clip when below certain value) + push!(numtokeep, findfirst(x -> abs(x) < 1e-5, coefs[0:2:end])) + end + + # Plot original function and reconstruction in each case + for (n, (maxN, coefs, a)) in enumerate(zip(numtokeep, stored_coefs, as)) + # Reconstruct functions from saved coefficients -- used saved number of coefficients + kernel = jackson_kernel(maxN) + rec = map( x -> cheb_eval(x, bounds, coefs; maxN), xs_rec) + jac = map( x -> cheb_eval(x, bounds, coefs; kernel, maxN), xs_rec) + + # Plot + ax = Axis(fig[2,n]; ylabel = n == 1 ? L"f(x)" : "", xlabel=L"x") + ref = map(x -> bose_reg(x, a), xs_ref) + lines!(ax, xs_ref, ref; color=(:black, 0.6), label="Reference") + scatter!(ax, xs_rec, rec; marker=:circle, markersize=5.0, color=:orange, label="Reconstruction") + scatter!(ax, xs_rec, jac; marker=:cross, markersize=5.0, color=:magenta, label="Reconstruction (Jackson)") + (n == 3) && Legend(fig[2,4], ax) + end + + fig +end \ No newline at end of file diff --git a/examples/extra/plotting2d.jl b/examples/extra/plotting2d.jl new file mode 100644 index 000000000..99a55102e --- /dev/null +++ b/examples/extra/plotting2d.jl @@ -0,0 +1,124 @@ +using LinearAlgebra +import ColorSchemes, ColorTypes + + +################################################# +# Functions for plotting on triangular plaquettes +################################################# + + +function plaquette_idcs(dims::Tuple{Int,Int,Int}) + dx, dy, dz = dims + (dz != 1) && println("Warning: Ignoring lattice c vector.") + Triple = Tuple{Int,Int,Int} + idcs = Array{Tuple{Triple,Triple,Triple},4}(undef, 2, dx + 1, dy + 1, 1) + for j ∈ 1:(dy+1) + for i ∈ 1:(dx+1) + idcs[1, i, j, 1] = ( + (mod1(i , dx), mod1(j, dy), 1), + (mod1(i+1, dx), mod1(j+1, dy), 1), + (mod1(i , dx), mod1(j+1, dy), 1), + ) + idcs[2, i, j, 1] = ( + (mod1(i, dx), mod1(j, dy), 1), + (mod1(i+1, dx), mod1(j, dy), 1), + (mod1(i+1, dx), mod1(j+1, dy), 1), + ) + end + end + return idcs +end + +plaquette_idcs(x::Array{T,4}) where T = plaquette_idcs(size(x)[1:3]) + +function plaquette_map(f::Function, x::Array{T,4}) where T + + dims = size(x) + @assert dims[3] == dims[4] == 1 "Multiple sites and multiple layers are not supported." + dims = dims[1:3] + out = zeros(Float64, 2, dims...) + idcs_all = plaquette_idcs(x) + for i ∈ CartesianIndices(dims) + idcs = idcs_all[1, i] + out[1, i] = f(x[idcs[1]..., 1], x[idcs[2]..., 1], x[idcs[3]..., 1]) + idcs = idcs_all[2, i] + out[2, i] = f(x[idcs[1]..., 1], x[idcs[2]..., 1], x[idcs[3]..., 1]) + end + out +end + + +################################################################################ +# Plotting functions +################################################################################ +function aspect_ratio(x_panel, y_panel, x_offset, y_offset, numrows, numcols; + adhoc_offset = (0.0, 0.0) +) + corners = [ + (0, 0), + numrows * y_panel + (numrows - 1) * y_offset, + numcols * x_panel + (numcols - 1) * x_offset, + numcols * x_panel + numrows * y_panel + (numrows - 1) * y_offset + (numcols - 1) * x_offset, + ] + xs = [c[1] for c ∈ corners] + ys = [c[2] for c ∈ corners] + x1, x2 = minimum(xs), maximum(xs) + adhoc_offset[1] + y1, y2 = minimum(ys), maximum(ys) + adhoc_offset[2] + return abs((x2 - x1) / (y2 - y1)) +end + + +function plot_triangular_plaquettes(f, frames; + colorscheme=ColorSchemes.RdBu, clims=(-0.5, 0.5), offset_spacing=1, + numcols=nothing, texts=nothing, force_aspect=true, text_offset = (0.0, 0.0), fig_kwargs... +) + # Consolidate lattice info and panel layout + numpanels = length(frames) + isnothing(numcols) && (numcols = numpanels) + numrows = floor(Int, (numpanels - 1) / numcols) + 1 + v₁, v₂ = [1, 0, 0], [-1/2, √3/2, 0] # Derives from lattice_vectors(a,a,c,90,90,120) + nx, ny = size(frames[1])[1:2] + v₁, v₂ = Point3f(v₁), Point3f(v₂) + x, y = [1.0, 0, 0], [0.0, 1, 0] + x_offset = offset_spacing * (x ⋅ v₁) * x + y_offset = -offset_spacing * (y ⋅ v₂) * y + x_panel = nx * v₁ + y_panel = -ny * v₂ + aspect = aspect_ratio(x_panel, y_panel, x_offset, y_offset, numrows, numcols; adhoc_offset=text_offset) + + # Set up figure + fig = Figure(; fig_kwargs...) + if force_aspect + ax = Axis(fig[1, 1:length(frames)]; aspect) + else + ax = Axis(fig[1, 1:length(frames)]; aspect=true) + end + hidespines!(ax) + hidedecorations!(ax) + + # Plot panels + plaq1(p) = Makie.Polygon(Point2f.([p, p+v₁+v₂, p+v₂])) + plaq2(p) = Makie.Polygon(Point2f.([p, p+v₁, p+v₂+v₁])) + for (i, frame) ∈ enumerate(frames) + r, c = fldmod1(i, numcols) + v₀ = (c - 1) * (x_panel + x_offset) + (r - 1) * (y_panel + y_offset) + + χ = plaquette_map(f, frame) + pgons = Makie.Polygon[] + colors = ColorTypes.RGB{Float64}[] + for r ∈ 1:nx + for c ∈ 1:ny + base = (r - 1) * v₁ + (c - 1) * v₂ + v₀ + push!(pgons, plaq1(base)) + push!(colors, get(colorscheme, χ[1, r, c, 1, 1], clims)) + push!(pgons, plaq2(base)) + push!(colors, get(colorscheme, χ[2, r, c, 1, 1], clims)) + end + end + poly!(ax, pgons; color=colors) + if !isnothing(texts) + text!(ax, v₀[1] - text_offset[1], v₀[2] - text_offset[2]; text=texts[i], fontsize=36) + end + end + return fig +end \ No newline at end of file diff --git a/examples/extra/plotting_2d.jl b/examples/extra/plotting_2d.jl new file mode 100644 index 000000000..e6d9ca966 --- /dev/null +++ b/examples/extra/plotting_2d.jl @@ -0,0 +1,145 @@ +using LinearAlgebra +import ColorSchemes, ColorTypes + + +################################################################################ +# Functions for plotting on triangular plaquettes, in particular Berry curvature +################################################################################ + +function vec3_triple_product(s₁, s₂, s₃) + s₁, s₂, s₃ = normalize.((s₁, s₂, s₃)) + return s₁ ⋅ (s₂ × s₃) +end + +function sun_berry_curvature(z₁, z₂, z₃) + z₁, z₂, z₃ = normalize.((z₁, z₂, z₃)) + n₁ = z₁ ⋅ z₂ + n₂ = z₂ ⋅ z₃ + n₃ = z₃ ⋅ z₁ + return angle(n₁ * n₂ * n₃) +end + + +function plaquette_idcs(dims::Tuple{Int,Int,Int}) + dx, dy, dz = dims + (dz != 1) && println("Warning: Ignoring lattice c vector.") + Triple = Tuple{Int,Int,Int} + idcs = Array{Tuple{Triple,Triple,Triple},4}(undef, 2, dx + 1, dy + 1, 1) + for j ∈ 1:(dy+1) + for i ∈ 1:(dx+1) + idcs[1, i, j, 1] = ( + (mod1(i , dx), mod1(j, dy), 1), + (mod1(i+1, dx), mod1(j+1, dy), 1), + (mod1(i , dx), mod1(j+1, dy), 1), + ) + idcs[2, i, j, 1] = ( + (mod1(i, dx), mod1(j, dy), 1), + (mod1(i+1, dx), mod1(j, dy), 1), + (mod1(i+1, dx), mod1(j+1, dy), 1), + ) + end + end + return idcs +end + +plaquette_idcs(x::Array{T,4}) where T = plaquette_idcs(size(x)[1:3]) + +function plaquette_map(f::Function, x::Array{T,4}) where T + + dims = size(x) + @assert dims[3] == dims[4] == 1 "Multiple sites and multiple layers are not supported." + dims = dims[1:3] + out = zeros(Float64, 2, dims...) + idcs_all = plaquette_idcs(x) + for i ∈ CartesianIndices(dims) + idcs = idcs_all[1, i] + out[1, i] = f(x[idcs[1]..., 1], x[idcs[2]..., 1], x[idcs[3]..., 1]) + idcs = idcs_all[2, i] + out[2, i] = f(x[idcs[1]..., 1], x[idcs[2]..., 1], x[idcs[3]..., 1]) + end + out +end + + +################################################################################ +# Plotting functions +################################################################################ +function aspect_ratio(x_panel, y_panel, x_offset, y_offset, numrows, numcols) + corners = [ + (0, 0), + numrows * y_panel + (numrows - 1) * y_offset, + numcols * x_panel + (numcols - 1) * x_offset, + numcols * x_panel + numrows * y_panel + (numrows - 1) * y_offset + (numcols - 1) * x_offset, + ] + xs = [c[1] for c ∈ corners] + ys = [c[2] for c ∈ corners] + x1, x2 = minimum(xs), maximum(xs) + y1, y2 = minimum(ys), maximum(ys) + return abs((x2 - x1) / (y2 - y1)) +end + + +function plot_chirality(xs, sys; + colorscheme=ColorSchemes.RdBu, clims=(-0.5, 0.5), offset_spacing=1, + numcols=nothing, texts=nothing, force_aspect=true, fig_kwargs... +) + # Consolidate lattice info and panel layout + numpanels = length(xs) + isnothing(numcols) && (numcols = numpanels) + numrows = floor(Int, (numpanels - 1) / numcols) + 1 + vecs = sys.crystal.latvecs + v₁, v₂ = vecs[:, 1], vecs[:, 2] + nx, ny = size(xs[1])[1:2] + v₁, v₂ = Point3f(v₁), Point3f(v₂) + x, y = [1.0, 0, 0], [0.0, 1, 0] + x_offset = offset_spacing * (x ⋅ v₁) * x + y_offset = -offset_spacing * (y ⋅ v₂) * y + x_panel = nx * v₁ + y_panel = -ny * v₂ + aspect = aspect_ratio(x_panel, y_panel, x_offset, y_offset, numrows, numcols) + + # Set up figure + fig = Figure(; fig_kwargs...) + if force_aspect + ax = Axis(fig[1, 1:length(xs)]; aspect) + else + ax = Axis(fig[1, 1:length(xs)]) + end + hidespines!(ax) + hidedecorations!(ax) + + # Plot panels + plaq1(p) = GLMakie.Polygon(Point2f.([p, p+v₁+v₂, p+v₂])) + plaq2(p) = GLMakie.Polygon(Point2f.([p, p+v₁, p+v₂+v₁])) + for (i, x) ∈ enumerate(xs) + r, c = fldmod1(i, numcols) + v₀ = (c - 1) * (x_panel + x_offset) + (r - 1) * (y_panel + y_offset) + + f = (eltype(x) == Sunny.Vec3) ? vec3_triple_product : sun_berry_curvature + χ = plaquette_map(f, x) + pgons = GLMakie.Polygon[] + colors = ColorTypes.RGB{Float64}[] + for r ∈ 1:nx + for c ∈ 1:ny + base = (r - 1) * v₁ + (c - 1) * v₂ + v₀ + push!(pgons, plaq1(base)) + push!(colors, get(colorscheme, χ[1, r, c, 1, 1], clims)) + push!(pgons, plaq2(base)) + push!(colors, get(colorscheme, χ[2, r, c, 1, 1], clims)) + end + end + poly!(ax, pgons; color=colors) + if !isnothing(texts) + text!(ax, v₀[1], v₀[2]; text=texts[i]) + end + end + return fig +end + +function plot_chirality(sys; + colorscheme=ColorSchemes.RdBu, clims=(-0.5, 0.5), offset=1, + numcols=nothing, texts=nothing, fig_kwargs... +) + x = (sys.mode == :dipole) ? sys.dipoles : sys.coherents + return plot_chirality([x], sys; colorscheme, clims, offset, numcols, texts, fig_kwargs...) +end \ No newline at end of file diff --git a/src/Intensities/ElementContraction.jl b/src/Intensities/ElementContraction.jl index 9f39d5099..27b3100e3 100644 --- a/src/Intensities/ElementContraction.jl +++ b/src/Intensities/ElementContraction.jl @@ -187,3 +187,15 @@ function intensity_formula(sc::SampledCorrelations, contractor::Contraction{T}; intensity = contract(correlations, k, contractor) end end + +function intensity_formula_kpm(swt::SpinWaveTheory, mode::Symbol; kwargs...) + contractor, string_formula = contractor_from_mode(swt, mode) + intensity_formula_kpm(swt, contractor; string_formula, kwargs...) +end + +function intensity_formula_kpm(swt::SpinWaveTheory, contractor::Contraction{T}; kwargs...) where T + intensity_formula_kpm(swt,required_correlations(contractor); return_type = T,kwargs...) do k,correlations + intensity = contract(correlations, k, contractor) + end +end + diff --git a/src/Intensities/LinearSpinWaveIntensities.jl b/src/Intensities/LinearSpinWaveIntensities.jl index bd5b4edb2..f5eabcf2d 100644 --- a/src/Intensities/LinearSpinWaveIntensities.jl +++ b/src/Intensities/LinearSpinWaveIntensities.jl @@ -65,8 +65,8 @@ function intensities_bands(swt::SpinWaveTheory, ks, formula::SpinWaveIntensityFo # Get the type parameter from the BandStructure return_type = typeof(formula).parameters[1].parameters[2] - band_dispersions = zeros(Float64,length(ks),nmodes) - band_intensities = zeros(return_type,length(ks),nmodes) + band_dispersions = zeros(Float64,length(ks),2nmodes) + band_intensities = zeros(return_type,length(ks),2nmodes) for kidx in CartesianIndices(ks) band_structure = formula.calc_intensity(swt, ks[kidx]) diff --git a/src/SpinWaveTheory/Chebyshev.jl b/src/SpinWaveTheory/Chebyshev.jl new file mode 100644 index 000000000..ea924619c --- /dev/null +++ b/src/SpinWaveTheory/Chebyshev.jl @@ -0,0 +1,83 @@ +""" + jackson_kernel(N) + +Generate a Jackson kernel for `N` Chebyshev coefficients. For use with +`cheb_eval`. +""" +function jackson_kernel(N) + Np = N + 1.0 + map(n -> (1/Np)*((Np-n)cos(n*π/Np) + sin(n*π/Np)/tan(π/Np)), 0:N-1) |> + A -> OffsetArray(A, 0:N-1) +end +@inline scale(x, bounds) = 2(x-bounds[1])/(bounds[2]-bounds[1]) - 1 +@inline unscale(x, bounds) = (x+1)*(bounds[2]-bounds[1])/2 + bounds[1] + + +# Add FFT planned version of following. Can make fully non-allocating version +# when it's time to optimize (e.g., remove internal buffer allocation) +""" + cheb_coefs(N, Nsamp, func, bounds) + +Generate `N` coefficients of the Chebyshev expansion using `Nsamp` samples of +the function `func`. Sample points are taken within the interval specified by +`bounds`, e.g., `bounds=(-1,1)`. +""" +function cheb_coefs(N, Nsamp, func, bounds) + @assert Nsamp >= N + buf = OffsetArray(zeros(Nsamp), 0:Nsamp-1) + out = OffsetArray(zeros(N), 0:N-1) + for i in 0:Nsamp-1 + x_i = cos((i+0.5)π / Nsamp) + buf[i] = func(unscale(x_i, bounds)) + end + FFTW.r2r!(buf.parent, FFTW.REDFT10) + for n in 0:N-1 + out[n] = (iszero(n) ? 0.5 : 1.0) * buf[n] / Nsamp + end + return out +end + +function view_cheb_samps(N, Nsamp, func, bounds) + @assert Nsamp >= N + xx = OffsetArray(zeros(Nsamp), 0:Nsamp-1) + buf = OffsetArray(zeros(Nsamp), 0:Nsamp-1) + out = OffsetArray(zeros(N), 0:N-1) + for i in 0:Nsamp-1 + x_i = cos((i+0.5)π / Nsamp) + xx[i] = unscale(x_i, bounds) + buf[i] = func(unscale(x_i, bounds)) + end + xx, buf +end + + + + +""" + cheb_eval(x, bounds, coefs::OffsetArray; kernel=nothing, maxN = nothing) + +Evaluate a function, specified in terms of the Chebyshev `coefs`, at point `x`. +`bounds` specifies the domain of the function. `kernel` is used to specify a +weighting of the `coefs`, using, for example, `jackson_kernel`. `maxN` specifies +how many coefficients (polynomials) to retain in the reconstruction. +""" +function cheb_eval(x, bounds, coefs::T; kernel = nothing, maxN = nothing) where T <: OffsetArray + ncoefs = length(coefs) + @assert ncoefs > 2 + maxN = isnothing(maxN) ? ncoefs : maxN + @assert maxN <= ncoefs + g = isnothing(kernel) ? _ -> 1.0 : n -> kernel[n] + + return cheb_eval_aux(scale(x, bounds), coefs, g, maxN) +end + +function cheb_eval_aux(x, c, g, maxN) + Tn2, Tn1 = 1, x + out = g(0)*c[0]*Tn2 + g(1)*c[1]*Tn1 + for n in 2:maxN-1 + Tn = 2x*Tn1 - Tn2 + out += g(n)*c[n]*Tn + Tn1, Tn2 = Tn, Tn1 + end + out +end diff --git a/src/SpinWaveTheory/DispersionAndIntensities.jl b/src/SpinWaveTheory/DispersionAndIntensities.jl index eebdcea18..3b99f8253 100644 --- a/src/SpinWaveTheory/DispersionAndIntensities.jl +++ b/src/SpinWaveTheory/DispersionAndIntensities.jl @@ -42,7 +42,8 @@ function bogoliubov!(V::Matrix{ComplexF64}, H::Matrix{ComplexF64}) # Inverse of λ gives eigenvalues of Ĩ H. We only care about the first L # eigenvalues, which are positive. A factor of 2 is needed to get the # physical quasiparticle energies. - disp = resize!(λ, L) + #disp = resize!(λ, L) + disp = λ @. disp = 2 / disp # In the special case that H(q) = H(-q) (i.e., a magnetic ordering with @@ -103,7 +104,7 @@ function dispersion(swt::SpinWaveTheory, qs) H = zeros(ComplexF64, 2nmodes, 2nmodes) V = zeros(ComplexF64, 2nmodes, 2nmodes) - disp = zeros(Float64, nmodes, length(qs)) + disp = zeros(Float64, 2nmodes, length(qs)) for (iq, q) in enumerate(qs) q_reshaped = to_reshaped_rlu(swt.sys, q) @@ -153,8 +154,8 @@ function dssf(swt::SpinWaveTheory, qs) qs = Vec3.(qs) nmodes = nbands(swt) - disp = zeros(Float64, nmodes, size(qs)...) - Sαβs = zeros(ComplexF64, 3, 3, nmodes, size(qs)...) + disp = zeros(Float64, 2nmodes, size(qs)...) + Sαβs = zeros(ComplexF64, 3, 3, 2nmodes, size(qs)...) # dssf(...) doesn't do any contraction, temperature correction, etc. # It simply returns the full Sαβ correlation matrix @@ -164,7 +165,7 @@ function dssf(swt::SpinWaveTheory, qs) for qidx in CartesianIndices(qs) q = qs[qidx] band_structure = formula.calc_intensity(swt,q) - for band = 1:nmodes + for band = 1:(2nmodes) disp[band,qidx] = band_structure.dispersion[band] Sαβs[:,:,band,qidx] .= band_structure.intensity[band] end @@ -247,7 +248,7 @@ function intensity_formula(f::Function,swt::SpinWaveTheory,corr_ix::AbstractVect H = zeros(ComplexF64, 2*nmodes, 2*nmodes) V = zeros(ComplexF64, 2*nmodes, 2*nmodes) Avec_pref = zeros(ComplexF64, Nm) - intensity = zeros(return_type, nmodes) + intensity = zeros(return_type, 2nmodes) # Expand formfactors for symmetry classes to formfactors for all atoms in # crystal @@ -363,7 +364,7 @@ function intensity_formula(f::Function,swt::SpinWaveTheory,corr_ix::AbstractVect # If there is no specified kernel, we are done: just return the # BandStructure - return BandStructure{nmodes,return_type}(disp, intensity) + return BandStructure{2*nmodes,return_type}(disp, intensity) else # Smooth kernel --> Intensity as a function of ω (or a list of ωs) return function(ω) @@ -375,4 +376,4 @@ function intensity_formula(f::Function,swt::SpinWaveTheory,corr_ix::AbstractVect end output_type = isnothing(kernel) ? BandStructure{nmodes,return_type} : return_type SpinWaveIntensityFormula{output_type}(string_formula, kernel_edep, calc_intensity) -end \ No newline at end of file +end diff --git a/src/SpinWaveTheory/HamiltonianSUN.jl b/src/SpinWaveTheory/HamiltonianSUN.jl index 4e3b476bd..cc7ab72e9 100644 --- a/src/SpinWaveTheory/HamiltonianSUN.jl +++ b/src/SpinWaveTheory/HamiltonianSUN.jl @@ -61,7 +61,7 @@ end # Set the dynamical quadratic Hamiltonian matrix in SU(N) mode. -function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_reshaped::Vec3) +function swt_hamiltonian_SUN!(H, swt::SpinWaveTheory, q_reshaped::Vec3) (; sys, data) = swt (; zeeman_operators) = data diff --git a/src/SpinWaveTheory/KPM.jl b/src/SpinWaveTheory/KPM.jl new file mode 100644 index 000000000..da0399bc4 --- /dev/null +++ b/src/SpinWaveTheory/KPM.jl @@ -0,0 +1,376 @@ +""" +This code will implement the Kernal Polynomial method to approximate the dynamical spin structure factor (DSSF). +The method takes advantage of the fact that the DSSF can be written as a matrix function and hence we can avoid +diagonalizing the Hamiltonian. Futher, implementing the KPM we do not even need to explicity construct the matrix +function, instead we only need to evaluate the application of the function to some vector. For sparce matrices this +is efficient making this approximate approach useful in systems with large unit cells. +""" + +""" + bose_function(kT, x) + +Returns the Bose occupation factor for energy, x, at temperature kT. +""" + +function bose_function(kT, x) + r = 1 / (exp(x/kT) - 1) + if isinf(r) + return 0. + end + return r +end + +""" + regularization_function(ω,σ) + +Returns a regularization factor to apply to the intensity at low energy according to a smooth approximation to a step function +with width, σ. + +""" + +function regularization_function(ω,σ) + if ω < 0 + return 0.0 + elseif 0 ≤ ω ≤ σ + return (4 - (3ω ./σ)) .*(ω.^3/σ.^3) + else + return 1.0 + end +end + +function srf(ω,σ) + if -σ ≤ ω ≤ σ + return (3(ω ./σ).^4 .- 2(ω ./σ).^6) + else + return 1.0 + end +end + + + +""" + get_all_coefficients(M, ωs, broadening, σ, kT,γ) + +Retrieves the Chebyshev coefficients up to index, M for a user-defined lineshape. A numerical regularization is applied to +treat the divergence of the dynamical correlation function at small energy. Here, σ² represents an energy cutoff scale +which should be related to the energy resolution. γ is the maximum eigenvalue used to rescale the spectrum to lie on the +interval [-1,1]. Regularization is treated using a cubic cutoff function and the negative eigenvalues are zeroed out. + +""" +function get_all_coefficients(M, ωs, broadening, σ, kT,γ;η=1.0, regularization_style) + f = if regularization_style == :cubic + (ω,x) -> regularization_function(γ*x,η*σ) * broadening(ω, x*γ, σ) * (1 + bose_function(kT, x*γ)) + elseif regularization_style == :tanh + (ω,x) -> tanh((γ*x/(η*σ))^2) * broadening(ω, x*γ, σ) * (1 + bose_function(kT, x*γ)) + elseif regularization_style == :none + (ω,x) -> broadening(ω, x*γ, σ) * (1 + bose_function(kT, x*γ)) + elseif regularization_style == :susceptibility + (ω,x) -> broadening(ω, x*γ, σ) + elseif regularization_style == :srf + (ω,x) -> srf(γ*x,η*σ) * broadening(ω, x*γ, σ) * (1 + bose_function(kT, x*γ)) + end + output = OffsetArray(zeros(M, length(ωs)), 0:M-1, 1:length(ωs)) + for i in eachindex(ωs) + output[:, i] = cheb_coefs(M, 2M, x -> f(ωs[i], x), (-1, 1)) + end + return output +end + + +""" + get_all_coefficients_legacy(M, ωs, broadening, σ, kT,γ) + +Retrieves the Chebyshev coefficients up to index, M for a user-defined lineshape. A numerical regularization is applied to +treat the divergence of the dynamical correlation function at small energy. Here, σ² represents an energy cutoff scale +which should be related to the energy resolution. γ is the maximum eigenvalue used to rescale the spectrum to lie on the +interval [-1,1]. Regularization is treated using a tanh cutoff function and the negative eigenvalues are not zeroed out. + +""" +function get_all_coefficients_legacy(M, ωs, broadening, σ, kT,γ) + f(ω, x) = tanh((x/σ)^2) * broadening(ω, x*γ, σ) * (1 + bose_function(kT, x)) + output = OffsetArray(zeros(M, length(ωs)), 0:M-1, 1:length(ωs)) + for i in eachindex(ωs) + output[:, i] = cheb_coefs(M, 2M, x -> f(ωs[i], x), (-1, 1)) + end + return output +end + +""" + apply_kernel(offset_array, kernel, M) + +Applies the Jackson kernel (defined in Chebyshev.jl) to an OffsetArray. The Jackson kernel damps the coefficients of the +Chebyshev expansion to reduce "Gibbs oscillations" or "ringing" -- numerical artifacts present due to the truncation offset_array +the Chebyshev series. It should be noted that employing the Jackson kernel comes at a cost to the energy resolution. +""" + +# +function apply_kernel(offset_array, kernel, M) + kernel === "jackson" ? offset_array .= offset_array .* jackson_kernel(M) : + return offset_array +end + +""" + kpm_dssf(swt::SpinWaveTheory, qs,ωlist,P::Int64,kT,σ,broadening; kernel) + +Calculated the Dynamical Spin Structure Factor (DSSF) using the Kernel Polynomial Method (KPM). Requires input in the form of a +SpinWaveTheory which contains System information and the rotated spin matrices. The calculation is carried out at each wavevectors +in qs for all energies appearing in ωlist. The Chebyshev expansion is taken to P terms and the lineshape is specified by the user- +defined function, broadening. kT is required for the calcualation of the bose function and σ is the width of the lineshape and +defines the low energy cutoff σ². There is a keyword argument, kernel, which speficies a damping kernel. +""" + +function kpm_dssf(swt::SpinWaveTheory, qs,ωlist,P::Int64,kT,σ,broadening; kernel, regularization_style) + # P is the max Chebyshyev coefficient + (; sys, data) = swt + qs = Sunny.Vec3.(qs) + Nm, Ns = length(sys.dipoles), sys.Ns[1] # number of magnetic atoms and dimension of Hilbert space + Nf = sys.mode == :SUN ? Ns-1 : 1 + N=Nf+1 + nmodes = Nf*Nm + M = sys.mode == :SUN ? 1 : (Ns-1) # scaling factor (=1) if in the fundamental representation + sqrt_M = √M #define prefactors + sqrt_Nm_inv = 1.0 / √Nm #define prefactors + S = (Ns-1) / 2 + sqrt_halfS = √(S/2) #define prefactors + Ĩ = spdiagm([ones(nmodes); -ones(nmodes)]) + n_iters = 50 + Hmat = zeros(ComplexF64, 2*nmodes, 2*nmodes) + Avec_pref = zeros(ComplexF64, Nm) # initialize array of some prefactors + chebyshev_moments = OffsetArray(zeros(ComplexF64,3,3,length(qs),P),1:3,1:3,1:length(qs),0:P-1) + Sαβs = zeros(ComplexF64,3,3,length(qs),length(ωlist)) + for qidx in CartesianIndices(qs) + q = qs[qidx] + q_reshaped = Sunny.to_reshaped_rlu(swt.sys, q) + u = zeros(ComplexF64,3,2*nmodes) + if sys.mode == :SUN + swt_hamiltonian_SUN!(Hmat, swt, q_reshaped) + else + swt_hamiltonian_dipole!(Hmat, swt, q_reshaped) + end + D = 2.0*sparse(Hmat) # calculate D (factor of 2 for correspondence) + lo,hi = Sunny.eigbounds(Ĩ*D,n_iters; extend=0.25) # calculate bounds + γ=max(abs(lo),abs(hi)) # select upper bound (combine with the preceeding line later) + A = Ĩ*D / γ + # u(q) calculation) + for site = 1:Nm + # note that d is the chemical coordinates + chemical_coor = sys.crystal.positions[site] # find chemical coords + phase = exp(2*im * π * dot(q_reshaped, chemical_coor)) # calculate phase + Avec_pref[site] = sqrt_Nm_inv * phase # define the prefactor of the tS matrices + end + # calculate u(q) + if sys.mode == :SUN + for site=1:Nm + @views tS_μ = data.dipole_operators[:, :, :, site]*Avec_pref[site] + for μ=1:3 + for j=2:N + u[μ,(j-1)+(site-1)*(N-1) ]=tS_μ[j,1,μ] + u[μ,(N-1)*Nm+(j-1)+(site-1)*(N-1) ]=tS_μ[1,j,μ] + end + end + end + elseif sys.mode == :dipole + for site = 1:Nm + R=data.R_mat[site] + u[1,site]= Avec_pref[site] * sqrt_halfS * (R[1,1] + 1im * R[1,2]) + u[1,site+nmodes] = Avec_pref[site] * sqrt_halfS * (R[1,1] - 1im * R[1,2]) + u[2,site]= Avec_pref[site] * sqrt_halfS * (R[2,1] + 1im * R[2,2]) + u[2,site+nmodes] = Avec_pref[site] * sqrt_halfS * (R[2,1] - 1im * R[2,2]) + u[3,site]= Avec_pref[site] * sqrt_halfS * (R[3,1] + 1im * R[3,2]) + u[3,site+nmodes] = Avec_pref[site] * sqrt_halfS * (R[3,1] - 1im * R[3,2]) + end + + end + for β=1:3 + α0 = zeros(ComplexF64,2*nmodes) + α1 = zeros(ComplexF64,2*nmodes) + mul!(α0,Ĩ,u[β,:]) # calculate α0 + mul!(α1,A,α0) # calculate α1 + for α=1:3 + chebyshev_moments[α,β,qidx,0] = (dot(u[α,:],α0)) #removed symmetrization + chebyshev_moments[α,β,qidx,1] = (dot(u[α,:],α1)) #removed symmetrization + end + for m=2:P-1 + αnew = zeros(ComplexF64,2*nmodes) + mul!(αnew,A,α1) + @. αnew = 2*αnew - α0 + for α=1:3 + chebyshev_moments[α,β,qidx,m] = (dot(u[α,:],αnew)) #removed symmetrization + end + (α1, α0) = (αnew, α1) + end + end + ωdep = get_all_coefficients(P,ωlist,broadening,σ,kT,γ;regularization_style) + apply_kernel(ωdep,kernel,P) + for w=1:length(ωlist) + for α=1:3 + for β=1:3 + Sαβs[α,β,qidx,w] = sum(chebyshev_moments[α,β,qidx,:] .* ωdep[:,w]) + end + end + end + end + return Sαβs +end + +""" + kpm_intensities(swt::SpinWaveTheory, qs,ωlist,P::Int64,kT,σ,broadening; kernel) + +Calculated the neutron scattering intensity using the Kernel Polynomial Method (KPM). Calls KPMddsf and so takes the same parameters. +Requires input in the form of a SpinWaveTheory which contains System information and the rotated spin matrices. The calculation is +carried out at each wavevectors in qs for all energies appearing in ωlist. The Chebyshev expansion is taken to P terms and the +lineshape is specified by the user-defined function, broadening. kT is required for the calcualation of the bose function and σ is +the width of the lineshape and defines the low energy cutoff σ². There is an optional keyword argument, kernel, which speficies a +damping kernel. The default is to include no damping. +""" +function kpm_intensities(swt::SpinWaveTheory, qs, ωvals,P::Int64,kT,σ,broadening; kernel = nothing,regularization_style) + (; sys) = swt + qs = Sunny.Vec3.(qs) + Sαβs = kpm_dssf(swt, qs,ωvals,P,kT,σ,broadening;kernel,regularization_style) + num_ω = length(ωvals) + is = zeros(Float64, size(qs)..., num_ω) + for qidx in CartesianIndices(qs) + q_reshaped = Sunny.to_reshaped_rlu(sys, qs[qidx]) + q_absolute = sys.crystal.recipvecs * q_reshaped + polar_mat = Sunny.polarization_matrix(q_absolute) + is[qidx, :] = real(sum(polar_mat .* Sαβs[:,:,qidx,:],dims=(1,2))) + end + return is +end + + +struct KPMIntensityFormula{T} + P :: Int64 + kT :: Float64 + σ :: Float64 + broadening + kernel + string_formula :: String + calc_intensity :: Function +end + +function Base.show(io::IO, ::MIME"text/plain", formula::KPMIntensityFormula{T}) where T + printstyled(io, "Quantum Scattering Intensity Formula (KPM Method)\n";bold=true, color=:underline) + + formula_lines = split(formula.string_formula,'\n') + + intensity_equals = " Intensity(Q,ω) = " + println(io,"At any (Q,ω), with S = ...:") + println(io) + println(io,intensity_equals,formula_lines[1]) + for i = 2:length(formula_lines) + precursor = repeat(' ', textwidth(intensity_equals)) + println(io,precursor,formula_lines[i]) + end + println(io,"P = $(formula.P), kT = $(formula.kT), σ = $(formula.σ)") +end + +function intensity_formula_kpm(f,swt::SpinWaveTheory,corr_ix::AbstractVector{Int64}; P =50, kT=Inf,σ=0.1,broadening, kernel=nothing , return_type = Float64, string_formula = "f(Q,ω,S{α,β}[ix_q,ix_ω])",regularization_style) + # P is the max Chebyshyev coefficient + (; sys, data) = swt + Nm, Ns = length(sys.dipoles), sys.Ns[1] # number of magnetic atoms and dimension of Hilbert space + Nf = sys.mode == :SUN ? Ns-1 : 1 + N=Nf+1 + nmodes = Nf*Nm + M = sys.mode == :SUN ? 1 : (Ns-1) # scaling factor (=1) if in the fundamental representation + sqrt_M = √M #define prefactors + sqrt_Nm_inv = 1.0 / √Nm #define prefactors + S = (Ns-1) / 2 + sqrt_halfS = √(S/2) #define prefactors + sqrt_Nm_inv = 1.0 / √Nm #define prefactors + S = (Ns-1) / 2 + sqrt_halfS = √(S/2) #define prefactors + Ĩ = spdiagm([ones(nmodes); -ones(nmodes)]) + n_iters = 50 + Hmat = zeros(ComplexF64, 2*nmodes, 2*nmodes) + Avec_pref = zeros(ComplexF64, Nm) # initialize array of some prefactors + chebyshev_moments = OffsetArray(zeros(ComplexF64,3,3,P),1:3,1:3,0:P-1) + + calc_intensity = function(swt::SpinWaveTheory,q::Vec3) + q_reshaped = Sunny.to_reshaped_rlu(sys, q) + q_absolute = swt.sys.crystal.recipvecs * q_reshaped + u = zeros(ComplexF64,3,2*nmodes) + if sys.mode == :SUN + swt_hamiltonian_SUN!(Hmat, swt, q_reshaped) + else + swt_hamiltonian_dipole!(Hmat, swt, q_reshaped) + end + D = 2.0*sparse(Hmat) # calculate D (factor of 2 for correspondence) + lo,hi = Sunny.eigbounds(Ĩ*D,n_iters; extend=0.25) # calculate bounds + + γ=max(lo,hi) # select upper bound (combine with the preceeding line later) + A = Ĩ*D / γ + # u(q) calculation + for site = 1:Nm + # note that d is the chemical coordinates + chemical_coor = sys.crystal.positions[site] # find chemical coords + phase = exp(2*im * π * dot(q_reshaped, chemical_coor)) # calculate phase + Avec_pref[site] = sqrt_Nm_inv * phase # define the prefactor of the tS matrices + end + # calculate u(q) + if sys.mode == :SUN + for site=1:Nm + @views tS_μ = data.dipole_operators[:, :, :, site]*Avec_pref[site] + for μ=1:3 + for j=2:N + u[μ,(j-1)+(site-1)*(N-1) ]=tS_μ[j,1,μ] + u[μ,(N-1)*Nm+(j-1)+(site-1)*(N-1) ]=tS_μ[1,j,μ] + end + end + end + elseif sys.mode == :dipole + for site = 1:Nm + R=data.local_rotations[site] + u[1,site]= Avec_pref[site] * sqrt_halfS * (R[1,1] + 1im * R[1,2]) + u[1,site+nmodes] = Avec_pref[site] * sqrt_halfS * (R[1,1] - 1im * R[1,2]) + u[2,site]= Avec_pref[site] * sqrt_halfS * (R[2,1] + 1im * R[2,2]) + u[2,site+nmodes] = Avec_pref[site] * sqrt_halfS * (R[2,1] - 1im * R[2,2]) + u[3,site]= Avec_pref[site] * sqrt_halfS * (R[3,1] + 1im * R[3,2]) + u[3,site+nmodes] = Avec_pref[site] * sqrt_halfS * (R[3,1] - 1im * R[3,2]) + end + + end + for β=1:3 + α0 = zeros(ComplexF64,2*nmodes) + α1 = zeros(ComplexF64,2*nmodes) + mul!(α0,Ĩ,u[β,:]) # calculate α0 + mul!(α1,A,α0) # calculate α1 + for α=1:3 + chebyshev_moments[α,β,0] = (dot(u[α,:],α0)) #removed symmetrization + chebyshev_moments[α,β,1] = (dot(u[α,:],α1)) #removed symmetrization + end + for m=2:P-1 + αnew = zeros(ComplexF64,2*nmodes) + mul!(αnew,A,α1) + @. αnew = 2*αnew - α0 + for α=1:3 + chebyshev_moments[α,β,m] = (dot(u[α,:],αnew)) #removed symmetrization + end + (α1, α0) = (αnew, α1) + end + end + + return function(ωlist) + intensity = zeros(Float64,length(ωlist)) + ωdep = get_all_coefficients(P,ωlist,broadening,σ,kT,γ;regularization_style) + apply_kernel(ωdep,kernel,P) + Sαβ = Matrix{ComplexF64}(undef,3,3) + corrs = Vector{ComplexF64}(undef,num_correlations(swt.observables)) + for (iω,ω) = enumerate(ωlist) + for α=1:3 + for β=1:3 + Sαβ[α,β] = sum(chebyshev_moments[α,β,:] .* ωdep[:,iω]) + end + end + + for (ci,i) in swt.observables.correlations + (α,β) = ci.I + corrs[i] = Sαβ[α,β] + end + intensity[iω] = f(q_absolute,corrs[corr_ix]) + end + return intensity + end + end + KPMIntensityFormula{return_type}(P,kT,σ,broadening,kernel,string_formula,calc_intensity) +end diff --git a/src/SpinWaveTheory/Lanczos.jl b/src/SpinWaveTheory/Lanczos.jl index 8d9fa083d..8ee0088a8 100644 --- a/src/SpinWaveTheory/Lanczos.jl +++ b/src/SpinWaveTheory/Lanczos.jl @@ -6,7 +6,7 @@ Lanczos algorithm. `niters` should be given a value smaller than the dimension of `A`. `extend` specifies how much to shift the upper and lower bounds as a percentage of the total scale of the estimated eigenvalues. """ -function eigbounds(A, niters; extend=0.0) +function eigbounds(A, niters; extend=0.5) lo, hi = lanczos(A, niters) |> eigvals |> xs -> (first(xs), last(xs)) slack = extend*(hi-lo) return lo-slack, hi+slack @@ -16,43 +16,81 @@ end """ lanczos(A, niters) -Takes an `N`x`N`` Hermitian matrix `A` and returns a real, symmetric, +Takes an `N`x`N`` quasi-Hermitian matrix `A` and returns a real, symmetric, tridiagonal matrix `T`. `niters` sets the number of iterations used by the Lanczos algorithm. """ + function lanczos(A, niters) - @assert ishermitian(A) "Matrix must be Hermitian" N = size(A, 1) - + nmodes= Int(N/2) + Ĩ = diagm([ones(nmodes); -ones(nmodes)]) + #@assert Ĩ * A == A' * Ĩ "Matrix must be quasi-Hermitian" αs = zeros(Float64, niters) # Main diagonal βs = zeros(Float64, niters-1) # Off diagonal - buf = zeros(ComplexF64, N, 3) # Vector buffer -- don't technically need 3 columns, but more legible this way - - lanczos_aux!(αs, βs, buf, A, niters) # Call non-allocating internal Lanczos function - + buf = zeros(ComplexF64, N, 3) # Vector buffer -- don't technically need 3 columns, but more legible this way + lanczos_QH_aux!(αs, βs, buf, A, niters) # Call non-allocating internal Lanczos function return SymTridiagonal(αs, βs) end - function lanczos_aux!(αs, βs, buf, A, niters) v, vprev, w = view(buf,:,1), view(buf,:,2), view(buf,:,3) - - @label initialize_lanczos - randn!(v) - normalize!(v) - mul!(w, A, v) - αs[1] = real(w⋅v) - @. w = w - αs[1]*v + randn!(vprev) + normalize!(vprev) + mul!(w, A, vprev) + αs[1] = real(w⋅vprev) + @. w = w - αs[1]*vprev for j in 2:niters - v, vprev = vprev, v βs[j-1] = norm(w) - (abs(βs[j-1]) < 1e-12) && (@goto initialize_lanczos) + iszero(βs[j-1]) && return @. v = w/βs[j-1] mul!(w, A, v) αs[j] = real(w⋅v) @. w = w - αs[j]*v - βs[j-1]*vprev + v, vprev = vprev, v end return nothing -end \ No newline at end of file +end +function Ĩ_dot_3(w, v, n) + @views real(dot(w[1:n], v[1:n]) - dot(w[n+1:2n], v[n+1:2n])) +end +function lanczos_QH_aux!(αs, βs, buf, A, niters) + v, vprev, w = view(buf,:,1), view(buf,:,2), view(buf,:,3) + randn!(vprev) + mat_size=Int(size(A,2)/2) + mul!(w, A, vprev) + b0 = sqrt(Ĩ_dot_3(vprev, w,mat_size)) + vprev = vprev/b0 + w = w/b0 + αs[1] = Ĩ_dot_3(w, w,mat_size) + @. w = w - αs[1]*vprev + for j in 2:niters + @. v = w + mul!(w, A, v) + βs[j-1] = sqrt(Ĩ_dot_3(v, w,mat_size)) # maybe v? + iszero(βs[j-1]) && return + @. v = v/βs[j-1] + @. w = w/βs[j-1] + αs[j] = Ĩ_dot_3(w, w,mat_size) + @. w = w - αs[j]*v - βs[j-1]*vprev + v, vprev = vprev, v + end + return nothing +end + + + + +function lanczos_legacy(A, niters) + @assert ishermitian(A) "Matrix must be Hermitian" + N = size(A, 1) + + αs = zeros(Float64, niters) # Main diagonal + βs = zeros(Float64, niters-1) # Off diagonal + buf = zeros(ComplexF64, N, 3) # Vector buffer -- don't technically need 3 columns, but more legible this way + + lanczos_aux!(αs, βs, buf, A, niters) # Call non-allocating internal Lanczos function + return SymTridiagonal(αs, βs) +end diff --git a/src/SpinWaveTheory/kpm.jl b/src/SpinWaveTheory/kpm.jl deleted file mode 100644 index 8bfa19838..000000000 --- a/src/SpinWaveTheory/kpm.jl +++ /dev/null @@ -1,44 +0,0 @@ -struct KPMIntensityFormula{T} - P :: Int64 - kT :: Float64 - σ :: Float64 - broadening - kernel - string_formula :: String - calc_intensity :: Function -end - -function Base.show(io::IO, formula::KPMIntensityFormula{T}) where T - print(io,"KPMIntensityFormula{$T}") -end - -function Base.show(io::IO, ::MIME"text/plain", formula::KPMIntensityFormula{T}) where T - printstyled(io, "Quantum Scattering Intensity Formula (KPM Method)\n";bold=true, color=:underline) - - formula_lines = split(formula.string_formula,'\n') - - intensity_equals = " Intensity(Q,ω) = " - println(io,"At any (Q,ω), with S = ...:") - println(io) - println(io,intensity_equals,formula_lines[1]) - for i = 2:length(formula_lines) - precursor = repeat(' ', textwidth(intensity_equals)) - println(io,precursor,formula_lines[i]) - end - println(io,"P = $(formula.P), kT = $(formula.kT), σ = $(formula.σ)") -end - - -function intensity_formula_kpm(f::Function,swt::SpinWaveTheory,corr_ix::AbstractVector{Int64}; P =, return_type = Float64, string_formula = "f(Q,ω,S{α,β}[ix_q,ix_ω])") - - error("KPM not yet implemented") - - stuff = setup_stuff(swt) - formula = function(swt::SpinWaveTheory,q::Vec3,ω::Float64) - Sαβ = do_KPM(swt,stuff) - return f(q,ω,Sαβ[corr_ix]) - end - KPMIntensityFormula{return_type}(P,kT,σ,broadening,kernel,string_formula,formula) -end - - diff --git a/src/Sunny.jl b/src/Sunny.jl index a56e21ef1..5e2d142ea 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -19,6 +19,15 @@ import CrystalInfoFramework as CIF import Spglib import RowEchelon: rref! +# Specific to SunnyGfx +import JSON +import Colors: distinguishable_colors, RGB, Colors +import Inflate: inflate_gzip +import Random: randstring, RandomDevice + +# Specific to KPM +import SparseArrays: spzeros, sparse, spdiagm + const Vec3 = SVector{3, Float64} const Vec5 = SVector{5, Float64} const Mat3 = SMatrix{3, 3, Float64, 9} @@ -96,7 +105,9 @@ include("SpinWaveTheory/HamiltonianDipole.jl") include("SpinWaveTheory/HamiltonianSUN.jl") include("SpinWaveTheory/DispersionAndIntensities.jl") include("SpinWaveTheory/Lanczos.jl") -export SpinWaveTheory, dispersion, dssf, delta_function_kernel +include("SpinWaveTheory/Chebyshev.jl") +include("SpinWaveTheory/KPM.jl") +export SpinWaveTheory, dispersion, intensities, dssf, kpm_dssf, kpm_intensities include("SampledCorrelations/SampledCorrelations.jl") include("SampledCorrelations/CorrelationUtils.jl") diff --git a/test/test_lswt.jl b/test/test_lswt.jl index b0e9ea23a..1758e7933 100644 --- a/test/test_lswt.jl +++ b/test/test_lswt.jl @@ -81,18 +81,74 @@ end -@testitem "Lanczos Bounds" begin - using LinearAlgebra, Random - - Random.seed!(100) - A = randn(ComplexF64, 500, 500) - A = 0.5(A' + A) - lo, hi = Sunny.eigbounds(A, 20) - vals = eigvals(A) +""" +@testitem "Lanczos Bounds" begin + using LinearAlgebra + n=10 + A = rand(ComplexF64,n,n) + A = 0.5(A+A') + B = rand(ComplexF64,n,n) + B = 0.5(B+transpose(B)) + D = [A B; conj(B) conj(A)] + ϵ = minimum(eigvals(D)) + D = D-(ϵ-0.1)*I(2n) + Ĩ = diagm([ones(n); -ones(n)]) + niters=40 + H = Ĩ*D + lo, hi = Sunny.eigbounds(H, niters; extend=0.) + vals = eigvals(H) @test (abs(lo/vals[1] - 1) < 0.025) && (abs(hi/vals[end] - 1) < 0.025) end +@testitem "Chebyshev Approximations" begin + using LinearAlgebra + + function cheb_coefs_manual(M, Mq, func, bounds) + function cheb_array!(out, x) + len = length(out) + @assert len > 2 + out[0], out[1] = 1.0, x + for m = 2:len-1 + out[m] = 2x*out[m-1] - out[m-2] + end + return nothing + end + + @assert Mq > M + out = Sunny.OffsetArray(zeros(M), 0:M-1) # To avoid bringing OffsetArrays into testing environment + fp = Sunny.OffsetArray(zeros(M), 0:M-1) + T = Sunny.OffsetArray(zeros(M), 0:M-1) + for i in 0:Mq-1 + x_i = cos((i+0.5)π/Mq) + f_i = func(Sunny.unscale(x_i, bounds)) + cheb_array!(T, x_i) + @. fp += f_i * T + end + for m in 0:M-1 + out[m] = (iszero(m) ? 1.0 : 2.0 ) * fp[m] / Mq + end + return out + end + + # Set up test function + bounds = (-1, 1) + σ = 0.1 + func(x) = (1/√(2π*σ^2))*exp(-0.5*(x/σ)^2) + + # Check coefficients from Sunny function match manual calculation + M = 100 + Mq = 400 + coefs = Sunny.cheb_coefs(M, Mq, func, bounds) + coefs_ref = cheb_coefs_manual(M, Mq, func, bounds) + @test norm(coefs - coefs_ref) < 1e-12 + + # Check resulting approximation sufficiently close to ground truth + xs = -0.9:0.1:0.9 + ref = func.(xs) + rec = map(x -> Sunny.cheb_eval(x, bounds, coefs), xs) + @test norm(rec - ref) < 1e-12 +end @testitem "Single Ion" begin # Tetragonal crystal From a77f930250a71ef050ca9c8f5f0adaf785c6c8af Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 13 Dec 2023 10:01:17 -0500 Subject: [PATCH 02/11] Remove unneeded dependencies --- src/Sunny.jl | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/Sunny.jl b/src/Sunny.jl index 5e2d142ea..55c3de9dc 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -19,12 +19,6 @@ import CrystalInfoFramework as CIF import Spglib import RowEchelon: rref! -# Specific to SunnyGfx -import JSON -import Colors: distinguishable_colors, RGB, Colors -import Inflate: inflate_gzip -import Random: randstring, RandomDevice - # Specific to KPM import SparseArrays: spzeros, sparse, spdiagm From 38f02c5e61f2e6fd20adb52a3dec7b1546b33e36 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 13 Dec 2023 10:01:17 -0500 Subject: [PATCH 03/11] Update [compat] and remove export `intensities` --- Project.toml | 1 + src/Sunny.jl | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index c14717117..c3b3569d1 100644 --- a/Project.toml +++ b/Project.toml @@ -45,6 +45,7 @@ Optim = "1.7.6" Printf = "1.9" Random = "1.9" RowEchelon = "0.2.1" +SparseArrays = "1.9" SpecialFunctions = "2.1.0" Spglib = "0.8.0" StaticArrays = "1.3.4" diff --git a/src/Sunny.jl b/src/Sunny.jl index 55c3de9dc..568132a76 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -101,7 +101,7 @@ include("SpinWaveTheory/DispersionAndIntensities.jl") include("SpinWaveTheory/Lanczos.jl") include("SpinWaveTheory/Chebyshev.jl") include("SpinWaveTheory/KPM.jl") -export SpinWaveTheory, dispersion, intensities, dssf, kpm_dssf, kpm_intensities +export SpinWaveTheory, dispersion, dssf, kpm_dssf, kpm_intensities include("SampledCorrelations/SampledCorrelations.jl") include("SampledCorrelations/CorrelationUtils.jl") From 6a479ec73c0d324743a2cf16cc1c12094e47509f Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 13 Dec 2023 10:01:17 -0500 Subject: [PATCH 04/11] Undo changes to LSWT --- src/Intensities/LinearSpinWaveIntensities.jl | 4 ++-- src/SpinWaveTheory/DispersionAndIntensities.jl | 15 +++++++-------- src/Sunny.jl | 2 +- test/test_lswt.jl | 4 +--- 4 files changed, 11 insertions(+), 14 deletions(-) diff --git a/src/Intensities/LinearSpinWaveIntensities.jl b/src/Intensities/LinearSpinWaveIntensities.jl index f5eabcf2d..bd5b4edb2 100644 --- a/src/Intensities/LinearSpinWaveIntensities.jl +++ b/src/Intensities/LinearSpinWaveIntensities.jl @@ -65,8 +65,8 @@ function intensities_bands(swt::SpinWaveTheory, ks, formula::SpinWaveIntensityFo # Get the type parameter from the BandStructure return_type = typeof(formula).parameters[1].parameters[2] - band_dispersions = zeros(Float64,length(ks),2nmodes) - band_intensities = zeros(return_type,length(ks),2nmodes) + band_dispersions = zeros(Float64,length(ks),nmodes) + band_intensities = zeros(return_type,length(ks),nmodes) for kidx in CartesianIndices(ks) band_structure = formula.calc_intensity(swt, ks[kidx]) diff --git a/src/SpinWaveTheory/DispersionAndIntensities.jl b/src/SpinWaveTheory/DispersionAndIntensities.jl index 3b99f8253..08c7e4f32 100644 --- a/src/SpinWaveTheory/DispersionAndIntensities.jl +++ b/src/SpinWaveTheory/DispersionAndIntensities.jl @@ -42,8 +42,7 @@ function bogoliubov!(V::Matrix{ComplexF64}, H::Matrix{ComplexF64}) # Inverse of λ gives eigenvalues of Ĩ H. We only care about the first L # eigenvalues, which are positive. A factor of 2 is needed to get the # physical quasiparticle energies. - #disp = resize!(λ, L) - disp = λ + disp = resize!(λ, L) @. disp = 2 / disp # In the special case that H(q) = H(-q) (i.e., a magnetic ordering with @@ -104,7 +103,7 @@ function dispersion(swt::SpinWaveTheory, qs) H = zeros(ComplexF64, 2nmodes, 2nmodes) V = zeros(ComplexF64, 2nmodes, 2nmodes) - disp = zeros(Float64, 2nmodes, length(qs)) + disp = zeros(Float64, nmodes, length(qs)) for (iq, q) in enumerate(qs) q_reshaped = to_reshaped_rlu(swt.sys, q) @@ -154,8 +153,8 @@ function dssf(swt::SpinWaveTheory, qs) qs = Vec3.(qs) nmodes = nbands(swt) - disp = zeros(Float64, 2nmodes, size(qs)...) - Sαβs = zeros(ComplexF64, 3, 3, 2nmodes, size(qs)...) + disp = zeros(Float64, nmodes, size(qs)...) + Sαβs = zeros(ComplexF64, 3, 3, nmodes, size(qs)...) # dssf(...) doesn't do any contraction, temperature correction, etc. # It simply returns the full Sαβ correlation matrix @@ -165,7 +164,7 @@ function dssf(swt::SpinWaveTheory, qs) for qidx in CartesianIndices(qs) q = qs[qidx] band_structure = formula.calc_intensity(swt,q) - for band = 1:(2nmodes) + for band = 1:nmodes disp[band,qidx] = band_structure.dispersion[band] Sαβs[:,:,band,qidx] .= band_structure.intensity[band] end @@ -248,7 +247,7 @@ function intensity_formula(f::Function,swt::SpinWaveTheory,corr_ix::AbstractVect H = zeros(ComplexF64, 2*nmodes, 2*nmodes) V = zeros(ComplexF64, 2*nmodes, 2*nmodes) Avec_pref = zeros(ComplexF64, Nm) - intensity = zeros(return_type, 2nmodes) + intensity = zeros(return_type, nmodes) # Expand formfactors for symmetry classes to formfactors for all atoms in # crystal @@ -364,7 +363,7 @@ function intensity_formula(f::Function,swt::SpinWaveTheory,corr_ix::AbstractVect # If there is no specified kernel, we are done: just return the # BandStructure - return BandStructure{2*nmodes,return_type}(disp, intensity) + return BandStructure{nmodes,return_type}(disp, intensity) else # Smooth kernel --> Intensity as a function of ω (or a list of ωs) return function(ω) diff --git a/src/Sunny.jl b/src/Sunny.jl index 568132a76..7305c216f 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -101,7 +101,7 @@ include("SpinWaveTheory/DispersionAndIntensities.jl") include("SpinWaveTheory/Lanczos.jl") include("SpinWaveTheory/Chebyshev.jl") include("SpinWaveTheory/KPM.jl") -export SpinWaveTheory, dispersion, dssf, kpm_dssf, kpm_intensities +export SpinWaveTheory, dispersion, dssf, kpm_dssf, kpm_intensities, delta_function_kernel include("SampledCorrelations/SampledCorrelations.jl") include("SampledCorrelations/CorrelationUtils.jl") diff --git a/test/test_lswt.jl b/test/test_lswt.jl index 1758e7933..d1bc692ac 100644 --- a/test/test_lswt.jl +++ b/test/test_lswt.jl @@ -81,8 +81,6 @@ end -""" - @testitem "Lanczos Bounds" begin using LinearAlgebra n=10 @@ -644,4 +642,4 @@ end is_ref = [0.042551643530614226 0.05061618774203612 0.06118972834822155 0.07541981401191966 0.09518339306332815 0.12371078181463013 0.1669136576105103 0.23644634156679015 0.3574143016120729 0.589319004913063 1.0795750388213625 2.0570916947398903 2.6569442073245604 1.6749367428248116 0.8709898932656854 0.4927038808971025 0.30849624153704924 0.20903618563778026 0.1502266820186879 0.11286975075833217 0.08777050885259896 0.07013929292833648 0.057300827672198455 0.047672208009393216 0.04027090049655778 0.0344619791000266 0.02982086993597283 0.02605519995256082 0.022958441606058772 0.020381422739929135 0.01821426297223774 0.016374601990583864 0.014799738769001156 0.013441266708211906 0.01226133976759786 0.01123002734145714 0.010323410070056573 0.009522188815697288 0.008810654796662341 0.008175917662895354 0.007607320304517221 0.007095990541228231 0.006634494316422953 0.006216564975067639 0.005836890143659928 0.00549094262866035 0.005174845247781923 0.0048852620341421574 0.004619310095626147 0.004374487768740016 0.00414861571474443; 0.14853118271457438 0.19360218524421224 0.2621797495216911 0.3732135013522173 0.5678610986270936 0.944407657666259 1.7450675448755357 3.2945578443571577 3.9947874291087895 2.418790606759786 1.2612861617795654 0.7201697967624341 0.45442367736983386 0.30970051154679945 0.22353509755861065 0.1685055354270234 0.1313752454162916 0.10520387618639207 0.0860938574888819 0.0717287252840356 0.060665219355115055 0.05196772774655005 0.045008912709399315 0.03935576182418519 0.03470177881993486 0.030825189142763776 0.027562374631854493 0.024790523018178495 0.02241601889071837 0.020366506674885078 0.018585357752204528 0.017027745220739847 0.015657814448345492 0.014446613655329999 0.013370560099471452 0.01241028925551664 0.011549781566933818 0.010775692876605122 0.010076836041215703 0.009443775967574439 0.008868510590520351 0.008344217576743833 0.007865051732026552 0.007425981842385972 0.007022658419610481 0.006651305841367626 0.006308633878300496 0.005991764727412878 0.005698172523177506 0.0054256329470820505 0.005172181054614204; 0.042551643530614205 0.0506161877420361 0.06118972834822153 0.07541981401191963 0.0951833930633281 0.12371078181463005 0.16691365761051016 0.23644634156678992 0.3574143016120724 0.5893190049130621 1.0795750388213603 2.0570916947398867 2.656944207324563 1.6749367428248163 0.8709898932656877 0.4927038808971036 0.30849624153704985 0.20903618563778062 0.15022668201868813 0.1128697507583323 0.08777050885259907 0.07013929292833657 0.05730082767219852 0.04767220800939327 0.04027090049655782 0.03446197910002663 0.02982086993597286 0.026055199952560844 0.02295844160605879 0.020381422739929156 0.018214262972237757 0.016374601990583878 0.014799738769001168 0.013441266708211913 0.01226133976759787 0.011230027341457147 0.010323410070056582 0.009522188815697295 0.008810654796662348 0.00817591766289536 0.007607320304517226 0.007095990541228234 0.006634494316422957 0.006216564975067644 0.005836890143659932 0.005490942628660353 0.005174845247781927 0.00488526203414216 0.00461931009562615 0.0043744877687400185 0.0041486157147444325] @test is ≈ is_ref -end \ No newline at end of file +end From 869eaf9e56c6698401d0b1ab4fc1b8652e4b193d Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 13 Dec 2023 10:01:17 -0500 Subject: [PATCH 05/11] Remove plotting2d duplicates --- examples/extra/plotting2d.jl | 124 ----------------------------- examples/extra/plotting_2d.jl | 145 ---------------------------------- 2 files changed, 269 deletions(-) delete mode 100644 examples/extra/plotting2d.jl delete mode 100644 examples/extra/plotting_2d.jl diff --git a/examples/extra/plotting2d.jl b/examples/extra/plotting2d.jl deleted file mode 100644 index 99a55102e..000000000 --- a/examples/extra/plotting2d.jl +++ /dev/null @@ -1,124 +0,0 @@ -using LinearAlgebra -import ColorSchemes, ColorTypes - - -################################################# -# Functions for plotting on triangular plaquettes -################################################# - - -function plaquette_idcs(dims::Tuple{Int,Int,Int}) - dx, dy, dz = dims - (dz != 1) && println("Warning: Ignoring lattice c vector.") - Triple = Tuple{Int,Int,Int} - idcs = Array{Tuple{Triple,Triple,Triple},4}(undef, 2, dx + 1, dy + 1, 1) - for j ∈ 1:(dy+1) - for i ∈ 1:(dx+1) - idcs[1, i, j, 1] = ( - (mod1(i , dx), mod1(j, dy), 1), - (mod1(i+1, dx), mod1(j+1, dy), 1), - (mod1(i , dx), mod1(j+1, dy), 1), - ) - idcs[2, i, j, 1] = ( - (mod1(i, dx), mod1(j, dy), 1), - (mod1(i+1, dx), mod1(j, dy), 1), - (mod1(i+1, dx), mod1(j+1, dy), 1), - ) - end - end - return idcs -end - -plaquette_idcs(x::Array{T,4}) where T = plaquette_idcs(size(x)[1:3]) - -function plaquette_map(f::Function, x::Array{T,4}) where T - - dims = size(x) - @assert dims[3] == dims[4] == 1 "Multiple sites and multiple layers are not supported." - dims = dims[1:3] - out = zeros(Float64, 2, dims...) - idcs_all = plaquette_idcs(x) - for i ∈ CartesianIndices(dims) - idcs = idcs_all[1, i] - out[1, i] = f(x[idcs[1]..., 1], x[idcs[2]..., 1], x[idcs[3]..., 1]) - idcs = idcs_all[2, i] - out[2, i] = f(x[idcs[1]..., 1], x[idcs[2]..., 1], x[idcs[3]..., 1]) - end - out -end - - -################################################################################ -# Plotting functions -################################################################################ -function aspect_ratio(x_panel, y_panel, x_offset, y_offset, numrows, numcols; - adhoc_offset = (0.0, 0.0) -) - corners = [ - (0, 0), - numrows * y_panel + (numrows - 1) * y_offset, - numcols * x_panel + (numcols - 1) * x_offset, - numcols * x_panel + numrows * y_panel + (numrows - 1) * y_offset + (numcols - 1) * x_offset, - ] - xs = [c[1] for c ∈ corners] - ys = [c[2] for c ∈ corners] - x1, x2 = minimum(xs), maximum(xs) + adhoc_offset[1] - y1, y2 = minimum(ys), maximum(ys) + adhoc_offset[2] - return abs((x2 - x1) / (y2 - y1)) -end - - -function plot_triangular_plaquettes(f, frames; - colorscheme=ColorSchemes.RdBu, clims=(-0.5, 0.5), offset_spacing=1, - numcols=nothing, texts=nothing, force_aspect=true, text_offset = (0.0, 0.0), fig_kwargs... -) - # Consolidate lattice info and panel layout - numpanels = length(frames) - isnothing(numcols) && (numcols = numpanels) - numrows = floor(Int, (numpanels - 1) / numcols) + 1 - v₁, v₂ = [1, 0, 0], [-1/2, √3/2, 0] # Derives from lattice_vectors(a,a,c,90,90,120) - nx, ny = size(frames[1])[1:2] - v₁, v₂ = Point3f(v₁), Point3f(v₂) - x, y = [1.0, 0, 0], [0.0, 1, 0] - x_offset = offset_spacing * (x ⋅ v₁) * x - y_offset = -offset_spacing * (y ⋅ v₂) * y - x_panel = nx * v₁ - y_panel = -ny * v₂ - aspect = aspect_ratio(x_panel, y_panel, x_offset, y_offset, numrows, numcols; adhoc_offset=text_offset) - - # Set up figure - fig = Figure(; fig_kwargs...) - if force_aspect - ax = Axis(fig[1, 1:length(frames)]; aspect) - else - ax = Axis(fig[1, 1:length(frames)]; aspect=true) - end - hidespines!(ax) - hidedecorations!(ax) - - # Plot panels - plaq1(p) = Makie.Polygon(Point2f.([p, p+v₁+v₂, p+v₂])) - plaq2(p) = Makie.Polygon(Point2f.([p, p+v₁, p+v₂+v₁])) - for (i, frame) ∈ enumerate(frames) - r, c = fldmod1(i, numcols) - v₀ = (c - 1) * (x_panel + x_offset) + (r - 1) * (y_panel + y_offset) - - χ = plaquette_map(f, frame) - pgons = Makie.Polygon[] - colors = ColorTypes.RGB{Float64}[] - for r ∈ 1:nx - for c ∈ 1:ny - base = (r - 1) * v₁ + (c - 1) * v₂ + v₀ - push!(pgons, plaq1(base)) - push!(colors, get(colorscheme, χ[1, r, c, 1, 1], clims)) - push!(pgons, plaq2(base)) - push!(colors, get(colorscheme, χ[2, r, c, 1, 1], clims)) - end - end - poly!(ax, pgons; color=colors) - if !isnothing(texts) - text!(ax, v₀[1] - text_offset[1], v₀[2] - text_offset[2]; text=texts[i], fontsize=36) - end - end - return fig -end \ No newline at end of file diff --git a/examples/extra/plotting_2d.jl b/examples/extra/plotting_2d.jl deleted file mode 100644 index e6d9ca966..000000000 --- a/examples/extra/plotting_2d.jl +++ /dev/null @@ -1,145 +0,0 @@ -using LinearAlgebra -import ColorSchemes, ColorTypes - - -################################################################################ -# Functions for plotting on triangular plaquettes, in particular Berry curvature -################################################################################ - -function vec3_triple_product(s₁, s₂, s₃) - s₁, s₂, s₃ = normalize.((s₁, s₂, s₃)) - return s₁ ⋅ (s₂ × s₃) -end - -function sun_berry_curvature(z₁, z₂, z₃) - z₁, z₂, z₃ = normalize.((z₁, z₂, z₃)) - n₁ = z₁ ⋅ z₂ - n₂ = z₂ ⋅ z₃ - n₃ = z₃ ⋅ z₁ - return angle(n₁ * n₂ * n₃) -end - - -function plaquette_idcs(dims::Tuple{Int,Int,Int}) - dx, dy, dz = dims - (dz != 1) && println("Warning: Ignoring lattice c vector.") - Triple = Tuple{Int,Int,Int} - idcs = Array{Tuple{Triple,Triple,Triple},4}(undef, 2, dx + 1, dy + 1, 1) - for j ∈ 1:(dy+1) - for i ∈ 1:(dx+1) - idcs[1, i, j, 1] = ( - (mod1(i , dx), mod1(j, dy), 1), - (mod1(i+1, dx), mod1(j+1, dy), 1), - (mod1(i , dx), mod1(j+1, dy), 1), - ) - idcs[2, i, j, 1] = ( - (mod1(i, dx), mod1(j, dy), 1), - (mod1(i+1, dx), mod1(j, dy), 1), - (mod1(i+1, dx), mod1(j+1, dy), 1), - ) - end - end - return idcs -end - -plaquette_idcs(x::Array{T,4}) where T = plaquette_idcs(size(x)[1:3]) - -function plaquette_map(f::Function, x::Array{T,4}) where T - - dims = size(x) - @assert dims[3] == dims[4] == 1 "Multiple sites and multiple layers are not supported." - dims = dims[1:3] - out = zeros(Float64, 2, dims...) - idcs_all = plaquette_idcs(x) - for i ∈ CartesianIndices(dims) - idcs = idcs_all[1, i] - out[1, i] = f(x[idcs[1]..., 1], x[idcs[2]..., 1], x[idcs[3]..., 1]) - idcs = idcs_all[2, i] - out[2, i] = f(x[idcs[1]..., 1], x[idcs[2]..., 1], x[idcs[3]..., 1]) - end - out -end - - -################################################################################ -# Plotting functions -################################################################################ -function aspect_ratio(x_panel, y_panel, x_offset, y_offset, numrows, numcols) - corners = [ - (0, 0), - numrows * y_panel + (numrows - 1) * y_offset, - numcols * x_panel + (numcols - 1) * x_offset, - numcols * x_panel + numrows * y_panel + (numrows - 1) * y_offset + (numcols - 1) * x_offset, - ] - xs = [c[1] for c ∈ corners] - ys = [c[2] for c ∈ corners] - x1, x2 = minimum(xs), maximum(xs) - y1, y2 = minimum(ys), maximum(ys) - return abs((x2 - x1) / (y2 - y1)) -end - - -function plot_chirality(xs, sys; - colorscheme=ColorSchemes.RdBu, clims=(-0.5, 0.5), offset_spacing=1, - numcols=nothing, texts=nothing, force_aspect=true, fig_kwargs... -) - # Consolidate lattice info and panel layout - numpanels = length(xs) - isnothing(numcols) && (numcols = numpanels) - numrows = floor(Int, (numpanels - 1) / numcols) + 1 - vecs = sys.crystal.latvecs - v₁, v₂ = vecs[:, 1], vecs[:, 2] - nx, ny = size(xs[1])[1:2] - v₁, v₂ = Point3f(v₁), Point3f(v₂) - x, y = [1.0, 0, 0], [0.0, 1, 0] - x_offset = offset_spacing * (x ⋅ v₁) * x - y_offset = -offset_spacing * (y ⋅ v₂) * y - x_panel = nx * v₁ - y_panel = -ny * v₂ - aspect = aspect_ratio(x_panel, y_panel, x_offset, y_offset, numrows, numcols) - - # Set up figure - fig = Figure(; fig_kwargs...) - if force_aspect - ax = Axis(fig[1, 1:length(xs)]; aspect) - else - ax = Axis(fig[1, 1:length(xs)]) - end - hidespines!(ax) - hidedecorations!(ax) - - # Plot panels - plaq1(p) = GLMakie.Polygon(Point2f.([p, p+v₁+v₂, p+v₂])) - plaq2(p) = GLMakie.Polygon(Point2f.([p, p+v₁, p+v₂+v₁])) - for (i, x) ∈ enumerate(xs) - r, c = fldmod1(i, numcols) - v₀ = (c - 1) * (x_panel + x_offset) + (r - 1) * (y_panel + y_offset) - - f = (eltype(x) == Sunny.Vec3) ? vec3_triple_product : sun_berry_curvature - χ = plaquette_map(f, x) - pgons = GLMakie.Polygon[] - colors = ColorTypes.RGB{Float64}[] - for r ∈ 1:nx - for c ∈ 1:ny - base = (r - 1) * v₁ + (c - 1) * v₂ + v₀ - push!(pgons, plaq1(base)) - push!(colors, get(colorscheme, χ[1, r, c, 1, 1], clims)) - push!(pgons, plaq2(base)) - push!(colors, get(colorscheme, χ[2, r, c, 1, 1], clims)) - end - end - poly!(ax, pgons; color=colors) - if !isnothing(texts) - text!(ax, v₀[1], v₀[2]; text=texts[i]) - end - end - return fig -end - -function plot_chirality(sys; - colorscheme=ColorSchemes.RdBu, clims=(-0.5, 0.5), offset=1, - numcols=nothing, texts=nothing, fig_kwargs... -) - x = (sys.mode == :dipole) ? sys.dipoles : sys.coherents - return plot_chirality([x], sys; colorscheme, clims, offset, numcols, texts, fig_kwargs...) -end \ No newline at end of file From e69b660525ce8311ca2fb06540a785ebf8669700 Mon Sep 17 00:00:00 2001 From: hlane33 <119619760+hlane33@users.noreply.github.com> Date: Wed, 13 Dec 2023 10:01:17 -0500 Subject: [PATCH 06/11] Implement matrix free KPM. --- src/SpinWaveTheory/KPM.jl | 37 ++++++++++------------------ src/SpinWaveTheory/Lanczos.jl | 45 +++++++++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 24 deletions(-) diff --git a/src/SpinWaveTheory/KPM.jl b/src/SpinWaveTheory/KPM.jl index da0399bc4..f599dfd61 100644 --- a/src/SpinWaveTheory/KPM.jl +++ b/src/SpinWaveTheory/KPM.jl @@ -134,7 +134,6 @@ function kpm_dssf(swt::SpinWaveTheory, qs,ωlist,P::Int64,kT,σ,broadening; kern sqrt_halfS = √(S/2) #define prefactors Ĩ = spdiagm([ones(nmodes); -ones(nmodes)]) n_iters = 50 - Hmat = zeros(ComplexF64, 2*nmodes, 2*nmodes) Avec_pref = zeros(ComplexF64, Nm) # initialize array of some prefactors chebyshev_moments = OffsetArray(zeros(ComplexF64,3,3,length(qs),P),1:3,1:3,1:length(qs),0:P-1) Sαβs = zeros(ComplexF64,3,3,length(qs),length(ωlist)) @@ -142,15 +141,9 @@ function kpm_dssf(swt::SpinWaveTheory, qs,ωlist,P::Int64,kT,σ,broadening; kern q = qs[qidx] q_reshaped = Sunny.to_reshaped_rlu(swt.sys, q) u = zeros(ComplexF64,3,2*nmodes) - if sys.mode == :SUN - swt_hamiltonian_SUN!(Hmat, swt, q_reshaped) - else - swt_hamiltonian_dipole!(Hmat, swt, q_reshaped) - end - D = 2.0*sparse(Hmat) # calculate D (factor of 2 for correspondence) - lo,hi = Sunny.eigbounds(Ĩ*D,n_iters; extend=0.25) # calculate bounds + multiply_by_hamiltonian!(y, x) = (swt.sys.mode == :SUN) ? multiply_by_hamiltonian_SUN!(y, x, swt, q_reshaped) : multiply_by_hamiltonian_dipole!(y, x, swt, q_reshaped) + lo,hi = Sunny.eigbounds_MF(swt,q_reshaped,n_iters; extend=0.25) # calculate bounds γ=max(abs(lo),abs(hi)) # select upper bound (combine with the preceeding line later) - A = Ĩ*D / γ # u(q) calculation) for site = 1:Nm # note that d is the chemical coordinates @@ -185,14 +178,16 @@ function kpm_dssf(swt::SpinWaveTheory, qs,ωlist,P::Int64,kT,σ,broadening; kern α0 = zeros(ComplexF64,2*nmodes) α1 = zeros(ComplexF64,2*nmodes) mul!(α0,Ĩ,u[β,:]) # calculate α0 - mul!(α1,A,α0) # calculate α1 + multiply_by_hamiltonian!(α1,α0) + mul!(α1,Ĩ,2α1/γ) for α=1:3 chebyshev_moments[α,β,qidx,0] = (dot(u[α,:],α0)) #removed symmetrization chebyshev_moments[α,β,qidx,1] = (dot(u[α,:],α1)) #removed symmetrization end for m=2:P-1 αnew = zeros(ComplexF64,2*nmodes) - mul!(αnew,A,α1) + multiply_by_hamiltonian!(αnew,a1) + mul!(αnew,Ĩ,2αnew/γ) @. αnew = 2*αnew - α0 for α=1:3 chebyshev_moments[α,β,qidx,m] = (dot(u[α,:],αnew)) #removed symmetrization @@ -282,7 +277,6 @@ function intensity_formula_kpm(f,swt::SpinWaveTheory,corr_ix::AbstractVector{Int sqrt_halfS = √(S/2) #define prefactors Ĩ = spdiagm([ones(nmodes); -ones(nmodes)]) n_iters = 50 - Hmat = zeros(ComplexF64, 2*nmodes, 2*nmodes) Avec_pref = zeros(ComplexF64, Nm) # initialize array of some prefactors chebyshev_moments = OffsetArray(zeros(ComplexF64,3,3,P),1:3,1:3,0:P-1) @@ -290,16 +284,9 @@ function intensity_formula_kpm(f,swt::SpinWaveTheory,corr_ix::AbstractVector{Int q_reshaped = Sunny.to_reshaped_rlu(sys, q) q_absolute = swt.sys.crystal.recipvecs * q_reshaped u = zeros(ComplexF64,3,2*nmodes) - if sys.mode == :SUN - swt_hamiltonian_SUN!(Hmat, swt, q_reshaped) - else - swt_hamiltonian_dipole!(Hmat, swt, q_reshaped) - end - D = 2.0*sparse(Hmat) # calculate D (factor of 2 for correspondence) - lo,hi = Sunny.eigbounds(Ĩ*D,n_iters; extend=0.25) # calculate bounds - - γ=max(lo,hi) # select upper bound (combine with the preceeding line later) - A = Ĩ*D / γ + multiply_by_hamiltonian!(y, x) = (swt.sys.mode == :SUN) ? multiply_by_hamiltonian_SUN!(y, x, swt, q_reshaped) : multiply_by_hamiltonian_dipole!(y, x, swt, q_reshaped) + lo,hi = Sunny.eigbounds_MF(swt,q_reshaped,n_iters; extend=0.25) # calculate bounds + γ=max(abs(lo),abs(hi)) # select upper bound (combine with the preceeding line later) # u(q) calculation for site = 1:Nm # note that d is the chemical coordinates @@ -334,14 +321,16 @@ function intensity_formula_kpm(f,swt::SpinWaveTheory,corr_ix::AbstractVector{Int α0 = zeros(ComplexF64,2*nmodes) α1 = zeros(ComplexF64,2*nmodes) mul!(α0,Ĩ,u[β,:]) # calculate α0 - mul!(α1,A,α0) # calculate α1 + multiply_by_hamiltonian!(α1,α0) + mul!(α1,Ĩ,2α1/γ) for α=1:3 chebyshev_moments[α,β,0] = (dot(u[α,:],α0)) #removed symmetrization chebyshev_moments[α,β,1] = (dot(u[α,:],α1)) #removed symmetrization end for m=2:P-1 αnew = zeros(ComplexF64,2*nmodes) - mul!(αnew,A,α1) + multiply_by_hamiltonian!(αnew,α1) + mul!(αnew,Ĩ,2αnew/γ) @. αnew = 2*αnew - α0 for α=1:3 chebyshev_moments[α,β,m] = (dot(u[α,:],αnew)) #removed symmetrization diff --git a/src/SpinWaveTheory/Lanczos.jl b/src/SpinWaveTheory/Lanczos.jl index 8ee0088a8..630bec8fa 100644 --- a/src/SpinWaveTheory/Lanczos.jl +++ b/src/SpinWaveTheory/Lanczos.jl @@ -81,6 +81,51 @@ function lanczos_QH_aux!(αs, βs, buf, A, niters) end +# Matrix free implementation of QH Lanczos on Hamiltonian +function eigbounds_MF(swt,q_reshaped, niters; extend=0.5) + lo, hi = lanczos_MF(swt,q_reshaped, niters) |> eigvals |> xs -> (first(xs), last(xs)) + slack = extend*(hi-lo) + return lo-slack, hi+slack +end + + +function lanczos_MF(swt,q_reshaped, niters) + nmodes=Sunny.num_bands(swt) + αs = zeros(Float64, niters) # Main diagonal + βs = zeros(Float64, niters-1) # Off diagonal + buf = zeros(ComplexF64, 2nmodes, 3) # Vector buffer -- don't technically need 3 columns, but more legible this way + lanczos_QH_MF_aux!(αs, βs, buf, swt ,q_reshaped, niters) + return SymTridiagonal(αs, βs) +end + +function lanczos_QH_MF_aux!(αs, βs, buf, swt ,q_reshaped, niters) + multiply_by_hamiltonian!(y, x) = (swt.sys.mode == :SUN) ? multiply_by_hamiltonian_SUN!(y,x, swt, q_reshaped) : multiply_by_hamiltonian_dipole!(y, x, swt, q_reshaped) + v, vprev, w = view(buf,:,1), view(buf,:,2), view(buf,:,3) + randn!(vprev) + mat_size=num_bands(swt) + Ĩ = diagm([ones(mat_size); -ones(mat_size)]) + multiply_by_hamiltonian!(w,vprev) + w = mul!(w,Ĩ,1w) + b0 = sqrt(Ĩ_dot_3(vprev, w,mat_size)) + vprev = vprev/b0 + w = w/b0 + αs[1] = Ĩ_dot_3(w, w,mat_size) + @. w = w - αs[1]*vprev + for j in 2:niters + @. v = w + multiply_by_hamiltonian!(w,v) + mul!(w,Ĩ,1w) + βs[j-1] = sqrt(Ĩ_dot_3(v, w,mat_size)) # maybe v? + iszero(βs[j-1]) && return + @. v = v/βs[j-1] + @. w = w/βs[j-1] + αs[j] = Ĩ_dot_3(w, w,mat_size) + @. w = w - αs[j]*v - βs[j-1]*vprev + v, vprev = vprev, v + end + return nothing +end + function lanczos_legacy(A, niters) From 167f57c15bca276db4af347440837d8662d913df Mon Sep 17 00:00:00 2001 From: hlane33 <119619760+hlane33@users.noreply.github.com> Date: Wed, 13 Dec 2023 10:01:17 -0500 Subject: [PATCH 07/11] Extend Lanczos buffer. --- src/SpinWaveTheory/KPM.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/SpinWaveTheory/KPM.jl b/src/SpinWaveTheory/KPM.jl index f599dfd61..f804a746e 100644 --- a/src/SpinWaveTheory/KPM.jl +++ b/src/SpinWaveTheory/KPM.jl @@ -272,9 +272,6 @@ function intensity_formula_kpm(f,swt::SpinWaveTheory,corr_ix::AbstractVector{Int sqrt_Nm_inv = 1.0 / √Nm #define prefactors S = (Ns-1) / 2 sqrt_halfS = √(S/2) #define prefactors - sqrt_Nm_inv = 1.0 / √Nm #define prefactors - S = (Ns-1) / 2 - sqrt_halfS = √(S/2) #define prefactors Ĩ = spdiagm([ones(nmodes); -ones(nmodes)]) n_iters = 50 Avec_pref = zeros(ComplexF64, Nm) # initialize array of some prefactors @@ -285,7 +282,7 @@ function intensity_formula_kpm(f,swt::SpinWaveTheory,corr_ix::AbstractVector{Int q_absolute = swt.sys.crystal.recipvecs * q_reshaped u = zeros(ComplexF64,3,2*nmodes) multiply_by_hamiltonian!(y, x) = (swt.sys.mode == :SUN) ? multiply_by_hamiltonian_SUN!(y, x, swt, q_reshaped) : multiply_by_hamiltonian_dipole!(y, x, swt, q_reshaped) - lo,hi = Sunny.eigbounds_MF(swt,q_reshaped,n_iters; extend=0.25) # calculate bounds + lo,hi = Sunny.eigbounds_MF(swt,q_reshaped,n_iters; extend=0.5) # calculate bounds γ=max(abs(lo),abs(hi)) # select upper bound (combine with the preceeding line later) # u(q) calculation for site = 1:Nm From d8eafdecb293fc26958a1476ba44d44ef7089ae9 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Thu, 14 Dec 2023 08:53:04 -0700 Subject: [PATCH 08/11] Alloc-free Chebyshev recursion --- src/SpinWaveTheory/KPM.jl | 42 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/src/SpinWaveTheory/KPM.jl b/src/SpinWaveTheory/KPM.jl index f804a746e..065ac3c4b 100644 --- a/src/SpinWaveTheory/KPM.jl +++ b/src/SpinWaveTheory/KPM.jl @@ -360,3 +360,45 @@ function intensity_formula_kpm(f,swt::SpinWaveTheory,corr_ix::AbstractVector{Int end KPMIntensityFormula{return_type}(P,kT,σ,broadening,kernel,string_formula,calc_intensity) end + +#= +function Itilde!(α, n) + view(α, n+1:2n) .*= -1 +end + +function Run_Recurrence_fast(swt,q_reshaped,γ,u,nmodes,chebyshev_moments,M) + α0 = zeros(ComplexF64,2*nmodes) + α1 = zeros(ComplexF64,2*nmodes) + α2 = zeros(ComplexF64,2*nmodes) + + chebyshev_moments .= 0 + for β=1:3 + # α0 = ̃I u + α0 .= view(u, β, :) + Itilde!(α0, nmodes) + + # α1 = ̃I A α0 + Sunny.multiply_by_hamiltonian_dipole!(α1,α0,swt,q_reshaped) + @. α1 = α1 * (2/γ) + Itilde!(α1, nmodes) + + for α=1:3 + chebyshev_moments[α,β,0] = dot(view(u, α,:), α0) + chebyshev_moments[α,β,1] = dot(view(u, α,:), α1) + end + + for m=2:M-1 + # α2 = ̃2 I A α1 - α0 + Sunny.multiply_by_hamiltonian_dipole!(α2,α1,swt,q_reshaped) + @. α2 = α2 * (2/γ) + Itilde!(α2, nmodes) + @. α2 = 2*α2 - α0 + + for α=1:3 + chebyshev_moments[α,β,m] = dot(view(u, α,:),α2) + end + (α1, α0, α2) = (α2, α1, α0) + end + end +end +=# \ No newline at end of file From c79bf1c158552563f998ed6cb91d9224ec778ed0 Mon Sep 17 00:00:00 2001 From: hlane33 <119619760+hlane33@users.noreply.github.com> Date: Wed, 7 Feb 2024 09:53:47 +0000 Subject: [PATCH 09/11] Tidy up some code. --- src/SpinWaveTheory/KPM.jl | 2 +- src/SpinWaveTheory/Lanczos.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/SpinWaveTheory/KPM.jl b/src/SpinWaveTheory/KPM.jl index 065ac3c4b..4a860bf89 100644 --- a/src/SpinWaveTheory/KPM.jl +++ b/src/SpinWaveTheory/KPM.jl @@ -186,7 +186,7 @@ function kpm_dssf(swt::SpinWaveTheory, qs,ωlist,P::Int64,kT,σ,broadening; kern end for m=2:P-1 αnew = zeros(ComplexF64,2*nmodes) - multiply_by_hamiltonian!(αnew,a1) + multiply_by_hamiltonian!(αnew,α1) mul!(αnew,Ĩ,2αnew/γ) @. αnew = 2*αnew - α0 for α=1:3 diff --git a/src/SpinWaveTheory/Lanczos.jl b/src/SpinWaveTheory/Lanczos.jl index 630bec8fa..ccd6ac41a 100644 --- a/src/SpinWaveTheory/Lanczos.jl +++ b/src/SpinWaveTheory/Lanczos.jl @@ -105,7 +105,7 @@ function lanczos_QH_MF_aux!(αs, βs, buf, swt ,q_reshaped, niters) mat_size=num_bands(swt) Ĩ = diagm([ones(mat_size); -ones(mat_size)]) multiply_by_hamiltonian!(w,vprev) - w = mul!(w,Ĩ,1w) + mul!(w,Ĩ,1w) b0 = sqrt(Ĩ_dot_3(vprev, w,mat_size)) vprev = vprev/b0 w = w/b0 From b7b2eaa8017ed1189ebd3f6589bca7920785d907 Mon Sep 17 00:00:00 2001 From: hlane33 <119619760+hlane33@users.noreply.github.com> Date: Wed, 7 Feb 2024 15:44:02 +0000 Subject: [PATCH 10/11] Add in-place functions as stopgap. --- src/SpinWaveTheory/HamiltonianDipole.jl | 106 +++++++++++++++++++++++- src/SpinWaveTheory/HamiltonianSUN.jl | 40 +++++++++ src/SpinWaveTheory/Lanczos.jl | 4 +- 3 files changed, 147 insertions(+), 3 deletions(-) diff --git a/src/SpinWaveTheory/HamiltonianDipole.jl b/src/SpinWaveTheory/HamiltonianDipole.jl index e0f1c59d6..83b3c44f0 100644 --- a/src/SpinWaveTheory/HamiltonianDipole.jl +++ b/src/SpinWaveTheory/HamiltonianDipole.jl @@ -242,4 +242,108 @@ function multiply_by_hamiltonian_dipole_aux!(y, x, phasebuf, qphase, swt) @inbounds @. y += swt.energy_ϵ * x nothing -end \ No newline at end of file +end + +# This is a stopgap measure to avoid unnecessary reshaping. Revisit on merge. + +# Calculate y = H*x, where H is the quadratic Hamiltonian matrix (dynamical matrix). +function multiply_by_hamiltonian_dipole!(y, x, swt::SpinWaveTheory, q_reshaped::Vec3) + (; sys, data) = swt + (; local_rotations, stevens_coefs) = data + + N = swt.sys.Ns[1] + S = (N-1)/2 + L = nbands(swt) + + x = reshape(x, natoms(sys.crystal), 2) + y = reshape(y, natoms(sys.crystal), 2) + y .= 0 + + # Add Zeeman term + (; extfield, gs, units) = sys + for i in 1:L + B = units.μB * (gs[1, 1, 1, i]' * extfield[1, 1, 1, i]) + B′ = dot(B, local_rotations[i][:, 3]) / 2 + + y[i, 1] += B′ * x[i, 1] + y[i, 2] += B′ * x[i, 2] + end + + # Add pairwise terms + for ints in sys.interactions_union + + # Bilinear exchange + for coupling in ints.pair + (; isculled, bond) = coupling + isculled && break + i, j = bond.i, bond.j + phase = exp(2π*im * dot(q_reshaped, bond.n)) # Phase associated with periodic wrapping + + if !iszero(coupling.bilin) + J = coupling.bilin # This is Rij in previous notation (transformed exchange matrix) + + P = 0.25 * (J[1, 1] - J[2, 2] - im*J[1, 2] - im*J[2, 1]) + Q = 0.25 * (J[1, 1] + J[2, 2] - im*J[1, 2] + im*J[2, 1]) + + y[i, 1] += Q * phase * x[j, 1] + y[j, 1] += conj(Q) * conj(phase) * x[i, 1] + y[i, 2] += conj(Q) * phase * x[j, 2] + y[j, 2] += Q * conj(phase) * x[i, 2] + + y[i, 2] += P * phase * x[j, 1] + y[j, 2] += P * conj(phase) * x[i, 1] + y[i, 1] += conj(P) * phase * x[j, 2] + y[j, 1] += conj(P) * conj(phase) * x[i, 2] + + y[i, 1] -= 0.5 * J[3, 3] * x[i, 1] + y[j, 1] -= 0.5 * J[3, 3] * x[j, 1] + y[i, 2] -= 0.5 * J[3, 3] * x[i, 2] + y[j, 2] -= 0.5 * J[3, 3] * x[j, 2] + end + + # Biquadratic exchange + if !iszero(coupling.biquad) + J = coupling.biquad # Transformed quadrupole exchange matrix + + y[i, 1] += -6J[3, 3] * x[i, 1] + y[j, 1] += -6J[3, 3] * x[j, 1] + + y[i, 2] += -6J[3, 3] * x[i, 2] + y[j, 2] += -6J[3, 3] * x[j, 2] + + y[i, 2] += 12*(J[1, 3] - im*J[5, 3]) * x[i, 1] + y[i, 1] += 12*(J[1, 3] + im*J[5, 3]) * x[i, 2] + y[j, 2] += 12*(J[3, 1] - im*J[3, 5]) * x[j, 1] + y[j, 1] += 12*(J[3, 1] + im*J[3, 5]) * x[j, 2] + + P = 0.25 * (-J[4, 4]+J[2, 2] - im*( J[4, 2]+J[2, 4])) + Q = 0.25 * ( J[4, 4]+J[2, 2] - im*(-J[4, 2]+J[2, 4])) + + y[i, 1] += Q * phase * x[j, 1] + y[j, 1] += conj(Q) * conj(phase) * x[i, 1] + y[i, 2] += conj(Q) * phase * x[j, 2] + y[j, 2] += Q * conj(phase) * x[i, 2] + + y[i, 2] += P * phase * x[j, 1] + y[j, 2] += P * conj(phase) * x[i, 1] + y[i, 1] += conj(P) * phase * x[j, 2] + y[j, 1] += conj(P) * conj(phase) * x[i, 2] + end + end + end + + # Add single-ion anisotropy + for i in 1:L + (; c2, c4, c6) = stevens_coefs[i] + y[i, 1] += (-3S*c2[3] - 40*S^3*c4[5] - 168*S^5*c6[7]) * x[i, 1] + y[i, 2] += (-3S*c2[3] - 40*S^3*c4[5] - 168*S^5*c6[7]) * x[i, 2] + y[i, 1] += (-im*(S*c2[5] + 6S^3*c4[7] + 16S^5*c6[9]) + (S*c2[1] + 6S^3*c4[3] + 16S^5*c6[5])) * x[i, 2] + y[i, 2] += (im*(S*c2[5] + 6S^3*c4[7] + 16S^5*c6[9]) + (S*c2[1] + 6S^3*c4[3] + 16S^5*c6[5])) * x[i, 1] + end + + # Add small constant shift for positive-definiteness + @. y += swt.energy_ϵ * x + + nothing +end + diff --git a/src/SpinWaveTheory/HamiltonianSUN.jl b/src/SpinWaveTheory/HamiltonianSUN.jl index cc7ab72e9..53f12a5a5 100644 --- a/src/SpinWaveTheory/HamiltonianSUN.jl +++ b/src/SpinWaveTheory/HamiltonianSUN.jl @@ -234,5 +234,45 @@ function multiply_by_hamiltonian_SUN_aux!(y, x, phasebuf, qphase, swt) # Add small constant shift for positive-definiteness @inbounds @. y += swt.energy_ϵ * x + nothing +end + +# This is a stopgap measure to avoid unnecessary reshaping. Revisit on merge. + +# Calculate y = H*x, where H is the quadratic Hamiltonian matrix (dynamical matrix). +function multiply_by_hamiltonian_SUN!(y, x, swt, q_reshaped) + (; sys, data) = swt + (; zeeman_operators) = data + y .= 0 + + # Add single-site terms (single-site anisotropy and external field) + # Couple percent speedup if this is removed and accumulated into onsite term + # (not pursuing for now to maintain parallelism with dipole mode). + for atom in 1:natoms(sys.crystal) + zeeman = view(zeeman_operators, :, :, atom) + multiply_by_onsite_coupling_SUN!(y, x, zeeman, swt, atom) + end + + # Add pair interactions that use explicit bases + for (atom, int) in enumerate(sys.interactions_union) + + # Set the onsite term + multiply_by_onsite_coupling_SUN!(y, x, int.onsite, swt, atom) + + for coupling in int.pair + # Extract information common to bond + (; isculled, bond) = coupling + isculled && break + + phase = exp(2π*im * dot(q_reshaped, bond.n)) # Phase associated with periodic wrapping + for (A, B) in coupling.general.data + multiply_by_pair_coupling_SUN!(y, x, A, B, swt, phase, bond) + end + end + end + + # # Add small constant shift for positive-definiteness + @. y += swt.energy_ϵ * x + nothing end \ No newline at end of file diff --git a/src/SpinWaveTheory/Lanczos.jl b/src/SpinWaveTheory/Lanczos.jl index ccd6ac41a..f6f071a64 100644 --- a/src/SpinWaveTheory/Lanczos.jl +++ b/src/SpinWaveTheory/Lanczos.jl @@ -90,7 +90,7 @@ end function lanczos_MF(swt,q_reshaped, niters) - nmodes=Sunny.num_bands(swt) + nmodes=Sunny.nbands(swt) αs = zeros(Float64, niters) # Main diagonal βs = zeros(Float64, niters-1) # Off diagonal buf = zeros(ComplexF64, 2nmodes, 3) # Vector buffer -- don't technically need 3 columns, but more legible this way @@ -102,7 +102,7 @@ function lanczos_QH_MF_aux!(αs, βs, buf, swt ,q_reshaped, niters) multiply_by_hamiltonian!(y, x) = (swt.sys.mode == :SUN) ? multiply_by_hamiltonian_SUN!(y,x, swt, q_reshaped) : multiply_by_hamiltonian_dipole!(y, x, swt, q_reshaped) v, vprev, w = view(buf,:,1), view(buf,:,2), view(buf,:,3) randn!(vprev) - mat_size=num_bands(swt) + mat_size=nbands(swt) Ĩ = diagm([ones(mat_size); -ones(mat_size)]) multiply_by_hamiltonian!(w,vprev) mul!(w,Ĩ,1w) From 6b0682c2ed40662eefd46c3dd5d2bf0576580190 Mon Sep 17 00:00:00 2001 From: hlane33 <119619760+hlane33@users.noreply.github.com> Date: Wed, 7 Feb 2024 17:34:09 +0000 Subject: [PATCH 11/11] Fix updated :SUN variable names. --- src/SpinWaveTheory/HamiltonianSUN.jl | 63 ++++++++++++++++++++++++++-- src/SpinWaveTheory/KPM.jl | 4 +- 2 files changed, 62 insertions(+), 5 deletions(-) diff --git a/src/SpinWaveTheory/HamiltonianSUN.jl b/src/SpinWaveTheory/HamiltonianSUN.jl index 53f12a5a5..9e8b747f2 100644 --- a/src/SpinWaveTheory/HamiltonianSUN.jl +++ b/src/SpinWaveTheory/HamiltonianSUN.jl @@ -250,14 +250,14 @@ function multiply_by_hamiltonian_SUN!(y, x, swt, q_reshaped) # (not pursuing for now to maintain parallelism with dipole mode). for atom in 1:natoms(sys.crystal) zeeman = view(zeeman_operators, :, :, atom) - multiply_by_onsite_coupling_SUN!(y, x, zeeman, swt, atom) + multiply_by_onsite_coupling_SUN_legacy!(y, x, zeeman, swt, atom) end # Add pair interactions that use explicit bases for (atom, int) in enumerate(sys.interactions_union) # Set the onsite term - multiply_by_onsite_coupling_SUN!(y, x, int.onsite, swt, atom) + multiply_by_onsite_coupling_SUN_legacy!(y, x, int.onsite, swt, atom) for coupling in int.pair # Extract information common to bond @@ -266,7 +266,7 @@ function multiply_by_hamiltonian_SUN!(y, x, swt, q_reshaped) phase = exp(2π*im * dot(q_reshaped, bond.n)) # Phase associated with periodic wrapping for (A, B) in coupling.general.data - multiply_by_pair_coupling_SUN!(y, x, A, B, swt, phase, bond) + multiply_by_pair_coupling_SUN_legacy!(y, x, A, B, swt, phase, bond) end end end @@ -275,4 +275,61 @@ function multiply_by_hamiltonian_SUN!(y, x, swt, q_reshaped) @. y += swt.energy_ϵ * x nothing +end + +# Calculate y = H_{onsite}*x, where H_{onsite} is the portion of the quadratic +# Hamiltonian matrix (dynamical matrix) due to onsite terms (other than Zeeman). +function multiply_by_onsite_coupling_SUN_legacy!(y, x, op, swt, atom) + sys = swt.sys + N = sys.Ns[1] + nflavors = N - 1 + + x = reshape(x, nflavors, natoms(sys.crystal), 2) + y = reshape(y, nflavors, natoms(sys.crystal), 2) + + for m in 1:N-1 + for n in 1:N-1 + c = 0.5 * (op[m, n] - δ(m, n) * op[N, N]) + y[m, atom, 1] += c * x[n, atom, 1] + y[n, atom, 2] += c * x[m, atom, 2] + end + end +end + +# Calculate y = H_{pair}*x, where H_{pair} is the portion of the quadratic +# Hamiltonian matrix (dynamical matrix) due to pair-wise interactions. +function multiply_by_pair_coupling_SUN_legacy!(y, x, Ti, Tj, swt, phase, bond) + (; i, j) = bond + sys = swt.sys + N = sys.Ns[1] + nflavors = N - 1 + + x = reshape(x, nflavors, natoms(sys.crystal), 2) + y = reshape(y, nflavors, natoms(sys.crystal), 2) + + for m in 1:N-1 + for n in 1:N-1 + c = 0.5 * (Ti[m,n] - δ(m,n)*Ti[N,N]) * Tj[N,N] + y[m, i, 1] += c * x[n, i, 1] + y[n, i, 2] += c * x[m, i, 2] + + c = 0.5 * Ti[N,N] * (Tj[m,n] - δ(m,n)*Tj[N,N]) + y[m, j, 1] += c * x[n, j, 1] + y[n, j, 2] += c * x[m, j, 2] + + c = 0.5 * Ti[m,N] * Tj[N,n] + y[m, i, 1] += c * phase * x[n, j, 1] + y[n, j, 2] += c * conj(phase) * x[m, i, 2] + + c = 0.5 * Ti[N,m] * Tj[n,N] + y[n, j, 1] += c * conj(phase) * x[m, i, 1] + y[m, i, 2] += c * phase * x[n, j, 2] + + c = 0.5 * Ti[m,N] * Tj[n,N] + y[m, i, 1] += c * phase * x[n, j, 2] + y[n, j, 1] += c * conj(phase) * x[m, i, 2] + y[m, i, 2] += conj(c * phase) * x[n, j, 1] + y[n, j, 2] += conj(c) *phase * x[m, i, 1] + end + end end \ No newline at end of file diff --git a/src/SpinWaveTheory/KPM.jl b/src/SpinWaveTheory/KPM.jl index 4a860bf89..47d269b5f 100644 --- a/src/SpinWaveTheory/KPM.jl +++ b/src/SpinWaveTheory/KPM.jl @@ -154,7 +154,7 @@ function kpm_dssf(swt::SpinWaveTheory, qs,ωlist,P::Int64,kT,σ,broadening; kern # calculate u(q) if sys.mode == :SUN for site=1:Nm - @views tS_μ = data.dipole_operators[:, :, :, site]*Avec_pref[site] + @views tS_μ = data.observable_operators[:, :, :, site]*Avec_pref[site] for μ=1:3 for j=2:N u[μ,(j-1)+(site-1)*(N-1) ]=tS_μ[j,1,μ] @@ -294,7 +294,7 @@ function intensity_formula_kpm(f,swt::SpinWaveTheory,corr_ix::AbstractVector{Int # calculate u(q) if sys.mode == :SUN for site=1:Nm - @views tS_μ = data.dipole_operators[:, :, :, site]*Avec_pref[site] + @views tS_μ = data.observable_operators[:, :, :, site]*Avec_pref[site] for μ=1:3 for j=2:N u[μ,(j-1)+(site-1)*(N-1) ]=tS_μ[j,1,μ]