From f9055cfd9d1047627ced8217f871f3b838aeae23 Mon Sep 17 00:00:00 2001 From: Kipton Barros Date: Sat, 15 Jul 2023 17:29:28 -0600 Subject: [PATCH] Absorb `set_biquadratic!` into `set_exchange!` Also some internal refactors that unify pair couplings. --- docs/src/versions.md | 5 +- src/Integrators.jl | 8 +- src/Reshaping.jl | 24 ++- src/SpinWaveTheory/SWTCalculations.jl | 63 +------- src/SpinWaveTheory/SpinWaveTheory.jl | 6 +- src/Sunny.jl | 4 +- src/Symmetry/AllowedCouplings.jl | 13 +- src/System/Interactions.jl | 206 ++++++++++++-------------- src/System/PairExchange.jl | 194 ++++++------------------ src/System/SingleIonAnisotropy.jl | 16 +- src/System/System.jl | 3 +- src/System/Types.jl | 16 +- test/shared.jl | 4 +- test/test_energy_consistency.jl | 4 +- test/test_lswt.jl | 3 +- 15 files changed, 206 insertions(+), 363 deletions(-) diff --git a/docs/src/versions.md b/docs/src/versions.md index 718f664e4..a8d774688 100644 --- a/docs/src/versions.md +++ b/docs/src/versions.md @@ -10,6 +10,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 @@ -60,7 +63,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 diff --git a/src/Integrators.jl b/src/Integrators.jl index ef81f4b1b..745f17245 100644 --- a/src/Integrators.jl +++ b/src/Integrators.jl @@ -174,7 +174,7 @@ 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]) @@ -182,7 +182,7 @@ function rhs_langevin!(ΔZ::Array{CVec{N}, 4}, Z::Array{CVec{N}, 4}, ξ::Array{C 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 diff --git a/src/Reshaping.jl b/src/Reshaping.jl index 537dbe5b4..fe4bbbb8b 100644 --- a/src/Reshaping.jl +++ b/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 diff --git a/src/SpinWaveTheory/SWTCalculations.jl b/src/SpinWaveTheory/SWTCalculations.jl index 1a205dfd0..9a3c7e8a7 100644 --- a/src/SpinWaveTheory/SWTCalculations.jl +++ b/src/SpinWaveTheory/SWTCalculations.jl @@ -21,7 +21,7 @@ function swt_hamiltonian!(swt::SpinWaveTheory, k̃ :: Vector{Float64}, Hmat::Mat @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].aniso.matrep) + 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 @@ -59,60 +59,13 @@ 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 + for (; isculled, bond, bilin, biquad) in ints.pair + isculled && break - 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 + J = bilin - # Quadratic exchange - for (; isculled, bond, J) in ints.exchange - isculled && break sub_i, sub_j, ΔRδ = bond.i, bond.j, bond.n tTi_μ = s̃_mat[:, :, :, sub_i] @@ -159,10 +112,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 + + J = biquad sub_i, sub_j, ΔRδ = bond.i, bond.j, bond.n tTi_μ = zeros(ComplexF64, N, N, 8) diff --git a/src/SpinWaveTheory/SpinWaveTheory.jl b/src/SpinWaveTheory/SpinWaveTheory.jl index ee12e600b..140fe98dc 100644 --- a/src/SpinWaveTheory/SpinWaveTheory.jl +++ b/src/SpinWaveTheory/SpinWaveTheory.jl @@ -94,7 +94,7 @@ function generate_local_sun_gens(sys :: System) 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].aniso.matrep * U_mat) + T̃_mat[:, :, atom] = Hermitian(U_mat' * sys.interactions_union[atom].onsite.matrep * U_mat) end elseif sys.mode == :dipole @@ -102,7 +102,7 @@ function generate_local_sun_gens(sys :: System) s̃_mat = Array{ComplexF64, 4}(undef, 2, 2, 3, Nₘ) - no_single_ion = isempty(sys.interactions_union[1].aniso.matrep) + 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ₘ) @@ -123,7 +123,7 @@ function generate_local_sun_gens(sys :: System) end if !no_single_ion - T̃_mat[:, :, atom] = Hermitian(U_mat_N' * sys.interactions_union[atom].aniso.matrep * U_mat_N)[1:2, 1:2] + T̃_mat[:, :, atom] = Hermitian(U_mat_N' * sys.interactions_union[atom].onsite.matrep * U_mat_N)[1:2, 1:2] end end end diff --git a/src/Sunny.jl b/src/Sunny.jl index 1d95e5c5b..0ce271714 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -62,9 +62,9 @@ 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, diff --git a/src/Symmetry/AllowedCouplings.jl b/src/Symmetry/AllowedCouplings.jl index 5f84e3ec5..839806394 100644 --- a/src/Symmetry/AllowedCouplings.jl +++ b/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) diff --git a/src/System/Interactions.jl b/src/System/Interactions.jl index 334b928dc..57a810644 100644 --- a/src/System/Interactions.jl +++ b/src/System/Interactions.jl @@ -1,17 +1,14 @@ function empty_interactions(na, N) return map(1:na) do _ - Interactions(empty_anisotropy(N), - Coupling{Float64}[], - Coupling{Mat3}[], - Coupling{Float64}[]) + Interactions(empty_anisotropy(N), PairCoupling[]) end end # Creates a clone of the lists of exchange interactions, which can be mutably # updated. function clone_interactions(ints::Interactions) - (; aniso, heisen, exchange, biquad) = ints - return Interactions(aniso, copy(heisen), copy(exchange), copy(biquad)) + (; onsite, pair) = ints + return Interactions(onsite, copy(pair)) end function interactions_homog(sys::System{N}) where N @@ -31,7 +28,7 @@ end Returns a copy of the system that allows for inhomogeneous interactions, which can be set using [`set_onsite_coupling_at!`](@ref), [`set_exchange_at!`](@ref), -[`set_biquadratic_at!`](@ref), and [`set_vacancy_at!`](@ref). +and [`set_vacancy_at!`](@ref). Inhomogeneous systems do not support symmetry-propagation of interactions or system reshaping. @@ -116,9 +113,9 @@ function local_energy_change(sys::System{N}, site, state::SpinState) where N (; latsize, extfield, dipoles, coherents, ewald) = sys if is_homogeneous(sys) - (; aniso, heisen, exchange, biquad) = interactions_homog(sys)[to_atom(site)] + (; onsite, pair) = interactions_homog(sys)[to_atom(site)] else - (; aniso, heisen, exchange, biquad) = interactions_inhomog(sys)[site] + (; onsite, pair) = interactions_inhomog(sys)[site] end s₀ = dipoles[site] @@ -126,48 +123,47 @@ function local_energy_change(sys::System{N}, site, state::SpinState) where N Δs = s - s₀ ΔE = 0.0 - cell = to_cell(site) - # Zeeman coupling to external field - ΔE -= sys.units.μB * extfield[site] ⋅ (sys.gs[site] * Δs) + ΔE -= sys.units.μB * dot(extfield[site], sys.gs[site], Δs) # Single-ion anisotropy, dipole or SUN mode if N == 0 - E_new, _ = energy_and_gradient_for_classical_anisotropy(s, aniso.stvexp) - E_old, _ = energy_and_gradient_for_classical_anisotropy(s₀, aniso.stvexp) + E_new, _ = energy_and_gradient_for_classical_anisotropy(s, onsite.stvexp) + E_old, _ = energy_and_gradient_for_classical_anisotropy(s₀, onsite.stvexp) ΔE += E_new - E_old else - Λ = aniso.matrep + Λ = onsite.matrep ΔE += real(dot(Z, Λ, Z) - dot(Z₀, Λ, Z₀)) end - # Heisenberg exchange - for (; bond, J) in heisen - sⱼ = dipoles[offsetc(cell, bond.n, latsize), bond.j] - ΔE += J * (Δs ⋅ sⱼ) - end - # Quadratic exchange matrix - for (; bond, J) in exchange - sⱼ = dipoles[offsetc(cell, bond.n, latsize), bond.j] - ΔE += dot(Δs, J, sⱼ) - end - - # Scalar biquadratic exchange - for (; bond, J) in biquad - cellⱼ = offsetc(cell, bond.n, latsize) + for coupling in pair + (; bond) = coupling + cellⱼ = offsetc(to_cell(site), bond.n, latsize) sⱼ = dipoles[cellⱼ, bond.j] - if sys.mode == :dipole - # Renormalization introduces a factor r and a Heisenberg term - Sᵢ = (sys.Ns[site]-1)/2 - Sⱼ = (sys.Ns[cellⱼ, bond.j]-1)/2 - S = √(Sᵢ*Sⱼ) - r = (1 - 1/S + 1/4S^2) - ΔE += J * (r*((s⋅sⱼ)^2 - (s₀⋅sⱼ)^2) - (Δs⋅sⱼ)/2) - elseif sys.mode == :large_S - ΔE += J * ((s⋅sⱼ)^2 - (s₀⋅sⱼ)^2) - elseif sys.mode == :SUN - error("Biquadratic currently unsupported in SU(N) mode.") + + if coupling.bilin isa Float64 + J = coupling.bilin::Float64 + ΔE += J * dot(Δs, sⱼ) + else + J = coupling.bilin::Mat3 + ΔE += dot(Δs, J, sⱼ) + end + + if !iszero(coupling.biquad) + J = coupling.biquad + if sys.mode == :dipole + # Renormalization introduces a factor r and a Heisenberg term + Sᵢ = (sys.Ns[site]-1)/2 + Sⱼ = (sys.Ns[cellⱼ, bond.j]-1)/2 + S = √(Sᵢ*Sⱼ) + r = (1 - 1/S + 1/4S^2) + ΔE += J * (r*((s⋅sⱼ)^2 - (s₀⋅sⱼ)^2) - (Δs⋅sⱼ)/2) + elseif sys.mode == :large_S + ΔE += J * ((s⋅sⱼ)^2 - (s₀⋅sⱼ)^2) + elseif sys.mode == :SUN + error("Biquadratic currently unsupported in SU(N) mode.") + end end end @@ -227,45 +223,42 @@ function energy_aux(sys::System{N}, ints::Interactions, i::Int, cells, foreachbo if N == 0 # Dipole mode for cell in cells s = dipoles[cell, i] - E += energy_and_gradient_for_classical_anisotropy(s, ints.aniso.stvexp)[1] + E += energy_and_gradient_for_classical_anisotropy(s, ints.onsite.stvexp)[1] end else # SU(N) mode for cell in cells - Λ = ints.aniso.matrep + Λ = ints.onsite.matrep Z = coherents[cell, i] E += real(dot(Z, Λ, Z)) end end - # Heisenberg exchange - foreachbond(ints.heisen) do J, site1, site2 + foreachbond(ints.pair) do coupling, site1, site2 sᵢ = dipoles[site1] sⱼ = dipoles[site2] - E += J * dot(sᵢ, sⱼ) - end - # Quadratic exchange matrix - foreachbond(ints.exchange) do J, site1, site2 - sᵢ = dipoles[site1] - sⱼ = dipoles[site2] - E += dot(sᵢ, J, sⱼ) - end + if coupling.bilin isa Float64 + J = coupling.bilin::Float64 + E += J * dot(sᵢ, sⱼ) + else + J = coupling.bilin::Mat3 + E += dot(sᵢ, J, sⱼ) + end - # Scalar biquadratic exchange - foreachbond(ints.biquad) do J, site1, site2 - sᵢ = dipoles[site1] - sⱼ = dipoles[site2] - if sys.mode == :dipole - # Renormalization introduces a factor r and a Heisenberg term - Sᵢ = (sys.Ns[site1]-1)/2 - Sⱼ = (sys.Ns[site2]-1)/2 - S = √(Sᵢ*Sⱼ) - r = (1 - 1/S + 1/4S^2) - E += J * (r*(sᵢ⋅sⱼ)^2 - (sᵢ⋅sⱼ)/2 + S^3 + S^2/4) - elseif sys.mode == :large_S - E += J * (sᵢ⋅sⱼ)^2 - elseif sys.mode == :SUN - error("Biquadratic currently unsupported in SU(N) mode.") + if !iszero(coupling.biquad) + J = coupling.biquad + if sys.mode == :dipole + # Renormalization introduces a factor r and a Heisenberg term + Sᵢ = (sys.Ns[site1]-1)/2 + Sⱼ = (sys.Ns[site2]-1)/2 + S = √(Sᵢ*Sⱼ) + r = (1 - 1/S + 1/4S^2) + E += J * (r*(sᵢ⋅sⱼ)^2 - (sᵢ⋅sⱼ)/2 + S^3 + S^2/4) + elseif sys.mode == :large_S + E += J * (sᵢ⋅sⱼ)^2 + elseif sys.mode == :SUN + error("Biquadratic currently unsupported in SU(N) mode.") + end end end @@ -312,58 +305,55 @@ function set_forces_aux!(B, dipoles::Array{Vec3, 4}, ints::Interactions, sys::Sy if N == 0 for cell in cells s = dipoles[cell, i] - B[cell, i] -= energy_and_gradient_for_classical_anisotropy(s, ints.aniso.stvexp)[2] + B[cell, i] -= energy_and_gradient_for_classical_anisotropy(s, ints.onsite.stvexp)[2] end end - # Heisenberg exchange - foreachbond(ints.heisen) do J, site1, site2 - sᵢ = dipoles[site1] - sⱼ = dipoles[site2] - B[site1] -= J * sⱼ - B[site2] -= J' * sᵢ - end - - # Quadratic exchange matrix - foreachbond(ints.exchange) do J, site1, site2 + foreachbond(ints.pair) do coupling, site1, site2 sᵢ = dipoles[site1] sⱼ = dipoles[site2] - B[site1] -= J * sⱼ - B[site2] -= J' * sᵢ - end - # Scalar biquadratic exchange - foreachbond(ints.biquad) do J, site1, site2 - sᵢ = dipoles[site1] - sⱼ = dipoles[site2] + if coupling.bilin isa Float64 + J = coupling.bilin::Float64 + B[site1] -= J * sⱼ + B[site2] -= J' * sᵢ + else + J = coupling.bilin::Mat3 + B[site1] -= J * sⱼ + B[site2] -= J' * sᵢ + end - if sys.mode == :dipole - Sᵢ = (sys.Ns[site1]-1)/2 - Sⱼ = (sys.Ns[site2]-1)/2 - S = √(Sᵢ*Sⱼ) - # Renormalization introduces a factor r and a Heisenberg term - r = (1 - 1/S + 1/4S^2) - B[site1] -= J * (2r*sⱼ*(sᵢ⋅sⱼ) - sⱼ/2) - B[site2] -= J * (2r*sᵢ*(sᵢ⋅sⱼ) - sᵢ/2) - elseif sys.mode == :large_S - B[site1] -= J * 2sⱼ*(sᵢ⋅sⱼ) - B[site2] -= J * 2sᵢ*(sᵢ⋅sⱼ) - elseif sys.mode == :SUN - error("Biquadratic currently unsupported in SU(N) mode.") + if !iszero(coupling.biquad) + J = coupling.biquad + if sys.mode == :dipole + Sᵢ = (sys.Ns[site1]-1)/2 + Sⱼ = (sys.Ns[site2]-1)/2 + S = √(Sᵢ*Sⱼ) + # Renormalization introduces a factor r and a Heisenberg term + r = (1 - 1/S + 1/4S^2) + B[site1] -= J * (2r*sⱼ*(sᵢ⋅sⱼ) - sⱼ/2) + B[site2] -= J * (2r*sᵢ*(sᵢ⋅sⱼ) - sᵢ/2) + elseif sys.mode == :large_S + B[site1] -= J * 2sⱼ*(sᵢ⋅sⱼ) + B[site2] -= J * 2sᵢ*(sᵢ⋅sⱼ) + elseif sys.mode == :SUN + error("Biquadratic currently unsupported in SU(N) mode.") + end end end end + # Produces a function that iterates over a list interactions for a given cell function inhomog_bond_iterator(latsize, cell) - return function foreachbond(f, ints) - for (; isculled, bond, J) in ints + return function foreachbond(f, pcs) + for pc in pcs # Early return to avoid double-counting a bond - isculled && break + pc.isculled && break # Neighboring cell may wrap the system - cell′ = offsetc(cell, bond.n, latsize) - f(J, CartesianIndex(cell, bond.i), CartesianIndex(cell′, bond.j)) + cell′ = offsetc(cell, pc.bond.n, latsize) + f(pc, CartesianIndex(cell, pc.bond.i), CartesianIndex(cell′, pc.bond.j)) end end end @@ -371,14 +361,14 @@ 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, ints) - for (; isculled, bond, J) in ints + return function foreachbond(f, pcs) + for pc in pcs # Early return to avoid double-counting a bond - isculled && break + pc.isculled && break # Iterate over all cells and periodically shifted neighbors - for (ci, cj) in zip(CartesianIndices(latsize), CartesianIndicesShifted(latsize, Tuple(bond.n))) - f(J, CartesianIndex(ci, bond.i), CartesianIndex(cj, bond.j)) + 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 diff --git a/src/System/PairExchange.jl b/src/System/PairExchange.jl index a544b9c6a..27f6eb910 100644 --- a/src/System/PairExchange.jl +++ b/src/System/PairExchange.jl @@ -14,68 +14,27 @@ function bond_parity(bond) return bond_delta > (0, 0, 0, 0) end - -""" - set_biquadratic!(sys::System, J, bond::Bond) - -Sets a scalar biquadratic interaction along `bond`, yielding a pairwise energy -``J (𝐒_i⋅𝐒_j)²``. This interaction will be propagated to equivalent bonds in -consistency with crystal symmetry. Any previous biquadratic exchange -interactions on these bonds will be overwritten. - -For systems restricted to dipoles, the biquadratic interactions will -automatically be renormalized to achieve maximum consistency with the more -variationally accurate SU(_N_) mode. This renormalization introduces also a -correction to the quadratic part of the exchange. - -See also [`set_exchange!`](@ref). -""" -function set_biquadratic!(sys::System{N}, J, bond::Bond) where N - is_homogeneous(sys) || error("Use `set_biquadratic_at!` for an inhomogeneous system.") - - # If `sys` has been reshaped, then operate first on `sys.origin`, which - # contains full symmetry information. - if !isnothing(sys.origin) - set_biquadratic!(sys.origin, J, bond) - set_interactions_from_origin!(sys) - return - end - - validate_bond(sys.crystal, bond) - - ints = interactions_homog(sys) - - - # Print a warning if an interaction already exists for bond - if any(x -> x.bond == bond, ints[bond.i].biquad) - println("Warning: Overriding biquadratic interaction for bond $bond.") - end - - for i in 1:natoms(sys.crystal) - bonds = all_symmetry_related_bonds_for_atom(sys.crystal, i, bond) - isempty(bonds) && continue - - for bond′ in bonds - # Remove any existing exchange for bond′ - matches_bond(c) = c.bond == bond′ - filter!(!matches_bond, ints[i].biquad) - - # The energy or force calculation only needs to see each bond once - isculled = bond_parity(bond′) - - # Add to list - coupling = Coupling(isculled, bond′, J) - push!(ints[i].biquad, coupling) - end - - # Sort interactions so that non-culled bonds appear first - sort!(ints[i].biquad, by=c->c.isculled) +# Convert J to Union{Float64, Mat3} +function to_float_or_mat3(J) + if J isa Number || isapprox(J, J[1] * I, atol=1e-12) + J = Float64(first(J)) + else + J = Mat3(J) end + return J end +# Internal function only +function push_coupling!(couplings, bond, bilin, biquad) + isculled = bond_parity(bond) + filter!(c -> c.bond != bond, couplings) + push!(couplings, PairCoupling(isculled, bond, bilin, biquad)) + sort!(couplings, by=c->c.isculled) + return +end """ - set_exchange!(sys::System, J, bond::Bond) + set_exchange!(sys::System, J, bond::Bond; biquad=0.) Sets a 3×3 spin-exchange matrix `J` along `bond`, yielding a pairwise interaction energy ``𝐒_i⋅J 𝐒_j``. This interaction will be propagated to @@ -84,11 +43,16 @@ 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. -Scalar `J` implies a pure Heisenberg exchange. +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 +the Dzyaloshinskii-Moriya pseudo-vector. The resulting interaction will be +``𝐃⋅(𝐒_i×𝐒_j)``. -As a convenience, `dmvec(D)` can be used to construct the antisymmetric part of -the exchange, where `D` is 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 +to dipoles, ``b`` will be automatically renormalized for maximum consistency +with the more variationally accurate SU(_N_) mode. This renormalization +introduces also a correction to the quadratic part of the exchange. # Examples ```julia @@ -104,26 +68,22 @@ set_exchange!(sys, J1, bond) J2 = 2*I + dmvec([0,0,3]) set_exchange!(sys, J2, bond) ``` - -See also [`set_biquadratic!`](@ref), [`dmvec`](@ref). """ -function set_exchange!(sys::System{N}, J, bond::Bond) where N - is_homogeneous(sys) || error("Use `set_exchange_at!` for an inhomogeneous system.") - +function set_exchange!(sys::System{N}, J, bond::Bond; biquad=0.) where N # 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) + set_exchange!(sys.origin, J, bond; biquad) set_interactions_from_origin!(sys) return end validate_bond(sys.crystal, bond) + + is_homogeneous(sys) || error("Use `set_exchange_at!` for an inhomogeneous system.") ints = interactions_homog(sys) - # Convert J to Mat3 - J = Mat3(J isa Number ? J*I : J) - is_heisenberg = isapprox(diagm([J[1,1],J[1,1],J[1,1]]), J; atol=1e-12) + J = to_float_or_mat3(J) # Verify that exchange is symmetry-consistent if !is_coupling_valid(sys.crystal, bond, J) @@ -133,36 +93,15 @@ function set_exchange!(sys::System{N}, J, bond::Bond) where N end # Print a warning if an interaction already exists for bond - if any(x -> x.bond == bond, vcat(ints[bond.i].heisen, ints[bond.i].exchange)) + if any(x -> x.bond == bond, ints[bond.i].pair) println("Warning: Overriding exchange for bond $bond.") end for i in 1:natoms(sys.crystal) bonds, Js = all_symmetry_related_couplings_for_atom(sys.crystal, i, bond, J) - isempty(bonds) && continue - for (bond′, J′) in zip(bonds, Js) - # Remove any existing exchange for bond′ - matches_bond(c) = c.bond == bond′ - filter!(!matches_bond, ints[i].heisen) - filter!(!matches_bond, ints[i].exchange) - - # The energy or force calculation only needs to see each bond once - isculled = bond_parity(bond′) - - # Add to list - if is_heisenberg - coupling = Coupling(isculled, bond′, J′[1,1]) - push!(ints[i].heisen, coupling) - else - coupling = Coupling(isculled, bond′, J′) - push!(ints[i].exchange, coupling) - end + push_coupling!(ints[i].pair, bond′, J′, biquad) end - - # Sort interactions so that non-culled bonds appear first - sort!(ints[i].heisen, by=c->c.isculled) - sort!(ints[i].exchange, by=c->c.isculled) end end @@ -211,44 +150,9 @@ function sites_to_internal_bond(sys::System{N}, site1::CartesianIndex{4}, site2: end end -# Internal function only -function push_coupling!(couplings, bond, J) - isculled = bond_parity(bond) - filter!(c -> c.bond != bond, couplings) - push!(couplings, Coupling(isculled, bond, J)) - sort!(couplings, by=c->c.isculled) - return -end """ - set_biquadratic_at!(sys::System, J, site1::Site, site2::Site; offset=nothing) - -Sets the scalar biquadratic interaction along the single bond connecting two -[`Site`](@ref)s, ignoring crystal symmetry. The system must support -inhomogeneous interactions via [`to_inhomogeneous`](@ref). - -If the system is relatively small, the direction of the bond can be ambiguous -due to possible periodic wrapping. Resolve this ambiguity by passing an explicit -`offset` vector, in multiples of unit cells. - -See also [`set_biquadratic!`](@ref). -""" -function set_biquadratic_at!(sys::System{N}, J, site1, site2; offset=nothing) where N - site1 = to_cartesian(site1) - site2 = to_cartesian(site2) - bond = sites_to_internal_bond(sys, site1, site2, offset) - - is_homogeneous(sys) && error("Use `to_inhomogeneous` first.") - ints = interactions_inhomog(sys) - - push_coupling!(ints[site1].biquad, bond, J) - push_coupling!(ints[site2].biquad, reverse(bond), J') - return -end - - -""" - set_exchange_at!(sys::System, J, site1::Site, site2::Site; offset=nothing) + set_exchange_at!(sys::System, J, site1::Site, site2::Site; biquad=0., offset=nothing) Sets the exchange interaction along the single bond connecting two [`Site`](@ref)s, ignoring crystal symmetry. The system must support @@ -260,7 +164,7 @@ 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; offset=nothing) where N +function set_exchange_at!(sys::System{N}, J, site1::Site, site2::Site; biquad=0., offset=nothing) where N site1 = to_cartesian(site1) site2 = to_cartesian(site2) bond = sites_to_internal_bond(sys, site1, site2, offset) @@ -268,17 +172,10 @@ function set_exchange_at!(sys::System{N}, J, site1::Site, site2::Site; offset=no is_homogeneous(sys) && error("Use `to_inhomogeneous` first.") ints = interactions_inhomog(sys) - # Convert J to Mat3 - J = Mat3(J isa Number ? J*I : J) - is_heisenberg = isapprox(diagm([J[1,1],J[1,1],J[1,1]]), J; atol=1e-12) - - if is_heisenberg - push_coupling!(ints[site1].heisen, bond, J[1,1]) - push_coupling!(ints[site2].heisen, reverse(bond), J[1,1]') - else - push_coupling!(ints[site1].exchange, bond, J) - push_coupling!(ints[site2].exchange, reverse(bond), J') - end + J = to_float_or_mat3(J) + push_coupling!(ints[site1].pair, bond, J, biquad) + push_coupling!(ints[site2].pair, reverse(bond), J', biquad') + return end @@ -300,16 +197,13 @@ function remove_periodicity!(sys::System{N}, dims) where N is_homogeneous(sys) && error("Use `to_inhomogeneous` first.") for site in all_sites(sys) - (; heisen, exchange, biquad) = interactions_inhomog(sys)[site] - - for ints in (heisen, exchange, biquad) - filter!(ints) do (; bond) - offset_cell = Tuple(to_cell(site)) .+ bond.n + ints = interactions_inhomog(sys)[site] + filter!(ints.pair) do (; bond) + offset_cell = Tuple(to_cell(site)) .+ bond.n - # keep bond if it is acceptable along every dimension (either - # `!dims` or if each cell component is within bounds) - all(@. !dims || 1 <= offset_cell <= sys.latsize) - end + # keep bond if it is acceptable along every dimension (either + # `!dims` or if each cell component is within bounds) + all(@. !dims || 1 <= offset_cell <= sys.latsize) end end end diff --git a/src/System/SingleIonAnisotropy.jl b/src/System/SingleIonAnisotropy.jl index 49348326b..450a40bd7 100644 --- a/src/System/SingleIonAnisotropy.jl +++ b/src/System/SingleIonAnisotropy.jl @@ -65,8 +65,8 @@ function empty_anisotropy(N) return SingleIonAnisotropy(matrep, stvexp) end -function Base.iszero(aniso::SingleIonAnisotropy) - return iszero(aniso.matrep) && iszero(aniso.stvexp.kmax) +function Base.iszero(onsite::SingleIonAnisotropy) + return iszero(onsite.matrep) && iszero(onsite.stvexp.kmax) end function rotate_operator(stevens::StevensExpansion, R) @@ -149,11 +149,11 @@ function set_onsite_coupling!(sys::System{N}, op::Matrix{ComplexF64}, i::Int) wh (1 <= i <= natoms(sys.crystal)) || error("Atom index $i is out of range.") - if !iszero(ints[i].aniso) + if !iszero(ints[i].onsite) println("Warning: Overriding anisotropy for atom $i.") end - aniso = SingleIonAnisotropy(sys, op, sys.Ns[1,1,1,i]) + onsite = SingleIonAnisotropy(sys, op, sys.Ns[1,1,1,i]) cryst = sys.crystal for j in all_symmetry_related_atoms(cryst, i) @@ -170,9 +170,9 @@ function set_onsite_coupling!(sys::System{N}, op::Matrix{ComplexF64}, i::Int) wh # In moving from site i to j, a spin S rotates to Q S. Transform the # anisotropy operator using the inverse rotation Q' so that the energy # remains invariant when applied to the transformed spins. - matrep′ = rotate_operator(aniso.matrep, Q') - stvexp′ = rotate_operator(aniso.stvexp, Q') - ints[j].aniso = SingleIonAnisotropy(matrep′, stvexp′) + matrep′ = rotate_operator(onsite.matrep, Q') + stvexp′ = rotate_operator(onsite.stvexp, Q') + ints[j].onsite = SingleIonAnisotropy(matrep′, stvexp′) end end @@ -190,7 +190,7 @@ function set_onsite_coupling_at!(sys::System{N}, op::Matrix{ComplexF64}, site::S is_homogeneous(sys) && error("Use `to_inhomogeneous` first.") ints = interactions_inhomog(sys) site = to_cartesian(site) - ints[site].aniso = SingleIonAnisotropy(sys, op, sys.Ns[site]) + ints[site].onsite = SingleIonAnisotropy(sys, op, sys.Ns[site]) end diff --git a/src/System/System.jl b/src/System/System.jl index 45840a7d2..620df9b76 100644 --- a/src/System/System.jl +++ b/src/System/System.jl @@ -264,8 +264,7 @@ end Given a [`Bond`](@ref) for the original (unreshaped) crystal, return all symmetry equivalent bonds in the [`System`](@ref). Each returned bond is represented as a pair of [`Site`](@ref)s, which may be used as input to -[`set_exchange_at!`](@ref) or [`set_biquadratic_at!`](@ref). Reverse bonds are -not included (no double counting). +[`set_exchange_at!`](@ref). Reverse bonds are not included (no double counting). # Example ```julia diff --git a/src/System/Types.jl b/src/System/Types.jl index 717628dd4..30d825bfd 100644 --- a/src/System/Types.jl +++ b/src/System/Types.jl @@ -16,17 +16,21 @@ struct SingleIonAnisotropy stvexp :: StevensExpansion # Renormalized coefficients for Stevens functions end -struct Coupling{T} +struct PairCoupling isculled :: Bool bond :: Bond - J :: T + + bilin :: Union{Float64, Mat3} # Bilinear exchange as 3×3 matrix + biquad :: Float64 # Scalar biquadratic, only valid in dipole mode + + # General pair interactions, only valid in SU(N) mode + # general :: Vector{Tuple{Hermitian{ComplexF64}, Hermitian{ComplexF64}}} + # TODO: update clone_interactions(), set_interactions_from_origin! end mutable struct Interactions - aniso :: SingleIonAnisotropy - heisen :: Vector{Coupling{Float64}} - exchange :: Vector{Coupling{Mat3}} - biquad :: Vector{Coupling{Float64}} + onsite :: SingleIonAnisotropy + pair :: Vector{PairCoupling} end const rFTPlan = FFTW.rFFTWPlan{Float64, -1, false, 5, UnitRange{Int64}} diff --git a/test/shared.jl b/test/shared.jl index 95b6f22de..90d2b8062 100644 --- a/test/shared.jl +++ b/test/shared.jl @@ -41,7 +41,9 @@ function add_quartic_interactions!(sys, mode) set_onsite_coupling!(sys, 0.2*(S[1]^4+S[2]^4+S[3]^4), 1) end if mode == :large_S - set_biquadratic!(sys, 0.2, Bond(1, 3, [0, 0, 0])) + # We must exclude :dipole because the renormalization will introduce a + # quadratic Heisenberg interaction + set_exchange!(sys, 0.0, Bond(1, 3, [0, 0, 0]); biquad=0.2) end end diff --git a/test/test_energy_consistency.jl b/test/test_energy_consistency.jl index 19aa911a2..54eed1c3a 100644 --- a/test/test_energy_consistency.jl +++ b/test/test_energy_consistency.jl @@ -27,8 +27,8 @@ @test energy(sys2) ≈ energy(sys) set_vacancy_at!(sys2, (1,1,1,1)) set_exchange_at!(sys2, 0.5, (1,1,1,1), (2,1,1,2); offset=(1, 0, 0)) - if mode != :SUN - set_biquadratic_at!(sys2, 0.7, (3,2,1,2), (3,1,1,3); offset=(0,-1,0)) + if mode != :SUN # kbtodo + set_exchange_at!(sys2, 0.0, (3,2,1,2), (3,1,1,3); biquad=0.7, offset=(0,-1,0)) end S = spin_operators(sys2, 4) diff --git a/test/test_lswt.jl b/test/test_lswt.jl index 9901225fd..195f407f9 100644 --- a/test/test_lswt.jl +++ b/test/test_lswt.jl @@ -136,8 +136,7 @@ end α = -0.4 * π J = 1.0 JL, JQ = J * cos(α), J * sin(α) / S^2 - set_exchange!(sys, JL, Bond(1, 1, [1, 0, 0])) - set_biquadratic!(sys, JQ, Bond(1, 1, [1, 0, 0])) + set_exchange!(sys, JL, Bond(1, 1, [1, 0, 0]); biquad=JQ) sys_swt = reshape_geometry(sys, [1 1 1; -1 1 0; 0 0 1]) polarize_spin!(sys_swt, (1, 0, 0), position_to_site(sys_swt, (0, 0, 0)))