Code Deduplication II: Progress

src/SpinWaveTheory/HamiltonianSUN.jl

@@ -1,8 +1,7 @@
 # Construct portion of Hamiltonian due to onsite terms (single-site anisotropy
 # or external field).
-function swt_onsite_coupling!(H, swt, op, site)
-    (; sys) = swt
-
+function swt_onsite_coupling!(H, op, swt, atom)
+    sys = swt.sys
     N = sys.Ns[1] 
     nflavors = N - 1 
     L = nflavors * natoms(sys.crystal)   
@@ -11,22 +10,25 @@ function swt_onsite_coupling!(H, swt, op, site)
 
     for m = 2:N
         for n = 2:N
-            ix_m = (site-1)*nflavors+m-1
-            ix_n = (site-1)*nflavors+n-1
+            m_idx = (atom-1)*nflavors+m-1
+            n_idx = (atom-1)*nflavors+n-1
             c = 0.5 * (op[m, n] - δ(m, n) * op[1, 1])
-            H11[ix_m, ix_n] += c
-            H22[ix_n, ix_m] += c
+            H11[m_idx, n_idx] += c
+            H22[n_idx, m_idx] += c
         end
     end
 end
 
-# Construct portion of Hamiltonian due to pair couplings (bilinear interactions
-# in spins or quadrupoles)
-function swt_pair_coupling!(H, swt, phase, metric, bond, Ti, Tj)
-    # Get views of Hamiltonian submatrices
-    N = swt.sys.Ns[1] 
+# Adds any interaction of the form,
+#
+# J_{ij}^{αβ} T_i^α T_j^β
+#
+# to the spin wave Hamiltonian H. TODO: Complete description.
+function swt_pair_coupling!(H, J, Ti, Tj, swt, phase, bond)
+    sys = swt.sys
+    N = sys.Ns[1] 
     nflavors = N - 1 
-    L = nflavors * natoms(swt.sys.crystal)   
+    L = nflavors * natoms(sys.crystal)   
     H11 = view(H, 1:L, 1:L)
     H12 = view(H, 1:L, L+1:2L)
     H22 = view(H, L+1:2L, L+1:2L)
@@ -41,84 +43,30 @@ function swt_pair_coupling!(H, swt, phase, metric, bond, Ti, Tj)
             i_n = sub_i_M1*nflavors+nM1
             j_n = sub_j_M1*nflavors+nM1
 
-            c = 0.5 * dot_no_conj(Ti[m,n] - δ(m,n)*Ti[1,1], metric, Tj[1,1])
+            c = 0.5 * dot_no_conj(Ti[m,n] - δ(m,n)*Ti[1,1], J, Tj[1,1])
             H11[i_m, i_n] += c
             H22[i_n, i_m] += c
 
-            c = 0.5 * dot_no_conj(Ti[1,1], metric, Tj[m,n] - δ(m,n)*Tj[1,1])
+            c = 0.5 * dot_no_conj(Ti[1,1], J, Tj[m,n] - δ(m,n)*Tj[1,1])
             H11[j_m, j_n] += c
             H22[j_n, j_m] += c
 
-            c = 0.5 * dot_no_conj(Ti[m,1], metric, Tj[1,n])
+            c = 0.5 * dot_no_conj(Ti[m,1], J, Tj[1,n])
             H11[i_m, j_n] += c * phase
             H22[j_n, i_m] += c * conj(phase)
 
-            c = 0.5 * dot_no_conj(Ti[1,m], metric, Tj[n,1])
+            c = 0.5 * dot_no_conj(Ti[1,m], J, Tj[n,1])
             H11[j_n, i_m] += c * conj(phase)
             H22[i_m, j_n] += c * phase
 
-            c = 0.5 * dot_no_conj(Ti[m,1], metric, Tj[n,1])
+            c = 0.5 * dot_no_conj(Ti[m,1], J, Tj[n,1])
             H12[i_m, j_n] += c * phase
             H12[j_n, i_m] += c * conj(phase)
         end
     end
 end
 
 
-# Add generalized couplings to Hamiltonian.
-# DD TODO: Reorganize data for general_pair_operators so can eliminate this
-# function.
-function swt_general_couplings!(H, swt, q)
-    (; sys, data) = swt
-    (; general_pair_operators) = data
-
-    # Get views of Hamiltonian submatrices
-    N = sys.Ns[1] 
-    nflavors = N - 1 
-    L = nflavors * natoms(sys.crystal)   
-    H11 = view(H, 1:L, 1:L)
-    H12 = view(H, 1:L, L+1:2L)
-    H22 = view(H, L+1:2L, L+1:2L)
-
-    for general_pair in general_pair_operators
-        ((A, B), bond) = general_pair
-        phase = exp(2π*im * dot(q, bond.n)) 
-
-        sub_i_M1, sub_j_M1 = bond.i - 1, bond.j - 1
-        for m = 2:N
-            mM1 = m - 1
-            i_m = (sub_i_M1 * nflavors) + mM1
-            j_m = (sub_j_M1 * nflavors) + mM1
-
-            for n = 2:N
-                nM1 = n - 1
-                i_n = (sub_i_M1 * nflavors) + nM1
-                j_n = (sub_j_M1 * nflavors) + nM1
-
-                c = 0.5 * (A[m,n] - δ(m, n)*A[1,1])*B[1,1]
-                H11[i_m, i_n] += c
-                H22[i_n, i_m] += c
-
-                c = 0.5 * A[1,1] * (B[m,n] - δ(m, n)*B[1,1]) 
-                H11[j_m, j_n] += c
-                H22[j_n, j_m] += c
-
-                c = 0.5 * A[m,1] * B[1,n] 
-                H11[i_m, j_n] += c * phase
-                H22[j_n, i_m] += c * conj(phase)
-
-                c = 0.5 * A[1,m] * B[n,1] 
-                H11[j_n, i_m] += c * conj(phase)
-                H22[i_m, j_n] += c * phase
-
-                c = 0.5 * A[m,1] * B[n,1]
-                H12[i_m, j_n] += c * phase
-                H12[j_n, i_m] += c * conj(phase)
-            end
-        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, data) = swt
@@ -135,39 +83,47 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh
     # Add single-site terms (single-site anisotropy and external field)
     for atom = 1:natoms(sys.crystal)
         site_aniso = view(onsite_operator, :, :, atom)
-        swt_onsite_coupling!(H, swt, site_aniso, atom)
+        swt_onsite_coupling!(H, site_aniso, swt, atom)
     end
 
     # Add pair interactions that use explicit bases
     for ints in sys.interactions_union
         for coupling in ints.pair
             # Extract information common to bond
-            (; isculled, bond, bilin, biquad, general) = coupling
+            (; isculled, bond, bilin, biquad) = coupling
             isculled && break
             phase = exp(2π*im * dot(q_reshaped, bond.n)) # Phase associated with periodic wrapping
 
             # Add bilinear interactions
             if !all(iszero, bilin)
                 Si = reinterpret(reshape, Sunny.CVec{3}, view(dipole_operators, :, :, :, bond.i))
                 Sj = reinterpret(reshape, Sunny.CVec{3}, view(dipole_operators, :, :, :, bond.j))
-                J = Mat3(bilin*I)  
+                # J = Mat3(bilin*I)  
 
-                swt_pair_coupling!(H, swt, phase, J, bond, Si, Sj)
+                swt_pair_coupling!(H, bilin, Si, Sj, swt, phase, bond)
             end
 
             # Add biquadratic interactions
             if !all(iszero, biquad)
                 Qi = reinterpret(reshape, Sunny.CVec{5}, view(quadrupole_operators, :, :, :, bond.i))
                 Qj = reinterpret(reshape, Sunny.CVec{5}, view(quadrupole_operators, :, :, :, bond.j))
-                J =  isa(biquad, Float64) ? biquad * scalar_biquad_metric_mat : biquad
+                J =  isa(biquad, Float64) ? biquad * scalar_biquad_metric : biquad
 
-                swt_pair_coupling!(H, swt, phase, J, bond, Qi, Qj)
+                swt_pair_coupling!(H, J, Qi, Qj, swt, phase, bond)
             end
         end
     end
 
     # Add generalized pair interactions
-    swt_general_couplings!(H, swt, q_reshaped)
+    for (bond, operator_pair) in data.bond_operator_pairs 
+        phase = exp(2π*im * dot(q_reshaped, bond.n))
+        veclen = size(operator_pair, 1) 
+        A = reinterpret(reshape, CVec{veclen}, view(operator_pair, :, :, :, 1))
+        B = reinterpret(reshape, CVec{veclen}, view(operator_pair, :, :, :, 2))
+        J = I(veclen)
+
+        swt_pair_coupling!(H, 1.0, A, B, swt, phase, bond)
+    end
 
     # Infer H21 by H=H'
     set_H21!(H)

src/SpinWaveTheory/SpinWaveTheory.jl

@@ -7,11 +7,11 @@ struct SWTDataDipole
 end
 
 struct SWTDataSUN
-    dipole_operators        :: Array{ComplexF64, 4}
-    quadrupole_operators    :: Array{ComplexF64, 4}
-    onsite_operator         :: Array{ComplexF64, 3}  
-    general_pair_operators  :: Vector{Tuple{Tuple{HermitianC64, HermitianC64}, Bond}}
-    observable_operators    :: Array{ComplexF64, 4}
+    dipole_operators      :: Array{ComplexF64, 4}
+    quadrupole_operators  :: Array{ComplexF64, 4}
+    onsite_operator       :: Array{ComplexF64, 3}  
+    bond_operator_pairs   :: Vector{Tuple{Bond, Array{ComplexF64, 4}}}
+    observable_operators  :: Array{ComplexF64, 4}
 end
 
 """
@@ -143,25 +143,29 @@ function swt_data_sun(sys::System{N}, obs) where N
 
     # Rotate operators generated by the tensor decomposition of generalized
     # interactions.
-    general_pair_operators_localized = Tuple{Tuple{HermitianC64, HermitianC64}, Bond}[]
+    bond_operator_pairs_localized = Tuple{Bond, Array{ComplexF64, 4}}[]
     for int in sys.interactions_union
         for pc in int.pair
             (; isculled, bond, general) = pc
-            isculled && continue 
-            atomi, atomj = bond.i, bond.j
-            Ui, Uj = local_quantization_bases[atomi], local_quantization_bases[atomj] 
-            for (A, B) in general.data
+            (isculled || length(general.data) == 0) && continue 
+
+            Ui, Uj = local_quantization_bases[bond.i], local_quantization_bases[bond.j] 
+            nops = length(general.data)
+            bond_operators = zeros(ComplexF64, nops, N, N, 2)
+            for (n, (A, B)) in enumerate(general.data)
                 A′, B′ = Hermitian.((Ui' * A * Ui, Uj' * B * Uj)) 
-                push!(general_pair_operators_localized, ((A′, B′), bond))
+                bond_operators[n,:,:,1] = A′
+                bond_operators[n,:,:,2] = B′
             end
+            push!(bond_operator_pairs_localized, (bond, bond_operators))
         end
     end
 
     return SWTDataSUN(
         dipole_operators_localized, 
         quadrupole_operators_localized, 
         onsite_operator_localized,
-        general_pair_operators_localized,
+        bond_operator_pairs_localized,
         observables_localized,
     )
 end

src/SpinWaveTheory/Util.jl

@@ -34,7 +34,7 @@ function hermiticity_norm(H)
 end
 
 # Modified from LinearAlgebra.jl to not perform any conjugation
-function dot_no_conj(x,A,y)
+function dot_no_conj(x, A, y)
     (axes(x)..., axes(y)...) == axes(A) || throw(DimensionMismatch())
     T = typeof(dot(first(x), first(A), first(y)))
     s = zero(T)
@@ -52,3 +52,23 @@ function dot_no_conj(x,A,y)
     end
     return s
 end
+
+# Constant "metric"
+function dot_no_conj(x, A::Float64, y)
+    s = zero(eltype(x))
+    axes(x) == axes(y) || throw(DimensionMismatch())
+    for j in eachindex(x)
+        s += x[j]*y[j]
+    end
+    return A*s
+end
+
+# Diagonal "metric"
+function dot_no_conj(x, A::SVector{N, Float64}, y) where N
+    s = eltype(x)
+    axes(x) == axes(y)  == A || throw(DimensionMismatch())
+    for j in eachindex(x)
+        s += A[j]*x[i]*y[j]
+    end
+    return s
+end