Skip to content

Commit

Permalink
Fix scalar biquadratic interaction (#313)
Browse files Browse the repository at this point in the history 
Relevant when passed with `biquad` argument to `set_exchange!`.
  • Loading branch information
@kbarros
kbarros authored Sep 3, 2024
1 parent 216cc5d commit 53105eb
Show file tree
Hide file tree
Showing 3 changed files with 81 additions and 28 deletions.
7 changes: 6 additions & 1 deletion docs/src/versions.md
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,12 @@
# Version History

## v0.7.1
(In development)
(Sep 3, 2024)

* Fix crash in `plot_intensities!` when all data is uniform.
* Correctness fix for scalar biquadratic interactions specified with option
`biquad` to [`set_exchange!`](@ref).
* Prototype implementation of entangled units.

## v0.7.0
(Aug 30, 2024)
Expand Down
63 changes: 37 additions & 26 deletions src/System/PairExchange.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -244,9 +244,9 @@ function set_pair_coupling_aux!(sys::System, scalar::Float64, bilin::Union{Float

# Renormalize biquadratic interactions (from rcs_factors with k=2)
if sys.mode == :dipole
S1 = spin_label(sys, bond.i)
S2 = spin_label(sys, bond.j)
biquad *= (1 - 1/2S1) * (1 - 1/2S2)
s1 = spin_label(sys, bond.i)
s2 = spin_label(sys, bond.j)
biquad *= (1 - 1/2s1) * (1 - 1/2s2)
end

# Propagate all couplings by symmetry
Expand Down Expand Up @@ -315,6 +315,33 @@ function set_pair_coupling!(sys::System{N}, fn::Function, bond; extract_parts=tr
end


# Use the operator identity Qᵢ⋅g Qⱼ = (Sᵢ⋅Sⱼ)² + Sᵢ⋅Sⱼ/2 - Sᵢ²Sⱼ²/3, where Qᵢ
# are the five Stevens quadrupoles, and g is the `scalar_biquad_metric`. The
# parameter `biquad` is accepted as the coefficient to (Sᵢ⋅Sⱼ)², but is returned
# as the coefficient to Qᵢ⋅g Qⱼ. This is achieved via a shift of the bilinear
# and scalar parts. In the special case of :dipole_large_s, the limiting
# behavior is Sᵢ²Sⱼ² → sᵢ²sⱼ² (just the spin labels squared), and 𝒪(s²) → 0
# (homogeneous in quartic order of spin).
function adapt_for_biquad(scalar, bilin, biquad, sys, site1, site2)
bilin = to_float_or_mat3(bilin)
biquad = Float64(biquad)

if !iszero(biquad)
if sys.mode in (:SUN, :dipole)
s1 = spin_label(sys, to_atom(site1))
s2 = spin_label(sys, to_atom(site2))
bilin -= (bilin isa Number) ? biquad/2 : (biquad/2)*I
scalar += biquad * s1*(s1+1) * s2*(s2+1) / 3
else
@assert sys.mode == :dipole_large_s
s1 = sys.κs[to_cartesian(site1)]
s2 = sys.κs[to_cartesian(site2)]
scalar += biquad * s1^2 * s2^2 / 3
end
end
return scalar, bilin, biquad
end

"""
set_exchange!(sys::System, J, bond::Bond; biquad=0)
Expand All @@ -330,9 +357,9 @@ antisymmetric part of the exchange, where `D` is the Dzyaloshinskii-Moriya
pseudo-vector. The resulting interaction will be ``𝐃⋅(𝐒_i×𝐒_j)``.
The optional numeric parameter `biquad` multiplies a scalar biquadratic
interaction, ``(𝐒_i⋅𝐒_j)^2``, with appropriate [Interaction
Renormalization](@ref). For more general interactions, use
[`set_pair_coupling!`](@ref) instead.
interaction, ``(𝐒_i⋅𝐒_j)^2``, with [Interaction Renormalization](@ref) if
appropriate. For more general interactions, use [`set_pair_coupling!`](@ref)
instead.
# Examples
```julia
Expand All @@ -350,17 +377,9 @@ set_exchange!(sys, J, bond)
```
"""
function set_exchange!(sys::System{N}, J, bond::Bond; biquad=0.0) where N
if !iszero(biquad)
# Reinterpret `biquad (Sᵢ⋅Sⱼ)²` by shifting its bilinear part into the
# usual 3×3 exchange J. What remains in `biquad` is a coupling between
# quadratic Stevens operators O[2,q] via `scalar_biquad_metric`.
biquad = Float64(biquad)
J -= (J isa Number) ? biquad/2 : (biquad/2)*I
end

is_homogeneous(sys) || error("Use `set_exchange_at!` for an inhomogeneous system.")
bilin = to_float_or_mat3(J)
set_pair_coupling_aux!(sys, 0.0, bilin, biquad, zero(TensorDecomposition), bond)
scalar, bilin, biquad = adapt_for_biquad(0.0, J, biquad, sys, (1, 1, 1, bond.i), (1, 1, 1, bond.j))
set_pair_coupling_aux!(sys, scalar, bilin, biquad, zero(TensorDecomposition), bond)
return
end

Expand Down Expand Up @@ -448,16 +467,8 @@ See also [`set_exchange!`](@ref) for more details on specifying `J` and
instead.
"""
function set_exchange_at!(sys::System{N}, J, site1::Site, site2::Site; biquad::Number=0.0, offset=nothing) where N
if !iszero(biquad)
# Reinterpret `biquad (Sᵢ⋅Sⱼ)²` by shifting its bilinear part into the
# usual 3×3 exchange J. What remains in `biquad` is a coupling between
# quadratic Stevens operators O[2,q] via `scalar_biquad_metric`.
biquad = Float64(biquad)
J -= (J isa Number) ? biquad/2 : (biquad/2)*I
end

bilin = to_float_or_mat3(J)
set_pair_coupling_at_aux!(sys, 0.0, bilin, biquad, zero(TensorDecomposition), site1, site2, offset)
scalar, bilin, biquad = adapt_for_biquad(0.0, J, biquad, sys, site1, site2)
set_pair_coupling_at_aux!(sys, scalar, bilin, biquad, zero(TensorDecomposition), site1, site2, offset)
return
end

Expand Down
39 changes: 38 additions & 1 deletion test/test_rescaling.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@ end

# Reference scalar biquadratic without renormalization
set_exchange!(sys1, 0, Bond(1, 1, [1,0,0]); biquad=1)
@test sys1.interactions_union[1].pair[1].bilin ≈ -1/2
@test sys1.interactions_union[1].pair[1].bilin ≈ 0
@test sys1.interactions_union[1].pair[1].biquad ≈ 1

# Same thing, but with renormalization
Expand All @@ -137,3 +137,40 @@ end
@test sys2.interactions_union[1].pair[1].bilin ≈ -1/2
@test sys2.interactions_union[1].pair[1].biquad ≈ 1 * rcs
end

@testitem "Biquadratic renormalization 2" begin
# Simple dimer model
latvecs = lattice_vectors(1, 1, 1, 90, 90, 90)
cryst = Crystal(latvecs, [[0, 0, 0], [0.3, 0, 0]]; types=["A", "B"])
s1 = 3/2
s2 = 2
v1 = randn(3)
v2 = randn(3)
biquad = 1.2
bond = Bond(1, 2, [0, 0, 0])

# Biquadratic energy in dipole mode with RCS
sys = System(cryst, [1 => Moment(s=s1, g=-1), 2 => Moment(s=s1, g=-1)], :dipole)
set_dipole!(sys, v1, (1, 1, 1, 1))
set_dipole!(sys, v2, (1, 1, 1, 2))
set_exchange!(sys, 0.0, bond; biquad)
E_dipole = energy(sys)

# Biquadratic energy in SU(N) mode is same
sys = System(cryst, [1 => Moment(s=s1, g=-1), 2 => Moment(s=s1, g=-1)], :SUN)
set_dipole!(sys, v1, (1, 1, 1, 1))
set_dipole!(sys, v2, (1, 1, 1, 2))
set_exchange!(sys, 0.0, bond; biquad)
E_SUN_1 = energy(sys)
set_pair_coupling!(sys, (Si, Sj) -> biquad * (Si'*Sj)^2, bond)
E_SUN_2 = energy(sys)
@test E_dipole ≈ E_SUN_1 ≈ E_SUN_2

# Biquadratic energy in large-s limit is classical formula
sys = System(cryst, [1 => Moment(s=s1, g=-1), 2 => Moment(s=s2, g=-1)], :dipole_large_s)
set_dipole!(sys, v1, (1, 1, 1, 1))
set_dipole!(sys, v2, (1, 1, 1, 2))
set_exchange!(sys, 0.0, bond; biquad)
E_large_s = energy(sys)
@test E_large_s ≈ biquad * (sys.dipoles[1]' * sys.dipoles[2])^2
end

2 comments on commit 53105eb

@kbarros
Copy link
Member Author

@kbarros kbarros commented on 53105eb Sep 3, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/114477

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.7.1 -m "<description of version>" 53105ebb670b60887bc185dee678889b6d8a8d05
git push origin v0.7.1

Please sign in to comment.