Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change convention for LSWT dynamical matrix #349

Merged
merged 2 commits into from
Jan 12, 2025
Merged

Change convention for LSWT dynamical matrix #349

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f3affd6
Enlarge dynamical matrix by factor of 2
kbarros Dec 17, 2024
0f70257
Reduce LSWT regularization by factor 1/2
kbarros Jan 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions docs/src/versions.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,9 @@
[`print_site`](@ref) accepts an optional reference atom `i_ref`, with default
of `i`. The optional reference bond `b_ref` of [`print_bond`](@ref) now
defaults to `b`.
* The `regularization` parameter in [`SpinWaveTheory`](@ref) is reduced by half,
and now corresponds to an effective energy shift. This may affect intensities,
especially at small excitation energies.

## v0.7.4
(Dec 6, 2024)
Expand Down
24 changes: 12 additions & 12 deletions src/EntangledUnits/EntangledSpinWaveTheory.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -288,7 +288,7 @@ function energy_per_site_lswt_correction(swt::EntangledSpinWaveTheory; opts...)
# The uniform correction to the classical energy (trace of the (1,1)-block
# of the spin-wave Hamiltonian)
dynamical_matrix!(H, swt, zero(Vec3))
δE₁ = -real(tr(view(H, 1:L, 1:L))) / Natoms
δE₁ = -real(tr(view(H, 1:L, 1:L))) / 2Natoms

# Integrate zero-point energy over the first Brillouin zone 𝐪 ∈ [0, 1]³ for
# magnetic cell in reshaped RLU
Expand Down Expand Up @@ -326,7 +326,7 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::EntangledSpinWaveTheor
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])
c = op[m, n] - δ(m, n) * op[N, N]
H11[m, i, n, i] += c
H22[n, i, m, i] += c
end
Expand All @@ -343,23 +343,23 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::EntangledSpinWaveTheor
# 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])
c = (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])

c = 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]

c = 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]

c = 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]

c = 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)
Expand All @@ -379,4 +379,4 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::EntangledSpinWaveTheor
for i in 1:2L
H[i,i] += swt.regularization
end
end
end
14 changes: 7 additions & 7 deletions src/SpinWaveTheory/DispersionAndIntensities.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,22 +28,22 @@ function bogoliubov!(T::Matrix{ComplexF64}, H::Matrix{ComplexF64})
view(T, :, j) .*= c
end

# Inverse of λ are eigenvalues of Ĩ H. A factor of 2 yields the physical
# quasiparticle energies.
# Inverse of λ are eigenvalues of Ĩ H, or equivalently, of √H Ĩ √H.
energies = λ # reuse storage
@. energies = 2 / λ
@. energies = 1 / λ

# The first L elements are positive, while the next L energies are negative.
# Their absolute values are excitation energies for the wavevectors q and
# -q, respectively.
# By Sylvester's theorem, "inertia" (sign signature) is invariant under a
# congruence transform Ĩ → √H Ĩ √H. The first L elements are positive,
# while the next L elements are negative. Their absolute values are
# excitation energies for the wavevectors q and -q, respectively.
@assert all(>(0), view(energies, 1:L)) && all(<(0), view(energies, L+1:2L))

# Disable tests below for speed. Note that the data in H has been
# overwritten by eigen!, so H0 should refer to an original copy of H.
#=
Ĩ = Diagonal([ones(L); -ones(L)])
@assert T' * Ĩ * T ≈ Ĩ
@assert diag(T' * H0 * T) ≈ Ĩ * energies / 2
@assert diag(T' * H0 * T) ≈ Ĩ * energies
# Reflection symmetry H(q) = H(-q) is identified as H11 = conj(H22). In this
# case, eigenvalues come in pairs.
if H0[1:L, 1:L] ≈ conj(H0[L+1:2L, L+1:2L])
Expand Down
103 changes: 49 additions & 54 deletions src/SpinWaveTheory/HamiltonianDipole.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,19 +18,19 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r
for (i, int) in enumerate(sys.interactions_union)
# Zeeman term
B = gs[1, 1, 1, i]' * extfield[1, 1, 1, i]
B′ = - dot(B, local_rotations[i][:, 3]) / 2
B′ = - dot(B, local_rotations[i][:, 3])
H11[i, i] += B′
H22[i, i] += B′

# Single-ion anisotropy
(; c2, c4, c6) = stevens_coefs[i]
s = sqrtS[i]^2
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])
A1 = -6s*c2[3] - 80*s^3*c4[5] - 336*s^5*c6[7]
A2 = 2s*(c2[1]+im*c2[5]) + 12s^3*(c4[3]+im*c4[7]) + 32s^5*(c6[5]+im*c6[9])
H11[i, i] += A1
H22[i, i] += A1
H12[i, i] += A2
H21[i, i] += conj(A2)
H21[i, i] += conj(A2)

# Pair interactions
for coupling in int.pair
Expand All @@ -47,22 +47,22 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r
if !iszero(coupling.bilin)
J = coupling.bilin # Transformed exchange matrix

Q = 0.25 * sij * (J[1, 1] + J[2, 2] - im*(J[1, 2] - J[2, 1]))
Q = 0.5 * sij * (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 * sij * (J[1, 1] - J[2, 2] - im*(J[1, 2] + J[2, 1]))
P = 0.5 * sij * (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
H12[j, i] += conj(P) * conj(phase)

H11[i, i] -= 0.5 * sj * J[3, 3]
H11[j, j] -= 0.5 * si * J[3, 3]
H22[i, i] -= 0.5 * sj * J[3, 3]
H22[j, j] -= 0.5 * si * J[3, 3]
H11[i, i] -= sj * J[3, 3]
H11[j, j] -= si * J[3, 3]
H22[i, i] -= sj * J[3, 3]
H22[j, j] -= si * J[3, 3]
end

# Biquadratic exchange
Expand All @@ -71,22 +71,22 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r

Sj2Si = sj^2 * si
Si2Sj = si^2 * sj
H11[i, i] += -6 * Sj2Si * K[3, 3]
H22[i, i] += -6 * Sj2Si * K[3, 3]
H11[j, j] += -6 * Si2Sj * K[3, 3]
H22[j, j] += -6 * Si2Sj * K[3, 3]
H21[i, i] += 2 * Sj2Si * (K[1, 3] - im*K[5, 3])
H12[i, i] += 2 * Sj2Si * (K[1, 3] + im*K[5, 3])
H21[j, j] += 2 * Si2Sj * (K[3, 1] - im*K[3, 5])
H12[j, j] += 2 * Si2Sj * (K[3, 1] + im*K[3, 5])

Q = 0.25 * sij^3 * ( K[4, 4]+K[2, 2] - im*(-K[4, 2]+K[2, 4]))
H11[i, i] += -12 * Sj2Si * K[3, 3]
H22[i, i] += -12 * Sj2Si * K[3, 3]
H11[j, j] += -12 * Si2Sj * K[3, 3]
H22[j, j] += -12 * Si2Sj * K[3, 3]
H21[i, i] += 4 * Sj2Si * (K[1, 3] - im*K[5, 3])
H12[i, i] += 4 * Sj2Si * (K[1, 3] + im*K[5, 3])
H21[j, j] += 4 * Si2Sj * (K[3, 1] - im*K[3, 5])
H12[j, j] += 4 * Si2Sj * (K[3, 1] + im*K[3, 5])

Q = 0.5 * sij^3 * ( K[4, 4]+K[2, 2] - im*(-K[4, 2]+K[2, 4]))
H11[i, j] += Q * phase
H11[j, i] += conj(Q * phase)
H22[i, j] += conj(Q) * phase
H22[j, i] += Q * conj(phase)

P = 0.25 * sij^3 * (-K[4, 4]+K[2, 2] - im*( K[4, 2]+K[2, 4]))
P = 0.5 * sij^3 * (-K[4, 4]+K[2, 2] - im*( K[4, 2]+K[2, 4]))
H21[i, j] += P * phase
H12[j, i] += conj(P * phase)
H21[j, i] += P * conj(phase)
Expand Down Expand Up @@ -122,31 +122,31 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r
J0 = sqrtS[i]*sqrtS[j] * Rs[i]' * J0 * Rs[j]

# 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]))
Q⁻ = 0.5 * (J[1, 1] + J[2, 2] - im*(J[1, 2] - J[2, 1]))
Q⁺ = 0.5 * (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]))
P⁻ = 0.5 * (J[1, 1] - J[2, 2] - im*(J[1, 2] + J[2, 1]))
P⁺ = 0.5 * (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⁻)

# 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]
H11[i, i] -= J0[3, 3]
H11[j, j] -= J0[3, 3]
H22[i, i] -= J0[3, 3]
H22[j, j] -= J0[3, 3]
end
end

# H must be hermitian up to round-off errors
@assert diffnorm2(H, H') < 1e-12

# Make H exactly hermitian
hermitianpart!(H)

Expand All @@ -157,11 +157,6 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r
end


function multiply_by_hamiltonian_dipole(x::AbstractMatrix{ComplexF64}, 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!(y::AbstractMatrix{ComplexF64}, x::AbstractMatrix{ComplexF64}, swt::SpinWaveTheory, qs_reshaped::Vector{Vec3};
phases=zeros(ComplexF64, size(qs_reshaped)))
Expand All @@ -182,11 +177,11 @@ function multiply_by_hamiltonian_dipole!(y::AbstractMatrix{ComplexF64}, x::Abstr
for i in 1:L
(; c2, c4, c6) = stevens_coefs[i]
s = sqrtS[i]^2
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])
A1 = -6s*c2[3] - 80*s^3*c4[5] - 336*s^5*c6[7]
A2 = 2s*(c2[1]+im*c2[5]) + 12s^3*(c4[3]+im*c4[7]) + 32s^5*(c6[5]+im*c6[9])

B = gs[1, 1, 1, i]' * extfield[1, 1, 1, i]
B′ = - dot(B, view(local_rotations[i], :, 3)) / 2
B′ = - dot(B, view(local_rotations[i], :, 3))

# Seems to be no benefit to breaking this into two loops acting on
# different final indices, presumably because memory access patterns for
Expand Down Expand Up @@ -217,28 +212,28 @@ function multiply_by_hamiltonian_dipole!(y::AbstractMatrix{ComplexF64}, x::Abstr
if !iszero(coupling.bilin)
J = coupling.bilin # This is Rij in previous notation (transformed exchange matrix)

P = 0.25 * sij * (J[1, 1] - J[2, 2] - im*J[1, 2] - im*J[2, 1])
Q = 0.25 * sij * (J[1, 1] + J[2, 2] - im*J[1, 2] + im*J[2, 1])
P = 0.5 * sij * (J[1, 1] - J[2, 2] - im*J[1, 2] - im*J[2, 1])
Q = 0.5 * sij * (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 * phases[q] * X[q, j, 1]
Y[q, i, 1] += conj(P) * phases[q] * X[q, j, 2]
Y[q, i, 1] -= 0.5 * sj * J[3, 3] * X[q, i, 1]
Y[q, i, 1] -= sj * J[3, 3] * X[q, i, 1]
end
@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 * sj * J[3, 3] * X[q, i, 2]
Y[q, i, 2] -= sj * J[3, 3] * X[q, i, 2]
end
@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 * si * J[3, 3] * X[q, j, 1]
Y[q, j, 1] -= si * J[3, 3] * X[q, j, 1]
end
@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 * si * J[3, 3] * X[q, j, 2]
Y[q, j, 2] -= si * J[3, 3] * X[q, j, 2]
end
end

Expand All @@ -248,30 +243,30 @@ function multiply_by_hamiltonian_dipole!(y::AbstractMatrix{ComplexF64}, x::Abstr

Sj2Si = sj^2 * si
Si2Sj = si^2 * sj
Q = 0.25 * sij^3 * ( J[4, 4]+J[2, 2] - im*(-J[4, 2]+J[2, 4]))
P = 0.25 * sij^3 * (-J[4, 4]+J[2, 2] - im*( J[4, 2]+J[2, 4]))
Q = 0.5 * sij^3 * ( J[4, 4]+J[2, 2] - im*(-J[4, 2]+J[2, 4]))
P = 0.5 * sij^3 * (-J[4, 4]+J[2, 2] - im*( J[4, 2]+J[2, 4]))

@inbounds for q in 1:Nq
Y[q, i, 1] += -6 * Sj2Si * J[3, 3] * X[q, i, 1]
Y[q, i, 1] += 2 * Sj2Si * (J[1, 3] + im*J[5, 3]) * X[q, i, 2]
Y[q, i, 1] += -12 * Sj2Si * J[3, 3] * X[q, i, 1]
Y[q, i, 1] += 4 * Sj2Si * (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] += -6 * Sj2Si * J[3, 3] * X[q, i, 2]
Y[q, i, 2] += 2 * Sj2Si * (J[1, 3] - im*J[5, 3]) * X[q, i, 1]
Y[q, i, 2] += -12 * Sj2Si * J[3, 3] * X[q, i, 2]
Y[q, i, 2] += 4 * Sj2Si * (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] += -6 * Si2Sj * J[3, 3] * X[q, j, 1]
Y[q, j, 1] += 2 * Si2Sj * (J[3, 1] + im*J[3, 5]) * X[q, j, 2]
Y[q, j, 1] += -12 * Si2Sj * J[3, 3] * X[q, j, 1]
Y[q, j, 1] += 4 * Si2Sj * (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] += -6 * Si2Sj * J[3, 3] * X[q, j, 2]
Y[q, j, 2] += 2 * Si2Sj * (J[3, 1] - im*J[3, 5]) * X[q, j, 1]
Y[q, j, 2] += -12 * Si2Sj * J[3, 3] * X[q, j, 2]
Y[q, j, 2] += 4 * Si2Sj * (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
Expand Down
Loading
Loading