From af1cf197a4c0dcb1bc7ed025bbcf76acb9620799 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Sun, 8 Oct 2023 07:37:14 -0600 Subject: [PATCH] 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 | 2 +- docs/src/versions.md | 4 + src/Integrators.jl | 9 +- src/Operators/Rotation.jl | 37 +++- src/Operators/Spin.jl | 14 +- src/Operators/Stevens.jl | 12 +- src/Operators/TensorOperators.jl | 92 +++++---- src/Optimization.jl | 2 +- src/Reshaping.jl | 12 +- src/SpinWaveTheory/SpinWaveTheory.jl | 8 + src/Sunny.jl | 12 +- src/Symmetry/AllowedAnisotropy.jl | 4 +- src/Symmetry/AllowedCouplings.jl | 48 +---- src/Symmetry/Crystal.jl | 14 ++ src/System/Interactions.jl | 267 +++++++++++++++------------ src/System/OnsiteCoupling.jl | 11 +- src/System/PairExchange.jl | 236 +++++++++++++++++------ src/System/System.jl | 7 - src/System/Types.jl | 19 +- src/Util/CartesianIndicesShifted.jl | 66 ------- test/shared.jl | 23 ++- test/test_jet.jl | 11 +- test/test_operators.jl | 7 +- test/test_symmetry.jl | 11 +- test/test_tensors.jl | 59 ++++++ 25 files changed, 578 insertions(+), 409 deletions(-) delete mode 100644 src/Util/CartesianIndicesShifted.jl create mode 100644 test/test_tensors.jl diff --git a/Project.toml b/Project.toml index e3f4ebffd..3a1d81af3 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/versions.md b/docs/src/versions.md index 0aaed6dd9..a5b090224 100644 --- a/docs/src/versions.md +++ b/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 diff --git a/src/Integrators.jl b/src/Integrators.jl index 1173d1a9f..f32e9e1f8 100644 --- a/src/Integrators.jl +++ b/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 diff --git a/src/Operators/Rotation.jl b/src/Operators/Rotation.jl index ce3e594ce..11fc6026e 100644 --- a/src/Operators/Rotation.jl +++ b/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 diff --git a/src/Operators/Spin.jl b/src/Operators/Spin.jl index 3c9df48a4..382c843d6 100644 --- a/src/Operators/Spin.jl +++ b/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 diff --git a/src/Operators/Stevens.jl b/src/Operators/Stevens.jl index aac49e693..b0af0e16a 100644 --- a/src/Operators/Stevens.jl +++ b/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 diff --git a/src/Operators/TensorOperators.jl b/src/Operators/TensorOperators.jl index f5368fc2c..eb3bad57e 100644 --- a/src/Operators/TensorOperators.jl +++ b/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 - -=# \ No newline at end of file diff --git a/src/Optimization.jl b/src/Optimization.jl index dd6abca63..1d10de262 100644 --- a/src/Optimization.jl +++ 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) diff --git a/src/Reshaping.jl b/src/Reshaping.jl index b21cb4c35..03002eae0 100644 --- a/src/Reshaping.jl +++ b/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 diff --git a/src/SpinWaveTheory/SpinWaveTheory.jl b/src/SpinWaveTheory/SpinWaveTheory.jl index 6cde296ac..8c648597e 100644 --- a/src/SpinWaveTheory/SpinWaveTheory.jl +++ b/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)) diff --git a/src/Sunny.jl b/src/Sunny.jl index bb740dd8e..395352441 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -29,9 +29,9 @@ 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") @@ -39,7 +39,7 @@ 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! diff --git a/src/Symmetry/AllowedAnisotropy.jl b/src/Symmetry/AllowedAnisotropy.jl index 14f856711..a40c6393e 100644 --- a/src/Symmetry/AllowedAnisotropy.jl +++ b/src/Symmetry/AllowedAnisotropy.jl @@ -42,7 +42,7 @@ function basis_for_symmetry_allowed_anisotropies(cryst::Crystal, i::Int; k::Int, # The Wigner D matrix, whose action on a spherical tensor corresponds to # the 3x3 rotation Q (see more below). - return unitary_for_rotation(Q; N=2k+1) + return unitary_irrep_for_rotation(Q; N=2k+1) end # A general operator in the spin-k representation can be decomposed in the @@ -63,7 +63,7 @@ function basis_for_symmetry_allowed_anisotropies(cryst::Crystal, i::Int; k::Int, # c′ satisfying c′ᵀ T′ = c T. Recall c′ᵀ T′ = (c′ᵀ D*) T = (D† c′)ᵀ T. The # constraint becomes D† c′ = c. Since D is unitary, we have c′ = D c. We # apply this transformation to each column c of C. - D = unitary_for_rotation(R; N=2k+1) + D = unitary_irrep_for_rotation(R; N=2k+1) C = D * C # Find an orthonormal basis for the columns of A, discarding linearly diff --git a/src/Symmetry/AllowedCouplings.jl b/src/Symmetry/AllowedCouplings.jl index 839806394..82d5ef603 100644 --- a/src/Symmetry/AllowedCouplings.jl +++ b/src/Symmetry/AllowedCouplings.jl @@ -118,7 +118,7 @@ function symmetry_allowed_couplings_operator(cryst::Crystal, b::BondPos) return P end -function transform_coupling_by_symmetry(cryst, J, symop, parity) +function transform_coupling_by_symmetry(cryst, J::Mat3, symop, parity) R = cryst.latvecs * symop.R * inv(cryst.latvecs) return R * (parity ? J : J') * R' end @@ -127,8 +127,8 @@ end 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...) + for (symop, parity) in symmetries_between_bonds(cryst, b, b) + J′ = transform_coupling_by_symmetry(cryst, J, symop, parity) # TODO use symprec to handle case where symmetry is inexact if !isapprox(J, J′; atol = 1e-12) return false @@ -270,45 +270,3 @@ function transform_coupling_for_bonds(cryst, b, b_ref, J_ref) return transform_coupling_by_symmetry(cryst, J_ref, first(syms)...) end -""" - all_symmetry_related_couplings_for_atom(cryst::Crystal, i::Int, b::Bond, J) - -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::T) where T - @assert is_coupling_valid(cryst, b_ref, J_ref) - - bs = Bond[] - Js = T[] - - for b in all_symmetry_related_bonds_for_atom(cryst, i, b_ref) - push!(bs, b) - push!(Js, transform_coupling_for_bonds(cryst, b, b_ref, J_ref)) - end - - return (bs, Js) -end - -""" - all_symmetry_related_couplings(cryst::Crystal, b::Bond, J) - -Given a reference bond `b` and coupling matrix `J` on that bond, return a list -of symmetry-equivalent bonds and a corresponding list of symmetry-transformed -coupling matrices. -""" -function all_symmetry_related_couplings(cryst::Crystal, b_ref::Bond, J_ref) - J_ref = Mat3(J_ref) - - bs = Bond[] - Js = Mat3[] - - for i in eachindex(cryst.positions) - (bs_i, Js_i) = all_symmetry_related_couplings_for_atom(cryst, i, b_ref, J_ref) - append!(bs, bs_i) - append!(Js, Js_i) - end - - return (bs, Js) -end diff --git a/src/Symmetry/Crystal.jl b/src/Symmetry/Crystal.jl index e9bae79cd..4533b5dc3 100644 --- a/src/Symmetry/Crystal.jl +++ b/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]] diff --git a/src/System/Interactions.jl b/src/System/Interactions.jl index bdd9a624d..7f1748d1b 100644 --- a/src/System/Interactions.jl +++ b/src/System/Interactions.jl @@ -17,8 +17,8 @@ function warn_coupling_override(str) end -# Creates a clone of the lists of exchange interactions, which can be mutably -# updated. +# Creates a copy of the Vector of PairCouplings. This is useful when cloning a +# system; mutable updates to one clone should not affect the other. function clone_interactions(ints::Interactions) (; onsite, pair) = ints return Interactions(onsite, copy(pair)) @@ -146,29 +146,38 @@ function local_energy_change(sys::System{N}, site, state::SpinState) where N E_old, _ = energy_and_gradient_for_classical_anisotropy(s₀, stvexp) ΔE += E_new - E_old else - Λ = onsite :: Matrix + Λ = onsite :: HermitianC64 ΔE += real(dot(Z, Λ, Z) - dot(Z₀, Λ, Z₀)) end # Quadratic exchange matrix - for coupling in pair - (; bond) = coupling - cellⱼ = offsetc(to_cell(site), bond.n, latsize) - sⱼ = dipoles[cellⱼ, bond.j] + for pc in pair + cellⱼ = offsetc(to_cell(site), pc.bond.n, latsize) + sⱼ = dipoles[cellⱼ, pc.bond.j] + Zⱼ = coherents[cellⱼ, pc.bond.j] # Bilinear - J = coupling.bilin + J = pc.bilin ΔE += dot(Δs, J, sⱼ) # Biquadratic - if !iszero(coupling.biquad) - J = coupling.biquad + if !iszero(pc.biquad) + J = pc.biquad if sys.mode == :dipole ΔE += J * ((s⋅sⱼ)^2 - (s₀⋅sⱼ)^2) elseif sys.mode == :SUN error("Biquadratic currently unsupported in SU(N) mode.") end end + + # General + if N > 0 + for (A, B) in pc.general.data + ΔĀ = real(dot(Z, A, Z) - dot(Z₀, A, Z₀)) + B̄ = real(dot(Zⱼ, B, Zⱼ)) + ΔE += ΔĀ * B̄ + end + end end # Long-range dipole-dipole @@ -194,30 +203,29 @@ end The total system energy. See also [`energy_per_site`](@ref). """ function energy(sys::System{N}) where N - (; crystal, latsize, dipoles, extfield, ewald) = sys - E = 0.0 # Zeeman coupling to external field for site in eachsite(sys) - E -= sys.units.μB * extfield[site] ⋅ (sys.gs[site] * dipoles[site]) + E -= sys.units.μB * sys.extfield[site] ⋅ (sys.gs[site] * sys.dipoles[site]) end # Anisotropies and exchange interactions - for i in 1:natoms(crystal) + for i in 1:natoms(sys.crystal) if is_homogeneous(sys) + # Interactions for sublattice i (same for every cell) interactions = sys.interactions_union[i] - E += energy_aux(sys, interactions, i, eachcell(sys), homog_bond_iterator(latsize)) + E += energy_aux(interactions, sys, i, eachcell(sys)) else for cell in eachcell(sys) interactions = sys.interactions_union[cell, i] - E += energy_aux(sys, interactions, i, (cell,), inhomog_bond_iterator(latsize, cell)) + E += energy_aux(interactions, sys, i, (cell,)) end end end # Long-range dipole-dipole - if !isnothing(ewald) + if !isnothing(sys.ewald) E += ewald_energy(sys) end @@ -225,42 +233,63 @@ function energy(sys::System{N}) where N end # Total energy contributed by sublattice `i`, summed over the list of `cells`. -# The function `foreachbond` enables efficient iteration over neighboring cell -# pairs (without double counting). -function energy_aux(sys::System{N}, ints::Interactions, i::Int, cells, foreachbond) where N - (; dipoles, coherents) = sys +function energy_aux(ints::Interactions, sys::System{N}, i::Int, cells) where N E = 0.0 # Single-ion anisotropy if N == 0 # Dipole mode stvexp = ints.onsite :: StevensExpansion for cell in cells - s = dipoles[cell, i] + s = sys.dipoles[cell, i] E += energy_and_gradient_for_classical_anisotropy(s, stvexp)[1] end else # SU(N) mode - Λ = ints.onsite :: Matrix + Λ = ints.onsite :: HermitianC64 for cell in cells - Z = coherents[cell, i] + Z = sys.coherents[cell, i] E += real(dot(Z, Λ, Z)) end end - foreachbond(ints.pair) do coupling, site1, site2 - sᵢ = dipoles[site1] - sⱼ = dipoles[site2] + for pc in ints.pair + (; bond, isculled) = pc + isculled && break - # Bilinear - J = coupling.bilin - E += dot(sᵢ, J, sⱼ) + for cellᵢ in cells + cellⱼ = offsetc(cellᵢ, bond.n, sys.latsize) + sᵢ = sys.dipoles[cellᵢ, bond.i] + sⱼ = sys.dipoles[cellⱼ, bond.j] - # Biquadratic - if !iszero(coupling.biquad) - J = coupling.biquad - if sys.mode == :dipole - E += J * (sᵢ⋅sⱼ)^2 - elseif sys.mode == :SUN - error("Biquadratic currently unsupported in SU(N) mode.") + # Scalar + if N > 0 + E += pc.scalar * sys.κs[cellᵢ, bond.i] * sys.κs[cellⱼ, bond.j] + else + E += pc.scalar + end + + # Bilinear + J = pc.bilin :: Union{Float64, Mat3} + E += dot(sᵢ, J, sⱼ) + + # Biquadratic + if !iszero(pc.biquad) + J = pc.biquad + if sys.mode == :dipole + E += J * (sᵢ⋅sⱼ)^2 + elseif sys.mode == :SUN + error("Biquadratic currently unsupported in SU(N) mode.") + end + end + + # General + if N > 0 + Zᵢ = sys.coherents[cellᵢ, bond.i] + Zⱼ = sys.coherents[cellⱼ, bond.j] + for (A, B) in pc.general.data + Ā = real(dot(Zᵢ, A, Zᵢ)) + B̄ = real(dot(Zⱼ, B, Zⱼ)) + E += Ā * B̄ + end end end end @@ -268,47 +297,45 @@ function energy_aux(sys::System{N}, ints::Interactions, i::Int, cells, foreachbo return E end + # Updates ∇E in-place to hold energy gradient, dE/ds, for each spin. In the case # of :SUN mode, s is interpreted as expected spin, and dE/ds only includes # contributions from Zeeman coupling, bilinear exchange, and long-range # dipole-dipole. Excluded terms include onsite coupling, and general pair # coupling (biquadratic and beyond). function set_energy_grad_dipoles!(∇E, dipoles::Array{Vec3, 4}, sys::System{N}) where N - (; crystal, latsize, extfield, ewald) = sys - fill!(∇E, zero(Vec3)) # Zeeman coupling for site in eachsite(sys) - ∇E[site] -= sys.units.μB * (sys.gs[site]' * extfield[site]) + ∇E[site] -= sys.units.μB * (sys.gs[site]' * sys.extfield[site]) end # Anisotropies and exchange interactions - for i in 1:natoms(crystal) + for i in 1:natoms(sys.crystal) if is_homogeneous(sys) - # Interaction is the same at every cell + # Interactions for sublattice i (same for every cell) interactions = sys.interactions_union[i] - set_energy_grad_dipoles_aux!(∇E, dipoles, interactions, sys, i, eachcell(sys), homog_bond_iterator(latsize)) + set_energy_grad_dipoles_aux!(∇E, dipoles, interactions, sys, i, eachcell(sys)) else for cell in eachcell(sys) - # There is a different interaction at every cell - interactions = sys.interactions_union[cell,i] - set_energy_grad_dipoles_aux!(∇E, dipoles, interactions, sys, i, (cell,), inhomog_bond_iterator(latsize, cell)) + # Interactions for sublattice i and a specific cell + interactions = sys.interactions_union[cell, i] + set_energy_grad_dipoles_aux!(∇E, dipoles, interactions, sys, i, (cell,)) end end end - if !isnothing(ewald) + if !isnothing(sys.ewald) accum_ewald_grad!(∇E, dipoles, sys) end end # Calculate the energy gradient `∇E' for the sublattice `i' at all elements of -# `cells`. The function `foreachbond` enables efficient iteration over -# neighboring cell pairs (without double counting). -function set_energy_grad_dipoles_aux!(∇E, dipoles::Array{Vec3, 4}, ints::Interactions, sys::System{N}, i::Int, cells, foreachbond) where N +# `cells`. +function set_energy_grad_dipoles_aux!(∇E, dipoles::Array{Vec3, 4}, ints::Interactions, sys::System{N}, i::Int, cells) where N # Single-ion anisotropy only contributes in dipole mode. In SU(N) mode, the - # anisotropy matrix will be incorporated directly into ℌ. + # anisotropy matrix will be incorporated directly into local H matrix. if N == 0 stvexp = ints.onsite :: StevensExpansion for cell in cells @@ -317,23 +344,29 @@ function set_energy_grad_dipoles_aux!(∇E, dipoles::Array{Vec3, 4}, ints::Inter end end - foreachbond(ints.pair) do coupling, site1, site2 - sᵢ = dipoles[site1] - sⱼ = dipoles[site2] - - # Bilinear - J = coupling.bilin - ∇E[site1] += J * sⱼ - ∇E[site2] += J' * sᵢ - - # Biquadratic - if !iszero(coupling.biquad) - J = coupling.biquad - if sys.mode == :dipole - ∇E[site1] += J * 2sⱼ*(sᵢ⋅sⱼ) - ∇E[site2] += J * 2sᵢ*(sᵢ⋅sⱼ) - elseif sys.mode == :SUN - error("Biquadratic currently unsupported in SU(N) mode.") + for pc in ints.pair + (; bond, isculled) = pc + isculled && break + + for cellᵢ in cells + cellⱼ = offsetc(cellᵢ, bond.n, sys.latsize) + sᵢ = dipoles[cellᵢ, bond.i] + sⱼ = dipoles[cellⱼ, bond.j] + + # Bilinear + J = pc.bilin + ∇E[cellᵢ, bond.i] += J * sⱼ + ∇E[cellⱼ, bond.j] += J' * sᵢ + + # Biquadratic + if !iszero(pc.biquad) + J = pc.biquad + if sys.mode == :dipole + ∇E[cellᵢ, bond.i] += J * 2sⱼ*(sᵢ⋅sⱼ) + ∇E[cellⱼ, bond.j] += J * 2sᵢ*(sᵢ⋅sⱼ) + elseif sys.mode == :SUN + error("Biquadratic currently unsupported in SU(N) mode.") + end end end end @@ -341,35 +374,66 @@ end # Updates `HZ` in-place to hold `dE/dZ̄`, which is the Schrödinger analog to the # quantity `dE/ds`. **Overwrites the first two dipole buffers in `sys`.** -function set_energy_grad_coherents!(HZ, Z, sys::System{N}) where N +function set_energy_grad_coherents!(HZ, Z::Array{CVec{N}, 4}, sys::System{N}) where N @assert N > 0 - # For efficiency, pre-calculate some of the terms associated with dE/ds, - # where s is the expected spin associated with Z. Note that dE_ds does _not_ - # include anything about the onsite coupling, the biquadratic interactions, - # or the general pair couplings, which must be handled in a more general - # way. + fill!(HZ, zero(CVec{N})) + + # Accumulate Zeeman, Ewald interactions, and spin-bilinear exchange + # interactions into dE/ds, where s is the expected spin associated with Z. + # Note that dE_ds does _not_ include the onsite coupling or biquadratic + # couplings, which must be handled differently. dE_ds, dipoles = get_dipole_buffers(sys, 2) @. dipoles = expected_spin(Z) set_energy_grad_dipoles!(dE_ds, dipoles, sys) - if is_homogeneous(sys) - ints = interactions_homog(sys) - for site in eachsite(sys) - Λ = ints[to_atom(site)].onsite :: Matrix - HZ[site] = mul_spin_matrices(Λ, dE_ds[site], Z[site]) + # Accumulate anisotropies and exchange interactions. + for i in 1:natoms(sys.crystal) + if is_homogeneous(sys) + # Interactions for sublattice i (same for every cell) + interactions = sys.interactions_union[i] + set_energy_grad_coherents_aux!(HZ, Z, dE_ds, interactions, sys, i, eachcell(sys)) + else + for cell in eachcell(sys) + # Interactions for sublattice i and a specific cell + interactions = sys.interactions_union[cell, i] + set_energy_grad_coherents_aux!(HZ, Z, dE_ds, interactions, sys, i, (cell,)) + end end - else - ints = interactions_inhomog(sys) - for site in eachsite(sys) - Λ = ints[site].onsite :: Matrix - HZ[site] = mul_spin_matrices(Λ, dE_ds[site], Z[site]) - end end - @. dE_ds = dipoles = Vec3(0,0,0) + fill!(dE_ds, zero(Vec3)) + fill!(dipoles, zero(Vec3)) end +function set_energy_grad_coherents_aux!(HZ, Z::Array{CVec{N}, 4}, dE_ds::Array{Vec3, 4}, ints::Interactions, sys::System{N}, i, cells) where N + for cell in cells + # HZ += (Λ + dE/ds S) Z + Λ = ints.onsite :: HermitianC64 + HZ[cell, i] += mul_spin_matrices(Λ, dE_ds[cell, i], Z[cell, i]) + end + + for pc in ints.pair + (; bond, isculled) = pc + isculled && break + + for (A, B) in pc.general.data + A = SMatrix{N, N}(A) + B = SMatrix{N, N}(B) + for cellᵢ in cells + cellⱼ = offsetc(cellᵢ, bond.n, sys.latsize) + Zᵢ = Z[cellᵢ, bond.i] + Zⱼ = Z[cellⱼ, bond.j] + Ā = real(dot(Zᵢ, A, Zᵢ)) + B̄ = real(dot(Zⱼ, B, Zⱼ)) + HZ[cellᵢ, bond.i] += (A * Zᵢ) * B̄ + HZ[cellⱼ, bond.j] += Ā * (B * Zⱼ) + end + end + end +end + + # Internal testing functions function energy_grad_dipoles(sys::System{N}) where N ∇E = zero(sys.dipoles) @@ -401,34 +465,3 @@ end end return :(CVec{$N}($(out...))) end - - -# Produces a function that iterates over a list interactions for a given cell -function inhomog_bond_iterator(latsize, cell) - return function foreachbond(f, pcs) - for pc in pcs - # Early return to avoid double-counting a bond - pc.isculled && break - - # Neighboring cell may wrap the system - cell′ = offsetc(cell, pc.bond.n, latsize) - f(pc, CartesianIndex(cell, pc.bond.i), CartesianIndex(cell′, pc.bond.j)) - end - end -end - -# Produces a function that iterates over a list of interactions, involving all -# pairs of cells in a homogeneous system -function homog_bond_iterator(latsize) - return function foreachbond(f, pcs) - for pc in pcs - # Early return to avoid double-counting a bond - pc.isculled && break - - # Iterate over all cells and periodically shifted neighbors - for (ci, cj) in zip(CartesianIndices(latsize), CartesianIndicesShifted(latsize, Tuple(pc.bond.n))) - f(pc, CartesianIndex(ci, pc.bond.i), CartesianIndex(cj, pc.bond.j)) - end - end - end -end diff --git a/src/System/OnsiteCoupling.jl b/src/System/OnsiteCoupling.jl index de0508a3f..fe028897f 100644 --- a/src/System/OnsiteCoupling.jl +++ b/src/System/OnsiteCoupling.jl @@ -1,13 +1,14 @@ -function onsite_coupling(sys, site, matrep::Matrix{ComplexF64}) +function onsite_coupling(sys, site, matrep::AbstractMatrix) N = sys.Ns[site] size(matrep) == (N, N) || error("Invalid matrix size.") + matrep ≈ matrep' || error("Requires Hermitian operator") if sys.mode == :SUN - return matrep + return Hermitian(matrep) elseif sys.mode == :dipole S = sys.κs[site] λ = anisotropy_renormalization(S) - c = matrix_to_stevens_coefficients(matrep) + c = matrix_to_stevens_coefficients(hermitianpart(matrep)) return StevensExpansion(λ .* c) end end @@ -41,7 +42,7 @@ function empty_anisotropy(mode, N) c = map(k -> zeros(2k+1), OffsetArray(0:6, 0:6)) return StevensExpansion(c) elseif mode == :SUN - return zeros(ComplexF64, N, N) + return Hermitian(zeros(ComplexF64, N, N)) end end @@ -91,7 +92,7 @@ function Base.getindex(this::StevensMatrices, k::Int, q::Int) k > 6 && error("Stevens operators 𝒪[k,q] currently require k <= 6.") !(-k <= q <= k) && error("Stevens operators 𝒪[k,q] require -k <= q <= k.") if k == 0 - return Matrix{ComplexF64}(I, this.N, this.N) + return HermitianC64(I, this.N, this.N) else # Stevens operators are stored in descending order: k, k-1, ... -k. return stevens_matrices(k; this.N)[k - q + 1] diff --git a/src/System/PairExchange.jl b/src/System/PairExchange.jl index cd454307c..db1df070f 100644 --- a/src/System/PairExchange.jl +++ b/src/System/PairExchange.jl @@ -29,47 +29,189 @@ function to_float_or_mat3(J) return J::Union{Float64, Mat3} end -# Internal function only -function push_coupling!(sys, couplings, bond, bilin, biquad, large_S) - # Perform renormalization derived in https://arxiv.org/abs/2304.03874 - if sys.mode == :dipole && !large_S - # Multiplicative average is the correct result for two sites with - # different S (derivation omitted). - S1 = (sys.Ns[bond.i]-1)/2 - S2 = (sys.Ns[bond.j]-1)/2 - S = sqrt(S1*S2) - r = (1 - 1/S + 1/4S^2) +# Perform renormalization derived in https://arxiv.org/abs/2304.03874 +function renormalize_biquad(sys, bond, bilin, biquad) + @assert sys.mode == :dipole + + # Multiplicative average is the correct result for two sites with + # different S (derivation omitted). + S1 = (sys.Ns[bond.i]-1)/2 + S2 = (sys.Ns[bond.j]-1)/2 + S = sqrt(S1*S2) + r = (1 - 1/S + 1/4S^2) + + # biquad (s⋅sⱼ)^2 -> biquad (r (sᵢ⋅sⱼ)^2 - (sᵢ⋅sⱼ)/2 + S^3 + S^2/4) + bilin = bilin - (biquad/2) * one(bilin) + biquad = r * biquad - # J_bq (s⋅sⱼ)^2 -> J_bq (r (sᵢ⋅sⱼ)^2 - (sᵢ⋅sⱼ)/2 + S^3 + S^2/4) - bilin = bilin - (biquad/2) * one(bilin) - biquad = r * biquad - - # Drop the constant shift, J_bq (S^3 + S^2/4). - end + # Drop the constant shift, `biquad (S^3 + S^2/4)`. + return (bilin, biquad) +end + +# Internal function only +function push_coupling!(couplings, bond::Bond, scalar::Float64, bilin::Union{Float64, Mat3}, biquad::Float64, tensordec::TensorDecomposition) # Remove previous coupling on this bond filter!(c -> c.bond != bond, couplings) # If the new coupling is exactly zero, return early - iszero(bilin) && iszero(biquad) && return + iszero(bilin) && iszero(biquad) && isempty(tensordec.data) && return # Otherwise, add the new coupling to the list isculled = bond_parity(bond) - push!(couplings, PairCoupling(isculled, bond, bilin, biquad)) + push!(couplings, PairCoupling(isculled, bond, scalar, bilin, biquad, tensordec)) + + # Sorting after each insertion will introduce quadratic scaling in length of + # `couplings`. In typical usage, the `couplings` list will be short. sort!(couplings, by=c->c.isculled) return end + +function decompose_general_coupling(op, gen1, gen2; fast) + N1 = size(gen1[1], 1) + N2 = size(gen2[1], 1) + + if fast + # Remove scalar part + scalar = real(tr(op)) + op = op - (scalar/size(op, 1))*I + + # Remove bilinear part + bilin = zeros(3, 3) + for α in 1:3, β in 1:3 + v = kron(gen1[α], gen2[β]) + J = tr(v' * op) / tr(v' * v) + @assert imag(J) < 1e-12 + bilin[α, β] = real(J) + op = op - v * bilin[α, β] + end + bilin = Mat3(bilin) + + # TODO: Remove biquadratic part + u = sum(kron(gen1[α], gen2[α]) for α in 1:3) + if norm(op) > 1e-12 && normalize(op) ≈ normalize(u^2 + u/2) + @info "Detected scalar biquadratic. Not yet optimized." + end + else + scalar = bilin = 0.0 + end + + return scalar, bilin, TensorDecomposition(gen1, gen2, svd_tensor_expansion(op, N1, N2)) +end + +function Base.zero(::Type{TensorDecomposition}) + gen = spin_matrices(; N=0) + return TensorDecomposition(gen, gen, []) +end + +function transform_coupling_by_symmetry(cryst, tensordec::TensorDecomposition, symop, parity) + (; gen1, gen2, data) = tensordec + isempty(data) && return tensordec + + if !parity + data = [(B, A) for (A, B) in data] + gen2, gen1 = (gen1, gen2) + end + R = cryst.latvecs * symop.R * inv(cryst.latvecs) + Q = R * det(R) + U1 = unitary_for_rotation(Q, gen1) + U2 = unitary_for_rotation(Q, gen2) + # 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 + +function Base.isapprox(op1::TensorDecomposition, op2::TensorDecomposition; kwargs...) + isempty(op1.data) == isempty(op2.data) && return true + op1′ = sum(kron(A, B) for (A, B) in op1.data) + op2′ = sum(kron(A, B) for (A, B) in op2.data) + return isapprox(op1′, op2′; kwargs...) +end + +function set_pair_coupling_aux!(sys::System, scalar::Float64, bilin::Union{Float64, Mat3}, biquad::Float64, tensordec::TensorDecomposition, bond::Bond) + # If `sys` has been reshaped, then operate first on `sys.origin`, which + # contains full symmetry information. + if !isnothing(sys.origin) + set_pair_coupling_aux!(sys.origin, scalar, bilin, biquad, tensordec, bond) + set_interactions_from_origin!(sys) + return + end + + # Simple checks on bond indices + validate_bond(sys.crystal, bond) + + # Verify that couplings are symmetry-consistent + if !is_coupling_valid(sys.crystal, bond, bilin) + effective = isempty(tensordec.data) ? "" : " effective" + @error """Symmetry-violating$effective exchange $bilin. + Use `print_bond(crystal, $bond)` for more information.""" + end + if !is_coupling_valid(sys.crystal, bond, tensordec) + @error """Symmetry-violating coupling. Use `print_bond(crystal, $bond)` for more information.""" + error("Interaction violates symmetry.") + end + + # Print a warning if an interaction already exists for bond + ints = interactions_homog(sys) + if any(x -> x.bond == bond, ints[bond.i].pair) + warn_coupling_override("Overriding coupling for $bond.") + end + + # Propagate all couplings by symmetry + for i in 1:natoms(sys.crystal) + for bond′ in all_symmetry_related_bonds_for_atom(sys.crystal, i, bond) + bilin′ = transform_coupling_for_bonds(sys.crystal, bond′, bond, bilin) + tensordec′ = transform_coupling_for_bonds(sys.crystal, bond′, bond, tensordec) + push_coupling!(ints[i].pair, bond′, scalar, bilin′, biquad, tensordec′) + end + end +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` is a represented as a +matrix acting in the tensor product space of the two sites, and typically +originates from [`to_product_space`](@ref). + +# Examples +```julia +# Add a bilinear and biquadratic exchange +S = spin_matrices(1/2) +Si, Sj = to_product_space(S, S) +set_pair_coupling!(sys, Si'*J1*Sj + (Si'*J2*Sj)^2, bond) +``` +""" +function set_pair_coupling!(sys::System{N}, tensordec::Matrix{ComplexF64}, bond; fast=true) where N + is_homogeneous(sys) || error("Use `set_pair_coupling_at!` for an inhomogeneous system.") + + gen1 = spin_operators(sys, bond.i) + gen2 = spin_operators(sys, bond.j) + scalar, bilin, tensordec = decompose_general_coupling(tensordec, gen1, gen2; fast) + biquad = 0.0 + + set_pair_coupling_aux!(sys, scalar, bilin, biquad, tensordec, bond) +end + """ set_exchange!(sys::System, J, bond::Bond; biquad=0, large_S=false) Sets a 3×3 spin-exchange matrix `J` along `bond`, yielding a pairwise interaction energy ``𝐒_i⋅J 𝐒_j``. This interaction will be propagated to -equivalent bonds in consistency with crystal symmetry. Any previous exchange -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. +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 parameter `J` may be scalar or matrix-valued. As a convenience, `dmvec(D)` can be used to construct the antisymmetric part of the exchange, where `D` is @@ -77,15 +219,13 @@ the Dzyaloshinskii-Moriya pseudo-vector. The resulting interaction will be ``𝐃⋅(𝐒_i×𝐒_j)``. The optional parameter `biquad` defines the strength ``b`` for scalar -biquadratic interactions of the form ``b (𝐒_i⋅𝐒_j)²`` For systems restricted +biquadratic interactions of the form ``b (𝐒_i⋅𝐒_j)²``. For systems restricted to dipoles, ``b`` will be automatically renormalized for maximum consistency with the more variationally accurate SU(_N_) mode. Set `large_S=true` to work in the large-``S`` limit and disable this renormalization. # Examples ```julia -using Sunny, LinearAlgebra - # An explicit exchange matrix J1 = [2 3 0; -3 2 0; @@ -99,40 +239,14 @@ set_exchange!(sys, J2, bond) """ function set_exchange!(sys::System{N}, J, bond::Bond; biquad=0, large_S=false) where N is_homogeneous(sys) || error("Use `set_exchange_at!` for an inhomogeneous system.") - ints = interactions_homog(sys) - - # If `sys` has been reshaped, then operate first on `sys.origin`, which - # contains full symmetry information. - if !isnothing(sys.origin) - set_exchange!(sys.origin, J, bond; biquad) - set_interactions_from_origin!(sys) - return - end - - validate_bond(sys.crystal, bond) - - J = to_float_or_mat3(J) - - # Verify that exchange is symmetry-consistent - if !is_coupling_valid(sys.crystal, bond, J) - @error """Symmetry-violating exchange: $J. - Use `print_bond(crystal, $bond)` for more information.""" - error("Interaction violates symmetry.") - end - - # Print a warning if an interaction already exists for bond - if any(x -> x.bond == bond, ints[bond.i].pair) - warn_coupling_override("Overriding coupling for $bond.") - end - - for i in 1:natoms(sys.crystal) - bonds, Js = all_symmetry_related_couplings_for_atom(sys.crystal, i, bond, J) - for (bond′, J′) in zip(bonds, Js) - push_coupling!(sys, ints[i].pair, bond′, J′, biquad, large_S) - end + bilin = to_float_or_mat3(J) + if sys.mode == :dipole && !large_S + bilin, biquad = renormalize_biquad(sys, bond, bilin, biquad) end + set_pair_coupling_aux!(sys, 0.0, bilin, Float64(biquad), zero(TensorDecomposition), bond) end + # Converts two sites to a bond with indices for possibly reshaped unit cell. For # internal use only. function sites_to_internal_bond(sys::System{N}, site1::CartesianIndex{4}, site2::CartesianIndex{4}, n_ref) where N @@ -199,9 +313,13 @@ function set_exchange_at!(sys::System{N}, J, site1::Site, site2::Site; biquad=0, is_homogeneous(sys) && error("Use `to_inhomogeneous` first.") ints = interactions_inhomog(sys) - J = to_float_or_mat3(J) - push_coupling!(sys, ints[site1].pair, bond, J, biquad, large_S) - push_coupling!(sys, ints[site2].pair, reverse(bond), J', biquad', large_S) + bilin = to_float_or_mat3(J) + if sys.mode == :dipole && !large_S + bilin, biquad = renormalize_biquad(sys, bond, bilin, biquad) + end + + push_coupling!(ints[site1].pair, bond, 0.0, bilin, Float64(biquad), zero(TensorDecomposition)) + push_coupling!(ints[site2].pair, reverse(bond), 0.0, bilin', Float64(biquad), zero(TensorDecomposition)) return end diff --git a/src/System/System.jl b/src/System/System.jl index 03fd3126a..ea2978540 100644 --- a/src/System/System.jl +++ b/src/System/System.jl @@ -445,10 +445,3 @@ function get_coherent_buffers(sys::System{N}, numrequested) where N end return view(sys.coherent_buffers, 1:numrequested) end - - -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 diff --git a/src/System/Types.jl b/src/System/Types.jl index f4d3daaee..7fefb8d0e 100644 --- a/src/System/Types.jl +++ b/src/System/Types.jl @@ -16,25 +16,34 @@ function StevensExpansion(c) return StevensExpansion(kmax, c[0], c[2], c[4], c[6]) end +struct TensorDecomposition + # Generators of rotations for (Aᵢ, Bᵢ) respectively + gen1 :: SVector{3, HermitianC64} + gen2 :: SVector{3, HermitianC64} + + # [(A₁, B₁), (A₂, B₂), ...] interpreted as ∑ₖ Aₖ ⊗ Bₖ + data :: Vector{Tuple{HermitianC64, HermitianC64}} +end + # Pair couplings are counted only once per bond struct PairCoupling isculled :: Bool # Bond directionality is used to avoid double counting bond :: Bond - # In :dipole mode, these will be renormalized couplings following - # the procedure in https://arxiv.org/abs/2304.03874 + # In :dipole mode, biquad couplings will be renormalized following the + # procedure in https://arxiv.org/abs/2304.03874 + scalar :: Float64 # Constant shift bilin :: Union{Float64, Mat3} # Bilinear exchange biquad :: Float64 # Scalar biquadratic # General pair interactions, only valid in SU(N) mode - # general :: Vector{Tuple{Hermitian{ComplexF64}, Hermitian{ComplexF64}}} - # TODO: update clone_interactions(), set_interactions_from_origin! + general :: TensorDecomposition end mutable struct Interactions # Onsite coupling is either an N×N Hermitian matrix or possibly renormalized # Stevens coefficients, depending on the mode :SUN or :dipole. - onsite :: Union{Matrix{ComplexF64}, StevensExpansion} + onsite :: Union{HermitianC64, StevensExpansion} # Pair couplings for every bond that starts at the given atom pair :: Vector{PairCoupling} end diff --git a/src/Util/CartesianIndicesShifted.jl b/src/Util/CartesianIndicesShifted.jl deleted file mode 100644 index 739e619d5..000000000 --- a/src/Util/CartesianIndicesShifted.jl +++ /dev/null @@ -1,66 +0,0 @@ -struct CartesianIndicesShifted{N} # <: AbstractArray{CartesianIndex{N},N} - limit::NTuple{N,Int} - shift::NTuple{N,Int} - - function CartesianIndicesShifted(limit, shift) - any(iszero, limit) && error("Empty iterator.") - N = length(limit) - return new{N}(limit, mod.(shift, limit)) - end -end - -CartesianIndicesShifted(a::AbstractArray, shift) = CartesianIndicesShifted(size(a), shift) - - -function Base.first(iter::CartesianIndicesShifted) - (; limit, shift) = iter - return CartesianIndex(mod1.(shift .+ 1, limit)) -end - -function Base.last(iter::CartesianIndicesShifted) - (; limit, shift) = iter - return CartesianIndex(mod1.(shift, limit)) -end - -@inline function Base.iterate(iter::CartesianIndicesShifted) - iterfirst = first(iter) - iterfirst, iterfirst -end - -@inline function Base.iterate(iter::CartesianIndicesShifted, state) - valid, I = inc_shifted(state.I, iter.limit, iter.shift) - valid || return nothing - return CartesianIndex(I), CartesianIndex(I) -end - - -# << Adapted from multidimensional.jl in Julia Base stdlib >> -# CartesianIndicesShifted continues the iteration in the next column when the -# current column is consumed. The implementation is written recursively to -# achieve this. `iterate` returns `Union{Nothing, Tuple}`, we explicitly pass a -# `valid` flag to eliminate the type instability inside the core `inc_shifted` -# logic, and this gives better runtime performance. -@inline inc_shifted(::Tuple{}, ::Tuple{}, ::Tuple{}) = false, () - -@inline function inc_shifted(state::Tuple{Int}, limit::Tuple{Int}, shift::Tuple{Int}) - i = state[1] - ip = i < limit[1] ? i+1 : 1 - if ip != shift[1]+1 - return true, (ip,) - else - # wrapped to starting index so we're finished - return false, (0,) - end -end - -@inline function inc_shifted(state::Tuple{Int,Int,Vararg{Int}}, limit::Tuple{Int,Int,Vararg{Int}}, shift::Tuple{Int,Int,Vararg{Int}}) - i = state[1] - ip = i < limit[1] ? i+1 : 1 - if ip != shift[1]+1 - return true, (ip, Base.tail(state)...) - else - valid, I = inc_shifted(Base.tail(state), Base.tail(limit), Base.tail(shift)) - return valid, (shift[1]+1, I...) - end -end - diff --git a/test/shared.jl b/test/shared.jl index 7422d393f..bbbe407e7 100644 --- a/test/shared.jl +++ b/test/shared.jl @@ -29,9 +29,28 @@ function add_exchange_interactions!(sys, _) end function add_quadratic_interactions!(sys, mode) - add_exchange_interactions!(sys, mode) + if mode == :dipole + add_exchange_interactions!(sys, mode) + else + add_exchange_interactions!(sys, mode) + + #= + # 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 + D = 0.4 # DM interaction + J_exch = [J Γ -D; + Γ J -D; + D D J+K] - # TODO: Include biquadratic in SU(N) mode + bond = Bond(1, 2, [0, 0, 0]) + S1 = spin_matrices(; N=sys.Ns[bond.i]) + S2 = spin_matrices(; N=sys.Ns[bond.j]) + Si, Sj = to_product_space(S1, S2) + set_pair_coupling!(sys, Si'*J_exch*Sj + 0.01(Si'*Sj)^2, bond) + =# + end end function add_quartic_interactions!(sys, mode) diff --git a/test/test_jet.jl b/test/test_jet.jl index 6f4f66679..e76e1ad05 100644 --- a/test/test_jet.jl +++ b/test/test_jet.jl @@ -9,20 +9,18 @@ set_exchange!(sys, -1.0, Bond(1,1,(1,0,0))) polarize_spins!(sys, (0,0,1)) - # Test type stability of LocalSampler + @test_opt energy(sys) + sampler = LocalSampler(kT=0.2, propose=propose_flip) @test_opt step!(sys, sampler) - # Test type stability with mixed proposals propose = @mix_proposals 0.5 propose_flip 0.5 propose_delta(0.2) sampler = LocalSampler(kT=0.2; propose) @test_opt step!(sys, sampler) - # Test type stability of Langevin langevin = Langevin(0.01, kT=0.2, λ=0.1) @test_opt step!(sys, langevin) - # Test type stability of ImplicitMidpoint integrator = ImplicitMidpoint(0.01) @test_opt step!(sys, integrator) end @@ -42,6 +40,11 @@ end set_exchange!(sys, -1.0, Bond(1,1,(1,0,0))) polarize_spins!(sys, (0,0,1)) + # TODO: Diagonose possible @allocated bug. BenchmarkTools.@btime + # suggests that this is actually zero-allocation. + energy(sys) + @test 16 >= @allocated energy(sys) + propose = @mix_proposals 0.5 propose_flip 0.5 propose_delta(0.2) sampler = LocalSampler(kT=0.2; propose) step!(sys, sampler) diff --git a/test/test_operators.jl b/test/test_operators.jl index d8dcdaac3..34b3d704a 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -147,9 +147,7 @@ end # Test decomposition of a random Hermitian matrix into Stevens coefficients let N = 7 # big enough to yield contributions at k=6 - A = randn(ComplexF64, N, N) - A = A + A' - + A = Hermitian(randn(ComplexF64, N, N)) c = Sunny.matrix_to_stevens_coefficients(A) acc = zeros(ComplexF64, N, N) @@ -186,8 +184,7 @@ end # Test that Stevens coefficients rotate properly let - A = randn(ComplexF64, N, N) - A = A + A' + A = Hermitian(randn(ComplexF64, N, N)) c = Sunny.matrix_to_stevens_coefficients(A) # Rotate coefficients directly diff --git a/test/test_symmetry.jl b/test/test_symmetry.jl index 346c15235..b49274e8f 100644 --- a/test/test_symmetry.jl +++ b/test/test_symmetry.jl @@ -53,11 +53,10 @@ # Calculate interaction table ref_bonds = reference_bonds(cryst, 2.) - b = ref_bonds[2] - basis = Sunny.basis_for_symmetry_allowed_couplings(cryst, b) - J = basis' * randn(length(basis)) - (bs, Js) = Sunny.all_symmetry_related_couplings_for_atom(cryst, b.i, b, J) - @test length(Js) == Sunny.coordination_number(cryst, b.i, b) + bond = ref_bonds[2] + basis = Sunny.basis_for_symmetry_allowed_couplings(cryst, bond) + bs = Sunny.all_symmetry_related_bonds_for_atom(cryst, bond.i, bond) + @test length(bs) == Sunny.coordination_number(cryst, bond.i, bond) ### Triangular lattice, primitive unit cell @@ -202,7 +201,7 @@ end # If coherents are present, perform same operation if mode == :SUN - U = Sunny.unitary_for_rotation(R; N=sys.Ns[1]) + U = Sunny.unitary_irrep_for_rotation(R; N=sys.Ns[1]) sys.coherents .= circshift(sys.coherents, (0,0,0,1)) sys.coherents .= [U*z for z in sys.coherents] end diff --git a/test/test_tensors.jl b/test/test_tensors.jl new file mode 100644 index 000000000..c29e6b721 --- /dev/null +++ b/test/test_tensors.jl @@ -0,0 +1,59 @@ + +@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) + + S = spin_matrices(; N=5) + Si, Sj = to_product_space(S, S) + 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