From 5d451200c459dff2547e9a23bdb1fdfc9eb91e71 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 01/29] Progress --- src/System/Ewald.jl | 37 +++++++++++++++++++++++-------------- 1 file changed, 23 insertions(+), 14 deletions(-) diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl index f6fe2e0c2..352f3d8c6 100644 --- a/src/System/Ewald.jl +++ b/src/System/Ewald.jl @@ -36,11 +36,20 @@ end # Tensor product of 3-vectors (⊗)(a::Vec3,b::Vec3) = reshape(kron(a,b), 3, 3) -# Precompute the dipole-dipole interaction matrix A and its Fourier transform -# F[A] -function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int}) :: Array{Mat3, 5} + +function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int}) + precompute_dipole_ewald_aux(cryst, latsize, Vec3(0,0,0), cos, Val{Float64}()) +end + +# Precompute the dipole-dipole interaction matrix A. If q_reshaped is zero then +# all calculated values will be real. For example, `exp(i (q+k)⋅r) → cos(k⋅r)` +# because the imaginary part cancels in the symmetric sum over ±k. Knowing this, +# one can replace `cis(x) ≡ exp(i x) = cos(x) + i sin(x)` with just `cos(x)` for +# efficiency. The parameter `T ∈ {Float64, ComplexF64}` controls the return type +# in a type-stable way. +function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, q_reshaped, cis, ::Val{T}) where T na = natoms(cryst) - A = zeros(Mat3, latsize..., na, na) + A = zeros(SMatrix{3, 3, T, 9}, latsize..., na, na) # Superlattice vectors and reciprocals for the full system volume sys_size = diagm(Vec3(latsize)) @@ -51,7 +60,7 @@ function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int}) :: Arra I₃ = Mat3(I) V = det(latvecs) L = cbrt(V) - # Roughly balances the real and Fourier space costs + # Roughly balances the real and Fourier space costs. Note that σ = 1/(√2 λ) σ = L/3 σ² = σ*σ σ³ = σ^3 @@ -70,14 +79,16 @@ function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int}) :: Arra # println("nmax $nmax mmax $mmax") for cell in CartesianIndices(latsize), j in 1:na, i in 1:na - acc = zero(Mat3) + acc = zero(eltype(A)) cell_offset = Vec3(cell[1]-1, cell[2]-1, cell[3]-1) Δr = cryst.latvecs * (cell_offset + cryst.positions[j] - cryst.positions[i]) - + ##################################################### ## Real space part for n1 = -nmax[1]:nmax[1], n2 = -nmax[2]:nmax[2], n3 = -nmax[3]:nmax[3] - rvec = Δr + latvecs * Vec3(n1, n2, n3) + n = Vec3(n1, n2, n3) + phase = cis(2π * dot(q_reshaped, n)) + rvec = Δr + latvecs * n r² = rvec⋅rvec if 0 < r² <= rmax*rmax r = √r² @@ -85,20 +96,18 @@ function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int}) :: Arra rhat = rvec/r erfc0 = erfc(r/(√2*σ)) gauss0 = √(2/π) * (r/σ) * exp(-r²/2σ²) - acc += (1/2) * ((I₃/r³) * (erfc0 + gauss0) - (3(rhat⊗rhat)/r³) * (erfc0 + (1+r²/3σ²) * gauss0)) + acc += (phase/2) * ((I₃/r³) * (erfc0 + gauss0) - (3(rhat⊗rhat)/r³) * (erfc0 + (1+r²/3σ²) * gauss0)) end end ##################################################### ## Fourier space part for m1 = -mmax[1]:mmax[1], m2 = -mmax[2]:mmax[2], m3 = -mmax[3]:mmax[3] - k = recipvecs * Vec3(m1, m2, m3) + m = Vec3(m1, m2, m3) + k = recipvecs * (m + q_reshaped) k² = k⋅k if 0 < k² <= kmax*kmax - # Replace exp(-ikr) -> cos(kr). It's valid to drop the imaginary - # component because it will cancel in the interaction that exchanges - # i ↔ j. - acc += (4π/2V) * (exp(-σ²*k²/2) / k²) * (k⊗k) * cos(k⋅Δr) + acc += (4π/2V) * (exp(-σ²*k²/2) / k²) * (k⊗k) * cis(k⋅Δr) end end From cd6239aaa57985d232a7e94abf47aee1ca25b26d Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 02/29] Interaction conventions --- src/System/Ewald.jl | 55 ++++++++++++++++++++++++--------------------- 1 file changed, 29 insertions(+), 26 deletions(-) diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl index 352f3d8c6..cf179caed 100644 --- a/src/System/Ewald.jl +++ b/src/System/Ewald.jl @@ -2,7 +2,7 @@ function Ewald(sys::System{N}) where N (; crystal, latsize, units) = sys - A = (units.μ0/4π) * precompute_dipole_ewald(crystal, latsize) + A = precompute_dipole_ewald(crystal, latsize, units.μ0) na = natoms(crystal) μ = zeros(Vec3, latsize..., na) @@ -37,17 +37,24 @@ end (⊗)(a::Vec3,b::Vec3) = reshape(kron(a,b), 3, 3) -function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int}) - precompute_dipole_ewald_aux(cryst, latsize, Vec3(0,0,0), cos, Val{Float64}()) +function precompute_dipole_ewald(cryst::Crystal, latsize::NTuple{3,Int}, μ0) + precompute_dipole_ewald_aux(cryst, latsize, μ0, Vec3(0,0,0), cos, Val{Float64}()) end -# Precompute the dipole-dipole interaction matrix A. If q_reshaped is zero then -# all calculated values will be real. For example, `exp(i (q+k)⋅r) → cos(k⋅r)` -# because the imaginary part cancels in the symmetric sum over ±k. Knowing this, -# one can replace `cis(x) ≡ exp(i x) = cos(x) + i sin(x)` with just `cos(x)` for -# efficiency. The parameter `T ∈ {Float64, ComplexF64}` controls the return type -# in a type-stable way. -function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, q_reshaped, cis, ::Val{T}) where T +# Precompute the pairwise interaction matrix A beween magnetic moments μ. For +# q_reshaped = 0, this yields the usual Ewald energy, E = μᵢ Aᵢⱼ μⱼ / 2. Nonzero +# q_reshaped is useful in spin wave theory. Physically, this amounts to a +# modification of the periodic boundary conditions, such that μ(q) can be +# incommensurate with the magnetic cell. In all cases, the energy is E = μᵢ(-q) +# Aᵢⱼ(-q) μⱼ(q) / 2 in Fourier space, where q should be interpreted as a Fourier +# transform of the cell offset. +# +# As an optimization, this function returns real values when q_reshaped is zero. +# Effectively, one can replace `exp(i (q+k)⋅r) → cos(k⋅r)` because the imaginary +# part cancels in the symmetric sum over ±k. Specifically, replace `cis(x) ≡ +# exp(i x) = cos(x) + i sin(x)` with just `cos(x)` for efficiency. The parameter +# `T ∈ {Float64, ComplexF64}` controls the return type in a type-stable way. +function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, μ0, q_reshaped, cis, ::Val{T}) where T na = natoms(cryst) A = zeros(SMatrix{3, 3, T, 9}, latsize..., na, na) @@ -96,7 +103,7 @@ function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, q_r rhat = rvec/r erfc0 = erfc(r/(√2*σ)) gauss0 = √(2/π) * (r/σ) * exp(-r²/2σ²) - acc += (phase/2) * ((I₃/r³) * (erfc0 + gauss0) - (3(rhat⊗rhat)/r³) * (erfc0 + (1+r²/3σ²) * gauss0)) + acc += phase * (μ0/4π) * ((I₃/r³) * (erfc0 + gauss0) - (3(rhat⊗rhat)/r³) * (erfc0 + (1+r²/3σ²) * gauss0)) end end @@ -107,14 +114,14 @@ function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, q_r k = recipvecs * (m + q_reshaped) k² = k⋅k if 0 < k² <= kmax*kmax - acc += (4π/2V) * (exp(-σ²*k²/2) / k²) * (k⊗k) * cis(k⋅Δr) + acc += (μ0/V) * (exp(-σ²*k²/2) / k²) * (k⊗k) * cis(k⋅Δr) end end ##################################################### ## Remove self energies if iszero(Δr) - acc += - I₃/(3√(2π)*σ³) + acc += - μ0*I₃/(3(2π)^(3/2)*σ³) end # For sites site1=(cell1, i) and site2=(cell2, j) offset by an amount @@ -147,15 +154,13 @@ function ewald_energy(sys::System{N}) where N else @views Fμ[:, 2:end, :, :, :] .*= √2 end - # * The field is a cross correlation, h(r) = - 2 A(Δr) s(r+Δr) = - 2A⋆s, or - # in Fourier space, F[h] = - 2 conj(F[A]) F[s] - # * The energy is an inner product, E = - (1/2)s⋅h, or using Parseval's - # theorem, E = - (1/2) conj(F[s]) F[h] / N - # * Combined, the result is: E = conj(F[s]) conj(F[A]) F[s] / N + + # In real space, E = μ (A ⋆ μ) / 2. In Fourier space, the convolution + # becomes an ordinary product using Parseval's theorem. (_, m1, m2, m3, na) = size(Fμ) ms = CartesianIndices((m1, m2, m3)) @inbounds for j in 1:na, i in 1:na, m in ms, α in 1:3, β in 1:3 - E += real(conj(Fμ[α, m, i]) * conj(FA[α, β, m, i, j]) * Fμ[β, m, j]) + E += (1/2) * real(conj(Fμ[α, m, i]) * conj(FA[α, β, m, i, j]) * Fμ[β, m, j]) end return E / prod(latsize) end @@ -186,7 +191,7 @@ function accum_ewald_grad!(∇E, dipoles, sys::System{N}) where N mul!(ϕr, ift_plan, Fϕ) for site in eachsite(sys) - ∇E[site] += 2 * units.μB * (gs[site]' * ϕ[site]) + ∇E[site] += units.μB * (gs[site]' * ϕ[site]) end end @@ -197,11 +202,9 @@ function ewald_pairwise_grad_at(sys::System{N}, site1, site2) where N cell_offset = mod.(Tuple(to_cell(site2)-to_cell(site1)), latsize) cell = CartesianIndex(cell_offset .+ (1,1,1)) - # A prefactor of 2 is always appropriate here. If site1 == site2, it - # accounts for the quadratic dependence on the dipole. If site1 != site2, it - # accounts for energy contributions from both ordered pairs (site1, site2) - # and (site2, site1). - return 2 * units.μB^2 * gs[site1]' * ewald.A[cell, to_atom(site1), to_atom(site2)] * gs[site2] * sys.dipoles[site2] + # The factor of 1/2 in the energy formula `E = μ (A ⋆ μ) / 2` disappears due + # to quadratic appearance of μ = - μB g S. + return units.μB^2 * gs[site1]' * ewald.A[cell, to_atom(site1), to_atom(site2)] * gs[site2] * sys.dipoles[site2] end # Calculate the field dE/ds at `site` generated by all `dipoles`. @@ -221,5 +224,5 @@ function ewald_energy_delta(sys::System{N}, site, s::Vec3) where N Δμ = units.μB * (sys.gs[site] * Δs) i = to_atom(site) ∇E = ewald_grad_at(sys, site) - return Δs⋅∇E + dot(Δμ, ewald.A[1, 1, 1, i, i], Δμ) + return Δs⋅∇E + dot(Δμ, ewald.A[1, 1, 1, i, i], Δμ) / 2 end From 57fd03838da0ecdc96f97ee4f3a0da50b8df61f2 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 03/29] Fix sign and rounding --- src/System/Ewald.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl index cf179caed..d37246778 100644 --- a/src/System/Ewald.jl +++ b/src/System/Ewald.jl @@ -94,7 +94,6 @@ function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, μ0 ## Real space part for n1 = -nmax[1]:nmax[1], n2 = -nmax[2]:nmax[2], n3 = -nmax[3]:nmax[3] n = Vec3(n1, n2, n3) - phase = cis(2π * dot(q_reshaped, n)) rvec = Δr + latvecs * n r² = rvec⋅rvec if 0 < r² <= rmax*rmax @@ -103,6 +102,7 @@ function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, μ0 rhat = rvec/r erfc0 = erfc(r/(√2*σ)) gauss0 = √(2/π) * (r/σ) * exp(-r²/2σ²) + phase = cis(2π * dot(q_reshaped, n)) acc += phase * (μ0/4π) * ((I₃/r³) * (erfc0 + gauss0) - (3(rhat⊗rhat)/r³) * (erfc0 + (1+r²/3σ²) * gauss0)) end end @@ -111,10 +111,11 @@ function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, μ0 ## Fourier space part for m1 = -mmax[1]:mmax[1], m2 = -mmax[2]:mmax[2], m3 = -mmax[3]:mmax[3] m = Vec3(m1, m2, m3) - k = recipvecs * (m + q_reshaped) + k = recipvecs * (m + q_reshaped - round.(q_reshaped)) k² = k⋅k if 0 < k² <= kmax*kmax - acc += (μ0/V) * (exp(-σ²*k²/2) / k²) * (k⊗k) * cis(k⋅Δr) + phase = cis(-k⋅Δr) + acc += phase * (μ0/V) * (exp(-σ²*k²/2) / k²) * (k⊗k) end end From 550e676bc2a2f046a97129a4dd132de02b244f84 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 04/29] Implementation attempt 1 --- src/SpinWaveTheory/HamiltonianDipole.jl | 61 +++++++++++++++++++++++++ src/SpinWaveTheory/SpinWaveTheory.jl | 4 -- src/System/Ewald.jl | 8 ++++ 3 files changed, 69 insertions(+), 4 deletions(-) diff --git a/src/SpinWaveTheory/HamiltonianDipole.jl b/src/SpinWaveTheory/HamiltonianDipole.jl index 86c72e463..2c05788b6 100644 --- a/src/SpinWaveTheory/HamiltonianDipole.jl +++ b/src/SpinWaveTheory/HamiltonianDipole.jl @@ -97,6 +97,67 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r H21[i, i] += conj(A2) end + # Add long-range dipole-dipole + if !isnothing(sys.ewald) + N = sys.Ns[1] + S = (N-1)/2 + μB² = units.μB^2 + Rs = local_rotations + A = precompute_dipole_ewald_aux(sys.crystal, (1,1,1), units.μ0, q_reshaped, cis, Val{ComplexF64}()) + A = reshape(A, L, L) + + # Loop over all unique sublattice pairs (i < j) + for i in 1:L, j in 1:L # i:L + + # DEBUG + i == j && continue + + # An ordered pair of magnetic moments contribute (μᵢ A μⱼ)/2 to the + # energy. Note that μ = -μB g S + J = μB² * gs[i]' * A[i, j] * gs[j] + + # Account for contributions from both (i, j) and (j, i) + #= + if i != j + J *= 2 + end + =# + + # Perform same transformation as appears in usual bilinear exchange. + # R denotes rotation from lab frame into global frame. + J = S * Rs[i]' * J * Rs[j] + + #= + 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]) + + H11[i, j] += Q * phase + H11[j, i] += conj(Q) * conj(phase) + H22[i, j] += conj(Q) * phase + H22[j, i] += Q * conj(phase) + + H21[i, j] += P * phase + H21[j, i] += P * conj(phase) + H12[i, j] += conj(P) * phase + H12[j, i] += conj(P) * conj(phase) + =# + + Q⁻ = 0.25 * (J[1, 1] + J[2, 2] - im*(J[1, 2] - J[2, 1])) + Q⁺ = 0.25 * (J[1, 1] + J[2, 2] + im*(J[1, 2] - J[2, 1])) + P⁻ = 0.25 * (J[1, 1] - J[2, 2] - im*(J[1, 2] + J[2, 1])) + P⁺ = 0.25 * (J[1, 1] - J[2, 2] + im*(J[1, 2] + J[2, 1])) + + H11[i, j] += Q⁻ + H11[j, i] += conj(Q⁻) + H22[i, j] += Q⁺ + H22[j, i] += conj(Q⁺) + H21[i, j] += P⁻ + H12[j, i] += conj(P⁻) + H21[j, i] += P⁺ + H12[i, j] += conj(P⁺) + end + end + # H must be hermitian up to round-off errors @assert diffnorm2(H, H') < 1e-12 diff --git a/src/SpinWaveTheory/SpinWaveTheory.jl b/src/SpinWaveTheory/SpinWaveTheory.jl index 0d905abce..cce60fbee 100644 --- a/src/SpinWaveTheory/SpinWaveTheory.jl +++ b/src/SpinWaveTheory/SpinWaveTheory.jl @@ -27,10 +27,6 @@ struct SpinWaveTheory end function SpinWaveTheory(sys::System{N}; energy_ϵ::Float64=1e-8, observables=nothing, correlations=nothing, apply_g=true) where N - if !isnothing(sys.ewald) - error("SpinWaveTheory does not yet support long-range dipole-dipole interactions.") - end - # Reshape into single unit cell. A clone will always be performed, even if # no reshaping happens. cellsize_mag = cell_shape(sys) * diagm(Vec3(sys.latsize)) diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl index d37246778..e9cb94fd8 100644 --- a/src/System/Ewald.jl +++ b/src/System/Ewald.jl @@ -131,6 +131,14 @@ function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, μ0 A[cell, i, j] = acc end + #= + for j in 1:na, i in 1:na + X = (A[1,1,1, i,j] + A[1,1,1, j,i]')/2 + A[1,1,1, i,j] = X + A[1,1,1, j,i] = X' + end + =# + return A end From 2f91631d5182228906735fbd9e630c41898fc1e8 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 05/29] Debugging --- src/SpinWaveTheory/HamiltonianDipole.jl | 29 ++++++++++++++++--------- 1 file changed, 19 insertions(+), 10 deletions(-) diff --git a/src/SpinWaveTheory/HamiltonianDipole.jl b/src/SpinWaveTheory/HamiltonianDipole.jl index 2c05788b6..971edd868 100644 --- a/src/SpinWaveTheory/HamiltonianDipole.jl +++ b/src/SpinWaveTheory/HamiltonianDipole.jl @@ -106,20 +106,16 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r A = precompute_dipole_ewald_aux(sys.crystal, (1,1,1), units.μ0, q_reshaped, cis, Val{ComplexF64}()) A = reshape(A, L, L) - # Loop over all unique sublattice pairs (i < j) - for i in 1:L, j in 1:L # i:L - - # DEBUG - i == j && continue - + # Loop over sublattice pairs + for i in 1:L, j in 1:L # An ordered pair of magnetic moments contribute (μᵢ A μⱼ)/2 to the # energy. Note that μ = -μB g S J = μB² * gs[i]' * A[i, j] * gs[j] - - # Account for contributions from both (i, j) and (j, i) + J *= 1/2 # TODO: save factor of two by taking (j in i:L) and using: #= - if i != j - J *= 2 + # Undo double counting in case where (i == j) + if i == j + J *= 1/2 end =# @@ -140,6 +136,11 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r H21[j, i] += P * conj(phase) H12[i, j] += conj(P) * phase H12[j, i] += conj(P) * conj(phase) + + H11[i, i] -= 0.5 * J[3, 3] + H11[j, j] -= 0.5 * J[3, 3] + H22[i, i] -= 0.5 * J[3, 3] + H22[j, j] -= 0.5 * J[3, 3] =# Q⁻ = 0.25 * (J[1, 1] + J[2, 2] - im*(J[1, 2] - J[2, 1])) @@ -155,6 +156,11 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r H12[j, i] += conj(P⁻) H21[j, i] += P⁺ H12[i, j] += conj(P⁺) + + H11[i, i] -= 0.5 * J[3, 3] + H11[j, j] -= 0.5 * J[3, 3] + H22[i, i] -= 0.5 * J[3, 3] + H22[j, j] -= 0.5 * J[3, 3] end end @@ -164,6 +170,9 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r # Make H exactly hermitian hermitianpart!(H) + display(H) + println(eigvals(H)) + # Add small constant shift for positive-definiteness for i in 1:2L H[i, i] += swt.energy_ϵ From c22fbf4344838e593feb6dadf4e75f810bdbfa7a Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 06/29] Progress --- src/SpinWaveTheory/HamiltonianDipole.jl | 64 ++++++++++--------------- src/System/Ewald.jl | 10 +--- 2 files changed, 26 insertions(+), 48 deletions(-) diff --git a/src/SpinWaveTheory/HamiltonianDipole.jl b/src/SpinWaveTheory/HamiltonianDipole.jl index 971edd868..044ba53c6 100644 --- a/src/SpinWaveTheory/HamiltonianDipole.jl +++ b/src/SpinWaveTheory/HamiltonianDipole.jl @@ -103,51 +103,33 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r S = (N-1)/2 μB² = units.μB^2 Rs = local_rotations + + # Interaction matrix for wavevector q A = precompute_dipole_ewald_aux(sys.crystal, (1,1,1), units.μ0, q_reshaped, cis, Val{ComplexF64}()) A = reshape(A, L, L) + # Interaction matrix for wavevector (0,0,0) + A0 = precompute_dipole_ewald(sys.crystal, (1,1,1), units.μ0) + A0 = reshape(A0, L, L) + # Loop over sublattice pairs for i in 1:L, j in 1:L # An ordered pair of magnetic moments contribute (μᵢ A μⱼ)/2 to the - # energy. Note that μ = -μB g S - J = μB² * gs[i]' * A[i, j] * gs[j] - J *= 1/2 # TODO: save factor of two by taking (j in i:L) and using: - #= - # Undo double counting in case where (i == j) - if i == j - J *= 1/2 - end - =# + # energy. A symmetric contribution will appear for the bond reversal + # (i, j) → (j, i). Note that μ = -μB g S. + J = μB² * gs[i]' * A[i, j] * gs[j] / 2 + J0 = μB² * gs[i]' * A0[i, j] * gs[j] / 2 # Perform same transformation as appears in usual bilinear exchange. - # R denotes rotation from lab frame into global frame. + # Rⱼ denotes a rotation from ẑ to the ground state dipole Sⱼ. J = S * Rs[i]' * J * Rs[j] + J0 = S * Rs[i]' * J0 * Rs[j] - #= - 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]) - - H11[i, j] += Q * phase - H11[j, i] += conj(Q) * conj(phase) - H22[i, j] += conj(Q) * phase - H22[j, i] += Q * conj(phase) - - H21[i, j] += P * phase - H21[j, i] += P * conj(phase) - H12[i, j] += conj(P) * phase - H12[j, i] += conj(P) * conj(phase) - - H11[i, i] -= 0.5 * J[3, 3] - H11[j, j] -= 0.5 * J[3, 3] - H22[i, i] -= 0.5 * J[3, 3] - H22[j, j] -= 0.5 * J[3, 3] - =# - + # Interactions for Jˣˣ, Jʸʸ, Jˣʸ, and Jʸˣ at wavevector q. Q⁻ = 0.25 * (J[1, 1] + J[2, 2] - im*(J[1, 2] - J[2, 1])) Q⁺ = 0.25 * (J[1, 1] + J[2, 2] + im*(J[1, 2] - J[2, 1])) P⁻ = 0.25 * (J[1, 1] - J[2, 2] - im*(J[1, 2] + J[2, 1])) P⁺ = 0.25 * (J[1, 1] - J[2, 2] + im*(J[1, 2] + J[2, 1])) - H11[i, j] += Q⁻ H11[j, i] += conj(Q⁻) H22[i, j] += Q⁺ @@ -157,10 +139,11 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r H21[j, i] += P⁺ H12[i, j] += conj(P⁺) - H11[i, i] -= 0.5 * J[3, 3] - H11[j, j] -= 0.5 * J[3, 3] - H22[i, i] -= 0.5 * J[3, 3] - H22[j, j] -= 0.5 * J[3, 3] + # Interactions for Jᶻᶻ at wavevector (0,0,0). + H11[i, i] -= 0.5 * J0[3, 3] + H11[j, j] -= 0.5 * J0[3, 3] + H22[i, i] -= 0.5 * J0[3, 3] + H22[j, j] -= 0.5 * J0[3, 3] end end @@ -168,10 +151,7 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r @assert diffnorm2(H, H') < 1e-12 # Make H exactly hermitian - hermitianpart!(H) - - display(H) - println(eigvals(H)) + hermitianpart!(H) # Add small constant shift for positive-definiteness for i in 1:2L @@ -305,6 +285,10 @@ function multiply_by_hamiltonian_dipole_aux!(y, x, phasebuf, qphase, swt) y[q, j, 2] += P * conj(phasebuf[q]) * x[q, i, 1] end end + + if !isnothing(sys.ewald) + error("Ewald not supported") + end end end @@ -312,4 +296,4 @@ 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 diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl index e9cb94fd8..06ae79654 100644 --- a/src/System/Ewald.jl +++ b/src/System/Ewald.jl @@ -131,14 +131,8 @@ function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, μ0 A[cell, i, j] = acc end - #= - for j in 1:na, i in 1:na - X = (A[1,1,1, i,j] + A[1,1,1, j,i]')/2 - A[1,1,1, i,j] = X - A[1,1,1, j,i] = X' - end - =# - + # TODO: Verify that A[off, i, j] ≈ A[-off, j, i]' + return A end From 5220ac967da2f0277195199fa6d02a0a3689c527 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 07/29] Error for SU(N) mode --- src/SpinWaveTheory/HamiltonianSUN.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/SpinWaveTheory/HamiltonianSUN.jl b/src/SpinWaveTheory/HamiltonianSUN.jl index be514be70..4b20578b0 100644 --- a/src/SpinWaveTheory/HamiltonianSUN.jl +++ b/src/SpinWaveTheory/HamiltonianSUN.jl @@ -65,7 +65,7 @@ end # Set the dynamical quadratic Hamiltonian matrix in SU(N) mode. function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_reshaped::Vec3) - (; sys, data) = swt + (; sys) = swt L = nbands(swt) # Number of quasiparticle bands @assert size(H) == (2L, 2L) @@ -73,12 +73,11 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh # Clear the Hamiltonian H .= 0 - # Add pair interactions that use explicit bases for (atom, int) in enumerate(sys.interactions_union) - - # Set the onsite term + # Set the onsite term for atom swt_onsite_coupling!(H, int.onsite, swt, atom) + # Set pair interactions that originate for atom for coupling in int.pair # Extract information common to bond (; isculled, bond) = coupling @@ -91,6 +90,10 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh end end + if !isnothing(sys.ewald) + error("Ewald not yet supported in LSWT for :SUN mode") + end + # H must be hermitian up to round-off errors @assert diffnorm2(H, H') < 1e-12 From 2519d4e0e7e6c8d37bb59ddbc16d9ba0bbfd730c Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 08/29] Dipole-dipole not broken? --- test/test_lswt.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_lswt.jl b/test/test_lswt.jl index e52742b4c..3b0bb2277 100644 --- a/test/test_lswt.jl +++ b/test/test_lswt.jl @@ -332,12 +332,13 @@ end end +#= @testitem "Dipole-dipole unimplemented" begin sys = System(Sunny.diamond_crystal(),(1,1,1),[SpinInfo(1,S=1/2,g=2)],:SUN;seed = 0) enable_dipole_dipole!(sys) @test_throws "SpinWaveTheory does not yet support long-range dipole-dipole interactions." SpinWaveTheory(sys) end - +=# @testitem "SW15-Langasite" begin # Ba3NbFe3Si2O14 From c2376b3ea7177988becb463d033c26ea06f4e6f8 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 09/29] Avoid recomputing Ewald --- src/SpinWaveTheory/HamiltonianDipole.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/SpinWaveTheory/HamiltonianDipole.jl b/src/SpinWaveTheory/HamiltonianDipole.jl index 044ba53c6..ef9dea865 100644 --- a/src/SpinWaveTheory/HamiltonianDipole.jl +++ b/src/SpinWaveTheory/HamiltonianDipole.jl @@ -109,7 +109,8 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r A = reshape(A, L, L) # Interaction matrix for wavevector (0,0,0) - A0 = precompute_dipole_ewald(sys.crystal, (1,1,1), units.μ0) + A0 = sys.ewald.A + # @assert A0 ≈ precompute_dipole_ewald(sys.crystal, (1,1,1), units.μ0) A0 = reshape(A0, L, L) # Loop over sublattice pairs From 36cbc5eb9594e3e09d525b8fd4a69a32829bdfd3 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 10/29] Refactors --- src/SpinWaveTheory/HamiltonianDipole.jl | 172 +++++++------- src/SpinWaveTheory/HamiltonianSUN.jl | 290 ++++++++++-------------- test/test_lswt.jl | 25 +- 3 files changed, 217 insertions(+), 270 deletions(-) diff --git a/src/SpinWaveTheory/HamiltonianDipole.jl b/src/SpinWaveTheory/HamiltonianDipole.jl index ef9dea865..3c0930955 100644 --- a/src/SpinWaveTheory/HamiltonianDipole.jl +++ b/src/SpinWaveTheory/HamiltonianDipole.jl @@ -2,6 +2,7 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_reshaped::Vec3) (; sys, data) = swt (; local_rotations, stevens_coefs) = data + (; extfield, gs, units) = sys N = swt.sys.Ns[1] S = (N-1)/2 @@ -16,25 +17,30 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r H21 = view(H, L+1:2L, 1:L) H22 = view(H, L+1:2L, L+1:2L) - # Add Zeeman term - (; extfield, gs, units) = sys - for i in 1:L + for (i, int) in enumerate(sys.interactions_union) + # Zeeman term B = units.μB * (gs[1, 1, 1, i]' * extfield[1, 1, 1, i]) B′ = dot(B, local_rotations[i][:, 3]) / 2 H11[i, i] += B′ H22[i, i] += B′ - end - # Add pairwise terms - for ints in sys.interactions_union + # Single-ion anisotropy + (; c2, c4, c6) = stevens_coefs[i] + A1 = -3S*c2[3] - 40*S^3*c4[5] - 168*S^5*c6[7] + A2 = im*(S*c2[5] + 6S^3*c4[7] + 16S^5*c6[9]) + (S*c2[1] + 6S^3*c4[3] + 16S^5*c6[5]) + H11[i, i] += A1 + H22[i, i] += A1 + H12[i, i] += A2 + H21[i, i] += conj(A2) - # Bilinear exchange - for coupling in ints.pair + # Pair interactions + for coupling in int.pair (; isculled, bond) = coupling isculled && break - i, j = bond.i, bond.j + (; i, j) = bond phase = exp(2π*im * dot(q_reshaped, bond.n)) # Phase associated with periodic wrapping + # Bilinear exchange if !iszero(coupling.bilin) J = coupling.bilin # This is Rij in previous notation (transformed exchange matrix) @@ -86,17 +92,6 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r end end - # Add single-ion anisotropy - for i in 1:L - (; c2, c4, c6) = stevens_coefs[i] - A1 = -3S*c2[3] - 40*S^3*c4[5] - 168*S^5*c6[7] - A2 = im*(S*c2[5] + 6S^3*c4[7] + 16S^5*c6[9]) + (S*c2[1] + 6S^3*c4[3] + 16S^5*c6[5]) - H11[i, i] += A1 - H22[i, i] += A1 - H12[i, i] += A2 - H21[i, i] += conj(A2) - end - # Add long-range dipole-dipole if !isnothing(sys.ewald) N = sys.Ns[1] @@ -160,26 +155,15 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r end end -# Calculate H*x for many qs simultaneously. -function multiply_by_hamiltonian_dipole(x::Array{ComplexF64, 2}, swt::SpinWaveTheory, qs_reshaped::Array{Vec3}) - # Preallocate buffers. - L = nbands(swt) - @assert size(x, 2) == 2L - y = zeros(ComplexF64, (size(qs_reshaped)..., 2L)) - phasebuf = zeros(ComplexF64, length(qs_reshaped)) - - # Precompute e^{2πq_α} components. - qphase = map(qs_reshaped) do q - (exp(2π*im*q[1]), exp(2π*im*q[2]), exp(2π*im*q[3])) - end - - # Perform batched matrix-vector multiply. - multiply_by_hamiltonian_dipole_aux!(reshape(y, (length(qs_reshaped), 2L)), x, phasebuf, qphase, swt) - return y +function multiply_by_hamiltonian_dipole(x::Array{ComplexF64, 2}, swt::SpinWaveTheory, qs_reshaped::Array{Vec3}) + y = zero(x) + multiply_by_hamiltonian_dipole!(y, x, swt, qs_reshaped) + return y end -function multiply_by_hamiltonian_dipole_aux!(y, x, phasebuf, qphase, swt) +function multiply_by_hamiltonian_dipole!(y::Array{ComplexF64, 2}, x::Array{ComplexF64, 2}, swt::SpinWaveTheory, qs_reshaped::Array{Vec3}; + phases=zeros(ComplexF64, size(qs_reshaped))) (; sys, data) = swt (; stevens_coefs, local_rotations) = data @@ -187,14 +171,14 @@ function multiply_by_hamiltonian_dipole_aux!(y, x, phasebuf, qphase, swt) S = (N-1)/2 L = natoms(sys.crystal) - y .= 0 - nq = size(y, 1) - x = Base.ReshapedArray(x, (nq, L, 2), ()) - y = Base.ReshapedArray(y, (nq, L, 2), ()) + Nq = length(qs_reshaped) + @assert size(x) == size(y) == (Nq, 2L) + X = reshape(x, (Nq, L, 2)) + Y = reshape(y, (Nq, L, 2)) + Y .= 0 - # Add single-site terms, both Zeeman and single-ion anisotropy. These - # entries are q-independent and can be precomputed, but there appears to be - # little or no advantage to this. + # Add Zeeman and single-ion anisotropy. These entries are q-independent and + # could be precomputed. (; units, extfield, gs) = sys for i in 1:L (; c2, c4, c6) = stevens_coefs[i] @@ -206,25 +190,25 @@ function multiply_by_hamiltonian_dipole_aux!(y, x, phasebuf, qphase, swt) # Seems to be no benefit to breaking this into two loops acting on # different final indices, presumably because memory access patterns for - # x and y cannot be simultaneously optimized. - @inbounds for q in 1:nq - y[q, i, 1] += (B′ + A1) * x[q, i, 1] + A2 * x[q, i, 2] - y[q, i, 2] += (B′ + A1) * x[q, i, 2] + conj(A2) * x[q, i, 1] + # X and Y cannot be simultaneously optimized. + @inbounds for q in 1:Nq + Y[q, i, 1] += (B′ + A1) * X[q, i, 1] + A2 * X[q, i, 2] + Y[q, i, 2] += (B′ + A1) * X[q, i, 2] + conj(A2) * X[q, i, 1] end end - # Add pairwise terms + # Pair interactions for ints in sys.interactions_union # Bilinear exchange for coupling in ints.pair (; isculled, bond) = coupling isculled && break - i, j = bond.i, bond.j + (; i, j) = bond - # Calculate phase associated with periodic wrapping - n1, n2, n3 = bond.n - map!(qp -> (qp[1]^n1)*(qp[2]^n2)*(qp[3]^n3), phasebuf, qphase) + map!(phases, qs_reshaped) do q + cis(2π*dot(q, bond.n)) + end if !iszero(coupling.bilin) J = coupling.bilin # This is Rij in previous notation (transformed exchange matrix) @@ -232,25 +216,25 @@ function multiply_by_hamiltonian_dipole_aux!(y, x, phasebuf, qphase, swt) 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]) - @inbounds for q in axes(y, 1) - y[q, i, 1] += Q * phasebuf[q] * x[q, j, 1] - y[q, i, 1] += conj(P) * phasebuf[q] * x[q, j, 2] - y[q, i, 1] -= 0.5 * J[3, 3] * x[q, i, 1] + @inbounds for q in axes(Y, 1) + Y[q, i, 1] += Q * phases[q] * X[q, j, 1] + Y[q, i, 1] += conj(P) * phases[q] * X[q, j, 2] + Y[q, i, 1] -= 0.5 * J[3, 3] * X[q, i, 1] end - @inbounds for q in axes(y, 1) - y[q, i, 2] += conj(Q) * phasebuf[q] * x[q, j, 2] - y[q, i, 2] += P * phasebuf[q] * x[q, j, 1] - y[q, i, 2] -= 0.5 * J[3, 3] * x[q, i, 2] + @inbounds for q in axes(Y, 1) + Y[q, i, 2] += conj(Q) * phases[q] * X[q, j, 2] + Y[q, i, 2] += P * phases[q] * X[q, j, 1] + Y[q, i, 2] -= 0.5 * J[3, 3] * X[q, i, 2] end - @inbounds for q in axes(y, 1) - y[q, j, 1] += conj(P) * conj(phasebuf[q]) * x[q, i, 2] - y[q, j, 1] += conj(Q) * conj(phasebuf[q]) * x[q, i, 1] - y[q, j, 1] -= 0.5 * J[3, 3] * x[q, j, 1] + @inbounds for q in axes(Y, 1) + Y[q, j, 1] += conj(P) * conj(phases[q]) * X[q, i, 2] + Y[q, j, 1] += conj(Q) * conj(phases[q]) * X[q, i, 1] + Y[q, j, 1] -= 0.5 * J[3, 3] * X[q, j, 1] end - @inbounds for q in axes(y, 1) - y[q, j, 2] += Q * conj(phasebuf[q]) * x[q, i, 2] - y[q, j, 2] += P * conj(phasebuf[q]) * x[q, i, 1] - y[q, j, 2] -= 0.5 * J[3, 3] * x[q, j, 2] + @inbounds for q in axes(Y, 1) + Y[q, j, 2] += Q * conj(phases[q]) * X[q, i, 2] + Y[q, j, 2] += P * conj(phases[q]) * X[q, i, 1] + Y[q, j, 2] -= 0.5 * J[3, 3] * X[q, j, 2] end end @@ -261,40 +245,40 @@ function multiply_by_hamiltonian_dipole_aux!(y, x, phasebuf, qphase, swt) 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])) - @inbounds for q in 1:nq - y[q, i, 1] += -6J[3, 3] * x[q, i, 1] - y[q, i, 1] += 12*(J[1, 3] + im*J[5, 3]) * x[q, i, 2] - y[q, i, 1] += Q * phasebuf[q] * x[q, j, 1] - y[q, i, 1] += conj(P) * phasebuf[q] * x[q, j, 2] + @inbounds for q in 1:Nq + Y[q, i, 1] += -6J[3, 3] * X[q, i, 1] + Y[q, i, 1] += 12*(J[1, 3] + im*J[5, 3]) * X[q, i, 2] + Y[q, i, 1] += Q * phases[q] * X[q, j, 1] + Y[q, i, 1] += conj(P) * phases[q] * X[q, j, 2] end - @inbounds for q in 1:nq - y[q, i, 2] += -6J[3, 3] * x[q, i, 2] - y[q, i, 2] += 12*(J[1, 3] - im*J[5, 3]) * x[q, i, 1] - y[q, i, 2] += conj(Q) * phasebuf[q] * x[q, j, 2] - y[q, i, 2] += P * phasebuf[q] * x[q, j, 1] + @inbounds for q in 1:Nq + Y[q, i, 2] += -6J[3, 3] * X[q, i, 2] + Y[q, i, 2] += 12*(J[1, 3] - im*J[5, 3]) * X[q, i, 1] + Y[q, i, 2] += conj(Q) * phases[q] * X[q, j, 2] + Y[q, i, 2] += P * phases[q] * X[q, j, 1] end - @inbounds for q in 1:nq - y[q, j, 1] += -6J[3, 3] * x[q, j, 1] - y[q, j, 1] += 12*(J[3, 1] + im*J[3, 5]) * x[q, j, 2] - y[q, j, 1] += conj(Q) * conj(phasebuf[q]) * x[q, i, 1] - y[q, j, 1] += conj(P) * conj(phasebuf[q]) * x[q, i, 2] + @inbounds for q in 1:Nq + Y[q, j, 1] += -6J[3, 3] * X[q, j, 1] + Y[q, j, 1] += 12*(J[3, 1] + im*J[3, 5]) * X[q, j, 2] + Y[q, j, 1] += conj(Q) * conj(phases[q]) * X[q, i, 1] + Y[q, j, 1] += conj(P) * conj(phases[q]) * X[q, i, 2] end - @inbounds for q in 1:nq - y[q, j, 2] += -6J[3, 3] * x[q, j, 2] - y[q, j, 2] += 12*(J[3, 1] - im*J[3, 5]) * x[q, j, 1] - y[q, j, 2] += Q * conj(phasebuf[q]) * x[q, i, 2] - y[q, j, 2] += P * conj(phasebuf[q]) * x[q, i, 1] + @inbounds for q in 1:Nq + Y[q, j, 2] += -6J[3, 3] * X[q, j, 2] + Y[q, j, 2] += 12*(J[3, 1] - im*J[3, 5]) * X[q, j, 1] + Y[q, j, 2] += Q * conj(phases[q]) * X[q, i, 2] + Y[q, j, 2] += P * conj(phases[q]) * X[q, i, 1] end end - - if !isnothing(sys.ewald) - error("Ewald not supported") - end end end + if !isnothing(sys.ewald) + error("Ewald not supported") + end + # Add small constant shift for positive-definiteness. - @inbounds @. y += swt.energy_ϵ * x + @inbounds @. Y += swt.energy_ϵ * X nothing end diff --git a/src/SpinWaveTheory/HamiltonianSUN.jl b/src/SpinWaveTheory/HamiltonianSUN.jl index 4b20578b0..00fe7645a 100644 --- a/src/SpinWaveTheory/HamiltonianSUN.jl +++ b/src/SpinWaveTheory/HamiltonianSUN.jl @@ -1,91 +1,67 @@ -# Construct portion of Hamiltonian due to onsite terms (single-site anisotropy -# or external field). -function swt_onsite_coupling!(H, op, swt, atom) - sys = swt.sys - N = sys.Ns[1] - nflavors = N - 1 - L = nflavors * natoms(sys.crystal) - newdims = (nflavors, natoms(sys.crystal), nflavors, natoms(sys.crystal)) - - H11 = reshape(view(H, 1:L, 1:L), newdims) - H22 = reshape(view(H, L+1:2L, L+1:2L), newdims) - - for m in 1:N-1 - for n in 1:N-1 - c = 0.5 * (op[m, n] - δ(m, n) * op[N, N]) - H11[m, atom, n, atom] += c - H22[n, atom, m, atom] += c - end - end -end - -# Adds contribution of a bilinear pair of operators Ai * Bj, where Ai is on site -# i, and Bj is on site j. Typically these will be generated by the sparse tensor -# decomposition stored in the `general` field of a `PairCoupling`. -function swt_pair_coupling!(H, Ai, Bj, swt, phase, bond) - (; i, j) = bond - sys = swt.sys - N = sys.Ns[1] - nflavors = N - 1 - L = nflavors * natoms(sys.crystal) - newdims = (nflavors, natoms(sys.crystal), nflavors, natoms(sys.crystal)) - - H11 = reshape(view(H, 1:L, 1:L), newdims) - H12 = reshape(view(H, 1:L, L+1:2L), newdims) - H21 = reshape(view(H, L+1:2L, 1:L), newdims) - H22 = reshape(view(H, L+1:2L, L+1:2L), newdims) - - for m in 1:N-1 - for n in 1:N-1 - c = 0.5 * (Ai[m,n] - δ(m,n)*Ai[N,N]) * (Bj[N,N]) - H11[m, i, n, i] += c - H22[n, i, m, i] += c - - c = 0.5 * Ai[N,N] * (Bj[m,n] - δ(m,n)*Bj[N,N]) - H11[m, j, n, j] += c - H22[n, j, m, j] += c - - c = 0.5 * Ai[m,N] * Bj[N,n] - H11[m, i, n, j] += c * phase - H22[n, j, m, i] += c * conj(phase) - - c = 0.5 * Ai[N,m] * Bj[n,N] - H11[n, j, m, i] += c * conj(phase) - H22[m, i, n, j] += c * phase - - c = 0.5 * Ai[m,N] * Bj[n,N] - H12[m, i, n, j] += c * phase - H12[n, j, m, i] += c * conj(phase) - H21[n, j, m, i] += conj(c) * conj(phase) - H21[m, i, n, j] += conj(c) * phase - end - end -end - # Set the dynamical quadratic Hamiltonian matrix in SU(N) mode. function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_reshaped::Vec3) (; sys) = swt - L = nbands(swt) # Number of quasiparticle bands - @assert size(H) == (2L, 2L) + N = sys.Ns[1] + Na = natoms(sys.crystal) + L = (N-1) * Na # Clear the Hamiltonian + @assert size(H) == (2L, 2L) H .= 0 + blockdims = (N-1, Na, N-1, Na) + H11 = reshape(view(H, 1:L, 1:L), blockdims) + H12 = reshape(view(H, 1:L, L+1:2L), blockdims) + H21 = reshape(view(H, L+1:2L, 1:L), blockdims) + H22 = reshape(view(H, L+1:2L, L+1:2L), blockdims) + + for (i, int) in enumerate(sys.interactions_union) + + # Onsite coupling, including Zeeman. Note that op has already been + # transformed according to the local frame of sublattice i. + op = int.onsite + for m in 1:N-1 + for n in 1:N-1 + c = 0.5 * (op[m, n] - δ(m, n) * op[N, N]) + H11[m, i, n, i] += c + H22[n, i, m, i] += c + end + end - for (atom, int) in enumerate(sys.interactions_union) - # Set the onsite term for atom - swt_onsite_coupling!(H, int.onsite, swt, atom) - - # Set pair interactions that originate for atom for coupling in int.pair - # Extract information common to bond (; isculled, bond) = coupling isculled && break - + (; i, j) = bond phase = exp(2π*im * dot(q_reshaped, bond.n)) # Phase associated with periodic wrapping - for (A, B) in coupling.general.data - swt_pair_coupling!(H, A, B, swt, phase, bond) + + # Set "general" pair interactions of the form Aᵢ⊗Bⱼ. Note that Aᵢ + # and Bᵢ have already been transformed according to the local frames + # of sublattice i and j, respectively. + for (Ai, Bj) in coupling.general.data + for m in 1:N-1, n in 1:N-1 + c = 0.5 * (Ai[m,n] - δ(m,n)*Ai[N,N]) * (Bj[N,N]) + H11[m, i, n, i] += c + H22[n, i, m, i] += c + + c = 0.5 * Ai[N,N] * (Bj[m,n] - δ(m,n)*Bj[N,N]) + H11[m, j, n, j] += c + H22[n, j, m, j] += c + + c = 0.5 * Ai[m,N] * Bj[N,n] + H11[m, i, n, j] += c * phase + H22[n, j, m, i] += c * conj(phase) + + c = 0.5 * Ai[N,m] * Bj[n,N] + H11[n, j, m, i] += c * conj(phase) + H22[m, i, n, j] += c * phase + + c = 0.5 * Ai[m,N] * Bj[n,N] + H12[m, i, n, j] += c * phase + H12[n, j, m, i] += c * conj(phase) + H21[n, j, m, i] += conj(c) * conj(phase) + H21[m, i, n, j] += conj(c) * phase + end end end end @@ -107,117 +83,93 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh 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!(y, x, op, swt, atom) - sys = swt.sys - N = sys.Ns[1] - nflavors = N - 1 - - nq = size(y, 1) - x = Base.ReshapedArray(x, (nq, nflavors, natoms(sys.crystal), 2), ()) - y = Base.ReshapedArray(y, (nq, 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]) - @inbounds for q in 1:nq - y[q, m, atom, 1] += c * x[q, n, atom, 1] - y[q, n, atom, 2] += c * x[q, m, atom, 2] - end - 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!(y, x, Ti, Tj, swt, phase, bond) - (; i, j) = bond - sys = swt.sys - N = sys.Ns[1] - nflavors = N - 1 - - nq = size(y, 1) - x = Base.ReshapedArray(x, (nq, nflavors, natoms(sys.crystal), 2), ()) - y = Base.ReshapedArray(y, (nq, nflavors, natoms(sys.crystal), 2), ()) - - for m in 1:N-1 - for n in 1:N-1 - c1 = 0.5 * (Ti[m,n] - δ(m,n)*Ti[N,N]) * Tj[N,N] - c2 = 0.5 * Ti[N,N] * (Tj[m,n] - δ(m,n)*Tj[N,N]) - c3 = 0.5 * Ti[m,N] * Tj[N,n] - c4 = 0.5 * Ti[N,m] * Tj[n,N] - c5 = 0.5 * Ti[m,N] * Tj[n,N] - - @inbounds for q in axes(y, 1) - y[q, m, i, 1] += c1 * x[q, n, i, 1] - y[q, n, i, 2] += c1 * x[q, m, i, 2] - - y[q, m, j, 1] += c2 * x[q, n, j, 1] - y[q, n, j, 2] += c2 * x[q, m, j, 2] - - y[q, m, i, 1] += c3 * phase[q] * x[q, n, j, 1] - y[q, n, j, 2] += c3 * conj(phase[q]) * x[q, m, i, 2] - - y[q, n, j, 1] += c4 * conj(phase[q]) * x[q, m, i, 1] - y[q, m, i, 2] += c4 * phase[q] * x[q, n, j, 2] - - y[q, m, i, 1] += c5 * phase[q] * x[q, n, j, 2] - y[q, n, j, 1] += c5 * conj(phase[q]) * x[q, m, i, 2] - y[q, m, i, 2] += conj(c5 * phase[q]) * x[q, n, j, 1] - y[q, n, j, 2] += conj(c5) * phase[q] * x[q, m, i, 1] - end - end - end +function multiply_by_hamiltonian_SUN(x::Array{ComplexF64, 2}, swt::SpinWaveTheory, qs_reshaped::Array{Vec3}) + y = zero(x) + multiply_by_hamiltonian_SUN!(y, x, swt, qs_reshaped) + return y end # Calculate y = H*x, where H is the quadratic Hamiltonian matrix (dynamical # matrix). Note that x is assumed to be a 2D array with first index # corresponding to q. -function multiply_by_hamiltonian_SUN(x::Array{ComplexF64, 2}, swt::SpinWaveTheory, qs_reshaped::Array{Vec3}) - # Preallocate buffers - L = nbands(swt) - @assert size(x, 2) == 2L - y = zeros(ComplexF64, (size(qs_reshaped)..., 2L)) - phasebuf = zeros(ComplexF64, length(qs_reshaped)) - - # Precompute e^{2πq_α} components - qphase = map(qs_reshaped) do q - (exp(2π*im*q[1]), exp(2π*im*q[2]), exp(2π*im*q[3])) - end - - # Perform batched matrix-vector multiply - multiply_by_hamiltonian_SUN_aux!(reshape(y, (length(qs_reshaped), 2L)), x, phasebuf, qphase, swt) - - return y -end - -function multiply_by_hamiltonian_SUN_aux!(y, x, phasebuf, qphase, swt) +function multiply_by_hamiltonian_SUN!(y::Array{ComplexF64, 2}, x::Array{ComplexF64, 2}, swt::SpinWaveTheory, qs_reshaped::Array{Vec3}; + phases=zeros(ComplexF64, size(qs_reshaped))) (; sys) = swt - y .= 0 - - # 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) + N = sys.Ns[1] + Na = natoms(sys.crystal) + L = (N-1) * Na + + Nq = length(qs_reshaped) + @assert size(x) == size(y) == (Nq, 2L) + X = reshape(x, (Nq, N-1, Na, 2)) + Y = reshape(y, (Nq, N-1, Na, 2)) + Y .= 0 + + # All operators appearing in interactions have been pre-rotated to local + # frame + for (i, int) in enumerate(sys.interactions_union) + + # Onsite coupling, including Zeeman + op = int.onsite + for m in 1:N-1, n in 1:N-1 + c = 0.5 * (op[m, n] - δ(m, n) * op[N, N]) + @inbounds for q in 1:Nq + Y[q, m, i, 1] += c * X[q, n, i, 1] + Y[q, n, i, 2] += c * X[q, m, i, 2] + end + end + # General pair interactions for coupling in int.pair # Extract information common to bond (; isculled, bond) = coupling isculled && break + (; i, j) = bond - # phase = exp(2π*im * dot(q_reshaped, bond.n)) # Phase associated with periodic wrapping - n1, n2, n3 = bond.n - map!(qp -> (qp[1]^n1)*(qp[2]^n2)*(qp[3]^n3), phasebuf, qphase) - for (A, B) in coupling.general.data - multiply_by_pair_coupling_SUN!(y, x, A, B, swt, phasebuf, bond) + map!(phases, qs_reshaped) do q + cis(2π*dot(q, bond.n)) + end + + # General pair interactions + for (A, B) in coupling.general.data + + for m in 1:N-1, n in 1:N-1 + c1 = 0.5 * (A[m,n] - δ(m,n)*A[N,N]) * B[N,N] + c2 = 0.5 * A[N,N] * (B[m,n] - δ(m,n)*B[N,N]) + c3 = 0.5 * A[m,N] * B[N,n] + c4 = 0.5 * A[N,m] * B[n,N] + c5 = 0.5 * A[m,N] * B[n,N] + + @inbounds for q in axes(Y, 1) + Y[q, m, i, 1] += c1 * X[q, n, i, 1] + Y[q, n, i, 2] += c1 * X[q, m, i, 2] + + Y[q, m, j, 1] += c2 * X[q, n, j, 1] + Y[q, n, j, 2] += c2 * X[q, m, j, 2] + + Y[q, m, i, 1] += c3 * phases[q] * X[q, n, j, 1] + Y[q, n, j, 2] += c3 * conj(phases[q]) * X[q, m, i, 2] + + Y[q, n, j, 1] += c4 * conj(phases[q]) * X[q, m, i, 1] + Y[q, m, i, 2] += c4 * phases[q] * X[q, n, j, 2] + + Y[q, m, i, 1] += c5 * phases[q] * X[q, n, j, 2] + Y[q, n, j, 1] += c5 * conj(phases[q]) * X[q, m, i, 2] + Y[q, m, i, 2] += conj(c5 * phases[q]) * X[q, n, j, 1] + Y[q, n, j, 2] += conj(c5) * phases[q] * X[q, m, i, 1] + end + end end end end + if !isnothing(sys.ewald) + error("Ewald not supported") + end + # Add small constant shift for positive-definiteness - @inbounds @. y += swt.energy_ϵ * x + @inbounds @. Y += swt.energy_ϵ * X nothing -end \ No newline at end of file +end diff --git a/test/test_lswt.jl b/test/test_lswt.jl index 3b0bb2277..0e4685e1d 100644 --- a/test/test_lswt.jl +++ b/test/test_lswt.jl @@ -138,9 +138,8 @@ end # Crystal a = 8.289 latvecs = lattice_vectors(a, a, a, 90, 90, 90) - types = ["MzR1"] positions = [[0, 0, 0]] - fcc = Crystal(latvecs, positions, 225; types) + fcc = Crystal(latvecs, positions, 225) S = 5/2 g = 2 @@ -332,13 +331,25 @@ end end -#= -@testitem "Dipole-dipole unimplemented" begin - sys = System(Sunny.diamond_crystal(),(1,1,1),[SpinInfo(1,S=1/2,g=2)],:SUN;seed = 0) +@testitem "Dipole-dipole" begin + latvecs = lattice_vectors(10, 10, 1, 90, 90, 90) + cryst = Crystal(latvecs, [[0,0,0]]) + sys = System(cryst, (1,1,1), [SpinInfo(1; S=1, g=1)], :dipole; units=Units.theory) enable_dipole_dipole!(sys) - @test_throws "SpinWaveTheory does not yet support long-range dipole-dipole interactions." SpinWaveTheory(sys) + + randomize_spins!(sys) + minimize_energy!(sys) # Polarized chain + @test energy_per_site(sys) ≈ -0.1913132980155851 + + swt = SpinWaveTheory(sys) + formula = intensity_formula(swt, :perp; kernel=delta_function_kernel) + + qpoints = [[0, 0, 0], [0, 0, 1/2], [0, 1/2, 1/2], [0, 0, 0]] + disps, is = intensities_bands(swt, qpoints, formula) + @test disps ≈ [0.5689399140467553, 0.23914164251944922, 0.23914164251948083, 0.5689399140467553] + @test is ≈ [1, 1, 201/202, 1] end -=# + @testitem "SW15-Langasite" begin # Ba3NbFe3Si2O14 From 543c5e2697c46ff350e5277d12b11b30ec33fc8b Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 11/29] SU(N) Ewald LSWT --- .../DispersionAndIntensities.jl | 6 +-- src/SpinWaveTheory/HamiltonianSUN.jl | 52 ++++++++++++++++++- src/SpinWaveTheory/SpinWaveTheory.jl | 14 +++-- test/test_lswt.jl | 28 +++++----- 4 files changed, 79 insertions(+), 21 deletions(-) diff --git a/src/SpinWaveTheory/DispersionAndIntensities.jl b/src/SpinWaveTheory/DispersionAndIntensities.jl index 821045a29..c979e38a5 100644 --- a/src/SpinWaveTheory/DispersionAndIntensities.jl +++ b/src/SpinWaveTheory/DispersionAndIntensities.jl @@ -236,7 +236,7 @@ The integral of a properly normalized kernel function over all `Δω` is one. function intensity_formula(f::Function, swt::SpinWaveTheory, corr_ix::AbstractVector{Int64}; kernel::Union{Nothing,Function}, return_type=Float64, string_formula="f(Q,ω,S{α,β}[ix_q,ix_ω])", formfactors=nothing) (; sys, data, observables) = swt - (; observable_operators) = data + (; observables_localized) = data # Number of atoms in magnetic cell Nm = length(sys.dipoles) @@ -326,7 +326,7 @@ function intensity_formula(f::Function, swt::SpinWaveTheory, corr_ix::AbstractVe v = reshape(view(V, :, band), N-1, Nm, 2) for i = 1:Nm for μ = 1:num_observables(observables) - @views O = observable_operators[:, :, μ, i] + @views O = observables_localized[:, :, μ, i] for α = 1:N-1 Avec[μ] += Avec_pref[i] * (O[α, N] * v[α, i, 2] + O[N, α] * v[α, i, 1]) end @@ -349,7 +349,7 @@ function intensity_formula(f::Function, swt::SpinWaveTheory, corr_ix::AbstractVe # (; local_rotations) = data # displacement_global_frame = local_rotations[i] * displacement_local_frame - @views O_local_frame = observable_operators[:,:,μ,i] + @views O_local_frame = observables_localized[:,:,μ,i] Avec[μ] += Avec_pref[i] * sqrt_halfS * (O_local_frame * displacement_local_frame)[1] end end diff --git a/src/SpinWaveTheory/HamiltonianSUN.jl b/src/SpinWaveTheory/HamiltonianSUN.jl index 00fe7645a..4e5447814 100644 --- a/src/SpinWaveTheory/HamiltonianSUN.jl +++ b/src/SpinWaveTheory/HamiltonianSUN.jl @@ -1,7 +1,9 @@ # Set the dynamical quadratic Hamiltonian matrix in SU(N) mode. function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_reshaped::Vec3) - (; sys) = swt + (; sys, data) = swt + (; spins_localized) = data + (; gs, units) = sys N = sys.Ns[1] Na = natoms(sys.crystal) @@ -67,7 +69,53 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh end if !isnothing(sys.ewald) - error("Ewald not yet supported in LSWT for :SUN mode") + N = sys.Ns[1] + μB² = units.μB^2 + + # Interaction matrix for wavevector q + A = precompute_dipole_ewald_aux(sys.crystal, (1,1,1), units.μ0, q_reshaped, cis, Val{ComplexF64}()) + A = reshape(A, Na, Na) + + # Interaction matrix for wavevector (0,0,0) + A0 = sys.ewald.A + A0 = reshape(A0, Na, Na) + + for i in 1:Na, j in 1:Na + # An ordered pair of magnetic moments contribute (μᵢ A μⱼ)/2 to the + # energy. A symmetric contribution will appear for the bond reversal + # (i, j) → (j, i). Note that μ = -μB g S. + J = μB² * gs[i]' * A[i, j] * gs[j] / 2 + J0 = μB² * gs[i]' * A0[i, j] * gs[j] / 2 + + for α in 1:3, β in 1:3 + Ai = view(spins_localized, :, :, α, i) + Bj = view(spins_localized, :, :, β, j) + + for m in 1:N-1, n in 1:N-1 + c = 0.5 * (Ai[m,n] - δ(m,n)*Ai[N,N]) * (Bj[N,N]) + H11[m, i, n, i] += c * J0[α, β] + H22[n, i, m, i] += c * J0[α, β] + + c = 0.5 * Ai[N,N] * (Bj[m,n] - δ(m,n)*Bj[N,N]) + H11[m, j, n, j] += c * J0[α, β] + H22[n, j, m, j] += c * J0[α, β] + + c = 0.5 * Ai[m,N] * Bj[N,n] + H11[m, i, n, j] += c * J[α, β] + H22[n, j, m, i] += c * conj(J[α, β]) + + c = 0.5 * Ai[N,m] * Bj[n,N] + H11[n, j, m, i] += c * conj(J[α, β]) + H22[m, i, n, j] += c * J[α, β] + + c = 0.5 * Ai[m,N] * Bj[n,N] + H12[m, i, n, j] += c * J[α, β] + H12[n, j, m, i] += c * conj(J[α, β]) + H21[n, j, m, i] += conj(c) * conj(J[α, β]) + H21[m, i, n, j] += conj(c) * J[α, β] + end + end + end end # H must be hermitian up to round-off errors diff --git a/src/SpinWaveTheory/SpinWaveTheory.jl b/src/SpinWaveTheory/SpinWaveTheory.jl index cce60fbee..2b8a58ded 100644 --- a/src/SpinWaveTheory/SpinWaveTheory.jl +++ b/src/SpinWaveTheory/SpinWaveTheory.jl @@ -1,12 +1,13 @@ struct SWTDataDipole local_rotations :: Vector{Mat3} stevens_coefs :: Vector{StevensExpansion} - observable_operators :: Array{ComplexF64, 4} # Observables in local frame (for intensity calcs) + observables_localized :: Array{ComplexF64, 4} # Observables in local frame (for intensity calcs) end struct SWTDataSUN local_unitaries :: Array{ComplexF64, 3} # Aligns to quantization axis on each site - observable_operators :: Array{ComplexF64, 4} # Observables in local frame (for intensity calcs) + spins_localized :: Array{ComplexF64, 4} # Spins in local frame + observables_localized :: Array{ComplexF64, 4} # Observables in local frame (for intensity calcs) end """ @@ -118,6 +119,7 @@ function swt_data_sun(sys::System{N}, obs) where N # Preallocate buffers for local unitaries and observables. local_unitaries = zeros(ComplexF64, N, N, n_magnetic_atoms) + spins_localized = zeros(ComplexF64, N, N, 3, n_magnetic_atoms) observables_localized = zeros(ComplexF64, N, N, num_observables(obs), n_magnetic_atoms) for atom in 1:n_magnetic_atoms @@ -132,7 +134,12 @@ function swt_data_sun(sys::System{N}, obs) where N U = view(local_unitaries, :, :, atom) # Rotate observables into local reference frames - for k = 1:num_observables(obs) + S = spin_matrices_of_dim(; N) + for k in 1:3 + spins_localized[:, :, k, atom] = Hermitian(U' * S[k] * U) + end + + for k in 1:num_observables(obs) A = observable_at_site(obs.observables[k],CartesianIndex(1,1,1,atom)) observables_localized[:, :, k, atom] = Hermitian(U' * convert(Matrix, A) * U) end @@ -167,6 +174,7 @@ function swt_data_sun(sys::System{N}, obs) where N return SWTDataSUN( local_unitaries, + spins_localized, observables_localized ) end diff --git a/test/test_lswt.jl b/test/test_lswt.jl index 0e4685e1d..73ff870af 100644 --- a/test/test_lswt.jl +++ b/test/test_lswt.jl @@ -334,20 +334,22 @@ end @testitem "Dipole-dipole" begin latvecs = lattice_vectors(10, 10, 1, 90, 90, 90) cryst = Crystal(latvecs, [[0,0,0]]) - sys = System(cryst, (1,1,1), [SpinInfo(1; S=1, g=1)], :dipole; units=Units.theory) - enable_dipole_dipole!(sys) - randomize_spins!(sys) - minimize_energy!(sys) # Polarized chain - @test energy_per_site(sys) ≈ -0.1913132980155851 - - swt = SpinWaveTheory(sys) - formula = intensity_formula(swt, :perp; kernel=delta_function_kernel) - - qpoints = [[0, 0, 0], [0, 0, 1/2], [0, 1/2, 1/2], [0, 0, 0]] - disps, is = intensities_bands(swt, qpoints, formula) - @test disps ≈ [0.5689399140467553, 0.23914164251944922, 0.23914164251948083, 0.5689399140467553] - @test is ≈ [1, 1, 201/202, 1] + for mode in (:dipole, :SUN) + sys = System(cryst, (1,1,1), [SpinInfo(1; S=1, g=1)], mode; units=Units.theory) + enable_dipole_dipole!(sys) + + polarize_spins!(sys, (0,0,1)) + @test energy_per_site(sys) ≈ -0.1913132980155851 + + swt = SpinWaveTheory(sys) + formula = intensity_formula(swt, :perp; kernel=delta_function_kernel) + + qpoints = [[0, 0, 0], [0, 0, 1/2], [0, 1/2, 1/2], [0, 0, 0]] + disps, is = intensities_bands(swt, qpoints, formula) + @test disps[:,end] ≈ [0.5689399140467553, 0.23914164251944922, 0.23914164251948083, 0.5689399140467553] + @test is[:,end] ≈ [1, 1, 201/202, 1] + end end From f92f6fe9e47303f63998edf165e7de2825f8b3c6 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 12/29] Fix painful bug violating U(1) gauge symmetry --- ext/PlottingExt.jl | 2 -- src/SpinWaveTheory/HamiltonianDipole.jl | 18 +++++++++--------- test/test_ewald.jl | 2 +- 3 files changed, 10 insertions(+), 12 deletions(-) diff --git a/ext/PlottingExt.jl b/ext/PlottingExt.jl index 6a0e24bcb..2fd984b16 100644 --- a/ext/PlottingExt.jl +++ b/ext/PlottingExt.jl @@ -881,8 +881,6 @@ function plot_spins!(ax, sys::System; notifier=Makie.Observable(nothing), arrows sys.latsize[[2,3]] == [1,1] || error("System not one-dimensional in (a₁)") end - supervecs = sys.crystal.latvecs * diagm(Vec3(sys.latsize)) - # Show bounding box of magnetic supercell in gray (this needs to come first # to set a scale for the scene in case there is only one atom). supervecs = sys.crystal.latvecs * diagm(Vec3(sys.latsize)) diff --git a/src/SpinWaveTheory/HamiltonianDipole.jl b/src/SpinWaveTheory/HamiltonianDipole.jl index 3c0930955..b1d927992 100644 --- a/src/SpinWaveTheory/HamiltonianDipole.jl +++ b/src/SpinWaveTheory/HamiltonianDipole.jl @@ -44,14 +44,13 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r 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]) - + Q = 0.25 * (J[1, 1] + J[2, 2] - im*(J[1, 2] - J[2, 1])) H11[i, j] += Q * phase H11[j, i] += conj(Q) * conj(phase) H22[i, j] += conj(Q) * phase H22[j, i] += Q * conj(phase) + P = 0.25 * (J[1, 1] - J[2, 2] - im*(J[1, 2] + J[2, 1])) H21[i, j] += P * phase H21[j, i] += P * conj(phase) H12[i, j] += conj(P) * phase @@ -76,14 +75,13 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r H21[j, j] += 12*(J[3, 1] - im*J[3, 5]) H12[j, j] += 12*(J[3, 1] + im*J[3, 5]) - 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])) - H11[i, j] += Q * phase H11[j, i] += conj(Q) * conj(phase) H22[i, j] += conj(Q) * phase H22[j, i] += Q * conj(phase) + P = 0.25 * (-J[4, 4]+J[2, 2] - im*( J[4, 2]+J[2, 4])) H21[i, j] += P * phase H21[j, i] += P * conj(phase) H12[i, j] += conj(P) * phase @@ -110,6 +108,7 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r # Loop over sublattice pairs for i in 1:L, j in 1:L + # An ordered pair of magnetic moments contribute (μᵢ A μⱼ)/2 to the # energy. A symmetric contribution will appear for the bond reversal # (i, j) → (j, i). Note that μ = -μB g S. @@ -124,16 +123,17 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r # Interactions for Jˣˣ, Jʸʸ, Jˣʸ, and Jʸˣ at wavevector q. Q⁻ = 0.25 * (J[1, 1] + J[2, 2] - im*(J[1, 2] - J[2, 1])) Q⁺ = 0.25 * (J[1, 1] + J[2, 2] + im*(J[1, 2] - J[2, 1])) - P⁻ = 0.25 * (J[1, 1] - J[2, 2] - im*(J[1, 2] + J[2, 1])) - P⁺ = 0.25 * (J[1, 1] - J[2, 2] + im*(J[1, 2] + J[2, 1])) H11[i, j] += Q⁻ H11[j, i] += conj(Q⁻) H22[i, j] += Q⁺ H22[j, i] += conj(Q⁺) + + P⁻ = 0.25 * (J[1, 1] - J[2, 2] - im*(J[1, 2] + J[2, 1])) + P⁺ = 0.25 * (J[1, 1] - J[2, 2] + im*(J[1, 2] + J[2, 1])) H21[i, j] += P⁻ + H21[j, i] += conj(P⁺) + H12[i, j] += P⁺ H12[j, i] += conj(P⁻) - H21[j, i] += P⁺ - H12[i, j] += conj(P⁺) # Interactions for Jᶻᶻ at wavevector (0,0,0). H11[i, i] -= 0.5 * J0[3, 3] diff --git a/test/test_ewald.jl b/test/test_ewald.jl index ca30dd9b6..420262122 100644 --- a/test/test_ewald.jl +++ b/test/test_ewald.jl @@ -26,7 +26,7 @@ # Same thing, with multiple unit cells sys = System(cryst, (2,3,4), infos, :dipole; units=Units.theory) enable_dipole_dipole!(sys) - @test isapprox(energy(sys), -(1/6)prod(sys.latsize); atol=1e-13) + @test isapprox(energy_per_site(sys), -1/6; atol=1e-13) # Create a random box latvecs = lattice_vectors(1.1,0.9,0.8,92,85,95) From 8420b38f8a98e44f5fc21a2bff9e3bff28066e09 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 13/29] Another ewald test --- test/test_lswt.jl | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/test/test_lswt.jl b/test/test_lswt.jl index 73ff870af..23574b470 100644 --- a/test/test_lswt.jl +++ b/test/test_lswt.jl @@ -350,6 +350,26 @@ end @test disps[:,end] ≈ [0.5689399140467553, 0.23914164251944922, 0.23914164251948083, 0.5689399140467553] @test is[:,end] ≈ [1, 1, 201/202, 1] end + + begin + cryst = Sunny.bcc_crystal() + sys = System(cryst, (1, 1, 1), [SpinInfo(1, S=1, g=2)], :dipole, seed=2) + enable_dipole_dipole!(sys) + polarize_spins!(sys, (1,2,3)) # arbitrary direction + + R = hcat([1,1,-1], [-1,1,1], [1,-1,1]) / 2 + sys_reshape = reshape_supercell(sys, R) + @test energy_per_site(sys_reshape) ≈ energy_per_site(sys) ≈ -0.89944235377 + + swt1 = SpinWaveTheory(sys, energy_ϵ=1e-8) + swt2 = SpinWaveTheory(sys_reshape, energy_ϵ=1e-8) + path = [[0.5, -0.1, 0.3]] + disp1 = dispersion(swt1, path) + disp2 = dispersion(swt2, path) + + @test disp1 ≈ [1.3236778213351734 0.9206655419611791] + @test disp2 ≈ [0.9206655419611772] + end end From d45b992d21e64588a8b904c71b2f3ebb17235da7 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 14/29] Update docs and add tutorial --- ...bosons.lyx => bogoliubov_perturbation.lyx} | 0 docs/lyx_notes/fourier_lswt.lyx | 1086 +++++++++++++++++ docs/lyx_notes/phase_averaging.lyx | 459 ------- examples/07_Dipole_Dipole.jl | 63 + 4 files changed, 1149 insertions(+), 459 deletions(-) rename docs/lyx_notes/{bosons.lyx => bogoliubov_perturbation.lyx} (100%) create mode 100644 docs/lyx_notes/fourier_lswt.lyx delete mode 100644 docs/lyx_notes/phase_averaging.lyx create mode 100644 examples/07_Dipole_Dipole.jl diff --git a/docs/lyx_notes/bosons.lyx b/docs/lyx_notes/bogoliubov_perturbation.lyx similarity index 100% rename from docs/lyx_notes/bosons.lyx rename to docs/lyx_notes/bogoliubov_perturbation.lyx diff --git a/docs/lyx_notes/fourier_lswt.lyx b/docs/lyx_notes/fourier_lswt.lyx new file mode 100644 index 000000000..fe09d7c3e --- /dev/null +++ b/docs/lyx_notes/fourier_lswt.lyx @@ -0,0 +1,1086 @@ +#LyX 2.3 created this file. For more info see http://www.lyx.org/ +\lyxformat 544 +\begin_document +\begin_header +\save_transient_properties true +\origin unavailable +\textclass article +\use_default_options true +\maintain_unincluded_children false +\language english +\language_package default +\inputencoding auto +\fontencoding global +\font_roman "default" "default" +\font_sans "default" "default" +\font_typewriter "default" "default" +\font_math "auto" "auto" +\font_default_family default +\use_non_tex_fonts false +\font_sc false +\font_osf false +\font_sf_scale 100 100 +\font_tt_scale 100 100 +\use_microtype false +\use_dash_ligatures true +\graphics default +\default_output_format default +\output_sync 0 +\bibtex_command default +\index_command default +\paperfontsize default +\use_hyperref false +\papersize default +\use_geometry false +\use_package amsmath 1 +\use_package amssymb 1 +\use_package cancel 1 +\use_package esint 1 +\use_package mathdots 1 +\use_package mathtools 1 +\use_package mhchem 1 +\use_package stackrel 1 +\use_package stmaryrd 1 +\use_package undertilde 1 +\cite_engine basic +\cite_engine_type default +\use_bibtopic false +\use_indices false +\paperorientation portrait +\suppress_date false +\justification true +\use_refstyle 1 +\use_minted 0 +\index Index +\shortcut idx +\color #008000 +\end_index +\secnumdepth 3 +\tocdepth 3 +\paragraph_separation indent +\paragraph_indentation default +\is_math_indent 0 +\math_numbering_side default +\quotes_style english +\dynamic_quotes 0 +\papercolumns 1 +\papersides 1 +\paperpagestyle default +\tracking_changes false +\output_changes false +\html_math_output 0 +\html_css_as_file 0 +\html_be_strict false +\end_header + +\begin_body + +\begin_layout Title +Fourier space exchange interactions for spin wave theory +\end_layout + +\begin_layout Author +Kipton Barros +\end_layout + +\begin_layout Section +Classical bilinear interactions expressed in Fourier space +\end_layout + +\begin_layout Standard +Exchange energy for one unit cell is a sum over all couplings +\begin_inset Formula $c$ +\end_inset + + contained in the interaction list, +\begin_inset Formula +\begin{equation} +E=\frac{1}{2}\sum_{\mathbf{n}}\sum_{c}\delta_{c=(ij\mathbf{n})}J_{\alpha\beta}^{(c)}S_{\alpha i\mathbf{0}}S_{\beta j\mathbf{n}}.\label{eq:E} +\end{equation} + +\end_inset + +In what follows, there is implicit summation over any repeated, unbound + indices – this includes spin indices +\begin_inset Formula $(\alpha\beta)$ +\end_inset + + and sublattice indices +\begin_inset Formula $(ij$ +\end_inset + +). + The three indices in +\begin_inset Formula $\mathbf{n}$ +\end_inset + + denote a cell offset. + The notation +\begin_inset Formula $\delta_{c=(ij\mathbf{n})}$ +\end_inset + + is the indicator function – it is one if the +\begin_inset Formula $c$ +\end_inset + +th coupling is along a bond from atom +\begin_inset Formula $i$ +\end_inset + + to +\begin_inset Formula $j$ +\end_inset + +, displaced by +\begin_inset Formula $\mathbf{n}$ +\end_inset + + unit cells; otherwise, it is zero. + Given the supercell vectors +\begin_inset Formula $\mathbf{A}_{\{1,2,3\}}$ +\end_inset + +, the real-space displacement corresponding to +\begin_inset Formula $\mathbf{n}$ +\end_inset + + is +\begin_inset Formula +\begin{equation} +\mathbf{r}_{\mathbf{n}}=n_{1}\mathbf{A}_{1}+n_{2}\mathbf{A}_{2}+n_{3}\mathbf{A}_{3}, +\end{equation} + +\end_inset + + +\end_layout + +\begin_layout Standard +Define Fourier transforms from cells +\begin_inset Formula $\mathbf{n}\in\mathbb{Z}^{3}$ +\end_inset + + into momenta +\begin_inset Formula $\mathbf{q}$ +\end_inset + + defined in reciprocal lattice units (i.e., periodic in the cube +\begin_inset Formula $[0,1]^{3}$ +\end_inset + +), +\begin_inset Formula +\begin{align} +S_{\alpha i}(\mathbf{q}) & =\sum_{\mathbf{n}}e^{2\pi i\mathbf{q}\cdot\mathbf{n}}S_{\alpha i\mathbf{n}}\\ +J_{(\alpha i),(\beta j)}(\mathbf{q}) & =\sum_{\mathbf{n}}e^{2\pi i\mathbf{q}\cdot\mathbf{n}}\sum_{c}\delta_{c=(ij\mathbf{n})}J_{\alpha\beta}^{(c)}.\label{eq:J_q} +\end{align} + +\end_inset + +By construction, +\begin_inset Formula $S_{\alpha i}(\mathbf{q})$ +\end_inset + + is extensive, and scales like the number of cells in the macroscopic volume, +\end_layout + +\begin_layout Standard +\begin_inset Formula +\begin{equation} +N_{\mathrm{cells}}=\sum_{\mathbf{n}}1. +\end{equation} + +\end_inset + +Since +\begin_inset Formula $J_{\alpha\beta}^{(c)}$ +\end_inset + + is real, taking the complex-conjugate of both sides yields +\begin_inset Formula +\begin{equation} +J_{(\alpha i),(\beta j)}^{\ast}(\mathbf{q})=J_{(\alpha i),(\beta j)}(-\mathbf{q})=\sum_{\mathbf{n}}e^{2\pi i\mathbf{q}\cdot\mathbf{n}}\delta_{c=(ij,-\mathbf{n})}J_{\alpha\beta}^{(c)}. +\end{equation} + +\end_inset + +Without loss of generality, we can assume the pair interactions +\begin_inset Formula $\delta_{c=(ij\mathbf{n})}J_{\alpha\beta}^{(c)}$ +\end_inset + + and +\begin_inset Formula $\delta_{c=(ji,-\mathbf{n})}J_{\beta\alpha}^{(c)}$ +\end_inset + + contribute equal parts to the energy in Eq. +\begin_inset space ~ +\end_inset + + +\begin_inset CommandInset ref +LatexCommand eqref +reference "eq:E" +plural "false" +caps "false" +noprefix "false" + +\end_inset + +. + This allows to effectively swap indices +\begin_inset Formula $(\alpha i)\longleftrightarrow(\beta j)$ +\end_inset + + and +\series bold + +\begin_inset Formula $\mathbf{n}\longleftrightarrow-\mathbf{n}$ +\end_inset + + +\series default +. + The result is, +\begin_inset Formula +\begin{equation} +J_{(\alpha i),(\beta j)}^{\ast}(\mathbf{q})=J_{(\beta j),(\alpha i)}(\mathbf{q}).\label{eq:J_conj} +\end{equation} + +\end_inset + +In other words, +\begin_inset Formula $J(\mathbf{q})$ +\end_inset + + viewed as a matrix is Hermitian. +\end_layout + +\begin_layout Standard +We will demonstrate that +\begin_inset Formula $E$ +\end_inset + + can be equivalently calculated as, +\begin_inset Formula +\begin{equation} +\mathcal{E}=\frac{1}{2N_{\mathrm{cells}}}\int_{[0,1]^{3}}d\mathbf{q}\,S_{\alpha i}(\mathbf{q})J_{(\alpha i)(\beta j)}(\mathbf{q})S_{\beta j}^{*}(\mathbf{q}),\label{eq:Energy_fourier} +\end{equation} + +\end_inset + +There is implicit summation over indices +\begin_inset Formula $\alpha\beta ij$ +\end_inset + + and each component of +\begin_inset Formula $\mathbf{q}$ +\end_inset + + should be integrated from 0 to 1. + +\end_layout + +\begin_layout Standard +Substitution yields, +\begin_inset Formula +\begin{align} +\mathcal{E} & =\frac{1}{2N_{\mathrm{cells}}}\int_{[0,1]^{3}}d\mathbf{q}\,\sum_{\mathbf{n},\mathbf{n}',\mathbf{n}''}e^{i2\pi\mathbf{q}\cdot(\mathbf{n}-\mathbf{n}'+\mathbf{n}'')}S_{\alpha i\mathbf{n}}S_{\beta j\mathbf{n}'}\sum_{c}\delta_{c=(ij\mathbf{n}'')}J_{\alpha\beta}^{(c)}, +\end{align} + +\end_inset + +The integral over +\begin_inset Formula $\mathbf{q}$ +\end_inset + + yields a factor of +\begin_inset Formula $\delta_{\mathbf{n}-\mathbf{n}'+\mathbf{n}''}$ +\end_inset + +. + Replacing +\begin_inset Formula $\mathbf{n}'\to\mathbf{n}+\mathbf{n}''$ +\end_inset + + results in, +\begin_inset Formula +\begin{equation} +\mathcal{E}=\frac{1}{2N_{\mathrm{cells}}}\sum_{\mathbf{n},\mathbf{n}''}S_{\alpha i\mathbf{n}}S_{\beta j\mathbf{n}+\mathbf{n}''}\sum_{c}\delta_{c=(ij\mathbf{n}'')}J_{\alpha\beta}^{(c)}. +\end{equation} + +\end_inset + +Since the interactions are the same for every unit cell, the sum over cell + index +\begin_inset Formula $\mathbf{n}$ +\end_inset + + cancels the denominator +\begin_inset Formula $N_{\mathrm{cells}}$ +\end_inset + +, yielding +\begin_inset Formula +\begin{equation} +\mathcal{E}=\frac{1}{2}\sum_{\mathbf{n}''}S_{\alpha i\mathbf{0}}S_{\beta j\mathbf{n}''}\sum_{c}\delta_{c=(ij\mathbf{n}'')}J_{\alpha\beta}^{(c)}. +\end{equation} + +\end_inset + +Upon relabeling +\begin_inset Formula $\mathbf{n}''\to\mathbf{n}$ +\end_inset + +, and keeping in mind the implicit summation over +\begin_inset Formula $\alpha\beta ij$ +\end_inset + + indices, we find +\begin_inset Formula $\mathcal{E}=E$ +\end_inset + + as expected. +\end_layout + +\begin_layout Standard +A lower bound on the possible energy is given by the smallest eigenvalue + of +\begin_inset Formula $J(\mathbf{q}),$ +\end_inset + +viewed as a Hermitian matrix. + The corresponding eigenvector might be used to define the spin field at + the two momenta +\begin_inset Formula $S_{\alpha i}(\mathbf{q})=S_{\alpha i}(-\mathbf{q})^{\ast}$ +\end_inset + +, guaranteeing that the inverse Fourier transform is real. + There is no general guarantee, however, that all the spins in real-space + have the correct normalization. +\end_layout + +\begin_layout Section +Fourier-transform of dipole-dipole interactions +\end_layout + +\begin_layout Standard +The pairwise interaction between two point magnetic dipoles is given by + the tensor +\begin_inset Formula +\begin{align*} +\mathcal{D}_{\alpha\beta}(\mathbf{R}_{ij\mathbf{n}}) & =-\frac{\mu_{0}}{4\pi}\nabla_{\alpha}\nabla_{\beta}\frac{1}{|\mathbf{R}_{ij\mathbf{n}}|}\\ + & =\frac{\mu_{0}}{4\pi}\left(\delta_{\alpha\beta}\frac{1}{|\mathbf{R}_{ij\mathbf{n}}|^{3}}-3\frac{R_{ij\mathbf{n}}^{\alpha}R_{ij\mathbf{n}}^{\beta}}{|\mathbf{R}_{ij\mathbf{n}}|^{5}}\right). +\end{align*} + +\end_inset + +Note +\begin_inset Formula $\mathbf{R}_{ij\mathbf{n}}=\mathbf{r}_{ij}+\mathbf{r}_{\mathbf{n}}$ +\end_inset + +, with +\begin_inset Formula $\mathbf{r}_{\mathbf{n}}=\mathbf{A}_{\mu}\mathbf{n}_{\mu}$ +\end_inset + + and +\begin_inset Formula $\mathbf{r}_{ij}=\mathbf{r}_{j}-\mathbf{r}_{i}$ +\end_inset + +. + Here +\begin_inset Formula $\mathbf{A_{\mu}}$ +\end_inset + + denote lattice vectors for the magnetic cell. +\end_layout + +\begin_layout Standard +Fourier transform yields a dipole-dipole equivalent to +\begin_inset Formula $J_{(\alpha i)(\beta j)}(\mathbf{q})$ +\end_inset + + +\begin_inset Formula +\begin{equation} +\mathcal{D}_{(\alpha i)(\beta j)}(\mathbf{q})=\sum_{\mathbf{n}}'e^{i2\pi\mathbf{q}\cdot\mathbf{n}}\mathcal{D}_{\alpha\beta}(\mathbf{R}_{ij\mathbf{n}}).\label{eq:dipole_sum} +\end{equation} + +\end_inset + +The notation +\begin_inset Formula $\sum'$ +\end_inset + + indicates that divergent self-interactions at +\begin_inset Formula $\mathbf{R}_{ij\mathbf{n}}=0$ +\end_inset + + are excluded from the sum (i.e., the term with +\begin_inset Formula $\mathbf{n}=\mathbf{0}$ +\end_inset + + is removed when +\begin_inset Formula $i=j$ +\end_inset + +). + The remaining sum is only conditionally convergent—different results can + be obtained by summing +\begin_inset Formula $\mathbf{n}$ +\end_inset + + in different orders. +\end_layout + +\begin_layout Standard +Following the standard Ewald procedure, the infinite summation over cells + +\begin_inset Formula $\mathbf{n}$ +\end_inset + + can be rearranged into three contributions, +\end_layout + +\begin_layout Standard +\begin_inset Formula +\begin{equation} +\mathcal{D}_{(\alpha i)(\beta j)}(\mathbf{q})=\mathcal{D}_{(\alpha i)(\beta j)}^{\mathrm{LR}}(\mathbf{q},\lambda)+\mathcal{D}_{(\alpha i)(\beta j)}^{\mathrm{SR}}(\mathbf{q},\lambda)-\mathcal{D}_{(\alpha i)(\beta j)}^{\mathrm{self}}(\lambda). +\end{equation} + +\end_inset + +This decomposition introduces an arbitrary parameter +\begin_inset Formula $\lambda$ +\end_inset + +, which can be selected for numerical convenience. + Note, however, that the left-hand side is independent of +\begin_inset Formula $\lambda$ +\end_inset + +. +\end_layout + +\begin_layout Standard +The interactions at long-range are best handled in Fourier-space, +\begin_inset Formula +\begin{equation} +\mathcal{D}_{(\alpha i)(\beta j)}^{\mathrm{LR}}(\mathbf{q},\lambda)=\frac{\mu_{0}}{V}\sum_{\mathbf{m}}\mathcal{K}_{\alpha\beta}(\mathbf{k}_{\mathbf{m}+\mathbf{q}})e^{-|\mathbf{k}_{\mathbf{m}+\mathbf{q}}|^{2}/4\lambda^{2}}e^{-i\mathbf{k}_{\mathbf{m}+\mathbf{q}}\cdot\mathbf{r}_{ij}},\label{eq:D_LR} +\end{equation} + +\end_inset + +where the sum runs over +\begin_inset Formula $\mathbf{m}\in\mathbb{Z}^{3}$ +\end_inset + + and +\begin_inset Formula +\begin{equation} +\mathbf{k}_{\mathbf{q}}=q_{1}\mathbf{B}_{1}+q_{2}\mathbf{B}_{2}+q_{3}\mathbf{B}_{3}, +\end{equation} + +\end_inset + +is a dimensionful moment vector, with reciprocal vectors +\begin_inset Formula $\mathbf{B}_{\mu}$ +\end_inset + + that satisfy +\begin_inset Formula $\mathbf{A}_{\mu}\cdot\mathbf{B}_{\nu}\equiv2\pi\delta_{\mu\nu}$ +\end_inset + +. + The meaning of +\begin_inset Formula +\begin{equation} +\mathcal{K}_{\alpha\beta}(\mathbf{k})=\frac{k_{\alpha}k_{\beta}}{|\mathbf{k}|^{2}} +\end{equation} + +\end_inset + +is ambiguous the limit that +\begin_inset Formula $\mathbf{k}\to\mathbf{0}$ +\end_inset + +. + This ambiguity originates in the conditional convergence of Eq. +\begin_inset space ~ +\end_inset + + +\begin_inset CommandInset ref +LatexCommand eqref +reference "eq:dipole_sum" +plural "false" +caps "false" +noprefix "false" + +\end_inset + +. + In practice, we will take +\begin_inset Formula $\mathcal{K}_{\alpha\beta}(\mathbf{0})=0$ +\end_inset + +, which corresponds to a certain interpretation of the sample surface, and + corresponding demagnetization field. + For examples and discussion, see final appendix of +\begin_inset Flex URL +status open + +\begin_layout Plain Layout + +https://doi.org/10.1088/0953-8984/16/43/r02 +\end_layout + +\end_inset + +. +\end_layout + +\begin_layout Standard +The short-range part of the interactions can be handled directly in real-space, +\begin_inset Formula +\begin{equation} +\mathcal{D}_{(\alpha i)(\beta j)}^{\mathrm{SR}}(\mathbf{q},\lambda)=\frac{\mu_{0}}{4\pi}\sum_{\mathbf{n}}'\mathcal{U}_{\alpha\beta}(\mathbf{R}_{ij\mathbf{n}},\lambda)e^{i2\pi\mathbf{q}\cdot\mathbf{n}}, +\end{equation} + +\end_inset + +where +\begin_inset Formula +\begin{align} +\mathcal{U}_{\alpha\beta}(\mathbf{r},\lambda) & =-\nabla^{\alpha}\nabla^{\beta}\frac{\mathrm{erfc}\left(\lambda r\right)}{r}\nonumber \\ + & =\frac{\delta_{\alpha\beta}}{r^{3}}\left[\mathrm{erfc}\left(\lambda r\right)+\frac{2\lambda r}{\sqrt{\pi}}e^{-\lambda^{2}r^{2}}\right]\nonumber \\ + & \quad\frac{3\hat{r}^{\alpha}\hat{r}^{\alpha}}{r^{3}}\left[\mathrm{erfc}\left(\lambda r\right)+\left(1+\frac{2}{3}\lambda^{2}r^{2}\right)\frac{2\lambda r}{\sqrt{\pi}}e^{-\lambda^{2}r^{2}}\right]. +\end{align} + +\end_inset + +Finally, one must subtract the term, +\begin_inset Formula +\begin{align} +\mathcal{D}_{(\alpha i)(\beta j)}^{\mathrm{self}}(\lambda) & =-\frac{\mu_{0}}{2\pi}\delta_{ij}\nabla^{\alpha}\nabla^{\beta}\frac{1}{r}\mathrm{erf}\left(\lambda r\right)|_{r\to0}\nonumber \\ + & =\mu_{0}\frac{\lambda^{3}}{3\pi^{3/2}}\delta_{ij}\delta_{\alpha\beta}. +\end{align} + +\end_inset + +Physically, this can be interpreted as eliminating self-interactions that + were erroneously introduced in the Fourier space sum of Eq. +\begin_inset space ~ +\end_inset + + +\begin_inset CommandInset ref +LatexCommand eqref +reference "eq:D_LR" +plural "false" +caps "false" +noprefix "false" + +\end_inset + +. +\end_layout + +\begin_layout Section +Linear spin wave theory in Fourier space +\end_layout + +\begin_layout Subsection +Review of Holstein-Primakoff bosons +\end_layout + +\begin_layout Standard +Without loss of generality, assume we have rotated into a +\emph on +local frame +\emph default + such that the magnetic ground state is uniformly oriented in the +\begin_inset Formula $\hat{z}$ +\end_inset + + direction. + Introduce on each sublattice +\begin_inset Formula $j$ +\end_inset + + and magnetic cell +\begin_inset Formula $\mathbf{n}$ +\end_inset + + a Holstein-Primakoff boson +\begin_inset Formula $b_{j,\mathbf{n}}$ +\end_inset + +. + Spin operators satisfying SU(2) commutation relations can be constructed + as, +\begin_inset Formula +\begin{align} +\hat{S}_{j,\mathbf{n}}^{+} & \equiv\sqrt{2s}\sqrt{1-\frac{b_{j,\mathbf{n}}^{\dagger}b_{j,\mathbf{n}}}{2s}}b_{j,\mathbf{n}}\\ +\hat{S}_{j,\mathbf{n}}^{z} & \equiv s-b_{j,\mathbf{n}}^{\dagger}b_{j,\mathbf{n}}. +\end{align} + +\end_inset + +Expanding in small +\begin_inset Formula $1/s$ +\end_inset + + produces an effective expansion in powers of bosons, +\begin_inset Formula +\begin{equation} +\hat{S}_{j,\mathbf{n}}^{+}\approx\sqrt{2s}\,b_{j,\mathbf{n}},\quad\hat{S}_{j,\mathbf{n}}^{-}\approx\sqrt{2s}\,b_{j,\mathbf{n}}^{\dagger}. +\end{equation} + +\end_inset + +Using +\begin_inset Formula $\hat{S}^{\pm}=\hat{S}^{x}\pm i\hat{S}^{y}$ +\end_inset + +, the result to quadratic order is, +\begin_inset Formula +\begin{align} +\hat{S}_{j,\mathbf{n}}^{x} & =\frac{\sqrt{2s}}{2}(b_{j,\mathbf{n}}+b_{j,\mathbf{n}}^{\dagger})\\ +\hat{S}_{j,\mathbf{n}}^{y} & =\frac{\sqrt{2s}}{2i}(b_{j,\mathbf{n}}-b_{j,\mathbf{n}}^{\dagger}) +\end{align} + +\end_inset + +Now Fourier transform, +\begin_inset Formula +\begin{equation} +b_{j,\mathbf{q}}\equiv\sum_{\mathbf{n}}e^{2\pi i\mathbf{q}\cdot\mathbf{n}}b_{j,\mathbf{n}}, +\end{equation} + +\end_inset + +with, +\begin_inset Formula +\begin{equation} +b_{j,\mathbf{q}}^{\dagger}\equiv(b_{j,\mathbf{q}})^{\dagger}=\sum_{\mathbf{n}}e^{-2\pi i\mathbf{q}\cdot\mathbf{n}}b_{j,\mathbf{n}}^{\dagger}. +\end{equation} + +\end_inset + +The Fourier transform of spin components, +\begin_inset Formula +\begin{equation} +\hat{S}_{j,\mathbf{q}}^{\alpha}\equiv\sum_{\mathbf{n}}e^{2\pi i\mathbf{q}\cdot\mathbf{n}}\hat{S}_{j,\mathbf{n}}, +\end{equation} + +\end_inset + +is, up to quadratic order, +\begin_inset Formula +\begin{align} +\hat{S}_{j,\mathbf{q}}^{x} & \approx\frac{\sqrt{2s}}{2}(b_{j,\mathbf{q}}+b_{j,\mathbf{-q}}^{\dagger})\label{eq:Sx_q}\\ +\hat{S}_{j,\mathbf{q}}^{y} & \approx\frac{\sqrt{2s}}{2i}(b_{j,\mathbf{q}}-b_{j,\mathbf{-q}}^{\dagger})\\ +\hat{S}_{j,\mathbf{q}}^{z} & \approx s\delta(\mathbf{q})-\int_{[0,1]^{3}}d\mathbf{q}'\,b_{j,\mathbf{q}'}^{\dagger}b_{j,\mathbf{q}+\mathbf{q}'}.\label{eq:Sz_q} +\end{align} + +\end_inset + +The last equation derives from the +\begin_inset Quotes eld +\end_inset + +cross-correlation theorem +\begin_inset Quotes erd +\end_inset + + (direct analog of convolution theorem). +\end_layout + +\begin_layout Subsection +LSWT Hamiltonian +\end_layout + +\begin_layout Standard +Adapting Eq. +\begin_inset space ~ +\end_inset + + +\begin_inset CommandInset ref +LatexCommand eqref +reference "eq:Energy_fourier" +plural "false" +caps "false" +noprefix "false" + +\end_inset + + to the quantum setting, the quadratic Hamiltonian becomes +\begin_inset Formula +\begin{equation} +\mathcal{H}=\frac{1}{2N_{\mathrm{cells}}}\int_{[0,1]^{3}}d\mathbf{q}\,\hat{\mathbf{S}}_{i,\mathbf{q}}^{\dagger}J_{ij,\mathbf{q}}\hat{\mathbf{S}}_{j,\mathbf{q}}, +\end{equation} + +\end_inset + +with the same +\begin_inset Formula $3\times3$ +\end_inset + + coupling matrices +\begin_inset Formula $J_{ij,\mathbf{q}}$ +\end_inset + + as appear in the classical case, Eq. +\begin_inset space ~ +\end_inset + + +\begin_inset CommandInset ref +LatexCommand eqref +reference "eq:J_q" +plural "false" +caps "false" +noprefix "false" + +\end_inset + +. +\end_layout + +\begin_layout Standard +To perform an expansion in powers of +\begin_inset Formula $1/s$ +\end_inset + +, it is convenient to locally rotate the coordinate system on each sublattice + +\begin_inset Formula $j$ +\end_inset + +, +\begin_inset Formula +\begin{align} +\tilde{\mathbf{S}}_{j,\mathbf{q}} & \equiv R_{j}^{T}\hat{\mathbf{S}}_{j,\mathbf{q}}\\ +\tilde{J}_{ij,\mathbf{q}} & \equiv R_{i}^{T}J_{ij,\mathbf{q}}R_{j}. +\end{align} + +\end_inset + +The rotations +\begin_inset Formula $R_{j}$ +\end_inset + + should be constructed such that, in the new local frame, ground state is + ferromagnetically aligned in the +\begin_inset Formula $\hat{z}$ +\end_inset + +-direction: +\begin_inset Formula $\langle\tilde{\mathbf{S}}_{j,\mathbf{q}}\rangle\propto\hat{\mathbf{z}}$ +\end_inset + +. + The quantum Hamiltonian in the local frame has the same form, +\begin_inset Formula +\begin{equation} +\mathcal{H}=\frac{1}{2N_{\mathrm{cells}}}\int_{[0,1]^{3}}d\mathbf{q}\,\tilde{J}_{ij,\mathbf{q}}^{\alpha\beta}\tilde{S}_{i,\mathbf{q}}^{\alpha\dagger}\tilde{S}_{j,\mathbf{q}}^{\beta}.\label{eq:H_rotated} +\end{equation} + +\end_inset + +Note that the rotated dipoles +\begin_inset Formula $\tilde{\mathbf{S}}_{j,\mathbf{q}}$ +\end_inset + + are still quantum spin operators (i.e., satisfy the SU(2) commutation relations + when expressed in real-space). +\end_layout + +\begin_layout Standard +Informally, an approximation can be obtained by assuming that the deviations + from the ground state are small, +\begin_inset Formula +\begin{equation} +\tilde{\mathbf{S}}_{j,\mathbf{q}}\approx\left\langle \tilde{\mathbf{S}}_{j,\mathbf{q}}\right\rangle . +\end{equation} + +\end_inset + +Formally, this is achieved by substituting Eqs. +\begin_inset space ~ +\end_inset + + +\begin_inset CommandInset ref +LatexCommand eqref +reference "eq:Sx_q" +plural "false" +caps "false" +noprefix "false" + +\end_inset + +– +\begin_inset CommandInset ref +LatexCommand eqref +reference "eq:Sz_q" +plural "false" +caps "false" +noprefix "false" + +\end_inset + + for the locally rotated dipoles +\begin_inset Formula $\tilde{\mathbf{S}}_{j,\mathbf{q}}$ +\end_inset + + and expanding in powers of bosons (equivalently, in powers of +\begin_inset Formula $1/s$ +\end_inset + +). + Up to quadratic order in bosons, the result is, +\begin_inset Formula +\begin{equation} +\mathcal{H}\approx\mathcal{H}_{0}+\cancel{\mathcal{H}_{1}}+\mathcal{H}_{2}. +\end{equation} + +\end_inset + +At zeroth order, +\begin_inset Formula $\mathcal{H}_{0}=\langle\mathcal{H}\rangle_{\mathrm{classical}}$ +\end_inset + + represents the ground state energy of the classical Hamiltonian. + At first order, +\begin_inset Formula $\mathcal{H}_{1}$ +\end_inset + + vanishes, as it is proportional to the gradient of the classical Hamiltonian + in the ground state. +\end_layout + +\begin_layout Standard +The terms contributing to +\begin_inset Formula $\mathcal{H}_{2}$ +\end_inset + + are obtained by expanding +\begin_inset Formula $\tilde{S}_{i,\mathbf{q}}^{\alpha\dagger}\tilde{S}_{i,\mathbf{q}}^{\beta\dagger}$ +\end_inset + + for all spin-components +\begin_inset Formula $\alpha$ +\end_inset + + and +\begin_inset Formula $\beta$ +\end_inset + + and retaining only terms at quadratic order in bosons. + The steps are lengthy but mechanical. +\end_layout + +\begin_layout Standard +Within the oscillating plane, +\begin_inset Formula +\begin{align} +\tilde{S}_{i,\mathbf{q}}^{x\dagger}\tilde{S}_{j,\mathbf{q}}^{x} & \to\frac{s}{2}(b_{i,\mathbf{q}}^{\dagger}b_{j,\mathbf{q}}+b_{i,-\mathbf{q}}b_{j,-\mathbf{q}}^{\dagger}+b_{i,\mathbf{q}}^{\dagger}b_{j,-\mathbf{q}}^{\dagger}+b_{i,-\mathbf{q}}b_{j,\mathbf{q}})\\ +\tilde{S}_{i,\mathbf{q}}^{y\dagger}\tilde{S}_{j,\mathbf{q}}^{y} & \to\frac{s}{2}(b_{i,\mathbf{q}}^{\dagger}b_{j,\mathbf{q}}+b_{i,-\mathbf{q}}b_{j,-\mathbf{q}}^{\dagger}-b_{i,\mathbf{q}}^{\dagger}b_{j,-\mathbf{q}}^{\dagger}-b_{i,-\mathbf{q}}b_{j,\mathbf{q}})\\ +\tilde{S}_{i,\mathbf{q}}^{x\dagger}\tilde{S}_{j,\mathbf{q}}^{y} & \to\frac{s}{2}(-ib_{i,\mathbf{q}}^{\dagger}b_{j,\mathbf{q}}+ib_{i,-\mathbf{q}}b_{j,-\mathbf{q}}^{\dagger}+ib_{i,\mathbf{q}}^{\dagger}b_{j,-\mathbf{q}}^{\dagger}-ib_{i,-\mathbf{q}}b_{j,\mathbf{q}})\\ +\tilde{S}_{i,\mathbf{q}}^{y\dagger}\tilde{S}_{j,\mathbf{q}}^{x} & \to\frac{s}{2}(+ib_{i,\mathbf{q}}^{\dagger}b_{j,\mathbf{q}}-ib_{i,\mathbf{-q}}b_{j,\mathbf{-q}}^{\dagger}+ib_{i,\mathbf{q}}^{\dagger}b_{j,\mathbf{-q}}^{\dagger}-ib_{i,\mathbf{-q}}b_{j,\mathbf{q}}). +\end{align} + +\end_inset + +At quadratic order, terms linear in +\begin_inset Formula $\tilde{S}_{j,\mathbf{q}}^{z}$ +\end_inset + + vanish, +\begin_inset Formula +\begin{equation} +\tilde{S}_{i,\mathbf{q}}^{x\dagger}\tilde{S}_{j,\mathbf{q}}^{z}\to0,\quad\tilde{S}_{i,\mathbf{q}}^{y\dagger}\tilde{S}_{j,\mathbf{q}}^{z}\to0. +\end{equation} + +\end_inset + +Finally, +\begin_inset Formula +\begin{equation} +\tilde{S}_{i,\mathbf{q}}^{z\dagger}\tilde{S}_{j,\mathbf{q}}^{z}\to-s\,\delta(\mathbf{q})\int_{[0,1]^{3}}d\mathbf{q}'\,(b_{i,\mathbf{q}'}^{\dagger}b_{i,\mathbf{q}'}+b_{j,\mathbf{q}'}^{\dagger}b_{j,\mathbf{q}'}).\label{eq:SzSz} +\end{equation} + +\end_inset + +Substituting these expansions into Eq. +\begin_inset space ~ +\end_inset + + +\begin_inset CommandInset ref +LatexCommand eqref +reference "eq:H_rotated" +plural "false" +caps "false" +noprefix "false" + +\end_inset + + yields, +\end_layout + +\begin_layout Standard +\begin_inset Formula +\begin{align} +\mathcal{H}_{2}=\frac{s}{2N_{\mathrm{cells}}}\int_{[0,1]^{3}}d\mathbf{q}\, & \Bigl[\frac{b_{i,\mathbf{q}}^{\dagger}b_{j,\mathbf{q}}}{2}(\tilde{J}_{ij,\mathbf{q}}^{xx}+\tilde{J}_{ij,\mathbf{q}}^{yy}-i\tilde{J}_{ij,\mathbf{q}}^{xy}+i\tilde{J}_{ij,\mathbf{q}}^{yx})+\nonumber \\ + & \frac{b_{i,\mathbf{-q}}b_{j,\mathbf{-q}}^{\dagger}}{2}(\tilde{J}_{ij,\mathbf{q}}^{xx}+\tilde{J}_{ij,\mathbf{q}}^{yy}+i\tilde{J}_{ij,\mathbf{q}}^{xy}-i\tilde{J}_{ij,\mathbf{q}}^{yx})+\nonumber \\ + & \frac{b_{i,\mathbf{q}}^{\dagger}b_{j,-\mathbf{q}}^{\dagger}}{2}(\tilde{J}_{ij,\mathbf{q}}^{xx}-\tilde{J}_{ij,\mathbf{q}}^{yy}+i\tilde{J}_{ij,\mathbf{q}}^{xy}+i\tilde{J}_{ij,\mathbf{q}}^{yx})+\nonumber \\ + & \frac{b_{i,-\mathbf{q}}b_{j,\mathbf{q}}}{2}(\tilde{J}_{ij,\mathbf{q}}^{xx}-\tilde{J}_{ij,\mathbf{q}}^{yy}-i\tilde{J}_{ij,\mathbf{q}}^{xy}-i\tilde{J}_{ij,\mathbf{q}}^{yx})+\nonumber \\ + & -s\,\left(b_{i,\mathbf{q}}^{\dagger}b_{i,\mathbf{q}}+b_{j,\mathbf{q}}^{\dagger}b_{j,\mathbf{q}}\right)\tilde{J}_{ij,\mathbf{0}}^{zz}\Bigr]. +\end{align} + +\end_inset + +Observe that the final term in the integrand involves coefficients +\begin_inset Formula $\tilde{J}_{ij,\mathbf{0}}^{zz}$ +\end_inset + + evaluated at the +\emph on +zero +\emph default + wavevector, which originates from the Dirac- +\begin_inset Formula $\delta$ +\end_inset + + appearing in Eq. +\begin_inset space ~ +\end_inset + + +\begin_inset CommandInset ref +LatexCommand eqref +reference "eq:SzSz" +plural "false" +caps "false" +noprefix "false" + +\end_inset + +. +\end_layout + +\begin_layout Standard +It is conventional to express this bilinear form using bosons concatenated + into a +\begin_inset Quotes eld +\end_inset + +Nambu spinor, +\begin_inset Quotes erd +\end_inset + + +\begin_inset Formula +\begin{equation} +\Psi_{\mathbf{q}}^{\dagger}=[b_{1,\mathbf{q}}^{\dagger},\dots,b_{N,\mathbf{q}}^{\dagger},b_{1,-\mathbf{q}},\dots,b_{N,-\mathbf{q}}], +\end{equation} + +\end_inset + +Then, +\begin_inset Formula +\begin{equation} +\mathcal{H}_{2}=\frac{1}{2N_{\mathrm{cells}}}\int_{[0,1]^{3}}d\mathbf{q}\Psi_{\mathbf{q}}^{\dagger}D_{\mathbf{q}}\Psi_{\mathbf{q}}. +\end{equation} + +\end_inset + +where the +\begin_inset Formula $2N\times2N$ +\end_inset + + dynamical matrix takes the form, +\begin_inset Formula +\begin{equation} +D_{\mathbf{q}}=\left(\begin{array}{cc} +A_{\mathbf{q}} & B_{\mathbf{q}}\\ +B_{\mathbf{q}}^{\dagger} & A_{-\mathbf{q}} +\end{array}\right). +\end{equation} + +\end_inset + + +\end_layout + +\begin_layout Standard +\begin_inset Note Comment +status open + +\begin_layout Plain Layout +Because +\begin_inset Formula $J_{q}$ +\end_inset + + was obtained by Fourier transformation of +\emph on +real +\emph default +exchange matrices, it follows that +\begin_inset Formula $J_{q}^{*}=J_{-q}$ +\end_inset + +. + Also, per Eq. +\begin_inset space ~ +\end_inset + + +\begin_inset CommandInset ref +LatexCommand eqref +reference "eq:J_conj" +plural "false" +caps "false" +noprefix "false" + +\end_inset + +, it is valid to symmetrize all interaction matrices over the two bond orientiat +ions +\begin_inset Formula $(ij)$ +\end_inset + + and +\begin_inset Formula $(ji)$ +\end_inset + +. + Without loss of generality, one can assume, +\begin_inset Formula +\begin{equation} +(J_{q,ij}^{\alpha\beta})^{*}=J_{-q,ij}^{\alpha\beta}=J_{q,ji}^{\beta\alpha}. +\end{equation} + +\end_inset + +Spin operators are Hermitian, and this is also satisfied by the expansion, +\begin_inset Formula +\begin{equation} +\hat{S}_{q,i}^{\alpha\dagger}=\hat{S}_{-q,i}^{\alpha}. +\end{equation} + +\end_inset + + +\end_layout + +\end_inset + + +\end_layout + +\end_body +\end_document diff --git a/docs/lyx_notes/phase_averaging.lyx b/docs/lyx_notes/phase_averaging.lyx deleted file mode 100644 index 3b81af64b..000000000 --- a/docs/lyx_notes/phase_averaging.lyx +++ /dev/null @@ -1,459 +0,0 @@ -#LyX 2.3 created this file. For more info see http://www.lyx.org/ -\lyxformat 544 -\begin_document -\begin_header -\save_transient_properties true -\origin unavailable -\textclass article -\use_default_options true -\maintain_unincluded_children false -\language english -\language_package default -\inputencoding auto -\fontencoding global -\font_roman "default" "default" -\font_sans "default" "default" -\font_typewriter "default" "default" -\font_math "auto" "auto" -\font_default_family default -\use_non_tex_fonts false -\font_sc false -\font_osf false -\font_sf_scale 100 100 -\font_tt_scale 100 100 -\use_microtype false -\use_dash_ligatures true -\graphics default -\default_output_format default -\output_sync 0 -\bibtex_command default -\index_command default -\paperfontsize default -\use_hyperref false -\papersize default -\use_geometry false -\use_package amsmath 1 -\use_package amssymb 1 -\use_package cancel 1 -\use_package esint 1 -\use_package mathdots 1 -\use_package mathtools 1 -\use_package mhchem 1 -\use_package stackrel 1 -\use_package stmaryrd 1 -\use_package undertilde 1 -\cite_engine basic -\cite_engine_type default -\use_bibtopic false -\use_indices false -\paperorientation portrait -\suppress_date false -\justification true -\use_refstyle 1 -\use_minted 0 -\index Index -\shortcut idx -\color #008000 -\end_index -\secnumdepth 3 -\tocdepth 3 -\paragraph_separation indent -\paragraph_indentation default -\is_math_indent 0 -\math_numbering_side default -\quotes_style english -\dynamic_quotes 0 -\papercolumns 1 -\papersides 1 -\paperpagestyle default -\tracking_changes false -\output_changes false -\html_math_output 0 -\html_css_as_file 0 -\html_be_strict false -\end_header - -\begin_body - -\begin_layout Title -Phase averaging -\end_layout - -\begin_layout Standard -The Bravais lattice for the lattice vectors -\begin_inset Formula $\{\mathbf{a}_{1},\mathbf{a}_{2},\mathbf{a}_{3}\}$ -\end_inset - - is the set of positions -\begin_inset Formula -\begin{equation} -\mathbf{r}_{n_{1},n_{2},n_{3}}=n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}+n_{3}\mathbf{a}_{3}, -\end{equation} - -\end_inset - -where -\begin_inset Formula $\{n_{1},n_{2},n_{3}\}$ -\end_inset - - are integer. - Decorated crystals require an additional sublattice index -\begin_inset Formula $\ell=1,\dots N_{\mathrm{sublat.}}$ -\end_inset - -, and the associated atom positions are -\begin_inset Formula -\begin{equation} -\mathbf{r}_{\mathbf{n},\ell}=\mathbf{r}_{n_{1},n_{2},n_{3},\ell}=n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}+n_{3}\mathbf{a}_{3}+\mathbf{r}_{\ell}^{0}.\label{eq:r_position} -\end{equation} - -\end_inset - -The relative shift between sublattices -\begin_inset Formula $\mathbf{r}_{\ell}^{0}-\mathbf{r}_{\ell'}^{0}$ -\end_inset - - will play an important role in calculating structure factor data, as we - will demonstrate below. -\end_layout - -\begin_layout Standard -Consider a finite magnetic unit cell, which is defined by a set of spins - -\begin_inset Formula $\mathbf{s}_{\mathbf{n},\ell}$ -\end_inset - - satisfying some periodicity. - That is, we will assume existence of some integer lengths -\begin_inset Formula $\{L_{1},L_{2},L_{3}\}$ -\end_inset - - such that -\begin_inset Formula $\mathbf{s}_{\mathbf{n},\ell}$ -\end_inset - - is invariant under the substitution -\begin_inset Formula $n_{i}\rightarrow n_{i}+L_{i}$ -\end_inset - - for each -\begin_inset Formula $i\in\{1,2,3\}$ -\end_inset - -. -\end_layout - -\begin_layout Standard -In the continuum, each discrete spin -\begin_inset Formula $\mathbf{s}_{\mathbf{n},\ell}$ -\end_inset - - is promoted to a Dirac- -\begin_inset Formula $\delta$ -\end_inset - - peak. - The continuum -\begin_inset Quotes eld -\end_inset - -field -\begin_inset Quotes erd -\end_inset - - of spins is defined to be -\begin_inset Formula -\begin{equation} -\mathbf{s}(\mathbf{r})=\sum_{\mathbf{n}}\sum_{\ell=1}^{N_{\mathrm{sublat.}}}\mathbf{s}_{\mathbf{n},\ell}\delta(\mathbf{r}-\mathbf{r}_{\mathbf{n},\ell}), -\end{equation} - -\end_inset - -where the -\begin_inset Formula $n_{i}$ -\end_inset - - sums run over all integers. - Periodicity of -\begin_inset Formula $\mathbf{s}_{\mathbf{n},\ell}$ -\end_inset - - implies periodicity of -\begin_inset Formula $\mathbf{s}(\mathbf{r})$ -\end_inset - -. - Namely, -\begin_inset Formula $\mathbf{s}(\mathbf{r}+L_{i}\mathbf{a}_{i})=\mathbf{s}(\mathbf{r})$ -\end_inset - - for -\begin_inset Formula $i\in\{1,2,3\}$ -\end_inset - -. -\end_layout - -\begin_layout Standard -The Fourier transform of the continuum spin field is, -\begin_inset Formula -\begin{equation} -\hat{\mathbf{s}}(\mathbf{k})=\int_{\mathbb{R}^{3}}e^{-i\mathbf{k}\cdot\mathbf{r}}\mathbf{s}(\mathbf{r})\,d\mathbf{r}=\sum_{\mathbf{n}}\sum_{\ell=1}^{N_{\mathrm{sublat.}}}e^{-i\mathbf{k}\cdot\mathbf{r}_{\mathbf{n},\ell}}\mathbf{s}_{\mathbf{n},\ell}. -\end{equation} - -\end_inset - -Substituting Eq. -\begin_inset space ~ -\end_inset - - -\begin_inset CommandInset ref -LatexCommand eqref -reference "eq:r_position" -plural "false" -caps "false" -noprefix "false" - -\end_inset - - and rearranging terms, -\begin_inset Formula -\begin{equation} -\hat{\mathbf{s}}(\mathbf{k})=\sum_{\ell=1}^{N_{\mathrm{sublat.}}}e^{-i\mathbf{k}\cdot\mathbf{r}_{\ell}^{0}}\left[\sum_{\mathbf{n}}e^{-i\mathbf{k}\cdot(n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}+n_{3}\mathbf{a}_{3})}\mathbf{s}_{\mathbf{n},\ell}\right].\label{eq:shat_k_1} -\end{equation} - -\end_inset - -Due to the assumed periodicity of -\begin_inset Formula $\mathbf{s}_{\mathbf{n},\ell}$ -\end_inset - -, the sum in brackets is nonzero only for a regular grid of Fourier wave - vectors -\begin_inset Formula $\mathbf{k}$ -\end_inset - -, -\begin_inset Formula -\begin{equation} -\mathbf{k}=\frac{m_{1}}{L_{1}}\mathbf{b}_{1}+\frac{m_{2}}{L_{2}}\mathbf{b}_{2}+\frac{m_{3}}{L_{3}}\mathbf{b}_{3},\label{eq:k_to_m} -\end{equation} - -\end_inset - -where -\begin_inset Formula $m_{i}\in\mathbb{Z}$ -\end_inset - -, and the reciprocal lattice vectors -\begin_inset Formula $\mathbf{b}_{i}$ -\end_inset - - are defined to satisfy -\end_layout - -\begin_layout Standard -\begin_inset Formula -\begin{equation} -\mathbf{a}_{i}\cdot\mathbf{b}_{j}=2\pi\delta_{ij}. -\end{equation} - -\end_inset - -Equation -\begin_inset space ~ -\end_inset - - -\begin_inset CommandInset ref -LatexCommand eqref -reference "eq:k_to_m" -plural "false" -caps "false" -noprefix "false" - -\end_inset - - defines the -\emph on -reciprocal -\emph default -( -\begin_inset Formula $k$ -\end_inset - --space) lattice. -\end_layout - -\begin_layout Standard -The -\begin_inset Formula $\mathbf{n}$ -\end_inset - - summation may be written -\begin_inset Formula -\begin{equation} -\sum_{\mathbf{n}}e^{-i\mathbf{k}\cdot(n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}+n_{3}\mathbf{a}_{3})}\mathbf{s}_{\mathbf{n},\ell}=\mathcal{V}\hat{\mathbf{s}}_{\mathbf{m},\ell}, -\end{equation} - -\end_inset - -where -\begin_inset Formula $\hat{\mathbf{s}}_{\mathbf{m},\ell}$ -\end_inset - - is the usual 3D discrete Fourier transform of the magnetic unit cell data - -\begin_inset Formula $\mathbf{s}_{\mathbf{n},\ell}$ -\end_inset - -, and -\begin_inset Formula $\mathcal{V}$ -\end_inset - - is a measure of total system volume (specifically, the count of magnetic - unit cells). - The final result becomes -\end_layout - -\begin_layout Standard -\begin_inset Box Boxed -position "t" -hor_pos "c" -has_inner_box 1 -inner_pos "t" -use_parbox 0 -use_makebox 0 -width "100col%" -special "none" -height "1in" -height_special "totalheight" -thickness "0.4pt" -separation "3pt" -shadowsize "4pt" -framecolor "black" -backgroundcolor "none" -status open - -\begin_layout Plain Layout -\begin_inset Formula -\begin{equation} -\hat{\mathbf{s}}(\mathbf{k})=\mathcal{V}\sum_{\ell=1}^{N_{\mathrm{sublat.}}}e^{-i\mathbf{k}\cdot\mathbf{r}_{\ell}^{0}}\hat{\mathbf{s}}_{\mathbf{m},\ell}, -\end{equation} - -\end_inset - - -\end_layout - -\end_inset - - -\end_layout - -\begin_layout Standard -Although -\begin_inset Formula $\hat{\mathbf{s}}_{\mathbf{m},\ell}$ -\end_inset - - is periodic under -\begin_inset Formula $m_{i}\rightarrow m_{i}+L_{i}$ -\end_inset - -, the continuum-space Fourier transform -\begin_inset Formula $\hat{\mathbf{s}}(\mathbf{k})$ -\end_inset - - does -\emph on -not -\emph default -have this simple periodicity. -\end_layout - -\begin_layout Standard - -\series bold -Example: -\series default - Consider, for example, a shift of -\begin_inset Formula $\mathbf{k}$ -\end_inset - - along one of the reciprocal lattice vectors, -\begin_inset Formula -\[ -\mathbf{k}\rightarrow\mathbf{k}+\mathbf{b}_{1}. -\] - -\end_inset - -A Bravais lattice ( -\begin_inset Formula $N_{\mathrm{sublat.}}=1)$ -\end_inset - - does have periodicity, -\begin_inset Formula $\hat{\mathbf{s}}(\mathbf{k})=\hat{\mathbf{s}}(\mathbf{k}+\mathbf{b}_{1})$ -\end_inset - -. - More generally, for a decorated lattice, -\begin_inset Formula -\begin{align*} -\hat{\mathbf{s}}\left(\mathbf{k}+\mathbf{b}_{1}\right) & =\mathcal{V}\sum_{\ell=1}^{N_{\mathrm{sublat.}}}e^{-i(\mathbf{k}+\mathbf{b}_{1})\cdot\mathbf{r}_{\ell}^{0}}\hat{\mathbf{s}}_{\{m_{1}+L,m_{2},m_{3}\},\ell}\\ - & =\mathcal{V}\sum_{\ell=1}^{N_{\mathrm{sublat.}}}e^{-i(\mathbf{k}+\mathbf{b}_{1})\cdot\mathbf{r}_{\ell}^{0}}\hat{\mathbf{s}}_{\mathbf{m},\ell}\\ - & \neq\mathcal{V}\sum_{\ell=1}^{N_{\mathrm{sublat.}}}e^{-i\mathbf{k}\cdot\mathbf{r}_{\ell}^{0}}\hat{\mathbf{s}}_{\mathbf{m},\ell}\\ - & =\hat{\mathbf{s}}\left(\mathbf{k}\right) -\end{align*} - -\end_inset - -Observe that -\begin_inset Formula $\hat{\mathbf{s}}\left(\mathbf{k}+\mathbf{b}_{1}\right)$ -\end_inset - - cannot generally be expressed as a function of -\begin_inset Formula $\hat{\mathbf{s}}\left(\mathbf{k}\right)$ -\end_inset - - for fixed -\begin_inset Formula $\mathbf{k}$ -\end_inset - --vector. - That is, calculation of the value -\begin_inset Formula $\hat{\mathbf{s}}\left(\mathbf{k}+\mathbf{b}_{1}\right)$ -\end_inset - - requires new phase-averaging over the discrete data -\begin_inset Formula $\hat{\mathbf{s}}_{\mathbf{m},\ell}$ -\end_inset - -, where -\begin_inset Formula $\mathbf{m}$ -\end_inset - - and -\begin_inset Formula $\mathbf{k}$ -\end_inset - - are related via Eq. -\begin_inset space ~ -\end_inset - - -\begin_inset CommandInset ref -LatexCommand eqref -reference "eq:k_to_m" -plural "false" -caps "false" -noprefix "false" - -\end_inset - -. -\end_layout - -\end_body -\end_document diff --git a/examples/07_Dipole_Dipole.jl b/examples/07_Dipole_Dipole.jl new file mode 100644 index 000000000..07dc231fb --- /dev/null +++ b/examples/07_Dipole_Dipole.jl @@ -0,0 +1,63 @@ +# # 7. Long-range dipole interactions +# +# This example demonstrates Sunny's ability to incorporate long-range dipole +# interactions using Ewald summation. The calculation reproduces previous +# results in [Del Maestro and Gingras, J. Phys.: Cond. Matter, **16**, 3339 +# (2004)](https://arxiv.org/abs/cond-mat/0403494). + +using Sunny, GLMakie + +# Create a Pyrochlore lattice from spacegroup 227 in the standard setting. +latvecs = lattice_vectors(10.19, 10.19, 10.19, 90, 90, 90) +cryst = Crystal(latvecs, [[0,0,0]], 227, setting="2") + +sys = System(cryst, (1, 1, 1), [SpinInfo(1, S=7/2, g=2)], :dipole, seed=2) + +J1 = 0.304 * meV_per_K +set_exchange!(sys, J1, Bond(1,2,[0,0,0])) + +shape = [1/2 1/2 0; 0 1/2 1/2; 1/2 0 1/2] +sys_prim = reshape_supercell(sys, shape) + +set_dipole!(sys_prim, [+1, -1, 0], position_to_site(sys_prim, [0,0,0])) +set_dipole!(sys_prim, [-1, +1, 0], position_to_site(sys_prim, [1/4,1/4,0])) +set_dipole!(sys_prim, [+1, +1, 0], position_to_site(sys_prim, [1/4,0,1/4])) +set_dipole!(sys_prim, [-1, -1, 0], position_to_site(sys_prim, [0,1/4,1/4])) + +plot_spins(sys_prim) + +q_points = [[0,0,0], [0,1,0], [1,1/2,0], [1/2,1/2,1/2], [3/4,3/4,0], [0,0,0]] +q_labels = ["Γ","X","W","L","K","Γ"] +density = 150 +path, xticks = reciprocal_space_path(cryst, q_points, density); +xticks = (xticks[1], q_labels) + +swt = SpinWaveTheory(sys_prim, energy_ϵ=1e-6) +disp1 = dispersion(swt, path) + +enable_dipole_dipole!(sys_prim) +swt = SpinWaveTheory(sys_prim, energy_ϵ=1e-6) +disp2 = dispersion(swt, path) + +# Plot dispersions with and without long-range dipole interactions. To reproduce +# Fig. 2 of [Del Maestro and Gingras](https://arxiv.org/abs/cond-mat/0403494), +# an empirical rescaling of energy is necessary. +fudge_factor = 1/2 + +fig = Figure(size=(900,300)) + +ax = Axis(fig[1,1]; xlabel="", ylabel="Energy (K)", xticks, xticklabelrotation=0) +ylims!(ax, 0, 2) +xlims!(ax, 1, size(disp1, 1)) +for i in axes(disp1, 2) + lines!(ax, 1:length(disp1[:,i]), fudge_factor*disp1[:,i]/meV_per_K) +end + +ax = Axis(fig[1,2]; xlabel="", ylabel="Energy (K)", xticks, xticklabelrotation=0) +ylims!(ax, 0.0, 3) +xlims!(ax, 1, size(disp2, 1)) +for i in axes(disp2, 2) + lines!(ax, 1:length(disp2[:,i]), fudge_factor*disp2[:,i]/meV_per_K) +end + +fig From 54d36bac8cef60e59e927b9eaa5a5c608ce1dc7c Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 15/29] Update readme and versions --- README.md | 4 +--- docs/src/versions.md | 3 +++ 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index f59e6d1e3..b077742aa 100644 --- a/README.md +++ b/README.md @@ -50,11 +50,9 @@ Sunny is inspired by SpinW, especially regarding symmetry analysis, model specif | [Interaction renormalization for dipoles](https://arxiv.org/abs/2304.03874) | ❌ | ❌ | ✅ | | [Multi-flavor spin wave theory](https://arxiv.org/abs/1307.7731) | ✅ | ❌ | ✅ | | [Classical SU(_N_) spin dynamics](https://arxiv.org/abs/2209.01265) | ❌ | ❌ | ✅ | -| Ewald summation for dipole-dipole | ❌ | ❌ | 🟨⁽¹⁾ | +| Ewald summation for dipole-dipole | ❌ | ❌ | ✅ | | Programming language | C++ | Matlab | [Julia](https://julialang.org/) | -_Fine print: (1) Dipole-dipole interactions currently supported in classical dyamics but not LSWT._ - The classical SU(_N_) spin dynamics in Sunny generalizes the Landau-Lifshitz equation for $S > 1/2$ quantum spins. Linearizing and quantizing SU(_N_) dynamics yields a generalization of spin wave theory involving multi-flavor bosons. Codes like [Spirit](https://github.com/spirit-code/spirit) and [Vampire](https://vampire.york.ac.uk/) focus less on capturing quantum effects, but may be good options for classical dynamics of pure dipoles. diff --git a/docs/src/versions.md b/docs/src/versions.md index 8e82d2bb1..5f6dd87ac 100644 --- a/docs/src/versions.md +++ b/docs/src/versions.md @@ -5,6 +5,9 @@ * [`view_crystal`](@ref) called on a [`System`](@ref) now optionally shows spin or magnetic dipoles. +* Interactions for [`enable_dipole_dipole!`](@ref) are now supported in linear + spin wave theory, with proper Ewald summation. + ## v0.5.9 (Mar 25, 2024) From 2fb833c375ecac066fd2bc2b59129395d93007e1 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 16/29] Improve comment in example. --- examples/07_Dipole_Dipole.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/07_Dipole_Dipole.jl b/examples/07_Dipole_Dipole.jl index 07dc231fb..3eb62f7d0 100644 --- a/examples/07_Dipole_Dipole.jl +++ b/examples/07_Dipole_Dipole.jl @@ -37,12 +37,13 @@ disp1 = dispersion(swt, path) enable_dipole_dipole!(sys_prim) swt = SpinWaveTheory(sys_prim, energy_ϵ=1e-6) -disp2 = dispersion(swt, path) +disp2 = dispersion(swt, path); # Plot dispersions with and without long-range dipole interactions. To reproduce # Fig. 2 of [Del Maestro and Gingras](https://arxiv.org/abs/cond-mat/0403494), # an empirical rescaling of energy is necessary. -fudge_factor = 1/2 + +fudge_factor = 1/2 # For purposes of reproducing prior work fig = Figure(size=(900,300)) From 9e25b61cf9ec7d209063fbd2e1215b8712402aa0 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 17/29] Docs tweaks --- examples/07_Dipole_Dipole.jl | 22 ++++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/examples/07_Dipole_Dipole.jl b/examples/07_Dipole_Dipole.jl index 3eb62f7d0..3da64ac2b 100644 --- a/examples/07_Dipole_Dipole.jl +++ b/examples/07_Dipole_Dipole.jl @@ -7,15 +7,21 @@ using Sunny, GLMakie -# Create a Pyrochlore lattice from spacegroup 227 in the standard setting. +# Create a Pyrochlore crystal from spacegroup 227. + latvecs = lattice_vectors(10.19, 10.19, 10.19, 90, 90, 90) -cryst = Crystal(latvecs, [[0,0,0]], 227, setting="2") +positions = [[0, 0, 0]] +cryst = Crystal(latvecs, positions, 227, setting="2") sys = System(cryst, (1, 1, 1), [SpinInfo(1, S=7/2, g=2)], :dipole, seed=2) J1 = 0.304 * meV_per_K set_exchange!(sys, J1, Bond(1,2,[0,0,0])) +# Reshape to the primitive cell with four atoms. To facilitate indexing, the +# function [`position_to_site`](@ref) accepts positions with respect to the +# original (cubic) cell. + shape = [1/2 1/2 0; 0 1/2 1/2; 1/2 0 1/2] sys_prim = reshape_supercell(sys, shape) @@ -24,7 +30,11 @@ set_dipole!(sys_prim, [-1, +1, 0], position_to_site(sys_prim, [1/4,1/4,0])) set_dipole!(sys_prim, [+1, +1, 0], position_to_site(sys_prim, [1/4,0,1/4])) set_dipole!(sys_prim, [-1, -1, 0], position_to_site(sys_prim, [0,1/4,1/4])) -plot_spins(sys_prim) +plot_spins(sys_prim; ghost_radius=8) + +# Calculate dispersions with and without long-range dipole interactions. The +# high-symmetry k-points are specified with respect to the conventional cubic +# cell. q_points = [[0,0,0], [0,1,0], [1,1/2,0], [1/2,1/2,1/2], [3/4,3/4,0], [0,0,0]] q_labels = ["Γ","X","W","L","K","Γ"] @@ -39,9 +49,9 @@ enable_dipole_dipole!(sys_prim) swt = SpinWaveTheory(sys_prim, energy_ϵ=1e-6) disp2 = dispersion(swt, path); -# Plot dispersions with and without long-range dipole interactions. To reproduce -# Fig. 2 of [Del Maestro and Gingras](https://arxiv.org/abs/cond-mat/0403494), -# an empirical rescaling of energy is necessary. +# To reproduce Fig. 2 of [Del Maestro and +# Gingras](https://arxiv.org/abs/cond-mat/0403494), an empirical rescaling of +# energy was found to be necessary. fudge_factor = 1/2 # For purposes of reproducing prior work From 985f4a3544d4eb5c5aedc6df7e9d509d3319850d Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 18/29] Update note --- docs/lyx_notes/fourier_lswt.lyx | 79 ++++++++++++++++++++++++--------- 1 file changed, 59 insertions(+), 20 deletions(-) diff --git a/docs/lyx_notes/fourier_lswt.lyx b/docs/lyx_notes/fourier_lswt.lyx index fe09d7c3e..192f3ee1e 100644 --- a/docs/lyx_notes/fourier_lswt.lyx +++ b/docs/lyx_notes/fourier_lswt.lyx @@ -364,17 +364,38 @@ Fourier-transform of dipole-dipole interactions \end_layout \begin_layout Standard -The pairwise interaction between two point magnetic dipoles is given by - the tensor +Consider two magnetic moment dipoles +\begin_inset Formula $\boldsymbol{\mu}_{1}$ +\end_inset + + and +\begin_inset Formula $\boldsymbol{\mu}_{2}$ +\end_inset + + separated by a displacement vector +\begin_inset Formula $\mathbf{r}$ +\end_inset + +. + Their pairwise interaction energy is +\begin_inset Formula $\boldsymbol{\mu}_{1}\cdot\mathcal{A}(\mathbf{r})\boldsymbol{\mu}_{2}$ +\end_inset + +, where \begin_inset Formula -\begin{align*} -\mathcal{D}_{\alpha\beta}(\mathbf{R}_{ij\mathbf{n}}) & =-\frac{\mu_{0}}{4\pi}\nabla_{\alpha}\nabla_{\beta}\frac{1}{|\mathbf{R}_{ij\mathbf{n}}|}\\ - & =\frac{\mu_{0}}{4\pi}\left(\delta_{\alpha\beta}\frac{1}{|\mathbf{R}_{ij\mathbf{n}}|^{3}}-3\frac{R_{ij\mathbf{n}}^{\alpha}R_{ij\mathbf{n}}^{\beta}}{|\mathbf{R}_{ij\mathbf{n}}|^{5}}\right). -\end{align*} +\begin{align} +\mathcal{A}_{\alpha\beta}(\mathbf{r}) & =-\frac{\mu_{0}}{4\pi}\nabla_{\alpha}\nabla_{\beta}\frac{1}{|\mathbf{r}|}\nonumber \\ + & =\frac{\mu_{0}}{4\pi}\left(\delta_{\alpha\beta}\frac{1}{|\mathbf{r}|^{3}}-3\frac{r^{\alpha}r^{\beta}}{|\mathbf{r}|^{5}}\right). +\end{align} \end_inset -Note + +\end_layout + +\begin_layout Standard +On the lattice, we will be interested in a discrete set of displacement + vectors \begin_inset Formula $\mathbf{R}_{ij\mathbf{n}}=\mathbf{r}_{ij}+\mathbf{r}_{\mathbf{n}}$ \end_inset @@ -392,17 +413,10 @@ Note \end_inset denote lattice vectors for the magnetic cell. -\end_layout - -\begin_layout Standard -Fourier transform yields a dipole-dipole equivalent to -\begin_inset Formula $J_{(\alpha i)(\beta j)}(\mathbf{q})$ -\end_inset - - + Fourier transform of the interaction tensor yields, \begin_inset Formula \begin{equation} -\mathcal{D}_{(\alpha i)(\beta j)}(\mathbf{q})=\sum_{\mathbf{n}}'e^{i2\pi\mathbf{q}\cdot\mathbf{n}}\mathcal{D}_{\alpha\beta}(\mathbf{R}_{ij\mathbf{n}}).\label{eq:dipole_sum} +\mathcal{A}_{(\alpha i)(\beta j)}(\mathbf{q})=\sum_{\mathbf{n}}'e^{i2\pi\mathbf{q}\cdot\mathbf{n}}\mathcal{A}_{\alpha\beta}(\mathbf{R}_{ij\mathbf{n}}).\label{eq:dipole_sum} \end{equation} \end_inset @@ -432,6 +446,31 @@ The notation in different orders. \end_layout +\begin_layout Standard +Magnetic moments are linearly related to spin angular momenta, +\begin_inset Formula $\boldsymbol{\mu}=-\mu_{B}g\mathbf{S}$ +\end_inset + +. + This implies an interaction tensor for spins, +\end_layout + +\begin_layout Standard +\begin_inset Formula +\begin{equation} +J_{ij}(\mathbf{q})=\mu_{B}^{2}g_{i}^{T}\mathcal{A}_{ij}(\mathbf{q})g_{j},\label{eq:dipole_sum-1} +\end{equation} + +\end_inset + +with matrix-multiplication acting as contraction on the implicit vector + indices +\begin_inset Formula $\alpha\beta$ +\end_inset + +. +\end_layout + \begin_layout Standard Following the standard Ewald procedure, the infinite summation over cells @@ -444,7 +483,7 @@ Following the standard Ewald procedure, the infinite summation over cells \begin_layout Standard \begin_inset Formula \begin{equation} -\mathcal{D}_{(\alpha i)(\beta j)}(\mathbf{q})=\mathcal{D}_{(\alpha i)(\beta j)}^{\mathrm{LR}}(\mathbf{q},\lambda)+\mathcal{D}_{(\alpha i)(\beta j)}^{\mathrm{SR}}(\mathbf{q},\lambda)-\mathcal{D}_{(\alpha i)(\beta j)}^{\mathrm{self}}(\lambda). +\mathcal{A}_{(\alpha i)(\beta j)}(\mathbf{q})=\mathcal{A}_{(\alpha i)(\beta j)}^{\mathrm{LR}}(\mathbf{q},\lambda)+\mathcal{A}_{(\alpha i)(\beta j)}^{\mathrm{SR}}(\mathbf{q},\lambda)-\mathcal{A}_{(\alpha i)(\beta j)}^{\mathrm{self}}(\lambda). \end{equation} \end_inset @@ -465,7 +504,7 @@ This decomposition introduces an arbitrary parameter The interactions at long-range are best handled in Fourier-space, \begin_inset Formula \begin{equation} -\mathcal{D}_{(\alpha i)(\beta j)}^{\mathrm{LR}}(\mathbf{q},\lambda)=\frac{\mu_{0}}{V}\sum_{\mathbf{m}}\mathcal{K}_{\alpha\beta}(\mathbf{k}_{\mathbf{m}+\mathbf{q}})e^{-|\mathbf{k}_{\mathbf{m}+\mathbf{q}}|^{2}/4\lambda^{2}}e^{-i\mathbf{k}_{\mathbf{m}+\mathbf{q}}\cdot\mathbf{r}_{ij}},\label{eq:D_LR} +\mathcal{A}_{(\alpha i)(\beta j)}^{\mathrm{LR}}(\mathbf{q},\lambda)=\frac{\mu_{0}}{V}\sum_{\mathbf{m}}\mathcal{K}_{\alpha\beta}(\mathbf{k}_{\mathbf{m}+\mathbf{q}})e^{-|\mathbf{k}_{\mathbf{m}+\mathbf{q}}|^{2}/4\lambda^{2}}e^{-i\mathbf{k}_{\mathbf{m}+\mathbf{q}}\cdot\mathbf{r}_{ij}},\label{eq:D_LR} \end{equation} \end_inset @@ -543,7 +582,7 @@ https://doi.org/10.1088/0953-8984/16/43/r02 The short-range part of the interactions can be handled directly in real-space, \begin_inset Formula \begin{equation} -\mathcal{D}_{(\alpha i)(\beta j)}^{\mathrm{SR}}(\mathbf{q},\lambda)=\frac{\mu_{0}}{4\pi}\sum_{\mathbf{n}}'\mathcal{U}_{\alpha\beta}(\mathbf{R}_{ij\mathbf{n}},\lambda)e^{i2\pi\mathbf{q}\cdot\mathbf{n}}, +\mathcal{A}_{(\alpha i)(\beta j)}^{\mathrm{SR}}(\mathbf{q},\lambda)=\frac{\mu_{0}}{4\pi}\sum_{\mathbf{n}}'\mathcal{U}_{\alpha\beta}(\mathbf{R}_{ij\mathbf{n}},\lambda)e^{i2\pi\mathbf{q}\cdot\mathbf{n}}, \end{equation} \end_inset @@ -561,7 +600,7 @@ where Finally, one must subtract the term, \begin_inset Formula \begin{align} -\mathcal{D}_{(\alpha i)(\beta j)}^{\mathrm{self}}(\lambda) & =-\frac{\mu_{0}}{2\pi}\delta_{ij}\nabla^{\alpha}\nabla^{\beta}\frac{1}{r}\mathrm{erf}\left(\lambda r\right)|_{r\to0}\nonumber \\ +\mathcal{A}_{(\alpha i)(\beta j)}^{\mathrm{self}}(\lambda) & =-\frac{\mu_{0}}{2\pi}\delta_{ij}\nabla^{\alpha}\nabla^{\beta}\frac{1}{r}\mathrm{erf}\left(\lambda r\right)|_{r\to0}\nonumber \\ & =\mu_{0}\frac{\lambda^{3}}{3\pi^{3/2}}\delta_{ij}\delta_{\alpha\beta}. \end{align} From f2e501268f4982ed1eea3979dced1cd7369f457c Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 19/29] accum_dipole_dipole_locally! --- src/System/Ewald.jl | 29 ++++++++++++++++++ src/System/PairExchange.jl | 62 ++++++++++++++++++++++++++++++++------ 2 files changed, 82 insertions(+), 9 deletions(-) diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl index 06ae79654..42871d112 100644 --- a/src/System/Ewald.jl +++ b/src/System/Ewald.jl @@ -229,3 +229,32 @@ function ewald_energy_delta(sys::System{N}, site, s::Vec3) where N ∇E = ewald_grad_at(sys, site) return Δs⋅∇E + dot(Δμ, ewald.A[1, 1, 1, i, i], Δμ) / 2 end + + +# Adds dipole-dipole interactions on top of existing exchange interactions. +function accum_dipole_dipole_locally!(sys::System{N}, cutoff) where N + (; gs) = sys + (; μB, μ0) = sys.units + + if !isnothing(sys.origin) + accum_dipole_dipole_locally!(sys.origin, cutoff) + transfer_interactions!(sys, sys.origin) + return + end + + is_homogeneous(sys) || error("Currently requires homogeneous system") + ints = interactions_homog(sys) + + for bond in reference_bonds(sys.crystal, cutoff) + for i in 1:natoms(sys.crystal) + for bond′ in all_symmetry_related_bonds_for_atom(sys.crystal, i, bond) + (; j) = bond′ + r = global_displacement(sys.crystal, bond′) + iszero(r) && continue + r̂ = normalize(r) + bilin = (μ0/4π) * μB^2 * gs[i]' * ((I - 3r̂⊗r̂) / norm(r)^3) * gs[j] + push_coupling!(ints[i].pair, bond′, 0.0, Mat3(bilin), 0.0, zero(TensorDecomposition); accum=true) + end + end + end +end diff --git a/src/System/PairExchange.jl b/src/System/PairExchange.jl index 3f0f441cd..e777871d6 100644 --- a/src/System/PairExchange.jl +++ b/src/System/PairExchange.jl @@ -30,20 +30,54 @@ function to_float_or_mat3(J) end +function Base.:+(c1::PairCoupling, c2::PairCoupling) + @assert c1.isculled == c2.isculled + @assert c1.bond == c2.bond + + scalar = c1.scalar + c2.scalar + + bilin = if (typeof(c1.bilin) == typeof(c2.bilin)) + c1.bilin + c2.bilin + else + Mat3(c1.bilin*I + c2.bilin*I) + end + + biquad = if (typeof(c1.biquad) == typeof(c2.biquad)) + c1.biquad + c2.biquad + else + Mat5(c1.biquad*I + c2.biquad*I) + end + + tensordec = c1.general + c2.general + + PairCoupling(c1.isculled, c1.bond, scalar, bilin, biquad, tensordec) +end + # Internal function only -function push_coupling!(couplings, bond::Bond, scalar::Float64, bilin::Union{Float64, Mat3}, biquad::Union{Float64, Mat5}, tensordec::TensorDecomposition) - # Remove previous coupling on this bond - filter!(c -> c.bond != bond, couplings) +function push_coupling!(couplings, bond::Bond, scalar::Float64, bilin::Union{Float64, Mat3}, biquad::Union{Float64, Mat5}, tensordec::TensorDecomposition; accum=false) + # Find and remove existing couplings for this bond + idxs = findall(c -> c.bond == bond, couplings) + existing = couplings[idxs] + deleteat!(couplings, idxs) - # If the new coupling is exactly zero, return early - iszero(bilin) && iszero(biquad) && isempty(tensordec.data) && return + # If the new coupling is exactly zero, and we're not accumulating, then + # return early + iszero(bilin) && iszero(biquad) && isempty(tensordec.data) && !accum && return - # Otherwise, add the new coupling to the list + # Create a new PairCoupling isculled = bond_parity(bond) - push!(couplings, PairCoupling(isculled, bond, scalar, bilin, biquad, tensordec)) + coupling = PairCoupling(isculled, bond, scalar, bilin, biquad, tensordec) + + # Optionally accumulate to an existing PairCoupling + if accum && !isempty(existing) + coupling += only(existing) + end - # Sorting after each insertion will introduce quadratic scaling in length of - # `couplings`. In typical usage, the `couplings` list will be short. + # Add to the list and sort by isculled. Sorting after each insertion will + # introduce quadratic scaling in length of `couplings`. If this becomes + # slow, we could swap two PairCouplings instead of performing a full sort. + isculled = bond_parity(bond) + push!(couplings, coupling) sort!(couplings, by=c->c.isculled) return @@ -116,6 +150,16 @@ function Base.zero(::Type{TensorDecomposition}) return TensorDecomposition(gen, gen, []) end +function Base.:+(op1::TensorDecomposition, op2::TensorDecomposition) + @assert op1.gen1 ≈ op2.gen1 + @assert op1.gen2 ≈ op2.gen2 + isempty(op2.data) && return op1 + total = sum(kron(A, B) for (A, B) in vcat(op1.data, op2.data)) + N1 = size(op1.gen1, 1) + N2 = size(op1.gen2, 1) + TensorDecomposition(op1.gen1, op1.gen2, svd_tensor_expansion(total, N1, N2)) +end + function Base.isapprox(op1::TensorDecomposition, op2::TensorDecomposition; kwargs...) isempty(op1.data) == isempty(op2.data) && return true op1′ = sum(kron(A, B) for (A, B) in op1.data) From b877249f685b921c838eba7d69bb1baabcc4d77c Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 20/29] Update example --- examples/07_Dipole_Dipole.jl | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/examples/07_Dipole_Dipole.jl b/examples/07_Dipole_Dipole.jl index 3da64ac2b..a76ae76d0 100644 --- a/examples/07_Dipole_Dipole.jl +++ b/examples/07_Dipole_Dipole.jl @@ -45,15 +45,21 @@ xticks = (xticks[1], q_labels) swt = SpinWaveTheory(sys_prim, energy_ϵ=1e-6) disp1 = dispersion(swt, path) -enable_dipole_dipole!(sys_prim) -swt = SpinWaveTheory(sys_prim, energy_ϵ=1e-6) +sys_prim′ = Sunny.clone_system(sys_prim) +enable_dipole_dipole!(sys_prim′) +swt = SpinWaveTheory(sys_prim′, energy_ϵ=1e-6) disp2 = dispersion(swt, path); +sys_prim′ = Sunny.clone_system(sys_prim) +Sunny.accum_dipole_dipole_locally!(sys_prim′, 5.0) +swt = SpinWaveTheory(sys_prim′, energy_ϵ=1e-6) +disp3 = dispersion(swt, path); + # To reproduce Fig. 2 of [Del Maestro and # Gingras](https://arxiv.org/abs/cond-mat/0403494), an empirical rescaling of -# energy was found to be necessary. +# energy is necessary. -fudge_factor = 1/2 # For purposes of reproducing prior work +fudge_factor = 1/2 # To reproduce prior work fig = Figure(size=(900,300)) @@ -71,4 +77,8 @@ for i in axes(disp2, 2) lines!(ax, 1:length(disp2[:,i]), fudge_factor*disp2[:,i]/meV_per_K) end +for i in axes(disp3, 2) + lines!(ax, 1:length(disp3[:,i]), fudge_factor*disp3[:,i]/meV_per_K; color=:gray, linestyle=:dash) +end + fig From a1c7ce13a87b5d6328205fad467eb8071782de9b Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 21/29] Better example --- examples/07_Dipole_Dipole.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/07_Dipole_Dipole.jl b/examples/07_Dipole_Dipole.jl index a76ae76d0..52c3d3c27 100644 --- a/examples/07_Dipole_Dipole.jl +++ b/examples/07_Dipole_Dipole.jl @@ -30,7 +30,7 @@ set_dipole!(sys_prim, [-1, +1, 0], position_to_site(sys_prim, [1/4,1/4,0])) set_dipole!(sys_prim, [+1, +1, 0], position_to_site(sys_prim, [1/4,0,1/4])) set_dipole!(sys_prim, [-1, -1, 0], position_to_site(sys_prim, [0,1/4,1/4])) -plot_spins(sys_prim; ghost_radius=8) +plot_spins(sys_prim; ghost_radius=8, color=[:red,:blue,:yellow,:purple]) # Calculate dispersions with and without long-range dipole interactions. The # high-symmetry k-points are specified with respect to the conventional cubic @@ -39,7 +39,7 @@ plot_spins(sys_prim; ghost_radius=8) q_points = [[0,0,0], [0,1,0], [1,1/2,0], [1/2,1/2,1/2], [3/4,3/4,0], [0,0,0]] q_labels = ["Γ","X","W","L","K","Γ"] density = 150 -path, xticks = reciprocal_space_path(cryst, q_points, density); +path, xticks = reciprocal_space_path(cryst, q_points, density) xticks = (xticks[1], q_labels) swt = SpinWaveTheory(sys_prim, energy_ϵ=1e-6) From 37315ca8d560361b16b684eebc64db5d74d87fd1 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 22/29] Tweak --- src/System/PairExchange.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/System/PairExchange.jl b/src/System/PairExchange.jl index e777871d6..0aca25988 100644 --- a/src/System/PairExchange.jl +++ b/src/System/PairExchange.jl @@ -36,13 +36,13 @@ function Base.:+(c1::PairCoupling, c2::PairCoupling) scalar = c1.scalar + c2.scalar - bilin = if (typeof(c1.bilin) == typeof(c2.bilin)) + bilin = if typeof(c1.bilin) == typeof(c2.bilin) c1.bilin + c2.bilin else Mat3(c1.bilin*I + c2.bilin*I) end - biquad = if (typeof(c1.biquad) == typeof(c2.biquad)) + biquad = if typeof(c1.biquad) == typeof(c2.biquad) c1.biquad + c2.biquad else Mat5(c1.biquad*I + c2.biquad*I) From f952ef15411ffacb348c9cbb804c9398bb3617ea Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 23/29] Comment on surface dipole term --- src/System/Ewald.jl | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl index 42871d112..aef527f01 100644 --- a/src/System/Ewald.jl +++ b/src/System/Ewald.jl @@ -113,7 +113,14 @@ function precompute_dipole_ewald_aux(cryst::Crystal, latsize::NTuple{3,Int}, μ0 m = Vec3(m1, m2, m3) k = recipvecs * (m + q_reshaped - round.(q_reshaped)) k² = k⋅k - if 0 < k² <= kmax*kmax + + ϵ² = 1e-16 + if k² <= ϵ² + # Consider including a surface dipole term as in S. W. DeLeeuw, + # J. W. Perram, and E. R. Smith, Proc. R. Soc. Lond. A 373, + # 27-56 (1980). For a spherical geometry, this term might be: + # acc += (μ0/2V) * I₃ + elseif ϵ² < k² <= kmax*kmax phase = cis(-k⋅Δr) acc += phase * (μ0/V) * (exp(-σ²*k²/2) / k²) * (k⊗k) end From 4a521aef52587e9e7502689a13766f8b79cdfec0 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 24/29] Remove explicit LSWT shift --- examples/07_Dipole_Dipole.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/07_Dipole_Dipole.jl b/examples/07_Dipole_Dipole.jl index 52c3d3c27..4bde74133 100644 --- a/examples/07_Dipole_Dipole.jl +++ b/examples/07_Dipole_Dipole.jl @@ -42,17 +42,17 @@ density = 150 path, xticks = reciprocal_space_path(cryst, q_points, density) xticks = (xticks[1], q_labels) -swt = SpinWaveTheory(sys_prim, energy_ϵ=1e-6) +swt = SpinWaveTheory(sys_prim) disp1 = dispersion(swt, path) sys_prim′ = Sunny.clone_system(sys_prim) enable_dipole_dipole!(sys_prim′) -swt = SpinWaveTheory(sys_prim′, energy_ϵ=1e-6) +swt = SpinWaveTheory(sys_prim′) disp2 = dispersion(swt, path); sys_prim′ = Sunny.clone_system(sys_prim) -Sunny.accum_dipole_dipole_locally!(sys_prim′, 5.0) -swt = SpinWaveTheory(sys_prim′, energy_ϵ=1e-6) +Sunny.accum_dipole_dipole_locally!(sys_prim′, 5.0) # Experimental function +swt = SpinWaveTheory(sys_prim′) disp3 = dispersion(swt, path); # To reproduce Fig. 2 of [Del Maestro and From 855926e34cf9f8a960d5eeab79151114678491ba Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 25/29] Spaces --- examples/07_Dipole_Dipole.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/07_Dipole_Dipole.jl b/examples/07_Dipole_Dipole.jl index 4bde74133..a7b8f2ce6 100644 --- a/examples/07_Dipole_Dipole.jl +++ b/examples/07_Dipole_Dipole.jl @@ -16,7 +16,7 @@ cryst = Crystal(latvecs, positions, 227, setting="2") sys = System(cryst, (1, 1, 1), [SpinInfo(1, S=7/2, g=2)], :dipole, seed=2) J1 = 0.304 * meV_per_K -set_exchange!(sys, J1, Bond(1,2,[0,0,0])) +set_exchange!(sys, J1, Bond(1, 2, [0,0,0])) # Reshape to the primitive cell with four atoms. To facilitate indexing, the # function [`position_to_site`](@ref) accepts positions with respect to the @@ -25,10 +25,10 @@ set_exchange!(sys, J1, Bond(1,2,[0,0,0])) shape = [1/2 1/2 0; 0 1/2 1/2; 1/2 0 1/2] sys_prim = reshape_supercell(sys, shape) -set_dipole!(sys_prim, [+1, -1, 0], position_to_site(sys_prim, [0,0,0])) -set_dipole!(sys_prim, [-1, +1, 0], position_to_site(sys_prim, [1/4,1/4,0])) -set_dipole!(sys_prim, [+1, +1, 0], position_to_site(sys_prim, [1/4,0,1/4])) -set_dipole!(sys_prim, [-1, -1, 0], position_to_site(sys_prim, [0,1/4,1/4])) +set_dipole!(sys_prim, [+1, -1, 0], position_to_site(sys_prim, [0, 0, 0])) +set_dipole!(sys_prim, [-1, +1, 0], position_to_site(sys_prim, [1/4, 1/4, 0])) +set_dipole!(sys_prim, [+1, +1, 0], position_to_site(sys_prim, [1/4, 0, 1/4])) +set_dipole!(sys_prim, [-1, -1, 0], position_to_site(sys_prim, [0, 1/4, 1/4])) plot_spins(sys_prim; ghost_radius=8, color=[:red,:blue,:yellow,:purple]) From f258b3b50c2fa57f8fa105caa1e37ebbda64b403 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 17:53:39 -0600 Subject: [PATCH 26/29] Update docs --- docs/src/versions.md | 4 +++- examples/07_Dipole_Dipole.jl | 6 +++--- src/Sunny.jl | 3 ++- src/System/Ewald.jl | 25 +++++++++++++++++++++---- 4 files changed, 29 insertions(+), 9 deletions(-) diff --git a/docs/src/versions.md b/docs/src/versions.md index 5f6dd87ac..e37d27f76 100644 --- a/docs/src/versions.md +++ b/docs/src/versions.md @@ -6,7 +6,9 @@ * [`view_crystal`](@ref) called on a [`System`](@ref) now optionally shows spin or magnetic dipoles. * Interactions for [`enable_dipole_dipole!`](@ref) are now supported in linear - spin wave theory, with proper Ewald summation. + spin wave theory, with proper Ewald summation. For a faster alternative, the + experimental function [`modify_exchange_with_truncated_dipole_dipole!`](@ref) + will accept a real-space cutoff. ## v0.5.9 diff --git a/examples/07_Dipole_Dipole.jl b/examples/07_Dipole_Dipole.jl index a7b8f2ce6..28fc6fda8 100644 --- a/examples/07_Dipole_Dipole.jl +++ b/examples/07_Dipole_Dipole.jl @@ -45,13 +45,13 @@ xticks = (xticks[1], q_labels) swt = SpinWaveTheory(sys_prim) disp1 = dispersion(swt, path) -sys_prim′ = Sunny.clone_system(sys_prim) +sys_prim′ = clone_system(sys_prim) enable_dipole_dipole!(sys_prim′) swt = SpinWaveTheory(sys_prim′) disp2 = dispersion(swt, path); -sys_prim′ = Sunny.clone_system(sys_prim) -Sunny.accum_dipole_dipole_locally!(sys_prim′, 5.0) # Experimental function +sys_prim′ = clone_system(sys_prim) +modify_exchange_with_truncated_dipole_dipole!(sys_prim′, 5.0) swt = SpinWaveTheory(sys_prim′) disp3 = dispersion(swt, path); diff --git a/src/Sunny.jl b/src/Sunny.jl index b6c5ef402..f325ca738 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -56,7 +56,8 @@ export SpinInfo, System, Site, clone_system, eachsite, position_to_site, global_ set_coherent!, set_dipole!, polarize_spins!, randomize_spins!, set_spin_rescaling!, energy, energy_per_site, spin_label, set_onsite_coupling!, set_pair_coupling!, set_exchange!, dmvec, enable_dipole_dipole!, set_external_field!, to_inhomogeneous, set_external_field_at!, set_vacancy_at!, set_onsite_coupling_at!, - set_exchange_at!, set_pair_coupling_at!, symmetry_equivalent_bonds, remove_periodicity! + set_exchange_at!, set_pair_coupling_at!, symmetry_equivalent_bonds, remove_periodicity!, + modify_exchange_with_truncated_dipole_dipole! include("MagneticOrdering.jl") export print_wrapped_intensities, suggest_magnetic_supercell, set_spiral_order!, set_spiral_order_on_sublattice! diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl index aef527f01..aff6a35b8 100644 --- a/src/System/Ewald.jl +++ b/src/System/Ewald.jl @@ -237,14 +237,31 @@ function ewald_energy_delta(sys::System{N}, site, s::Vec3) where N return Δs⋅∇E + dot(Δμ, ewald.A[1, 1, 1, i, i], Δμ) / 2 end - -# Adds dipole-dipole interactions on top of existing exchange interactions. -function accum_dipole_dipole_locally!(sys::System{N}, cutoff) where N +""" + modify_exchange_with_truncated_dipole_dipole!(sys::System, cutoff) + +This *experimental* function is subject to change. + +Like [`enable_dipole_dipole!`](@ref), the purpose of this function is to +introduce long-range dipole-dipole interactions between magnetic moments. +Whereas `enable_dipole_dipole!` employs Ewald summation, this function instead +employs real-space pair couplings, with truncation at some user-specified +`cutoff` distance. If the cutoff is relatively small, then this function can be +faster than `enable_dipole_dipole!`, especially for spin wave theory +calculations. + +**Caution**. This function will modify existing bilinear couplings between spins +by adding dipole-dipole interactions. It must therefore be called _after_ all +other pair couplings have been specified. Conversely, any calls to +`set_exchange!`, `set_pair_coupling!`, etc. will irreversibly delete the +dipole-dipole interactions that have been introduced by this function. +""" +function modify_exchange_with_truncated_dipole_dipole!(sys::System{N}, cutoff) where N (; gs) = sys (; μB, μ0) = sys.units if !isnothing(sys.origin) - accum_dipole_dipole_locally!(sys.origin, cutoff) + modify_exchange_with_truncated_dipole_dipole!(sys.origin, cutoff) transfer_interactions!(sys, sys.origin) return end From de9047354b4920463818c3f3118d2876b8a1638e Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 18:17:57 -0600 Subject: [PATCH 27/29] Rename `push_coupling!` to `replace_coupling!` Also streamline calculation of `isculled`. --- src/Reshaping.jl | 3 +-- src/SpinWaveTheory/SpinWaveTheory.jl | 10 ++++---- src/System/Ewald.jl | 2 +- src/System/PairExchange.jl | 36 +++++++++++++--------------- src/System/Types.jl | 4 ++++ 5 files changed, 27 insertions(+), 28 deletions(-) diff --git a/src/Reshaping.jl b/src/Reshaping.jl index 105b5cabd..082c89b95 100644 --- a/src/Reshaping.jl +++ b/src/Reshaping.jl @@ -51,8 +51,7 @@ function transfer_interactions!(sys::System{N}, src::System{N}) where N new_pc = PairCoupling[] for pc in src_int.pair new_bond = map_bond_to_other_crystal(src.crystal, pc.bond, sys.crystal, new_i) - isculled = bond_parity(new_bond) - push!(new_pc, PairCoupling(isculled, new_bond, pc.scalar, pc.bilin, pc.biquad, pc.general)) + push!(new_pc, PairCoupling(new_bond, pc.scalar, pc.bilin, pc.biquad, pc.general)) end new_pc = sort!(new_pc, by=c->c.isculled) new_ints[new_i].pair = new_pc diff --git a/src/SpinWaveTheory/SpinWaveTheory.jl b/src/SpinWaveTheory/SpinWaveTheory.jl index 2b8a58ded..2f4e8b2bb 100644 --- a/src/SpinWaveTheory/SpinWaveTheory.jl +++ b/src/SpinWaveTheory/SpinWaveTheory.jl @@ -72,7 +72,7 @@ end # contains all information about the interaction in the `general` (tensor # decomposition) field. function as_general_pair_coupling(pc, sys) - (; isculled, bond, scalar, bilin, biquad, general) = pc + (; bond, scalar, bilin, biquad, general) = pc N1 = sys.Ns[bond.i] N2 = sys.Ns[bond.j] @@ -99,16 +99,16 @@ function as_general_pair_coupling(pc, sys) # Generate new interaction with extract_parts=false scalar, bilin, biquad, general = decompose_general_coupling(accum, N1, N2; extract_parts=false) - return PairCoupling(isculled, bond, scalar, bilin, biquad, general) + return PairCoupling(bond, scalar, bilin, biquad, general) end function rotate_general_coupling_into_local_frame(pc, U1, U2) - (; isculled, bond, scalar, bilin, biquad, general) = pc + (; bond, scalar, bilin, biquad, general) = pc data_new = map(general.data) do (A, B) (Hermitian(U1'*A*U1), Hermitian(U2'*B*U2)) end td = TensorDecomposition(general.gen1, general.gen2, data_new) - return PairCoupling(isculled, bond, scalar, bilin, biquad, td) + return PairCoupling(bond, scalar, bilin, biquad, td) end # Prepare local operators and observables for SU(N) spin wave calculation by @@ -243,7 +243,7 @@ function swt_data_dipole(sys::System{0}, obs) @assert isempty(general.data) - ints.pair[c] = PairCoupling(isculled, bond, scalar, bilin, biquad, general) + ints.pair[c] = PairCoupling(bond, scalar, bilin, biquad, general) end end diff --git a/src/System/Ewald.jl b/src/System/Ewald.jl index aff6a35b8..4ebd1b5f9 100644 --- a/src/System/Ewald.jl +++ b/src/System/Ewald.jl @@ -277,7 +277,7 @@ function modify_exchange_with_truncated_dipole_dipole!(sys::System{N}, cutoff) w iszero(r) && continue r̂ = normalize(r) bilin = (μ0/4π) * μB^2 * gs[i]' * ((I - 3r̂⊗r̂) / norm(r)^3) * gs[j] - push_coupling!(ints[i].pair, bond′, 0.0, Mat3(bilin), 0.0, zero(TensorDecomposition); accum=true) + replace_coupling!(ints[i].pair, PairCoupling(bond′, 0.0, Mat3(bilin), 0.0, zero(TensorDecomposition)); accum=true) end end end diff --git a/src/System/PairExchange.jl b/src/System/PairExchange.jl index 0aca25988..3ec55962c 100644 --- a/src/System/PairExchange.jl +++ b/src/System/PairExchange.jl @@ -29,6 +29,9 @@ function to_float_or_mat3(J) return J::Union{Float64, Mat3} end +function Base.iszero(c::PairCoupling) + return iszero(c.scalar) && iszero(c.bilin) && iszero(c.biquad) && isempty(c.general.data) +end function Base.:+(c1::PairCoupling, c2::PairCoupling) @assert c1.isculled == c2.isculled @@ -50,23 +53,21 @@ function Base.:+(c1::PairCoupling, c2::PairCoupling) tensordec = c1.general + c2.general - PairCoupling(c1.isculled, c1.bond, scalar, bilin, biquad, tensordec) + PairCoupling(c1.bond, scalar, bilin, biquad, tensordec) end # Internal function only -function push_coupling!(couplings, bond::Bond, scalar::Float64, bilin::Union{Float64, Mat3}, biquad::Union{Float64, Mat5}, tensordec::TensorDecomposition; accum=false) +function replace_coupling!(list, coupling::PairCoupling; accum=false) + (; bond) = coupling + # Find and remove existing couplings for this bond - idxs = findall(c -> c.bond == bond, couplings) - existing = couplings[idxs] - deleteat!(couplings, idxs) + idxs = findall(c -> c.bond == bond, list) + existing = list[idxs] + deleteat!(list, idxs) # If the new coupling is exactly zero, and we're not accumulating, then # return early - iszero(bilin) && iszero(biquad) && isempty(tensordec.data) && !accum && return - - # Create a new PairCoupling - isculled = bond_parity(bond) - coupling = PairCoupling(isculled, bond, scalar, bilin, biquad, tensordec) + iszero(coupling) && !accum && return # Optionally accumulate to an existing PairCoupling if accum && !isempty(existing) @@ -76,17 +77,12 @@ function push_coupling!(couplings, bond::Bond, scalar::Float64, bilin::Union{Flo # Add to the list and sort by isculled. Sorting after each insertion will # introduce quadratic scaling in length of `couplings`. If this becomes # slow, we could swap two PairCouplings instead of performing a full sort. - isculled = bond_parity(bond) - push!(couplings, coupling) - sort!(couplings, by=c->c.isculled) + push!(list, coupling) + sort!(list, by=c->c.isculled) return end -function allapproxequal(a; kwargs...) - mean = sum(a; init=0.0) / length(a) - all(x -> isapprox(mean, x), a) -end # If A ≈ α B, then return the scalar α. Otherwise, return A. function proportionality_factor(A, B; atol=1e-12) @@ -243,7 +239,7 @@ function set_pair_coupling_aux!(sys::System, scalar::Float64, bilin::Union{Float bilin′ = transform_coupling_for_bonds(sys.crystal, bond′, bond, bilin) biquad′ = transform_coupling_for_bonds(sys.crystal, bond′, bond, biquad) tensordec′ = transform_coupling_for_bonds(sys.crystal, bond′, bond, tensordec) - push_coupling!(ints[i].pair, bond′, scalar, bilin′, biquad′, tensordec′) + replace_coupling!(ints[i].pair, PairCoupling(bond′, scalar, bilin′, biquad′, tensordec′)) end end end @@ -424,8 +420,8 @@ function set_pair_coupling_at_aux!(sys::System, scalar::Float64, bilin::Union{Fl site2 = to_cartesian(site2) bond = sites_to_internal_bond(sys, site1, site2, offset) - push_coupling!(ints[site1].pair, bond, scalar, bilin, biquad, tensordec) - push_coupling!(ints[site2].pair, reverse(bond), scalar, bilin', biquad', reverse(tensordec)) + replace_coupling!(ints[site1].pair, PairCoupling(bond, scalar, bilin, biquad, tensordec)) + replace_coupling!(ints[site2].pair, PairCoupling(reverse(bond), scalar, bilin', biquad', reverse(tensordec))) end """ diff --git a/src/System/Types.jl b/src/System/Types.jl index 6567ba32a..2cbba9912 100644 --- a/src/System/Types.jl +++ b/src/System/Types.jl @@ -52,6 +52,10 @@ struct PairCoupling # General pair interactions, only valid in SU(N) mode general :: TensorDecomposition + + function PairCoupling(bond, scalar, bilin, biquad, general) + return new(bond_parity(bond), bond, scalar, bilin, biquad, general) + end end mutable struct Interactions From a6a86d7633c57d87d8a58cc22472a5ae81a82b5b Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 18:35:23 -0600 Subject: [PATCH 28/29] Tweak addition of tensor decomposition --- src/System/PairExchange.jl | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/System/PairExchange.jl b/src/System/PairExchange.jl index 3ec55962c..57d4e2eeb 100644 --- a/src/System/PairExchange.jl +++ b/src/System/PairExchange.jl @@ -51,9 +51,8 @@ function Base.:+(c1::PairCoupling, c2::PairCoupling) Mat5(c1.biquad*I + c2.biquad*I) end - tensordec = c1.general + c2.general - - PairCoupling(c1.bond, scalar, bilin, biquad, tensordec) + general = c1.general + c2.general + PairCoupling(c1.bond, scalar, bilin, biquad, general) end # Internal function only @@ -147,13 +146,21 @@ function Base.zero(::Type{TensorDecomposition}) end function Base.:+(op1::TensorDecomposition, op2::TensorDecomposition) + isempty(op2.data) && return op1 + isempty(op1.data) && return op2 @assert op1.gen1 ≈ op2.gen1 @assert op1.gen2 ≈ op2.gen2 - isempty(op2.data) && return op1 + + # We could re-optimize the SVD decomposition as below, but doing this + # unnecessarily would cost some floating point precision. + return TensorDecomposition(op1.gen1, op1.gen2, vcat(op1.data, op2.data)) + + #= total = sum(kron(A, B) for (A, B) in vcat(op1.data, op2.data)) N1 = size(op1.gen1, 1) N2 = size(op1.gen2, 1) - TensorDecomposition(op1.gen1, op1.gen2, svd_tensor_expansion(total, N1, N2)) + return TensorDecomposition(op1.gen1, op1.gen2, svd_tensor_expansion(total, N1, N2)) + =# end function Base.isapprox(op1::TensorDecomposition, op2::TensorDecomposition; kwargs...) From 507e864dc032c9dcbdb151955db9a6de7fd00027 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Wed, 1 May 2024 18:50:15 -0600 Subject: [PATCH 29/29] Style tweaks --- examples/07_Dipole_Dipole.jl | 16 ++++++++-------- src/SpinWaveTheory/HamiltonianDipole.jl | 4 ++-- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/examples/07_Dipole_Dipole.jl b/examples/07_Dipole_Dipole.jl index 28fc6fda8..ea3002395 100644 --- a/examples/07_Dipole_Dipole.jl +++ b/examples/07_Dipole_Dipole.jl @@ -45,14 +45,14 @@ xticks = (xticks[1], q_labels) swt = SpinWaveTheory(sys_prim) disp1 = dispersion(swt, path) -sys_prim′ = clone_system(sys_prim) -enable_dipole_dipole!(sys_prim′) -swt = SpinWaveTheory(sys_prim′) -disp2 = dispersion(swt, path); - -sys_prim′ = clone_system(sys_prim) -modify_exchange_with_truncated_dipole_dipole!(sys_prim′, 5.0) -swt = SpinWaveTheory(sys_prim′) +sys_prim_dd = clone_system(sys_prim) +enable_dipole_dipole!(sys_prim_dd) +swt = SpinWaveTheory(sys_prim_dd) +disp2 = dispersion(swt, path) + +sys_prim_tdd = clone_system(sys_prim) +modify_exchange_with_truncated_dipole_dipole!(sys_prim_tdd, 5.0) +swt = SpinWaveTheory(sys_prim_tdd) disp3 = dispersion(swt, path); # To reproduce Fig. 2 of [Del Maestro and diff --git a/src/SpinWaveTheory/HamiltonianDipole.jl b/src/SpinWaveTheory/HamiltonianDipole.jl index b1d927992..4ce619ee8 100644 --- a/src/SpinWaveTheory/HamiltonianDipole.jl +++ b/src/SpinWaveTheory/HamiltonianDipole.jl @@ -101,9 +101,9 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r A = precompute_dipole_ewald_aux(sys.crystal, (1,1,1), units.μ0, q_reshaped, cis, Val{ComplexF64}()) A = reshape(A, L, L) - # Interaction matrix for wavevector (0,0,0) + # Interaction matrix for wavevector (0,0,0). It could be recalculated as: + # precompute_dipole_ewald(sys.crystal, (1,1,1), units.μ0) A0 = sys.ewald.A - # @assert A0 ≈ precompute_dipole_ewald(sys.crystal, (1,1,1), units.μ0) A0 = reshape(A0, L, L) # Loop over sublattice pairs