diff --git a/Project.toml b/Project.toml index 1c8f91d23..c3b3569d1 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" @@ -44,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/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/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/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..08c7e4f32 100644 --- a/src/SpinWaveTheory/DispersionAndIntensities.jl +++ b/src/SpinWaveTheory/DispersionAndIntensities.jl @@ -375,4 +375,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/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 4e3b476bd..9e8b747f2 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 @@ -235,4 +235,101 @@ function multiply_by_hamiltonian_SUN_aux!(y, x, phasebuf, qphase, swt) @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_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_legacy!(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_legacy!(y, x, A, B, swt, phase, bond) + end + end + end + + # # Add small constant shift for positive-definiteness + @. 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 new file mode 100644 index 000000000..47d269b5f --- /dev/null +++ b/src/SpinWaveTheory/KPM.jl @@ -0,0 +1,404 @@ +""" +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 + 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) + 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 + 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.observable_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 + 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) + multiply_by_hamiltonian!(αnew,α1) + mul!(αnew,Ĩ,2αnew/γ) + @. α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 + Ĩ = spdiagm([ones(nmodes); -ones(nmodes)]) + n_iters = 50 + 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) + 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.5) # 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 + 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.observable_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 + 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) + multiply_by_hamiltonian!(αnew,α1) + mul!(αnew,Ĩ,2αnew/γ) + @. α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 + +#= +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 diff --git a/src/SpinWaveTheory/Lanczos.jl b/src/SpinWaveTheory/Lanczos.jl index 8d9fa083d..f6f071a64 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,126 @@ 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 + + +# 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.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 + 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=nbands(swt) + Ĩ = diagm([ones(mat_size); -ones(mat_size)]) + multiply_by_hamiltonian!(w,vprev) + 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) + @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..7305c216f 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -19,6 +19,9 @@ import CrystalInfoFramework as CIF import Spglib import RowEchelon: rref! +# 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 +99,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, 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 b0e9ea23a..d1bc692ac 100644 --- a/test/test_lswt.jl +++ b/test/test_lswt.jl @@ -82,17 +82,71 @@ 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) - + 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 @@ -588,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