PairCoupling refactors (#94)

Absorb `set_biquadratic!` into `set_exchange!`

System now stores a single list of PairCoupling objects. System now has `pair` and `onsite` fields, with types `PairCoupling` and `OnsiteCoupling`, respectively.

docs/src/versions.md

@@ -18,6 +18,9 @@ a polynomial of [`spin_operators`](@ref). To understand the mapping between
 these two, the new function [`print_stevens_expansion`](@ref) acts on an
 arbitrary local operator.
 
+Replace `set_biquadratic!` with an optional keyword argument `biquad` to
+[`set_exchange!`](@ref).
+
 Symbolic representations of operators are now hidden unless the package
 `DynamicPolynomials` is explicitly loaded by the user. The functionality of
 `print_anisotropy_as_stevens` has been replaced with
@@ -68,7 +71,7 @@ The function [`to_inhomogeneous`](@ref) creates a system that supports
 inhomogeneous interactions, which can be set using [`set_exchange_at!`](@ref),
 etc.
 
-[`set_biquadratic!`](@ref) replaces `set_exchange_with_biquadratic!`.
+`set_biquadratic!` replaces `set_exchange_with_biquadratic!`.
 
 
 # Version 0.4.0

src/Integrators.jl

@@ -174,15 +174,15 @@ function rhs_langevin!(ΔZ::Array{CVec{N}, 4}, Z::Array{CVec{N}, 4}, ξ::Array{C
     if is_homogeneous(sys)
         ints = interactions_homog(sys)
         for site in all_sites(sys)
-            Λ = ints[to_atom(site)].aniso.matrep
+            Λ = ints[to_atom(site)].onsite.matrep
             HZ = mul_spin_matrices(Λ, -B[site], Z[site]) # HZ = (Λ - B⋅s) Z
             ΔZ′ = -im*√(2*Δt*kT*λ)*ξ[site] - Δt*(im+λ)*HZ
             ΔZ[site] = proj(ΔZ′, Z[site])
         end 
     else
         ints = interactions_inhomog(sys)
         for site in all_sites(sys)
-            Λ = ints[site].aniso.matrep
+            Λ = ints[site].onsite.matrep
             HZ = mul_spin_matrices(Λ, -B[site], Z[site]) # HZ = (Λ - B⋅s) Z
             ΔZ′ = -im*√(2*Δt*kT*λ)*ξ[site] - Δt*(im+λ)*HZ
             ΔZ[site] = proj(ΔZ′, Z[site])
@@ -221,14 +221,14 @@ function rhs_ll!(ΔZ, Z, B, integrator, sys)
     if is_homogeneous(sys)
         ints = interactions_homog(sys)
         for site in all_sites(sys)
-            Λ = ints[to_atom(site)].aniso.matrep
+            Λ = ints[to_atom(site)].onsite.matrep
             HZ = mul_spin_matrices(Λ, -B[site], Z[site]) # HZ = (Λ - B⋅s) Z
             ΔZ[site] = - Δt*im*HZ
         end 
     else
         ints = interactions_inhomog(sys)
         for site in all_sites(sys)
-            Λ = ints[site].aniso.matrep
+            Λ = ints[site].onsite.matrep
             HZ = mul_spin_matrices(Λ, -B[site], Z[site]) # HZ = (Λ - B⋅s) Z
             ΔZ[site] = - Δt*im*HZ
         end 

src/Reshaping.jl

@@ -31,21 +31,17 @@ function set_interactions_from_origin!(sys::System{N}) where N
     orig_ints = interactions_homog(origin)
 
     for new_i in 1:natoms(sys.crystal)
-        function map_couplings(couplings::Vector{Coupling{T}}) where T
-            new_couplings = Coupling{T}[]
-            for (; bond, J) in couplings
-                new_bond = transform_bond(sys.crystal, new_i, origin.crystal, bond)
-                isculled = bond_parity(new_bond)
-                push!(new_couplings, Coupling(isculled, new_bond, J))
-            end
-            return sort!(new_couplings, by=c->c.isculled)
-        end
-
         i = map_atom_to_crystal(sys.crystal, new_i, origin.crystal)
-        new_ints[new_i].aniso    = orig_ints[i].aniso
-        new_ints[new_i].heisen   = map_couplings(orig_ints[i].heisen)
-        new_ints[new_i].exchange = map_couplings(orig_ints[i].exchange)
-        new_ints[new_i].biquad   = map_couplings(orig_ints[i].biquad)
+        new_ints[new_i].onsite = orig_ints[i].onsite
+
+        new_pair = PairCoupling[]
+        for (; bond, bilin, biquad) in orig_ints[i].pair
+            new_bond = transform_bond(sys.crystal, new_i, origin.crystal, bond)
+            isculled = bond_parity(new_bond)
+            push!(new_pair, PairCoupling(isculled, new_bond, bilin, biquad))
+        end
+        new_pair = sort!(new_pair, by=c->c.isculled)
+        new_ints[new_i].pair = new_pair
     end
 end
 

src/SpinWaveTheory/SWTCalculations.jl

@@ -26,6 +26,16 @@ function swt_hamiltonian!(swt::SpinWaveTheory, k̃ :: Vector{Float64}, Hmat::Mat
     N  = Nf + 1
     L  = Nf * Nm
     @assert size(Hmat) == (2*L, 2*L)
+    # scaling factor (=1) if in the fundamental representation
+    M = sys.mode == :SUN ? 1 : (Ns-1)
+    no_single_ion = isempty(swt.sys.interactions_union[1].onsite.matrep)
+
+    # the "metric" of scalar biquad interaction. Here we are using the following identity:
+    # (𝐒ᵢ⋅𝐒ⱼ)² = -(𝐒ᵢ⋅𝐒ⱼ)/2 + ∑ₐ (OᵢᵃOⱼᵃ)/2, a=4,…,8
+    # where the definition of Oᵢᵃ is given in Appendix B of *Phys. Rev. B 104, 104409*
+    # Note: this is only valid for the `:dipole` mode, for `:SUN` mode, we consider 
+    # different implementations
+    biquad_metric = 1/2 * diagm([-1, -1, -1, 1/M, 1/M, 1/M, 1/M, 1/M])
 
     for k̃ᵢ in k̃
         (k̃ᵢ < 0.0 || k̃ᵢ ≥ 1.0) && throw("k̃ outside [0, 1) range")
@@ -55,60 +65,14 @@ function swt_hamiltonian!(swt::SpinWaveTheory, k̃ :: Vector{Float64}, Hmat::Mat
     # pairexchange interactions
     for matom = 1:Nm
         ints = sys.interactions_union[matom]
-        # Heisenberg exchange
-        for (; isculled, bond, J) in ints.heisen
-            isculled && break
-            sub_i, sub_j, ΔRδ = bond.i, bond.j, bond.n
-
-            tTi_μ = s̃_mat[:, :, :, sub_i]
-            tTj_ν = s̃_mat[:, :, :, sub_j]
-            phase  = exp(2im * π * dot(k̃, ΔRδ))
-            cphase = conj(phase)
-            sub_i_M1, sub_j_M1 = sub_i - 1, sub_j - 1
-
-            for m = 2:N
-                mM1 = m - 1
-                T_μ_11 = conj(tTi_μ[1, 1, :])
-                T_μ_m1 = conj(tTi_μ[m, 1, :])
-                T_μ_1m = conj(tTi_μ[1, m, :])
-                T_ν_11 = tTj_ν[1, 1, :]
-
-                for n = 2:N
-                    nM1 = n - 1
-                    δmn = δ(m, n)
-                    T_μ_mn, T_ν_mn = conj(tTi_μ[m, n, :]), tTj_ν[m, n, :]
-                    T_ν_n1 = tTj_ν[n, 1, :]
-                    T_ν_1n = tTj_ν[1, n, :]
-
-                    c1 = J * dot(T_μ_mn - δmn * T_μ_11, T_ν_11)
-                    c2 = J * dot(T_μ_11, T_ν_mn - δmn * T_ν_11)
-                    c3 = J * dot(T_μ_m1, T_ν_1n)
-                    c4 = J * dot(T_μ_1m, T_ν_n1)
-                    c5 = J * dot(T_μ_m1, T_ν_n1)
-                    c6 = J * dot(T_μ_1m, T_ν_1n)
 
-                    Hmat11[sub_i_M1*Nf+mM1, sub_i_M1*Nf+nM1] += 0.5 * c1
-                    Hmat11[sub_j_M1*Nf+mM1, sub_j_M1*Nf+nM1] += 0.5 * c2
-                    Hmat22[sub_i_M1*Nf+nM1, sub_i_M1*Nf+mM1] += 0.5 * c1
-                    Hmat22[sub_j_M1*Nf+nM1, sub_j_M1*Nf+mM1] += 0.5 * c2
-
-                    Hmat11[sub_i_M1*Nf+mM1, sub_j_M1*Nf+nM1] += 0.5 * c3 * phase
-                    Hmat22[sub_j_M1*Nf+nM1, sub_i_M1*Nf+mM1] += 0.5 * c3 * cphase
-
-                    Hmat22[sub_i_M1*Nf+mM1, sub_j_M1*Nf+nM1] += 0.5 * c4 * phase
-                    Hmat11[sub_j_M1*Nf+nM1, sub_i_M1*Nf+mM1] += 0.5 * c4 * cphase
-
-                    Hmat12[sub_i_M1*Nf+mM1, sub_j_M1*Nf+nM1] += 0.5 * c5 * phase
-                    Hmat12[sub_j_M1*Nf+nM1, sub_i_M1*Nf+mM1] += 0.5 * c5 * cphase
-                    Hmat21[sub_i_M1*Nf+mM1, sub_j_M1*Nf+nM1] += 0.5 * c6 * phase
-                    Hmat21[sub_j_M1*Nf+nM1, sub_i_M1*Nf+mM1] += 0.5 * c6 * cphase
-                end
-            end
-        end
-
-        # Quadratic exchange
-        for (; isculled, bond, J) in ints.exchange
+        for coupling in ints.pair
+            (; isculled, bond) = coupling
             isculled && break
+
+            ### Bilinear exchange
+
+            J = Mat3(coupling.bilin*I)
             sub_i, sub_j, ΔRδ = bond.i, bond.j, bond.n
 
             tTi_μ = s̃_mat[:, :, :, sub_i]
@@ -155,10 +119,10 @@ function swt_hamiltonian!(swt::SpinWaveTheory, k̃ :: Vector{Float64}, Hmat::Mat
                     Hmat21[sub_j_M1*Nf+nM1, sub_i_M1*Nf+mM1] += 0.5 * c6 * cphase
                 end
             end
-        end
 
-        for (; isculled, bond, J) in ints.biquad
-            isculled && break
+            ### Biquadratic exchange
+
+            J = coupling.biquad
             sub_i, sub_j, ΔRδ = bond.i, bond.j, bond.n
             phase  = exp(2im * π * dot(k̃, ΔRδ))
             cphase = conj(phase)

src/SpinWaveTheory/SpinWaveTheory.jl

@@ -84,22 +84,58 @@ function generate_local_sun_gens(sys :: System)
     Q_mat_N[4] = s_mat_N[1] * s_mat_N[2] + s_mat_N[2] * s_mat_N[1]
     Q_mat_N[5] = √3 * s_mat_N[3] * s_mat_N[3] - 1/√3 * S * (S+1) * I
 
-    U_mat = Matrix{ComplexF64}(undef, N, N)
-
-    s̃_mat = Array{ComplexF64, 4}(undef, N, N, 3, Nₘ)
-    T̃_mat = Array{ComplexF64, 3}(undef, N, N, Nₘ)
-    Q̃_mat = zeros(ComplexF64, N, N, 5, Nₘ)
-
-    for atom = 1:Nₘ
-        U_mat[:, 1] = sys.coherents[1, 1, 1, atom]
-        U_mat[:, 2:N] = nullspace(U_mat[:, 1]')
-        for μ = 1:3
-            s̃_mat[:, :, μ, atom] = Hermitian(U_mat' * s_mat_N[μ] * U_mat)
+    if sys.mode == :SUN
+        s_mat_N = spin_matrices(; N)
+
+        s̃_mat = Array{ComplexF64, 4}(undef, N, N, 3, Nₘ)
+        T̃_mat = Array{ComplexF64, 3}(undef, N, N, Nₘ)
+        Q̃_mat = zeros(ComplexF64, N, N, 5, Nₘ)
+
+        U_mat = Matrix{ComplexF64}(undef, N, N)
+
+        for atom = 1:Nₘ
+            U_mat[:, 1] = sys.coherents[1, 1, 1, atom]
+            U_mat[:, 2:N] = nullspace(U_mat[:, 1]')
+            @assert isapprox(U_mat * U_mat', I) "rotation matrix from (global frame to local frame) not unitary"
+            for μ = 1:3
+                s̃_mat[:, :, μ, atom] = Hermitian(U_mat' * s_mat_N[μ] * U_mat)
+            end
+            for ν = 1:5
+                Q̃_mat[:, :, ν, atom] = Hermitian(U_mat' * Q_mat_N[ν] * U_mat)
+            end
+            T̃_mat[:, :, atom] = Hermitian(U_mat' * sys.interactions_union[atom].onsite.matrep * U_mat)
         end
-        for ν = 1:5
-            Q̃_mat[:, :, ν, atom] = Hermitian(U_mat' * Q_mat_N[ν] * U_mat)
+
+    elseif sys.mode == :dipole
+        s_mat_2 = spin_matrices(N=2)
+
+        s̃_mat = Array{ComplexF64, 4}(undef, 2, 2, 3, Nₘ)
+
+        no_single_ion = isempty(sys.interactions_union[1].onsite.matrep)
+        T̃_mat = no_single_ion ? zeros(ComplexF64, 0, 0, 0) : Array{ComplexF64, 3}(undef, 2, 2, Nₘ)
+        Q̃_mat = Array{ComplexF64, 4}(undef, 2, 2, 5, Nₘ)
+
+        U_mat_2 = Matrix{ComplexF64}(undef, 2, 2)
+        U_mat_N = Matrix{ComplexF64}(undef, N, N)
+
+        for atom = 1:Nₘ
+            θ, ϕ = dipole_to_angles(sys.dipoles[1, 1, 1, atom])
+            U_mat_N[:] = exp(-1im * ϕ * s_mat_N[3]) * exp(-1im * θ * s_mat_N[2])
+            U_mat_2[:] = exp(-1im * ϕ * s_mat_2[3]) * exp(-1im * θ * s_mat_2[2])
+            @assert isapprox(U_mat_N * U_mat_N', I) "rotation matrix from (global frame to local frame) not unitary"
+            @assert isapprox(U_mat_2 * U_mat_2', I) "rotation matrix from (global frame to local frame) not unitary"
+            for μ = 1:3
+                s̃_mat[:, :, μ, atom] = Hermitian(U_mat_2' * s_mat_2[μ] * U_mat_2)
+            end
+            for ν = 1:5
+                Q̃_mat[:, :, ν, atom] = Hermitian(U_mat_N' * Q_mat_N[ν] * U_mat_N)[1:2, 1:2]
+            end
+
+            if !no_single_ion
+                T̃_mat[:, :, atom] = Hermitian(U_mat_N' * sys.interactions_union[atom].onsite.matrep * U_mat_N)[1:2, 1:2]
+            end
         end
-        T̃_mat[:, :, atom] = Hermitian(U_mat' * sys.interactions_union[atom].aniso.matrep * U_mat)
+        T̃_mat[:, :, atom] = Hermitian(U_mat' * sys.interactions_union[atom].onsite.matrep * U_mat)
     end
 
     return s̃_mat, T̃_mat, Q̃_mat
@@ -128,7 +164,7 @@ function generate_local_stevens_coefs(sys :: System)
         R[:] = [-sin(ϕ) -cos(ϕ)*cos(θ) cos(ϕ)*sin(θ);
                     cos(ϕ) -sin(ϕ)*cos(θ) sin(ϕ)*sin(θ);
                     0.0     sin(θ)        cos(θ)]
-        (; c2, c4, c6) = sys.interactions_union[atom].aniso.stvexp
+        (; c2, c4, c6) = sys.interactions_union[atom].onsite.stvexp
         SR  = Mat3(R)
         # In Cristian's note, S̃ = R S, so here we should pass SR'
         push!(R_mat, SR')

src/Sunny.jl

@@ -60,14 +60,14 @@ include("System/SpinInfo.jl")
 include("System/Types.jl")
 include("System/System.jl")
 include("System/PairExchange.jl")
-include("System/SingleIonAnisotropy.jl")
+include("System/OnsiteCoupling.jl")
 include("System/Ewald.jl")
 include("System/Interactions.jl")
 export SpinInfo, System, Site, all_sites, position_to_site,
     global_position, magnetic_moment, polarize_spin!, polarize_spins!, randomize_spins!, energy, forces,
-    spin_operators, stevens_operators, set_external_field!, set_onsite_coupling!, set_exchange!, set_biquadratic!,
+    spin_operators, stevens_operators, set_external_field!, set_onsite_coupling!, set_exchange!,
     dmvec, enable_dipole_dipole!, to_inhomogeneous, set_external_field_at!, set_vacancy_at!, set_onsite_coupling_at!,
-    symmetry_equivalent_bonds, set_exchange_at!, set_biquadratic_at!, remove_periodicity!
+    symmetry_equivalent_bonds, set_exchange_at!, remove_periodicity!
 
 include("Reshaping.jl")
 export reshape_geometry, resize_periodically, repeat_periodically, 

src/Symmetry/AllowedCouplings.jl

@@ -124,7 +124,9 @@ function transform_coupling_by_symmetry(cryst, J, symop, parity)
 end
 
 # Check whether a coupling matrix J is consistent with symmetries of a bond
-function is_coupling_valid(cryst::Crystal, b::BondPos, J::Mat3)
+function is_coupling_valid(cryst::Crystal, b::BondPos, J)
+    J isa Number && return true
+
     for sym in symmetries_between_bonds(cryst, b, b)
         J′ = transform_coupling_by_symmetry(cryst, J, sym...)
         # TODO use symprec to handle case where symmetry is inexact
@@ -135,7 +137,7 @@ function is_coupling_valid(cryst::Crystal, b::BondPos, J::Mat3)
     return true
 end
 
-function is_coupling_valid(cryst::Crystal, b::Bond, J::Mat3)
+function is_coupling_valid(cryst::Crystal, b::Bond, J)
     return is_coupling_valid(cryst, BondPos(cryst, b), J)
 end
 
@@ -261,6 +263,8 @@ function basis_for_symmetry_allowed_couplings(cryst::Crystal, b::Bond)
 end
 
 function transform_coupling_for_bonds(cryst, b, b_ref, J_ref)
+    J_ref isa Number && return J_ref
+
     syms = symmetries_between_bonds(cryst, BondPos(cryst, b), BondPos(cryst, b_ref))
     isempty(syms) && error("Bonds $b and $b_ref are not symmetry equivalent.")
     return transform_coupling_by_symmetry(cryst, J_ref, first(syms)...)
@@ -273,12 +277,11 @@ Given a reference bond `b` and coupling matrix `J` on that bond, return a list
 of symmetry-equivalent bonds (constrained to start from atom `i`), and a
 corresponding list of symmetry-transformed coupling matrices.
 """
-function all_symmetry_related_couplings_for_atom(cryst::Crystal, i::Int, b_ref::Bond, J_ref)
-    J_ref = Mat3(J_ref)
+function all_symmetry_related_couplings_for_atom(cryst::Crystal, i::Int, b_ref::Bond, J_ref::T) where T
     @assert is_coupling_valid(cryst, b_ref, J_ref)
 
     bs = Bond[]
-    Js = Mat3[]
+    Js = T[]
 
     for b in all_symmetry_related_bonds_for_atom(cryst, i, b_ref)
         push!(bs, b)