From 6a89b0616c15e803d7a4844c95c8f382ee9bc68c Mon Sep 17 00:00:00 2001 From: Hao-Phys Date: Sat, 28 Oct 2023 08:30:24 -0600 Subject: [PATCH 1/2] Implemented gen-biquad interaction in dipole mode --- src/SpinWaveTheory/HamiltonianDipole.jl | 62 +++++++++++-------------- src/SpinWaveTheory/Util.jl | 15 ++++++ 2 files changed, 41 insertions(+), 36 deletions(-) diff --git a/src/SpinWaveTheory/HamiltonianDipole.jl b/src/SpinWaveTheory/HamiltonianDipole.jl index 01b9acc1a..5244607e5 100644 --- a/src/SpinWaveTheory/HamiltonianDipole.jl +++ b/src/SpinWaveTheory/HamiltonianDipole.jl @@ -60,42 +60,32 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r ### Biquadratic exchange - (coupling.biquad isa Number) || error("General biquadratic interactions not yet implemented in LSWT.") - J = coupling.biquad::Float64 - - # ⟨Ω₂, Ω₁|[(𝐒₁⋅𝐒₂)^2 + 𝐒₁⋅𝐒₂/2]|Ω₁, Ω₂⟩ = (Ω₁⋅Ω₂)^2 - # The biquadratic part - Ri = R_mat[sub_i] - Rj = R_mat[sub_j] - Rʳ = Ri' * Rj - C0 = Rʳ[3, 3]*S^2 - C1 = S*√S/2*(Rʳ[1, 3] + im*Rʳ[2, 3]) - C2 = S*√S/2*(Rʳ[3, 1] + im*Rʳ[3, 2]) - A11 = -Rʳ[3, 3]*S - A22 = -Rʳ[3, 3]*S - A21 = S/2*(Rʳ[1, 1] - im*Rʳ[1, 2] - im*Rʳ[2, 1] + Rʳ[2, 2]) - A12 = S/2*(Rʳ[1, 1] + im*Rʳ[1, 2] + im*Rʳ[2, 1] + Rʳ[2, 2]) - B21 = S/4*(Rʳ[1, 1] + im*Rʳ[1, 2] + im*Rʳ[2, 1] - Rʳ[2, 2]) - B12 = B21 - - H[sub_i, sub_i] += J* (C0*A11 + C1*conj(C1)) - H[sub_j, sub_j] += J* (C0*A22 + C2*conj(C2)) - H[sub_i, sub_j] += J* ((C0*A12 + C1*conj(C2)) * phase) - H[sub_j, sub_i] += J* ((C0*A21 + C2*conj(C1)) * conj(phase)) - H[sub_i+L, sub_i+L] += J* (C0*A11 + C1*conj(C1)) - H[sub_j+L, sub_j+L] += J* (C0*A22 + C2*conj(C2)) - H[sub_j+L, sub_i+L] += J* ((C0*A12 + C1*conj(C2)) * conj(phase)) - H[sub_i+L, sub_j+L] += J* ((C0*A21 + C2*conj(C1)) * phase) - - H[sub_i, sub_i+L] += J* (C1*conj(C1)) - H[sub_j, sub_j+L] += J* (C2*conj(C2)) - H[sub_i+L, sub_i] += J* (C1*conj(C1)) - H[sub_j+L, sub_j] += J* (C2*conj(C2)) - - H[sub_i, sub_j+L] += J* ((2C0*B12 + C1*C2) * phase) - H[sub_j, sub_i+L] += J* ((2C0*B21 + C2*C1) * conj(phase)) - H[sub_i+L, sub_j] += J* (conj(2C0*B12 + C1*C2) * phase) - H[sub_j+L, sub_i] += J* (conj(2C0*B21 + C2*C1) * conj(phase)) + biquad = coupling.biquad + J = biquad isa Number ? biquad * Mat5(diagm(scalar_biquad_metric)) : biquad::Mat5 + J = S^3 * transform_coupling_for_general_biquad(J, R_mat_i, R_mat_j) + + H[sub_i, sub_i] += -6J[3, 3] + H[sub_j, sub_j] += -6J[3, 3] + H[sub_i+L, sub_i+L] += -6J[3, 3] + H[sub_j+L, sub_j+L] += -6J[3, 3] + H[sub_i+L, sub_i] += 12*(J[1, 3] - 1im*J[5, 3]) + H[sub_i, sub_i+L] += 12*(J[1, 3] + 1im*J[5, 3]) + H[sub_j+L, sub_j] += 12*(J[3, 1] - 1im*J[3, 5]) + H[sub_j, sub_j+L] += 12*(J[3, 1] + 1im*J[3, 5]) + + P = 0.25 * (-J[4, 4]+J[2, 2] - 1im*( J[4, 2]+J[2, 4])) + Q = 0.25 * ( J[4, 4]+J[2, 2] - 1im*(-J[4, 2]+J[2, 4])) + + H[sub_i, sub_j] += Q * phase + H[sub_j, sub_i] += conj(Q) * conj(phase) + H[sub_i+L, sub_j+L] += conj(Q) * phase + H[sub_j+L, sub_i+L] += Q * conj(phase) + + H[sub_i+L, sub_j] += P * phase + H[sub_j+L, sub_i] += P * conj(phase) + H[sub_i, sub_j+L] += conj(P) * phase + H[sub_j, sub_i+L] += conj(P) * conj(phase) + end end diff --git a/src/SpinWaveTheory/Util.jl b/src/SpinWaveTheory/Util.jl index 321036309..9d86a20b1 100644 --- a/src/SpinWaveTheory/Util.jl +++ b/src/SpinWaveTheory/Util.jl @@ -59,4 +59,19 @@ function dot_no_conj(x, A::SVector{N, Float64}, y) where N s += A[j]*x[j]*y[j] end return s +end + +function transform_coupling_for_general_biquad(biquad::Mat5, Ri::Mat3, Rj::Mat3) + k = 2 + + # Notice that in our implementation of (renormalized) :dipole mode for the LSWT, we have the convention that S = RS̃. As a result, we should pass R' below. + Di = unitary_irrep_for_rotation(Ri'; N=2k+1) + Dj = unitary_irrep_for_rotation(Rj'; N=2k+1) + + Vi = stevens_α[k] * transpose(Di) * stevens_αinv[k] + Vj = stevens_α[k] * transpose(Dj) * stevens_αinv[k] + + ret = Vi' * biquad * Vj + @assert norm(imag(ret)) < 1e-12 + return Mat5(real(ret)) end \ No newline at end of file From 3d1a1af33fa69f57a0b9e2a60484d34e692c0f70 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Sat, 28 Oct 2023 11:40:11 -0600 Subject: [PATCH 2/2] Tweaks --- src/Operators/Stevens.jl | 26 +++--- src/SpinWaveTheory/HamiltonianDipole.jl | 102 +++++++++++++----------- src/SpinWaveTheory/Util.jl | 15 ---- src/Symmetry/AllowedCouplings.jl | 23 ++---- test/test_operators.jl | 8 ++ 5 files changed, 86 insertions(+), 88 deletions(-) diff --git a/src/Operators/Stevens.jl b/src/Operators/Stevens.jl index c260ee673..b84431e5b 100644 --- a/src/Operators/Stevens.jl +++ b/src/Operators/Stevens.jl @@ -152,19 +152,25 @@ function matrix_to_stevens_coefficients(A::HermitianC64) end # Spherical tensors T_q rotate as T_q -> D*_{q,q′} T_q′, where D = exp(-i θ n⋅S) -# in dimension 2k+1 irrep. The Stevens operators 𝒪_q are linearly related to -# T_q via 𝒪 = α T, and therefore rotate as 𝒪 -> α conj(D) α⁻¹ 𝒪. -# -# Consider now an operator expansion 𝒜 = cᵀ 𝒪. This operator rotates as 𝒜 -> -# cᵀ α conj(D) α⁻¹ 𝒪 = c′ᵀ 𝒪. The rotated Stevens coefficients must therefore -# satisfy c′ = α⁻ᵀ D† αᵀ c. +# in dimension 2k+1 irrep, for axis-angle (n, θ). The Stevens operators 𝒪_q are +# linearly related to T_q via 𝒪 = α T, and therefore rotate as 𝒪 -> V 𝒪, +# where V = α conj(D) α⁻¹. +function operator_for_stevens_rotation(k, R) + D = unitary_irrep_for_rotation(R; N=2k+1) + V = stevens_α[k] * conj(D) * stevens_αinv[k] + @assert norm(imag(V)) < 1e-12 + return real(V) +end + +# Let c denote coefficients of an operator expansion 𝒜 = c† 𝒪. Under the +# rotation R, Stevens operators transform as 𝒪 → V 𝒪. Alternatively, we can +# treat the Stevens operators as fixed, provided the coefficients transform as +# c† → c† V, or c → V† c. function rotate_stevens_coefficients(c, R::Mat3) N = length(c) k = Int((N-1)/2) - D = unitary_irrep_for_rotation(R; N) - c′ = transpose(stevens_αinv[k]) * D' * transpose(stevens_α[k]) * c - @assert norm(imag(c′)) < 1e-12 - return real(c′) + V = operator_for_stevens_rotation(k, R) + return V' * c end diff --git a/src/SpinWaveTheory/HamiltonianDipole.jl b/src/SpinWaveTheory/HamiltonianDipole.jl index 5244607e5..af591ed13 100644 --- a/src/SpinWaveTheory/HamiltonianDipole.jl +++ b/src/SpinWaveTheory/HamiltonianDipole.jl @@ -10,7 +10,7 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r N = sys.Ns[1] # Dimension of SU(N) coherent states S = (N-1)/2 # Spin magnitude - L = natoms(sys.crystal) # Number of quasiparticle bands + L = natoms(sys.crystal) # Number of quasiparticle bands @assert size(H) == (2L, 2L) @@ -32,60 +32,68 @@ function swt_hamiltonian_dipole!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_r (; isculled, bond) = coupling isculled && break - J = Mat3(coupling.bilin*I) - sub_i, sub_j = bond.i, bond.j - phase = exp(2π*im * dot(q_reshaped, bond.n)) # Phase associated with periodic wrapping + if !iszero(coupling.bilin) + J = Mat3(coupling.bilin*I) + sub_i, sub_j = bond.i, bond.j + phase = exp(2π*im * dot(q_reshaped, bond.n)) # Phase associated with periodic wrapping - R_mat_i = R_mat[sub_i] - R_mat_j = R_mat[sub_j] - Rij = S * (R_mat_i' * J * R_mat_j) + # Transform couplings according to the local quantization basis + R_mat_i = R_mat[sub_i] + R_mat_j = R_mat[sub_j] + Rij = S * (R_mat_i' * J * R_mat_j) - P = 0.25 * (Rij[1, 1] - Rij[2, 2] - im*Rij[1, 2] - im*Rij[2, 1]) - Q = 0.25 * (Rij[1, 1] + Rij[2, 2] - im*Rij[1, 2] + im*Rij[2, 1]) + P = 0.25 * (Rij[1, 1] - Rij[2, 2] - im*Rij[1, 2] - im*Rij[2, 1]) + Q = 0.25 * (Rij[1, 1] + Rij[2, 2] - im*Rij[1, 2] + im*Rij[2, 1]) - H[sub_i, sub_j] += Q * phase - H[sub_j, sub_i] += conj(Q) * conj(phase) - H[sub_i+L, sub_j+L] += conj(Q) * phase - H[sub_j+L, sub_i+L] += Q * conj(phase) + H[sub_i, sub_j] += Q * phase + H[sub_j, sub_i] += conj(Q) * conj(phase) + H[sub_i+L, sub_j+L] += conj(Q) * phase + H[sub_j+L, sub_i+L] += Q * conj(phase) - H[sub_i+L, sub_j] += P * phase - H[sub_j+L, sub_i] += P * conj(phase) - H[sub_i, sub_j+L] += conj(P) * phase - H[sub_j, sub_i+L] += conj(P) * conj(phase) + H[sub_i+L, sub_j] += P * phase + H[sub_j+L, sub_i] += P * conj(phase) + H[sub_i, sub_j+L] += conj(P) * phase + H[sub_j, sub_i+L] += conj(P) * conj(phase) - H[sub_i, sub_i] -= 0.5 * Rij[3, 3] - H[sub_j, sub_j] -= 0.5 * Rij[3, 3] - H[sub_i+L, sub_i+L] -= 0.5 * Rij[3, 3] - H[sub_j+L, sub_j+L] -= 0.5 * Rij[3, 3] + H[sub_i, sub_i] -= 0.5 * Rij[3, 3] + H[sub_j, sub_j] -= 0.5 * Rij[3, 3] + H[sub_i+L, sub_i+L] -= 0.5 * Rij[3, 3] + H[sub_j+L, sub_j+L] -= 0.5 * Rij[3, 3] + end ### Biquadratic exchange - biquad = coupling.biquad - J = biquad isa Number ? biquad * Mat5(diagm(scalar_biquad_metric)) : biquad::Mat5 - J = S^3 * transform_coupling_for_general_biquad(J, R_mat_i, R_mat_j) - - H[sub_i, sub_i] += -6J[3, 3] - H[sub_j, sub_j] += -6J[3, 3] - H[sub_i+L, sub_i+L] += -6J[3, 3] - H[sub_j+L, sub_j+L] += -6J[3, 3] - H[sub_i+L, sub_i] += 12*(J[1, 3] - 1im*J[5, 3]) - H[sub_i, sub_i+L] += 12*(J[1, 3] + 1im*J[5, 3]) - H[sub_j+L, sub_j] += 12*(J[3, 1] - 1im*J[3, 5]) - H[sub_j, sub_j+L] += 12*(J[3, 1] + 1im*J[3, 5]) - - P = 0.25 * (-J[4, 4]+J[2, 2] - 1im*( J[4, 2]+J[2, 4])) - Q = 0.25 * ( J[4, 4]+J[2, 2] - 1im*(-J[4, 2]+J[2, 4])) - - H[sub_i, sub_j] += Q * phase - H[sub_j, sub_i] += conj(Q) * conj(phase) - H[sub_i+L, sub_j+L] += conj(Q) * phase - H[sub_j+L, sub_i+L] += Q * conj(phase) - - H[sub_i+L, sub_j] += P * phase - H[sub_j+L, sub_i] += P * conj(phase) - H[sub_i, sub_j+L] += conj(P) * phase - H[sub_j, sub_i+L] += conj(P) * conj(phase) - + if !iszero(coupling.biquad) + J = coupling.biquad + J = Mat5(J isa Number ? J * diagm(scalar_biquad_metric) : J) + + # Transform couplings according to the local quantization basis + Vi = operator_for_stevens_rotation(2, R_mat_i) + Vj = operator_for_stevens_rotation(2, R_mat_j) + J = S^3 * Mat5(Vi' * J * Vj) + + H[sub_i, sub_i] += -6J[3, 3] + H[sub_j, sub_j] += -6J[3, 3] + H[sub_i+L, sub_i+L] += -6J[3, 3] + H[sub_j+L, sub_j+L] += -6J[3, 3] + H[sub_i+L, sub_i] += 12*(J[1, 3] - 1im*J[5, 3]) + H[sub_i, sub_i+L] += 12*(J[1, 3] + 1im*J[5, 3]) + H[sub_j+L, sub_j] += 12*(J[3, 1] - 1im*J[3, 5]) + H[sub_j, sub_j+L] += 12*(J[3, 1] + 1im*J[3, 5]) + + P = 0.25 * (-J[4, 4]+J[2, 2] - 1im*( J[4, 2]+J[2, 4])) + Q = 0.25 * ( J[4, 4]+J[2, 2] - 1im*(-J[4, 2]+J[2, 4])) + + H[sub_i, sub_j] += Q * phase + H[sub_j, sub_i] += conj(Q) * conj(phase) + H[sub_i+L, sub_j+L] += conj(Q) * phase + H[sub_j+L, sub_i+L] += Q * conj(phase) + + H[sub_i+L, sub_j] += P * phase + H[sub_j+L, sub_i] += P * conj(phase) + H[sub_i, sub_j+L] += conj(P) * phase + H[sub_j, sub_i+L] += conj(P) * conj(phase) + end end end diff --git a/src/SpinWaveTheory/Util.jl b/src/SpinWaveTheory/Util.jl index 9d86a20b1..1315b6e8b 100644 --- a/src/SpinWaveTheory/Util.jl +++ b/src/SpinWaveTheory/Util.jl @@ -60,18 +60,3 @@ function dot_no_conj(x, A::SVector{N, Float64}, y) where N end return s end - -function transform_coupling_for_general_biquad(biquad::Mat5, Ri::Mat3, Rj::Mat3) - k = 2 - - # Notice that in our implementation of (renormalized) :dipole mode for the LSWT, we have the convention that S = RS̃. As a result, we should pass R' below. - Di = unitary_irrep_for_rotation(Ri'; N=2k+1) - Dj = unitary_irrep_for_rotation(Rj'; N=2k+1) - - Vi = stevens_α[k] * transpose(Di) * stevens_αinv[k] - Vj = stevens_α[k] * transpose(Dj) * stevens_αinv[k] - - ret = Vi' * biquad * Vj - @assert norm(imag(ret)) < 1e-12 - return Mat5(real(ret)) -end \ No newline at end of file diff --git a/src/Symmetry/AllowedCouplings.jl b/src/Symmetry/AllowedCouplings.jl index c2a26696d..8f429df94 100644 --- a/src/Symmetry/AllowedCouplings.jl +++ b/src/Symmetry/AllowedCouplings.jl @@ -123,23 +123,14 @@ function transform_coupling_by_symmetry(J::Mat3, R::Mat3, parity) end function transform_coupling_by_symmetry(biquad::Mat5, R::Mat3, parity) + # Under a rotation R, Stevens operators transform as 𝒪 → V 𝒪. To maintain + # `𝒪† biquad 𝒪` as an invariant, the coupling coefficients must transform + # as `biquad -> inv(V)† biquad inv(V)`. As an optimization, note that + # `inv(V(R)) = V(inv(R)) = V(R†)`. k = 2 - D = unitary_irrep_for_rotation(R; N=2k+1) - - # Spherical tensors rotate as `T -> D* T`, involving the complex conjugate - # of the Wigner-D matrix defined above. Stevens operators are `𝒪 = α T`. - # Therefore Stevens operators rotate as `𝒪 -> α D* α⁻¹ 𝒪`. The coupling - # `𝒪† biquad 𝒪` is an invariant, so we impose the transformation rule - # `biquad -> V† biquad V`, where - # - # V = (α D* α⁻¹)⁻¹ = α Dᵀ α⁻¹ - # - # See also `transform_spherical_to_stevens_coefficients`. - - V = stevens_α[k] * transpose(D) * stevens_αinv[k] - ret = V' * (parity ? biquad : biquad') * V - @assert norm(imag(ret)) < 1e-12 - return Mat5(real(ret)) + inv_V = operator_for_stevens_rotation(k, R') + ret = inv_V' * (parity ? biquad : biquad') * inv_V + return Mat5(ret) end # Check whether a coupling matrix J is consistent with symmetries of a bond diff --git a/test/test_operators.jl b/test/test_operators.jl index ef4f63149..84ba2f415 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -199,6 +199,14 @@ end @test R * S ≈ rotate_operator.(S, Ref(R)) end + # Test that Stevens quadrupoles rotate correctly + let + O = stevens_matrices(S₀) + Q = [O[2, q] for q in 2:-1:-2] + V = Sunny.operator_for_stevens_rotation(2, R) + @test V * Q ≈ rotate_operator.(Q, Ref(R)) + end + # Test that Stevens coefficients rotate properly let A = Hermitian(randn(ComplexF64, N, N))