Generalized pair couplings (#176)

Initial support for generalized pair couplings.

* Changes to public interface are conservative. This commit adds behaviors but does not remove anything.
* SVD to optimally decompose a tensor-space operator as a sum of tensor products.
* Extract scalar and bilinear parts specially. No special handling of quadrupoles yet.

Project.toml

@@ -1,7 +1,7 @@
 name = "Sunny"
 uuid = "2b4a2ac8-8f8b-43e8-abf4-3cb0c45e8736"
 authors = ["The Sunny team"]
-version = "0.5.5"
+version = "0.5.6"
 
 [deps]
 CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193"

docs/src/versions.md

@@ -1,5 +1,9 @@
 # Version History
 
+## v0.5.6
+
+* General pair couplings are now supported in [`set_pair_coupling!`](@ref). 
+
 ## v0.5.5
 
 * [`reshape_supercell`](@ref) now allows reshaping to multiples of the primitive

src/Integrators.jl

@@ -121,16 +121,15 @@ function step!(sys::System{0}, integrator::ImplicitMidpoint)
 end
 
 function fast_isapprox(x, y; atol)
-    sqTol = atol^2
     acc = 0.
     for i in eachindex(x)
         diff = x[i] - y[i]
         acc += real(dot(diff,diff))
-        if acc > sqTol
+        if acc > atol^2
             return false
         end
     end
-    return true
+    return !isnan(acc)
 end
 
 
@@ -173,10 +172,10 @@ end
 
 
 # Implicit Midpoint Method applied to the nonlinear Schrödinger dynamics, as
-# proposed in Phys. Rev. B 106, 054423 (2022). Integrates dZ/dt = - i ℌ(Z) Z one
+# proposed in Phys. Rev. B 106, 054423 (2022). Integrates dZ/dt = - i H(Z) Z one
 # timestep Z → Z′ via the implicit equation
 #
-#   (Z′-Z)/Δt = - i ℌ(Z̄) Z, where Z̄ = (Z+Z′)/2
+#   (Z′-Z)/Δt = - i H(Z̄) Z, where Z̄ = (Z+Z′)/2
 #
 function step!(sys::System{N}, integrator::ImplicitMidpoint; max_iters=100) where N
     (; atol) = integrator

src/Operators/Rotation.jl

@@ -68,24 +68,51 @@ function random_orthogonal(rng, N::Int; special=false)
     return special ? O*det(O) : O
 end
 
-function unitary_for_rotation(R::Mat3; N::Int)
+# Unitary for a rotation matrix built from abstract generators.
+function unitary_for_rotation(R::Mat3, gen)
     !(R'*R ≈ I)   && error("Not an orthogonal matrix, R = $R.")
     !(det(R) ≈ 1) && error("Matrix includes a reflection, R = $R.")
-    S = spin_matrices(; N)
     n, θ = matrix_to_axis_angle(R)
-    return exp(-im*θ*(n'*S))
+    return exp(-im*θ*(n'*gen))
 end
 
+# Unitary for a rotation matrix in the N-dimensional irrep of SU(2).
+function unitary_irrep_for_rotation(R::Mat3; N::Int)
+    gen = spin_matrices(; N)
+    unitary_for_rotation(R, gen)
+end
+
+# Unitary for a rotation matrix in the (N1⊗N2⊗...)-dimensional irrep of SU(2).
+function unitary_tensor_for_rotation(R::Mat3; Ns)
+    Us = [unitary_irrep_for_rotation(R; N) for N in Ns]
+    if length(Ns) == 1
+        return Us[1]
+    elseif length(Ns) == 2
+        return kron(Us[1], Us[2])
+    else
+        error("Tensor products currently accept only two operators")
+    end
+end
+
+# TODO: Replace this with a function that takes generators.
 """
     rotate_operator(A, R)
 
 Rotates the local quantum operator `A` according to the ``3×3`` rotation matrix
 `R`.
 """
-function rotate_operator(A::Matrix, R) # TODO: Arbitrary generators of rotation
+function rotate_operator(A::HermitianC64, R)
     isempty(A) && return A
     R = convert(Mat3, R)
     N = size(A, 1)
-    U = unitary_for_rotation(R; N)
+    U = unitary_irrep_for_rotation(R; N)
+    return Hermitian(U'*A*U)
+end
+
+
+# Given an operator A in a tensor product representation labeled by (N1, N2,
+# ...), perform a physical rotation by the matrix R .
+function rotate_tensor_operator(A::Matrix, R; Ns)
+    U = unitary_tensor_for_rotation(R; Ns)
     return U'*A*U
 end

src/Operators/Spin.jl

@@ -7,17 +7,17 @@ appropriate value of `N` for a given site index.
 """
 function spin_matrices(; N::Int)
     if N == 0
-        return fill(zeros(ComplexF64,0,0), 3)
+        return fill(Hermitian(zeros(ComplexF64,0,0)), 3)
     end
 
-    S = (N-1)/2
+    S = (N-1)/2 + 0im
     j = 1:N-1
-    off = @. sqrt(2(S+1)*j - j*(j+1)) / 2
+    off = @. sqrt(2(S+1)*j - j*(j+1)) / 2 + 0im
 
-    Sx = diagm(1 => off, -1 => off)
-    Sy = diagm(1 => -im*off, -1 => +im*off)
-    Sz = diagm(S .- (0:N-1))
-    return [Sx, Sy, Sz]
+    Sx = Hermitian(diagm(1 => off, -1 => off))
+    Sy = Hermitian(diagm(1 => -im*off, -1 => +im*off))
+    Sz = Hermitian(diagm(S .- (0:N-1)))
+    return SVector(Sx, Sy, Sz)
 end
 
 

src/Operators/Stevens.jl

@@ -75,9 +75,9 @@ end
 # Construct Stevens operators as polynomials in the spin operators.
 function stevens_matrices(k::Int; N::Int)
     if k >= N
-        return fill(zeros(ComplexF64, N, N), 2k+1)
+        return fill(Hermitian(zeros(ComplexF64, N, N)), 2k+1)
     else
-        return stevens_abstract_polynomials(; J=spin_matrices(; N), k)
+        return Hermitian.(stevens_abstract_polynomials(; J=spin_matrices(; N), k))
     end
 end
 
@@ -134,7 +134,7 @@ end
 const stevens_αinv = map(inv, stevens_α)
 
 
-function matrix_to_stevens_coefficients(A::Matrix{ComplexF64})
+function matrix_to_stevens_coefficients(A::HermitianC64)
     N = size(A,1)
     @assert N == size(A,2)
 
@@ -161,7 +161,7 @@ end
 function rotate_stevens_coefficients(c, R::Mat3)
     N = length(c)
     k = Int((N-1)/2)
-    D = unitary_for_rotation(R; N)
+    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′)
@@ -187,12 +187,12 @@ print_stevens_expansion(S[1]^4 + S[2]^4 + S[3]^4)
 # Prints: (1/20)𝒪₄₀ + (1/4)𝒪₄₄ + (3/5)𝒮⁴
 ```
 """
-function print_stevens_expansion(op::Matrix{ComplexF64})
+function print_stevens_expansion(op::AbstractMatrix)
     op ≈ op' || error("Requires Hermitian operator")
     terms = String[]
 
     # Decompose op into Stevens coefficients
-    c = matrix_to_stevens_coefficients(op)
+    c = matrix_to_stevens_coefficients(hermitianpart(op))
     for k in 1:6
         for (c_km, m) in zip(reverse(c[k]), -k:k)
             abs(c_km) < 1e-12 && continue

src/Operators/TensorOperators.jl

@@ -1,16 +1,17 @@
-# Produces a matrix representation of a tensor product of operators, C = A⊗B.
-# Like built-in `kron` but with permutation. Returns C_{acbd} = A_{ab} B_{cd}.
-function kron_operator(A::AbstractMatrix, B::AbstractMatrix)
-    TS = promote_type(eltype(A), eltype(B))
-    C = zeros(TS, size(A,1), size(B,1), size(A,2), size(B,2))
-    for ci in CartesianIndices(C)
-        a, c, b, d = Tuple(ci)
-        C[ci] = A[a,b] * B[c,d]
-    end
-    return reshape(C, size(A,1)*size(B,1), size(A,2)*size(B,2))
+# Defined so that `reverse_kron(A⊗B) == B⊗A`.
+# 
+# To understand the implementation, note that:
+# - `A_ij B_kl = (A⊗B)_kilj`
+# - `(A⊗B)_ijkl = A_jl B_ik`
+function reverse_kron(C, N1, N2)
+    @assert length(C) == N2*N1*N2*N1
+    C = reshape(C, N2, N1, N2, N1)
+    C = permutedims(C, (2, 1, 4, 3))
+    return reshape(C, N1*N2, N1*N2)
 end
 
-
+# Return list of groups of indices. Within each group, the indexed values of `S`
+# are approximately equal. Assumes S is sorted.
 function degeneracy_groups(S, tol)
     acc = UnitRange{Int}[]
     isempty(S) && return acc
@@ -33,10 +34,10 @@ function svd_tensor_expansion(D::Matrix{T}, N1, N2) where T
     tol = 1e-12
 
     @assert size(D, 1) == size(D, 2) == N1*N2
-    D̃ = permutedims(reshape(D, N1, N2, N1, N2), (1,3,2,4))
+    D̃ = permutedims(reshape(D, N2, N1, N2, N1), (2,4,1,3))
     D̃ = reshape(D̃, N1*N1, N2*N2)
     (; S, U, V) = svd(D̃)
-    ret = []
+    ret = Tuple{HermitianC64, HermitianC64}[]
 
     # Rotate columns of U and V within each degenerate subspace so that all
     # columns (when reshaped) are Hermitian matrices
@@ -45,7 +46,7 @@ function svd_tensor_expansion(D::Matrix{T}, N1, N2) where T
 
         U_sub = view(U, :, range)
         n = length(range)
-        Q = zeros(n, n)
+        Q = zeros(ComplexF64, n, n)
         for k in 1:n, k′ in 1:n
             uk  = reshape(view(U_sub, :, k), N1, N1)
             uk′ = reshape(view(U_sub, :, k′), N1, N1)
@@ -69,15 +70,16 @@ function svd_tensor_expansion(D::Matrix{T}, N1, N2) where T
             # Check factors are really Hermitian
             @assert norm(u - u') < tol
             @assert norm(v - v') < tol
-            u = Hermitian(u+u')/2
-            v = Hermitian(v+v')/2
+            u = hermitianpart(u)
+            v = hermitianpart(v)
             push!(ret, (σ*u, conj(v)))
         end
     end
     return ret
 end
 
-function local_quantum_operators(A, B)
+
+function to_product_space_orig(A, B)
     (isempty(A) || isempty(B)) && error("Nonempty lists required")
 
     @assert allequal(size.(A))
@@ -88,39 +90,31 @@ function local_quantum_operators(A, B)
 
     I1 = Ref(Matrix(I, N1, N1))
     I2 = Ref(Matrix(I, N2, N2))
-    return (kron_operator.(A, I2), kron_operator.(I1, B))
+    return (kron.(A, I2), kron.(I1, B))
 end
 
+"""
+    to_product_space(A, B, Cs...)
+
+Given lists of operators acting on local Hilbert spaces individually, return the
+corresponding operators that act on the tensor product space. In typical usage,
+the inputs will represent local physical observables and the outputs will be
+used to define quantum couplings.
+"""
+function to_product_space(A, B, Cs...)
+    lists = [A, B, Cs...]
+
+    Ns = map(enumerate(lists)) do (i, list)
+        isempty(list) && error("Empty operator list in argument $i.")
+        allequal(size.(list)) || error("Unequal sized operators in argument $i.")
+        return size(first(list), 1)
+    end
 
-#=
-
-# Returns the spin operators for two sites, with Hilbert space dimensions N₁ and
-# N₂, respectively.
-function spin_pair(N1, N2)
-    S1 = spin_matrices(N=N1)
-    S2 = spin_matrices(N=N2)
-    return local_quantum_operators(S1, S2)
+    return map(enumerate(lists)) do (i, list)
+        I1 = I(prod(Ns[begin:i-1]))
+        I2 = I(prod(Ns[i+1:end]))
+        return map(list) do op
+            kron(I1, op, I2)
+        end
+    end
 end
-
-N1 = 3
-N2 = 3
-
-# Check Kronecker product of operators
-A1 = randn(N1,N1)
-A2 = randn(N1,N1)
-B1 = randn(N2,N2)
-B2 = randn(N2,N2)
-@assert kron_operator(A1, B1) * kron_operator(A2, B2) ≈ kron_operator(A1*A2, B1*B2)
-
-# Check SVD decomposition
-S1, S2 = spin_pair(N1, N2)
-B = (S1' * S2)^2                                # biquadratic interaction
-D = svd_tensor_expansion(B, N1, N2)             # a sum of 9 tensor products
-@assert sum(kron_operator(d...) for d in D) ≈ B # consistency check
-
-
-B = S1' * S2
-D = svd_tensor_expansion(B, N1, N2)             # a sum of 9 tensor products
-@assert sum(kron_operator(d...) for d in D) ≈ B
-
-=#

src/Optimization.jl

@@ -99,7 +99,7 @@ function minimize_energy!(sys::System{N}; maxiters=100, subiters=10, method=Opti
     # Functions to calculate energy and gradient for the state `αs`
     function f(αs)
         optim_set_spins!(sys, αs, ns)
-        return energy(sys) # TODO: `Sunny.energy` seems to allocate and is type-unstable (7/20/2023)
+        return energy(sys)
     end
     function g!(G, αs)
         optim_set_spins!(sys, αs, ns)

src/Reshaping.jl

@@ -42,14 +42,14 @@ function set_interactions_from_origin!(sys::System{N}) where N
         i = map_atom_to_crystal(sys.crystal, new_i, origin.crystal)
         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)
+        new_pc = PairCoupling[]
+        for pc in orig_ints[i].pair
+            new_bond = transform_bond(sys.crystal, new_i, origin.crystal, pc.bond)
             isculled = bond_parity(new_bond)
-            push!(new_pair, PairCoupling(isculled, new_bond, bilin, biquad))
+            push!(new_pc, PairCoupling(isculled, new_bond, pc.scalar, pc.bilin, pc.biquad, pc.general))
         end
-        new_pair = sort!(new_pair, by=c->c.isculled)
-        new_ints[new_i].pair = new_pair
+        new_pc = sort!(new_pc, by=c->c.isculled)
+        new_ints[new_i].pair = new_pc
     end
 end
 

src/SpinWaveTheory/SpinWaveTheory.jl

@@ -39,6 +39,14 @@ function SpinWaveTheory(sys::System{N}; energy_ϵ::Float64=1e-8, energy_tol::Flo
         error("SpinWaveTheory does not yet support long-range dipole-dipole interactions.")
     end
 
+    for int in sys.interactions_union
+        for pc in int.pair
+            if !isempty(pc.general.data)
+                error("SpinWaveTheory does not yet support general pair interactions.")
+            end
+        end
+    end
+
     # Reshape into single unit cell. A clone will always be performed, even if
     # no reshaping happens.
     cellsize_mag = cell_shape(sys) * diagm(SVector(sys.latsize))

src/Sunny.jl

@@ -29,17 +29,17 @@ import Random: randstring, RandomDevice
 const Vec3 = SVector{3, Float64}
 const Mat3 = SMatrix{3, 3, Float64, 9}
 const CVec{N} = SVector{N, ComplexF64}
+const HermitianC64 = Hermitian{ComplexF64, Matrix{ComplexF64}}
 
-
-include("Util/CartesianIndicesShifted.jl")
+@static if VERSION < v"1.10" hermitianpart(A) = Hermitian(A+A')/2 end
 
 include("Operators/Spin.jl")
 include("Operators/Rotation.jl")
 include("Operators/Stevens.jl")
 include("Operators/TensorOperators.jl")
 include("Operators/Symbolic.jl")
 include("Operators/Observables.jl")
-export spin_matrices, rotate_operator, print_stevens_expansion
+export spin_matrices, to_product_space, rotate_operator, print_stevens_expansion
 
 include("Symmetry/LatticeUtils.jl")
 include("Symmetry/SymOp.jl")
@@ -66,9 +66,9 @@ include("System/Interactions.jl")
 export SpinInfo, System, Site, eachsite, position_to_site, global_position, magnetic_moment, 
     set_coherent!, set_dipole!, polarize_spins!, randomize_spins!, energy, energy_per_site,
     spin_operators, stevens_operators, large_S_spin_operators, large_S_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!, remove_periodicity!
+    set_onsite_coupling!, set_pair_coupling!, set_exchange!, dmvec, enable_dipole_dipole!, set_external_field!,
+    to_inhomogeneous, set_external_field_at!, set_vacancy_at!, set_onsite_coupling_at!, set_exchange_at!,
+    symmetry_equivalent_bonds, remove_periodicity!
 
 include("MagneticOrdering.jl")
 export print_wrapped_intensities, suggest_magnetic_supercell, set_spiral_order!, set_spiral_order_on_sublattice!