From c31c611c1e55484d80508a6242979a0e2c4cdf51 Mon Sep 17 00:00:00 2001 From: ddahlbom Date: Thu, 19 Oct 2023 17:23:35 -0400 Subject: [PATCH] Code Deduplication II: Progress --- src/SpinWaveTheory/HamiltonianSUN.jl | 114 ++++++++------------------- src/SpinWaveTheory/SpinWaveTheory.jl | 28 ++++--- src/SpinWaveTheory/Util.jl | 22 +++++- 3 files changed, 72 insertions(+), 92 deletions(-) diff --git a/src/SpinWaveTheory/HamiltonianSUN.jl b/src/SpinWaveTheory/HamiltonianSUN.jl index 323ae555a..7976441dd 100644 --- a/src/SpinWaveTheory/HamiltonianSUN.jl +++ b/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,23 +43,23 @@ 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 @@ -65,60 +67,6 @@ function swt_pair_coupling!(H, swt, phase, metric, bond, Ti, Tj) 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,14 +83,14 @@ 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 @@ -150,24 +98,32 @@ function swt_hamiltonian_SUN!(H::Matrix{ComplexF64}, swt::SpinWaveTheory, q_resh 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) diff --git a/src/SpinWaveTheory/SpinWaveTheory.jl b/src/SpinWaveTheory/SpinWaveTheory.jl index a9d9bc828..f3109f6df 100644 --- a/src/SpinWaveTheory/SpinWaveTheory.jl +++ b/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,17 +143,21 @@ 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 @@ -161,7 +165,7 @@ function swt_data_sun(sys::System{N}, obs) where N dipole_operators_localized, quadrupole_operators_localized, onsite_operator_localized, - general_pair_operators_localized, + bond_operator_pairs_localized, observables_localized, ) end diff --git a/src/SpinWaveTheory/Util.jl b/src/SpinWaveTheory/Util.jl index a514746cf..309b74fc1 100644 --- a/src/SpinWaveTheory/Util.jl +++ b/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 \ No newline at end of file