From 2620243f12dcc3ad0f803bea6072f1e430113bbd Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Thu, 12 Oct 2023 10:54:20 -0600 Subject: [PATCH] Progress --- src/Reshaping.jl | 2 +- src/SpinWaveTheory/SWTCalculations.jl | 4 +- src/System/Interactions.jl | 12 ++-- src/System/PairExchange.jl | 91 +++++++++++---------------- src/System/Types.jl | 2 +- 5 files changed, 48 insertions(+), 63 deletions(-) diff --git a/src/Reshaping.jl b/src/Reshaping.jl index bbdf5313a..ea51cbfc0 100644 --- a/src/Reshaping.jl +++ b/src/Reshaping.jl @@ -46,7 +46,7 @@ function set_interactions_from_origin!(sys::System{N}) where N 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_pc, PairCoupling(isculled, new_bond, pc.scalar, pc.bilin, pc.biquad, pc.general)) + push!(new_pc, PairCoupling(isculled, new_bond, pc.scalar, pc.bilin, pc.scalar_biquad, pc.general)) end new_pc = sort!(new_pc, by=c->c.isculled) new_ints[new_i].pair = new_pc diff --git a/src/SpinWaveTheory/SWTCalculations.jl b/src/SpinWaveTheory/SWTCalculations.jl index c50e25e80..d13241e05 100644 --- a/src/SpinWaveTheory/SWTCalculations.jl +++ b/src/SpinWaveTheory/SWTCalculations.jl @@ -146,7 +146,7 @@ function swt_hamiltonian_SUN!(swt::SpinWaveTheory, q_reshaped::Vec3, Hmat::Matri ### Biquadratic exchange - J = coupling.biquad + J = coupling.scalar_biquad Ti_11 = view(sun_basis_i, 1, 1, :) @@ -313,7 +313,7 @@ function swt_hamiltonian_dipole!(swt::SpinWaveTheory, q_reshaped::Vec3, Hmat::Ma ### Biquadratic exchange - J = coupling.biquad + J = coupling.scalar_biquad # ⟨Ω₂, Ω₁|(𝐒₁⋅𝐒₂)^2|Ω₁, Ω₂⟩ = (1-1/S+1/(4S^2)) (Ω₁⋅Ω₂)^2 - 1/2 Ω₁⋅Ω₂ + const. # The biquadratic part including the biquadratic scaling factor. Ri = R_mat[sub_i] diff --git a/src/System/Interactions.jl b/src/System/Interactions.jl index 0644ab765..ad04ba66a 100644 --- a/src/System/Interactions.jl +++ b/src/System/Interactions.jl @@ -161,8 +161,8 @@ function local_energy_change(sys::System{N}, site, state::SpinState) where N ΔE += dot(Δs, J, sⱼ) # Biquadratic - if !iszero(pc.biquad) - J = pc.biquad + if !iszero(pc.scalar_biquad) + J = pc.scalar_biquad if sys.mode in (:dipole, :dipole_large_S) ΔE += J * ((s⋅sⱼ)^2 - (s₀⋅sⱼ)^2) else @@ -272,8 +272,8 @@ function energy_aux(ints::Interactions, sys::System{N}, i::Int, cells) where N E += dot(sᵢ, J, sⱼ) # Biquadratic - if !iszero(pc.biquad) - J = pc.biquad + if !iszero(pc.scalar_biquad) + J = pc.scalar_biquad if sys.mode in (:dipole, :dipole_large_S) E += J * (sᵢ⋅sⱼ)^2 else @@ -359,8 +359,8 @@ function set_energy_grad_dipoles_aux!(∇E, dipoles::Array{Vec3, 4}, ints::Inter ∇E[cellⱼ, bond.j] += J' * sᵢ # Biquadratic - if !iszero(pc.biquad) - J = pc.biquad + if !iszero(pc.scalar_biquad) + J = pc.scalar_biquad if sys.mode in (:dipole, :dipole_large_S) ∇E[cellᵢ, bond.i] += J * 2sⱼ*(sᵢ⋅sⱼ) ∇E[cellⱼ, bond.j] += J * 2sᵢ*(sᵢ⋅sⱼ) diff --git a/src/System/PairExchange.jl b/src/System/PairExchange.jl index 689fd31c7..da1679142 100644 --- a/src/System/PairExchange.jl +++ b/src/System/PairExchange.jl @@ -29,37 +29,18 @@ function to_float_or_mat3(J) return J::Union{Float64, Mat3} end -# 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 - - # 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) +function push_coupling!(couplings, bond::Bond, scalar::Float64, bilin::Union{Float64, Mat3}, scalar_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) && isempty(tensordec.data) && return + iszero(bilin) && iszero(scalar_biquad) && isempty(tensordec.data) && return # Otherwise, add the new coupling to the list isculled = bond_parity(bond) - push!(couplings, PairCoupling(isculled, bond, scalar, bilin, biquad, tensordec)) + push!(couplings, PairCoupling(isculled, bond, scalar, bilin, scalar_biquad, tensordec)) # Sorting after each insertion will introduce quadratic scaling in length of # `couplings`. In typical usage, the `couplings` list will be short. @@ -69,9 +50,8 @@ function push_coupling!(couplings, bond::Bond, scalar::Float64, bilin::Union{Flo end -function decompose_general_coupling(op, gen1, gen2; fast) - N1 = size(gen1[1], 1) - N2 = size(gen2[1], 1) +function decompose_general_coupling(op, N1, N2; fast) + @assert size(op) == (N1*N2, N1*N2) if fast # Remove scalar part @@ -79,6 +59,8 @@ function decompose_general_coupling(op, gen1, gen2; fast) op = op - (scalar/size(op, 1))*I # Remove bilinear part + gen1 = spin_matrices_of_dim(; N=N1) + gen2 = spin_matrices_of_dim(; N=N2) bilin = zeros(3, 3) for α in 1:3, β in 1:3 v = kron(gen1[α], gen2[β]) @@ -89,16 +71,24 @@ function decompose_general_coupling(op, gen1, gen2; fast) 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." + # Remove biquadratic part + biquad = zeros(5, 5) + Oi = stevens_matrices_of_dim(2; N=N1) + Oj = stevens_matrices_of_dim(2; N=N2) + for α in 1:5, β in 1:5 + v = kron(Oi[α], Oj[β]) + J = tr(v' * op) / tr(v' * v) + @assert imag(J) < 1e-12 + biquad[α, β] = real(J) + op = op - v * bilin[α, β] end + biquad = SMatrix{5, 5}(biquad) else scalar = bilin = 0.0 + biquad = nothing end - return scalar, bilin, TensorDecomposition(gen1, gen2, svd_tensor_expansion(op, N1, N2)) + return scalar, bilin, biquad, TensorDecomposition(gen1, gen2, svd_tensor_expansion(op, N1, N2)) end function Base.zero(::Type{TensorDecomposition}) @@ -133,11 +123,11 @@ function Base.isapprox(op1::TensorDecomposition, op2::TensorDecomposition; kwarg 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) +function set_pair_coupling_aux!(sys::System, scalar::Float64, bilin::Union{Float64, Mat3}, 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_pair_coupling_aux!(sys.origin, scalar, bilin, tensordec, bond) set_interactions_from_origin!(sys) return end @@ -167,7 +157,8 @@ function set_pair_coupling_aux!(sys::System, scalar::Float64, bilin::Union{Float 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′) + scalar_biquad = 0.0 + push_coupling!(ints[i].pair, bond′, scalar, bilin′, scalar_biquad, tensordec′) end end end @@ -199,12 +190,12 @@ function set_pair_coupling!(sys::System{N}, op::AbstractMatrix, bond; fast=true) error("System with mode `:dipole_large_S` requires a symbolic operator.") end - gen1 = spin_matrices(spin_irrep_label(sys, bond.i)) - gen2 = spin_matrices(spin_irrep_label(sys, bond.j)) - scalar, bilin, tensordec = decompose_general_coupling(op, gen1, gen2; fast) - biquad = 0.0 + N1 = Int(2spin_irrep_label(sys, bond.i)+1) + N2 = Int(2spin_irrep_label(sys, bond.j)+1) + scalar, bilin, biquad, tensordec = decompose_general_coupling(op, N1, N2; fast) + scalar_biquad = 0.0 - set_pair_coupling_aux!(sys, scalar, bilin, biquad, tensordec, bond) + set_pair_coupling_aux!(sys, scalar, bilin, scalar_biquad, tensordec, bond) end function set_pair_coupling!(sys::System{N}, fn::Function, bond) where N @@ -250,15 +241,12 @@ J2 = 2*I + dmvec([0,0,3]) set_exchange!(sys, J2, bond) ``` """ -function set_exchange!(sys::System{N}, J, bond::Bond; biquad=0, large_S=false) where N - large_S && @assert sys.mode == :dipole_large_S +function set_exchange!(sys::System{N}, J, bond::Bond; scalar_biquad=nothing) where N + !isnothing(scalar_biquad) && error("Use `set_pair_coupling!` to assign a biquadratic interaction.") is_homogeneous(sys) || error("Use `set_exchange_at!` for an inhomogeneous system.") bilin = to_float_or_mat3(J) - if sys.mode == :dipole - bilin, biquad = renormalize_biquad(sys, bond, bilin, biquad) - end - set_pair_coupling_aux!(sys, 0.0, bilin, Float64(biquad), zero(TensorDecomposition), bond) + set_pair_coupling_aux!(sys, 0.0, bilin, zero(TensorDecomposition), bond) end @@ -308,7 +296,7 @@ end """ - set_exchange_at!(sys::System, J, site1::Site, site2::Site; biquad=0, large_S=false, offset=nothing) + set_exchange_at!(sys::System, J, site1::Site, site2::Site; offset=nothing) Sets the exchange interaction along the single bond connecting two [`Site`](@ref)s, ignoring crystal symmetry. The system must support @@ -320,7 +308,9 @@ due to possible periodic wrapping. Resolve this ambiguity by passing an explicit See also [`set_exchange!`](@ref). """ -function set_exchange_at!(sys::System{N}, J, site1::Site, site2::Site; biquad=0, large_S=false, offset=nothing) where N +function set_exchange_at!(sys::System{N}, J, site1::Site, site2::Site; biquad=nothing, offset=nothing) where N + !isnothing(biquad) && error("Use `set_pair_coupling_at!` to assign a biquadratic interaction.") + site1 = to_cartesian(site1) site2 = to_cartesian(site2) bond = sites_to_internal_bond(sys, site1, site2, offset) @@ -328,13 +318,8 @@ 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) - 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)) + push_coupling!(ints[site1].pair, bond, 0.0, bilin, zero(TensorDecomposition)) + push_coupling!(ints[site2].pair, reverse(bond), 0.0, bilin', zero(TensorDecomposition)) return end diff --git a/src/System/Types.jl b/src/System/Types.jl index 5892e0544..9eda56203 100644 --- a/src/System/Types.jl +++ b/src/System/Types.jl @@ -34,7 +34,7 @@ struct PairCoupling # procedure in https://arxiv.org/abs/2304.03874 scalar :: Float64 # Constant shift bilin :: Union{Float64, Mat3} # Bilinear exchange - biquad :: Float64 # Scalar biquadratic + scalar_biquad :: Float64 # Scalar biquadratic # General pair interactions, only valid in SU(N) mode general :: TensorDecomposition