PairCoupling refactors #94

kbarros · 2023-07-16T13:11:44Z

Internal cleanups to unify treatment of pair couplings. In particular, each atom now stores only a single list of PairCouplings, which may contain bilinear (scalar or matrix) and biquadratic (scalar-only) exchange.

The function set_biquadratic! has been removed. Instead, call set_exchange! with the keyword argument biquad.

The changes to the internal datastructures will require updates to LSWT, which should also be incorporated into the branches of @Hao-Phys and @hlane33.

Lazersmoke · 2023-07-17T13:15:52Z

src/SpinWaveTheory/SWTCalculations.jl

-        # 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 coupling in ints.pair
+            (; isculled, bond) = coupling
+            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
+            ### Bilinear exchange
+            if coupling.bilin isa Number
+                J = Mat3(coupling.bilin*I)
+            else
+                J = coupling.bilin


It seems like a problem to me that these lines are just removed now? Is there supposed to be further changes re-introducing the behavior?

The behavior should be the same (if J is stored as a scalar, it has the meaning of a scalar times the 3x3 identity matrix). The priority here is code de-duplication. Performance optimization will be a separate pass.

Lazersmoke · 2023-07-17T13:18:51Z

src/System/Interactions.jl

-        sⱼ = dipoles[site2]
-        E += dot(sᵢ, J, sⱼ)
-    end
+        if coupling.bilin isa Float64


Dynamic dispatch? Benchmark. Also if you run this in @profview and see red in the top row on this line, that's an easy way to check if this is a meaningful problem

Actually scratch that, this branch can be removed entirely! dot(v,::Float64,v) is defined, let dot branch internally! :) @kbarros

There is an open Julia issue that Union{Float64, Mat3} cannot be stored on the stack, and would get heap allocated instead.
JuliaLang/julia#29289

That is the reason to never load coupling.bilin into a local variable, but instead put it behind the isa conditionals. Those statements are treated specially by the optimizer at inference time.

After some experimentation, the changes that Sam suggested seem to work well without heap allocation. Specifically, in this PR, it is possible to do:

J = coupling.bilin dot(v, J, v)

and the compiler can avoid heap allocation.

Lazersmoke · 2023-07-17T13:23:27Z

src/System/Interactions.jl

+        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


Ditto, these code are the same thing and is type stable, don't branch explicitly

OH, I see you're restricting the possible options it has to deal with. I think it's probably still not needed...

Another possibility: This is what UniformScaling was designed for, right? there's also dot(v,::UniformScaling,v)

This code was the result of some experiments and reading the docs. Maybe easiest to discuss on slack.

Following that discussion, I have updated the PR to use multiplication of J :: Union{Float64,Mat3}, which Julia seems to optimize well. Also, if there were type instability, the JET.jl unit test seem to be good about finding it.

Lazersmoke · 2023-07-17T13:29:24Z

src/System/PairExchange.jl

-        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)


Why is this not a floating point equality check? Is there some situation where we want to round a user's matrix?

You mean why not J isa Float64?

This is because users frequently pass us integers, and expect them to be treated as floats. For example, set_exchange!(sys, 2, bond) should work just like set_exchange!(sys, 2.0, bond).

Also, the code above can figure out that set_exchange!(sys, Mat3(2I), bond) should function like set_exchange!(sys, 2.0, bond). I.e. it will downcast from Mat3 to Float64 if it recognizes that the matrix is diagonal.

I will think about how to improve the comments.

No, I mean why isapprox instead of ==. I see now that this behavior was already here before the PR, but why?

Let's avoid floating point equality checks wherever possible.

As a concrete example, suppose we are generating J by a rotation: J = R * J0 * R', where J0 happens to be exactly diagonal. Mathematically, we should have R * R' = I, but numerically this might only be approximate. The approximate equality check will correctly infer that J is intended to be diagonal.

I agree that whether it is intended to be diagonal is what's important here. In that case, it seems to me that the user should explicitly round off the off-diagonal entries (equivalent of Chop in Mathematica, to signal the intent to be diagonal), and that Sunny rounding it for them would be surprising behavior. But not tooooo surprising and also somewhat helpful, so this convention is probably OK.

Also, from the user's perspective, the stakes are:

change in numerical result at the level of 1e-12, presumably not a big deal in most contexts.

possibly significant speedup

Its hard to imagine a circumstance where 1e-12 deviation would be significant, especially given that double precision floating point round off will already be present at this order. One interesting case is finite difference approximation to derivatives. (But in that context, the small perturbations should likely not be smaller than 1e-6 anyway).

The general principle I follow is that one should almost never check exact FP equality.

After thinking about it more, maybe we should log a warning if this situation is detected.

Hmm, I think I agree with you now 😆. The situation is rare, performance difference will likely be small, and using a full Mat3 would never be wrong.

I actually quite like the warning idea! Behavior is to use whatever they give us, but if it's numerically close to diagonal, a warning/suggestion that they could make it exactly diagonal for better performance may be warranted!

Lazersmoke · 2023-07-17T13:33:17Z

src/System/PairExchange.jl

@@ -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)


Feature request; flag to disable this warning. I use this to re-run the forwards problem with different J each time and currently no way to turn this off easily

Good idea. Per Slack discussion, let's address this in a follow-up PR using https://docs.julialang.org/en/v1/stdlib/Logging/

Lazersmoke · 2023-07-17T13:37:23Z

src/System/PairExchange.jl

+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.


Is this the (only) way people typically specify this? I can easily see dot(Si,B,Sj) where B is square root of b.

What you're suggesting is a general biquadratic interaction, which we don't support yet. There are currently no plans to support it in dipole mode because it would be fairly complicated (a generalization of the renormalization scheme here would need to be derived, https://arxiv.org/abs/2304.03874). The plan for these more complex interactions is to jump straight to SU(N) mode, which can be fully general. This way, people can build arbitrary spin polynomials, not restricted to just biquadratic.

Lazersmoke · 2023-07-17T13:39:58Z

Looks good! I think the type stuff (Float64 vs UniformScaling) can even be addressed at a later time, the interface probably shouldn't change based on that

Also some internal refactors that unify pair couplings.

A Union of isbits types cannot yet be stored on stack.

kbarros requested a review from Lazersmoke July 16, 2023 13:19

Lazersmoke reviewed Jul 17, 2023

View reviewed changes

kbarros added 4 commits July 18, 2023 09:22

Absorb set_biquadratic! into set_exchange!

59bed8b

Also some internal refactors that unify pair couplings.

Rename SingleIonAnisotropy to OnsiteCoupling

691e497

Avoid heap allocation of PairCoupling.bilin

eb76356

A Union of isbits types cannot yet be stored on stack.

Cleanup

d451774

kbarros force-pushed the general_interactions branch from c4c557d to d451774 Compare July 18, 2023 15:37

Update comment

40dbccc

kbarros merged commit 7865786 into main Jul 18, 2023
3 checks passed

kbarros deleted the general_interactions branch July 18, 2023 16:44

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

PairCoupling refactors #94

PairCoupling refactors #94

kbarros commented Jul 16, 2023

Lazersmoke Jul 17, 2023

kbarros Jul 17, 2023

Lazersmoke Jul 17, 2023

Lazersmoke Jul 17, 2023

kbarros Jul 17, 2023 •

edited

Loading

kbarros Jul 18, 2023

Lazersmoke Jul 17, 2023

Lazersmoke Jul 17, 2023

Lazersmoke Jul 17, 2023

kbarros Jul 17, 2023

kbarros Jul 18, 2023 •

edited

Loading

Lazersmoke Jul 17, 2023

kbarros Jul 17, 2023

Lazersmoke Jul 17, 2023

kbarros Jul 18, 2023

Lazersmoke Jul 18, 2023

kbarros Jul 18, 2023 •

edited

Loading

kbarros Jul 18, 2023

kbarros Jul 18, 2023

Lazersmoke Jul 18, 2023

Lazersmoke Jul 17, 2023

kbarros Jul 18, 2023

Lazersmoke Jul 17, 2023

kbarros Jul 17, 2023 •

edited

Loading

Lazersmoke commented Jul 17, 2023

PairCoupling refactors #94

PairCoupling refactors #94

Conversation

kbarros commented Jul 16, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kbarros Jul 17, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kbarros Jul 18, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kbarros Jul 18, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kbarros Jul 17, 2023 • edited Loading

Choose a reason for hiding this comment

Lazersmoke commented Jul 17, 2023

kbarros Jul 17, 2023 •

edited

Loading

kbarros Jul 18, 2023 •

edited

Loading

kbarros Jul 18, 2023 •

edited

Loading

kbarros Jul 17, 2023 •

edited

Loading