diff --git a/docs/src/versions.md b/docs/src/versions.md index 30ef33085..667eae311 100644 --- a/docs/src/versions.md +++ b/docs/src/versions.md @@ -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) diff --git a/src/EntangledUnits/EntangledSpinWaveTheory.jl b/src/EntangledUnits/EntangledSpinWaveTheory.jl index 78c1cb8e9..850b94449 100644 --- a/src/EntangledUnits/EntangledSpinWaveTheory.jl +++ b/src/EntangledUnits/EntangledSpinWaveTheory.jl @@ -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 @@ -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 @@ -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) @@ -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 \ No newline at end of file +end diff --git a/src/SpinWaveTheory/DispersionAndIntensities.jl b/src/SpinWaveTheory/DispersionAndIntensities.jl index 9e6a0c1e5..643a750c8 100644 --- a/src/SpinWaveTheory/DispersionAndIntensities.jl +++ b/src/SpinWaveTheory/DispersionAndIntensities.jl @@ -28,14 +28,14 @@ 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 @@ -43,7 +43,7 @@ function bogoliubov!(T::Matrix{ComplexF64}, H::Matrix{ComplexF64}) #= Ĩ = 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]) diff --git a/src/SpinWaveTheory/HamiltonianDipole.jl b/src/SpinWaveTheory/HamiltonianDipole.jl index 85c5903b3..76f0a1d2a 100644 --- a/src/SpinWaveTheory/HamiltonianDipole.jl +++ b/src/SpinWaveTheory/HamiltonianDipole.jl @@ -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 @@ -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 @@ -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) @@ -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) @@ -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))) @@ -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 @@ -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 @@ -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 diff --git a/src/SpinWaveTheory/HamiltonianSUN.jl b/src/SpinWaveTheory/HamiltonianSUN.jl index f6376de4c..11a14a3fb 100644 --- a/src/SpinWaveTheory/HamiltonianSUN.jl +++ b/src/SpinWaveTheory/HamiltonianSUN.jl @@ -7,7 +7,7 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh N = sys.Ns[1] Na = natoms(sys.crystal) - L = (N-1) * Na + L = (N-1) * Na # Clear the Hamiltonian @assert size(H) == (2L, 2L) @@ -25,7 +25,7 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh 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 @@ -42,23 +42,23 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh # 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) @@ -91,23 +91,23 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh Bj = 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]) + c = (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]) + c = 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] + c = 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] + c = 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] + c = 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[α, β]) @@ -130,12 +130,6 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh end -function multiply_by_hamiltonian_SUN(x::AbstractMatrix{ComplexF64}, 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. @@ -160,7 +154,7 @@ function multiply_by_hamiltonian_SUN!(y::AbstractMatrix{ComplexF64}, x::Abstract # 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]) + c = 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] @@ -180,13 +174,12 @@ function multiply_by_hamiltonian_SUN!(y::AbstractMatrix{ComplexF64}, x::Abstract # 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] + c1 = (A[m,n] - δ(m,n)*A[N,N]) * B[N,N] + c2 = A[N,N] * (B[m,n] - δ(m,n)*B[N,N]) + c3 = A[m,N] * B[N,n] + c4 = A[N,m] * B[n,N] + c5 = A[m,N] * B[n,N] @inbounds for q in axes(Y, 1) Y[q, m, i, 1] += c1 * X[q, n, i, 1] diff --git a/src/SpinWaveTheory/LSWTCorrections.jl b/src/SpinWaveTheory/LSWTCorrections.jl index 7d809fe8d..b74892b18 100644 --- a/src/SpinWaveTheory/LSWTCorrections.jl +++ b/src/SpinWaveTheory/LSWTCorrections.jl @@ -24,7 +24,7 @@ function energy_per_site_lswt_correction(swt::SpinWaveTheory; 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 diff --git a/src/SpinWaveTheory/SpinWaveTheory.jl b/src/SpinWaveTheory/SpinWaveTheory.jl index 1b51ee968..1c11c493d 100644 --- a/src/SpinWaveTheory/SpinWaveTheory.jl +++ b/src/SpinWaveTheory/SpinWaveTheory.jl @@ -95,6 +95,13 @@ function to_reshaped_rlu(sys::System{N}, q) where N return sys.crystal.recipvecs \ (orig_crystal(sys).recipvecs * q) end +function dynamical_matrix(swt::SpinWaveTheory, q_reshaped) + L = nbands(swt) + H = zeros(ComplexF64, 2L, 2L) + dynamical_matrix!(H, swt, q_reshaped) + return Hermitian(H) +end + function dynamical_matrix!(H, swt::SpinWaveTheory, q_reshaped) if swt.sys.mode == :SUN swt_hamiltonian_SUN!(H, swt, q_reshaped) @@ -104,14 +111,22 @@ function dynamical_matrix!(H, swt::SpinWaveTheory, q_reshaped) end end +function mul_dynamical_matrix(swt, x, qs_reshaped) + y = zero(x) + mul_dynamical_matrix!(swt, y, x, qs_reshaped) +end + function mul_dynamical_matrix!(swt, y, x, qs_reshaped) + @assert size(x) == size(y) "Dimensions of x and y do not match" + @assert size(x, 1) == length(qs_reshaped) + if swt.sys.mode == :SUN multiply_by_hamiltonian_SUN!(y, x, swt, qs_reshaped) else multiply_by_hamiltonian_dipole!(y, x, swt, qs_reshaped) end - # TODO: Incorporate this factor into Hamiltonian itself - y .*= 2 + + return y end diff --git a/src/Spiral/SpinWaveTheorySpiral.jl b/src/Spiral/SpinWaveTheorySpiral.jl index 3768c99b9..4e0989d80 100644 --- a/src/Spiral/SpinWaveTheorySpiral.jl +++ b/src/Spiral/SpinWaveTheorySpiral.jl @@ -51,7 +51,7 @@ function swt_hamiltonian_dipole_spiral!(H::Matrix{ComplexF64}, sswt::SpinWaveThe for ints in sys.interactions_union for c in ints.pair - (; i, j, n) = c.bond + (; i, j, n) = c.bond θ = (2*π * dot(k,n)) Rn = axis_angle_to_matrix(axis, θ) @@ -84,12 +84,10 @@ function swt_hamiltonian_dipole_spiral!(H::Matrix{ComplexF64}, sswt::SpinWaveThe end end - @. H /= 2 - # Add Zeeman term for i in 1:L B = sys.extfield[1, 1, 1, i]' * sys.gs[1, 1, 1, i] - B′ = - (B * local_rotations[i][:, 3]) / 2 + B′ = - B * local_rotations[i][:, 3] H[i, i] += B′ H[i+L, i+L] += conj(B′) end @@ -98,10 +96,10 @@ function swt_hamiltonian_dipole_spiral!(H::Matrix{ComplexF64}, sswt::SpinWaveThe for i in 1:L s = sqrtS[i]^2 (; c2, c4, c6) = stevens_coefs[i] - H[i, i] += -3s*c2[3] - 40*s^3*c4[5] - 168*s^5*c6[7] - H[i+L, i+L] += -3s*c2[3] - 40*s^3*c4[5] - 168*s^5*c6[7] - H[i, i+L] += +im*(s*c2[5] + 6s^3*c4[7] + 16s^5*c6[9]) + (s*c2[1] + 6s^3*c4[3] + 16s^5*c6[5]) - H[i+L, i] += -im*(s*c2[5] + 6s^3*c4[7] + 16s^5*c6[9]) + (s*c2[1] + 6s^3*c4[3] + 16s^5*c6[5]) + H[i, i] += -6s*c2[3] - 80*s^3*c4[5] - 336*s^5*c6[7] + H[i+L, i+L] += -6s*c2[3] - 80*s^3*c4[5] - 336*s^5*c6[7] + H[i, i+L] += 2s*(c2[1]+im*c2[5]) + 12s^3*(c4[3]+im*c4[7]) + 32s^5*(c6[5]+im*c6[9]) + H[i+L, i] += 2s*(c2[1]-im*c2[5]) + 12s^3*(c4[3]-im*c4[7]) + 32s^5*(c6[5]-im*c6[9]) end isnothing(sys.ewald) || error("Ewald interactions not yet supported") diff --git a/test/test_lswt.jl b/test/test_lswt.jl index c8fdb98ee..5d0cc48c8 100644 --- a/test/test_lswt.jl +++ b/test/test_lswt.jl @@ -50,12 +50,12 @@ swt = SpinWaveTheory(sys; measure=ssf_perp(sys; apply_g=false)) res = intensities_bands(swt, qs) - # println(round.(res.disp; digits=12)) - # println(round.(res.data; digits=12)) - disps_golden = [1394.440092589898 1393.728009961772 1393.008551261241 1392.919524984054 1279.23991907861 1279.094568482213 1278.22451852538 1277.691761492894 1194.366336265056 1193.750083635324 1191.583519669771 1189.794451350786 1131.422439597908 1131.202770084574 1065.242927860663 1065.095892456642 1026.649340932941 1024.028348568104 1022.830406309672 1020.767349649408 945.202397540707 944.795817861582 835.545028404179 832.001588705254 827.939501419137 827.307586957877 821.216582176438 820.430993577092 820.294548786087 818.594571008014 810.207001100209 808.553158283705 766.524411081004 766.516102760229 766.513825862473 766.508655568197 758.579854177684 754.683765895715 750.572578901226 750.471006262524 665.954573018229 662.421047663209 651.46556256004 651.417940134413 581.258189162714 568.105209810117 559.053702306466 558.493005833015 552.043762746846 550.131096080954 539.733572957862 530.698033203026 499.661483520139 494.928560833195 435.233706072008 427.70227707436 408.128705863823 399.856401759966 370.069343073402 369.845327696313 365.049514250289 363.639416679443 354.648012601512 346.609483937092 341.98916517756 339.373361078069 318.363717394716 276.219249213429 263.1610538422 257.409506256945 230.539454204132 229.778324183075 203.971681289205 197.504237163938 193.879371544746 189.866421885131 189.815806977935 167.944134441876 154.923566511974 146.21953885867] - data_golden = natoms * [9.6662565e-5 0.0 0.001807884262 0.0 0.002166250472 0.0 0.003835132693 0.0 0.0 0.013550066967 0.018327409336 0.0 0.001319003523 0.0 0.0 0.006677115951 0.0 0.015583193495 0.028006929969 0.0 0.007782414773 0.0 0.028892276057 0.0 0.0 0.009748743725 0.0 0.007788123595 0.007278285719 0.0 0.0 0.00116773942 0.000190731866 0.000235139044 0.0 0.0 0.002193317995 0.0 0.008315770402 0.0 0.016520647062 0.0 0.0 0.014989172648 0.08363192972 0.001902190356 0.0 0.0 0.0 0.017111545439 0.007341641476 0.0 0.0 0.151419922476 0.094666992635 0.0 0.214639628105 0.0 0.0 0.072514184691 0.087279329707 0.0 0.0 0.23891458934 0.0 0.401546184942 0.0 0.050427766439 0.04544662783 0.0 0.0 0.078427781361 0.192350495231 0.0 0.002594420277 0.0 0.037800182163 0.081921558936 0.0 0.0] - @test isapprox(res.disp, disps_golden'; atol=1e-9) - @test isapprox(res.data, data_golden'; atol=1e-9) + # println(round.(vec(res.disp); digits=12)) + # println(round.(vec(res.data); digits=12)) + disps_golden = [1394.440092579932, 1393.728009951748, 1393.008551251227, 1392.919524974106, 1279.239919068656, 1279.094568472208, 1278.224518515369, 1277.691761482934, 1194.366336255048, 1193.750083625344, 1191.583519659752, 1189.794451340786, 1131.422439587814, 1131.202770074483, 1065.242927850645, 1065.095892446623, 1026.649340922942, 1024.028348558074, 1022.830406299683, 1020.767349639389, 945.202397530637, 944.795817851532, 835.545028394177, 832.001588695239, 827.939501409159, 827.307586947865, 821.216582166477, 820.430993567146, 820.294548776124, 818.594570998049, 810.207001090183, 808.553158273676, 766.524411071072, 766.516102750272, 766.51382585253, 766.508655558251, 758.5798541677, 754.683765885704, 750.572578891224, 750.47100625256, 665.954573008178, 662.421047653195, 651.465562549975, 651.417940124351, 581.258189152584, 568.105209800095, 559.053702296455, 558.493005822971, 552.043762736843, 550.131096070956, 539.733572947825, 530.698033192909, 499.661483510074, 494.928560823138, 435.233706061892, 427.702277064325, 408.128705853663, 399.856401749648, 370.069343063308, 369.845327686246, 365.049514240266, 363.639416669427, 354.648012591404, 346.609483927002, 341.989165167266, 339.373361067994, 318.363717384335, 276.219249203178, 263.1610538318, 257.409506246766, 230.539454193854, 229.778324172693, 203.971681278992, 197.504237153931, 193.879371534726, 189.866421875022, 189.815806967694, 167.94413443161, 154.923566498384, 146.21953884777] + data_golden = [0.00038665026, 0.0, 0.007231537047, 0.0, 0.008665001888, 0.0, 0.015340530772, 0.0, 0.0, 0.054200267869, 0.073309637341, 0.0, 0.005276014091, 0.0, 0.0, 0.026708463804, 0.0, 0.062332773978, 0.112027719874, 0.0, 0.031129659091, 0.0, 0.115569104228, 0.0, 0.0, 0.038994974901, 0.0, 0.031152494381, 0.029113142876, 0.0, 0.0, 0.00467095768, 0.000762927464, 0.000940556177, 0.0, 0.0, 0.00877327198, 0.0, 0.033263081608, 0.0, 0.066082588247, 0.0, 0.0, 0.059956690592, 0.334527718881, 0.007608761425, 0.0, 0.0, 0.0, 0.068446181755, 0.029366565904, 0.0, 0.0, 0.605679689903, 0.378667970539, 0.0, 0.858558512424, 0.0, 0.0, 0.290056738762, 0.349117318827, 0.0, 0.0, 0.955658357376, 0.0, 1.606184739756, 0.0, 0.201711065754, 0.181786511323, 0.0, 0.0, 0.313711125445, 0.769401980923, 0.0, 0.010377681107, 0.0, 0.151200728665, 0.327686235743, 0.0, 0.0] + @test isapprox(res.disp, disps_golden; atol=1e-9) + @test isapprox(res.data, data_golden; atol=1e-9) # Test first 5 output matrices formfactors = [1 => FormFactor("Fe2")] @@ -260,7 +260,7 @@ end # P1 point group for most general single-ion anisotropy cryst = Crystal(latvecs, positions, 1) - s = 2 + s = 3 sys_dip = System(cryst, [1 => Moment(; s, g=-1)], :dipole) sys_SUN = System(cryst, [1 => Moment(; s, g=-1)], :SUN) @@ -329,8 +329,8 @@ end disp1 = dispersion(swt1, [q]) disp2 = dispersion(swt2, [q]) - @test disp1 ≈ [1.3236778213351734, 0.9206655419611791] - @test disp2 ≈ [0.9206655419611772] + @test disp1 ≈ [1.3236778041378718, 0.9206655245623444] + @test disp2 ≈ [0.9206655245623366] end end @@ -533,13 +533,7 @@ end @testitem "Equivalence of dense and sparse Hamiltonian constructions" begin - import LinearAlgebra: diagm - - function onehot(i, n) - out = zeros(ComplexF64, 1, n) - out[1, i] = 1.0 - return out - end + using LinearAlgebra # Build System and SpinWaveTheory with exchange, field and single-site anisotropy function simple_swt(mode) @@ -558,35 +552,20 @@ end return SpinWaveTheory(sys; measure=nothing) end - # Construct dipole Hamiltonian standard way - swt = simple_swt(:dipole) - H1 = zeros(ComplexF64, 16, 16) - q = Sunny.Vec3([0.5,0,0]) - Sunny.swt_hamiltonian_dipole!(H1, swt, q) - - # Construct dipole Hamiltonian from sparse matrix-vector multiplies - H2 = zero(H1) - L = Sunny.natoms(swt.sys.crystal) - for i in 1:2L - H2[:, i] = Sunny.multiply_by_hamiltonian_dipole(onehot(i, 2L), swt, [q]) - end + q = Sunny.Vec3(0.5, 0, 0) - @test isapprox(H1, H2; atol=1e-12) + for mode = (:dipole, :SUN) + # Construct Hamiltonian directly + swt = simple_swt(mode) + H1 = Sunny.dynamical_matrix(swt, q) - # Construct SU(N) Hamiltonian standard way - swt = simple_swt(:SUN) - N = swt.sys.Ns[1] - L = (N-1) * Sunny.natoms(swt.sys.crystal) - H1 = zeros(ComplexF64, 2L, 2L) - Sunny.swt_hamiltonian_SUN!(H1, swt, q) + # Construct Hamiltonian by sparse matrix-vector multiplies + L = Sunny.nbands(swt) + x = Matrix{ComplexF64}(I, 2L, 2L) + H2 = transpose(Sunny.mul_dynamical_matrix(swt, x, fill(q, 2L))) - # Construct SU(N) Hamiltonian from sparse matrix-vector multiplies - H2 = zero(H1) - for i in 1:2L - H2[:, i] = Sunny.multiply_by_hamiltonian_SUN(onehot(i, 2L), swt, [q]) + @test H1 ≈ H2 end - - @test isapprox(H1, H2; atol=1e-12) end @@ -630,12 +609,12 @@ end swt = SpinWaveTheory(sys; measure) qs = [[0,0,0], [0.5,0,0], [1,0,0]] - energies = 0:0.1:5 + energies = 0:0.5:5 kernel = lorentzian(fwhm=0.3) res = intensities(swt, qs; energies, kernel) - data_ref = [0.042551643530614226 0.05061618774203612 0.06118972834822155 0.07541981401191966 0.09518339306332815 0.12371078181463013 0.1669136576105103 0.23644634156679015 0.3574143016120729 0.589319004913063 1.0795750388213625 2.0570916947398903 2.6569442073245604 1.6749367428248116 0.8709898932656854 0.4927038808971025 0.30849624153704924 0.20903618563778026 0.1502266820186879 0.11286975075833217 0.08777050885259896 0.07013929292833648 0.057300827672198455 0.047672208009393216 0.04027090049655778 0.0344619791000266 0.02982086993597283 0.02605519995256082 0.022958441606058772 0.020381422739929135 0.01821426297223774 0.016374601990583864 0.014799738769001156 0.013441266708211906 0.01226133976759786 0.01123002734145714 0.010323410070056573 0.009522188815697288 0.008810654796662341 0.008175917662895354 0.007607320304517221 0.007095990541228231 0.006634494316422953 0.006216564975067639 0.005836890143659928 0.00549094262866035 0.005174845247781923 0.0048852620341421574 0.004619310095626147 0.004374487768740016 0.00414861571474443; 0.14853118271457438 0.19360218524421224 0.2621797495216911 0.3732135013522173 0.5678610986270936 0.944407657666259 1.7450675448755357 3.2945578443571577 3.9947874291087895 2.418790606759786 1.2612861617795654 0.7201697967624341 0.45442367736983386 0.30970051154679945 0.22353509755861065 0.1685055354270234 0.1313752454162916 0.10520387618639207 0.0860938574888819 0.0717287252840356 0.060665219355115055 0.05196772774655005 0.045008912709399315 0.03935576182418519 0.03470177881993486 0.030825189142763776 0.027562374631854493 0.024790523018178495 0.02241601889071837 0.020366506674885078 0.018585357752204528 0.017027745220739847 0.015657814448345492 0.014446613655329999 0.013370560099471452 0.01241028925551664 0.011549781566933818 0.010775692876605122 0.010076836041215703 0.009443775967574439 0.008868510590520351 0.008344217576743833 0.007865051732026552 0.007425981842385972 0.007022658419610481 0.006651305841367626 0.006308633878300496 0.005991764727412878 0.005698172523177506 0.0054256329470820505 0.005172181054614204; 0.042551643530614205 0.0506161877420361 0.06118972834822153 0.07541981401191963 0.0951833930633281 0.12371078181463005 0.16691365761051016 0.23644634156678992 0.3574143016120724 0.5893190049130621 1.0795750388213603 2.0570916947398867 2.656944207324563 1.6749367428248163 0.8709898932656877 0.4927038808971036 0.30849624153704985 0.20903618563778062 0.15022668201868813 0.1128697507583323 0.08777050885259907 0.07013929292833657 0.05730082767219852 0.04767220800939327 0.04027090049655782 0.03446197910002663 0.02982086993597286 0.026055199952560844 0.02295844160605879 0.020381422739929156 0.018214262972237757 0.016374601990583878 0.014799738769001168 0.013441266708211913 0.01226133976759787 0.011230027341457147 0.010323410070056582 0.009522188815697295 0.008810654796662348 0.00817591766289536 0.007607320304517226 0.007095990541228234 0.006634494316422957 0.006216564975067644 0.005836890143659932 0.005490942628660353 0.005174845247781927 0.00488526203414216 0.00461931009562615 0.0043744877687400185 0.0041486157147444325] - - @test res.data ≈ data_ref' + # println(round.(res.data; digits=10)) + data_ref = [0.0425516442 0.148531187 0.0425516442; 0.1237107851 0.9444077155 0.1237107851; 1.0795751088 1.2612860859 1.0795751088; 0.4927038544 0.1685055314 0.4927038544; 0.0877705066 0.0606652185 0.0877705066; 0.0344619785 0.0308251889 0.0344619785; 0.0182142627 0.0185853576 0.0182142627; 0.0112300272 0.0124102892 0.0112300272; 0.0076073202 0.0088685106 0.0076073202; 0.0054909426 0.0066513058 0.0054909426; 0.0041486157 0.005172181 0.0041486157] + @test res.data ≈ data_ref end