Stabilize KPM

docs/src/versions.md

@@ -5,6 +5,8 @@
 
 * Fix error in `print_symmetry_table` for slightly-distorted crystal cells ([PR
   #317](https://github.com/SunnySuite/Sunny.jl/pull/317)).
+* Stabilize [`SpinWaveTheoryKPM`](@ref). It now automatically selects the
+  polynomial order according to an error tolerance.
 
 ## v0.7.2
 (Sep 11, 2024)

examples/01_LSWT_CoRh2O4.jl

@@ -23,11 +23,11 @@ using Sunny, GLMakie
 # GLMakie from the [built-in package
 # manager](https://github.com/SunnySuite/Sunny.jl/wiki/Getting-started-with-Julia#the-built-in-julia-package-manager).
 
-# ### Units
+# ### Units system
 
-# The [`Units`](@ref Sunny.Units) object selects reference energy and length
-# scales, and uses these to provide physical constants. For example, `units.K`
-# returns one kelvin as 0.086 meV, where the Boltzmann constant is implicit.
+# The [`Units`](@ref Units) object selects reference energy and length scales,
+# and uses these to provide physical constants. For example, `units.K` returns
+# one kelvin as 0.086 meV, where the Boltzmann constant is implicit.
 
 units = Units(:meV, :angstrom);
 

examples/09_Disorder_KPM.jl

@@ -0,0 +1,121 @@
+# # 9. Disordered system with KPM
+#
+# This example uses the kernel polynomial method (KPM) to efficiently calculate
+# the neutron scattering spectrum of a disordered triangular antiferromagnet.
+# The model is inspired by YbMgGaO4, as studied in [Paddison et al, Nature
+# Phys., **13**, 117–122 (2017)](https://doi.org/10.1038/nphys3971) and [Zhu et
+# al, Phys. Rev. Lett. **119**, 157201
+# (2017)](https://doi.org/10.1103/PhysRevLett.119.157201). Disordered occupancy
+# of non-magnetic Mg/Ga sites can be modeled as a stochastic distribution of
+# exchange constants and ``g``-factors. Including this disorder introduces
+# broadening of the spin wave spectrum.
+
+using Sunny, GLMakie
+
+# Set up minimal triangular lattice system. Include antiferromagnetic exchange
+# interactions between nearest neighbor bonds. Energy minimization yields the
+# magnetic ground state with 120° angles between spins in triangular plaquettes.
+
+latvecs = lattice_vectors(1, 1, 10, 90, 90, 120)
+cryst = Crystal(latvecs, [[0, 0, 0]])
+sys = System(cryst, [1 => Moment(s=1/2, g=1)], :dipole; dims=(3, 3, 1))
+set_exchange!(sys, +1.0, Bond(1, 1, [1,0,0]))
+
+randomize_spins!(sys)
+minimize_energy!(sys)
+plot_spins(sys; color=[S[3] for S in sys.dipoles], ndims=2)
+
+# Select a ``𝐪``-space path for the spin wave calculations.
+
+qs = [[0, 0, 0], [1/3, 1/3, 0], [1/2, 0, 0], [0, 0, 0]]
+labels = ["Γ", "K", "M", "Γ"]
+path = q_space_path(cryst, qs, 150; labels)
+kernel = lorentzian(fwhm=0.4);
+
+# Perform a traditional spin wave calculation. The spectrum shows sharp modes
+# associated with coherent excitations about the K-point ordering wavevector,
+# ``𝐪 = [1/3, 1/3, 0]``.
+
+energies = range(0.0, 3.0, 150)
+swt = SpinWaveTheory(sys; measure=ssf_perp(sys))
+res = intensities(swt, path; energies, kernel)
+plot_intensities(res)
+
+# Use [`repeat_periodically`](@ref) to enlarge the system by a factor of 10 in
+# each dimension. Use [`to_inhomogeneous`](@ref) to disable symmetry
+# constraints, and allow for the addition of disordered interactions.
+
+sys_inhom = to_inhomogeneous(repeat_periodically(sys, (10, 10, 1)))
+
+# Use [`symmetry_equivalent_bonds`](@ref) to iterate over all nearest neighbor
+# bonds of the inhomogeneous system. Modify each AFM exchange with a noise term
+# that has variance of 1/3. The newly minimized energy configuration allows for
+# long wavelength modulations on top of the original 120° order.
+
+for (site1, site2, offset) in symmetry_equivalent_bonds(sys_inhom, Bond(1,1,[1,0,0]))
+    noise = randn()/3
+    set_exchange_at!(sys_inhom, 1.0 + noise, site1, site2; offset)
+end
+
+minimize_energy!(sys_inhom, maxiters=5_000)
+plot_spins(sys_inhom; color=[S[3] for S in sys_inhom.dipoles], ndims=2)
+
+# Traditional spin wave theory calculations become impractical for large system
+# sizes. Significant acceleration is possible with the [kernel polynomial
+# method](https://arxiv.org/abs/2312.08349). Enable it by selecting
+# [`SpinWaveTheoryKPM`](@ref) in place of the traditional
+# [`SpinWaveTheory`](@ref). Using KPM, the cost of an [`intensities`](@ref)
+# calculation becomes linear in system size and scales inversely with the width
+# of the line broadening `kernel`. Error tolerance is controlled through the
+# dimensionless `tol` parameter. A relatively small value, `tol = 0.01`, helps
+# to resolve the large intensities near the ordering wavevector. The alternative
+# choice `tol = 0.1` would be twice faster, but would introduce significant
+# numerical artifacts.
+#
+# Observe from the KPM calculation that disorder in the nearest-neighbor
+# exchange serves to broaden the discrete excitation bands into a continuum.
+
+swt = SpinWaveTheoryKPM(sys_inhom; measure=ssf_perp(sys_inhom), tol=0.01)
+res = intensities(swt, path; energies, kernel)
+plot_intensities(res)
+
+# Now apply a magnetic field of magnitude 7.5 (energy units) along the global
+# ``ẑ`` axis. This field fully polarizes the spins. Because gap opens, a larger
+# tolerance of `tol = 0.1` can be used to accelerate the KPM calculation without
+# sacrificing much accuracy. The resulting spin wave spectrum shows a sharp mode
+# at the Γ-point (zone center) that broadens into a continuum along the K and M
+# points (zone boundary).
+
+set_field!(sys_inhom, [0, 0, 7.5])
+randomize_spins!(sys_inhom)
+minimize_energy!(sys_inhom)
+
+energies = range(0.0, 9.0, 150)
+swt = SpinWaveTheoryKPM(sys_inhom; measure=ssf_perp(sys_inhom), tol=0.1)
+res = intensities(swt, path; energies, kernel)
+plot_intensities(res)
+
+# Add disorder to the ``z``-component of each magnetic moment ``g``-tensor. This
+# further broadens intensities, now across the entire path. Some intensity
+# modulation within the continuum is also apparent. This modulation is a
+# finite-size effect, and would be mitigated by enlarging the system beyond
+# 30×30 chemical cells.
+
+for site in eachsite(sys_inhom)
+    noise = randn()/6
+    sys_inhom.gs[site] = [1 0 0; 0 1 0; 0 0 1+noise]
+end
+
+swt = SpinWaveTheoryKPM(sys_inhom; measure=ssf_perp(sys_inhom), tol=0.1)
+res = intensities(swt, path; energies, kernel)
+plot_intensities(res)
+
+# For reference, the equivalent non-disordered system shows a single coherent
+# mode.
+
+set_field!(sys, [0, 0, 7.5])
+randomize_spins!(sys)
+minimize_energy!(sys)
+swt = SpinWaveTheory(sys; measure=ssf_perp(sys))
+res = intensities(swt, path; energies, kernel)
+plot_intensities(res)

examples/spinw_tutorials/SW02_AFM_Heisenberg_chain.jl

@@ -37,7 +37,7 @@ plot_spins(sys; ndims=2, ghost_radius=8)
 
 swt = SpinWaveTheory(sys; measure=ssf_perp(sys))
 qs = [[0,0,0], [1,0,0]]
-path = q_space_path(cryst, qs, 400)
+path = q_space_path(cryst, qs, 401)
 res = intensities_bands(swt, path)
 plot_intensities(res; units)
 

examples/spinw_tutorials/SW03_Frustrated_chain.jl

@@ -53,6 +53,6 @@ plot_spins(sys_enlarged; ndims=2)
 
 swt = SpinWaveTheorySpiral(sys; measure=ssf_perp(sys), k, axis)
 qs = [[0,0,0], [1,0,0]]
-path = q_space_path(cryst, qs, 400)
+path = q_space_path(cryst, qs, 401)
 res = intensities_bands(swt, path)
 plot_intensities(res; units)

examples/spinw_tutorials/SW04_Frustrated_square.jl

@@ -34,6 +34,6 @@ plot_spins(sys)
 
 swt = SpinWaveTheory(sys; measure=ssf_perp(sys))
 qs = [[0, 0, 0], [1/2, 0, 0], [1/2, 1/2, 0], [0, 0, 0]]
-path = q_space_path(cryst, qs, 400);
+path = q_space_path(cryst, qs, 400)
 res = intensities_bands(swt, path)
 plot_intensities(res; units)

src/KPM/Chebyshev.jl

@@ -56,3 +56,33 @@ function cheb_eval(x, bounds, coefs)
     end
     return ret
 end
+
+"""
+    cheb_moments_to_density(μs, N; γ=1)
+
+Transform Chebyshev expansion moments μ_m to densities ρ_n for discrete points
+x_n = cos[(π/N)(n+1/2)], where n = 0 … N-1. ⟦UNTESTED⟧.
+"""
+function cheb_moments_to_density(μs, N; γ=1)
+    M = length(μs)
+    @assert N >= M
+    xs = zeros(N)
+    ρs = zeros(N)
+    copy!(ρs, μs)
+    apply_jackson_kernel!(ρs)
+    plan = FFTW.plan_r2r!(ρs, FFTW.REDFT01)
+    mul!(ρs, plan, ρs)
+
+    for i in 1:N
+        n = i-1
+        x = cos((π/N) * (n+1/2))
+        push!(xs, x)
+        w = 1 / sqrt(1 - x^2)
+        ρs[i] *= w / π
+    end
+
+    xs .*= γ
+    ρs ./= γ
+
+    return (xs, ρs)
+end

src/KPM/Lanczos.jl

@@ -1,79 +1,112 @@
+# Reference implementation, without consideration of speed.
+function lanczos_ref(A, S, v; niters)
+    αs = Float64[]
+    βs = Float64[]
 
-"""
-    parallel_lanczos!(buf1, buf2, buf3; niters, mat_vec_mul!, inner_product!)
+    v /= sqrt(real(dot(v, S, v)))
+    w = A * S * v
+    α = real(dot(w, S, v))
+    w = w - α * v
+    push!(αs, α)
 
-Parallelized Lanczos. Takes functions for matrix-vector multiplication and inner
-products. Returns a list of tridiagonal matrices. `niters` sets the number of
-iterations used by the Lanczos algorithm.
+    for _ in 2:niters
+        β = sqrt(real(dot(w, S, w)))
+        iszero(β) && break
+        vp = w / β
+        w = A * S * vp
+        α = real(dot(w, S, vp))
+        w = w - α * vp - β * v
+        v = vp
+        push!(αs, α)
+        push!(βs, β)
+    end
 
-References:
-- [H. A. van der Vorst, Math. Comp. 39, 559-561
-  (1982)](https://doi.org/10.1090/s0025-5718-1982-0669648-0).
-- [M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156
-  (2011)](https://doi.org/10.1016/j.commatsci.2011.02.021).
-"""
-function parallel_lanczos!(vprev, v, w;  niters, mat_vec_mul!, inner_product)
-    # ...
+    return SymTridiagonal(αs, βs)
 end
 
 
-function lanczos(swt, q_reshaped, niters)
-    L = nbands(swt)
 
-    function inner_product(w, v)
-        return @views real(dot(w[1:L], v[1:L]) - dot(w[L+1:2L], v[L+1:2L]))
-    end
+# Use Lanczos iterations to find a truncated tridiagonal representation of A S,
+# up to similarity transformation. Here, A is any Hermitian matrix, while S is
+# both Hermitian and positive definite. Traditional Lanczos uses the identity
+# matrix, S = I. The extension to non-identity matrices S is as follows: Each
+# matrix-vector product A v becomes A S v, and each vector inner product w† v
+# becomes w† S v. The implementation below follows the notation of Wikipedia,
+# and is the most stable of the four variants considered by Paige [1]. Each
+# Lanczos iteration requires only one matrix-vector multiplication for A and S,
+# respectively.
 
-    function mat_vec_mul!(w, v, q_reshaped)
-        mul_dynamical_matrix!(swt, reshape(w, 1, :), reshape(v, 1, :), [q_reshaped])
-        view(w, L+1:2L) .*= -1
-    end
+# Similar generalizations of Lanczos have been considered in [2] and [3].
+#
+# 1. C. C. Paige, IMA J. Appl. Math., 373-381 (1972),
+# https://doi.org/10.1093%2Fimamat%2F10.3.373.
+# 2. H. A. van der Vorst, Math. Comp. 39, 559-561 (1982),
+# https://doi.org/10.1090/s0025-5718-1982-0669648-0
+# 3. M. Grüning, A. Marini, X. Gonze, Comput. Mater. Sci. 50, 2148-2156 (2011),
+# https://doi.org/10.1016/j.commatsci.2011.02.021.
+function lanczos(mulA!, mulS!, v; niters)
+    αs = Float64[]
+    βs = Float64[]
+
+    vp  = zero(v)
+    Sv  = zero(v)
+    Svp = zero(v)
+    w   = zero(v)
+    Sw  = zero(v)
 
-    αs = zeros(Float64, niters)   # Diagonal part
-    βs = zeros(Float64, niters-1) # Off-diagonal part
-
-    vprev = randn(ComplexF64, 2L)
-    v = zeros(ComplexF64, 2L)
-    w = zeros(ComplexF64, 2L)
-
-    mat_vec_mul!(w, vprev, q_reshaped)
-    b0 = sqrt(inner_product(vprev, w))
-    vprev = vprev / b0
-    w = w / b0
-    αs[1] = inner_product(w, w)
-    @. w = w - αs[1]*vprev
-    for j in 2:niters
-        @. v = w
-        mat_vec_mul!(w, v, q_reshaped)
-        βs[j-1] = sqrt(inner_product(v, w)) # maybe v?
-        if abs(βs[j-1]) < 1e-12
-            # TODO: Terminate early if all β are small
-            βs[j-1] = 0
-            @. v = w = 0
-        else
-            @. v = v / βs[j-1]
-            @. w = w / βs[j-1]
-        end
-        αs[j] = inner_product(w, w)
-        @. w = w - αs[j]*v - βs[j-1]*vprev
-        v, vprev = vprev, v
+    mulS!(Sv, v)
+    c = sqrt(real(dot(v, Sv)))
+    @. v /= c
+    @. Sv /= c
+
+    mulA!(w, Sv)
+    α = real(dot(w, Sv))
+    @. w = w - α * v
+    mulS!(Sw, w)
+    push!(αs, α)
+
+    for _ in 2:niters
+        β = sqrt(real(dot(w, Sw)))
+        iszero(β) && break
+        @. vp = w / β
+        @. Svp = Sw / β
+        mulA!(w, Svp)
+        α = real(dot(w, Svp))
+        @. w = w - α * vp - β * v
+        mulS!(Sw, w)
+        @. v = vp
+        push!(αs, α)
+        push!(βs, β)
     end
 
     return SymTridiagonal(αs, βs)
 end
 
 
 """
-    eigbounds(A, niters; extend=0.0)
+    eigbounds(swt, niters)
 
-Returns estimates of the extremal eigenvalues of Hermitian matrix `A` using the
-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.
+Returns estimates of the extremal eigenvalues of the matrix (Ĩ D)^2 using the
+Lanczos algorithm. Here Ĩ is the Bogoliubov identity and D is the dynamical
+matrix from LSWT for the wavevector q_reshaped. `niters` should be given a value
+smaller than the dimension of `A`.
 """
-function eigbounds(swt, q_reshaped, niters; extend)
-    tridiag = lanczos(swt, Vec3(q_reshaped), niters)
-    lo, hi = extrema(eigvals!(tridiag)) # TODO: pass irange=1 and irange=2L
-    slack = extend*(hi-lo)
-    return lo-slack, hi+slack
+function eigbounds(swt, q_reshaped, niters)
+    q_reshaped = Vec3(q_reshaped)
+    L = nbands(swt)
+
+    # w = Ĩ v
+    function mulA!(w, v)
+        @views w[1:L]    = +v[1:L]
+        @views w[L+1:2L] = -v[L+1:2L]
+    end
+
+    # w = D v
+    function mulS!(w, v)
+        mul_dynamical_matrix!(swt, reshape(w, 1, :), reshape(v, 1, :), [q_reshaped])
+    end
+
+    v = randn(ComplexF64, 2L)
+    tridiag = lanczos(mulA!, mulS!, v; niters)
+    return eigmin(tridiag), eigmax(tridiag)
 end