Skip to content

Commit

Permalink
Enlarge dynamical matrix by factor of 2
Browse files Browse the repository at this point in the history 
Simplify interface for multiplication by dynamical matrix.
  • Loading branch information
@kbarros
kbarros committed Jan 9, 2025
1 parent c7c3155 commit f3affd6
Show file tree
Hide file tree
Showing 8 changed files with 132 additions and 152 deletions.
26 changes: 13 additions & 13 deletions src/EntangledUnits/EntangledSpinWaveTheory.jl
View file
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 @@ -377,6 +377,6 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::EntangledSpinWaveTheor

# Add small constant shift for positive-definiteness
for i in 1:2L
H[i,i] += swt.regularization
H[i,i] += 2 * swt.regularization
end
end
end
14 changes: 7 additions & 7 deletions src/SpinWaveTheory/DispersionAndIntensities.jl
View file
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
107 changes: 51 additions & 56 deletions src/SpinWaveTheory/HamiltonianDipole.jl
View file
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,46 +122,41 @@ 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)

# Add small constant shift for positive-definiteness
for i in 1:2L
H[i, i] += swt.regularization
H[i, i] += 2 * swt.regularization
end
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 All @@ -284,7 +279,7 @@ function multiply_by_hamiltonian_dipole!(y::AbstractMatrix{ComplexF64}, x::Abstr
end

# Add small constant shift for positive-definiteness.
@inbounds @. Y += swt.regularization * X
@inbounds @. Y += 2 * swt.regularization * X

nothing
end
Loading

0 comments on commit f3affd6

Please sign in to comment.