Tests and bug fixes.

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 interactions are partially implemented in [`set_pair_coupling!`](@ref). 
+
 ## v0.5.5
 
 * [`reshape_supercell`](@ref) now allows reshaping to multiples of the primitive

src/Operators/TensorOperators.jl

@@ -78,9 +78,13 @@ function svd_tensor_expansion(D::Matrix{T}, N1, N2) where T
     return ret
 end
 
-# Given two lists of local operators acting on sites i and j individually,
-# return the corresponding local operators that act on the tensor product space.
-function local_quantum_operators(A, B)
+"""
+    to_product_space(A, B)
+
+Given two lists of local operators acting on sites i and j individually, return
+the corresponding local operators that act on the tensor product space.
+"""
+function to_product_space(A, B)
     (isempty(A) || isempty(B)) && error("Nonempty lists required")
 
     @assert allequal(size.(A))

src/Symmetry/Crystal.jl

@@ -610,6 +610,20 @@ end
 
 #= Definitions of common crystals =#
 
+function square_crystal(; a=1.0, c=10a)
+    latvecs = lattice_vectors(a, a, c, 90, 90, 90)
+    positions = [[0, 0, 0]]
+    cryst = Crystal(latvecs, positions)
+    return cryst
+end
+
+function triangular_crystal(; a=1.0, c=10a)
+    latvecs = lattice_vectors(a, a, c, 90, 90, 120)
+    positions = [[0, 0, 0]]
+    cryst = Crystal(latvecs, positions)
+    return cryst
+end
+
 function kagome_crystal(; a=1.0, c=10a)
     latvecs = lattice_vectors(a, a, c, 90, 90, 120)
     positions = [[0, 0, 0], [0.5, 0, 0], [0, 0.5, 0]]

src/System/PairExchange.jl

@@ -89,7 +89,7 @@ function decompose_general_coupling(op, gen1, gen2; fast)
             @info "Detected scalar biquadratic. Not yet optimized."
         end
     else
-        bilin = 0
+        bilin = 0.0
     end
 
     return bilin, TensorDecomposition(gen1, gen2, svd_tensor_expansion(op, N1, N2))
@@ -112,7 +112,11 @@ function transform_coupling_by_symmetry(cryst, tensordec::TensorDecomposition, s
     Q = R * det(R)
     U1 = unitary_for_rotation(Q, gen1)
     U2 = unitary_for_rotation(Q, gen2)
-    data = [(Hermitian(U1'*A*U1), Hermitian(U2'*B*U2)) for (A, B) in data]
+    # Under the symop, coherents transform as `Z -> U Z`. Then couplings must
+    # transform as `A -> U A U'` so that the expected energy on a bond `⟨A⟩⟨B⟩`
+    # is invariant. By analogy, spin rotates as `S -> R S` and the 3×3 exchange
+    # matrix transforms as `J -> R J Rᵀ` to preserve `Sᵀ J S`.
+    data = [(Hermitian(U1*A*U1'), Hermitian(U2*B*U2')) for (A, B) in data]
     return TensorDecomposition(gen1, gen2, data)
 end
 
@@ -163,20 +167,34 @@ function set_pair_coupling_aux!(sys::System, bilin::Union{Float64, Mat3}, biquad
 end
 
 
+"""
+    spin_operators_pair(sys::System, i::Int, j::Int)
+
+Returns a pair of spin dipoles `(Si, Sj)` associated with atoms `i` and `j`,
+respectively. The components of these return values are operators that act in
+the tensor product space of the two sites, which makes them useful for defining
+interactions as input to [`set_pair_coupling!`](@ref).
+"""
+function spin_operators_pair(sys::System{N}, i::Int, j::Int) where N
+    Si = spin_matrices(N=sys.Ns[i])
+    Sj = spin_matrices(N=sys.Ns[j])
+    return to_product_space(Si, Sj)
+end
+
 """
     set_pair_coupling!(sys::System, coupling, bond)
 
 Sets an arbitrary `coupling` along `bond`. This coupling will be propagated to
 equivalent bonds in consistency with crystal symmetry. Any previous interactions
 on these bonds will be overwritten. The parameter `bond` has the form `Bond(i,
 j, offset)`, where `i` and `j` are atom indices within the unit cell, and
-`offset` is a displacement in unit cells. The `coupling` will typically be
-formed as a polynomial of operators obtained from [`spin_operators_pair`](@ref).
+`offset` is a displacement in unit cells. The `coupling` may be formed as a
+polynomial of operators obtained from [`spin_operators_pair`](@ref).
 
 # Examples
 ```julia
 # Add a bilinear and biquadratic exchange
-Si, Sj = spin_operators_pair(sys, bond)
+S = spin_operators_pair(sys, bond.i, bond.j)
 set_pair_coupling!(sys, Si'*J1*Sj + (Si'*J2*Sj)^2, bond)
 ```
 """

src/System/System.jl

@@ -445,18 +445,3 @@ function get_coherent_buffers(sys::System{N}, numrequested) where N
     end
     return view(sys.coherent_buffers, 1:numrequested)
 end
-
-
-"""
-    spin_operators_pair(sys::System, b::Bond)
-
-Returns a pair of spin dipoles `(Si, Sj)` associated with atoms `b.i` and `b.j`,
-respectively. The components of these return values are operators that act in
-the tensor product space of the two sites, which makes them useful for defining
-interactions as input to [`set_pair_coupling!`](@ref).
-"""
-function spin_operators_pair(sys::System{N}, b::Bond) where N
-    Si = spin_matrices(N=sys.Ns[b.i])
-    Sj = spin_matrices(N=sys.Ns[b.j])
-    return local_quantum_operators(Si, Sj)
-end

test/shared.jl

@@ -34,8 +34,8 @@ function add_quadratic_interactions!(sys, mode)
     else
         add_exchange_interactions!(sys, mode)
 
-        # This alternative should work, but it is too slow to enable by default.        
         #=
+        # This alternative must work, but is too slow to enable by default.
         J  = 0.5   # Anti-ferro nearest neighbor
         K  = 1.0   # Scale of Kitaev term
         Γ  = 0.2   # Off-diagonal exchange
@@ -44,9 +44,8 @@ function add_quadratic_interactions!(sys, mode)
                   Γ   J   -D;
                   D   D  J+K]
 
-        bond = Bond(1, 2, [0, 0, 0])
-        Si, Sj = spin_operators_pair(sys, bond)
-        set_pair_coupling!(sys, Si'*J_exch*Sj + 0.01(Si'*Sj)^2, bond)
+        Si, Sj = spin_operators_pair(sys, 1, 2)
+        set_pair_coupling!(sys, Si'*J_exch*Sj + 0.01(Si'*Sj)^2, Bond(1, 2, [0, 0, 0]))
         =#
     end
 end

test/test_operators.jl

@@ -277,33 +277,3 @@ end
     c = Sunny.operator_to_stevens_coefficients(p, 1.0)
     @test c ≈ c_ref
 end
-
-
-@testitem "Tensor Operators" begin
-    Ni, Nj = (3, 4)
-    Si0 = spin_matrices(N=Ni)
-    Sj0 = spin_matrices(N=Nj)
-    Si, Sj = Sunny.local_quantum_operators(Si0, Sj0)
-
-    # Basic property of Kronecker product
-    A1 = randn(ComplexF64, Ni, Ni)
-    A2 = randn(ComplexF64, Ni, Ni)
-    B1 = randn(ComplexF64, Nj, Nj)
-    B2 = randn(ComplexF64, Nj, Nj)
-    @test kron(A1, B1) * kron(A2, B2) ≈ kron(A1*A2, B1*B2)
-
-    # Check transpose
-    @test Sunny.reverse_kron(kron(A1, B1), Ni, Nj) ≈ kron(B1, A1)
-
-    # Check factorization: S₁ˣ⊗S₂ˣ + S₁ˣ⊗S₂ʸ == S₁ˣ⊗(S₂ˣ+S₂ʸ)
-    B = Si[1] * Sj[1] + Si[1] * Sj[2]
-    D = Sunny.svd_tensor_expansion(B, Ni, Nj)
-    @test length(D) == 1
-    @test sum(kron(d...) for d in D) ≈ B
-
-    # Check complicated SVD decomposition
-    B = Si' * randn(3, 3) * Sj + (Si' * randn(3, 3) * Sj)^2
-    D = Sunny.svd_tensor_expansion(B, Ni, Nj)
-    @test length(D) == 9 # a nice factorization is always possible
-    @test sum(kron(d...) for d in D) ≈ B
-end

test/test_tensors.jl

@@ -0,0 +1,58 @@
+
+@testitem "Tensors basic" begin
+    Ni, Nj = (3, 4)
+    Si0 = spin_matrices(N=Ni)
+    Sj0 = spin_matrices(N=Nj)
+    Si, Sj = Sunny.to_product_space(Si0, Sj0)
+
+    # Basic properties of Kronecker product
+    A1 = randn(ComplexF64, Ni, Ni)
+    A2 = randn(ComplexF64, Ni, Ni)
+    B1 = randn(ComplexF64, Nj, Nj)
+    B2 = randn(ComplexF64, Nj, Nj)
+    @test kron(A1, B1) * kron(A2, B2) ≈ kron(A1*A2, B1*B2)
+    @test kron(A1, A2, B1, B2) ≈ kron(kron(A1, A2), kron(B1, B2))
+
+    # Check transpose
+    @test Sunny.reverse_kron(kron(A1, B1), Ni, Nj) ≈ kron(B1, A1)
+
+    # Check factorization: S₁ˣ⊗S₂ˣ + S₁ˣ⊗S₂ʸ == S₁ˣ⊗(S₂ˣ+S₂ʸ)
+    B = Si[1] * Sj[1] + Si[1] * Sj[2]
+    D = Sunny.svd_tensor_expansion(B, Ni, Nj)
+    @test length(D) == 1
+    @test sum(kron(d...) for d in D) ≈ B
+
+    # Check complicated SVD decomposition
+    B = Si' * randn(3, 3) * Sj + (Si' * randn(3, 3) * Sj)^2
+    D = Sunny.svd_tensor_expansion(B, Ni, Nj)
+    @test length(D) == 9 # a nice factorization is always possible
+    @test sum(kron(d...) for d in D) ≈ B
+end
+
+
+@testitem "General interactions" begin
+    cryst = Sunny.diamond_crystal()
+    sys = System(cryst, (2, 2, 2), [SpinInfo(1; S=2, g=2)], :SUN)
+    randomize_spins!(sys)
+
+    J = 0.5
+    K = 1.0
+    Γ = 0.2
+    D = 0.4
+    J_exch = [J   Γ   -D;
+              Γ   J   -D;
+              D   D  J+K]
+    bond = Bond(1, 2, [0, 0, 0])
+
+    set_exchange!(sys, J_exch, bond)
+    E = energy(sys)
+    dE_dZ = Sunny.energy_grad_coherents(sys)
+
+    Si, Sj = spin_operators_pair(sys, bond.i, bond.j)
+    set_pair_coupling!(sys, Si'*J_exch*Sj, bond; fast=false)
+    E′ = energy(sys)
+    dE_dZ′ = Sunny.energy_grad_coherents(sys)
+
+    @test E ≈ E′
+    @test dE_dZ ≈ dE_dZ′        
+end